用 SVD 计算 spinv
Calculating spinv with SVD
背景
我正在从事一个涉及求解大型欠定方程组的项目。
我当前的算法计算表示给定系统的矩阵的 SVD (numpy.linalg.svd)，然后使用其结果计算 Moore-Penrose 伪逆和矩阵的右零空间。我使用零空间来查找具有唯一解的所有变量，并使用伪逆来找出它的值。
但是，MPP（Moore Penrose 伪逆）非常密集，对于我的服务器来说有点太大，无法处理。
问题
我发现了以下 paper，它详细描述了一个更稀疏的伪逆，它保持了 MPP 的大部分基本属性。这显然对我很感兴趣，但我只是没有数学背景来理解他是如何计算伪逆的。可以用SVD计算吗？如果没有，最好的方法是什么？
详情
这些是我认为可能相关的论文中的几行，但我还不够陈旧，无法理解
spinv(A) = arg min ||B||受制于 BA = 其中 ||B||表示 B
这通常是一个不易处理的问题，因此我们使用标准线性松弛和 l1 范数
sspinv(A) = ητ {[spinv(A)]}，其中 ητ (u) = u1|u|≥τ
编辑
找到我的代码和实际实现的更多细节here
据我了解，论文中是这样描述稀疏伪逆的：
它说
We aim at minimizing the number of non-zeros in spinv(A)
这意味着您应该采用 L0 范数（参见 David Donoho 的定义 here：非零项的数量），这使得问题变得棘手。
spinv(A) = argmin ||B||_0 subject to B.A = I
所以他们求助于这个问题的凸松弛，这样就可以用线性规划来解决了。
This is in general a non-tractable problem, so we use the standard
linear relaxation with the `1 norm.
那么轻松的问题就是
spinv(A) = argmin ||B||_1 subject to B.A = I (6)
这有时称为 Basis pursuit 并倾向于产生稀疏解（请参阅 Boyd 和 Vandenberghe 的 凸优化，第 6.2 节最小范数问题).
所以，解决这个轻松的问题。
The linear program (6) is separable and can be solved by computing one
row of B at a time
所以，你可以解决下面表格的一系列问题来获得解决方案。
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。老实说，我不知道两者之间是否有区别。这里我使用了 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]]
To further sparsify the matrix, we can apply entrywise hard
thresholding, thus sacrificing the inverting property and computing an
approximate sparse pseudoinverse
如果你不想在 sparse-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]]
希望我回答了您的问题并提供了足够的参考资料（如果您需要更多详细信息）。如果您有任何问题，请告诉我。我不是专家，所以如果您对我的说法有任何疑问，您可以随时询问数学堆栈交换方面的专家（当然没有任何代码），请让我知道。
这是一个有趣的问题。它使我复习了我的线性代数和我知道的很少的优化，所以谢谢你:)
背景
我正在从事一个涉及求解大型欠定方程组的项目。
我当前的算法计算表示给定系统的矩阵的 SVD (numpy.linalg.svd)，然后使用其结果计算 Moore-Penrose 伪逆和矩阵的右零空间。我使用零空间来查找具有唯一解的所有变量，并使用伪逆来找出它的值。
但是，MPP（Moore Penrose 伪逆）非常密集，对于我的服务器来说有点太大，无法处理。
问题
我发现了以下 paper，它详细描述了一个更稀疏的伪逆，它保持了 MPP 的大部分基本属性。这显然对我很感兴趣，但我只是没有数学背景来理解他是如何计算伪逆的。可以用SVD计算吗？如果没有，最好的方法是什么？
详情
这些是我认为可能相关的论文中的几行，但我还不够陈旧，无法理解
spinv(A) = arg min ||B||受制于 BA = 其中 ||B||表示 B
的 entrywise l1 范数
这通常是一个不易处理的问题，因此我们使用标准线性松弛和 l1 范数
sspinv(A) = ητ {[spinv(A)]}，其中 ητ (u) = u1|u|≥τ
编辑
找到我的代码和实际实现的更多细节here
据我了解，论文中是这样描述稀疏伪逆的：
它说
We aim at minimizing the number of non-zeros in spinv(A)
这意味着您应该采用 L0 范数（参见 David Donoho 的定义 here：非零项的数量），这使得问题变得棘手。
spinv(A) = argmin ||B||_0 subject to B.A = I
所以他们求助于这个问题的凸松弛，这样就可以用线性规划来解决了。
This is in general a non-tractable problem, so we use the standard linear relaxation with the `1 norm.
那么轻松的问题就是
spinv(A) = argmin ||B||_1 subject to B.A = I (6)
这有时称为 Basis pursuit 并倾向于产生稀疏解（请参阅 Boyd 和 Vandenberghe 的 凸优化，第 6.2 节最小范数问题).
所以，解决这个轻松的问题。
The linear program (6) is separable and can be solved by computing one row of B at a time
所以，你可以解决下面表格的一系列问题来获得解决方案。
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。老实说，我不知道两者之间是否有区别。这里我使用了 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]]
To further sparsify the matrix, we can apply entrywise hard thresholding, thus sacrificing the inverting property and computing an approximate sparse pseudoinverse
如果你不想在 sparse-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]]
希望我回答了您的问题并提供了足够的参考资料（如果您需要更多详细信息）。如果您有任何问题，请告诉我。我不是专家，所以如果您对我的说法有任何疑问，您可以随时询问数学堆栈交换方面的专家（当然没有任何代码），请让我知道。
这是一个有趣的问题。它使我复习了我的线性代数和我知道的很少的优化，所以谢谢你:)