Descomposición colérica

Si una matriz hermitiana A es solo semidefinida positiva, en lugar de definida positiva, entonces aún tiene una descomposición de la forma A = LL * donde las entradas diagonales de L pueden ser cero. ^{[ 7 ]} La descomposición no tiene por qué ser única, por ejemplo: $${\displaystyle {\begin{bmatrix}0&0\\0&1\end{bmatrix}}=\mathbf {L} \mathbf {L} ^{*},\quad \quad \mathbf {L} ={\begin{bmatrix}0&0\\\cos \theta &\sin \theta \end{bmatrix}},}$$ para cualquier θ . Sin embargo, si el rango de A es r , entonces existe un único triángulo inferior L con exactamente r elementos diagonales positivos y n − r columnas que contienen todos ceros. ^{[ 8 ]}

Alternativamente, la descomposición puede hacerse única cuando se fija una elección pivotante. Formalmente, si A es una matriz semidefinida positiva n × n de rango r , entonces existe al menos una matriz de permutación P tal que PAP ^T tiene una descomposición única de la forma PAP ^T = LL ^* con ${\textstyle \mathbf {L} ={\begin{bmatrix}\mathbf {L} _{1}&0\\\mathbf {L} _{2}&0\end{bmatrix}}}$, donde L ₁ es una matriz triangular inferior r × r con diagonal positiva. ^{[ 9 ]}

La descomposición de Cholesky se utiliza principalmente para la solución numérica de ecuaciones lineales.${\textstyle \mathbf {Ax} =\mathbf {b} }$Si A es simétrica y definida positiva, entonces ${\textstyle \mathbf {Ax} =\mathbf {b} }$se puede resolver calculando primero la descomposición de Cholesky ${\textstyle \mathbf {A} =\mathbf {LL} ^{\mathrm {*} }}$, luego resolviendo${\textstyle \mathbf {Ly} =\mathbf {b} }$para y mediante sustitución hacia adelante y finalmente resolviendo${\textstyle \mathbf {L^{*}x} =\mathbf {y} }$para x por sustitución hacia atrás .

Una forma alternativa de eliminar el cálculo de raíces cuadradas en el${\textstyle \mathbf {LL} ^{\mathrm {*} }}$La descomposición consiste en calcular la descomposición LDL.${\textstyle \mathbf {A} =\mathbf {LDL} ^{\mathrm {*} }}$, luego resolviendo${\textstyle \mathbf {Ly} =\mathbf {b} }$para y , y finalmente resolviendo${\textstyle \mathbf {DL} ^{\mathrm {*} }\mathbf {x} =\mathbf {y} }$.

Para sistemas lineales que pueden expresarse en forma simétrica, la descomposición de Cholesky (o su variante LDL) es el método preferido, debido a su mayor eficiencia y estabilidad numérica . En comparación con la descomposición LU , es aproximadamente el doble de eficiente. ^{[ 2 ]}

En el problema de mínimos cuadrados lineales se busca una solución x de un sistema sobredeterminado Ax = l , tal que la norma cuadrática del vector residual Ax-l sea mínima. Esto se puede lograr resolviendo mediante la descomposición de Cholesky de ecuaciones normales.${\displaystyle \mathbf {Nx} =\mathbf {A} ^{\mathsf {T}}\mathbf {l} }$, dónde${\displaystyle \mathbf {N} =\mathbf {A} ^{\mathsf {T}}\mathbf {A} }$es simétrica definida positiva. La matriz de la ecuación simétrica también puede provenir de un funcional de energía, que debe ser positivo por consideraciones físicas; esto sucede frecuentemente en la solución numérica de ecuaciones diferenciales parciales .

Este método es económico y funciona bien en muchas aplicaciones; sin embargo, falla para valores de N cercanos a singular . Esto se ilustra mejor en el caso patológico de un cuadrado.${\displaystyle \mathbf {A} }$donde el determinante de N es el cuadrado del del sistema original Ax = l . Entonces es mejor aplicar la descomposición SVD o QR. La descomposición QR tiene la ventaja de que, de forma similar a las ecuaciones normales, no es necesario conservar toda la matriz A, ya que es posible actualizar el factor de Cholesky con filas consecutivas de A.

