SOLVED: scipy.optimize.least_squares problem

[Skytter13]

Hi,
I'm trying to make a code to solve a overdetermined system of non-linear equations.
It is to be used to locate positions of instruments on sea bottom: for each instrument there are several surface sources with known coordinates (Xs, Ys, Zs) and the time the signal propagated to the instrument Ts, so I want to find location of the instrument (Xi, Yi, Zi) and water velocity (Vw).
I have a large number of points for which I need to run it, but I'm stuck at modifying the code from single to multiple points. Basically I want to be able to set equation coefficients from the array and then loop solution over it.
Any advice will be appreciated.

here is my code for single point solution:

import numpy as np
import scipy.optimize

# observed values
Ts = np.array([0.0025, .002, .0039, .0041, .0036, .0015])


def f6(p):  # 6 equations to solve with hard coded sources coordinates (Xs, Ys, Zs=0) ; p[0:2] - Xi, Yi, Zi, p[3] - Vw
     return np.array([
        ((p[0] - 0) ** 2 + (p[1] - 0) ** 2 + (p[2] - 0) ** 2) / p[3]**2,
        ((p[0] - 0) ** 2 + (p[1] - 5) ** 2 + (p[2] - 0) ** 2) / p[3]**2,
        ((p[0] - 6) ** 2 + (p[1] - 5) ** 2 + (p[2] - 0) ** 2) / p[3]**2,
        ((p[0] - 6) ** 2 + (p[1] - 0) ** 2 + (p[2] - 0) ** 2) / p[3]**2,
        ((p[0] - 6) ** 2 + (p[1] - 3) ** 2 + (p[2] - 0) ** 2) / p[3]**2,
        ((p[0] - 0) ** 2 + (p[1] - 3) ** 2 + (p[2] - 0) ** 2) / p[3]**2])


def system(p):  # Returns the residuals
    return f6(p) - Ts**2


res = scipy.optimize.least_squares(system, [1, 1, 1, 1450],
                                   jac='2-point', loss='soft_l1', gtol=5e-16, verbose=2)
 

This code works as expected. But when I try to move to parameter definition of system of equations is fails.

import numpy as np
import scipy.optimize

# observed values
Ts = np.array([0.0025, .002, .0039, .0041, .0036, .0015])

# array with coefficients for equations  scoor[0] - Xs, scoor[1] - Ys, Zs=0
scoor = np.array([[0, 0, 6, 6, 6, 0],
                 [0, 5, 5, 0, 3, 3]])

def f6s(p):  # making 6 equations to solve, using scoor coefficients; p[0:2] - Xi, Yi, Zi, p[3] - Vw
        return np.asarray([((p[0] - scoor[0]) ** 2 + (p[1] - scoor[1]) ** 2 + (p[2] - 0) ** 2) / p[3]**2])

def system1(p):  Returns the residuals
    return f6s(p) - Ts**2

res2 = scipy.optimize.least_squares(system1, [1, 1, 1, 1450],
                                   jac='2-point', loss='soft_l1', gtol=5e-16, verbose=2)

This does not work (( and I can't figure out what is wrong ...

Error:
Traceback (most recent call last):
  File "C:\ProgramData\Anaconda3\envs\tomo\lib\site-packages\IPython\core\interactiveshell.py", line 3444, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-19-979c61a9e4d5>", line 29, in <module>
    res2 = scipy.optimize.least_squares(system1, [1, 1, 1, 1450],
  File "C:\ProgramData\Anaconda3\envs\tomo\lib\site-packages\scipy\optimize\_lsq\least_squares.py", line 823, in least_squares
    raise ValueError("`fun` must return at most 1-d array_like. "
ValueError: `fun` must return at most 1-d array_like. f0.shape: (1, 6)

What do I have wrong? Thanks in advance

[Gribouillis]

The error message is clear. The 'least_squares' function expects 'system1' to return a 1-d array_like. Looking in the documentation, it is said that

The argument x passed to this function is an ndarray of shape (n,) (never a scalar, even for n=1). It must allocate and return a 1-D array_like of shape (m,) or a scalar.

I suggest that you check the shape of the value returned by system1() and raise an exception if it is not a tuple with one element like (m,) to understand the problem.

[Skytter13]

The error message is clear. The 'least_squares' function expects 'system1' to return a 1-d array_like. Looking in the documentation, it is said that

The argument x passed to this function is an ndarray of shape (n,) (never a scalar, even for n=1). It must allocate and return a 1-D array_like of shape (m,) or a scalar.

I suggest that you check the shape of the value returned by system1() and raise an exception if it is not a tuple with one element like (m,) to understand the problem.

Thanks for the prompt reply. Solved, the solution was to flatten equation array:

def f6s(p):  # making 6 equations to solve, using scoor coefficients; p[0:2] - Xi, Yi, Zi, p[3] - Vw
        return np.asarray([((p[0] - scoor[0]) ** 2 + (p[1] - scoor[1]) ** 2 + (p[2] - 0) ** 2) / p[3]**2]).flatten()