Factorización de Cholesky

A factorización ou descomposición de Cholesky é un método alxébrico de descomposición de matrices descuberto polo oficial militar francés e matemático André-Louis Cholesky. Cholesky faleceu polas feridas recibidas no campo de batalla o 31 de agosto de 1918 no norte de Francia. Despois da súa morte, un dos seus compañeiros oficiais, o comandante Benoit, publicou o método de Cholesky para calcular solucións ás ecuacións normais para algúns problemas de axuste de datos de mínimos cadrados en Note sur une méthode de résolution des équations normales provenant de l'application de la méthode des moindres carrés. à un système d'équations lineaires en nombre inferior à aquel des inconnues. Aplicación do método á resolución dun sistema definido de ecuacións lineais (Procede du Commandant Cholesky) Ⓣ, publicado no Bulletin Géodesique en 1924.

Para introducir a factorización de Cholesky é preciso contextualizar que ao resolver un sistema de ecuacións lineais, que defina procesos físicos ou de enxeñaría, é necesario seleccionar o método de resolución, ben sexa directo ou iterativo. A principal desvantaxe dos métodos directos, como o método de eliminación de Gauss,é a súa complexidade computacional. Para unha matriz de tamaño N × N, o método de eliminación gaussiano é de orde $O(n^{3})$ que para matrices moi grandes é un custo computacional caro. Unha solución a este problema é proporcionada polos métodos iterativos como: Jacobi, Gauss-Seidel, SOR ou gradiente conxugado. Estes métodos teñen menor complexidade, o que os fai máis eficientes. Ademais, os erros de redondeo, que poden ocorrer ao traballar con precisión finita na aplicación de métodos directos, fai os métodos iterativos máis apropiados. Na práctica tamén se pode considerar non necesario ter varios díxitos na solución, polo que os métodos das probas iterativas son máis que suficientes para ter unha boa aproximación. No caso específico das matrices simétricas definidas positivas, a factorización de Cholesky demostrou ser un dos condicionadores máis comúns, xa que só require almacenar unha matriz triangular inferior na memoria. Na práctica, moitos servizos de enxeñaría só están interesados en resolver unha matriz de coeficientes simétrica definida positiva, como por exemplo as matrices que xorden da discretización de problemas de dispersión electromagnética.

Definición

Nunha definición xeral cando nunha factorización LU é certo que L=U a factorización chámase factorización de Cholesky. Conceptualmente, obsérvase que esta factorización realiza unha función análoga á raíz cadrada dentro do corpo dos escalares reais positivos. Do mesmo xeito que a raíz cadrada ten limitacións para a súa aplicación dentro dos números reais (o número a tratar debe ser non negativo) a factorización de Cholesky tamén ten certas limitacións en canto á súa aplicación. Diremos que unha matriz A de orde n, definida positiva simétrica (spd), ten factorización de Cholesky se existe unha matriz triangular inferior L con elementos positivos na diagonal, tal que $A=LL^{t}$. Diremos que A é definida positiva se A é simétrica e $x^{t}Ax>0$ para todo x diferente de 0. Doutra maneira podemos dicir que a matriz A é definida positiva se todos os elementos da diagonal son positivos. Entón, a factorización de Cholesky existe e é única.

Método

O factorización de Cholesky dunha matriz $A$ $n\times n$ pódese calcular co seguinte algoritmo recursivo en $n$ pasos. Cada paso produce unha columna de $L$.

A ecuación $A=LL^{T}$ escríbese como:

${\begin{pmatrix}a_{1,1}&A_{2,1}^{T}\\A_{2,1}&A_{2,2}\end{pmatrix}}={\begin{pmatrix}l_{1,1}&0\\L_{2,1}&L_{2,2}\end{pmatrix}}{\begin{pmatrix}l_{1,1}&L_{2,1}^{T}\\0&L_{2,2}^{T}\end{pmatrix}}={\begin{pmatrix}l_{1,1}^{2}&l_{1,1}L_{2,1}^{T}\\l_{1,1}L_{2,1}&L_{2,2}L_{2,2}^{T}+L_{2,1}L_{2,1}^{T}\end{pmatrix}}$

