Tópicos selectos de computación científica I
Cálculo de matriz inversa con factorización cholesky y factorización LU Williams Antonio Pantoja Laces 1 , Angélica Alejandra Serrano Rubio 2 , Javier García Vieyra 3 1−2−3
Dpto. de Computación Computación,,
CINVESTAV-IPN, México, D.F, México 1
{ wpantoja,
2
aserrano,
3
jgarcia}@computacion.cs.cinvestav.mx
2 de julio de 2014 Resumen
Invertir matrices es un problema clásico y que puede ser muy complicado para matrices muy grandes. Hay muchas formas de simplificar esta tarea para tipos especiales de matrices.Entre estas maneras existe la de transformar la matriz en un conjunto de matrices superiores o inferiores, en este documento documento se presentan presentan dos de estos métodos, métodos, específicamen específicamente te utilizando la factorización de cholesky y la otra utilizando la factorización LU , los cuales pueden ayudarnos a ahorrar tiempo en la búsqueda de la matriz inversa.
1. INTR INTROD ODUC UCCI CIÓN ÓN En algebra lineal la factorización de una matriz es la descomposición de la misma como producto de dos o mas matrices según una forma canónica. Mientras que la multiplicación de matrices implica una síntesis de datos(combinando el efecto de dos o mas transformaciones lineales en una sola matriz), la factorización de matrices es un análisis de datos. En el lenguaje de las ciencias de la computación, la expresión de A como un producto equivale a un preprocesamiento de los datos de A, A , el cual organiza esos datos en dos o más partes cuyas 1
estructuras son más útiles de algún modo, quizá por ser más accesibles para realizar cálculos con ellas. La factorización de matrices es el centro de varios programas de computadora usados ampliamente en aplicaciones, las factorizaciones de matrices son un punto clave en el desarrollo de este documento, ya que el documento se enfoca en encontrar la inversa de una matriz utilizando dos de los métodos existentes para la factorización de matrices, el primero es la factorización LU , y el segundo la factorización de cholesky, en las siguientes secciones se presenta el marco teorico con los principales conceptos y características de estos métodos, una sección de experimentos con la implementación de los algoritmos necesarios para dichos métodos y algunas conclusiones obtenidas después de realizar los experimentos.
2. MATERIALES Y MÉTODOS En esta sección se describen los materiales y los métodos que se utilizaron para la implementación de los algoritmos necesarios para el cálculo de la matriz inversa por medio de la factorización LU y la factorización de cholesky.
2.1. Lenguaje C C es un lenguaje de programación creado en 1972 por Dennis M. Ritchie en los Laboratorios Bell como evolución del anterior lenguaje B, a su vez basado en BCPL. Al igual que B, es un lenguaje orientado a la implementación de Sistemas Operativos, concretamente Unix. C es apreciado por la eficiencia del código que produce y es el lenguaje de programación más popular para crear software de sistemas, aunque también se utiliza para crear aplicaciones. Se trata de un lenguaje de tipos de datos estáticos, débilmente tipificado, de medio nivel pero con muchas características de bajo nivel. Dispone de las estructuras típicas de los lenguajes de alto nivel pero, a su vez, dispone de construcciones del lenguaje que permiten un control a muy bajo nivel. Los compiladores suelen ofrecer extensiones al lenguaje que posibilitan mezclar código en ensamblador con código C o acceder directamente a memoria o dispositivos periféricos.
2.2. Matrices Una matriz es un arreglo bidimensional de números (llamados entradas de la matriz) ordenados en filas (o renglones) y columnas, donde una fila es cada una de las líneas horizontales de la matriz y una columna es cada una de las líneas verticales. A una matriz con n filas y m columnas se le denomina matriz n-por-m (escrito n × m) donde n, m ∈ N − {0}. El conjunto de las matrices de tamaño n × m se representa como M n m (K), donde K es el campo al cual pertenecen las entradas. El tamaño de una matriz siempre se da con el número de filas primero y el número de columnas después. Dos matrices se dice que son iguales si tienen el mismo tamaño y los mismos elementos en las mismas posiciones. ×
2
A la entrada de una matriz que se encuentra en la fila i-ésima y la columna j-ésima se le llama entrada i, j o entrada (i, j)-ésimo de la matriz. En estas expresiones también se consideran primero las filas y después las columnas. Casi siempre se denotan a las matrices con letras mayúsculas mientras que se utilizan las correspondientes letras en minúsculas para denotar las entradas de las mismas. Por ejemplo, al elemento de una matriz A de tamaño n × m que se encuentra en la fila i-ésima y la columna j -ésima se le denota como a ij , donde 1 ≤ i ≤ ny1 ≤ j ≤ m. Cuando se va a representar explícitamente una entrada la cuál está indexada con un i o un j con dos cifras se introduce una coma entre el índice de filas y de columnas. Así por ejemplo, la entrada que está en la primera fila y la segunda columna de la matriz A de tamaño 50 × 100 se representa como a 1,2 mientras que la entrada que está en la fila número 23 y la columna 100 se representa como a 23,100 .
2.3. Factorización de Cholesky En matemáticas, la factorización o descomposición de Cholesky toma su nombre del matemático André-Louis Cholesky, quien encontró que una matriz simétrica definida positiva puede ser descompuesta como el producto de una matriz triangular inferior y la traspuesta de la matriz triangular inferior. La matriz triangular inferior es el triángulo de Cholesky de la matriz original positiva definida. El resultado de Cholesky ha sido extendido a matrices con entradas complejas. Es una manera de resolver sistemas de ecuaciones matriciales y se deriva de la factorización LU con una pequeña variación. Cualquier matriz cuadrada A con pivotes no nulos puede ser escrita como el producto de una matriz triangular inferior L y una matriz triangular superior U ; esto recibe el nombre de factorización LU . Sin embargo, si A es simétrica y definida positiva, se pueden escoger los factores tales que U es la transpuesta de L, y esto se llama la descomposición o factorización de Cholesky. Tanto la descomposición LU como la descomposición de Cholesky son usadas para resolver sistemas de ecuaciones lineales. Cuando es aplicable, la descomposición de Cholesky es dos veces más eficiente que la descomposición LU . En general, si A es Hermitiana y definida positiva, entonces A puede ser descompuesta como ∗
A = LL , donde L es una matriz triangular inferior con entradas diagonales estrictamente positivas y L∗ representa la conjugada traspuesta de L. Esta es la descomposición de Cholesky. Para encontrar la factorización , bastaría ver la forma de L y observar las ecuaciones que el producto derecho nos conduce al igualar elementos:
3
a11 a21
a12 a22
a31 .. .
a32 .. .
an1
an2
··· ··· .. .
a1n a2n a3n .. .
..
. ···
ann
así obtendríamos que:
=
l11 l21
0 l22
l31 .. .
l32 .. .
ln1
ln2
··· ··· .. .
0 0 0 .. .
..
. ···
0
l11 0
l21 l22
0 .. .
0 .. .
0
0
··· ··· .. . ..
. ···
ln1 ln2 ln3 .. . lnn
√
2 −→ l11 = 2 a11 = l 11 a21 = l 21 l11 ←− a21/l11, ...l n1 = a n1 /l11 2 21
2 22
a22 = l + l a32 = l 31 l21
←− l 2 = + l l ←− l 2
(a22 − a221)
32
32 22
= (a32 − l31l21)l22, etc
y de manera general, para i = 1, ..., n y j = i + 1, ...n :
− − i−1
lii =
l2i k
aii
k=1
i−1
l ji =
aji
l jk lik
/lii
k=1
Ahora bien, ya que A es simétrica y definida positiva, podemos asegurar que los elementos sobre la diagonal de L son positivos y los restantes elementos reales desde luego.
2.4. Factorización LU En el álgebra lineal, la factorización o descomposición LU (del inglés LowerUpper) es una forma de factorización de una matriz como el producto de una matriz triangular inferior y una superior. Debido a la inestabilidad de este método, deben tenerse en cuenta algunos casos especiales, por ejemplo, si uno o varios elemento de la diagonal principal de la matriz a factorizar es cero, es necesario premultiplicar la matriz por una o varias matrices elementales de permutación. Método llamado factorización P A = LU o LU con pivote. Esta descomposición se usa en el análisis numérico para resolver sistemas de ecuaciones (más eficientemente) o encontrar las matrices inversas. Sea A una matriz no singular (si lo fuera, entonces la descomposición podría no ser única) A = LU, donde L y U son matrices inferiores y superiores triangulares respectivamente. Para matrices 3 × 3, esto es: 4
a11 a21 a31
a12 a22 a32
a13 a23 a33
=
1
0 1
l21 l31
l32
0 0 1
u11 0 0
u12 u22 0
u13 u23 u33
2.5. Matriz Inversa Se dice que una matriz cuadrada A es inversible, si existe una matriz B con la propiedad de que: A ∗ B = B ∗ A = I siendo I la matriz identidad. Se denomina a B la matriz inversa de A y se denota por A 1 . Una matriz se dice que es inversible o regular si posee inversa, en caso contrario se dice que es ssingular o degenerada. Una matriz es singular si y solo si su determinante es nulo. −
2.6. Sustitución hacia atrás(Back-substitution) 2.7. Algoritmos En esta sección se presentan los algoritmos utilizados para la implementación del algoritmo para el cáclulo de la matriz inversa utilizando la factorización LU y factorización de cholesky.
2.7.1.
Algoritmo de Tridiagonalización
Para la parte de tridiagonalización de matrices se utilizó el algoritmo que se presenta en el Algoritmo 1, el cual recibe como entrada la matriz con la cual se operara A de tamaño N × N y un valor h el cual, es el valor de discretización y determina la granularidad.
2.7.2.
Algoritmo Factorización LU
También se utilizó un algoritmo para la factorización LU después de tener la matriz tridiagonal, este algoritmo es el que se muestra en el Algoritmo 2, el cual recibe como entrada la matriz obtenida el algoritmo de tridiagonalización de tamaño N y como salida devuelve las matrices L y U
5
Algoritmo 1 Tridiagonalizar matriz 1: para i = 1 hasta N hacer 2: para j = 1 hasta N hacer si i == j && i == 0 entonces 3:
13:
A[i][ j] = 2/h2 A[i][ j + 1] = −1/h2 j = j + 1 si no, si i == j && i == N − 1 entonces A[i][ j] = 2/h2 A[i][ j − 1] = −1/h2 si no, si i == j entonces A[i][ j] = 2/h2 A[i][ j − 1] = A[i][ j + 1] = −1/h2 j + +
14:
si no
4: 5: 6: 7: 8: 9: 10: 11: 12:
15:
A[i][ j] = 0
fin si fin para 18: fin para 16: 17:
Algoritmo 2 Factorización LU 1: para k = 1 hasta N hacer 2: 3: 4: 5: 6: 7: 8: 9:
Especificar un valor para l ii o u kk Calcular el otro termino mediante 1 lkk ukk = a kk − ks=1 lks usk para j = k + 1 hasta N hacer 1 ukj ← akj − ks=1 lks usk /lkk
−
fin para para i = k + 1 hasta N hacer
ljk ← aik −
fin para 11: fin para 10:
−
k−1 s=1
lis usk /lkk
6
2.7.3.
Algoritmo de multiplicación de matrices
2.7.4.
Principal
En algoritmo principal es el que se muestra en el Algoritmo 4, en donde primero recibe como entrada la matriz que se quiere diagonalizar, y de la cual se tiene que encontrar el espectro (los valores propios), primero se tridiagonaliza la matriz, para después hacer la factorización LU , teniendo la factorización LU se hace la multiplicación UL, se verifica si la matriz es diagonal, de ser asi, finaliza el algoritmo, en caso contrario se sigue iterando hasta conseguir la diagonalización. Algoritmo 3 Diagonalización iterativa LU
Entrada: A ∈ Rn n Salida: D ∈ Rn 2: D ← A 3: mientras D no sea diagonal hacer D → LU 4: D ← UL 5: ×
1:
6: fin 7:
n
×
mientras
Regresa D
3. RESULTADOS 4. CONCLUSIONES 5. REFERENCIAS 1. Evguenii Kurmyshev, “Fundamentos de métodos matemáticos para física e ingeniería”, 4ta Edición, Editorial Limusa, 2009. 2. C.Cerjan,“Numerical Grid Methods and Their Application to Schrödinger’s Equation”,Ed. Springer, 1993. 3. Valor y Vector Propio. En Wikipedia, recuperado el 12 de junio de 2014. http://es.wikipedia.org/wiki/Vector_propio_y_valor_propio 4. Matriz Diagonal. En Wikipedia recuperado el 12 de junio de 2014. http://es.wikipedia.org/wiki/Matriz_diagonal 5. Matriz Diagonalizable. En Wikipedia recuperado el 12 de junio de 2014. http://es.wikipedia.org/wiki/Matriz_diagonalizable
7
