Eigen

Question

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

检查矩阵是否为正定 (PD) 就足够了，因为可以在特征值中看到 "semi-" 部分。矩阵非常大（nxn，其中 n 是几千的数量级）所以特征分析很昂贵。

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

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

Answer 1

您可以使用 Cholesky 分解 (LLT)，其中 returns 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!");
}

这项检查很重要，例如一些 Eigen 求解器（如 LTDT) requires PSD(or NSD) matrix input. In fact, there exists non-symmetric and hence non-PSD matrix A that passes the A_llt.info() != Eigen::NumericalIssue test. Consider the following example (numbers taken from 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。但是正如问题中提到的，这可能是计算上的贵。

Answer 3

您必须测试矩阵是否对称 (A.isApprox(A.transpose()))，然后创建 LDLT（而不是 LLT，因为 LDLT 处理特征值之一为 0 的情况，即不是严格正向的），然后测试数值问题和正向性：

template <class MatrixT>
bool isPsd(const MatrixT& A) {
  if (!A.isApprox(A.transpose())) {
    return false;
  }
  const auto ldlt = A.template selfadjointView<Eigen::Upper>().ldlt();
  if (ldlt.info() == Eigen::NumericalIssue || !ldlt.isPositive()) {
    return false;
  }
  return true;
}

我在

上测试过这个

1 2
2 3

具有负特征值（因此不是 PSD）。如果没有 isPositive() 测试，isPsd() 错误地 returns 为真。

以后

1 2
2 4

具有空特征值（因此 PSD 但不是 PD）。

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

Eigen - Check if matrix is Positive (Semi-)Definite

c++

eigen3