Help with solving the system of N linear equations

[rc]

Hi,

I have a simple problem which I can approach with scipy minimize, but I would like to use lmfit due to the ease of use and additional options (bounds and constraints on parameters, and most important getting the error bars on estimated parameters).

Something which was also mentioned in M Newville's reply here:

https://stackoverflow.com/questions/43593592/errors-to-fit-parameters-of-scipy-optimize

I have a set of N linear equations (i.e., 3), which I would like to solve, for example:

sua = [0.02, 0.05, 55.0]
sub = [99.7, 99.5, 32.0]
suc = [0.02, 0.03, 0.4]
analysis_results = [34.14, 60.8, 5.0]

A = np.array([sua, sub, suc])
B = np.array(analysis_results)

With scipy, I would do this using:

fun = lambda x: np.linalg.norm(np.dot(A,x)-B)

scipy.optimize.minimize(fun, np.zeros(len(B)), method='L-BFGS-B', bounds=[(0.,100.) for x in range(len(B))])

However, I am having troubles transforming this into lmfit syntax. I tried to google some example with a simple N linear equations solving using lmfit, but could not find any.

Currently, I have the following code: 

def fcn2min(params, A, B):
    u1 = params['u1']
    u2 = params['u2']
    u3 = params['u3']
    u = np.array([u1, u2, u3])
    model = np.linalg.norm(np.dot(A,u)-B)
    return model

params = Parameters()

params.add('u1', value=45.)

params.add('u2', value=10.)

params.add('u3', value=45.)

minner = Minimizer(fcn2min, params, fcn_args=(A, B))

result = minner.minimize()

but it does not work properly, returning: 

TypeError: Improper input: N=3 must not exceed M=1

I have a feeling that my fcn2min is not properly defined for this problem. I have used lmfit in the past for curve fitting and ODE solving, but never for solving the system of N linear equations. 

Can you please suggest a solution ?

Many thanks

[Matt Newville]

That looks mostly OK to me, though it should be stated clearly that you don't really need an iterative method to solve a linear problem.

Your use of lmfit and scipy.optimize are pretty different though.  With scipy.optimize, you give initial values for parameters as 0 and set bounds for all of them between 0 and 100.  As a side note, it is a really bad idea to put an initial value at one of the bounds, and that is especially true when that value is 0 (ie, "how zero do you mean":  0.0001 or 1e-120?).  The code will not know what order of magnitude you mean.

With lmfit, you give non-zero initial values (and very far from "0"), and no bounds at all.

With scipy.optimize you specify using the  'lbfgs-b' solver.  With lmfit, you do not,  which tells it to use the default solver, 'leastsq'.  That default solver (sensibly!) requires there to be more observations than variables.  But your objective function is returning a scalar value -- the norm of  (A.x-B).   That is the meaning behind the error message you are getting: N (the number of variables) must not exceed the number of observations (M) used to optimize the variables.  

You could try using a scalar solver, including the 'lbfgs-b' solver, in lmfit.  Or you could solve 'Ax-B' in the least-squares sense, changing your objective function return just `np.dot(A, u) - B`.

--Matt