Factorizaciones LU con R, Python y Octave

3.3 Factorizaciones LU con `R`

Para realizar una factorización LU con R, podemos utilizar la función LU() introduciendo por parámetro una matriz cuadrada.

La función devolverá una lista con tres componentes: \(P\), \(L\) y \(U\)

Veámoslo con un ejemplo

Ejemplo 1

Encontremos la factorización LU de la siguiente matriz \[A = \begin{pmatrix}1 & 3 & 0 & -1\\ 2 & 1 & -1& 5\\ 0 & -2& 3 & -1\\ 1 & 1 & 3 & 1\end{pmatrix}\]

A = rbind(c(1,3,0,-1), c(2,1,-1,5), c(0,-2,3,-1), c(1,1,3,1))
luA = LU(A)

luA$P

     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    1    0    0
[3,]    0    0    1    0
[4,]    0    0    0    1

En este caso, como no se han permutado filas, la matriz \(P\) es la matriz identidad

luA$L

     [,1] [,2] [,3] [,4]
[1,]    1  0.0    0    0
[2,]    2  1.0    0    0
[3,]    0  0.4    1    0
[4,]    1  0.4    1    1

luA$U

     [,1] [,2] [,3] [,4]
[1,]    1    3  0.0 -1.0
[2,]    0   -5 -1.0  7.0
[3,]    0    0  3.4 -3.8
[4,]    0    0  0.0  3.0

Ejemplo 2

Encontremos ahora la factorización LU de la matriz \[A = \begin{pmatrix}0 & 1 & 3\\ 1 & 3& -2\\ -3 & -2 & -1 \end{pmatrix}\]

A = matrix(c(0,1,3,1,3,-2,-3,-2,-1), byrow = T, nrow = 3, ncol = 3)
luA = LU(A)

luA$P

     [,1] [,2] [,3]
[1,]    0    1    0
[2,]    1    0    0
[3,]    0    0    1

En este caso, podemos ver como sí se han permutado filas, ya que la matriz \(P\) no es la matriz identidad

luA$L

     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]   -3    7    1

luA$U

     [,1] [,2] [,3]
[1,]    1    3   -2
[2,]    0    1    3
[3,]    0    0  -28

Finalmente, también podemos resolver sistemas de ecuaciones lineales aplicando antes la factorización LU a la matriz de coeficientes.

Esto se lleva a cabo con la misma función utilizada hasta el momento: LU(), pero añadiendo un parámetro más. Lo que implica que además de las matrices \(P\),\(L\) y \(U\), la función nos devuelve dos vectores más: \(d\) y \(x\)

\(x\) es la solución del sistema
\(d\) es el vector solución del sistema \(Ld = b\), que luego nos sirve para resolver \(Ux = d\)

Ejemplo 3

Consideremos el sistema

\[\left\{\begin{matrix} &&x_2&+&3x_3&=&1\\ x_1&+&3x_2&-&2x_3 &=& 3\\ -3x_1&-&2x_2&-&x_3&=&-2 \end{matrix}\right.\]

A = rbind(c(0,1,3), c(1,3,-2), c(-3,-2,-1))
b = c(1,3,-2)
sistema = LU(A,b)

Warning in if (!backword) 1L:len else len:1L: the condition has length > 1 and
only the first element will be used

sistema$P

     [,1] [,2] [,3]
[1,]    0    1    0
[2,]    1    0    0
[3,]    0    0    1

sistema$L

     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]   -3    7    1

sistema$U

     [,1] [,2] [,3]
[1,]    1    3   -2
[2,]    0    1    3
[3,]    0    0  -28

sistema$d

     [,1]
[1,]    3
[2,]    1
[3,]    0

sistema$x

     [,1]
[1,]    0
[2,]    1
[3,]    0

3.4 Factorizaciones LU con `Python`

Para realizar una factorización LU con Python, podemos utilizar la función scipy.linalg.lu() introduciendo por parámetro una matriz cuadrada. Para ello, habrá que instalar la librería scipy mediante py_install("scipy").

La función scipy.linalg.lu() devolverá tres matrices: \(P\), \(L\) y \(U\)

Veámoslo con un ejemplo

Ejemplo 1

Encontremos la factorización LU de la siguiente matriz \[A = \begin{pmatrix}1 & 3 & 0 & -1\\ 2 & 1 & -1& 5\\ 0 & -2& 3 & -1\\ 1 & 1 & 3 & 1\end{pmatrix}\]

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

P