Aquí determínase a primeira columna de $L$:

$l_{1,1}={\sqrt {a_{1,1}}}$ (o elemento da diagonal)

$L_{2,1}={\frac {A_{2,1}}{l_{1,1}}}$ (isto é o resto da columna).

A matriz descoñecida $L_{2,2}$ sae da ecuación:

$A_{2,2}-L_{2,1}L_{2,1}^{T}=L_{2,2}L_{2,2}^{T}$

($L_{2,1}L_{2,1}^{T}$ é o produto vectorial do vector columna $L_{2,1}$ e o vector fila $L_{2,1}^{T}$)

Esta é de novo unha ecuación da forma $A=LL^{T}$, mais agora cunha matriz que ten unha columna e unha fila menos. Agora pódese repetir o anterior para atopar a primeira columna de $L_{2,2}$, que é a segunda columna de $L$, e así sucesivamente.

O método de Cholesky é numéricamente estábel. Na descomposición de Cholesky non é necesario trocar filas ou columnas para evitar obter un número de pivote cero. Mais a orde na que se eliminan as columnas e as filas pode ter un gran impacto no tempo de execución do algoritmo. Este é certamente o caso se a matriz está ocupada de raro en raro (unha matriz dispersa con moitos ceros). Para manter a simetría, no método de Cholesky hai que trocar sempre as filas coas columnas correspondentes. Para iso, mantense unha matriz de permutación $P$, de xeito que a descomposición se expresa como:

$PAP^{T}=LL^{T}$.

Escollendo coidadosamente a orde de eliminación, pódese obter unha factorización de Cholesky $L$ cun cálculo rápido.

Exemplo

Procuramos a descomposición da matriz

$A={\begin{pmatrix}4&12&-16\\12&37&-43\\-16&-43&98\end{pmatrix}}$

Primeiro paso

${\begin{pmatrix}4&12&-16\\12&37&-43\\-16&-43&98\end{pmatrix}}={\begin{pmatrix}l_{1,1}&0&0\\l_{2,1}&l_{2,2}&0\\l_{3,1}&l_{3,2}&l_{3,3}\end{pmatrix}}{\begin{pmatrix}l_{1,1}&l_{2,1}&l_{3,1}\\0&l_{2,2}&l_{3,2}\\0&0&l_{3,3}\end{pmatrix}}$

A primeira columna de $L$ pasa a ser:

$l_{1,1}={\sqrt {4}}=2$

$l_{2,1}=12/2=6$

$l_{3,1}=-16/2=-8$

A matriz 2x2 para o seguinte paso é

${\begin{pmatrix}37&-43\\-43&98\end{pmatrix}}-{\begin{pmatrix}6\\-8\end{pmatrix}}{\begin{pmatrix}6&-8\end{pmatrix}}={\begin{pmatrix}1&5\\5&34\end{pmatrix}}$

Segundo paso

O segundo paso é o mesmo que o primeiro paso:

$l_{2,2}={\sqrt {1}}=1$

$l_{3,2}=5$

O que fica é a ecuación matricial con matrices 1×1:

${\begin{pmatrix}34\end{pmatrix}}={\begin{pmatrix}5&l_{3,3}\end{pmatrix}}{\begin{pmatrix}5\\l_{3,3}\end{pmatrix}}={\begin{pmatrix}25+l_{3,3}^{2}\end{pmatrix}}$

Terceiro e último paso

Da ecuación anterior simplemente despréndese que

$l_{3,3}={\sqrt {9}}=3$

Así que a descomposición de Cholesky de $A$ é

