特征 - 检查矩阵是否为正(半)确定

Question

我正在实现一个谱聚类算法,我必须确保矩阵(laplacian)是正半正定的.

检查矩阵是否为正定(PD)就足够了,因为可以在特征值中看到"半"部分.矩阵非常大(nxn,其中n约为数千),因此特征分析是昂贵的.

在Eigen中是否有任何检查可以在运行时产生bool结果？

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

Answer 1

您可以使用Cholesky分解(LLT),Eigen::NumericalIssue如果矩阵为负数则返回,请参阅文档.

示例如下:

#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）矩阵输入。事实上，存在通过测试的非对称非PSD矩阵。考虑以下示例（数字取自《九章算术》，第 8 章，第 1题）：AA_llt.info() != Eigen::NumericalIssue

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 的特征值 >= 0来应用PSD的定义。A但正如问题中提到的，这可能在计算上很昂贵。