Descomposición de Cholesky

format_list_bulleted Contenido keyboard_arrow_down

ImprimirCitar

Método de descomposición de matriz

En álgebra lineal, la descomposición de Cholesky o factorización de Cholesky (pronunciado shə-LES-kee) es una descomposición de una matriz definida positiva hermitiana en el producto de una matriz triangular inferior y su transpuesta conjugada, que es útil para soluciones numéricas eficientes, por ejemplo, simulaciones de Monte Carlo. Fue descubierto por André-Louis Cholesky para matrices reales y publicado póstumamente en 1924. Cuando es aplicable, la descomposición de Cholesky es aproximadamente el doble de eficiente que la descomposición LU para resolver sistemas de ecuaciones lineales.

Declaración

La descomposición de Cholesky de una matriz definida positiva hermítica A, es una descomposición de la forma

{displaystyle mathbf {A} =mathbf {LL} ^{*},}

donde L es una matriz triangular inferior con entradas diagonales reales y positivas, y L* denota la transpuesta conjugada de L. Cada matriz definida positiva hermitiana (y, por lo tanto, también cada matriz definida positiva simétrica de valor real) tiene una descomposición de Cholesky única.

Lo contrario es trivial: si A se puede escribir como LL* para alguna L invertible, triangular inferior o de otro tipo, entonces A es hermítica y definida positiva.

Cuando A es una matriz real (por lo tanto simétrica definida positiva), la factorización se puede escribir

{displaystyle mathbf {A} =mathbf {LL} ^{mathsf {T}},}

donde L es una matriz triangular inferior real con entradas diagonales positivas.

Matrices semidefinidas positivas

Si una matriz hermitiana A es solo semidefinida positiva, en lugar de definida positiva, entonces todavía tiene una descomposición de la forma A = LL* donde las entradas diagonales de L pueden ser cero. La descomposición no necesita ser única, por ejemplo:

{displaystyle {begin{bmatrix}0&0\0&1end{bmatrix}}=mathbf {L} mathbf {L} ^{*},quad quad mathbf {L} ={begin{bmatrix}0&0\cos theta &sin theta end{bmatrix}}.}

Sin embargo, si el rango de A es r, entonces hay un L triangular inferior único con exactamente r elementos diagonales positivos y columnas n−r que contienen todos ceros.

Alternativamente, la descomposición se puede hacer única cuando se fija una opción pivotante. Formalmente, si A es un n × n matriz semidefinida positiva de rango r, entonces hay al menos una matriz de permutación P tales que P A P^T tiene una descomposición única de la forma P A P^T = L^* con ${displaystyle mathbf {L} ={begin{bmatrix}mathbf {L} _{1}&0\mathbf {L} _{2}&0end{bmatrix}}}$ , Donde L₁ es un r × r matriz triangular inferior con diagonal positivo.

Descomposición de LDL

Una variante estrechamente relacionada de la descomposición clásica de Cholesky es la descomposición de LDL,

{displaystyle mathbf {A} =mathbf {LDL} ^{*},}

donde L es una matriz triangular (unitriangular) de unidad inferior y D es una matriz diagonal. Es decir, se requiere que los elementos diagonales de L sean 1 a costa de introducir una matriz diagonal adicional D en la descomposición. La principal ventaja es que la descomposición de LDL se puede calcular y utilizar esencialmente con los mismos algoritmos, pero evita la extracción de raíces cuadradas.

Por esta razón, la descomposición de LDL a menudo se denomina descomposición de Cholesky sin raíces cuadradas. Para matrices reales, la factorización tiene la forma A = LDL^T y a menudo se denomina como descomposición de LDLT (o descomposición de LDL^T, o LDL′). Es una reminiscencia de la descomposición propia de matrices simétricas reales, A = QΛQ^T, pero es bastante diferente en la práctica porque Λ y D no son matrices similares.

La descomposición LDL está relacionada con la descomposición clásica de Cholesky de la forma LL* de la siguiente manera:

{displaystyle mathbf {A} =mathbf {LDL} ^{*}=mathbf {L} mathbf {D} ^{1/2}left(mathbf {D} ^{1/2}right)^{*}mathbf {L} ^{*}=mathbf {L} mathbf {D} ^{1/2}left(mathbf {L} mathbf {D} ^{1/2}right)^{*}.}

