我怎样才能获得与Matlab的'A\b`(mldivide)算子使用numpy/scipy返回的欠定线性系统相同的"特殊"解？

Question

我找到了一个链接,其中显示了一个示例,当线性方程组具有无限多个解时,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中尝试了相同的东西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' ₀,x ' ₁,...,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)