3

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!

Rusty Rob
  • 16,489
  • 8
  • 100
  • 116
  • I don't think that linear algebra is the right tool for this job - it will only tell you that there is an infinite number of solutions and the area containing these solutions. Linear algebra is pretty useless if your condition is "the solution must contain only integer numbers". – Wladimir Palant Aug 16 '11 at 10:34
  • yes I'm expecting infinite solutions. I expect that there will be a unique point where the infinite solutions cross the unit-ball defined by the one norm in the 1st quadrant which is the solution i'm looking for. http://en.wikipedia.org/wiki/Unit_sphere#Unit_balls_in_normed_vector_spaces – Rusty Rob Aug 16 '11 at 10:44
  • I can very easily write my own script to get this matrix into row echelon form and then do back substitution but I figured this must already be out there. – Rusty Rob Aug 16 '11 at 10:47
  • http://mail.scipy.org/pipermail/numpy-discussion/2008-November/038705.html and http://stackoverflow.com/questions/5889142/python-numpy-scipy-finding-the-null-space-of-a-matrix I'll check this out in the morning – Rusty Rob Aug 16 '11 at 10:50

3 Answers3

3

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()
eat
  • 7,440
  • 1
  • 19
  • 27
  • Thank you, this is extremely powerful and is exactly what I was looking for. The only problem is it didn't quite find the right solution on a large system http://codepad.org/6DUD7Wrf (see bottom of code comment) Which hard coded thresholds do I need to change to get the solution for this larger system. – Rusty Rob Aug 17 '11 at 11:46
  • @robert king: Well, first you may like to increase the maximum number iterations (in `xrange(.)`). But eventually you may need to twist all of the thresholds and norm specification (`p`) :(. Please note that `irls(.)` is by no mean a 'silver bullet' kind of solution for this kind of problems. Anyway, for much more advanced details, please study the link I provided in my answer. References therein may also give you more specific guidance how to 'fine tune' `irls(.)` proper manner, in practice. Thanks – eat Aug 17 '11 at 18:30
  • Thanks, I've marked your solution as correct and ill fine tune it when I have time to read the article cheers. – Rusty Rob Aug 17 '11 at 20:16
2

Depending on the input you expect you might be better off with a simple tree search algorithm. Your result vector contains relatively low numbers which allows cutting off most tree branches early. My attempt at implementing this algorithm produces the expected result after 0.2 seconds:

def testSolution(a, b, x):
  result = 0
  for i in range(len(b)):
    n = 0
    for j in range(len(a[i])):
      n += a[i][j] * x[j]
    if n < b[i]:
      result = -1
    elif n > b[i]:
      return 1
  return result

def solve(a, b):
  def solveStep(a, b, result, step):
    if step >= len(result):
      return False

    result[step] = 1
    areWeThere = testSolution(a, b, result)
    if areWeThere == 0:
      return True
    elif areWeThere < 0 and solveStep(a, b, result, step + 1):
      return True
    result[step] = 0
    return solveStep(a, b, result, step + 1)

  result = map(lambda x: 0, range(len(a[0])))
  if solveStep(a, b, result, 0):
    return result
  else:
    return None

matrix_a = [[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]]

vector_b = [2, 0, 2, 1, 2, 1, 1, 1, 1, 1, 1]

print solve(matrix_a, vector_b)

This had to test only 1325 possible vectors with your input which is a lot less than all the possible results (67 million). Worst-case scenario is still 67 million tests of course.

Wladimir Palant
  • 56,865
  • 12
  • 98
  • 126
  • Thanks. I was going to try something like this. I'll give this a go when I get home (albeit on a much larger matrix). If that works I'll mark your answer as correct. – Rusty Rob Aug 16 '11 at 22:28
  • Thanks. I forgot to mention I have larger systems that this is too slow on. Thanks anyway, this works on the smaller ones so +1 – Rusty Rob Aug 17 '11 at 11:08
1

Let Ax = b be the system and A|b be augmented matrix then

There are 3 possibilities

  1. No solutions iff rank(A) < rank(A|b)
  2. Only one solution iff rank(A) = rank(A|b) = n
  3. Infinite number of solutions iff rank(A) = rank(A|b) < n

where n is number of unknowns.

Pratik Deoghare
  • 35,497
  • 30
  • 100
  • 146
  • Thanks. I therefore expect infinite solutions. However the restriction is that of these solutions I want one that is distance one away from the origin using the one norm. – Rusty Rob Aug 16 '11 at 22:24
  • 1
    @robert: I don't think that the one norm will help you define your solution space. ||x|| is actually 5 in your example when using that norm. – Wladimir Palant Aug 17 '11 at 05:40
  • yes sorry, i think i meant the inf norm. basically I just wanted x to be a vector with 1's and 0's in it that satisfied the matrix – Rusty Rob Aug 17 '11 at 11:51