Por el contrario, dada la descomposición clásica de Cholesky ${displaystyle mathbf {A} =mathbf {C} mathbf {C} ^{*}}$ de una matriz definida positiva, si S es una matriz diagonal que contiene la diagonal principal $mathbf C$ , entonces un A puede ser descompuesto ${displaystyle mathbf {L} mathbf {D} mathbf {L} ^{*}}$ Donde

{displaystyle mathbf {L} =mathbf {C} mathbf {S} ^{-1}}

(esto reescala cada columna para hacer elementos diagonales 1),

{displaystyle mathbf {D} =mathbf {S} ^{2}.}

Si A es positivo definido entonces los elementos diagonales de D son todos positivos. Para semidefinido positivo A, un ${displaystyle mathbf {L} mathbf {D} mathbf {L} ^{*}}$ descomposición existe donde el número de elementos no cero en la diagonal D es exactamente el rango de A. Algunas matrices indefinidas para las cuales no existe descomposición de Cholesky tienen una descomposición LDL con entradas negativas D: basta que el primero n−1 principales menores de edad A son no fijos.

Ejemplo

Aquí está la descomposición de Cholesky de una matriz real simétrica:

{displaystyle {begin{aligned}left({begin{array}{*{3}{r}}4&12&-16\12&37&-43\-16&-43&98\end{array}}right)=left({begin{array}{*{3}{r}}2&0&0\6&1&0\-8&5&3\end{array}}right)left({begin{array}{*{3}{r}}2&6&-8\0&1&5\0&0&3\end{array}}right).end{aligned}}}

Y aquí está su descomposición de LDL^T:

{displaystyle {begin{aligned}left({begin{array}{*{3}{r}}4&12&-16\12&37&-43\-16&-43&98\end{array}}right)&=left({begin{array}{*{3}{r}}1&0&0\3&1&0\-4&5&1\end{array}}right)left({begin{array}{*{3}{r}}4&0&0\0&1&0\0&0&9\end{array}}right)left({begin{array}{*{3}{r}}1&3&-4\0&1&5\0&0&1\end{array}}right).end{aligned}}}

Aplicaciones

La descomposición de Cholesky se utiliza principalmente para la solución numérica de ecuaciones lineales $mathbf {Ax} =mathbf {b}$ . Si A es simétrico y positivo definido, entonces podemos resolver $mathbf {Ax} =mathbf {b}$ por primera vez computar la descomposición de Cholesky ${displaystyle mathbf {A} =mathbf {LL} ^{mathrm {*} }}$ , luego resolver ${displaystyle mathbf {Ly} =mathbf {b} }$ para Sí. por sustitución anticipada, y finalmente solución ${displaystyle mathbf {L^{*}x} =mathbf {y} }$ para x por sustitución trasera.

Una manera alternativa de eliminar la toma de raíces cuadradas ${displaystyle mathbf {LL} ^{mathrm {*} }}$ descomposición es computar la descomposición LDL ${displaystyle mathbf {A} =mathbf {LDL} ^{mathrm {*} }}$ , luego resolver ${displaystyle mathbf {Ly} =mathbf {b} }$ para Sí., y finalmente resolver ${displaystyle mathbf {DL} ^{mathrm {*} }mathbf {x} =mathbf {y} }$ .

Para sistemas lineales que se pueden poner en forma simétrica, la descomposición de Cholesky (o su variante LDL) es el método de elección, para una eficiencia superior y estabilidad numérica. En comparación con la descomposición LU, es aproximadamente el doble de eficiente.

Mínimos cuadrados lineales

Los sistemas de la forma Ax = b con A simétrica y definida positiva surgen con bastante frecuencia en las aplicaciones. Por ejemplo, las ecuaciones normales en problemas de mínimos cuadrados lineales tienen esta forma. También puede ocurrir que la matriz A provenga de un funcional energético, que debe ser positivo por consideraciones físicas; esto sucede con frecuencia en la solución numérica de ecuaciones diferenciales parciales.

Optimización no lineal

