计算科学 - 计算大型稀疏矩阵的对数行列式 - 吾爱随笔录

计算大型稀疏矩阵的对数行列式

计算科学矩阵稀疏矩阵计算物理学

2021-12-10 10:42:36

我需要计算其中是大型稀疏矩阵，它们是实数、对称和半正定矩阵。我希望有到个。物理上，每个是一个大型系统的刚度矩阵，如果有帮助，我知道如何将其写为（其中不是方阵） . $\log(\det (\mathbf M_i))$ $\mathbf M_i$ $10$ $100$ $\mathbf M_i$ $\mathbf M_i=\mathbf B_i^T\mathbf B_i$ $\mathbf B_i$

我不是计算科学家，所以我更喜欢使用书面包而不是实现我自己的算法。我目前在 MATLAB 和/或 Mathematica 工作，但我也精通 C++，如果有帮助，我愿意在任何其他环境中工作。

到目前为止，我尝试的是计算特征值（使用 MATLAB 的eig）并总结它们的日志，但我怀疑还有其他更有效的方法可以做到这一点。例如，我偶然发现了这篇论文，但我不确定它的算法是否适合我的情况，也找不到实现的代码。

另一个复杂因素是它可能会发生 $\mathbf M_i$ 不可逆，即 $\det\mathbf M_i=0$ ，在这种情况下，我对所有非零特征值的乘积感兴趣（在这些情况下，我确切地知道如何有许多零特征值）。

任何帮助将非常感激。

3个回答

就像 Brian Borchers 在他的评论中所说，行列式可以从 Cholesky 因子的对角线条目中确定。

理论：

考虑 Cholesky 分解其中

R^{T} R = S^{T} M S,

$R^TR = S^TMS,$

$R$ 是一个上三角矩阵，并且
$S$ 是一个置换矩阵。

然后

\begin{aligned} det (M) & = det (S R^{T} R S^{T}) \\ = det (S) det (R^{T}) det (R) det (S^{T}) \\ = det (R)^{2} . \end{aligned}

$\begin{align} \text{det}(M) &= \text{det}(SR^TRS^T) \\ &= \text{det}(S) \text{det}(R^T) \text{det}(R) \text{det}(S^T) \\ &= \text{det}(R)^2. \end{align}$

由于是三角形的，它的特征值是对角线条目，因此 $R$

\begin{aligned} log det (M) & = 2 log det (R) \\ = 2 \sum_{i} log (R_{i i}) . \end{aligned}

$\begin{align} \text{log}~\text{det}(M) &= 2 \text{log}~\text{det}(R) \\ &= 2 \sum_i \text{log}(R_{ii}). \end{align}$

代码：

执行此操作的 Matlab 代码是，例如，

[R,~,S] = chol(M); 
log_det_chol = sum(2*log(diag(R)));

在 python 中，可以使用 scikits.sparse.cholmod 包从稀疏 Cholesky 分解计算对数行列式。从链接页面解释，执行此操作的代码是：

from scikits.sparse.cholmod import cholesky
R = cholesky(M)
log_det_chol = R.logdet()

如果您只想考虑非零特征值，您可以对的对角线的非零条目求和。 $R$

测试用例：

我在使用 Dirichlet BC 的二维有限差分拉普拉斯算子上针对 Matlab 内置的 det() 函数进行了测试，得到了高达 400×400 左右的矩阵的精确一致性，之后使用 matlab 的 det() 得到 Inf，而cholesky 方法在我的笔记本电脑上达到百万倍的规模没有问题。下面是代码，供参考：

n = 20; % M is n^2-by-n^2. Matlab's det() returns Inf for n=25+
M = gallery('poisson',n); 
[R,~,S] = chol(M); 
log_det_chol = sum(2*log(diag(R))), 
log_det = log(det(M))

内存注意事项：

由于此方法基于 Cholesky 分解，因此它继承了速度和稳定性。

对于大问题，使用这种方法的瓶颈可能是内存。特别是，使用具有良好旋转（Matlab 很好）的 Cholesky 分解，您可以预期需要内存，其中是划分自由分成不同的部分。 $O(s^2)$ $s$

例如，在具有 n×n 网格的 2D 问题中，将域分成两部分的 1D 分离器具有的自由度。因此所需的内存量为，与存储右侧的数量级相同。 $s = O(n)$ $O(n^2)$

一个 3D n×n×n 问题有一个 2D 分离器前端，其自由度为内存。 $s = O(n^2)$ $O(n^4)$

在算法上，一个快速、稳定的行列式方法是首先使用 Householder 反射转换为对称三对角线形式（这是 eig 在对称矩阵上使用的第一步；这个子操作在 matlab 中可能是可见的）。

对称三对角矩阵的行列式可以通过以相当简单的方式沿对角线工作来找到。但是，它需要在每个步骤中进行乘法和加法，因此如果最终值（或中间值）太大或太小而无法以对数形式表示，则需要防止该过程发生上溢或下溢。

对于特征值为零的情况 - 我不确定这里会发生什么，如果你在对角线中出现“间隙”（零对角线元素和 4 个相邻的非对角线元素也为零），那么你可以处理子矩阵分别在每一侧。但它可能不会那么干净。对于这种情况，您可能只想使用“没有特征向量的 eig”操作。

我在玩弄 K_n（三对角线 -1 2 -1）的稀疏表示，理论上它的行列式是 n+1。我遇到了这个旧帖子。一旦超过 n=1e6，SuperLU 和 CHOLMOD 似乎都没有获得正确行列式的精度。

from scipy.sparse import diags
from scipy.sparse.linalg import splu
from sksparse.cholmod import cholesky
from math import exp

n=int(5e6)
K = diags([-1.],-1,shape=(n,n)) + diags([2.],shape=(n,n)) + diags([-1.],1,shape=(n,n))
lu = splu(K.tocsc())
diagL = lu.L.diagonal()
diagU = lu.U.diagonal()
det=diagL.prod()*diagU.prod()
print(det)

factor = cholesky(K.tocsc())
ld = factor.logdet()
print(exp(ld))

4999993.625461911

4999993.625461119

其它你可能感兴趣的问题

上一篇能不能把这个简单的二次优化问题变成一个简单的特征值问题？下一篇贝塞尔函数计算的最快开源实现是什么？