用python求解矩形矩阵，得到任意参数的解

Question

我想解一个矩形系统(解中有任意参数)。如果失败，我想添加行到我的矩阵，直到它是正方形。

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

给:

 [[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

我的完整代码在http://codepad.org/tSgstpYe

你可以看到我有一个系统Ax=b。我知道每个x值x1 x2。必须是1或0，并且我期望在这个限制下，系统应该只有一个解。

事实上我希望答案是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)

我看了看http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html#numpy.linalg.solve但它似乎只取方阵

任何帮助解决这个系统将是伟大的!

Answer 1

这是一个简单的实现(带有硬编码阈值)，但它提供了您正在寻找的测试数据的解决方案。

它基于迭代加权最小二乘。

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 2

根据您期望的输入，您可能更适合使用简单的树搜索算法。你的结果向量包含相对较低的数字，这允许早期切断大多数树枝。我尝试实现这个算法，在0.2秒后产生预期的结果:

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)

这只需要用你的输入测试1325个可能的向量，这比所有可能的结果(6700万)少得多。最坏的情况当然还是6700万次测试

Answer 3

设Ax = b为系统，A|b为增广矩阵，则

有三种可能

1. rank(A) < rank(A|b)
无解
2. rank(A) = rank(A|b) = n
只有一个解
3. rank(A) = rank(A|b) < n无限解

其中n为未知数