Factorización de Cholesky incompleta

En el análisis numérico , una factorización de Cholesky incompleta de una matriz definida positiva simétrica es una aproximación dispersa de la factorización de Cholesky . Una factorización de Cholesky incompleta se utiliza a menudo como preacondicionador para algoritmos como el método del gradiente conjugado .

La factorización de Cholesky de una matriz definida positiva A es A = LL * donde L es una matriz triangular inferior . Una factorización de Cholesky incompleta viene dada por una matriz triangular inferior dispersa K que en cierto sentido está cerca de L. El precondicionador correspondiente es KK *.

Una forma popular de encontrar dicha matriz K es usar el algoritmo para encontrar la descomposición de Cholesky exacta en la que K tiene el mismo patrón de escasez que A (cualquier entrada de K se establece en cero si la entrada correspondiente en A también es cero). Esto da una factorización de Cholesky incompleta que es tan escasa como la matriz A.

Motivación

Consideremos la siguiente matriz como ejemplo:

${\displaystyle \mathbf {A} ={\begin{bmatrix}5&-2&0&-2&-2\\-2&5&-2&0&0\\0&-2&5&-2&0\\-2&0&-2&5&-2\\-2&0&0&-2&5\\\end{bmatrix}}}$

Si aplicamos la descomposición regular completa de Cholesky , obtenemos:

${\displaystyle \mathbf {L} ={\begin{bmatrix}2.24&0&0&0&0\\-0.89&2.05&0&0&0\\0&-0.98&2.02&0&0\\-0.89&-0.39&-1.18&1.63&0\\-0.89&-0.39&-0.19&-1.95&0.45\\\end{bmatrix}}}$

Y, por definición:

