python - 如何获得与 Matlab 的 'special' (mldivide) 运算符使用 numpy/scipy 返回的欠定线性系统相同的 `A\b` 解？

Question

我找到了一个 link其中显示了一个示例，当线性方程组有无限多个解时，Matlab mldivide 运算符 (\) 给出“特殊”解。

例如:

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

给出输出:

c_mldivide =
                 0
                 4
  0.66666666666667


c_pinv =

  0.918032786885245
  3.54098360655738
  1.27868852459016

在解决方案 c_mldivide 中非零项的数量等于 rank(A) 的意义上，该解决方案是“特殊的”(在本例中为 2 ).我使用 numpy.linalg.lstsq 在 numpy 中尝试了同样的事情，它给出了与 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

非负最小二乘法 ( scipy.optimize.nnls ) 不是此问题的通用解决方案。如果所有可能的解决方案都包含负系数，那么它会失败的一个简单情况是:

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 未定的情况下，

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

将选择一个解决方案 x' 最小化 ||x||_L2 受 ||A·x - b||_L2 = 0。这恰好不是我们正在寻找的特定解决方案，但我们可以对其进行线性变换以获得我们想要的结果。为此，我们将首先计算 right null space的 A，表征 A·x = 的所有可能解的空间b。我们可以使用 rank-revealing QR decomposition 得到它:

from scipy.linalg import qr

def qr_null(A, tol=None):
    Q, R, P = qr(A.T, mode='full', pivoting=True)
    tol = np.finfo(R.dtype).eps if tol is None else tol
    rnk = min(A.shape) - np.abs(np.diag(R))[::-1].searchsorted(tol)
    return Q[:, rnk:].conj()

Z = qr_null(A)

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

如果 x' 中的最后一个 n - rnk(A) 系数为零，则此强加于此

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

因此我们可以求解 C(准确地说，因为我们知道 Z[rnk:] 必须是满秩):

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)