I wish to use the scipy.optimize.leastsq() method from scipy to optimize three parameters a,b,c. What I do have are these two equations.
1*a+2*b+3*c = x1
4*a+5*b+6*c = x2
Analytically this set of equations is underdetermined, but numerically I am trying to find a,b,c to minimize the error to given results from a measurement [2,2]:
1*a+2*b+3*c - 2 = 0
4*a+5*b+6*c - 2 = 0
Therefore I have written some code.
def function(a,b,c,t):
return np.array([1*a+2*b+3*c+t[1],4*a+5*b+6*c+t[1]])
a0 = 1
b0 = 1
c0 = 1
measdata = np.array([2,2])
t = [1,2]
def residual(x0,measdata,t):
return measdata - function(x0[0],x0[1],x0[2],t)
erg = optimize.leastsq(func=residual,x0=(a0,b0,c0),args=(measdata,t))
It always results in:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-296-ab0fc90a2253> in <module>()
14 return result - function(x0[0],x0[1],x0[2],t)
15
---> 16 erg = optimize.leastsq(func = residual, x0 = (a0,b0,c0) , args=(result,t), maxfev=10000)
17
18 function(erg[0][0],erg[0][1])
//anaconda/lib/python3.5/site-packages/scipy/optimize/minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
378 m = shape[0]
379 if n > m:
--> 380 raise TypeError('Improper input: N=%s must not exceed M=%s' % (n, m))
381 if epsfcn is None:
382 epsfcn = finfo(dtype).eps
TypeError: Improper input: N=3 must not exceed M=2
How do I get it to find a minimum? I am aware it's only a local minimum, but I would be happy with it.
The error is telling you what you already know, i.e., the system is underdetermined, with n being the number of parameters and m the number of constraints.
If you fix one of the parameters so that n > m is False, the code will stop complaining. For example, change
def residual(x0,measdata,t):
return measdata - function(x0[0],x0[1],x0[2],t)
erg = optimize.leastsq(func=residual,x0=(a0,b0,c0),args=(measdata,t))
to
def residual(x0,measdata,t):
# we fix the value of `c` here
return measdata - function(x0[0],x0[1],5,t)
# only two parameters for `x0`
erg = optimize.leastsq(func=residual,x0=(a0,b0),args=(measdata,t))
To answer the question how you can do what you want, I'm not sure it can be done with scipy. I found this issue saying that scipy can't handle underdetermined systems:
Interesting, I assumed the MINPACK routines would handle also m < n, but apparently not. The reason why they don't is probably that for m < n the minimum is some manifold of points, which causes problems in the termination conditions.
So would be after all some interest in adding also a small-scale solver for underdetermined least-square problems.
Even though that post was from 3 years ago, I still can't find any evidence in the docs that scipy can do what you want. However, I found a SO answer that claims you can solve for an underdetermined matrix, but I haven't fully grasped the mathematics yet to make sure whether it'd be applicable for your case. Since I find it hard to summarise the post, I'll just quote what seems to be the most important section.
In the case where A·x = b is underdetermined,
x1, res, rnk, s = np.linalg.lstsq(A, b)will pick a solution x' that minimizes ||x||L2 subject to ||A·x - b||L2 = 0. This happens not to be the particular solution we are looking for, but we can linearly transform it to get what we want. In order to do that, we'll first compute the right null space of A, which characterizes the space of all possible solutions to A·x = b. We can get this using a rank-revealing QR decomposition