Una matriz simétrica definida positiva pode ser descomposta como o produto de uma matriz triangular inferior e a transposta da matriz triangular inferior. Uma condição necessária e suficiente para que uma matriz \(A\) qualquer admita feitorização de Cholesky é que seja simétrica e definida positiva. Então, \(A\) pode ser decomposta como
\[A=LL^T,\]
onde \(L\) é una matriz triangular inferior com entradas diagonais estritamente positivas e \(L^T\) representa la conjugada transposta de \(L\).
Algoritmo para a decomposição de Cholesky
Existem vários métodos para gerar a decomposição de Cholesky, uns deles é o seguinte.
Particionar a matriz \(A=LL^T\) da seguinte forma
\[\left[\begin{array}{cc}
a_{11} & A_{21}^T \\
A_{21} & A_{22} \\
\end{array}\right]
=
\left[\begin{array}{cc}
l_{11} & 0 \\
L_{21}^T & L_{22}^T \\
\end{array}\right]
\left[\begin{array}{cc}
l_{11} & L_{21}^T \\
0 & L_{22}^T \\
\end{array}\right]
=
\left[\begin{array}{cc}
l_{11}^2 & l_{11}L_{21}^T \\
l_{11}L_{21} & L_{21} L_{21}^T + L_{22} L_{22}^T\\
\end{array}\right]
\]
1) Determine \(l_{11}\) e \(L_{21}\):
\[l_{11} = \sqrt{a_{11}}, \ L_{21}= \frac{A_{21}}{l_{11}}. \]
2) A partir da igualdade acima,
\[ A_{22} - L_{21} L_{21}^T = L_{22} L_{22}^T, \]
determine \(L_{22} \).
Note que esta é a decomposição de Cholesky de ordem \(n-1\); o tempo de execução do algoritmo é \(\mathcal{O}(n^3)\).
Vamos demonstrar que este algoritmo funciona com matrizes \(A\) positiva definida de ordem \(n\).
Se a é positiva definida, logo
1) \( \ a_{11}>0\)
2) \( \ A_{22} - L_{21} L_{21}^T = A_{22} - \frac{1}{a_{11}} A_{21} A_{2_1}^T \) é positiva definida pelo complemento de Schur.
Se o algoritmo funciona para \(n=m\), funciona para \(n=m-1\); claramente, funciona para \(n=1\).
O complemento de Schur estabelece que se \(A\) é positiva definida, então o complemento de Schur \(S = A_{22} - \frac{1}{a_{11}} A_{21} A_{21}^T\) é positiva definida, onde
\[
A = \left[\begin{array}{cc}
a_{11}^2 & A_{21}^T \\
A_{21} & A_{22} ^T\\
\end{array}\right]
\]
Exemplo
Vamos utilizar a decomposição de Cholesky para criar uma distribuição normal multivariada definido por um vetor de médias \(\mu\) e uma matriz de covariância \(\Sigma\). Obter \( L\) a partir de \(\Sigma\) usando algum algoritmo de decomposição de Cholesky. Para gerar amostras de esta distribuição \(K\)-variada, geramos uma amostra de k observações independentes a partir de uma distribuição normal padrão. Suponha que armazenamos a amostra num vetor \(z\), logo o resultado desejado será
\[ x = \mu + Cz.\]
Por que o vetor \(x\) tem correlação \(\Sigma\)?
Observe que,
\[ E[(x-u)(x-u)^T] = E[ Cz (Cz)^T] = CE[zz^T]C^T = CIC^T = \Sigma,\]
a penúltima igualdade ocorre por que as observações no vetor \(z\) são independentes; por tanto, a covariância é zero, e a variância (autocovariância) é um.
Algoritmo
O exemplo vai ser implementado em Python.
import numpy as np
from scipy.linalg import cholesky
from scipy.stats import norm
tamano_amostra = 10000
# matriz de covariancia desejada
Sigma = np.array([
[ 1.00, 0.75, 0.25],
[ 0.75, 1.00, 0.50],
[ 0.25, 0.50, 1.00]
])
#vetor de medias desejado
u =np.array([[1.],[2.],[3.] ])
#gerar um amostragem a partir de uma variavel normalmente distribuid
#padronizada(media = 0, desvio padrao = 1)
z = norm.rvs(size=(3, tamano_amostra))
# calcular a decomposicao de Cholesky
c = cholesky(Sigma, lower=True)
#convertir a amostra numa distribuicao correlacionada segundo a definicao de sigma
#e media segundo o vetor u
x = np.dot(c, z)
x += u
print('covariancia')
print(np.cov(x))
print('media')
print(np.mean(x, axis=1))
Saida
covariancia
[[ 1.02262798 0.7712547 0.25287991]
[ 0.7712547 1.02521253 0.51647226]
[ 0.25287991 0.51647226 1.01653229]]
media
[ 0.99318368 1.99785094 3.01186588]