在Python中求解矩形矩阵以获得任意参数的解

问题描述 投票：0回答：4

我想求解一个矩形系统（解中具有任意参数）。如果失败，我想向矩阵添加行，直到它成为正方形。

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 但它似乎只需要方阵。

任何解决这个系统的帮助都会很棒！

python matrix numpy scipy linear-algebra
4个回答
4
投票

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

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

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()
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 万次测试。

1
投票

设 

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

有3种可能

  1. 没有解决方案，如果
    rank(A) < rank(A|b)
  2. 只有一种解决方案当且仅当
    rank(A)  = rank(A|b) = n
  3. 无限多个解决方案 iff 
    rank(A) = rank(A|b) < n

其中 

n
是未知数。

0
投票

对于矩形矩阵，最好使用最小二乘法 - IRLS 看起来像正态方程 - 具有矩阵梯度下降：

import numpy as np
from numpy.linalg import lstsq, inv, norm

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

结果仍然

[ 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.]
与你的不同。

也许，原因是你的线性系统不确定（许多x，更少的方程）-没有单一的解决方案，因此我认为不同的方法会导致不同的解决方案

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.