pythonscipymathematical-optimization

minimizing a multidimentional solution over a dataset


I have a rectangle within a 2d space and a set of points within the rectangle. The rectangle is moving in the following manner: the center moves a value of u on the x-axis, v on the y-axis and everything is scaled in a factor of sx and sy respectively. I only get the location of the points after the motion. My goal is to estimate the (u, v, sx, sy) vector.

I have the following data at runtime:

The equation to calculate the position of the point in the current frame given the (u, v, sx, sy) vector is given by:

x_currentFrameCalaulated = sx * px + u + x0 * (1 - sx)
y_currentFrameCalaulated = sy * py + v + y0 * (1 - sy)

Note that the axes are non-dependent. I have defined the following function to be minimized:

def estimateCurr(vec, previous, current, center):
    return np.abs(current - (vec[1]* previous + vec[0]+ center * (1 - vec[1])))

Here, vec[0] represents the motion concerning the axis, and vec[1] represents the scale. I am setting the bound according to my problem in the following manner:

h = 270
w = 480
boundsX = Bounds([-w, -1], [w, w - 1])
boundsY = Bounds([-h, -1], [h, h - 1])

Initializing a guess to [0, 0]:

init = np.zeros((2 , ))

I am resuming to try and find the optimal solution by:

res = minimize(estimateCurr, x0 = init, args=(px, cx, centerX), method='trust-constr',  bounds=boundsX, options={'verbose': 1})
print(res)
res = minimize(estimateCurr, x0 = init, args=(py, cy, centerY), method='trust-constr',  bounds=boundsY, options={'verbose': 1})
print(res)

I am getting:

ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 2 dimension(s) and the array at index 1 has 1 dimension(s)

This is obvious since I have 17 points, and the sizes of printing: print(px.shape) print(cx.shape) print(centerX.shape) print(init.shape)

(17,) (17,) (17,) (2,)

However, I am not sure how to set the correct sizes for the problem, or even if I am using the correct solver. I tried multiplying the (u, v, sx, sy) to fit the size of the data to no evil. How do I approach this problem?

My formal question is, given a set of measured datapoints, how do I fit a multidimensional solution that minimizes the error using Python?


Adding a dummy example concerning comments request. I have shortened the dataset to 6 points:

from scipy.optimize import Bounds
from scipy.optimize import minimize
import numpy as np

def estimateCurr(vec, previous, current, center):
    # np.full(())
    return np.abs(current - (vec[1] * previous + vec[0] + center * (1 - vec[1])))


h = 270
w = 480
x0 = 90.8021
y0 = -20.8282

px = np.array([86.7581, 74.5433, 85.0012, 84.348, 83.704, 91.6176])
py = np.array([-19.5163, -17.3714, -3.39899, -4.83069, -1.97073, -2.20099])
cx = np.array([89.7436, 75.8955, 87.5827, 87.1492, 86.0817, 92.6683])
cy = np.array([-19.2132, -16.3913, -2.9177, -4.81898, -1.49321, -2.43572])

numPoints = px.shape[0]

boundsX = Bounds([-w, -1], [w, w - 1])
boundsY = Bounds([-h, -1], [h, h - 1])
centerX = np.full((numPoints,), x0)
centerY = np.full((numPoints,), y0)
init = np.zeros((2 ,))
print(px.shape)
print(cx.shape)
print(centerX.shape)
print(init.shape)
res = minimize(estimateCurr, x0 = init, args=(px, cx, centerX), method='trust-constr',  bounds=boundsX, options={'verbose': 1})
print(res)
res = minimize(estimateCurr, x0 = init, args=(py, cy, centerY), method='trust-constr',  bounds=boundsY, options={'verbose': 1})
print(res)

Adding the full traceback:

--------------------------------------------------------------------------- ValueError                                Traceback (most recent call last) /tmp/ipykernel_6249/3484902536.py in <module>
     25 print(centerX.shape)
     26 print(init.shape)
---> 27 res = minimize(estimateCurr, x0 = init, args=(px, cx, centerX), method='trust-constr',  bounds=boundsX, options={'verbose': 1})
     28 print(res)
     29 res = minimize(estimateCurr, x0 = init, args=(py, cy, centerY), method='trust-constr',  bounds=boundsY, options={'verbose': 1})

~/.conda/envs/oxip-new6/lib/python3.7/site-packages/scipy/optimize/_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    634         return _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
    635                                             bounds, constraints,
--> 636                                             callback=callback, **options)
    637     elif meth == 'dogleg':
    638         return _minimize_dogleg(fun, x0, args, jac, hess,

~/.conda/envs/oxip-new6/lib/python3.7/site-packages/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py in _minimize_trustregion_constr(fun, x0, args, grad, hess, hessp, bounds, constraints, xtol, gtol, barrier_tol, sparse_jacobian, callback, maxiter, verbose, finite_diff_rel_step, initial_constr_penalty, initial_tr_radius, initial_barrier_parameter, initial_barrier_tolerance, factorization_method, disp)
    518             initial_barrier_tolerance,
    519             initial_constr_penalty, initial_tr_radius,
--> 520             factorization_method)
    521 
    522     # Status 3 occurs when the callback function requests termination,

~/.conda/envs/oxip-new6/lib/python3.7/site-packages/scipy/optimize/_trustregion_constr/tr_interior_point.py in tr_interior_point(fun, grad, lagr_hess, n_vars, n_ineq, n_eq, constr, jac, x0, fun0, grad0, constr_ineq0, jac_ineq0, constr_eq0, jac_eq0, stop_criteria, enforce_feasibility, xtol, state, initial_barrier_parameter, initial_tolerance, initial_penalty, initial_trust_radius, factorization_method)
    306         barrier_parameter, tolerance, enforce_feasibility, ...
    347 

<__array_function__ internals> in concatenate(*args, **kwargs)

ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 2 dimension(s) and the array at index 1 has 1 dimension(s)

Solution

  • The objective function must be a scalar function. In your definition the function returns array. So you need to aggregate it in some way. The most obvious choice is generally the sum.

    For instance the following objective function works and converges:

    def estimateCurr(vec, previous, current, center):
        return np.sum(np.power(current - (vec[1] * previous + vec[0] + center * (1 - vec[1])), 2))
    
    #     message: `xtol` termination condition is satisfied.
    #     success: True
    #      status: 2
    #         fun: 3.1694871819514345
    #           x: [ 2.274e+00  1.013e+00]
    # ...