Las funciones multivariadas no lineales pueden minimizarse sobre sus parámetros utilizando variantes del método de Newton llamado quasi-Newton métodos. En la iteración k, los pasos de búsqueda en una dirección ${displaystyle p_{k}}$ definida por la resolución ${displaystyle B_{k}p_{k}=-g_{k}}$ para ${displaystyle p_{k}}$ , donde ${displaystyle p_{k}}$ es la dirección paso, ${displaystyle g_{k}}$ es el gradiente, y $B_{k}$ es una aproximación a la matriz hesiana formada por la repetición de actualizaciones de rango-1 en cada iteración. Dos fórmulas de actualización bien conocidas se llaman Davidon–Fletcher–Powell (DFP) y Broyden–Fletcher–Goldfarb–Shanno (BFGS). La pérdida de la condición positiva-definida a través del error redondeado se evita si en lugar de actualizar una aproximación a la inversa del Hessian, se actualiza la descomposición Cholesky de una aproximación de la matriz hesiana misma .

Simulación de Montecarlo

La descomposición de Cholesky se usa comúnmente en el método Monte Carlo para simular sistemas con múltiples variables correlacionadas. La matriz de covarianza se descompone para dar el L triangular inferior. Aplicar esto a un vector de muestras no correlacionadas u produce un vector de muestra Lu con las propiedades de covarianza del sistema que se está modelando.

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 correlativas $x_{1}$ y $x_{2}$ con coeficiente de correlación dado $rho$ . Para lograrlo, es necesario primero generar dos variables aleatorias de Gauss $z_{1}$ y $z_{2}$ , que se puede hacer utilizando una transformación Box-Muller. Dado el coeficiente de correlación requerido $rho$ , las variables normales correlativas se pueden obtener a través de las transformaciones $x_{1}=z_{1}$ y ${textstyle x_{2}=rho z_{1}+{sqrt {1-rho ^{2}}}z_{2}}$ .

Filtros Kalman

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

Inversión de matriz

La inversa explícita de una matriz hermitiana puede ser computada por la descomposición de Cholesky, de una manera similar a resolver sistemas lineales, utilizando $n^{3}$ operaciones ${displaystyle {tfrac {1}{2}}n^{3}}$ multiplicaciones). Toda la inversión puede incluso ser realizada eficientemente en el lugar.

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

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

Cálculo

Existen varios métodos para calcular la descomposición de Cholesky. La complejidad computacional de los algoritmos de uso común es O(n³) en general. Los algoritmos que se describen a continuación involucran aproximadamente (1/3)n³ FLOP (n³/6 multiplicaciones y el mismo número de adiciones) para sabores reales y (4/3)n³ FLOP para sabores complejos, donde n es el tamaño del matriz A. Por lo tanto, tienen la mitad del costo de la descomposición LU, que usa 2n³/3 FLOP (ver Trefethen y Bau 1997).

Cuál de los siguientes algoritmos es más rápido depende de los detalles de la implementación. Generalmente, el primer algoritmo será un poco más lento porque accede a los datos de una manera menos regular.

El algoritmo de Cholesky

El algoritmo 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⁽ⁱ⁾ tiene la siguiente forma:

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 ahora definimos la matriz L_i por

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}},

(tenga en cuenta que a_i,i > 0 desde A^{( i)} es positivo definido), entonces podemos escribir A⁽ⁱ⁾ como

mathbf {A} ^{(i)}=mathbf {L} _{i}mathbf {A} ^{(i+1)}mathbf {L} _{i}^{*}

dónde

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 externo, por lo que este algoritmo se denomina versión del producto externo en (Golub & Van Loan).

Repetimos esto para i del 1 al n. Después de n pasos, obtenemos A⁽ⁿ⁺¹⁾ = I. Por lo tanto, la matriz triangular inferior L que estamos buscando se calcula como

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

Algoritmos Cholesky-Banachiewicz y Cholesky-Crout

Patrón de acceso (blanco) y patrón de escritura (amarillo) para el algoritmo en el lugar Cholesky—Banachiewicz en una matriz 5×5

Si escribimos 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}}}

obtenemos 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}}},}

