Help with solving the system of N linear equations

[Roo Riop]

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 error message is telling you that you have 3 variables and only one value in the residual.  That's an ill-posed problem for non-linear least squares.  The basic answer is to not take the norm, but return the residual array, `np.dot(A, u) - B`.

With 3 parameters and 3 data values, the result might be somewhat suspicious, but at least it runs to completion:

import numpy as np
from lmfit import Parameters, Minimizer, fit_report

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)

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 np.dot(A,u)-B   # return the residual array, not the norm

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()

print(fit_report(result))