Eigen - проверьте, является ли матрица положительной (полу) определенной

Я реализую алгоритм спектральной кластеризации, и я должен убедиться, что матрица (лапласиан) является положительно полуопределенной.

Достаточно проверить, является ли матрица положительно определенной (ПД), поскольку в собственных значениях можно увидеть «полу-» часть. Матрица довольно большая (nxn, где n порядка нескольких тысяч), поэтому собственный анализ стоит дорого.

Есть ли в Eigen какая-либо проверка, которая дает логический результат во время выполнения?

Matlab может дать результат с помощью метода chol(), выдав исключение, если матрица не является PD. Следуя этой идее, Eigen возвращает результат, не жалуясь на LLL.llt().matrixL(), хотя я ожидал некоторого предупреждения/ошибки. У Eigen также есть метод isPositive, но из-за ошибки его нельзя использовать. для систем со старой версией Eigen.


c++ eigen eigen3
person dim_tz    schedule 05.02.2016    source источник
comment
Разве вы не можете сначала проверить, является ли он эрмитовым, а затем посмотреть на собственные значения? Проверить герметичность несложно.   -  person vsoftco    schedule 05.02.2016
comment
Вы правы насчет эрмитовой части, но в идеале я хотел бы избежать многократного вычисления собственных значений для огромной матрицы, поскольку это мой желаемый результат, поэтому я бы хотел, чтобы это произошло только один раз, если это возможно.   -  person dim_tz    schedule 05.02.2016
comment
Возможно, вы можете попробовать декомпозицию Холецкого из Eigen, которая возвращает NumericalIssue, если матрица отрицательна, см. eigen. tuxfamily.org/dox/classEigenNumericalIssue1LLT.html   -  person vsoftco    schedule 05.02.2016
comment
Возможно проблема с версией, потому что: error: 'class Eigen::LDLT has no member named 'info'   -  person dim_tz    schedule 05.02.2016
comment
О, на самом деле это работает для LLT, но не для LDLT, спасибо за указатель! Если вы хотите, вы можете написать ответ, чтобы я принял это, в противном случае я опубликую фрагмент в качестве ответа позже.   -  person dim_tz    schedule 05.02.2016
comment
Рад, что это помогло, я написал ответ. Я не уверен, что это лучший способ, тем более что вы действительно заинтересованы в том, чтобы наименьшее собственное значение было больше или равно нулю (и эрмитичность, конечно), но это делает работу.   -  person vsoftco    schedule 05.02.2016

Ответы (2)


arrow_upward
16
arrow_downward

Вы можете использовать разложение Холецкого (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!");
    }    
}
person vsoftco    schedule 05.02.2016
comment
Сбивает с толку, потому что в этой документации также говорится: Помните, что разложения Холецкого не раскрывают ранг. Это разложение LLT устойчиво только на положительно определенных матрицах, вместо этого используйте LDLT для полуопределенного случая. Я думаю, это зависит от того, как вы определяете раскрытие ранга. - person Taylor; 14.06.2019

arrow_upward
5
arrow_downward

В дополнение к ответу @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 (например, 1LDLT.html" rel="noreferrer">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

Примечания: чтобы быть более точным, можно применить определение PSD проверив, что A является симметричным, а все собственные значения A >= 0. Но, как уже упоминалось в вопросе, это может быть дорогостоящим в вычислительном отношении.

person Yixing    schedule 07.02.2019