j.}" xmlns="http://www.w3.org/1998/Math/MathML">Li,j=1Lj,j()Ai,j− − .. k=1j− − 1Li,kLj,k)parai■j.{displaystyle L_{i,j}={frac {1}{L_{j,j}}left(A_{i,j}-sum ¿Por qué? {texto {fnMicrosoft Sans Serif}j.}" aria-hidden="true" class="mwe-math-fallback-image-inline" src="https://wikimedia.org/api/rest_v1/media/math/render/svg/89f8de93310247c939102367ecc284d1f4beff9d" style="vertical-align: -3.171ex; width:44.493ex; height:7.509ex;"/>

Para matrices complejas y reales, se permiten cambios de signo arbitrarios intrascendentes de elementos diagonales y asociados fuera de la diagonal. La expresión debajo de la raíz cuadrada siempre es positiva si A es real y definida positiva.

Para matriz hermítica compleja, 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}^{*}}},}

j.}" xmlns="http://www.w3.org/1998/Math/MathML">Li,j=1Lj,j()Ai,j− − .. k=1j− − 1Li,kLj,kAlternativa Alternativa )parai■j.{displaystyle L_{i,j}={frac {1}{L_{j,j}}left(A_{i,j}-sum ¿Por qué? {texto {fnMicrosoft Sans Serif}j.}" aria-hidden="true" class="mwe-math-fallback-image-inline" src="https://wikimedia.org/api/rest_v1/media/math/render/svg/bb5adc5916e0762b2eca921de3e70ccae9bd2999" style="vertical-align: -3.171ex; width:44.493ex; height:7.509ex;"/>

Entonces podemos calcular la entrada (i, j) si conocemos las entradas a la izquierda y arriba. El cálculo generalmente se organiza en cualquiera de los siguientes órdenes:

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

para ()i = 0; i . dimension; i++) {} para ()j = 0; j . i; j++) {} flotador suma = 0; para ()k = 0; k . j; k++) suma += L[i[ ]k] * L[j[ ]k]; si ()i == j) L[i[ ]j] = Sqrt()A[i[ ]i] - suma); más L[i[ ]j] = ()1.0 / L[j[ ]j] * ()A[i[ ]j] - suma)); }}

El Algoritmo de Cholesky-Crout comienza desde la esquina superior izquierda de la matriz L y procede a calcular la columna matriz por columna.

para ()j = 0; j . dimension; j++) {} flotador suma = 0; para ()k = 0; k . j; k++) {} suma += L[j[ ]k] * L[j[ ]k]; } L[j[ ]j] = Sqrt()A[j[ ]j] - suma); para ()i = j + 1; i . dimension; 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)); }}

Cualquier patrón de acceso permite que todo el cálculo se realice en el lugar si se desea.

Estabilidad del cómputo

Supongamos que queremos resolver un sistema de ecuaciones lineales bien condicionado. Si se usa la descomposición LU, entonces el algoritmo es inestable a menos que usemos algún tipo de estrategia de pivote. En este último caso, el error depende del llamado factor de crecimiento de la matriz, que suele ser (pero no siempre) pequeño.

Ahora, suponga que se aplica la descomposición de Cholesky. 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 queremos resolver Ax = b, y y denota la solución calculada, entonces y resuelve el sistema perturbado (A + E)y = b, donde

|mathbf {E} |_{2}leq c_{n}varepsilon |mathbf {A} |_{2}.

Aquí ||·||₂ es la matriz 2-norma, c_n es una pequeña constante que depende de n , y ε denota el redondeo de unidades.

Una preocupación con la descomposición de Cholesky a tener en cuenta es el uso de raíces cuadradas. Si la matriz que se factoriza es definida positiva como se requiere, los números debajo de 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 suceder si la matriz está muy mal acondicionada. Una forma de abordar esto es agregar una matriz de corrección diagonal a la matriz que se descompone en un intento de promover la definición positiva. 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 la optimización, agregar una matriz diagonal puede mejorar la estabilidad cuando está lejos del óptimo.

Descomposición de LDL

Una forma alternativa, que elimina la necesidad de sacar raíces cuadradas cuando A es simétrica, es la factorización indefinida simétrica

{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},}

j.}" xmlns="http://www.w3.org/1998/Math/MathML">Lij=1Dj()Aij− − .. k=1j− − 1LikLjkDk)parai■j.{displaystyle L_{ij}={frac {1}{d_{j}}left(A_{ij}-sum ¿Por qué? {texto {fnMicrosoft Sans Serif}j.}" aria-hidden="true" class="mwe-math-fallback-image-inline" src="https://wikimedia.org/api/rest_v1/media/math/render/svg/1335e3e24c2341fe4138935f261bf2751eb7bdb8" style="vertical-align: -3.171ex; width:44.883ex; height:7.509ex;"/>

