Si B
est clairsemée, il peut être efficace (à savoir O (n), en supposant un bon nombre de condition de B
) à résoudre pour x_i
dans
B x_i = a_i
(échantillon Conjugate Gradient code est donnée sur Wikipedia). Prenant a_i
pour être les vecteurs colonne de A
, vous obtenez la matrice B^{-1} A
dans O (n^2). Ensuite, vous pouvez additionner les éléments diagonaux pour obtenir la trace. Généralement, il est plus facile de faire cette multiplication inverse éparse que d'obtenir l'ensemble complet des valeurs propres. Pour comparaison, Cholesky decomposition est O (n^3). (voir le commentaire de Darren Engwirda ci-dessous sur Cholesky).
Si vous avez seulement besoin d'une approximation à la trace, vous pouvez réellement réduire le coût à O (q n) en faisant la moyenne
r^T (A B^{-1}) r
sur q
vecteurs aléatoires r
. Habituellement q << n
. Ceci est une estimation non biaisée à condition que les composantes du vecteur aléatoire r
satisfont
< r_i r_j > = \delta_{ij}
où <...>
indique une moyenne sur la répartition des r
. Par exemple, les composants r_i
peuvent être des distributions gaussiennes indépendantes avec variance unitaire. Ou ils pourraient être sélectionnés uniformément à partir de + -1. Typiquement les échelles de trace comme O (n) et l'erreur dans l'estimation de trace échelles comme O (sqrt (n/q)), de sorte que l'erreur relative échelles comme O (sqrt (1/nq)).
Je pense que la balise C++ appartient réellement ici, car la question concerne une implémentation utilisant Eigen, une bibliothèque de manipulation de matrice C++. –
Est-ce qu'un semi-défini positif ou positif défini? –
@DavidZaslavsky J'ai retiré l'étiquette – yannick