对于欠定线性系统，如何获得与Matlab的“A\b”（mldivide）运算符使用numpy/scipy返回的“特殊”解？

Question

我发现了一个link，其中举例说明了当线性方程组有无穷多个解时，Matlab mldivide算子（\）给出了“特殊”解。在

例如：

A = [1 2 0; 0 4 3];
b = [8; 18];
c_mldivide = A \ b
c_pinv = pinv(A) * b

给出输出：

^{pr2}$

解决方案是“特殊”的，因为解决方案c_mldivide中的非零条目的数目等于rank(A)（在本例中为2）。我在numpy中使用numpy.linalg.lstsq尝试了相同的操作，这给了c_pinv相同的结果。在

有没有一种方法可以在Python中实现c_mldivide解决方案？在

还有一个非常相似的问题here，但我想“特殊”这个词的解释不够清楚。 Another question询问了mldivide运算符的内部工作方式，但接受的答案似乎没有解决这个问题。在

编辑1:numpy代码

In [149]: test_A = np.array([[1,2,0],[0,4,3]])
          test_b = np.array([[8],[18]])
          np.linalg.lstsq(test_A,test_b)

Out[149]:
(array([[ 0.918 ],
    [ 3.541 ],
    [ 1.2787]]), array([], dtype=float64), 2, array([ 5.2732,  1.4811]))

编辑2:使用scipy.optimize.nnls

In[189]:
from scipy.optimize import nnls
nnls(test_A,test_b)
Out[190]:
    ValueError                                Traceback (most recent call last)
<ipython-input-165-19ed603bd86c> in <module>()
      1 from scipy.optimize import nnls
      2 
----> 3 nnls(test_A,test_b)

C:\Users\abhishek\Anaconda\lib\site-packages\scipy\optimize\nnls.py in nnls(A, b)
     43         raise ValueError("expected matrix")
     44     if len(b.shape) != 1:
---> 45         raise ValueError("expected vector")
     46 
     47     m, n = A.shape

    ValueError: expected vector

Answer 1

非负最小二乘（^{}）是这个问题的一般解。如果所有可能的解决方案都包含负系数，则会失败：

import numpy as np
from scipy.optimize import nnls

A = np.array([[1, 2, 0],
              [0, 4, 3]])
b = np.array([-1, -2])

print(nnls(A, b))
# (array([ 0.,  0.,  0.]), 2.23606797749979)

---

在A·x=b的情况下

^{pr2}$

将选择一个最小化L2的解决方案。这恰好不是我们要寻找的特定解决方案，但我们可以线性变换它以得到我们想要的。为了做到这一点，我们首先计算A的right null space，它表征了A·x=b的所有可能解的空间。我们可以使用rank-revealing QR decomposition：

^{3}$

Z是一个向量（或者，在n-rnk（a）>；1的情况下，跨越a的子空间的一组基向量），使得a·Z=0：

print(A.dot(Z))
# [[  0.00000000e+00]
#  [  8.88178420e-16]]

换句话说，Z的列是与A中所有行正交的向量。这意味着对于任何解x'到A·x=b，那么x'=x+Z·c也必须是任意标度因子c的解。这意味着，通过选取适当的c值，我们可以将解中任何系数的n-rnk（A）设置为零。在

例如，假设要将最后一个系数的值设置为零：

c = -x1[-1] / Z[-1, 0]
x2 = x1 + Z * c
print(x2)
# [ -8.32667268e-17  -5.00000000e-01   0.00000000e+00]
print(A.dot(x2))
# [-1. -2.]

---

更一般的情况是n-rnk（A）≤1则稍微复杂一些：

A = np.array([[1, 4, 9, 6, 9, 2, 7],
              [6, 3, 8, 5, 2, 7, 6],
              [7, 4, 5, 7, 6, 3, 2],
              [5, 2, 7, 4, 7, 5, 4],
              [9, 3, 8, 6, 7, 3, 1]])
x_exact = np.array([ 1,  2, -1, -2,  5,  0,  0])
b = A.dot(x_exact)
print(b)
# [33,  4, 26, 29, 30]

我们得到x'和Z：

x1, res, rnk, s = np.linalg.lstsq(A, b)
Z = qr_null(A)

现在为了使解向量中的零值系数的数目最大化，我们要找到一个向量C，这样

x' = x + Z·C = [x'₀, x'₁, ..., x'_rnk(A)-1, 0, ..., 0]^T

如果n-rnk（A）中的最后一个n-rnk（A）系数为零，则这就要求

Z_{{rnk(A),...,n}}·C = -x_{{rnk(A),...,n}}

因此，我们必须清楚地知道（因为我们可以精确地求解）

C = np.linalg.solve(Z[rnk:], -x1[rnk:])

并计算x'：

x2 = x1 + Z.dot(C)
print(x2)
# [  1.00000000e+00   2.00000000e+00  -1.00000000e+00  -2.00000000e+00
#    5.00000000e+00   5.55111512e-17   0.00000000e+00]
print(A.dot(x2))
# [ 33.   4.  26.  29.  30.]

---

要将它们组合成一个功能：

import numpy as np
from scipy.linalg import qr

def solve_minnonzero(A, b):
    x1, res, rnk, s = np.linalg.lstsq(A, b)
    if rnk == A.shape[1]:
        return x1   # nothing more to do if A is full-rank
    Q, R, P = qr(A.T, mode='full', pivoting=True)
    Z = Q[:, rnk:].conj()
    C = np.linalg.solve(Z[rnk:], -x1[rnk:])
    return x1 + Z.dot(C)

Answer 2

np.array([[8],[18]]).shape

是

^{pr2}$

但你想

^{3}$

输出：

expected vector
(array([ 0.        ,  4.        ,  0.66666667]), 0.0)