array([[0., 1., 0., 0.],
       [1., 0., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

En este caso, aunque no fuese necesaria la permutación, Python la ha realizado

array([[ 1. ,  0. ,  0. ,  0. ],
       [ 0.5,  1. ,  0. ,  0. ],
       [ 0. , -0.8,  1. ,  0. ],
       [ 0.5,  0.2,  1. ,  1. ]])

array([[ 2. ,  1. , -1. ,  5. ],
       [ 0. ,  2.5,  0.5, -3.5],
       [ 0. ,  0. ,  3.4, -3.8],
       [ 0. ,  0. ,  0. ,  3. ]])

Ejemplo 2

Encontremos ahora la factorización LU de la matriz \[A = \begin{pmatrix}0 & 1 & 3\\ 1 & 3& -2\\ -3 & -2 & -1 \end{pmatrix}\]

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

P

array([[0., 0., 1.],
       [0., 1., 0.],
       [1., 0., 0.]])

En este caso, podemos ver como también se han permutado filas, ya que la matriz \(P\) no es la matriz identidad

array([[ 1.        ,  0.        ,  0.        ],
       [-0.33333333,  1.        ,  0.        ],
       [-0.        ,  0.42857143,  1.        ]])

array([[-3.        , -2.        , -1.        ],
       [ 0.        ,  2.33333333, -2.33333333],
       [ 0.        ,  0.        ,  4.        ]])

Para realizar una factorización LU con Python, también podemos utilizar la función scipy.linalg.lu_factor() introduciendo por parámetro una matriz cuadrada.

La función devolverá dos elementos matrices:

Una matriz en cuya parte inferior se corresponde con la matriz triangular inferior \(L\) y cuya parte superior se corresponde con la matriz triangular superior \(U\)
Un vector de índices, piv que indica que la fila \(i\) se ha intercambiado con la fila piv[i]

Veámoslo con un ejemplo

Ejemplo 2

Encontremos ahora la factorización LU de la matriz \[A = \begin{pmatrix}0 & 1 & 3\\ 1 & 3& -2\\ -3 & -2 & -1 \end{pmatrix}\]

A = scipy.array([[0,1,3], [1,3,-2], [-3,-2,-1]])
LU, piv = scipy.linalg.lu_factor(A)

LU

array([[-3.        , -2.        , -1.        ],
       [-0.33333333,  2.33333333, -2.33333333],
       [-0.        ,  0.42857143,  4.        ]])

L #Resultado anterior

array([[ 1.        ,  0.        ,  0.        ],
       [-0.33333333,  1.        ,  0.        ],
       [-0.        ,  0.42857143,  1.        ]])

U #Resultado anterior

array([[-3.        , -2.        , -1.        ],
       [ 0.        ,  2.33333333, -2.33333333],
       [ 0.        ,  0.        ,  4.        ]])

piv

array([2, 1, 2], dtype=int32)

P #Resultado anterior

array([[0., 0., 1.],
       [0., 1., 0.],
       [1., 0., 0.]])

Aquí observamos que la primera fila (la 0), se ha cambiado con la tercera; la segunda se ha quedado tal cual; y la tercera, una vez realizado el primer intercambio, se ha quedado en el sitio.

Con lo visto anteriormente, ahora somos capaces de resolver sistemas utilizando factorización LU con Python. Esto lo podemos hacer con la función scipy.linalg.lu_solve() a la cual introducimos por parámetros la tupla (LU,piv) y el vector de términos independientes del sistema, b.

Esta función devuelve únicamente el vector solución:

Ejemplo 3

Consideremos el sistema

\[\left\{\begin{matrix} &&x_2&+&3x_3&=&1\\ x_1&+&3x_2&-&2x_3 &=& 3\\ -3x_1&-&2x_2&-&x_3&=&-2 \end{matrix}\right.\]

A = scipy.array([[0,1,3], [1,3,-2], [-3,-2,-1]])
LU, piv = scipy.linalg.lu_factor(A)
b = [1, 3, -2]
x = scipy.linalg.lu_solve((LU,piv),b) 
x

array([-5.28677631e-18,  1.00000000e+00, -1.98254112e-18])

3.5 Factorizaciones LU con `Octave`

A la hora de realizar factorizaciones LU con Octave, utilizamos la función lu() introduciéndole por parámetro una matriz cuadrada

Ejemplo 1

Encontremos la factorización LU de la siguiente matriz \[A = \begin{pmatrix}1 & 3 & 0 & -1\\ 2 & 1 & -1& 5\\ 0 & -2& 3 & -1\\ 1 & 1 & 3 & 1\end{pmatrix}\]

A = [1,3,0,-1; 2,1,-1,5; 0,-2,3,-1; 1,1,3,1];
[L,U,P] = lu(A);

P

P =

Permutation Matrix

   0   1   0   0
   1   0   0   0
   0   0   1   0
   0   0   0   1

De nuevo, aunque no era necesario, la función ha llevado a cabo una permutación de filas.

A = [1,3,0,-1; 2,1,-1,5; 0,-2,3,-1; 1,1,3,1]; [L,U,P] = lu(A);
L
U

L =

   1.00000   0.00000   0.00000   0.00000
   0.50000   1.00000   0.00000   0.00000
   0.00000  -0.80000   1.00000   0.00000
   0.50000   0.20000   1.00000   1.00000

U =

   2.00000   1.00000  -1.00000   5.00000
   0.00000   2.50000   0.50000  -3.50000
   0.00000   0.00000   3.40000  -3.80000
   0.00000   0.00000   0.00000   3.00000

Ejemplo 2

Encontremos ahora la factorización LU de la matriz

\[A = \begin{pmatrix}0 & 1 & 3\\ 1 & 3& -2\\ -3 & -2 & -1 \end{pmatrix}\]

A = [0,1,3; 1,3,-2; -3,-2,-1];
[L,U,P] = lu(A)

P

L =

   1.00000   0.00000   0.00000
  -0.33333   1.00000   0.00000
  -0.00000   0.42857   1.00000

U =

  -3.00000  -2.00000  -1.00000
   0.00000   2.33333  -2.33333
   0.00000   0.00000   4.00000

P =

Permutation Matrix

   0   0   1
   0   1   0
   1   0   0

P =

Permutation Matrix

   0   0   1
   0   1   0
   1   0   0

En este caso sí que era necesaria una permutación de filas.

A = [0,1,3; 1,3,-2; -3,-2,-1]; 
[L,U,P] = lu(A);
L
U

L =

   1.00000   0.00000   0.00000
  -0.33333   1.00000   0.00000
  -0.00000   0.42857   1.00000

U =

  -3.00000  -2.00000  -1.00000
   0.00000   2.33333  -2.33333
   0.00000   0.00000   4.00000

Factorizaciones LU con R, Python y Octave

3.3 Factorizaciones LU con R

3.4 Factorizaciones LU con Python

3.5 Factorizaciones LU con Octave

3.3 Factorizaciones LU con `R`

3.4 Factorizaciones LU con `Python`

3.5 Factorizaciones LU con `Octave`