Esto funciona siempre que los elementos diagonales generados en D no sean cero. La descomposición es entonces ú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},}

j.}" xmlns="http://www.w3.org/1998/Math/MathML">Lij=1Dj()Aij− − .. k=1j− − 1LikLjkAlternativa Alternativa Dk)parai■j.{displaystyle L_{ij}={frac {1}{d_{j}}left(A_{ij}-sum ¿Por qué? {texto {fnMicrosoft Sans Serif}j.}" aria-hidden="true" class="mwe-math-fallback-image-inline" src="https://wikimedia.org/api/rest_v1/media/math/render/svg/88fb8d13a2eb2c7442a570544cf0b8ebb370dc70" style="vertical-align: -3.171ex; width:44.883ex; height:7.509ex;"/>

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

Variante de bloque

Cuando se usa en matrices indefinidas, se sabe que la factorización LDL* es inestable sin un pivote cuidadoso; específicamente, los elementos de la factorización pueden crecer arbitrariamente. Una posible mejora es realizar la factorización en submatrices de bloque, comúnmente 2 × 2:

{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. De aquí se siguen estas 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 de matriz e inversión explícita, lo que limita el tamaño práctico del bloque.

Actualizando la descomposición

Una tarea que a menudo surge en la práctica es que hay que actualizar una descomposición de Cholesky. En más detalles, uno ya ha computado la descomposición Cholesky ${displaystyle mathbf {A} =mathbf {L} mathbf {L} ^{*}}$ de alguna matriz $mathbf {A}$ , entonces uno cambia la matriz $mathbf {A}$ de alguna manera en otra matriz, digamos ${tilde {mathbf {A} }}$ , y uno quiere computar la descomposición Cholesky de la matriz actualizada: ${tilde {mathbf {A} }}={tilde {mathbf {L} }}{tilde {mathbf {L} }}^{*}$ . La pregunta es ahora si se puede utilizar la descomposición de Cholesky $mathbf {A}$ que se computó antes para calcular la descomposición de Cholesky ${tilde {mathbf {A} }}$ .

Actualización de rango uno

El caso específico, donde la matriz actualizada ${tilde {mathbf {A} }}$ está relacionado con la matriz $mathbf {A}$ por ${tilde {mathbf {A} }}=mathbf {A} +mathbf {x} mathbf {x} ^{*}$ , es conocido como actualización de rango uno.

Aquí hay una función escrita en la sintaxis de Matlab que realiza una actualización de rango uno:

función [L] = cholupdate()L, x) n = longitud()x); para k = 1:n r = Sqrt()L()k, k)^2 + x()k)^2); c = r / L()k, k); s = x()k) / L()k, k); L()k, k) = r; si k . n L()k+1):n, k) = ()L()k+1):n, k) + s * x()k+1):n) / c; x()k+1):n) = c * x()k+1):n) - s * L()k+1):n, k); final finalfinal

A rank-n update es uno donde para una matriz $mathbf {M}$ uno actualiza la descomposición tal que ${displaystyle {tilde {mathbf {A} }}=mathbf {A} +mathbf {M} mathbf {M} ^{*}}$ . Esto se puede lograr realizando sucesivamente actualizaciones de rango uno para cada una de las columnas de $mathbf {M}$ .

Rebaja de rango uno

A Fecha baja de la categoría 1 es similar a una actualización de rango uno, excepto que la adición es reemplazada por la resta: ${tilde {mathbf {A} }}=mathbf {A} -mathbf {x} mathbf {x} ^{*}$ . Esto solo funciona si la nueva matriz ${tilde {mathbf {A} }}$ sigue siendo positivo.

El código para la actualización de rango uno que se muestra arriba se puede adaptar fácilmente para hacer una actualización de rango uno: uno simplemente necesita reemplazar las dos adiciones en la asignación a r y L ((k+1):n, k) por restas.

Agregar y eliminar filas y columnas

