2011-09-22 9 views
6

Ho due matrici quadrate A e B. A è simmetrico, B è simmetrico positivo definito. Vorrei calcolare $ trace (A.B^{- 1}) $. Per ora, computo la scomposizione di Cholesky di B, risolvo per C nell'equazione $ A = C.B $ e riassumo gli elementi diagonali.calcolo efficiente di Trace (AB^{- 1}) dati A e B

Esiste un modo più efficiente di procedere?

Ho intenzione di utilizzare Eigen. Potresti fornire un'implementazione se le matrici sono sparse (A può essere spesso diagonale, B è spesso a banda diagonale)?

+0

Penso che il tag C++ appartenga effettivamente qui, poiché la domanda riguarda un'implementazione utilizzando Eigen, una libreria di manipolazione della matrice C++. –

+0

È una semidefinita positiva o definita positiva? –

+0

@DavidZaslavsky Ho rimosso il tag – yannick

risposta

5

Se B è scarsa, può essere efficace (cioè, O (n), assumendo buon numero condizioni di B) per risolvere per x_i in

B x_i = a_i 

(campione Conjugate Gradient codice è fornito su Wikipedia). Prendendo a_i per essere i vettori di colonna di A, si ottiene la matrice B^{-1} A in O (n^2). Quindi è possibile sommare gli elementi diagonali per ottenere la traccia. In genere, è più semplice eseguire questa moltiplicazione inversa sparsa che ottenere l'intero set di autovalori. Per il confronto, Cholesky decomposition è O (n^3). (vedere il commento di Darren Engwirda in basso su Cholesky).

Se avete solo bisogno un'approssimazione alla traccia, si può effettivamente ridurre il costo a O (q n) dalla media

r^T (A B^{-1}) r 

oltre q vettori aleatori r. Solitamente q << n. Questa è una stima imparziale purché le componenti del vettore casuale r soddisfano

< r_i r_j > = \delta_{ij} 

dove <...> indica una media per la distribuzione di r. Ad esempio, i componenti r_i potrebbero essere distribuiti gaussiani indipendenti con varianza unità. Oppure potrebbero essere selezionati uniformemente da + -1. In genere la traccia viene ridimensionata come O (n) e l'errore nella stima di traccia viene ridimensionato come O (sqrt (n/q)), pertanto l'errore relativo viene ridimensionato come O (sqrt (1/nq)).

+0

Grazie per la risposta. Come si fa la media con r? da quello che scrivi, sembra che tu abbia bisogno di calcolare A.B^{- 1} che probabilmente non è quello che volevi dire. – yannick

+0

Kipton probabilmente significa che dovresti calcolare r^T A B^{- 1} r risolvendo prima B x = r e poi calcolando r^T A x. Ma non vedo come ottenga un costo di O (n) per l'approccio probabilistico: risolvendo n sistemi con costo O (n) ognuno dà un costo di O (n^2). Forse il numero di vettori casuali può essere preso più piccolo di n = dimensione di A? –

+0

@ Jitse, sì, grazie per aver trovato l'errore di battitura. –

1

Se gli autovalori generalizzati sono più efficienti da calcolare, è possibile calcolare gli autovalori generalizzati, A*v = lambda* B *v e quindi sommare tutti i lambda.

Problemi correlati