2011-09-22 3 views
6

J'ai deux matrices carrées A et B. A est symétrique, B est symétrique positif défini. Je voudrais calculer $ trace (A.B^{- 1}) $. Pour l'instant, je calcule la décomposition de Cholesky de B, résous C dans l'équation $ A = C.B $ et résume les éléments diagonaux.calcul efficace de Trace (AB^{- 1}) donné A et B

Y a-t-il une façon plus efficace de procéder?

Je prévois utiliser Eigen. Pourriez-vous fournir une implémentation si les matrices sont éparses (A peut souvent être en diagonale, B est souvent en bande diagonale)?

+0

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++. –

+0

Est-ce qu'un semi-défini positif ou positif défini? –

+0

@DavidZaslavsky J'ai retiré l'étiquette – yannick

Répondre

5

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} 

<...> 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)).

+0

Merci pour votre réponse. Comment faites-vous la moyenne avec r? D'après ce que vous écrivez, il semble que vous ayez besoin de calculer A.B^{- 1} ce qui n'est probablement pas ce que vous vouliez dire. – yannick

+0

Kipton signifie probablement que vous devez calculer r^T A B^{- 1} r en résolvant d'abord B x = r et ensuite en calculant r^T A x. Mais je ne vois pas comment il obtient un coût de O (n) pour l'approche probabiliste: la résolution de n systèmes avec coût O (n) donne chacun un coût de O (n^2). Peut-être que le nombre de vecteurs aléatoires peut être pris plus petit que n = taille de A? –

+0

@Jitse, oui, merci de trouver la faute de frappe. –

1

Si les valeurs propres généralisées sont plus efficaces à calculer, vous pouvez calculer les valeurs propres généralisées, A*v = lambda* B *v, puis additionner tous les lambdas.

Questions connexes