Los mínimos cuadrados no lineales son un caso particular de optimización no lineal.${\textstyle \mathbf {f} (\mathbf {x} )=\mathbf {l} }$ser un sistema de ecuaciones sobredeterminado con una función no lineal${\displaystyle \mathbf {f} }$Devuelve resultados vectoriales. El objetivo es minimizar la norma cuadrática de los residuos.${\textstyle \mathbf {v} =\mathbf {f} (\mathbf {x} )-\mathbf {l} }$. Se obtiene una solución aproximada del método de Newton mediante la expansión${\displaystyle \mathbf {f} }$en la serie Taylor reducida${\displaystyle {\bf {f(x_{\rm {0}}+\delta x)\approx f(x_{\rm {0}})+(\partial f/\partial x)\delta x}}}$problema de mínimos cuadrados lineales que produce${\displaystyle {\bf {\delta x}}}$

${\displaystyle {\bf {(\partial f/\partial x)\delta x=l-f(x_{\rm {0}})=v,\;\;\min _{\delta x}=\|v\|^{2}}}.}$

Por supuesto, debido a la omisión de términos de Taylor superiores, dicha solución es solo aproximada, si es que existe. Ahora se podría actualizar el punto de expansión a${\displaystyle {\bf {x_{\rm {n+1}}=x_{\rm {n}}+\delta x}}}$ y repetir todo el procedimiento, esperando que (i) las iteraciones converjan a una solución y (ii) que la solución sea la necesaria. Desafortunadamente, ninguna de las dos cosas está garantizada y deben verificarse.

Los mínimos cuadrados no lineales también pueden aplicarse al problema de mínimos cuadrados lineales estableciendo ${\displaystyle {\bf {x_{\rm {0}}=0}}}$y${\displaystyle {\bf {f(x_{\rm {0}})=Ax}}}$Esto puede resultar útil si la descomposición de Cholesky produce una inversa inexacta.${\displaystyle {\bf {R^{\rm {-1}}}}}$para la matriz triangular donde ${\displaystyle {\bf {R^{\rm {T}}R=N}}}$Debido a errores de redondeo, este procedimiento se denomina corrección diferencial de la solución. Siempre que las iteraciones converjan, en virtud del teorema del punto fijo de Banach, proporcionan la solución con una precisión limitada únicamente por la precisión de los residuos calculados.${\displaystyle {\bf {v=Ax-l}}}$. La precisión es independiente de los errores de redondeo en${\displaystyle {\bf {R^{\rm {-1}}}}}$. Pobre${\displaystyle {\bf {R^{\rm {-1}}}}}$puede restringir la región inicial${\displaystyle {\bf {x_{\rm {0}}}}}$ lo que produce convergencia o la impide por completo. Normalmente la convergencia es más lenta, por ejemplo lineal, de modo que${\displaystyle {\bf {\|\delta x_{\rm {n+1}}\|\approx \|=\alpha \delta x_{\rm {n}}\|}}}$donde constante ${\displaystyle \alpha <1}$Aitken puede acelerar dicha convergencia lenta.${\displaystyle \delta ^{2}}$método. Si el cálculo de${\displaystyle {\bf {R^{\rm {-1}}}}}$Es muy costoso, pero es posible usarlo desde iteraciones anteriores siempre que se mantenga la convergencia. Este procedimiento de Cholesky puede funcionar incluso para matrices de Hilbert, que son notoriamente difíciles de invertir. ^{[ 15 ]}