Si tenemos una matriz definida simétrica y positiva ${displaystyle mathbf {A} }$ representado en forma de bloque como

{displaystyle mathbf {A} ={begin{pmatrix}mathbf {A} _{11}&mathbf {A} _{13}\mathbf {A} _{13}^{mathrm {T} }&mathbf {A} _{33}\end{pmatrix}}}

y su factor de Cholesky superior

{displaystyle mathbf {L} ={begin{pmatrix}mathbf {L} _{11}&mathbf {L} _{13}\0&mathbf {L} _{33}\end{pmatrix}},}

entonces para una nueva matriz ${tilde {mathbf {A} }}$ , que es el mismo ${displaystyle mathbf {A} }$ pero con la inserción de nuevas filas y columnas,

{displaystyle {begin{aligned}{tilde {mathbf {A} }}&={begin{pmatrix}mathbf {A} _{11}&mathbf {A} _{12}&mathbf {A} _{13}\mathbf {A} _{12}^{mathrm {T} }&mathbf {A} _{22}&mathbf {A} _{23}\mathbf {A} _{13}^{mathrm {T} }&mathbf {A} _{23}^{mathrm {T} }&mathbf {A} _{33}\end{pmatrix}}end{aligned}}}

estamos interesados en encontrar la factorización Cholesky de ${tilde {mathbf {A} }}$ , que llamamos ${displaystyle {tilde {mathbf {S} }}}$ , sin calcular directamente toda la descomposición.

{displaystyle {begin{aligned}{tilde {mathbf {S} }}&={begin{pmatrix}mathbf {S} _{11}&mathbf {S} _{12}&mathbf {S} _{13}\0&mathbf {S} _{22}&mathbf {S} _{23}\0&0&mathbf {S} _{33}\end{pmatrix}}.end{aligned}}}

Escritura ${displaystyle mathbf {A} setminus mathbf {b} }$ para la solución ${displaystyle mathbf {A} mathbf {x} =mathbf {b} }$ , que se puede encontrar fácilmente para matrices triangulares, y ${displaystyle {text{chol}}(mathbf {M})}$ para la descomposición de Cholesky ${mathbf M}$ , se pueden encontrar las siguientes relaciones:

{displaystyle {begin{aligned}mathbf {S} _{11}&=mathbf {L} _{11},\mathbf {S} _{12}&=mathbf {L} _{11}^{mathrm {T} }setminus mathbf {A} _{12},\mathbf {S} _{13}&=mathbf {L} _{13},\mathbf {S} _{22}&=mathrm {chol} left(mathbf {A} _{22}-mathbf {S} _{12}^{mathrm {T} }mathbf {S} _{12}right),\mathbf {S} _{23}&=mathbf {S} _{22}^{mathrm {T} }setminus left(mathbf {A} _{23}-mathbf {S} _{12}^{mathrm {T} }mathbf {S} _{13}right),\mathbf {S} _{33}&=mathrm {chol} left(mathbf {L} _{33}^{mathrm {T} }mathbf {L} _{33}-mathbf {S} _{23}^{mathrm {T} }mathbf {S} _{23}right).end{aligned}}}

Estas fórmulas se pueden utilizar para determinar el factor de Cholesky después de la inserción de filas o columnas en cualquier posición, si establecemos las dimensiones de fila y columna de manera adecuada (incluso a cero). El problema inverso, cuando tenemos

{displaystyle {begin{aligned}{tilde {mathbf {A} }}&={begin{pmatrix}mathbf {A} _{11}&mathbf {A} _{12}&mathbf {A} _{13}\mathbf {A} _{12}^{mathrm {T} }&mathbf {A} _{22}&mathbf {A} _{23}\mathbf {A} _{13}^{mathrm {T} }&mathbf {A} _{23}^{mathrm {T} }&mathbf {A} _{33}\end{pmatrix}}end{aligned}}}

con descomposición de Cholesky conocida

{displaystyle {begin{aligned}{tilde {mathbf {S} }}&={begin{pmatrix}mathbf {S} _{11}&mathbf {S} _{12}&mathbf {S} _{13}\0&mathbf {S} _{22}&mathbf {S} _{23}\0&0&mathbf {S} _{33}\end{pmatrix}}end{aligned}}}

y desea determinar el factor de Cholesky

