[Tutor] constrained least square fitting using python

[eryksun]

On Wed, Aug 7, 2013 at 5:37 PM, sudipta <sudipta.mml at gmail.com> wrote:
>
> 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. Further,
> there is a equality constraint on P which is Sum(P(i))=0.0. How do I proceed
> to solve that? 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 can use scipy.optimize.minimize() with method='SLSQP' (it uses
fmin_slsqp). You'll need an objective function that returns the sum of
the squared residuals; the algorithm minimizes this. You'll also need
an 'eq' constraint function; the algorithm tries to make this zero.

Here's a contrived example. Hopefully if there's a problem with this
example, someone else on the list (Oscar?) can set it straight. I
don't use scipy.optimize very much. I know a bit more about
scipy.signal.

This example uses a 100x3 matrix X that's set up to compute a
polynomial fit. Normally I'd use the pseudoinverse of the matrix, but
in this case there's an added constraint. Also, strictly speaking the
matrix isn't required when using minimize(); you could just evaluate
the polynomial to compute the residuals at each step, but this way
matches your problem description.

    import numpy as np
    from scipy import optimize
    import matplotlib.pyplot as plt

    np.random.seed(42)

    def peval(p, X):
        return np.dot(X, p)

    def objective(p, y, X):
        return np.sum((y - peval(p, X)) ** 2)

    # y = -2*x**2 + 2.5*x - 0.5
    t = np.reshape(np.arange(100) * 0.1, (100, 1))
    X = np.hstack([t**2, t**1, t**0])
    ptrue = np.array([-2.0, 2.5, -0.5])
    ytrue = peval(ptrue, X)
    # add measurement noise
    ymeas = ytrue + 5 * np.random.randn(len(ytrue))

    p0 = np.zeros(3)
    cons = {
        'type': 'eq',
        'fun': np.sum,
    }

    res = optimize.minimize(
        objective, p0, args=(ymeas, X),
        method='SLSQP', constraints=cons)

    pcons = res.x
    yfit = peval(pcons, X)

    ss_res = np.sum((ymeas - yfit) ** 2)
    ss_tot = np.var(ymeas) * len(ymeas)
    rsq = 1 - ss_res / ss_tot

    plt.plot(t, ymeas, 'o', t, ytrue, t, yfit)
    plt.title('Least Squares Fit')
    plt.legend(['Noisy Data', 'True', 'Fit'])
    plt.show()


Graph:
http://i.imgur.com/BHLPbCa.png

Estimated parameters and R**2 value:

    >>> pcons
    array([-1.94471575,  1.95757699, -0.01286124])
    >>> np.sum(pcons)
    -2.0816681711721685e-17
    >>> rsq
    0.99247909721611327

Pseudoinverse result for comparison:

    >>> plsq = np.dot(np.linalg.pinv(X), ymeas)
    >>> plsq
    array([-1.96260467,  2.19944918, -0.75938176])
    >>> np.sum(plsq)
    -0.52253724312209116

    >>> ss_res = np.sum((ymeas - np.dot(X, plsq)) ** 2)
    >>> 1 - ss_res / ss_tot
    0.99250546279261975