solve rectangular matrix in python to get solution with arbitrary parameters

Question

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!

Answer 1

as alternative to the accepted answer given as "implementation (with hard coded thresholds)" - solution in one line (without data):

if rounding to integer numbers is allowed (as OP require in comments) & inspired by comments about the best here in least-squares discussion, tried scipy.optimize.nnls (least squares with non-negativity constraint)

import numpy as np
from scipy.optimize import nnls

A = np.array([[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]])

b = np.array([2, 0, 2, 1, 2, 1, 1, 1, 1, 1, 1])    

##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]

x, rnorm = nnls(A,b)
print( x.round())
## OK: [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.]

P.S. regularization to Normal Equation is still under consideration (though according to the OP - only integer numbers accepted as answer)... or for nnls

P.S. if "solution must contain only integer numbers" - can see libs for Mixed-Integer Linear Programming problems - such as GEKKO, pulp, scipy.optimize.milp, or scipy.optimize.minimize with gurobi/mip solver [n.b. differences between MILP & MIP] & use objective fuction minimizing MSE, as for least-squares technique

p.s. or can try milp

Answer 2

for rectangular matrix better use least-squares method (1, ) - IRLS looks like this for Normal Equation - with matix gradient descent:

import numpy as np

A = np.array([[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]])

b = np.array([2, 0, 2, 1, 2, 1, 1, 1, 1, 1, 1])

######################
##https://people.duke.edu/~ccc14/sta-663-2020/notebooks/S09F_Least_Squares_Optimization.html

# Matrix form of gradient
def grad_m(X, y, b):
    return X.T@X@b- X.T@y

max_iter = 1000
a1 = 0.01 # learning rate
beta2 = np.zeros(26)
for i in range(max_iter):
    beta2 -= a1 * grad_m(A, b, beta2)
print(np.round(beta2,0))

###################
# same as
##from scipy.linalg import lstsq
##p, res, rnk, s = lstsq(A, b)
##print(p.round())

still the result [ 0. 1. 0. -0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 1. -0. 0. 0. 1. 0. -0. 0. 0. -0. 0. 0. 0.] differs from yours.

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]

Perhaps, the reason is that your Linear System is underdetermined (many x, less equations) - no single solution, thus different methods leads to different solutions, I suppose. OR perhaps your problem is Non-Convex? (the reason why GD is not working correct)

p.s. really, method with Regularization (recovered link), suggested above @eat, or here p.53, or here p29, regularizes for Matrix Rank Minimization for approximate solution of system, when matrix is not square or determinant=0

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

Answer 4

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.

Answer 5

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.