{displaystyle {begin{aligned}mathbf {L} &={begin{pmatrix}mathbf {L} _{11}&mathbf {L} _{13}\0&mathbf {L} _{33}\end{pmatrix}}end{aligned}}}

de la matriz ${displaystyle mathbf {A} }$ con filas y columnas removidas,

{displaystyle {begin{aligned}mathbf {A} &={begin{pmatrix}mathbf {A} _{11}&mathbf {A} _{13}\mathbf {A} _{13}^{mathrm {T} }&mathbf {A} _{33}\end{pmatrix}},end{aligned}}}

produce las siguientes reglas:

{displaystyle {begin{aligned}mathbf {L} _{11}&=mathbf {S} _{11},\mathbf {L} _{13}&=mathbf {S} _{13},\mathbf {L} _{33}&=mathrm {chol} left(mathbf {S} _{33}^{mathrm {T} }mathbf {S} _{33}+mathbf {S} _{23}^{mathrm {T} }mathbf {S} _{23}right).end{aligned}}}

Note que las ecuaciones anteriores que implican encontrar la descomposición Cholesky de una nueva matriz son todas de la forma ${displaystyle {tilde {mathbf {A} }}=mathbf {A} pm mathbf {x} mathbf {x} ^{*}}$ , lo que les permite ser eficientemente calculados utilizando los procedimientos de actualización y actualización detallados en la sección anterior.

Prueba para matrices semidefinidas positivas

Prueba mediante argumento limitante

Los algoritmos anteriores muestran que cada matriz definida positiva $mathbf {A}$ tiene una descomposición de Cholesky. Este resultado puede extenderse al caso semidefinido positivo por un argumento limitador. El argumento no es totalmente constructivo, es decir, no da ningún algoritmo numérico explícito para calcular los factores de Cholesky.

Si $mathbf {A}$ es un $ntimes n$ matriz semi-definida positiva, entonces la secuencia ${textstyle left(mathbf {A} _{k}right)_{k}:=left(mathbf {A} +{frac {1}{k}}mathbf {I} _{n}right)_{k}}$ consiste en matrices definidas positivas. (Esta es una consecuencia inmediata de, por ejemplo, el teorema de mapeo espectral para el cálculo funcional polinomio.) También,

{displaystyle mathbf {A} _{k}rightarrow mathbf {A} quad {text{for}}quad krightarrow infty }