Las funciones multivariables no lineales pueden minimizarse sobre sus parámetros utilizando variantes del método de Newton llamadas métodos cuasi-Newton . En la iteración k, la búsqueda avanza en una dirección${\textstyle p_{k}}$definido por la resolución${\textstyle B_{k}p_{k}=-g_{k}}$para${\textstyle p_{k}}$, dónde${\textstyle p_{k}}$es la dirección del paso,${\textstyle g_{k}}$es el gradiente y${\textstyle B_{k}}$es una aproximación a la matriz hessiana formada mediante actualizaciones repetidas de rango 1 en cada iteración. Dos fórmulas de actualización bien conocidas son Davidon-Fletcher-Powell (DFP) y Broyden-Fletcher-Goldfarb-Shanno (BFGS). La pérdida de la condición definida positiva debido al error de redondeo se evita si, en lugar de actualizar una aproximación a la inversa de la hessiana, se actualiza la descomposición de Cholesky de una aproximación de la propia matriz hessiana. ^{[ 16 ]}

La descomposición de Cholesky se utiliza comúnmente en el método de Monte Carlo para simular sistemas con múltiples variables correlacionadas. La matriz de covarianza se descompone para obtener la matriz triangular inferior L. Al aplicar esto a un vector de observaciones no correlacionadas en una muestra u, se obtiene un vector muestral Lu con las propiedades de covarianza del sistema que se está modelando. ^{[ 17 ]}

El siguiente ejemplo simplificado muestra la economía que se obtiene de la descomposición de Cholesky: supongamos que el objetivo es generar dos variables normales correlacionadas.${\textstyle x_{1}}$y${\textstyle x_{2}}$con el coeficiente de correlación dado${\textstyle \rho }$Para lograrlo, primero es necesario generar dos variables aleatorias gaussianas no correlacionadas.${\textstyle z_{1}}$y${\textstyle z_{2}}$(por ejemplo, mediante una transformación de Box-Muller ). Dado el coeficiente de correlación requerido${\textstyle \rho }$, las variables normales correlacionadas se pueden obtener mediante las transformaciones${\textstyle x_{1}=z_{1}}$y${\textstyle x_{2}=\rho z_{1}+{\sqrt {1-\rho ^{2}}}z_{2}}$.

Los filtros de Kalman sin aroma suelen utilizar la descomposición de Cholesky para seleccionar un conjunto de puntos sigma. El filtro de Kalman registra el estado promedio de un sistema como un vector x de longitud N y la covarianza como una matriz P de N × N. La matriz P es siempre semidefinida positiva y puede descomponerse en L L ^T. Las columnas de L se pueden sumar y restar de la media x para formar un conjunto de 2 N vectores llamados puntos sigma . Estos puntos sigma capturan completamente la media y la covarianza del estado del sistema.

La inversa explícita de una matriz hermitiana se puede calcular mediante la descomposición de Cholesky, de manera similar a la resolución de sistemas lineales, utilizando${\textstyle n^{3}}$operaciones (${\textstyle {\tfrac {1}{2}}n^{3}}$multiplicaciones). ^{[ 11 ]} La inversión completa incluso se puede realizar de manera eficiente en el mismo lugar.

Una matriz no hermitiana B también puede invertirse utilizando la siguiente identidad, donde BB * siempre será hermitiana:

$${\displaystyle \mathbf {B} ^{-1}=\mathbf {B} ^{*}(\mathbf {BB} ^{*})^{-1}.}$$

La descomposición de Cholesky también puede utilizarse para imputar datos. Variaciones del algoritmo de maximización de la esperanza, entre otros algoritmos de imputación de datos, utilizan la descomposición de Cholesky. ^{[ 18 ]}

El algoritmo de Cholesky , utilizado para calcular la matriz de descomposición L , es una versión modificada de la eliminación gaussiana .

El algoritmo recursivo comienza con i  := 1 y

A ⁽¹⁾  := A .

En el paso i , la matriz A ^{( i )} tiene la siguiente forma: $${\displaystyle \mathbf {A} ^{(i)}={\begin{pmatrix}\mathbf {I} _{i-1}&0&0\\0&a_{i,i}&\mathbf {b} _{i}^{*}\\0&\mathbf {b} _{i}&\mathbf {B} ^{(i)}\end{pmatrix}},}$$ donde I _{i −1} denota la matriz identidad de dimensión i − 1 .

