MAESTRÍA EN INVESTIGACIÓN E INNOVACIÓN DE TECNOLOGÍAS DE LA INFORMACIÓN Y COMUNICACIÓN LABORATORY 2: NUMERICAL AND DATA-INTENSIVE COMPUTING (COURSE 2012/2013) SYSTEM OF LINEAR EQUATIONS PROFESOR: DR. CARLOS SANTACRUZ José Luis Carrillo Medina Cesar Naranjo Hidalgo
1. Ejercicios 1. Dada la descomposición descomposición LU general general de una matriz cuadrada de dimensión n. A = LU Preguntas: a) Determine el valor valor de |A| El determinante es una función que le asigna a una matriz de orden n, un único número real llamado el determinante de la matriz. Si A es una matriz de orden n, el determinante de la matriz A lo denotaremos por det(A) o también por DEFINICIÓN Si
(las barras no significan valor absoluto).
2.1(Determinante
de
una
matriz
de
orden
1)
es una matriz de orden uno, entonces det(A)=a.
DEFINICIÓN 2.2(Menores y cofactores de una matriz de orden n) Sea A una matriz de orden elemento
, definimos el menor
de A como el determinante de la matriz que se obtiene al
eliminar la fila i y la columna j de la matriz A. El cofactor elemento
asociado al
de A esta dado por
asociado al
.
DEFINICIÓN 2.3(Determinante de una matriz de orden superior) Si A es una matriz de orden , entonces el determinante de la matriz A es la suma de los elementos de la primera fila f ila de A multiplicados por sus respectivos cofactores.
b) Describa un procedimiento para calcular el inverso inverso de A 1 −
La matriz inversa de A es otra matriz que representamos por A -1 y que verifica: A*A-1 = A-1*A = I
La matriz inversa de A es:
| |
)) -1
A·A
-1
= A · A = I
Propiedades (A * B) - 1 = B - 1 * A - 1 (A - 1) - 1 = A (k * A) - 1 = k - 1 * A - 1 (A t) - 1 = (A - 1) t
Cálculo por el método de Gauss Sea A una matriz cuadrada de orden n. Para calcular la matriz inversa de A, que denotaremos como A-1, seguiremos los siguientes pasos: 1º Construir una matriz del tipo M = (A | I), es decir, A está en la mitad izquierda de M y la matriz identidad I en la derecha. Consideremos una matriz 3x3 arbitraria
La ampliamos con la matriz identidad de orden 3.
2º Utilizando el método Gauss vamos a transformar la mitad izquierda, A, en la matriz identidad, que ahora está a la derecha, y la matriz que resulte en el lado derecho será la matriz inversa: A-1. F2 - F1
F3 + F2
F2 - F3
F1 + F2
(-1) F2
La matriz inversa es:
c) Basada en la descomposición de LU, discuta cuando la Matriz es invertible. Descomposición LU e Inversión de matrices Descomposición LU El principal atractivo de este método es que el paso de eliminación, que consume tiempo, se puede reformular de tal manera que involucre sólo operaciones sobre los elementos de la matriz de coeficientes, A De esta forma, es muy adecuado para aquellas situaciones donde se debe evaluar muchos vectores { B} El método de eliminación de Gauss puede implementarse como una descomposición LU La descomposición LU proporciona un medio eficaz para calcular la matriz inversa, la cual a su vez permite evaluar la condición de un sistema Partiendo de un sistema de ecuaciones lineales de la forma, •
•
•
•
•
A X B Este se puede ordenar como,
A X B 0 El primer paso de la eliminación de Gauss resulta en un sistema con una matriz triángular superior. u11 u12 u13 x1 d 1 0 u u23 x2 d 2 22 0 0 u33 x3 d 3
Que puede ser expresada como,
•
Ahora, suponga que existe una matriz triangular inferior con números 1 sobre la diagonal 1 L l 21 l 31
0 1
l 32
0
1 0
que tiene la siguiente propiedad
LU X D A X B
si esta propiedad se cumple, de las reglas de multiplicación de matrices se obtiene.
L U A L D B Estrategia para resolver el sistema 1. Paso de descomposición LU: la matriz [ A], se factoriza o descompone en matrices triangulares inferior [L] y superior [ U ] 2. Paso de sustitución: [L] y [U ] se usan par determinar una solución { X } para un vector {B}. Este paso consta de dos subpasos: Se determina el vector intermedio {D} resolviendo [L]{D}={B} por sustitución hacia delante, debido a que [ L] es una matriz triangular inferior Se determina { X } resolviendo [ U ]{ X }={D} por sustitución hacia atrás Descomposición LU con base en la eliminación de Gauss Partiendo de una matriz de coeficientes, se llega a una matriz triangular superior Para llegar a esta matriz [U ] •
•
Para llegar a esta matriz [U ] se multiplicó la fila 1 por el factor f 21 = a21/a11 y restando el resultado a la fila 2 se eliminó a21 se multiplicó la fila 1 por el factor f 31 = a31/a11 y restando el resultado a la fila 3 se eliminó a31 se multiplicó la fila 2 por el factor f 32 = a32’/a22’ y restando el resultado a la fila 3 se eliminó a32’ •
La matriz triangular inferior que tiene la propiedad requerida para la descomposición LU es 0 0 1 L f 21 1 0 f 31 f 32 1
Descomposición LU con base en la eliminación de Gauss 1. Paso de descomposición LU: la matriz [ A], se factoriza o descompone en matrices triangulares inferior [L] y superior [ U ] 2. Paso de sustitución: [L] y [U ] se usan par determinar una solución { X } para un vector {B}. Este paso consta de dos subpasos:
–
Se determina el vector intermedio {D} resolviendo [L]{D}={B} por sustitución hacia delante
d 1 b1 ;
d i bi
i 1
l b
para i 2, , n
ij j
j 1
–
xn
d n unn
Se determina { X } resolviendo [ U ]{ X }={D} por sustitución hacia atrás d i
;
xi
n
u d ij
j i 1
uii
j
para i n 1, n 1,,1
Matriz inversa •
•
•
•
•
Para una matriz cuadrada [ A], hay otra matriz [ A]-1 conocida como la inversa de [ A], para la cual se cumple, [A] [ A]-1 = [A]-1 [A] = I La matriz inversa se puede calcular en una forma de columna por columna a partir de vectores unitarios como vector de constantes del sistema de ecuaciones lineales algebraicas Por ejemplo, para determinar la primera columna de la matriz inversa se resuelve el sistema con el vector de constantes B=[1 0 0]T para determinar la segunda columna se usa B=[0 1 0]T y así sucesivamente La descomposición LU representa la mejor forma para implementar el cálculo de la matriz inversa, ya que una vez obtenida la descomposición LU de la matriz A se puede calcular su inversa resolviendo cada columna con los vectores unitarios como constantes Ejemplo: Determinar la inversa de
Análisis de error y condición del sistema La matriz inversa permite determinar si un sistema está mal condicionado, para esto existen 3 métodos: 1. Escalar la matriz de coeficientes [A ], de tal manera que el elemento más grande en cada fila sea 1. Si al invertir la matriz escalada existen elementos de la inversa [ A]-1 que sean varios ordenes de magnitud mayores que la unidad, es probable que el sistema esté mal condicionado 2. Multiplicar la inversa por la matriz de coeficientes original y verificar que [ A][ A]-1 I . Si no es así, indica que el sistema está mal condicionado 3. Invertir la matriz inversa y verificar que el resultado está lo suficientemente cercano a la matriz original. Si no es así, indica que el sistema está mal condicionado •
cond A A A1 •
• •
Este número mide la sensibilidad de la solución de un sistema de ecuaciones lineales a errores en los datos Valores cercanos a 1 indican que el sistema está bien condicionado Valores grandes indican que la matriz es casi singular
Número de condición de una matriz
• • •
Matrices banda Matrices simétricas Una matriz banda es una matriz cuadrada en la que todos sus elementos son cero, con excepción de una banda centrada sobre la diagonal principal
Las dimensiones de un sistema de banda se pueden cuantificar con dos parámetros: Ancho de banda, BW Ancho de media banda, HBW Estos se relacionan por BW = 2HBW + 1 En general un sistema de banda es aquel para el cual aij = 0 si |i – j | > HBW • •
Matrices especiales La eliminación de Gauss o la descomposición LU pueden emplearse para resolver sistemas de banda, pero si el pivoteo no es necesario resultan ineficientes, porque se utilizaría tiempo y espacio innecesario en el almacenamiento y manejo de ceros •
2. Regresión Polinómica: Dado un conjunto de datos (xi, yi) para i = 1, . . . , n, donde x y y representan las variables independiente y dependiente, respectivamente. Explique un procedimiento para calcular los parámetros de bi de los siguiente regresión polinomial.
La Regresión trata de buscar una función que permita explicar los valores de una variable en función de otra
Modelo de Regresión Siendo X la variable explicativa o independiente e Y la variable respuesta o dependiente, tendremos la Regresión Simple: Y = f(X) Si la variable respuesta, Y , depende de varias variables explicativas, X , tendremos la Regresión Múltiple: Y = f(X1 , X2 , . . . , Xn ) In statistics, regression analysis includes any techniques for modeling and analyzing several variables, when the focus is on the relationship between a dependent variable, Y , and one or more independent variables, X. Según la naturaleza de la función f podemos tener distintos tipos de Modelos de Regresión: Regresión Lineal Simple: Y=a+b·X Regresión Polinómica Simple: Regresión Lineal Múltiple:
Y=a+b·X+c·X
Y = a + b1· X1+ b2· X+ b3 · X3+ · · · + bn· Xn También hay Regresión: Logarítmica, Exponencial, Potencial, Hiperbólica, Trigonométrica, etc. Los valores desconocidos que caracterizan la función f se denominan Parámetros de Regresión, Regression Parameters. Qué criterio utilizar para escoger unos valores adecuados para los parámetros.
3. Regresión Múltiple: Un funcionario en una agencia gubernamental administra 5 nuevas pruebas de aptitud para cada uno de los 25 aspirantes para ciertas posiciones en la agencia. Para propósitos de estudio de los 25 aspirantes, los cuales fueron aceptados de acuerdo a su calificación. Después de período de prueba, cada aspirante fue valorado para competir por el trabajo. La calificación de las 5 pruebas (X1,X2,X3,X4,X5) y del score de competitividad (Y) fueron los siguientes: X1 X2 X3 X4 X5 Y ---------------------------------------------1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
86 62 110 101 100 78 120 105 112 120 87 133 140 84 106 109 104 150 98 120 74 96 104 94 91
110 97 107 117 101 85 77 122 119 89 81 120 121 113 102 129 83 118 125 94 121 114 73 121 129
100 99 103 93 95 95 80 116 106 105 90 113 96 98 109 102 100 107 108 95 91 114 93 115 97
95 99 101 91 102 94 89 112 110 87 90 101 100 85 99 101 93 108 100 96 95 91 80 85 105
87 100 103 95 88 84 74 102 105 97 88 108 89 78 109 108 102 110 95 90 85 103 80 104 83
88 80 96 76 80 73 58 116 104 99 64 126 94 71 111 109 100 127 99 82 67 109 78 115 83
Preguntas: a) Dado los datos anteriores calcular los parámetros, b, del siguiente modelo: Y = b0 + b1X1 + b2X2 + b3X3 + b4X4 + b5X5 Para el cálculo de los coeficientes de b, se implemento un programa en MatLab utilizando la multiplicación de matrices y la inversa
datos = load ('matrices.txt'); %Cargamos el archivo en la variable datos X = ones(25,6); %Llenamos de unos la variable X X(1:25,2:6)=datos(1:25,1:5); % Y = zeros(25,1); % Y(:)=datos(1:25,6); %La matriz Y llenamos los valores de la %matriz teniendo como filas del 1 al 25 %y como columna la 6. B = inv(X' * X) * (X'*Y)
B= -127.2370 0.2933 0.0372 1.3185 0.0474 0.5036 Puntuación del trabajo: Y = -127.2370 + 0.2933*X1 + 0.0372*X2 + 1.3185*X3 + 0.0474*X4 + 0.5036*X5
b) Cuál será la calificación de un nuevo aspirante, si los valores obtenidos en las pruebas son: (80, 90, 100, 110, 120) para X1 hasta X5 respectivamente? A (1,:) = [1]; A(2:6) = [80 90 100 110 120]; Y1 = B'*A' Y1 = 97.0694
c) Repita la pregunta anterior, si la relación entre la calificación de competitividad de trabajo Y y la calificación del test para X; esta dado por: Y = b0 + b1X1 + b2X2 + b3X3 + b4X4 + b5X5 + b6X11X2 + b7X1X3 + b8X1X4 + b9X1X5 clear close all clc datos = load ('matrices.txt');
%Cargamos el archivo en la variable datos
X=datos(:,1:5);
%Ponemos los datos de las columnas del 1 al 5 en X
Punos = ones(25,1); Y=datos(:,6);
%Ponemos una columna de unos en Punos %Ponemos la última columna de resultados en Y
X1=X(:,1); X6=X1.*X(:,2);
%Ponemos la primera columna en X %Ponemos la multiplicación de la primera columna por %la segunda columna de datos en X6 %Ponemos la multiplicación de la primera columna por %la tercera columna de datos en X7 %Ponemos la multiplicación de la primera columna por %la cuarta columna de datos en X8 %Ponemos la multiplicación de la primera columna por %la quinta columna de datos en X9
X7=X1.*X(:,3); X8=X1.*X(:,4); X9=X1.*X(:,5); nX=[Punos X X6 X7 X8 X9]; B=(inv(nX'*nX))*nX'*Y
%Armamos la matriz Nx con Punos, x1, X2, X3, X4, X5, %X6, X7, X8, X9 %Aplicamos la fórmula de cálculo de los coeficientes %de ecuaciones
X1
X2
X3
X4
X5
X1x2
X1x3
X1x4
X1x5
Y
1
86 110 100
95
87
9460
8600
8 170
7482
88
1
62
99 100
6014
6138
6138
6200
80
1 110 107 103 101 103 11770 11330 11110 11330
96
1 101 117
93
1 100 101
95 102
1
97
99
91
95
11817
9393
9191
9595
76
88 10100
9500 10200
8800
80
78
85
95
94
84
6630
7410
7332
6552
73
1 120
77
80
89
74
9240
9 600 10680
8880
58
1 105 122 116 112 102 12810 12180 11760 10710
116
1 112 119 106 110 105 13328 11872 12320 11760
104
1 120
89 105
87
97 10680 12600 10440 11640
99
1
81
90
88
7656
64
1 133 120 113 101 108 15960 15029 13433 14364
126
87
90
7047
7830
7830
1 140 121
96 100
89 16940 13440 14000 12460
94
1
98
78
6552
71
99 109 10812 11554 10494 11554
111
1 109 129 102 101 108 14061 11118 11009 11772
109
1 104
9672 10608
100
1 150 118 107 108 110 17700 16050 16200 16500
127
84 113
1 106 102 109
1
83 100
85
93 102
98 125 108 100
1 120
9492
8232
8632 10400
95 12250 10584
7140
9800
9310
99
94
95
96
90 11280 11400 11520 10800
82
1
74 121
91
95
85
1
96 114 114
1 104
73
93
1
94 121 115
1
91
129
97
6734
7030
6290
67
91 103 10944 10944
8736
9888
109
80
9672
8320
8320
78
85 104 11374 10810
7990
9776
115
9555
7553
83
105
80
83
8954
7592
11739
8827
Valores de los coeficientes: b0 = -238.9476221786654 b1 = 1.4102911447797 b2 = -0.4394416291384 b3 = 3.3907597939755 b4 = 1.7146905696797 b5 = -1.6673621662433 b6 = 0.0047034226729 b7 = -0.00.0207789758329 b8 = -0.0166652962838 b9 = 0.0218103251770
d) Discuta que podría pasar si el número de parámetros estimados son más grande que el numero de ejemplos. datos = load ('matrices.txt'); X = ones(25,6); X(1:25,2:6)=datos(1:25,1:5); X_N = X(1:25,:); i=25; while (i>2) XCU = X_N' * X_N; XAV = eig(XCU);
%Se obtiene los autovector
XAV_MIN = min(XAV); XAV_MAX = max(XAV);
%Se obtiene el mínimo autovalor %Se obtiene el máximo autovalor
FAC = XAV_MAX/XAV_MIN;
%Se obtiene la relación max/min de %autovalores
i = i - 1; X_N = X(1:i,:);
%Se escoge los nuevos valores para %X fprintf('I:%d FAC: %f\n',i, FAC); %Se imprime el FAC de cada %Iteración end I:24 I:23 I:22 I:21 I:20 I:19 I:18 I:17 I:16 I:15 I:14 I:13 I:12 I:11 I:10 I:9 I:8 I:7 I:6 I:5 I:4 I:3 I:2
FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC: FAC:
12584005.738765 12159397.250186 11780276.047055 12533095.584723 11944809.026838 12373783.803461 11808734.422008 11269030.036506 11464383.853942 10906782.625809 10187402.365430 9588712.605028 12359476.498200 11260243.821029 10950139.592433 10805641.607664 9713392.258705 12160359.020977 26838475.490894 55986551.891059 6546366786901714700000.000000 -16140180156686172000.000000 -53097578351227408.000000
De los datos obtenidos se puede indicar que el factor FAC = XAV_MAX/XAV_MIN nos indica la relación entre número de parámetros a estimar y el número de ecuaciones; cuando hay mas parámetros (incógnitas) que ecuaciones (datos) (como por ejemplo en las Iteraciones 4, 3, 2) la matriz es singular, es decir no tiene inversa, es inconsistente y no tiene solución. El número de ecuaciones debe ser igual o mayor que el número de variables de acuerdo a los valores obtenidos se puede indicar que a partir de la iteración 5 en adelante se obtienen soluciones estables y pueden tener matrices invertibles. En cualquier sistema de ecuaciones debe haber más ecuaciones que variables para que el modelo sea válido y sea consistente.
