Étant donné une matrice symétrique A et un vecteur x, j'ai souvent besoin de calculer x^T * A * x. Je peux le faire Eigen ++ 3.x avec x.transpose() * (A * x), mais cela n'exploite pas l'information selon laquelle x est la même des deux côtés et A est symétrique. Existe-t-il un moyen plus efficace de calculer cela?

1
quant_dev 3 sept. 2020 à 15:58

2 réponses

Meilleure réponse

Pour les petites matrices, écrire moi-même la boucle for semble être plus rapide que de s'appuyer sur le code d'Eigen. Pour les grandes matrices, j'ai obtenu de bons résultats en utilisant .selfAdjointView:

double xAx_symmetric(const Eigen::MatrixXd& A, const Eigen::VectorXd& x)
{
    const auto dim = A.rows();
    if (dim < 15) {
        double sum = 0;
        for (Eigen::Index i = 0; i < dim; ++i) {
            const auto x_i = x[i];
            sum += A(i, i) * x_i * x_i;
            for (Eigen::Index j = 0; j < i; ++j) {
                sum += 2 * A(j, i) * x_i * x[j];
            }
        }
        return sum;
    }
    else {
        return x.transpose() * A.selfadjointView<Eigen::Upper>() * x;
    }
}
2
quant_dev 5 sept. 2020 à 00:17

À quelle fréquence calculez-vous cela? Si très souvent pour différents x, alors cela pourrait donner un peu d'accélération pour calculer une décomposition Cholesky ou LDLT de A et utiliser que le produit d'une matrice triangulaire avec un vecteur n'a besoin que de la moitié du multiplications.

Ou peut-être même plus simple, si vous décomposez A=L+D+L.T, où L est strictement triangulaire inférieur et D est diagonal, alors

x.T*A*x = x.T*D*x + 2*x.T*(L*x)

Où le premier terme est la somme sur d[k]*x[k]**2. Le deuxième terme, s'il utilise soigneusement la structure triangulaire, utilise la moitié des opérations de l'expression d'origine.

Si le produit matrice-vecteur triangulaire doit être implémenté en dehors des procédures Eigen, cela pourrait détruire l'efficacité / les optimisations des opérations de bloc de type BLAS dans le produit vecteur matrice générique. En fin de compte, cette réduction du nombre d'opérations arithmétiques pourrait ne pas être améliorée.

2
Lutz Lehmann 3 sept. 2020 à 14:32