非负最小二乘(^{})是这个问题的一般解。如果所有可能的解决方案都包含负系数,则会失败: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的情况下
^{pr2}$
将选择一个最小化L2的解决方案。这恰好不是我们要寻找的特定解决方案,但我们可以线性变换它以得到我们想要的。为了做到这一点,我们首先计算A的right null space,它表征了A·x=b的所有可能解的空间。我们可以使用rank-revealing QR decomposition:
^{3}$
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'0, x'1, ..., x'rnk(A)-1, 0, ..., 0]T
如果n-rnk(A)中的最后一个n-rnk(A)系数为零,则这就要求Z{rnk(A),...,n}·C = -x{rnk(A),...,n}
因此,我们必须清楚地知道(因为我们可以精确地求解)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)