Si la matriz L _i se define por $${\displaystyle \mathbf {L} _{i}:={\begin{pmatrix}\mathbf {I} _{i-1}&0&0\\0&{\sqrt {a_{i,i}}}&0\\0&{\frac {1}{\sqrt {a_{i,i}}}}\mathbf {b} _{i}&\mathbf {I} _{n-i}\end{pmatrix}},}$$ (nótese que a _i,i > 0 ya que A ^{( i )} es definida positiva), entonces A ^{( i )} se puede escribir como $${\displaystyle \mathbf {A} ^{(i)}=\mathbf {L} _{i}\mathbf {A} ^{(i+1)}\mathbf {L} _{i}^{*}}$$ dónde $${\displaystyle \mathbf {A} ^{(i+1)}={\begin{pmatrix}\mathbf {I} _{i-1}&0&0\\0&1&0\\0&0&\mathbf {B} ^{(i)}-{\frac {1}{a_{i,i}}}\mathbf {b} _{i}\mathbf {b} _{i}^{*}\end{pmatrix}}.}$$ Tenga en cuenta que b _i b _i * es un producto exterior , por lo tanto, este algoritmo se denomina versión de producto exterior en (Golub & Van Loan).

Esto se repite para i desde 1 hasta n . Después de n pasos, se obtiene A ^{( n +1)} = I , y por lo tanto, la matriz triangular inferior L buscada se calcula como

$${\displaystyle \mathbf {L}$$ :=\mathbf {L} _{1}\mathbf {L} _{2}\dots \mathbf {L} _{n}.}

Si la ecuación $${\displaystyle {\begin{aligned}\mathbf {A} =\mathbf {LL} ^{T}&={\begin{pmatrix}L_{11}&0&0\\L_{21}&L_{22}&0\\L_{31}&L_{32}&L_{33}\\\end{pmatrix}}{\begin{pmatrix}L_{11}&L_{21}&L_{31}\\0&L_{22}&L_{32}\\0&0&L_{33}\end{pmatrix}}\\[8pt]&={\begin{pmatrix}L_{11}^{2}&&({\text{symmetric}})\\L_{21}L_{11}&L_{21}^{2}+L_{22}^{2}&\\L_{31}L_{11}&L_{31}L_{21}+L_{32}L_{22}&L_{31}^{2}+L_{32}^{2}+L_{33}^{2}\end{pmatrix}},\end{aligned}}}$$

Se escribe y se obtiene lo siguiente:

$${\displaystyle {\begin{aligned}\mathbf {L} ={\begin{pmatrix}{\sqrt {A_{11}}}&0&0\\A_{21}/L_{11}&{\sqrt {A_{22}-L_{21}^{2}}}&0\\A_{31}/L_{11}&\left(A_{32}-L_{31}L_{21}\right)/L_{22}&{\sqrt {A_{33}-L_{31}^{2}-L_{32}^{2}}}\end{pmatrix}}\end{aligned}}}$$

y por lo tanto las siguientes fórmulas para las entradas de L :

$${\displaystyle L_{j,j}=(\pm ){\sqrt {A_{j,j}-\sum _{k=1}^{j-1}L_{j,k}^{2}}},}$$$${\displaystyle L_{i,j}={\frac {1}{L_{j,j}}}\left(A_{i,j}-\sum _{k=1}^{j-1}L_{i,k}L_{j,k}\right)\quad {\text{for }}i>j.}$$

Para matrices complejas y reales, se permiten cambios de signo arbitrarios e intrascendentes en los elementos diagonales y sus elementos no diagonales asociados. La expresión bajo la raíz cuadrada siempre es positiva si A es real y definida positiva.

Para matrices hermitianas complejas, se aplica la siguiente fórmula:

$${\displaystyle L_{j,j}={\sqrt {A_{j,j}-\sum _{k=1}^{j-1}L_{j,k}^{*}L_{j,k}}},}$$$${\displaystyle L_{i,j}={\frac {1}{L_{j,j}}}\left(A_{i,j}-\sum _{k=1}^{j-1}L_{j,k}^{*}L_{i,k}\right)\quad {\text{for }}i>j.}$$

y se puede demostrar que${\displaystyle L_{j,j}}$siempre es real y positivo si A es definida positiva. ^{[ 21 ] : 49}

