Descomposición LU: ¿diferencia entre el cálculo manual y el solver?

Question

Tengo un $3 \times 3$ matriz $A$ y tienen que realizar el $LU$ Factorización (1)

$$ A = \begin{bmatrix} 1 & 2 & 3 \\ 1 & -1 & 3 \\ -2 & -10 & -2 \end{bmatrix}$$

Utilizando la reducción de filas, primero restaría la segunda fila por (1/1) la primera fila, lo que da (2):

$$ \begin{bmatrix} 1 & 2 & 3 \\ 0 & -3 & 0 \\ -2 & -10 & -2 \end{bmatrix}$$

Entonces restaría la tercera fila por (-2/1) la primera fila, lo que da (3)

$$ \begin{bmatrix} 1 & 2 & 3 \\ 0 & -3 & 0 \\ 0 & -6 & 4 \end{bmatrix}$$

Por último, restaría la tercera fila por (-6/2) la segunda fila, lo que daría (4):

$$ \begin{bmatrix} 1 & 2 & 3 \\ 0 & -3 & 0 \\ 0 & 0 & 4 \end{bmatrix} = U$$

con el factor encontrado en (2) (3) y (4) puedo obtener la matriz $U$

$$L = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ -2 & 2 & 1 \end{bmatrix}$$

Puedo verificar que $L \times U = A$

Sin embargo, cuando calculo [L, U] = lu(A) con Octave, obtengo $$ U = \begin{bmatrix} -2 & -10 & -2 \\ 0 & -6 & 2 \\ 0 & 0 & 1 \end{bmatrix}$$

y

$$ L= \begin{bmatrix} -0.5 & 0.5 & 1 \\ -0.5 & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix}$$

Aquí está el octave / matlab código:

A = [1 2 3; 1 -1 3; -2 -10 -2];
U_hand = [1 2 3; 0 -3 0; 0 0 4];
L_hand = [1 0 0; 1 1 0; -2 2 1];

L_hand * U_hand

[L, U] = lu(A)

¿Cómo se pueden explicar las diferencias? Probablemente me equivoque en alguna parte, pero ¿dónde?

Answer 1

No deberías haber conseguido eso para tu descomposición de la LU. Usé python que usa el mismo LAPACK

import scipy.linalg
import 

A = scipy.array([[1 ,2,3],[1, -1, 3 ] ,[-2,-10,-2]])
P,L,U = scipy.linalg.lu(A)

L
Out[3]: 
array([[ 1. ,  0. ,  0. ],
       [-0.5,  1. ,  0. ],
       [-0.5,  0.5,  1. ]])

U
Out[4]: 
array([[ -2., -10.,  -2.],
       [  0.,  -6.,   2.],
       [  0.,   0.,   1.]])

Tras una nueva inspección. La referencia dice

Es el producto de la matriz de permutación y la matriz L

Con dos argumentos de salida, devuelve las formas permutadas de las matrices triangulares superior e triangulares superiores e inferiores, de forma que A = L * U. Con un argumento de salida argumento y, entonces se devuelve la matriz devuelta por las rutinas LAPACK devuelve. Si la matriz de entrada es dispersa, entonces la matriz L se incrusta en U para dar un valor de retorno similar al caso completo. Tanto para el caso completo y matrices dispersas, lu pierde la información de permutación.

Lo que nos lleva a probarlo en Python.

Atest = np.dot(P,L)

Out[2]: 
array([[-0.5,  0.5,  1. ],
       [-0.5,  1. ,  0. ],
       [ 1. ,  0. ,  0. ]])