Pidiendo una solución al problema 21.5 en Numerical linear algebra por Lloyd N. Trefethen, David Bau III

Question

No consigo averiguar cómo es posible que las entradas de la matriz L no puedan exceder de 1, ya que el pivoteo simétrico sólo afecta a la matriz diagonal. Considere una matriz donde las entradas diagonales son muy pequeñas. ¿Cómo podría |l_{ij}|<=1 posible?

Instantánea del problema:

Por favor, que alguien me ayude...

Answer 1

Puedes tener una descomposición $PAP^T=LDL^T$ donde $D$ es una matriz tridiagonal y $L$ es una matriz triangular inferior con diagonal unitaria y todas las entradas no diagonales con valor absoluto $\leq 1$ . Encontré la idea y el algoritmo para hacer esto en este documento .

Como mencionas, al utilizar el pivote simétrico, sólo puedes intercambiar la posición de las entradas de la diagonal y no puedes poner una entrada fuera de la diagonal. Aasen propuesta para utilizar el pivote simétrico en el paso $i$ para poner el mayor elemento fuera de la diagonal de la columna $i$ en la posición $(i+1,i)$ y luego eliminar todas las entradas posteriores $(i+2:m,i)$ en la columna utilizando ese elemento como pivote.

Digamos que sí: $$A=\begin{pmatrix}1 & 2 & 3 & 4 \\ 2 & 3 & 5 & 6 \\ 3 & 5 & 4 & 9 \\ 4 & 6 & 9 & 7\end{pmatrix}$$ .

El primer paso es pivotar $A(4,1)$ a la posición $(2,1)$ utilizando la matriz de permutación $P_1=\begin{pmatrix}1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1\\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0\end{pmatrix}$ :

$$P_1AP_1=\begin{pmatrix}1 & 4 & 3 & 2\\4 & 7 & 9 & 6\\3 & 9 & 4 & 5\\2 & 6 & 5 & 3\end{pmatrix}$$ .

Ahora podemos construir una matriz diagonal unitaria $L_1$ que establecerá las entradas de la columna $1$ en filas $3$ y $4$ a $0$ :

