Eigen risultati errati del risolutore sparse

Aug 26 2020

Sto cercando di risolvere un sistema lineare sparso Ax=B con la libreria Eigen in C++, tuttavia il seguente esempio banale sembra fornire una soluzione errata:

#include <Eigen/SparseCholesky>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <iostream>
#include <vector>

using namespace std;
using namespace Eigen;

int main(){

    SimplicialLDLT<SparseMatrix<double>> solver;
    SparseMatrix<double> A(9,9);
    typedef Triplet<double> T;
    vector<T> triplet;
    VectorXd B(9);

    for(int i=0; i<4; i++){
        triplet.push_back(T(i,i,1));
        triplet.push_back(T(i+5,i+5,1));
    }

    triplet.push_back(T(4,1,-1));
    triplet.push_back(T(4,3,-1));
    triplet.push_back(T(4,5,-1));
    triplet.push_back(T(4,7,-1));
    triplet.push_back(T(4,4,4));

    A.setFromTriplets(triplet.begin(),triplet.end());
    B << 0,0,0,0,0.387049,0,0,0,0;

    solver.compute(A);
    VectorXd x = solver.solve(B);

    cout << "A\n" << A << "\n";
    cout << "B\n" << B << "\n";
    cout << "x\n" << x << "\n";

    return 0;
}

Non vedo errori, l'algoritmo restituisce "0" che significa "Successo", tuttavia la soluzione che ottengo è

x = 0 0.193524 0 0.193524 0.193524 0 0 0 0

che ovviamente non è la soluzione a questo sistema, quella corretta lo è

x = 0 0 0 0 0.0967621 0 0 0 0

Risposte

4 Soonts Aug 26 2020 at 04:24

Ecco la documentazione per il SimplicialLDLTrisolutore:

Questa classe fornisce fattorizzazioni di Cholesky LDL^T senza radice quadrata di matrici sparse autoaggiunte e definite positive .

Quando la matrice memorizza numeri reali negli elementi, autoaggiunto == simmetrico. La tua matrice chiaramente non è simmetrica. Inoltre, non tutte le matrici simmetriche sono definite positivamente, vedi esempi .

In breve, il risolutore che hai scelto è applicabile solo a classi di matrici molto ristrette. Come hai già scoperto, SparseLUil risolutore funziona per i tuoi dati di input.

ConjugateGradientnemmeno il risolutore funzionerà, non richiede che la matrice sia definita positiva ma richiede che sia autoaggiunta.