[FIXED] solve rectangular matrix in python to get solution with arbitrary parameters

Issue

I want to solve a rectangular system (with arbitrary parameters in the solution). Failing that I would like to add rows to my matrix until it is square.

print matrix_a

print vector_b

print len(matrix_a),len(matrix_a[0])

gives:

 [[0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0],
 [0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0],
 [0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0],
 [0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
 [0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0],
 [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0],
 [1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]]

[2, 0, 2, 1, 2, 1, 1, 1, 1, 1, 1]

11 26

my full code is at http://codepad.org/tSgstpYe

As you can see I have a system Ax=b. I know that each x value x1,x2.. has to be either 1 or 0 and I expect with this restriction that there should be only one solution to the system.

In fact the answer I expect is x=[0,1,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,0]

I looked at http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html#numpy.linalg.solve but it appears to only take square matrices.

Any help on solving this system would be great!

Solution

Here is a simple implementation (with hard coded thresholds), but it gives the solution you are looking for with the test data.

It's based on Iteratively Reweighted Least Squares.

from numpy import abs, diag, dot, zeros
from numpy.linalg import lstsq, inv, norm

def irls(A, b, p= 1):
    """Solve least squares problem min ||x||_p s.t. Ax= b."""
    x, p, e= zeros((A.shape[1], 1)), p/ 2.- 1, 1.
    xp= lstsq(A, b)[0]
    for k in xrange(1000):
        if e< 1e-6: break
        Q= dot(diag(1./ (xp** 2+ e)** p), A.T)
        x= dot(dot(Q, inv(dot(A, Q))), b)
        x[abs(x)< 1e-1]= 0
        if norm(x- xp)< 1e-2* e** .5: e*= 1e-1
        xp= x
    return k, x.round()

Answered By - eat