如何使用 numpy/scipy 获得与 Matlab 的 `A \ b` (mldivide) 运算符 returns 相同的 'special' 欠定线性系统解？
How can I obtain the same 'special' solutions to underdetermined linear systems that Matlab's `A \ b` (mldivide) operator returns using numpy/scipy?
我找到了一个link，其中显示了一个例子，当线性方程组有无限多个解决方案。
例如：
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
解决方案是 'special'，因为解决方案 c_mldivide 中的非零条目数等于 rank(A)（在本例中为 2）。我在 numpy 中使用 numpy.linalg.lstsq 尝试了同样的事情，它给出了与 c_pinv.
相同的结果
有没有办法在Python中实现c_mldivide的解决方案？
还有一个非常相似的问题here，但我想'special'这个词的解释不够清楚。
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
np.array([[8],[18]]).shape
是
(2,1)
但是你想要
(2,)
#!/usr/bin/env python3
import numpy as np
from scipy.optimize import nnls
test_A = np.array([[1,2,0],[0,4,3]])
try:
test_b = np.array([[8],[18]]) # wrong
print(nnls(test_A,test_b))
except Exception as e:
print(str(e))
test_b = np.array([8,18]) # sic!
print(nnls(test_A,test_b))
输出：
expected vector
(array([ 0. , 4. , 0.66666667]), 0.0)
非负最小二乘法 (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 of A, which characterizes the space of all possible solutions to A·x = b. We can get this using a 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的列是与[=76中的所有行正交的向量=]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'0, x'1, ..., x'rnk(A)-1, 0, ..., 0]T
如果n - rnk(A)系数在x' 为零，这意味着
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)
我找到了一个link，其中显示了一个例子，当线性方程组有无限多个解决方案。
例如：
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
解决方案是 'special'，因为解决方案 c_mldivide 中的非零条目数等于 rank(A)（在本例中为 2）。我在 numpy 中使用 numpy.linalg.lstsq 尝试了同样的事情，它给出了与 c_pinv.
有没有办法在Python中实现c_mldivide的解决方案？
还有一个非常相似的问题here，但我想'special'这个词的解释不够清楚。
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
np.array([[8],[18]]).shape
是
(2,1)
但是你想要
(2,)
#!/usr/bin/env python3
import numpy as np
from scipy.optimize import nnls
test_A = np.array([[1,2,0],[0,4,3]])
try:
test_b = np.array([[8],[18]]) # wrong
print(nnls(test_A,test_b))
except Exception as e:
print(str(e))
test_b = np.array([8,18]) # sic!
print(nnls(test_A,test_b))
输出：
expected vector
(array([ 0. , 4. , 0.66666667]), 0.0)
非负最小二乘法 (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 of A, which characterizes the space of all possible solutions to A·x = b. We can get this using a 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的列是与[=76中的所有行正交的向量=]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'0, x'1, ..., x'rnk(A)-1, 0, ..., 0]T
如果n - rnk(A)系数在x' 为零，这意味着
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)