Así pues, ahora es posible calcular la entrada ( i , j ) si se conocen las entradas de la izquierda y de arriba. El cálculo suele organizarse en alguno de los siguientes órdenes:

• El algoritmo de Cholesky-Banachiewicz comienza desde la esquina superior izquierda de la matriz L y procede a calcular la matriz fila por fila.

para ( i = 0 ; i < dimensionSize ; i ++ ) { para ( j = 0 ; j <= i ; j ++ ) { float sum = 0 ; para ( 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 )); } }

El algoritmo anterior se puede expresar sucintamente como la combinación de un producto escalar y una multiplicación de matrices en lenguajes de programación vectorizados como Fortran , de la siguiente manera:

hacer i = 1 , tamaño ( A , 1 ) L ( i , i ) = sqrt ( A ( i , i ) - producto_punto ( L ( i , 1 : i - 1 ), L ( i , 1 : i - 1 ))) L ( i + 1 :, i ) = ( A ( i + 1 :, i ) - matmul ( conjg ( L ( i , 1 : i - 1 )), L ( i + 1 :, 1 : i - 1 ))) / L ( i , i ) fin hacer

donde conjgse refiere al conjugado complejo de los elementos.

• El algoritmo de Cholesky-Crout comienza en la esquina superior izquierda de la matriz L y procede a calcular la matriz columna por columna.

  para ( j = 0 ; j < dimensionSize ; j ++ ) { float sum = 0 ; para ( k = 0 ; k < j ; k ++ ) { sum += L [ j ][ k ] * L [ j ][ k ]; } L [ j ][ j ] = sqrt ( A [ j ][ j ] - sum );para ( i = j + 1 ; i < dimensionSize ; i ++ ) { suma = 0 ; para ( k = 0 ; k < j ; k ++ ) { suma += L [ i ][ k ] * L [ j ][ k ]; } L [ i ][ j ] = ( 1.0 / L [ j ][ j ] * ( A [ i ][ j ] - suma )); } }

El algoritmo anterior se puede expresar sucintamente como la combinación de un producto escalar y una multiplicación de matrices en lenguajes de programación vectorizados como Fortran , de la siguiente manera:

hacer i = 1 , tamaño ( A , 1 ) L ( i , i ) = sqrt ( A ( i , i ) - producto_punto ( L ( 1 : i - 1 , i ), L ( 1 : i - 1 , i ))) L ( i , i + 1 :) = ( A ( i , i + 1 :) - matmul ( conjg ( L ( 1 : i - 1 , i )), L ( 1 : i - 1 , i + 1 :))) / L ( i , i ) fin hacer

donde conjgse refiere al conjugado complejo de los elementos.

Cualquiera de los dos patrones de acceso permite realizar todo el cálculo in situ si se desea.

Supongamos que se desea resolver un sistema de ecuaciones lineales bien condicionado . Si se utiliza la descomposición LU, el algoritmo resulta inestable a menos que se emplee alguna estrategia de pivoteo. En este último caso, el error depende del denominado factor de crecimiento de la matriz, que suele ser pequeño (aunque no siempre).

Ahora bien, supongamos que la descomposición de Cholesky es aplicable. Como se mencionó anteriormente, el algoritmo será el doble de rápido. Además, no es necesario pivotar y el error siempre será pequeño. Específicamente, si Ax = b y y denota la solución calculada, entonces y resuelve el sistema perturbado ( A + E ) y = b , donde $${\displaystyle \|\mathbf {E} \|_{2}\leq c_{n}\varepsilon \|\mathbf {A} \|_{2}.}$$ Aquí ||·|| ₂ es la norma 2 de la matriz , c _n es una pequeña constante que depende de n , y ε denota el redondeo unitario .

