[Tutor] constrained least square fitting using python

[Oscar Benjamin]

On 7 August 2013 22:37, sudipta <sudipta.mml at gmail.com> wrote:
> Hi All,
>
> I am facing a problem for constrained linear least square fitting. In my
> case the matrix equation looks like [Y]nX1=[X]nXm[P]mX1, where Y and P are
> vectors and X is a matrix and n, m are dimension of the matrix.

It took me a moment to understand that equation. I would have written it as

    y = X p

where I use capitals to distinguish matrices and vectors.

Okay so to change the terminology slightly if

    E[y] = X p

is the estimate of y given X and p. We want to choose p_lsq to minimise

    res = sum((E[y] - y)^2) = sum((X p - y)^2)

The optimal solution is given by the exact solution of the linear system

    X' X p_lsq = X' y

where X' means the transpose of X.

> Further,
> there is a equality constraint on P which is Sum(P(i))=0.0. How do I proceed
> to solve that?

The least-squares choice p_lsqcon that satisfies the linear constraint

    A p = a

is given by

    p_lsqcon = p_lsq - Z A' [A Z A']^(-1) (A p - a)

where Z is the inverse of X' X. In general A is a matrix but in your
case it is a 1xm vector of all 1s and a is just 0. The expressions
above may look complicated but when you actually work it out they're
not so bad.

> Which function of python is suitable for this? I saw few of
> discussion on scipy.optimize.fmin_slsqp() function but the implementation of
> this function is not very straightforward. Therefore, I need your help. I am
> new in SCIPY. Please help me out in this regard.

You should look at how to do matrix multiplication and how to solve a
linear system. Here's a quick ipython session that illustrates how to
do this:

In [1]: import numpy as np

In [2]: M = np.array([[1, 2], [3, 4]])

In [3]: M
Out[3]:
array([[1, 2],
       [3, 4]])

In [4]: x = np.array([[1], [2]])

In [5]: x
Out[5]:
array([[1],
       [2]])

In [6]: np.dot(M, x)  # matrix multiplication
Out[6]:
array([[ 5],
       [11]])

In [7]: b = np.dot(M, x)

In [8]: b
Out[8]:
array([[ 5],
       [11]])

In [9]: np.linalg.solve(M, b)
Out[9]:
array([[ 1.],
       [ 2.]])

I don't know whether using this analytic result will be more/less
accurate/efficient than the methods that Eryksun suggested but this
would be my first attempt.


Oscar