背景
我在做一个项目,涉及求解大型未充分确定的方程组。
我目前的算法计算了一个表示给定系统的矩阵的SVD(numpy.linalg.svd),然后用它的结果计算出该矩阵的Moore-Penrose伪逆和右零空间。我使用nullspace来查找具有唯一解的所有变量,使用伪逆来查找其值。
但是,MPP(Moore-Penrose伪逆)非常密集,对于我的服务器来说有点太大,无法处理。
问题
我发现下面的paper详细描述了一个稀疏伪逆,它保持了MPP的大部分基本性质。这显然是我非常感兴趣的,但我只是没有数学背景来理解他是如何计算伪逆的。可以用SVD计算吗?如果没有,最好的办法是什么?
细节
这些是我认为可能相关的文章,但我还没有过时到足以理解
自旋v(A)=arg min | | B | | |服从BA=其中| | B | |表示B的入口l1范数
这通常是一个不可处理的问题,所以我们使用l1范数的标准线性松弛
sspinv(A)=ητ{[spinv(A)]},其中ητ(u)=u1 | u |≥τ
编辑
查找我的代码和有关实际实现的更多详细信息here

最佳答案

据我所知,关于稀疏伪逆,这篇论文是这么说的:
上面写着
我们的目标是最小化自旋v(A)中的非零数
这意味着您应该采用L0范数(参见David Donoho的定义here:非零条目数),这使得问题难以解决。
spinv(A) = argmin ||B||_0 subject to B.A = I
所以他们把这个问题转化为凸松弛问题,这样就可以用线性规划来解决。
这通常是一个无法处理的问题,所以我们使用标准
线性松弛与'1范数。
放松的问题是
spinv(A) = argmin ||B||_1 subject to B.A = I (6)
这有时被称为Basis pursuit并倾向于产生稀疏解(参见Boyd和Vandenberghe的凸优化,第6.2节最小范数问题)。
所以,解决这个轻松的问题。
线性程序(6)是可分离的,可以通过计算一个
一次B行
因此,您可以通过解决以下形式的一系列问题来获得解决方案。
spinv(A)_i = argmin ||B_i||_1 subject to B_i.A = I_i
其中_i表示矩阵的第i行。
请参见here以了解如何将此绝对值问题转换为线性程序。
在下面的代码中,我将问题稍微更改为spinv(A)_i = argmin ||B_i||_1 subject to A.B_i = I_i,其中_i是矩阵的第I列,因此问题变成spinv(A) = argmin ||B||_1 subject to A.B = I。老实说,我不知道这两者之间有没有区别。这里我用的是scipylinprog单纯形法。我不知道simplex的内部是否使用SVD。

import numpy as np
from scipy import optimize

# argmin ||B_i||_1 stubect to A.B_i = I_i, where _i is the ith column
# let B_i = u_i - v_i where u_i >= 0 and v_i >= 0
# then ||B_i||_1 = [1' 1'][u_i;v_i] which is the objective function
# and A.B_i = I_i becomes
# A.[u_i - v_i] = I_i
# [A -A][u_i;v_i] = I_i which is the equality constraint
# and [u_i;v_i] >= 0 the bounds
# here A is n x m (as opposed to m x n in paper)

A = np.random.randn(4, 6)
n, m = A.shape
I = np.eye(n)

Aeq = np.hstack((A, -A))
# objective
c = np.ones((2*m))
# spinv
B = np.zeros((m, n))

for i in range(n):
    beq = I[:, i]
    result = optimize.linprog(c, A_eq=Aeq, b_eq=beq)
    x = result.x[0:m]-result.x[m:2*m]
    B[:, i] = x

print('spinv(A) = \n' + str(B))
print('pinv(A) = \n' + str(np.linalg.pinv(A)))
print('A.B = \n' + str(np.dot(A, B)))

这是一个输出。spinv(A)pinv(A)稀疏。
spinv(A) =
[[ 0.         -0.33361925  0.          0.        ]
 [ 0.04987467  0.          0.12741509  0.02897778]
 [ 0.          0.         -0.52306324  0.        ]
 [ 0.43848257  0.12114828  0.15678815 -0.19302049]
 [-0.16814546  0.02911103 -0.41089271  0.50785258]
 [-0.05696924  0.13391736  0.         -0.43858428]]
pinv(A) =
[[ 0.05626402 -0.1478497   0.19953692 -0.19719524]
 [ 0.04007696 -0.07330993  0.19903311  0.14704798]
 [ 0.01177361 -0.05761487 -0.23074996  0.15597663]
 [ 0.44471989  0.13849828  0.18733242 -0.20824972]
 [-0.1273604   0.15615595 -0.24647117  0.38047901]
 [-0.04638221  0.09879972  0.21951122 -0.33244635]]
A.B =
[[ 1.00000000e+00 -1.82225048e-17  6.73349443e-18 -2.39383542e-17]
 [-5.20584593e-18  1.00000000e+00 -3.70118759e-16 -1.62063433e-15]
 [-8.83342417e-18 -5.80049814e-16  1.00000000e+00  3.56175852e-15]
 [ 2.31629738e-17 -1.13459832e-15 -2.28503999e-16  1.00000000e+00]]

为了进一步稀疏矩阵,我们可以逐级应用
阈值化,从而牺牲反转属性并计算
近似稀疏伪逆
如果不想在稀疏pinv中保留小条目,可以这样删除它们:
Bt = B.copy()
Bt[np.abs(Bt) < 0.1] = 0
print('sspinv_0.1(A) = \n' + str(Bt))
print('A.Bt = \n' + str(np.dot(A, Bt)))

得到
sspinv_0.1(A) =
[[ 0.         -0.33361925  0.          0.        ]
 [ 0.          0.          0.12741509  0.        ]
 [ 0.          0.         -0.52306324  0.        ]
 [ 0.43848257  0.12114828  0.15678815 -0.19302049]
 [-0.16814546  0.         -0.41089271  0.50785258]
 [ 0.          0.13391736  0.         -0.43858428]]
A.Bt =
[[ 9.22717491e-01  1.17555372e-02  6.73349443e-18 -1.10993934e-03]
 [ 1.24361576e-01  9.41538212e-01 -3.70118759e-16  1.15028494e-02]
 [-8.76662313e-02 -1.36349311e-02  1.00000000e+00 -7.48302663e-02]
 [-1.54387852e-01 -3.27969169e-02 -2.28503999e-16  9.39161039e-01]]

希望我能回答你的问题,并提供足够的参考资料,如果你想进一步的细节。如果你有任何问题,请告诉我。我不是专家,所以如果你对我的说法有任何疑问,你可以问数学堆栈交换专家(当然没有任何代码),请告诉我。
这是个有趣的问题。它使我能够重新学习我的线性代数和我所知道的少量优化,所以谢谢:)

关于python - 使用SVD计算Spinv,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/51273038/

10-09 18:16