我如何才能获得与Matlab的A \ b(mldivide)运算符使用numpy/scipy返回的欠定线性系统相同的“特殊"解决方案? [英] How can I obtain the same 'special' solutions to underdetermined linear systems that Matlab's `A \ b` (mldivide) operator returns using numpy/scipy?
问题描述
I found a link where it is shown with an example that the Matlab mldivide operator (\) gives 'special' solutions when the system of linear equations has infinitely many solutions.
例如:
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给出了相同的结果.
The solution is 'special' in the sense that the number of non-zero entries in the solution c_mldivide is equal to rank(A) (in this case 2). I tried the same thing in numpy using numpy.linalg.lstsq, which gives an identical result to c_pinv.
有没有办法在Python中实现c_mldivide解决方案?
Is there a way to achieve the c_mldivide solution in Python?
在此处，还有另一个非常相似的问题，但我想对特殊"一词不够清楚.
另一个问题询问了内部mldivide运算符的工作原理，但可接受的答案似乎无法解决此问题.
There was another very similar question here, but I suppose the explanation of the word 'special' was not clear enough.
Another question asked about the internal workings of the mldivide operator, but the accepted answer doesn't seem to address this behavior.
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]))
使用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
推荐答案
Non-negative least squares (scipy.optimize.nnls) is not a general solution to this problem. A trivial case where it will fail is if all of the possible solutions contain negative coefficients:
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 的情况下，
In the case where A·x = b is underdetermined,
x1, res, rnk, s = np.linalg.lstsq(A, b)
将选择一个解决方案 x' ，该解决方案将 || x || L2 最小化em>符合 || A·x - b || L2 = 0 .碰巧这不是我们要寻找的特定解决方案，但是我们可以对其进行线性转换以得到我们想要的.为此，我们将首先计算< A ，它表示 A·x = b .我们可以使用等级显示QR分解:
will pick a solution x' that minimizes ||x||L2 subject to ||A·x - b||L2 = 0. This happens not to be the particular solution we are looking for, but we can linearly transform it to get what we want. In order to do that, we'll first compute the 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 :
Z is a vector (or, in case where n - rnk(A) > 1, a set of basis vectors spanning a subspace of A) such that 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 )设置为零
In other words, the column(s) of Z are vectors that are orthogonal to all of the rows in A. This means that for any solution x' to A·x = b, then x' = x + Z·c must also be a solution for any arbitrary scaling factor c. This means that by picking an appropriate value of c, we can set any n - rnk(A) of the coefficients in the solution to zero.
例如，假设我们要将最后一个系数的值设置为零:
For example, let's say we wanted to set the value of the last coefficient to zero:
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 有点复杂:
The more general case where n - rnk(A) ≤ 1 is a little bit more complicated:
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 :
We get x' and Z as before:
x1, res, rnk, s = np.linalg.lstsq(A, b)
Z = qr_null(A)
现在，为了最大化解向量中零值系数的数量，我们想找到一个向量 C
Now in order to maximise the number of zero-valued coefficients in the solution vector, we want to find a vector C such that
x' = x + Z·C = [x' 0 ，x' 1 ，...，x' rnk(A)-1 ，0，...，0] T
x' = x + Z·C = [x'0, x'1, ..., x'rnk(A)-1, 0, ..., 0]T
如果 x' 中最后一个 n-rnk( A )系数为零，则此施加
If the last n - rnk(A) coefficients in x' are to be zeros, this imposes that
Z {rnk(A)，...，n} ·C =- x {rnk(A)，... ，n}
Z{rnk(A),...,n}·C = -x{rnk(A),...,n}
因此，我们可以求解 C (正是因为我们知道Z[rnk:]必须是全等级):
We can thus solve for C (exactly, since we know that Z[rnk:] must be full-rank):
C = np.linalg.solve(Z[rnk:], -x1[rnk:])
并计算 x' :
and compute 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.]
将所有内容整合为一个函数:
To put it all together into a single function:
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)
这篇关于我如何才能获得与Matlab的A \ b(mldivide)运算符使用numpy/scipy返回的欠定线性系统相同的“特殊"解决方案?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
