艾根 - 检查矩阵为正（半）定

Question

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

如果矩阵是正定的（PD）的检查是足够的，因为“半”的一部分可以在特征值中可以看出。该矩阵是相当大的（n×n个，其中n是在几千量级），所以特征分析是昂贵的。

是否有任何的本征检查，让在运行时布尔结果呢？

Matlab的可通过如果一个矩阵不是PD抛出异常给与chol()方法的结果。依照该思路，艾根返回一个结果，而不抱怨的LLL.llt().matrixL()，但我期待一些警告/错误。本征也有方法isPositive，但是由于bug它是不能用于与旧版本的本征系统。

Answer 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!");
    }    
}

Answer 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!");
}    

这种检查是很重要的，例如一些本征求解器（如LTDT）要求PSD（或NSD）矩阵输入。事实上，有存在经过A测试非对称的，并且因此非PSD矩阵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

注：更准确地说，人们可以申请通过检查definition of PSD的A是对称的，所有A的特征值> = 0，但在问题中提到，这可能是计算昂贵的。