$$L_1=\begin{pmatrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & -3/4 & 1 & 0 \\ 0 & -1/2 & 0 & 1\end{pmatrix}$$

Utilizando $L_1$ nos encontramos con que:

$$L_1P_1AP_1L_1^T=\begin{pmatrix}1 & 4 & 0 & 0\\4 & 7 & 3.75 & 2.5\\0 & 3.75 & -5.5625 & -1.375 \\ 0 & 2.5 & -1.375 & -1.25\end{pmatrix}$$ .

Este paso se repite, en caso de una $m\times m$ matriz en conjunto $m-2$ tiempos.

Lo principal que hay que observar aquí es que cada $L_i$ tiene entradas subdiagonales no nulas en la columna $(i+1)$ sólo y que cada $P_i$ hace pivotar el mayor elemento subdiagonal de la columna $i$ a la posición $(i+1,i)$ y no $(i,i)$ .

A continuación se muestra el pseudocódigo del algoritmo descrito anteriormente que aplica la permutación y la eliminación simétrica en la matriz completa en cada paso. Para la implementación completa en Octave, véase aquí .

input: A, mxm matrix
output: L, D, P

L=I %identity matrix
D=tril(A)
for i=1:m-2
    [sigma, r] = max(abs(D(i+1:m,i)));

    if sigma == 0
        continue
    endif

    r=r+i
    Pi=permute(I,i+1,r)
    L=Pi*L*Pi

    %symmetric swap on a triangular matrix
    D=symmetricRowAndColumnSwap(D,i+1,r) 

    L(i+2:m,i+1)=D(i+2:m,i)/D(i+1,i)

    % eliminating the i-th column
    for j=m:-1:i+2
        D(j,i)=0
        D(j,i+1:j)=D(j,i+1:j)-D(i+1:j,i+1)'*Li(j,i+1)
    endfor

    % eliminating the i-th row
    for j=i+2:m
        D(j:m,j)=D(j:m,j)-Li(j,i+1)*D(j:m,i+1)
    endfor

    P=Pi*P
endfor

D=D+tril(D,-1)’

Se necesita $\approx 2/3m^3$ flops para ejecutar el algoritmo como se hace en el pseudocódigo anterior. Sin embargo, Aasen realmente implementa la actualización de manera diferente, actualizando una columna a la vez en cada paso. Esto requiere $\approx m^3/3$ fracasos. Además del papel original, también he utilizado esta presentación para ayudarme a implementar la actualización columna por columna y entender la complejidad del algoritmo.

La implementación de la actualización más rápida, columna por columna, es aquí .

Esta versión aprovecha el hecho de que en cada paso $i$ tenemos toda la información que necesitamos para actualizar la columna $i+1$ en la matriz reducida. Al explicar el algoritmo, ignoraré aquí las permutaciones para facilitar la explicación de la eliminación columna por columna:

Después de realizar el algoritmo para $i$ pasos que conocemos las matrices $L_i$ que hacen la eliminación en cada paso. Podemos escribir: $$\begin{eqnarray}\nonumber A=L_1^{-1} \cdots L_i^{-1} D^{(i)} (L_i^{-1})^T \cdots (L_{1}^{-1})^T\end{eqnarray}$$ $$\begin{eqnarray} \nonumber L^{(i)}_{m \times m}=L_1^{-1} \cdots L_i^{-1} =\begin{pmatrix} L_{11} & 0 \\ L_{21} & L_{22}^{(i)} \end{pmatrix}\end{eqnarray}$$

con $L_{22}^{(i)}= \begin{pmatrix} l_{i+1} & 0 \end{pmatrix}_{(m-i) \times (m-i)}$

$L^{(i)}$ es la matriz de eliminación $L$ después del paso $i$ . $l_{i+1}$ es un vector columna de dimensión $(m-i) \times 1$ . $L_{11; (i \times i)}$ es una matriz triangular unitaria inferior y $L_{21; (m-i) \times i}$ es una matriz rectangular con todas las entradas potencialmente no nulas (excepto las entradas de la primera columna).

Podemos factorizar A y $D^{(i)}$ de manera similar:

$$A=\begin{pmatrix} A_{11(i \times i)} & A_{21(i \times m-i)}^T \\ A_{21(m-i \times i)} & A_{22(m-i \times m-i)} \end{pmatrix}$$

$$D^{(i)}=\begin{pmatrix} D_{11(i \times i)} & D_{21(i \times m-i)}^T \\ D_{21(m-i \times i)} & D_{22(m-i \times m-i)}^{(i)} \end{pmatrix}$$

donde $D_{11}$ es un bloque tridiagonal de la matriz final, $D_{21}$ es también un bloque de la matriz final $D$ y sólo tiene un valor distinto de cero en la intersección de su primera fila y última columna y $D_{22}^{(i)}$ es la parte que aún debe reducirse más después del paso $i$ .

Lo tenemos:

$$\begin{pmatrix}A_{11} & A_{21}^T \\ A_{21} & A_{22}\end{pmatrix}=\begin{pmatrix} L_{11} & 0 \\L_{21} & L_{22}^{(i)} \end{pmatrix} \cdot \begin{pmatrix} D_{11} & D_{21}^T \\ D_{21} & D_{22}^{(i)} \end{pmatrix} \cdot \begin{pmatrix} L_{11}^T & L_{21}^T \\ 0 & L_{22}^{(i)T} \end{pmatrix} $$

De ahí sacamos:

$$A_{22}=L_{21}D_{11}L_{21}^T + L_{22}^{(i)}D_{21}L_{21}^T + L_{21}D_{21}^TL_{22}^{(i)T} + L_{22}^{(i)}D_{22}^{(i)}L_{22}^{(i)T}$$

y

$$L_{22}^{(i)}D_{22}^{(i)}L_{22}^{(i)T} = A_{22}-L_{21}D_{11}L_{21}^T - L_{22}^{(i)}D_{21}L_{21}^T - L_{21}D_{21}^TL_{22}^{(i)T}$$

De la ecuación anterior podemos obtener la primera columna de $L_{22}^{(i)}D_{22}^{(i)}L_{22}^{(i)T}$ . Lo que realmente queremos es la primera columna de $D_{22}^{(i)}$ . Por suerte, tenemos la información que necesitamos para calcular esa columna. $L_{22}^{(i)}$ sólo tiene una columna con entradas no-diagonales distintas de cero, su primera columna. Ésta resulta ser la parte no nula del $(i+1)$ -Columna de $L_i^{-1}$ la inversa de la matriz de reducción utilizada para crear ceros en el $i$ -columna de la matriz reducida en el paso $i$ . Podemos usar esto para obtener la primera columna de $D_{22}^{(i)}$ que es todo lo que necesitamos para construir $L_{i+1}$ y continuar el proceso descrito columna por columna.

Para las operaciones exactas, aquí está el fragmento de código que hace la actualización de la columna en el paso $i$ :

L(i+2:m,i+1)=D(i+2:m,i)/D(i+1,i);
D(i+2:m,i)=0;

v=zeros(i,1);
if i > 1
    v(1)=D(1,2)*L(i+1,2);
    for j=2:i-1
        v(j)=D(j,j-1:j)*L(i+1,j-1:j)'+D(j+1,j)*L(i+1,j+1);
    endfor
    v(i)=D(i,i-1:i)*L(i+1,i-1:i)';
    endif

D(i+1:m,i+1)=D(i+1:m,i+1)-L(i+1:m,1:i)*v-L(i+1:m,i+1)*D(i+1,i)*L(i+1,i)-L(i+1:m,i)*D(i+1,i);
D(i+2:m,i+1)=D(i+2:m,i+1)-L(i+2:m,i+1)*D(i+1,i+1);

(b) La forma que toma la factorización es $PAP^T=LDL^T$ donde $P$ es una matriz de permutación y L es una matriz triangular unitaria cuya primera columna son todas las $0$ s aparte del 1 de la diagonal. Además, todas las entradas subdiagonales de $L$ tienen valor absoluto menor que 1. D es una matriz simétrica tridiagonal.

(c) La primera implementación del algoritmo requiere $\approx 2/3m^3$ fracasos y la segunda mitad de eso.

Del pseudocódigo de la primera versión del algoritmo vemos que cada actualización de entradas requiere una multiplicación y una suma. Hay $\approx (m-i-1)*(m-i)/2$ entradas a actualizar al aprovechar la simetría que viene a $\approx (m-i)^2/2*2=(m-i)^2$ para la compensación de columnas y la misma cantidad para la compensación de filas. La actualización se produce $m-2$ veces por lo que tenemos:

$$\sum_{i=1}^{m-2}2(m-i)^2 \approx \sum_{k=1}^{m}2k^2\approx 2m^3/3$$

La actualización basada en columnas requiere $\approx m^3/3$ flops. El fragmento de código anterior calcula un vector $v$ que es la primera columna de $D_{11}L_{21}^T$ así que $v=D_{11}L_{21}(1,:)^T$ . Dado que cada fila de $D_{11}$ tiene como máximo $3$ valores distintos de cero, calculando $v$ toma $\approx 5i$ operaciones (como máximo 3 multiplicaciones y 2 sumas por cada entrada de $v$ ).

Sólo necesitamos la primera columna de $L_{21}T_{11}L_{21}^T$ que es igual a $L_{21}v$ . El cálculo de esta columna requiere $\approx 2i(m-i)$ fracasos.

También necesitamos las primeras columnas de $L_{22}^{(i)}D_{21}L_{21}^T$ y $L_{21}D_{21}^TL_{22}^{(i)T}$ . Desde $D_{21}$ sólo tiene una entrada no nula cada una de estas columnas toma $2(m-i)$ fracasos.

A esto hay que añadir las sustracciones del $3$ elementos de la primera columna de $A_{22}$ que viene a $\approx{3(m-i)}$ y las operaciones necesarias para obtener la primera columna de $D_{22}$ de $L_{22}^{(i)}D_{22}^{(i)}$ . Esto requiere $\approx 2(m-i)$ operaciones.

Hay $m-1$ pasos y todas las sumas de los otros factores excepto $2i(m-i)$ son del orden $m^2$ el número de flops está dominado por el término $2i(m-i)$ :

$$\sum_{i=1}^{m-1}=2m\sum_{i=1}^{m-1}i - 2\sum_{i=1}^{m-1}i^2 = m^3/3-m/3 \approx m^3/3$$

De esto se desprende que el número de flops para la eliminación columna por columna es de alrededor de $m^3/3$ .