Una preocupación con la descomposición de Cholesky que debe tenerse en cuenta es el uso de raíces cuadradas. Si la matriz que se está factorizando es definida positiva, como se requiere, los números bajo las raíces cuadradas siempre son positivos en aritmética exacta . Desafortunadamente, los números pueden volverse negativos debido a errores de redondeo , en cuyo caso el algoritmo no puede continuar. Sin embargo, esto solo puede ocurrir si la matriz está muy mal condicionada. Una forma de abordar esto es agregar una matriz de corrección diagonal a la matriz que se está descomponiendo para intentar promover la definición positiva. ^{[ 22 ]} Si bien esto podría disminuir la precisión de la descomposición, puede ser muy favorable por otras razones; por ejemplo, al realizar el método de Newton en optimización , agregar una matriz diagonal puede mejorar la estabilidad cuando está lejos del óptimo.

Una forma alternativa, que elimina la necesidad de tomar raíces cuadradas cuando A es simétrica, es la factorización indefinida simétrica ^{[ 21 ] : 84}$${\displaystyle {\begin{aligned}\mathbf {A} =\mathbf {LDL} ^{\mathrm {T} }&={\begin{pmatrix}1&0&0\\L_{21}&1&0\\L_{31}&L_{32}&1\\\end{pmatrix}}{\begin{pmatrix}D_{1}&0&0\\0&D_{2}&0\\0&0&D_{3}\\\end{pmatrix}}{\begin{pmatrix}1&L_{21}&L_{31}\\0&1&L_{32}\\0&0&1\\\end{pmatrix}}\\[8pt]&={\begin{pmatrix}D_{1}&&(\mathrm {symmetric} )\\L_{21}D_{1}&L_{21}^{2}D_{1}+D_{2}&\\L_{31}D_{1}&L_{31}L_{21}D_{1}+L_{32}D_{2}&L_{31}^{2}D_{1}+L_{32}^{2}D_{2}+D_{3}.\end{pmatrix}}.\end{aligned}}}$$

Las siguientes relaciones recursivas se aplican a las entradas de D y L : $${\displaystyle D_{j}=A_{jj}-\sum _{k=1}^{j-1}L_{jk}^{2}D_{k},}$$$${\displaystyle L_{ij}={\frac {1}{D_{j}}}\left(A_{ij}-\sum _{k=1}^{j-1}L_{ik}L_{jk}D_{k}\right)\quad {\text{for }}i>j.}$$

Esto funciona siempre que los elementos diagonales generados en D no sean cero. En ese caso, la descomposición es única. D y L son reales si A es real.

Para la matriz hermitiana compleja A , se aplica la siguiente fórmula:

$${\displaystyle D_{j}=A_{jj}-\sum _{k=1}^{j-1}L_{jk}L_{jk}^{*}D_{k},}$$$${\displaystyle L_{ij}={\frac {1}{D_{j}}}\left(A_{ij}-\sum _{k=1}^{j-1}L_{ik}L_{jk}^{*}D_{k}\right)\quad {\text{for }}i>j.}$$

Nuevamente, el patrón de acceso permite que, si se desea, todo el cálculo se realice in situ.

Cuando se utiliza en matrices indefinidas, se sabe que la factorización LDL * es inestable sin un pivoteo cuidadoso; ^{[ 23 ]} específicamente, los elementos de la factorización pueden crecer arbitrariamente. Una posible mejora es realizar la factorización en submatrices de bloques, comúnmente de 2 × 2: ^{[ 24 ]}

