Multiplication de matrice très lente en Eigen

Question

J'ai implémenté un processus d'optimisation de Gauss-Newton qui consiste à calculer l'incrément en résolvant un système linéarisé Hx = b. Le H matrx est calculé par H = J.transpose() * W * J et b est calculé à partir de b = J.transpose() * (W * e) où e est le vecteur d'erreur. Jacobien est une matrice n-by-6 où n est en milliers et reste inchangé à travers les itérations et W est une matrice de poids diagonale n-by-n qui changera à travers les itérations (certains éléments diagonaux seront mis à zéro). Cependant j'ai rencontré un problème de vitesse. Lorsque je n'ajoute pas la matrice de poids W, à savoir H = J.transpose()*J et b = J.transpose()*e, mon processus Gauss-Newton peut fonctionner très rapidement en 0.02 sec pour 30 itérations. Cependant quand j'ajoute la matrice W qui est définie en dehors de la boucle d'itération, elle devient si lente (0.3 ~ 0.7 sec pour 30 itérations) et je ne comprends pas si c'est mon problème de codage ou cela prend normalement beaucoup de temps.

Tout ici sont des matrices et des vecteurs propres.

J'ai défini ma matrice W en utilisant la fonction .asDiagonal() dans la bibliothèque Eigen à partir d'un vecteur de variances inverses. alors juste utilisé dans le calcul pour H ad b. Ensuite, il devient très lent. Je voudrais avoir quelques indices sur les raisons potentielles de ce ralentissement énorme.

EDIT:

Il n'y a que deux matrices. Jacobian est définitivement dense. La matrice de poids est générée à partir d'un vecteur par la fonction vec.asDiagonal() qui provient de la bibliothèque dense, donc je suppose qu'elle est également dense.

Le code est vraiment simple et la seule différence qui cause le changement d'heure est l'addition de la matrice de poids. Voici un extrait de code:

for (int iter=0; iter<max_iter; ++iter) { 
    // obtain error vector 
    error = ... 
    // calculate H and b - the fast one 
    Eigen::MatrixXf H = J.transpose() * J; 
    Eigen::VectorXf b = J.transpose() * error; 
    // calculate H and b - the slow one 
    Eigen::MatrixXf H = J.transpose() * weight_ * J; 
    Eigen::VectorXf b = J.transpose() * (weight_ * error); 
    // obtain delta and update state 
    del = H.ldlt().solve(b); 
    T <- T(del) // this is pseudo code, meaning update T with del 
} 

Il est en fonction dans une classe, et de la matrice de poids maintenant à des fins de débogage est défini comme une variable de classe qui peut être accessible par la fonction et est définie avant que la fonction est appelée .

Answer 1

Parce que la multiplication de la matrice est juste la diagonale, vous pouvez le modifier à utiliser le coefficient de multiplication sage comme ceci:

MatrixXd m; 
VectorXd w; 
w.setLinSpaced(5, 2, 6); 

m.setOnes(5,5); 

std::cout << (m.array().rowwise() * w.array().transpose()).matrix() << "\n"; 

De même, le produit vectoriel de la matrice peut être écrit comme:

(w.array() * error.array()).matrix() 

Cela évite les éléments zéro dans la matrice. Sans un MCVE pour me baser sur cela, YMMV ...

Answer 2

Je suppose que weight_ est déclaré comme MatrixXf dense? Si oui, alors le remplacer par w.asDiagonal() partout où vous utilisez weight_, ou faire le plus tard un alias à l'expression asDiagonal:

auto weight = w.asDiagonal(); 

De cette façon Eigen sait que weight est une matrice diagonale et les calculs seront optimisés comme prévu.