Résultats erronés du solveur clairsemé propre

Question

J'essaie de résoudre un système linéaire clairsemé Ax = B avec la bibliothèque Eigen en C ++, mais l'exemple trivial suivant semble donner une solution incorrecte:

#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;
}

Je ne vois aucune erreur, l'algorithme renvoie "0" signifiant "Succès", mais la solution que j'obtiens est

x = 0 0.193524 0 0.193524 0.193524 0 0 0 0

ce qui n'est évidemment pas la solution à ce système, la bonne est

x = 0 0 0 0 0.0967621 0 0 0 0

Answer 1

Voici la documentation pour le SimplicialLDLTsolveur:

Cette classe fournit une factorisation LDL ^ T Cholesky sans racine carrée de matrices clairsemées qui sont auto-jointes et définies positives .

Lorsque la matrice stocke des nombres réels dans les éléments, auto-adjoint == symétrique. Votre matrice n'est clairement pas symétrique. De plus, toutes les matrices symétriques ne sont pas définies positivement, voir les exemples .

En bref, le solveur que vous avez choisi ne s'applique qu'à une classe très étroite de matrices. Comme vous l'avez déjà découvert, le SparseLUsolveur fonctionne pour vos données d'entrée.

ConjugateGradientle solveur ne fonctionnera pas non plus, il ne nécessite pas que la matrice soit définie positivement, mais il nécessite qu'elle soit auto-adjointe.