使用SVD计算spinv

Question

使用SVD计算spinv

pythonsvdmatrix-inverse

4

背景

我正在处理一个涉及解决大型欠定方程组的项目。

我的当前算法计算表示给定系统的矩阵的SVD（numpy.linalg.svd），然后使用其结果计算Moore-Penrose伪逆和矩阵的右零空间。我使用零空间来找到所有具有唯一解的变量，使用伪逆来确定它的值。

然而，MPP（Moore Penrose伪逆）相当密集，对于我的服务器来说有点太大了。

问题

我发现了以下论文，详细介绍了一种保持MPP大部分基本属性的更稀疏的伪逆。这显然非常吸引我，但我简单地没有数学背景来理解他如何计算伪逆。是否可以使用SVD计算它？如果不能，最好的方法是什么？

细节

这些是我认为可能相关但我不够老练以理解的论文行

spinv（A）= arg min ||B|| subject to BA = In其中||B||表示B的逐项l1范数
这通常是一个非可计算问题，因此我们使用具有l1范数的标准线性松弛
sspinv（A）= ητ {[spinv（A）]}，其中ητ（u）= u1 | u |≥τ

编辑

在这里找到我的代码和更多实现细节

- Rushabh Mehta

1

这个问题太宽泛了。你有具体的问题和相关的代码吗？也许这更适合数学堆栈交换？ - erip

我可以给一个包含所有代码的问题打标签。 - Rushabh Mehta

1

这个问题因为与计算机科学相关而在math.stackexchange上被关闭，所以我被发送到这里。 - Rushabh Mehta

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- dhanushka · Accepted Answer

根据论文所说，稀疏伪逆的目标是最小化 spinv(A) 中非零元素的数量。这意味着应该采用 L0 范数（请参见 David Donoho 的定义 here：非零条目的数量），但这使问题难以处理。因此，他们转向对该问题进行凸松弛，以便通过线性规划来解决。

一般情况下，这是一个难以处理的问题，因此我们使用标准的线性松弛和“1范数”。经过松弛后的问题如下：spinv(A) = argmin ||B||_1 subject to B.A = I (6)。

有时被称为基 Pursuit，倾向于产生稀疏解决方案（参见Boyd和Vandenberghe的凸优化第6.2节最小规范问题）。

因此，解决这个放松的问题。

线性规划（6）是可分离的，并且可以通过逐行计算B来解决

因此，您可以解决一系列以下形式的问题以获得解决方案。 spinv(A)_i = argmin ||B_i||_1 subject to B_i.A = I_i 其中_i表示矩阵的第i行。

请参见此处，了解如何将此绝对值问题转换为线性规划。

在下面的代码中，我稍微修改了问题，使得 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。老实说，我不知道这两者之间是否有区别。这里我使用scipy的linprog单纯形法。我不知道单纯形法的内部是否使用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]]

为了进一步稀疏矩阵，我们可以应用逐元素硬阈值处理，从而牺牲反演属性并计算近似稀疏伪逆。

如果您不想在稀疏伪逆中保留小条目，可以像这样删除它们。

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]]

希望我回答了你的问题并提供了足够的参考资料，如果你需要进一步了解，请告诉我。我不是专家，如果你对我的观点有任何疑问，你可以在数学堆栈交换中询问专家（当然不要附上任何代码），并请告知我。

这是一个有趣的问题。它让我复习了线性代数和我所知道的优化知识的一些内容，所以谢谢你 :)