Kann ich eine Cholesky-Zerlegung aus der Umkehrung einer Matrix erhalten?

7

Ich habe die Umkehrung einer riesigen Kovarianzmatrix, aus der ich zufällige Instanzen ziehen möchte.

Ich weiß, wie das geht, indem ich eine Cholesky-Zerlegung in der Kovarianzmatrix durchführe und damit einen Vektor unabhängiger Gaußscher transformiere. Der einfache Prozess besteht also darin, die Riesenmatrix zu invertieren und dann die Cholesky-Zerlegung durchzuführen.

Da die Cholesky-Zerlegung in erster Linie verwendet werden kann, um Inversen zu erhalten, scheint es einen Weg zu geben, das zu bekommen, was ich brauche, ohne den Inversionsschritt (der lange dauert).

Frage: Gibt es eine Möglichkeit, die Cholesky-Zerlegung zu verwenden, ohne die Matrix zu invertieren?

Chuck37
quelle

Antworten:

2

Sie können das Invertieren der Matrix vermeiden, indem Sie mit der Eigendekompositionsmethode Zeichnungen zeichnen. Nach dieser Methode werden die Ziehungen mit diesem Produkt generiert:

(V.D.)X.,

Dabei ist die Eigenvektoren der Matrix, eine Diagonalmatrix mit den Quadratwurzeln der Eigenwerte und eine Matrix mit Zeichnungen aus der univariaten Standardverteilung .V.D.X.N.(0,1)

Es ist einfach, diese Methode anzupassen und zu vermeiden, die ursprüngliche Kovarianzmatrix wiederherzustellen, indem folgende Ergebnisse verwendet werden: 1) Die Eigenwerte einer Matrix sind der Kehrwert der Eigenwerte ihrer inversen ; 2) Die Eigenvektoren einer Matrix A sind auch Eigenvektoren ihrer inversen .EINEIN- -1EIN- -1


Beispiel:

Angenommen, die ursprüngliche Kovarianzmatrix lautet wie folgt:

EIN=[10,8- -0,40,820,3- -0,40,33]].
Aber Sie haben die Umkehrung dieser Matrix, : B.=EIN- -1
EIN- -1=B.=[1,699- -0,7250,299- -0,7250,817- -0,1780,299- -0,1780,391]].

Die auf der ursprünglichen Matrix basierende Eigendekompositionsmethode ergibt die folgende Kovarianzmatrix für eine Stichprobe von Zeichnungen:EIN

A <- rbind(c(1,0.8,-0.4), c(0.8,2,0.3), c(-0.4,0.3,3))
e1 <- eigen(A, symmetric=TRUE)
set.seed(1)
X <- matrix(rnorm(5000*ncol(A)), ncol=ncol(A))
draws1 <- t(e1$vectors %*% sqrt(diag(e1$values)) %*% t(X))
draws1.cov <- cov(draws1)
draws1.cov
#           [,1]      [,2]       [,3]
#[1,]  0.9765023 0.8030752 -0.3970233
#[2,]  0.8030752 1.9941052  0.3229827
#[3,] -0.3970233 0.3229827  3.1689348

Mit der Matrix, die Sie haben (die Umkehrung von A), müssen Sie nur die Eigenwerte invertieren:

B <- solve(A)
e2 <- eigen(B, symmetric=TRUE)
e2$values <- 1/e2$values
draws2 <- t(e2$vectors %*% sqrt(diag(e2$values)) %*% t(X))
draws2.cov <- cov(draws2)
draws2.cov
#           [,1]      [,2]       [,3]
#[1,]  0.9765023 0.8030752 -0.3970233
#[2,]  0.8030752 1.9941052  0.3229827
#[3,] -0.3970233 0.3229827  3.1689348

Man erhält eine Kovarianzmatrix, die eng mit der ursprünglichen übereinstimmt, und wir mussten nicht invertieren B, um die ursprüngliche Kovarianzmatrix wiederherzustellen A.


Eine kleine Simulation zur Überprüfung der Gültigkeit dieses Ansatzes:

set.seed(3)
niter <- 1000
m <- matrix(0, nrow=ncol(A), ncol=ncol(A))
for (i in seq_len(niter))
{
    X <- matrix(rnorm(5000*ncol(A)), ncol=ncol(A))
    draws2 <- t(e2$vectors %*% sqrt(diag(e2$values)) %*% t(X))
    m <- m + cov(draws2)
}
m/niter
# average covariance matrix
#           [,1]      [,2]       [,3]
#[1,]  1.0005129 0.7995872 -0.4005644
#[2,]  0.7995872 1.9993231  0.2990850
#[3,] -0.4005644 0.2990850  2.9957277
# original covariance matrix 'A'
#    [,1] [,2] [,3]
#[1,]  1.0  0.8 -0.4
#[2,]  0.8  2.0  0.3
#[3,] -0.4  0.3  3.0

Wir können sehen, dass die Kovarianzmatrix der Draws im Durchschnitt sehr nahe an der ursprünglichen Kovarianzmatrix liegt A.

javlacalle
quelle
Unabhängige Anfrage: Ich habe vorgeschlagen, arma zu einem Synonym für arima zu machen, und dies vor einiger Zeit hier auf Meta beworben (ohne großen Erfolg). Als einer der Topscorer unter dem Arima- Tag können Sie hier über den Vorschlag abstimmen .
Richard Hardy