${\displaystyle \mathbf {A} =\mathbf {L} \mathbf {L'} }$

Sin embargo, al aplicar la descomposición de Cholesky, observamos que algunos elementos cero en la matriz original terminan siendo elementos distintos de cero en la matriz descompuesta, como los elementos (4,2), (5,2) y (5,3) en este ejemplo. Estos elementos se conocen como "relleno".

Esto no es un problema en sí, pero es muy problemático cuando se trabaja con matrices dispersas , ya que la generación de rellenos es mayormente impredecible y reduce la escasez de la matriz, lo que afecta la eficiencia de los algoritmos de matrices dispersas.

Por lo tanto, dada la importancia de la descomposición de Cholesky en los cálculos matriciales, resulta sumamente relevante reutilizar el método regular, de modo de eliminar la generación de rellenos. La factorización de Cholesky incompleta hace exactamente eso, al generar una matriz L similar a la generada por el método regular, pero conservando los elementos cero en la matriz original.

Naturalmente:

${\displaystyle \mathbf {A} \neq \mathbf {L_{ichol}} \mathbf {L_{ichol}'} }$

La multiplicación de la matriz L generada por la factorización de Cholesky incompleta por su transpuesta no producirá la matriz original, sino una similar.

Algoritmo

Para desde hasta : ${\displaystyle i}$ ${\displaystyle 1}$ ${\displaystyle N}$

${\displaystyle L_{ii}=\left({a_{ii}-\sum \limits _{k=1}^{i-1}{L_{ik}^{2}}}\right)^{1 \over 2}}$

Para desde hasta : ${\displaystyle j}$ ${\displaystyle i+1}$ ${\displaystyle N}$

${\displaystyle L_{ji}={1 \over {L_{ii}}}\left({a_{ji}-\sum \limits _{k=1}^{i-1}{L_{ik}L_{jk}}}\right)}$

Implementación

Implementación de la factorización incompleta de Cholesky en el lenguaje GNU Octave . La factorización se almacena como una matriz triangular inferior, con los elementos del triángulo superior configurados a cero.

función  a = ichol ( a ) n = tamaño ( a , 1 );    para  k  =  1 : n a ( k , k )  =  sqrt ( a ( k , k )); para  i  =  ( k + 1 ): n  si  ( a ( i , k )  !=  0 )  a ( i , k )  =  a ( i , k ) / a ( k , k );  finsi finpara para  j  =  ( k + 1 ): n  para  i  =  j : n  si  ( a ( i , j )  !=  0 )  a ( i , j )  =  a ( i , j )  -  a ( i , k ) * a ( j , k );  finsi  finpara finpara finpara para  i  =  1 : n  para  j  =  i + 1 : n  a ( i , j )  =  0 ;  finpara  finpara  finfunción

Ejemplo escaso

Consideremos nuevamente la matriz mostrada al principio de este artículo. Como es simétrica y el método solo utiliza los elementos triangulares inferiores, podemos representarla de la siguiente manera:

${\displaystyle \mathbf {A_{tri}} ={\begin{bmatrix}5&0&0&0&0\\-2&5&0&0&0\\0&-2&5&0&0\\-2&0&-2&5&0\\-2&0&0&-2&5\end{bmatrix}}}$

Más específicamente, en su forma dispersa como una lista de coordenadas, barriendo las filas primero:

Valor 5 -2 -2 -2 5 -2 5 -2 5 -2 5Fila 1 2 4 5 2 3 3 4 4 5 5Columna 1 1 1 1 2 2 3 3 4 4 5

Luego, tomamos la raíz cuadrada de (1,1) y dividimos los otros (i,1) elementos por el resultado:

Valor 2,24 -0,89 -0,89 -0,89 | 5 -2 5 -2 5 -2 5Fila 1 2 4 5 | 2 3 3 4 4 5 5Columna 1 1 1 1 | 2 2 3 3 4 4 5

Después de eso, para todos los demás elementos con columna mayor que 1, calcula (i,j)=(i,j)-(i,1)*(j,1) si (i,1) y (j,1) existen. Por ejemplo: (5,4) = (5,4)-(5,1)*(4,1) = -2 -(-0,89*-0,89) = -2,8.

Valor 2,24 -0,89 -0,89 -0,89 | 4,2 -2 5 -2 4,2 -2,8 4,2Fila 1 2 4 5 | 2 3 3 4 4 5 5Columna 1 1 1 1 | 2 2 3 3 4 4 5 ↑ ↑ ↑ ↑

Los elementos (2,2), (4,4), (5,4) y (5,5) (marcados con una flecha) han sido recalculados, ya que cumplen esta regla. Por otro lado, los elementos (3,2), (3,3) y (4,3) no serán recalculados ya que el elemento (3,1) no existe, aunque los elementos (2,1) y (4,1) sí existan.

Ahora, repita el proceso, pero para (i,2). Tome la raíz cuadrada de (2,2) y divida los otros (i,2) elementos por el resultado:

Valor 2,24 -0,89 -0,89 -0,89 | 2,05 -0,98 | 5 -2 4,2 -2,8 4,2Fila 1 2 4 5 | 2 3 | 3 4 4 5 5Columna 1 1 1 1 | 2 2 | 3 3 4 4 5

Nuevamente, para elementos con columna mayor que 2, calcule (i,j)=(i,j)-(i,2)*(j,2) si (i,2) y (j,2) existen:

Valor 2,24 -0,89 -0,89 -0,89 | 2,05 -0,98 | 4,05 -2 4,2 -2,8 4,2Fila 1 2 4 5 | 2 3 | 3 4 4 5 5Columna 1 1 1 1 | 2 2 | 3 3 4 4 5 ↑

Repita para (i,3). Tome la raíz cuadrada de (3,3) y divida el otro (i,3):

Valor 2,24 -0,89 -0,89 -0,89 2,05 -0,98 | 2,01 -0,99 | 4,2 -2,8 4,2Fila 1 2 4 5 2 3 | 3 4 | 4 5 5Columna 1 1 1 1 2 2 | 3 3 | 4 4 5

Para elementos con columna mayor que 3, calcule (i,j)=(i,j)-(i,3)*(j,3) si existen (i,3) y (j,3):

Valor 2,24 -0,89 -0,89 -0,89 2,05 -0,98 | 2,01 -0,99 | 3,21 -2,8 4,2Fila 1 2 4 5 2 3 | 3 4 | 4 5 5Columna 1 1 1 1 2 2 | 3 3 | 4 4 5 ↑

Repita para (i,4). Tome la raíz cuadrada de (4,4) y divida el otro (i,4):

Valor 2,24 -0,89 -0,89 -0,89 2,05 -0,98 2,01 -0,99 | 1,79 -1,56 | 4,2Fila 1 2 4 5 2 3 3 4 | 4 5 | 5Columna 1 1 1 1 2 2 3 3 | 4 4 | 5

Para elementos con columnas mayores que 4, calcule (i,j)=(i,j)-(i,4)*(j,4) si existen (i,4) y (j,4):

Valor 2,24 -0,89 -0,89 -0,89 2,05 -0,98 2,01 -0,99 | 1,79 -1,56 | 1,76Fila 1 2 4 5 2 3 3 4 | 4 5 | 5Columna 1 1 1 1 2 2 3 3 | 4 4 | 5 ↑

Finalmente toma la raíz cuadrada de (5,5):

Valor 2,24 -0,89 -0,89 -0,89 2,05 -0,98 2,01 -0,99 1,79 -1,56 | 1,33Fila 1 2 4 5 2 3 3 4 4 5 | 5Columna 1 1 1 1 2 2 3 3 4 4 | 5

Ampliando la matriz a su forma completa:

${\displaystyle \mathbf {L} ={\begin{bmatrix}2.24&0&0&0&0\\-0.89&2.05&0&0&0\\0&-0.98&2.01&0&0\\-0.89&0&-0.99&1.79&0\\-0.89&0&0&-1.56&1.33\end{bmatrix}}}$

Obsérvese que en este caso no se generaron rellenos en comparación con la matriz original. Los elementos (4,2), (5,2) y (5,3) siguen siendo cero.

Sin embargo, si realizamos la multiplicación de L por su transpuesta:

${\displaystyle \mathbf {LL'} ={\begin{bmatrix}5&-2&0&-2&-2\\-2&5&-2&0.8&0.8\\0&-2&5&-2&0\\-2&0.8&-2&5&-2\\-2&0.8&0&-2&5\end{bmatrix}}}$

Obtenemos una matriz ligeramente diferente a la original, ya que la descomposición no tuvo en cuenta todos los elementos, para eliminar los rellenos.

Implementación dispersa

La versión dispersa de la factorización de Cholesky incompleta (el mismo procedimiento presentado anteriormente) implementada en MATLAB se puede ver a continuación. Naturalmente, MATLAB tiene sus propios medios para tratar con matrices dispersas, pero el código a continuación se hizo explícito con fines pedagógicos. Este algoritmo es eficiente, ya que trata la matriz como una matriz secuencial unidimensional, evitando automáticamente los elementos cero.

función  A = Sp_ichol ( A ) n = tamaño ( A , 1 ); ncols = A ( n ). col ; c_end = 0 ; para col = 1 : ncols is_next_col = 0 ; c_start = c_end + 1 ; para i = c_start : n si A ( i ). col == col % en la columna actual (col): si A ( i ). col == A ( i ). row A ( i ). val = sqrt ( A ( i ). val ); % toma la raíz cuadrada del elemento diagonal de la columna actual div = A ( i ). val ; de lo contrario A ( i ). val = A ( i ). val / div ; % divide los otros elementos de la columna actual por la raíz cuadrada del elemento diagonal end end si A ( i ). col > col % en las siguientes columnas (col+1 ... ncols): if is_next_col == 0 c_end = i - 1 ; is_next_col = 1 ; end v1 = 0 ; v2 = 0 ; for j = c_start : c_end if A ( j ). col == col if A ( j ). row == A ( i ). row % busca los elementos A(j) de la columna actual (col) cuyo índice de fila sea el mismo que el índice de fila A(i) del elemento actual v1 = A ( j ). val ;                                        fin si A ( j ). fila == A ( i ). col % busca los elementos de la columna actual (col) A(j) cuyo índice de fila es el mismo que el índice de la columna del elemento actual A(i) v2 = A ( j ). val ; fin si v1 ~= 0 && v2 ~= 0 % si estos elementos existen en la columna actual (col), recalcula el elemento actual A(i): A ( i ). val = A ( i ). val - v1 * v2 ; break ; fin fin fin fin fin fin fin fin                  

Referencias

• Factorización de Cholesky incompleta en la wiki de CFD Online
• Golub, Gene H. ; Van Loan, Charles F. (1996), Cálculos matriciales (3.ª ed.), Johns Hopkins, ISBN 978-0-8018-5414-9. Véase la Sección 10.3.2.