Eigen - 检查矩阵是否为正（半）定答案

【问题标题】：Eigen - Check if matrix is Positive (Semi-)DefiniteEigen - 检查矩阵是否为正（半）定
【发布时间】：2016-02-05 14:58:15
【问题描述】：

我正在实现一个谱聚类算法，我必须确保一个矩阵（拉普拉斯算子）是半正定的。

检查矩阵是否为正定 (PD) 就足够了，因为可以在特征值中看到“半”部分。矩阵非常大（nxn，其中 n 大约为数千），因此特征分析很昂贵。

在 Eigen 中是否有任何检查会在运行时给出 bool 结果？

如果矩阵不是 PD，Matlab 可以使用 chol() 方法通过抛出异常来给出结果。遵循这个想法，Eigen 返回一个结果而不抱怨LLL.llt().matrixL()，尽管我期待一些警告/错误。 Eigen 也有方法isPositive，但由于bug，它不能用于具有旧 Eigen 版本的系统。

【问题讨论】：

你不能先检查它是否是厄米特，然后再看特征值吗？检查封闭性很简单。
你对厄米特部分是正确的，但理想情况下，我想避免多次计算一个巨大矩阵的特征值，因为这是我想要的输出，所以我希望尽可能只发生一次。
也许您可以尝试从 Eigen 进行 Cholesky 分解，如果矩阵为负，则返回 NumericalIssue，请参阅 eigen.tuxfamily.org/dox/classEigen_1_1LLT.html
可能是版本问题，因为：error: 'class Eigen::LDLT has no member named 'info'
哦，实际上它适用于 LLT，虽然它不适用于 LLT，感谢您的指点！如果你愿意，你可以写一个答案让我接受，否则我稍后会发布一个 sn-p 作为答案。

标签： c++ eigen eigen3

【解决方案1】：

您可以使用 Cholesky 分解 (LLT)，如果矩阵为负，则返回 Eigen::NumericalIssue，请参阅 documentation。

下面的例子：

#include <Eigen/Dense>

#include <iostream>
#include <stdexcept>

int main()
{
    Eigen::MatrixXd A(2, 2);
    A << 1, 0 , 0, -1; // non semi-positive definitie matrix
    std::cout << "The matrix A is" << std::endl << A << std::endl;
    Eigen::LLT<Eigen::MatrixXd> lltOfA(A); // compute the Cholesky decomposition of A
    if(lltOfA.info() == Eigen::NumericalIssue)
    {
        throw std::runtime_error("Possibly non semi-positive definitie matrix!");
    }    
}

【讨论】：

令人困惑，因为该文档还说“请记住，Cholesky 分解不能揭示秩。这种 LLT 分解仅在正定矩阵上稳定，使用 LDLT 代替半定情况。”我想这取决于你如何定义“排名揭示”。

【解决方案2】：

除了@vsoftco 的回答，我们还要检查矩阵对称性，因为 PD/PSD 的定义需要对称矩阵。

Eigen::LLT<Eigen::MatrixXd> A_llt(A);
if (!A.isApprox(A.transpose()) || A_llt.info() == Eigen::NumericalIssue) {
    throw std::runtime_error("Possibly non semi-positive definitie matrix!");
}

这项检查很重要，例如一些 Eigen 求解器（如 LTDT）需要 PSD（或 NSD）矩阵输入。事实上，存在 非对称因此非 PSD 矩阵 A 通过 A_llt.info() != Eigen::NumericalIssue 测试。考虑以下示例（数字取自 Jiuzhang Suanshu，第 8 章，问题 1）：

Eigen::Matrix3d A;
Eigen::Vector3d b;
Eigen::Vector3d x;

// A is full rank and all its eigen values >= 0
// However A is not symmetric, thus not PSD
A << 3, 2, 1, 
     2, 3, 1, 
     1, 2, 3;
b << 39, 34, 26;

// This alone doesn't check matrix symmetry, so can't guarantee PSD
Eigen::LLT<Eigen::Matrix3d> A_llt(A);
std::cout << (A_llt.info() == Eigen::NumericalIssue) 
          << std::endl;  // false, no issue detected

// ldlt solver requires PSD, wrong answer
x = A.ldlt().solve(b);
std::cout << x << std::endl;  // Wrong solution [10.625, 1.5, 4.125]
std::cout << b.isApprox(A * x) << std::endl;  // false

// ColPivHouseholderQR doesn't assume PSD, right answer
x = A.colPivHouseholderQr().solve(b);
std::cout << x << std::endl;  // Correct solution [9.25, 4.25, 2.75]
std::cout << b.isApprox(A * x) << std::endl;  // true

注意：更准确地说，可以通过检查 A 是对称的并且 A 的所有特征值 >= 0 来应用 definition of PSD。但正如问题中所述，这可能在计算上很昂贵。

【讨论】：