en la norma del operador. Desde el caso definitivo positivo, cada uno ${mathbf {A}}_{k}$ tiene la descomposición de Cholesky ${displaystyle 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 $leq$ Oportunidad $M_{n}({mathbb {C}})$ equipado con la norma del operador es un álgebra C*. Así que... ${displaystyle left(mathbf {L} _{k}right)_{k}}$ es un conjunto consolidado en el espacio de Banach de los operadores, por lo tanto relativamente compacto (porque el espacio vectorial subyacente es finito-dimensional). En consecuencia, tiene una subsecuencia convergente, también denotada por ${displaystyle left(mathbf {L} _{k}right)_{k}}$ , con límite $mathbf{L}$ . Se puede comprobar fácilmente que esto $mathbf{L}$ tiene las propiedades deseadas, es decir. ${displaystyle mathbf {A} =mathbf {L} mathbf {L} ^{*}}$ , y $mathbf{L}$ es triangular inferior con entradas diagonales no negativas: para todos $x$ y $y$ ,

{displaystyle langle mathbf {A} x,yrangle =leftlangle lim mathbf {A} _{k}x,yrightrangle =langle lim mathbf {L} _{k}mathbf {L} _{k}^{*}x,yrangle =langle mathbf {L} mathbf {L} ^{*}x,yrangle ,.}

Por lo tanto, ${displaystyle mathbf {A} =mathbf {L} mathbf {L} ^{*}}$ . Debido a que el espacio vectorial subyacente es finito-dimensional, todas las topologías en el espacio de los operadores son equivalentes. Así que... ${displaystyle left(mathbf {L} _{k}right)_{k}}$ tiende a $mathbf{L}$ en la norma ${displaystyle left(mathbf {L} _{k}right)_{k}}$ tiende a $mathbf{L}$ A la entrada. Esto implica, a su vez, que ${displaystyle mathbf {L} _{k}}$ es triangular inferior con entradas diagonales no negativas, $mathbf{L}$ es también.

Prueba por descomposición QR

Vamos $mathbf {A}$ ser una matriz hermitiana semi-definida positiva. Entonces se puede escribir como producto de su matriz de raíz cuadrada, ${displaystyle mathbf {A} =mathbf {B} mathbf {B} ^{*}}$ . Ahora la descomposición QR se puede aplicar a ${displaystyle mathbf {B} ^{*}}$ , resultando en ${displaystyle mathbf {B} ^{*}=mathbf {Q} mathbf {R} }$ , donde $mathbf {Q}$ es unitario y $mathbf {R}$ es triangular superior. Inserción de la descomposición en los rendimientos originales de igualdad ${displaystyle A=mathbf {B} mathbf {B} ^{*}=(mathbf {QR})^{*}mathbf {QR} =mathbf {R} ^{*}mathbf {Q} ^{*}mathbf {QR} =mathbf {R} ^{*}mathbf {R} }$ . Ajuste ${displaystyle mathbf {L} =mathbf {R} ^{*}}$ completa la prueba.

Generalización

La factorización Cholesky se puede generalizar a matrices (no necesariamente finitas) con entradas del operador. Vamos ${{mathcal {H}}_{n}}$ ser una secuencia de espacios de Hilbert. Considere la matriz del operador

mathbf {A} ={begin{bmatrix}mathbf {A} _{11}&mathbf {A} _{12}&mathbf {A} _{13}&;\mathbf {A} _{12}^{*}&mathbf {A} _{22}&mathbf {A} _{23}&;\mathbf {A} _{13}^{*}&mathbf {A} _{23}^{*}&mathbf {A} _{33}&;\;&;&;&ddots end{bmatrix}}

actuando sobre la suma directa

{displaystyle {mathcal {H}}=bigoplus _{n}{mathcal {H}}_{n},}

donde cada

mathbf {A} _{ij}:{mathcal {H}}_{j}rightarrow {mathcal {H}}_{i}

es un operador acotado. Si A es positivo (semidefinido) en el sentido de que para todo k finito y para cualquier

{displaystyle hin bigoplus _{n=1}^{k}{mathcal {H}}_{k},}

tenemos $langle h,mathbf {A} hrangle geq 0$ , entonces existe una matriz de operador triangular inferior L tales que A = LL*. Uno también puede tomar las entradas diagonales L para ser positivo.

Implementaciones en bibliotecas de programación

Lenguaje de programación C: la Biblioteca Científica GNU proporciona varias implementaciones de la descomposición Cholesky.
Maxima sistema de álgebra de computadora: función cholesky computes Cholesky descomposición.
GNU Octave sistema de computaciones numéricas proporciona varias funciones para calcular, actualizar y aplicar una descomposición Cholesky.
La biblioteca LAPACK proporciona una implementación de alto rendimiento de la descomposición Cholesky que se puede acceder desde Fortran, C y la mayoría de los idiomas.
En Python, la función cholesky de la numpy.linalg módulo realiza Descomposición de Cholesky.
En Matlab, el chol función da la descomposición de Cholesky. Note que chol utiliza el factor triangular superior de la matriz de entrada por defecto, es decir, computes ${displaystyle A=R^{*}R}$ Donde $R$ es triangular superior. En cambio, se puede pasar una bandera para utilizar el factor triangular inferior.
En R, el chol función da la descomposición de Cholesky.
En Julia, la cholesky función de la LinearAlgebra biblioteca estándar da la descomposición de Cholesky.
En Mathematica, la función "CholeskyDecomposition" se puede aplicar a una matriz.
En C+++, varias bibliotecas lineales de álgebra soportan esta descomposición:
- El Armadillo (C++) suministra el comando chol para realizar la descomposición de Cholesky.
- La biblioteca Eigen suministra factorizaciones Cholesky tanto para matrices escasas como densas.
- En el paquete ROOT, TDecompChol la clase está disponible.

En Analytica, la función Decompose da la descomposición de Cholesky.
La biblioteca de Matemáticas de Apache Commons tiene una implementación que se puede utilizar en Java, Scala y cualquier otro idioma JVM.

Contenido relacionado

Más resultados...