我如何才能获得与Matlab的A \ b（mldivide）运算符使用numpy / scipy返回的欠定线性系统相同的“特殊”解决方案？

Question

我找到了一个链接，并在其中举例说明了当线性方程组具有无限多个解时，Matlabmldivide运算符（\）提供“特殊”解。

例如：

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中使用尝试了相同的操作numpy.linalg.lstsq，它给带来了相同的结果c_pinv。

有没有c_mldivide办法用Python实现解决方案？

还有一个非常类似的问题在这里，但我想这个词“特殊”的解释是不够清晰。
另一个问题询问了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_ 。碰巧这并不是我们要寻找的特定解决方案，但是我们可以对其进行线性转换以得到我们想要的。为了做到这一点，我们首先计算正确的零空间的
一个 ，表征所有可能的解决方案的空间 A·X = b
。我们可以通过显示排名的QR分解获得此信息： _
_ _
_
* __ __*
*

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 + Ž·C = [X’ 0，x ‘ 1，…，X’ RNK（A）-1，0，…，0] Ť

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