使用 svd 求解欠定的 scipy.sparse 矩阵

Question

问题

我有一组方程，其中变量用小写变量表示，常量用大写变量表示

A = a + b  
B = c + d  
C = a + b + c + d + e  

我在 pandas DataFrame 中提供了有关这些方程式结构的信息，其中包含两列：Constants 和 Variables

例如

df = pd.DataFrame([['A','a'],['A','b'],['B','c'],['B','d'],['C','a'],['C','b'], 
['C','c'],['C','d'],['C','e']],columns=['Constants','Variables'])

然后我使用 NetworkX 将其转换为稀疏 CSC 矩阵

table = nx.bipartite.biadjacency_matrix(nx.from_pandas_dataframe(df,'Constants','Variables')  
,df.Constants.unique(),df.Variables.unique(),format='csc')

当转换为密集矩阵时，表格如下所示

矩阵([[1, 1, 0, 0, 0],[0, 0, 1, 1, 0],[1, 1, 1, 1, 1]], dtype=int64)

我想从这里找出哪些变量是可解的（在这个例子中，只有 e 是可解的），对于每个可解变量，它的值依赖于哪些常量（在这种情况下，由于e = C-B-A，所以依赖于A， B 和 C)

解决方案的尝试

我首先尝试使用 rref 来求解可解变量。我使用了符号库 sympy 和函数 sympy.Matrix.rref，这正是我想要的，因为任何可解变量都有自己的行，几乎全是零和 1 个，我可以逐行检查。

但是，此解决方案并不稳定。首先，它非常慢，并且没有利用我的数据集可能非常稀疏的事实。此外， rref 在浮点方面做得不太好。所以我决定转向另一种受Removing unsolvable equations from an underdetermined system 启发的方法，它建议使用 svd

方便的是，scipy.sparse库中有一个svd函数，即scipy.sparse.linalg.svds。但是，由于我缺乏线性代数背景，我不明白在我的桌子上运行这个函数输出的结果，或者如何使用这些结果来得到我想要的。

问题的更多细节

1. 我的问题中每个变量的系数都是1。这就是前面显示的两列pandas DataFrame中数据的表达方式
2. 我的实际示例中的绝大多数变量都无法求解。目标是找到少数可解决的问题
3. 如果它符合此问题的限制条件，我非常愿意尝试另一种方法。

这是我第一次发布问题，所以如果这不完全符合准则，我深表歉意。请留下建设性的批评，但要温柔！

Answer 1

以 ewcz 的答案为基础，可以使用 numpy.linalg.svd 计算零空间和伪逆。请参阅以下链接：

pseudo-inverse

nullspace

Answer 2

你正在解决的系统有形式

[ 1 1 0 0 0 ] [a]   [A]
[ 0 0 1 1 0 ] [b] = [B]
[ 1 1 1 1 1 ] [c]   [C]
              [d]
              [e]

即五个变量a, b, c, d, e 的三个方程。正如您问题中链接的答案所提到的，可以使用pseudoinverse 来解决这种不确定的系统，Numpy 直接根据pinv 函数提供了该系统。

由于M 具有线性独立的行，因此在这种情况下，psudoinverse 具有M.pinv(M) = I 的属性，其中I 表示单位矩阵（在这种情况下为 3x3）。因此，形式上，我们可以将解决方案写成：

v = pinv(M) . b

其中v 是 5 分量解向量，b 表示右侧 3 分量向量 [A, B, C]。然而，这个解决方案并不是唯一的，因为可以从所谓的内核或矩阵M 的null space 中添加一个向量（即w 的向量M.w=0），它仍然是解决方案：

M.(v + w) = M.v + M.w = b + 0 = b

因此，唯一有唯一解的变量是那些来自M 的零空间的所有可能向量的对应分量为零的变量。换句话说，如果你将零空间的基组装成一个矩阵（每列一个基向量），那么“可解变量”将对应于该矩阵的零行（列的任何线性组合的相应分量将也为零）。

让我们将此应用于您的特定示例：

import numpy as np
from numpy.linalg import pinv

M = [
    [1, 1, 0, 0, 0],
    [0, 0, 1, 1, 0],
    [1, 1, 1, 1, 1]
]

print(pinv(M))

[[ 5.00000000e-01 -2.01966890e-16  1.54302378e-16]
 [ 5.00000000e-01  1.48779676e-16 -2.10806254e-16]
 [-8.76351626e-17  5.00000000e-01  8.66819360e-17]
 [-2.60659800e-17  5.00000000e-01  3.43000417e-17]
 [-1.00000000e+00 -1.00000000e+00  1.00000000e+00]]

从这个伪逆，我们看到变量e（最后一行）确实可以表示为- A - B + C。但是，它也“预测”了a=A/2 和b=A/2。为了消除这些非唯一的解决方案（例如，a=A 和 b=0 也同样有效），让我们借用 SciPy Cookbook 的函数来计算空空间：

print(nullspace(M))

[[ 5.00000000e-01 -5.00000000e-01]
 [-5.00000000e-01  5.00000000e-01]
 [-5.00000000e-01 -5.00000000e-01]
 [ 5.00000000e-01  5.00000000e-01]
 [-1.77302319e-16  2.22044605e-16]]

这个函数已经返回了组装成矩阵的零空间的基础（每列一个向量），我们看到，在合理的精度内，唯一的零行确实只是对应于变量e 的最后一行.

编辑：

对于方程组

A = a + b, B = b + c, C = a + c

对应的矩阵M是

[ 1 1 0 ]
[ 0 1 1 ]
[ 1 0 1 ]

这里我们看到矩阵实际上是正方形的，并且是可逆的（行列式是2）。因此，伪逆与“正常”逆相吻合：

[[ 0.5 -0.5  0.5]
 [ 0.5  0.5 -0.5]
 [-0.5  0.5  0.5]]

对应于解决方案a = (A - B + C)/2, ...。由于M 是可逆的，它的内核/空空间是空的，这就是为什么cookbook 函数只返回[]。为了看到这一点，让我们使用内核的定义——它由所有非零向量x 组成，这样M.x = 0。然而，由于M^{-1} 存在，x 被给出为x = M^{-1} . 0 = 0，这是一个矛盾。形式上，这意味着找到的解决方案是唯一的（或者所有变量都是“可解决的”）。