检查存储在 scipy 稀疏矩阵中的对称矩阵的半定的更好方法？

Question

我有一个非常大的对称矩阵要在 RAM 中存储和操作（约40,000 * 40,000），所以我使用 scispy.sparse 格式来存储它的一半，下面是我的代码

import numpy as np
from scipy.sparse import coo_matrix
def to_sparse(W):
    tmp = np.tril(W)
    del W
    return coo_matrix(tmp)

然后我想根据公式L = D - W，计算W的拉普拉斯矩阵（对称矩阵），（其中D是对角矩阵，其对角线是W 的列或行的总和，根据L 的定义，我需要检查它是否是正半正定（PSD）。 但是好像计算稀疏格式存储的'L'的特征值不等于原来的问题，所以我必须将L转换为密集格式，然后检查？ 我的问题：有没有更好的方法可以在不将 L 转换回密集矩阵的情况下实现这一目标？

下面是我生成 L 并检查 PSD 的代码：

from scipy.sparse import diags
from numpy import linalg
def gen_D(W):
    # W is sparse matrix
    tmp = W.sum(axis = 0)
    tmp2 = W.sum(axis = 1).T - W.diagonal()
    return np.squeeze(np.asarray(tmp + tmp2))
def gen_laplacian(W):
    D = gen_D(W)
    return diags(D, 0).tocoo() - W
def check_PSD(L):
    origin_L = L.todense() + L.T.todense() - np.diag(L.diagonal())
    E, V = linalg.eigh(origin_L)
    return np.all(E > -1e-10)

很抱歉我在上传之前没有检查代码示例，现在我已经修复了错误，下面是一个示例：

W = np.random.random_integers(1, 100, size = (100, 100))
W = .5 * (W + W.T)
W = to_sparse(W)
L = gen_laplacian(W)

Answer 1

原始数组是稀疏的（很多零），还是只是tril 的乘积？如果是后者，您可能不会通过使用稀疏代码来节省空间或时间。例如

def gen1(W):
    tmp = np.tril(W)
    d = tmp.sum(0)+tmp.sum(1)-tmp.diagonal()
    return np.diag(d) - tmp

在这个 100x100 测试中比 gen_laplacian(to_sparse(W)) 快 8 倍。

def check(L):
    oL = L+L.T-np.diag(L.diagonal())
    E, V = linalg.eigh(oL)
    return np.all(E > -1e-10)

即使 check(L.A) 也快 2 倍。

但是你需要构造orginal_L吗？看起来linalg.eigh(L.A) 产生了同样的东西，因为the calculation is done with the lower triangular part。并不是说构建orginal_L 需要很多时间。

但测试问题的症结在于，sparse.linalg 之类的 eigsh 函数是否会产生与 eigh 等价的东西。

我对那些 eigh 函数了解不多，无法提供帮助 - 没有深入研究文档。