$${\displaystyle {\begin{aligned}\mathbf {A} =\mathbf {LDL} ^{\mathrm {T} }&={\begin{pmatrix}\mathbf {I} &0&0\\\mathbf {L} _{21}&\mathbf {I} &0\\\mathbf {L} _{31}&\mathbf {L} _{32}&\mathbf {I} \\\end{pmatrix}}{\begin{pmatrix}\mathbf {D} _{1}&0&0\\0&\mathbf {D} _{2}&0\\0&0&\mathbf {D} _{3}\\\end{pmatrix}}{\begin{pmatrix}\mathbf {I} &\mathbf {L} _{21}^{\mathrm {T} }&\mathbf {L} _{31}^{\mathrm {T} }\\0&\mathbf {I} &\mathbf {L} _{32}^{\mathrm {T} }\\0&0&\mathbf {I} \\\end{pmatrix}}\\[8pt]&={\begin{pmatrix}\mathbf {D} _{1}&&(\mathrm {symmetric} )\\\mathbf {L} _{21}\mathbf {D} _{1}&\mathbf {L} _{21}\mathbf {D} _{1}\mathbf {L} _{21}^{\mathrm {T} }+\mathbf {D} _{2}&\\\mathbf {L} _{31}\mathbf {D} _{1}&\mathbf {L} _{31}\mathbf {D} _{1}\mathbf {L} _{21}^{\mathrm {T} }+\mathbf {L} _{32}\mathbf {D} _{2}&\mathbf {L} _{31}\mathbf {D} _{1}\mathbf {L} _{31}^{\mathrm {T} }+\mathbf {L} _{32}\mathbf {D} _{2}\mathbf {L} _{32}^{\mathrm {T} }+\mathbf {D} _{3}\end{pmatrix}},\end{aligned}}}$$

donde cada elemento de las matrices anteriores es una submatriz cuadrada. A partir de esto, se derivan las siguientes relaciones recursivas análogas:

$${\displaystyle \mathbf {D} _{j}=\mathbf {A} _{jj}-\sum _{k=1}^{j-1}\mathbf {L} _{jk}\mathbf {D} _{k}\mathbf {L} _{jk}^{\mathrm {T} },}$$$${\displaystyle \mathbf {L} _{ij}=\left(\mathbf {A} _{ij}-\sum _{k=1}^{j-1}\mathbf {L} _{ik}\mathbf {D} _{k}\mathbf {L} _{jk}^{\mathrm {T} }\right)\mathbf {D} _{j}^{-1}.}$$

Esto implica productos matriciales e inversión explícita, lo que limita el tamaño práctico del bloque.

Una tarea que surge con frecuencia en la práctica es la necesidad de actualizar una descomposición de Cholesky. En más detalle, ya se ha calculado la descomposición de Cholesky.${\textstyle \mathbf {A} =\mathbf {L} \mathbf {L} ^{*}}$de alguna matriz${\textstyle \mathbf {A} }$, entonces se cambia la matriz${\textstyle \mathbf {A} }$de alguna manera en otra matriz, por ejemplo${\textstyle {\tilde {\mathbf {A} }}}$y se desea calcular la descomposición de Cholesky de la matriz actualizada:${\textstyle {\tilde {\mathbf {A} }}={\tilde {\mathbf {L} }}{\tilde {\mathbf {L} }}^{*}}$. La pregunta ahora es si se puede utilizar la descomposición de Cholesky de${\textstyle \mathbf {A} }$que se calculó previamente para calcular la descomposición de Cholesky de${\textstyle {\tilde {\mathbf {A} }}}$.

function L = updateChol ( L,x,c ) % dada la descomposición de Cholesky L*L' de una matriz, calcula el factor actualizado % L para que tengamos la descomposición de Cholesky de L*L'+c*x*x'; n = length ( x ); for k = 1 : n - 1 l = L (:, k ); % valor antiguo de la k-ésima columna lk = l ( k ); xk = x ( k ); dk = sqrt ( lk ^ 2 + c * xk ^ 2 ); % nuevo valor diagonal L (:, k )=( lk / dk ) * l + ( c * xk / dk ) * x ; % nuevo valor de columna x = x - l * ( xk / lk ); c = c * ( lk / dk ) ^ 2 ; fin L ( n , n )= sqrt ( L ( n , n ) ^ 2 + c * x ( n ) ^ 2 ); fin

Los algoritmos anteriores muestran que toda matriz definida positiva${\textstyle \mathbf {A} }$posee una descomposición de Cholesky. Este resultado puede extenderse al caso semidefinido positivo mediante un argumento límite. El argumento no es completamente constructivo, es decir, no proporciona algoritmos numéricos explícitos para calcular los factores de Cholesky.