$A={\begin{pmatrix}2&0&0\\6&1&0\\-8&5&3\\\end{pmatrix}}{\begin{pmatrix}2&6&-8\\0&1&5\\0&0&3\\\end{pmatrix}}$

Algoritmo

Unha das formas máis sinxelas de implementar o código de factorización de Cholesky é a seguinte (implementación tendo en conta que A é a matriz asociada ao sistema de ecuacións lineares (SEL) con dimensión $n\times n$):

for (i = 0; i < n; i++) {
    for (j = 0; j <= i; j++) {
        float sum = 0;
        for (k = 0; k < j; k++)
            sum += L[i][k] * L[j][k];

        if (i == j)
            L[i][j] = sqrt(A[i][i] - sum);
        else
            L[i][j] = (1.0 / L[j][j] * (A[i][j] - sum));
    }
}

Observando o código é relativamente sinxelo calcular o custo computacional da factorización de Cholesky. Dado que A é unha matriz de orde n, e tendo en conta que os métodos de substitución cara adiante e cara atrás teñen un custo $n^{2}$, o custo de resolver o SEL mediante a factorización de Cholesky é do tipo:

$(1/3)n^{3}+2n^{2}$

Aplicacións

Se se buscan aplicacións concretas para todo o explicado ata agora, a aplicación máis obvia para a factorización de Cholesky é a resolución de SEL simétrico definido positivo, pero ten varias aplicacións concretas adicionais no mundo da ciencia, a estatística, as comunicacións e a informática:

A. Mínimos cadrados lineais. Este método é unha ferramenta estatística que se encarga de axustar as curvas a un determinado conxunto de datos, optimizando o mínimo do erro cadrado medio total. A factorización de Cholesky pódese usar para calcular os parámetros deste algoritmo. Neste caso é certo que para axustar as curvas hai que resolver o seguinte sistema

$A^{t}Ax=A^{t}b$ sendo $A^{t}A$ spd.

B. Método de Monte Carlo. A método de Monte Carlo é unha técnica cuantitativa que fai uso da estatística e ordenadores para imitar, mediante modelos matemáticos, o comportamento aleatorio de sistemas reais non dinámicos. Mediante a factorización de Cholesky é posíbel xerar ruído segundo a modelaxe dun sistema determinado a partir de ruído non correlacionado, permitindo adaptar a simulación de Monte Carlo a calquera problema específico.

C. Xeración de vectores aleatorios. Relacionado coa simulación (ou método) de Monte Carlo, obsérvase que a factorización de Cholesky pódese usar para xerar vectores aleatorios.

D. Filtro D. Kalman. Este tipo de filtro caracterízase por controlar o estado medio e a covarianza do sistema. Trátase dun tipo de filtro adaptativo que emprega a factorización de Cholesky para extraer datos estatísticos da matriz de covarianza (que sempre reúne as condicións para factorizar) para que os parámetros estatísticos do filtro se poidan axustar en todo momento.

E. Núcleo para o cálculo de factorizacións en matrices dispersas. As matrices raras (ou dispersas) son aquelas nas que un gran número dos seus elementos son nulos. Estas matrices adoitan seguir certos modelos de distribución de elementos distintos de cero e hai varios deles nos que a factorización de Cholesky pode ser útil. O máis claro destes modelos é o daquelas matrices dispersas nas que os elementos non nulos se concentran en bloques e arredor da diagonal principal. A factorización densa de Cholesky úsase nestes casos nos bloques que residen preto da diagonal que son spd.

F. Precondicionamento de métodos iterativos. Ás veces pódese usar a factorización de Cholesky para precondicionar métodos iterativos e mellorar a súa converxencia. O obxectivo final do precondicionamento é aproximar a matriz de precondicionamento o máis preto posible da matriz SEL orixinal para mellorar as súas propiedades numéricas e a converxencia do método numérico.

Véxase tamén

Ligazóns externas

• Implementación e evaluación da Factorización de Cholesky mediante TBB