c++ - 使用 Eigen::Cholmod Supernodal LLT 时如何提取 matrixyL() 和 matrixU()？

Question

我正在尝试使用 Eigen::CholmodSupernodalLLT然而，对于 Cholesky 分解，我似乎无法得到 matrixL()和 matrixU() .如何提取matrixL()和 matrixU()来自 Eigen::CholmodSupernodalLLT供将来使用？

Answer 1

整合其他人所说的部分答案。

考虑 Y ~ MultivariateNormal(0, A)。人们可能想要 (1) 评估 (log-) 似然性(多元正态密度)，(2)来自这种密度的样本。

对于(1)，需要求解 Ax = b 其中 A 是对称正定的，并计算其对数行列式。 (2) 要求 L 使得 A = L * L.transpose() 因为 Y ~ MultivariateNormal(0, A) 可以找到为 Y = L u 其中 u ~ MultivariateNormal(0, I)。

Cholesky LLT 或 LDLT 分解很有用，因为 chol(A)可用于这两个目的。求解Ax=b给定分解很容易，并且(log)行列式可以很容易地从 D 的(log-)分量的(和)乘积导出或 L 的对角线.根据定义，L 可用于采样。

所以，在 Eigen可以使用:

Eigen::SimplicialLDLT solver(A) (或 Eigen::SimplicialLLT)，当 solver.solve(b)并使用 solver.vectorD().diag() 计算行列式.很有用，因为如果 A是协方差矩阵，则 solver可用于似然评估，以及 matrixL()用于采样。

Eigen::CholmodDecomposition不允许访问 matrixL()或 vectorD()但暴露.logDeterminant()实现(1)目标但不实现(2)。

Eigen::PardisoLDLT不允许访问 matrixL()或 vectorD()并且没有公开获得行列式的方法。

在某些应用程序中，步骤 (2) - 采样 - 可以在稍后阶段完成，因此 Eigen::CholmodDecomposition足够。至少在我的配置中， Eigen::CholmodDecomposition比 Eigen::SimplicialLDLT 快 2 到 5 倍(我猜是因为在引擎盖下进行了排列以促进并行化)

示例:在贝叶斯空间高斯过程回归中，空间随机效应可以被积分出来，不需要采样。因此 MCMC 可以通过 Eigen::CholmodDecomposition 迅速进行以实现未知参数的收敛。然后可以使用 Eigen::SimplicialLDLT 并行恢复空间随机效应。 .通常这只是计算的一小部分，但具有 matrixL()直接来自 CholmodDecomposition将它们简化一点。