Si${\textstyle \mathbf {A} }$es un${\textstyle n\times n}$matriz semidefinida positiva , entonces la secuencia${\textstyle \left(\mathbf {A} _{k}\right)_{k}:=\left(\mathbf {A} +{\frac {1}{k}}\mathbf {I} _{n}\right)_{k}}$consta de matrices definidas positivas . (Esto es una consecuencia inmediata de, por ejemplo, el teorema de mapeo espectral para el cálculo funcional polinomial). Además, $${\displaystyle \mathbf {A} _{k}\rightarrow \mathbf {A} \quad {\text{for}}\quad k\rightarrow \infty }$$ en la norma del operador . Del caso definido positivo, cada${\textstyle \mathbf {A} _{k}}$tiene descomposición de Cholesky${\textstyle \mathbf {A} _{k}=\mathbf {L} _{k}\mathbf {L} _{k}^{*}}$. Por propiedad de la norma del operador,

$${\displaystyle \|\mathbf {L} _{k}\|^{2}\leq \|\mathbf {L} _{k}\mathbf {L} _{k}^{*}\|=\|\mathbf {A} _{k}\|\,.}$$

El${\textstyle \leq }$se sostiene porque${\textstyle M_{n}(\mathbb {C} )}$equipado con la norma del operador es un álgebra C*. Por lo tanto,${\textstyle \left(\mathbf {L} _{k}\right)_{k}}$es un conjunto acotado en el espacio de Banach de operadores, por lo tanto relativamente compacto (porque el espacio vectorial subyacente es de dimensión finita). En consecuencia, tiene una subsucesión convergente, también denotada por${\textstyle \left(\mathbf {L} _{k}\right)_{k}}$, con límite${\textstyle \mathbf {L} }$Se puede comprobar fácilmente que esto${\textstyle \mathbf {L} }$tiene las propiedades deseadas, es decir${\textstyle \mathbf {A} =\mathbf {L} \mathbf {L} ^{*}}$, y${\textstyle \mathbf {L} }$es triangular inferior con entradas diagonales no negativas: para todo${\textstyle x}$y${\textstyle y}$,

$${\displaystyle \langle \mathbf {A} x,y\rangle =\left\langle \lim \mathbf {A} _{k}x,y\right\rangle =\langle \lim \mathbf {L} _{k}\mathbf {L} _{k}^{*}x,y\rangle =\langle \mathbf {L} \mathbf {L} ^{*}x,y\rangle \,.}$$

Por lo tanto,${\textstyle \mathbf {A} =\mathbf {L} \mathbf {L} ^{*}}$Debido a que el espacio vectorial subyacente es de dimensión finita, todas las topologías en el espacio de operadores son equivalentes. Por lo tanto,${\textstyle \left(\mathbf {L} _{k}\right)_{k}}$tiende a${\textstyle \mathbf {L} }$en norma significa${\textstyle \left(\mathbf {L} _{k}\right)_{k}}$tiende a${\textstyle \mathbf {L} }$entrada por entrada. Esto a su vez implica que, dado que cada${\textstyle \mathbf {L} _{k}}$es triangular inferior con entradas diagonales no negativas,${\textstyle \mathbf {L} }$También lo es.

Dejar${\textstyle \mathbf {A} }$Sea una matriz hermitiana semidefinida positiva . Entonces se puede escribir como un producto de su matriz de raíz cuadrada ,${\textstyle \mathbf {A} =\mathbf {B} \mathbf {B} ^{*}}$Ahora se puede aplicar la descomposición QR a${\textstyle \mathbf {B} ^{*}}$, Resultando en${\textstyle \mathbf {B} ^{*}=\mathbf {Q} \mathbf {R} }$ , dónde${\textstyle \mathbf {Q} }$es unitario y${\textstyle \mathbf {R} }$es triangular superior. Al insertar la descomposición en la igualdad original se obtiene${\textstyle A=\mathbf {B} \mathbf {B} ^{*}=(\mathbf {QR} )^{*}\mathbf {QR} =\mathbf {R} ^{*}\mathbf {Q} ^{*}\mathbf {QR} =\mathbf {R} ^{*}\mathbf {R} }$. Configuración${\textstyle \mathbf {L} =\mathbf {R} ^{*}}$Completa la demostración.