=
=
Álgebra lineal numérica con Matlab =
Métodos Matemáticos de Especialidad (Mecánica-Máquinas) =
= Escuela Técnica Superior de Ingenieros Industriales Universidad Politécnica de Madrid
Javier García de Jalón de la Fuente Septiembre 2004
Álgebra lineal numérica con Matlab Métodos Matemáticos de Especialidad (Mecánica-Máquinas)
Escuela Técnica Superior de Ingenieros Industriales Universidad Politécnica de Madrid Javier García de Jalón de la Fuente Septiembre 2004
Índice
pág. i
ÁLGEBRA LINEAL NUMÉRICA CON MATLAB Índice 0. 1.
Prefacio.......................................................................................................................................v Introducción................................................................................................................................1 1.1 1.2
1.3 1.4 1.5 1.6
1.7
Tipos de matrices .......................................................................................................................... 1 Espacios vectoriales euclídeos ...................................................................................................... 1 1.2.1 1.2.2 1.2.3 1.2.4 1.2.5 1.2.6 1.2.7
Definición de producto escalar y de espacio vectorial euclídeo o hermítico............................... 1 Bases ortonormales en espacios euclídeos .................................................................................. 1 Coordenadas de un vector en una base ortonormal ..................................................................... 2 Existencia de bases ortonormales: método de Gram-Schmidt..................................................... 2 Interpretación gráfica del método de Gram-Schmidt .................................................................. 3 Matrices ortogonales ................................................................................................................... 4 Matrices de columnas ortogonales .............................................................................................. 4
1.6.1 1.6.2 1.6.3
Matrices de proyección y simetría en R 2 ..................................................................................... 7 Matriz de rotación en 2-D ........................................................................................................... 8 Matrices de rotación de Givens ................................................................................................... 8
Subespacios de una matriz A∈R m×n .............................................................................................. 4 Matrices de rango 1....................................................................................................................... 5 Dos formas de interpretar el producto de matrices ....................................................................... 6 Matrices de rotación, proyección y simetría ................................................................................. 7
Aproximación en norma cuadrática: Teorema de la proyección ortogonal .................................. 9 1.7.1 1.7.2 1.7.3 1.7.4 1.7.5 1.7.6
1.8
Normas de vectores y matrices ................................................................................................... 13 1.8.1 1.8.2 1.8.3 1.8.4 1.8.5 1.8.6 1.8.7
2.
Teorema de la proyección ortogonal ........................................................................................... 9 Matriz de proyección ortogonal sobre un subespacio................................................................ 10 Simetría ortogonal respecto de un subespacio........................................................................... 10 Matriz de Householder .............................................................................................................. 11 Aplicación de las matrices de Householder............................................................................... 11 Almacenamiento de las matrices de Householder:.................................................................... 13 Normas vectoriales en R n ó Cn .................................................................................................. 13 Norma de una matriz................................................................................................................. 15 Norma matricial inducida por una norma vectorial (norma natural):........................................ 15 Norma-1 matricial ..................................................................................................................... 16 Norma-∞ matricial .................................................................................................................... 16 Norma espectral......................................................................................................................... 17 Teoremas relacionados con las normas matriciales (sin demostración) .................................... 17
Sistemas de ecuaciones lineales ...............................................................................................19 2.1 2.2 2.3 2.4 2.5
Introducción a los sistemas de ecuaciones lineales..................................................................... 19 Interpretaciones del sistema Ax=b.............................................................................................. 19 Algunos casos posibles del sistema Ax=b en 2-D ...................................................................... 19 Sistema de m ecuaciones con n incógnitas.................................................................................. 20 El método de eliminación de Gauss básico................................................................................. 21 2.5.1 2.5.2 2.5.3 2.5.4 2.5.5
2.6
Método de eliminación de Gauss con pivotamiento ................................................................... 25 2.6.1 2.6.2 2.6.3
2.7
Necesidad del pivotamiento ...................................................................................................... 25 Método de Gauss con pivotamiento por columnas.................................................................... 26 Método de Gauss con pivotamiento total .................................................................................. 28
Resumen del método de eliminación de Gauss........................................................................... 29 2.7.1 2.7.2 2.7.3
2.8
Operaciones con filas y matrices elementales ........................................................................... 21 Factorización LU equivalente ................................................................................................... 22 Programa de Gauss básico n×n ................................................................................................. 23 Programa para resolver un sistema en la forma LUx=b ............................................................ 24 Programa de Gauss vectorizado ................................................................................................ 25
Reducción de una matriz a la forma de escalera........................................................................ 29 Conclusiones de la forma de escalera PA=LU.......................................................................... 29 La eliminación de Gauss y los cuatro subespacios de la matriz A ............................................ 30
Algunas funciones de Matlab en relación con el método de Gauss............................................ 30
Álgebra lineal numérica con Matlab
2.9
pág. ii
Errores en la resolución de sistemas de ecuaciones lineales....................................................... 31
2.9.1 2.9.2 2.9.3 2.9.4 2.9.5 2.9.6 2.9.7
Número de condición de una matriz.......................................................................................... 31 Casos particulares de propagación de errores en sistemas de ecuaciones lineales .................... 32 Fórmula general de propagación de errores en sistemas de ecuaciones lineales ....................... 32 Expresión general del número de condición.............................................................................. 33 Número de condición para la norma euclídea ........................................................................... 33 Inversas de matrices perturbadas............................................................................................... 34 Conclusiones sobre errores en sistemas de ecuaciones lineales ................................................ 35
2.10 Sistemas de ecuaciones lineales redundantes (m>n=r ) ............................................................... 35 2.10.1 2.10.2 2.10.3
Aplicación del teorema de la proyección ortogonal .................................................................. 35 Sistema de ecuaciones ampliado ............................................................................................... 36 Problema de mínimos cuadrados con matrices de columnas ortogonales ................................. 37
2.11 Sistemas de ecuaciones indeterminados (r =m
3.
Valores y vectores propios .......................................................................................................45 3.1 3.2 3.3
Definición del problema de valores y vectores propios.............................................................. 45 Interpretación geométrica de los valores y vectores propios ...................................................... 46 Propiedades de los valores y vectores propios............................................................................ 46 3.3.1 3.3.2 3.3.3 3.3.4
3.4
3.5
Definición y propiedades de las formas cuadráticas.................................................................. 52 Transformaciones de congruencia............................................................................................. 52 Ley de Inercia de Sylvester ....................................................................................................... 52 Matrices definidas-positivas...................................................................................................... 53 Matrices semi-definidas positivas e indefinidas........................................................................ 54
Cociente de Rayleigh .................................................................................................................. 55 3.7.1 3.7.2 3.7.3 3.7.4 3.7.5
3.8
Definición, casos particulares y propiedades ............................................................................ 49 Teorema espectral para matrices normales................................................................................ 50 Corolarios del teorema espectral ............................................................................................... 51 Descomposición espectral de matrices normales ...................................................................... 51
Formas cuadráticas y transformaciones de congruencia............................................................. 52 3.6.1 3.6.2 3.6.3 3.6.4 3.6.5
3.7
Matrices semejantes .................................................................................................................. 47 Diagonalización mediante transformaciones de semejanza....................................................... 47 Reducción a forma triangular mediante transformaciones de semejanza.................................. 48 Transformaciones de semejanza unitarias ................................................................................. 48 Lema de Schur........................................................................................................................... 49
Matrices normales....................................................................................................................... 49 3.5.1 3.5.2 3.5.3 3.5.4
3.6
Subespacios propios .................................................................................................................. 46 Relación de los valores propios con la traza y el determinante ................................................. 46 Propiedad de "desplazamiento" de los valores propios ............................................................. 46 Casos particulares del problema de valores propios.................................................................. 47
Transformaciones de semejanza ................................................................................................. 47 3.4.1 3.4.2 3.4.3 3.4.4 3.4.5
Definición y relación con los valores y vectores propios.......................................................... 55 Error en los valores propios estimados mediante el cociente de Rayleigh ................................ 55 Teorema mini-max y maxi-min (Teoremas de Courant-Fischer).............................................. 56 Interpretación geométrica de los teoremas mini-max y maxi-min ............................................ 57 Propiedad de "separación" de los valores propios..................................................................... 57
Valores y vectores propios generalizados................................................................................... 58 3.8.1 3.8.2 3.8.3 3.8.4 3.8.5
4.
Interpretación del sistema de ecuaciones general Ax=b............................................................ 42
Introducción a partir del problema estándar.............................................................................. 58 Planteamiento del problema generalizado de valores y vectores propios.................................. 58 Reducción del problema generalizado al problema estándar..................................................... 59 Cociente de Rayleigh para Ax=λ Bx.......................................................................................... 61 Convergencia a otros valores propios. Técnicas de deflacción matricial .................................. 61
Factorizaciones de una matriz ..................................................................................................63 4.1
Factorización LU ........................................................................................................................ 63 4.1.1 4.1.2 4.1.3 4.1.4
Ejemplo de factorización LU directa......................................................................................... 63 Fórmulas generales para la factorización LU............................................................................ 64 Factorización LU con matrices simétricas................................................................................. 64 Factorización LU vectorizada ................................................................................................... 66
Índice
pág. iii 4.1.5
4.2
Factorización QR ........................................................................................................................ 67 4.2.1 4.2.2 4.2.3 4.2.4
4.3 4.4
Factorización de Choleski ......................................................................................................... 67 Factorización QR por Gram-Schmidt........................................................................................ 67 Factorización QR de matrices rectangulares ............................................................................. 68 Factorización QR mediante matrices de Householder............................................................... 70 Factorización QR mediante rotaciones de Givens..................................................................... 71
Descomposición espectral de matrices normales........................................................................ 71 Descomposición de valores singulares (DVS)............................................................................ 72 4.4.1 4.4.2 4.4.3 4.4.4 4.4.5 4.4.6 4.4.7
Cálculo y existencia de la DVS................................................................................................. 73 Propiedades de las matrices U y V ............................................................................................ 74 Cálculo de la DVS..................................................................................................................... 74 Aplicaciones de la DVS ............................................................................................................ 75 Valores singulares y perturbaciones en una matriz ................................................................... 76 Problemas de reducción de la dimensionalidad......................................................................... 77 Aplicación de la DVS a la solución del sistema general Ax=b ................................................. 78
Prefacio
pág. v
0. Prefacio Los métodos numéricos para la resolución de problemas de Álgebra Lineal tienen una enorme im portancia práctica en la asignatura de Métodos Matemáticos de Especialidad, asignatura troncal que se imparte en el primer semestre del cuarto curso en la Escuela Técnica Superior de Ingenieros Industriales de la Universidad Politécnica de Madrid. Los métodos matriciales se aplican en multitud de áreas de la ingeniería, y además su resolución numérica es algo que se les da bastante bien a los computadores actuales. La experiencia del curso 2003-04 –primero en el que se impartió esta asignatura del Plan 2000– y de la asignatura Matemáticas de Especialidad del Plan 1976, aconseja no dar por supuesto que se mantienen los conocimientos adquiridos en las asignaturas de Álgebra I y II de primer curso. Tam poco es posible dedicar un gran número de horas de clase a presentar de nuevo los contenidos fundamentales de dichas asignaturas. El objetivo de estos apuntes es facilitar a loa alumnos de cuarto el repaso personal de los temas más importantes de cálculo vectorial y matricial, que vieron al comienzo de la carrera. La idea no es dar por sabida el Álgebra Lineal numérica en la asignatura de cuarto, sino liberarse de los aspectos más teóricos y poder dedicar las clases a los aspectos más algorítmicos y a las aplicaciones prácticas. Por tratarse de unos apuntes de repaso, el orden es diferente del que se utilizó o utilizaría en primero, conteniendo algunas referencias a temas que se tratarán posteriormente. Algunas demostraciones se han omitido, sobre todo cuando son complicadas y no tienen gran interés pedagógico. Muchas otras se han incluido, aunque pueden omitirse en una primera lectura. En muchos casos, se les ha restado importancia utilizando una letra más pequeña. Se han incluido tam bién algunos programas de Matlab, con objeto de ilustrar los algoritmos explicados en teoría. Estos algoritmos pueden servir de base para muchos otros, introduciendo en ellos las modificaciones oportunas. El autor de estos apuntes asume íntegramente la responsabilidad de sus limitaciones y deficiencias, pero no puede desaprovechar la oportunidad de agradecer a sus compañeros del Departamento de Matemática Aplicada a la Ingeniería Industrial, en particular a los profesores de Métodos Matemáticos de Especialidad y de Álgebra II, el apoyo y la ayuda prestadas en estas materias a lo largo de los últimos cursos. Madrid, 24 de septiembre de 2004,
Javier García de Jalón de la Fuente
Introducción
pág. 1
1. Introducción 1.1
Tipos de matrices
Una matriz A ∈ Cm×n es un conjunto de escalares –reales o complejos– dispuestos en m filas y n columnas. Si m = n la matriz se dice cuadrada y si m ≠ n la matriz se dice rectangular. Esta distinción es importante entre otras razones porque ciertas propiedades sólo tienen sentido si la matriz es cuadrada. Por ejemplo, sólo las matrices cuadradas tienen determinante, sólo ellas pueden ser invertibles, sólo ellas tienen valores y vectores propios, etc. Las matrices rectangulares tienen algunas propiedades que se pueden considerar como generalización de las anteriores. Así, los conceptos de matriz seudoinversa y de valores singulares, pueden ser considerados como una generalización de los conceptos de matriz inversa y de valores propios. Teniendo en cuenta los elementos distintos de cero, las matrices pueden ser diagonales, bidiagonales, triangulares superior o inferior, tridiagonales, de Hessenberg, etc.
1.2
Espacios vectoriales euclídeos
1.2.1 Definición de producto escalar y de espacio vectorial euclídeo o hermítico Un espacio vectorial euclídeo es un espacio vectorial real , de dimensión finita , dotado de un pro ducto escalar . Análogamente, un espacio vectorial hermítico es un espacio vectorial complejo, de dimensión finita, dotado de un producto escalar . Sea E un espacio vectorial euclídeo o hermítico. Se define el producto escalar en E como una aplicación: , : E × E → R con las tres propiedades siguientes: a ) Simetría:
u,v = v , u
b) Lineal en la 1ª variable:
α u+β v , w
c) Definida-positiva
= α u,w + β v, w u,u > 0, si u ≠ 0; u,u = 0, sólo si u = 0
(1)
La notación a indica el complejo conjugado de a . El producto escalar estándar (euclídeo en R n ó hermítico en Cn ), se expresa respectivamente de modo matricial en la forma: En R n : x, y = x T y = y T x = x1 y1 + x2 y2 + ... + xn yn
(2)
En Cn : x, y = x H y = x T y = 1 y1 + x2 y2 + ... + xn yn
(3)
El considerar el conjugado y transpuesto del primer factor en Cn es necesario para asegurar la pro piedad c) de (1). Los espacios euclídeos son un caso particular de los espacios herméticos en el que todos los vectores y escalares son reales. En lo sucesivo se considerarán sólo los espacios vectoriales euclídeos, salvo que se indique otra cosa. Prácticamente todo lo que se diga sobre ellos es válido para los espacios hermíticos cambiando la notación de transpuesta (T ) por la de conjugada y transpuesta ( H ).
1.2.2 Bases ortonormales en espacios euclídeos Se dice que un vector u ∈ R n está normalizado cuando cumple: 2
uT u = u 2 = 1
(4)
Álgebra lineal numérica con Matlab
pág. 2
Se dice que dos vectores u, v ∈ R n son ortogonales cuando su producto escalar es nulo:
uT v = 0
(5)
Se dice que una base ( q1 , q 2 ,..., q n ) en un espacio vectorial de dimensión n es ortonormal cuando está formada por n vectores de norma unidad y ortogonales entre sí :
qT iq j = δ ij
i, j = 1,2,..., n
(6)
1.2.3 Coordenadas de un vector en una base ortonormal Sea a un vector cualquiera en R n y ( q1 , q2 ,..., q n ) los vectores de una base ortonormal de R n . Se desean hallar las componentes α i de a expresado en dicha base ortonormal. El vector a se puede expresar como combinación lineal de los vectores qi con coeficientes por determinar α i :
a = α1q1 + α2q 2 + ... + α nq n
(7)
Para determinar la componente α i basta realizar el producto escalar por q i :
qTi a = α1qiT q 1 + α 2qiT q 2 + ... + αnq iT q n = 0 + ... + αiq iT qi + ... + 0 = αi
(8)
Despejando α i de esta ecuación y sustituyendo en (7):
a = ( q1T a ) q1 + (q T2 a ) q 2 + ... + (q T na ) q n
(9)
O bien, utilizando formulación matricial:
⎡α 1 ⎤ ⎢α ⎥ a = α1q1 + α2q2 + ... + α nq n = [q1 , q2 ,..., q n ] ⎢ 2 ⎥ =Qα ⎢⎥ ⎢α ⎥ ⎣ n⎦
⇒
α
= Q−1a
(10)
La matriz Q ≡ [q1 , q 2 ,..., q n ] es una matriz cuadrada y ortogonal cuyas columnas son los vectores de la base. Comparando las expresiones (9) y (10) se concluye que la inversa de Q es su transpuesta, lo que caracteriza a las matrices ortogonales:
⎡q1T a⎤ ⎢qT a⎥ T −1 α=⎢ 2 ⎥=Q a=Q a ⎢ ⎥ ⎢⎣qT 2 a⎥⎦
⇒
Q−1 = Q T
(11)
1.2.4 Existencia de bases ortonormales: método de Gram-Schmidt Se parte de n vectores ( a1 , a2 ,..., a n ) linealmente independientes y se trata de hallar a partir de ellos una base ortogonal ( q1 , q 2 ,..., q n ) . El primer vector q1 se puede obtener normalizando a1:
b1 ≡ a1
q1 = b1 b1
(12)
Introducción
pág. 3
Los vectores q2, q3, … se obtienen quitando a los vectores a2, a3, … las componentes según los " q" anteriores (calculadas según (9)) y normalizando (siempre existe bi ≠ 0 , por ser los ai linealmente independientes):
b 2 ≡ a 2 − ( q1T a2 ) q1
q 2 = b2 b2
(13)
b3 ≡ a3 − ( q1T a3 ) q1 − (q2T a3 ) q2
q3 = b3 b3
(14)
q n = bn b n
(15)
…
bn ≡ an − ( q1T a n ) q1 − − ( q T n a n ) q n
Existe otra forma de organizar los cálculos (método de Gram-Schmidt modificado) consistente en eliminar la componente según qi a todos los vectores a j ( j=i+1,...,n), tan pronto como qi ha sido calculado (se obtienen menores errores numéricos):
q1 = a1 a1
a j = a j − ( q1T a j ) q1
j = 2,3,..., n
q 2 = a 2 a2
a j = a j − ( q T 2 a j ) q 2
j = 3,4,..., n
a j = a j − ( q iT a j ) qi
j = i + 1,..., n
...
q i = ai ai
(16)
Obsérvese que en la expresión de cada qi (ai) sólo intervienen a1,..., ai (q1,..., qi). b3 = a3 − ( q1T a3 ) q1 − ( q2T a3 ) q 2
a3
q3 q1
q2
b 2 = a2 − ( q1T a2 ) q1
b1 = a1
a2 Figura 1. Interpretación gráfica del método de Gram-Schmidt.
1.2.5 Interpretación gráfica del método de Gram-Schmidt Se realizará para el espacio euclídeo tridimensional. Las expresiones (13)-(15) en este caso son: b1 ≡ a1
q1 = b1 b1
b2 ≡ a 2 − ( q1T a2 ) q1
q 2 = b2 b2
b3 ≡ a3 − ( q1T a3 ) q1 − (q2T a3 ) q2
q3 = b3 b3
(17)
Teniendo en cuenta que cada vector ai sólo tiene proyección no nula sobre los vectores ortonormales ( q1 , q2 ,..., qi ) , los subespacios generados por ( a1 , a2 ,..., ai ) y por ( q1 , q 2 ,..., q i ) coinciden:
Álgebra lineal numérica con Matlab
pág. 4
L [q1 ] = L [ a1 ] L [q1 , q 2 ] = L [ a1 , a2 ]
(18)
L [q1 , q 2 , q3 ] = L [ a1 , a2 , a3 ] La relación entre los vectores ( a1 , a2 , a3 ) , ( b1 , b2 , b3 ) y ( q1 , q2 , q3 ) se muestra en la Figura 1.
1.2.6 Matrices ortogonales Una matriz Q es ortogonal si es cuadrada y sus columnas q i son ortonormales entre sí:
( qiT qi = 1;
Q ∈ R n×n , qTi q j = δ ij
qiT q j = 0, i ≠ j )
(19)
Estas condiciones se expresan de modo gráfico y matricial tal como se muestra en la Figura 2 y en la ecuación (20). q1T
q
T 2
•
q1 q 2
qn
qT n
=
1 0 0 1
0 0
0 0
1
Figura 2. Ortogonalidad entre las columnas de la matriz Q.
QT Q = I
(20)
Propiedades de las matrices ortogonales (se pueden demostrar como ejercicio): T –1 − Se verifica que la inversa es igual a la traspuesta Q =Q , como se concluye de la ec. (20). −
Todos los valores propios de una matriz ortogonal tienen módulo unidad.
−
El determinante de una matriz ortogonal es +1 ó –1.
−
Las matrices ortogonales conservan el producto escalar, y por tanto distancias y ángulos. Casos particulares de las matrices ortogonales son las matrices de permutación, las de rotación y las de simetría.
−
1.2.7 Matrices de columnas ortogonales Se llama matriz de columnas ortogonales Q a una matriz rectangular m × n ( m > n ) , tal que sus columnas son un conjunto de vectores ortonormales. No son matrices ortogonales porque no son cuadradas. Las matrices rectangulares de columnas ortogonales, aunque cumplan QT Q=I, no son ortogonales pues QQT ≠I . En este caso QT no es la inversa de Q, sino una inversa por la izquierda. Se verá más adelante que QQT es una matriz de proyección sobre el subespacio de columnas de Q, Im(Q).
1.3
Subespacios de una matriz A R m×n
Un subespacio se puede determinar, entre otras formas, de los dos modos siguientes (representados respectivamente en la Figura 3 y en la Figura 4):
Introducción
pág. 5
1. Por un conjunto de vectores generador (ejemplo, un plano determinado por dos vectores no colineales). Estos vectores pueden constituir una base, o ser un sistema linealmente dependiente del que se puede extraer una base. 2. Por un conjunto de restricciones sobre un espacio vectorial que lo contiene (por ejemplo, un plano queda determinado por la recta a la que es perpendicular).
Figura 3. Base de un subespacio.
Figura 4. Complemento ortogonal.
En relación con una matriz rectangular A, de tamaño m×n y rango r , se pueden considerar los cuatro subespacios vectoriales fundamentales siguientes: 1. Subespacio de columnas Im( A). Es el subespacio de R m generado por las columnas de A. Tendrá dimensión r , pues hay r columnas independientes. Se verifica que Ax∈Im(A), ∀x. 2. Subespacio de filas Im( AT ). Es el subespacio de R n generado por las filas de A. También tendrá dimensión r , pues sólo hay r filas independientes. 3. Subespacio nulo Ker( A). Es subespacio de R n formado por todos los vectores x tales que Ax=0. Todos los vectores de Ker( A) son ortogonales a las filas de A (y a todas sus combinaciones lineales), por lo que Ker( A) e Im(AT ) son subespacios ortogonales y complementarios en R n. De aquí se concluye que la dimensión de Ker( A) es n – r. Se tiene: ⊥
Ker ( A ) = Im ( A T )
(21)
4. Subespacio nulo de la transpuesta Ker( AT ). Es el subespacio de R m formado por todos los vectores y que satisfacen yT A=0, es decir por todos los vectores ortogonales a las columnas de A (y a sus combinaciones lineales). Los subespacios Ker( AT ) e Im(A) son ortogonales y com plementarios, por lo que la dimensión de Ker( AT ) es m – r . En este caso, se tendrá: ⊥
Ker ( A T ) = Im ( A )
1.4
(22)
Matrices de rango 1
Las matrices de rango 1 se pueden obtener como producto de un vector columna por un vector fila: A=uvT . Las columnas de A son múltiplos de u y las filas de A son múltiplos de vT , por ejemplo:
⎡ u1 v T ⎤ ⎧ u1 ⎫ ⎡ u1v1 u1v2 u1v3 ⎤ ⎪ ⎪ A = uvT = ⎨u2 ⎬ {v1 v2 v3 } = ⎢⎢ u2 v1 u2 v2 u2v3 ⎥⎥ = [ uv1 uv2 uv3 ] = ⎢⎢ u2 v T ⎥⎥ (23) ⎪u ⎪ ⎢⎣ u3 vT ⎥⎦ ⎩ 3⎭ ⎣⎢ u3v1 u3v2 u3v3 ⎦⎥ Como sólo hay una fila y una columna linealmente independientes, el rango es uno . Como consecuencia de su rango, la matriz A=uvT tendrá un valor propio 0 con multiplicidad n−1, y un único valor propio distinto de cero. Es fácil determinar este valor propio y el vector propio asociado. La ecuación de valores y vectores propios de la matriz A es: (24) Ax = λ x
Álgebra lineal numérica con Matlab
pág. 6
En la ec. (24) el miembro izquierdo pertenece a Im( A) porque es una combinación de las columnas de A, luego para λ ≠ 0 también el miembro derecho, esto es x, pertenecerá a Im( A). Por ello, en el caso de la matriz de rango uno A=uvT , el vector propio asociado con el único valor propio distinto de cero será u, asociado con el valor propio vT u, ya que se cumple:
Au = ( uvT ) u = u ( v T u ) = ( v T u ) u = λu
⇒ λ = v T u
(25)
También se puede decir algo acerca de los n−1 vectores propios asociados con el valor propio nulo. Son todos los vectores del subespacio nulo o núcleo Ker( A), definido por la ecuación:
( uvT ) x = 0 ⇒
u ( vT x ) = 0
⇒
v T x = 0
(26)
El núcleo es pues el complemento ortogonal del vector v: ⊥
Ker ( A ) = L [ v ]
(27)
Una matriz de rango uno P = uvT , con vT u = 1 es una matriz de proyección, pues cumple:
P 2 = ( uv T )( uv T ) = u ( v T u ) v T = 1 ⋅ uv T = P
(28)
Sin embargo, no es una matriz de proyección ortogonal , pues no es simétrica. Se podría hablar de proyección oblicua . Si v = u , imponiendo la condición vT u = uT u = 1 se obtiene la expresión general de la matriz de proyección ortogonal sobre un subespacio de dimensión 1, que es:
uuT P = T = uuT uu
1.5
(29)
Dos formas de interpretar el producto de matrices
El producto de matrices C = AB se puede interpretar de varias formas diferentes, dos de las cuales están basadas respectivamente en el producto escalar o interior y en el producto exterior de vectores. Denotando respectivamente como a j y a jT a la columna j y a la fila j de la matriz A:
⎡ a1T ⎤ ⎢ aT ⎥ C = AB = ⎢ 2 ⎥ [b1 b2 ⎢⎥ ⎢ T⎥ ⎣⎢am ⎦⎥
C = AB = [a1 a2
⎡ a1T b1 a1T b2 ⎢ aT b a T b bn ] = ⎢ 2 1 2 2 ⎢ ⎢ T T ⎣⎢ amb1 a mb 2
a1T b n ⎤ ⎥ a2T b n ⎥ ⎥ ⎥ aTmb n ⎦⎥
⎡b1T ⎤ ⎢bT ⎥ l T T T 2⎥ ⎢ = a1b1 + a2b2 + ... + a l b l = ∑ aibiT al ] ⎢⎥ i =1 ⎢ T ⎥ ⎣⎢bl ⎦⎥
(30)
(31)
En la primera forma, el elemento (i, j) de la matriz producto es el producto escalar de la fila i de A por la columna j de B. En la segunda forma, el producto matricial se formula como una suma de l matrices de rango uno , formadas cada una por el producto de una columna de A por la correspondiente fila de B. Esta segunda forma es muy poco práctica desde el punto de vista del cálculo numérico del producto, pero puede arrojar luz sobre ciertos aspectos del producto de matrices. Por ejemplo, las columnas o las filas de la matriz producto se pueden expresar en la forma:
Introducción
pág. 7
⎡ c1T ⎤ ⎡ a1T ⎤ ⎢ cT ⎥ ⎢ aT ⎥ 2 ⎢ ⎥ = AB = ⎢ 2 ⎥ B ⎢⎥ ⎢⎥ ⎢ T⎥ ⎢ T ⎥ ⎢⎣c m ⎥⎦ ⎢⎣ am ⎥⎦
[c1 , c 2 ,..., c n ] = AB = A [b1 , b2 ,..., b n ] ,
1.6
(32)
Matrices de rotación, proyección y simetría
1.6.1 Matrices de proyección y simetría en R 2 Se comenzará considerando la forma de una matriz de proyección ortogonal P en 2-D. Se trata de determinar la proyección ortogonal de un punto b sobre el subespacio generado por el vector a.
b
b–p
a
p=Pb Sb Figura 5. Proyección y simetría sobre una recta en R 2.
Sea p el vector proyección de b sobre a. Se trata de hallar una matriz P tal que Pb=p. Si x es la magnitud de dicha proyección, la matriz de proyección P se puede determinar como sigue:
p = xa, ( b − p ) ⊥ a ⇒ aT ( b − xa ) =0 aT b aT b aaT aa T x = T ⇒ p = xa = T a = T b ≡ Pb P ≡ T aa aa aa aa
(33)
Esta expresión de P coincide con la ec. (29). Esta expresión de P es general, en el sentido de que sirve para R n cuando se proyecta sobre un subespacio de dimensión 1. Algunas propiedades de la matriz P son las siguientes: −
Es simétrica (sólo si la proyección es ortogonal). Se verifica que P2=PP=P (idempotente).
−
La matriz P tiene rango 1 (es singular) y su subespacio Im( A) está generado por el vector a.
−
El subespacio nulo Ker( A) es el complemento ortogonal a Im( A). − La proyección Pb pertenece a Im( A) (todas las columnas de P tienen la dirección de a). La Figura 5 muestra también el vector simétrico de b respecto de la recta a. La matriz de simetría S se puede determinar en función de la matriz de proyección P en la forma: −
Sb = b + 2 ( Pb − b ) = ( 2P − I ) b
⇒
S = 2P − I
(34)
La matriz de simetría es simétrica y ortogonal, es decir, su inversa es ella misma:
S = ST ,
ST S = S 2 = ( 2P − I )( 2P − I ) = 4P 2 − 4P + I = I
⇒
S −1 = S
(35)
Definiendo el vector a por medio de sus componentes c ≡ cos α y s ≡ sen α , la matriz de proyección P, correspondiente a la Figura 6, se puede calcular a partir de la expresión (33):
Álgebra lineal numérica con Matlab
pág. 8
⎧c ⎫ ⎨ ⎬ {c s} ⎡c 2 cs ⎤ aaT ⎩ s ⎭ = P = T = aa ⎧c ⎫ ⎢⎣ sc s 2 ⎥⎦ {c s} ⎨ ⎬ ⎩ s ⎭
(36)
Como era de esperar, la matriz P calculada es simétrica, singular e idempotente. Sustituyendo el resultado de (36) en (34) se obtiene una expresión para la matriz de simetría en R 2 mostrada en la Figura 7:
⎡ c 2 cs ⎤ ⎡1 0⎤ ⎡ c2 − s2 2 cs ⎤ S = 2P − I = 2 ⎢ −⎢ (37) ⎥ = ⎢ 2 sc s2 − c2 ⎥ 2⎥ 0 1 sc s ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ Más adelante se determinarán expresiones análogas y más generales para P y S en espacios R n.
b
b
b
a
a α
α
Pb
Figura 6. Proyección 2-D.
Sb
Figura 7. Simetría 2-D.
α
a a a x
Figura 8. Rotación 2-D.
1.6.2 Matriz de rotación en 2-D La Figura 8 muestra el resultado de rotar el vector a un ángulo α para convertirse en el vector b. En este caso, en función del coseno y el seno del ángulo girado, la matriz de rotación A es una matriz ortogonal dada por ( b=Aa; AT A= AAT =I):
⎧b x ⎫ ⎡c − s ⎤ ⎧ a x ⎫ ⎨b ⎬ = ⎢ ⎥⎨ ⎬ ⎩ y ⎭ ⎣ s c ⎦ ⎩a y ⎭
c ≡ cos α , s ≡ senα
(38)
1.6.3 Matrices de rotación de Givens Las matrices de rotación de Givens son rotaciones planas ó 2-D que se aplican a matrices n×n. Se utilizan para hacer cero un elemento determinado ( i, j) de una matriz A o como transformaciones de semejanza unitarias (como matrices de rotación, son ortogonales). Sea la matriz Gij una matriz n×n con todos sus elementos iguales a los de la matriz identidad I, excepto los cuatro elementos siguientes: gii = g jj = cos α ≡ c;
g ji = − gij = sen α ≡ s
(39)
Al pre-multiplicar la matriz A por Gij se combinan las filas i y j, y al post-multiplicar se combinan las columnas. En el primer caso se obtiene el siguiente resultado:
Introducción
⎡1 ⎢ ⎢ ⎢0 Gij A = ⎢⎢ ⎢0 ⎢ ⎢ ⎢⎣0
pág. 9
0
c
−s
s
c
0
0
0
⎡ a11 0⎤ ⎢ ⎥ ⎢ ⎥ ⎢ cai1 − sa j1 0⎥ ⎢ ⎥ A=⎢ ⎥ ⎢ sai1 + ca j1 0⎥ ⎢ ⎥ ⎥ ⎢ ⎢ 1⎥⎦ ⎢⎣ an1
a 1i
caii − sa ji
caij − sa jj
saii + ca ji
saij + ca jj
a nj
a1 j
a ni
a1n
⎤ ⎥ ⎥ cain − sa jn ⎥ ⎥ (40) ⎥ sain + ca jn ⎥ ⎥ ⎥ ⎥ ann ⎥⎦
Los valores de c y s se pueden determinar con la condición de que el nuevo elemento ( j,i) se anule: saii + ca ji = 0;
c=
−a ji aii s = , 2 2 2 2 aii + a ji aii + a ji
(41)
Las matrices de Givens se pueden utilizar para triangularizar una matriz A de modo semejante a como se hace en el método de eliminación de Gauss. Haciendo ceros por columnas, de izquierda a derecha, los ceros anteriores se conservan y la matriz A se triangulariza con un nº finito de operaciones aritméticas. Como se ha dicho anteriormente y es fácil de comprobar, las matrices de Givens son ortogonales y su producto, es decir la aplicación sucesiva de varias rotaciones de este tipo, también lo es. Por tanto, la factorización de A con matrices de Givens será del tipo A=QR (ver apartado 4.2). La factorización de Givens es más costosa en operaciones aritméticas que la descomposición LU.
1.7
Aproximación en norma cuadrática: Teorema de la proyección ortogonal
1.7.1 Teorema de la proyección ortogonal Enunciado: Sea E un espacio euclídeo (o hermítico), y F un subespacio de E . Sea v ∈ E un vector cualquiera de E . Las dos condiciones siguientes son equivalentes: w ∗ proyección ortogonal de b sobre F
b − w ∗ 2 = min b − w 2
⇔
w∈ F
(42)
La 2ª parte de la expresión (42) indica que el vector w*∈ F es una aproximación óptima en F del vector b∈ E . Además, dicha aproximación óptima es única. La Figura 9 ilustra este teorema. Demostración: Si w ∗ es la proyección ortogonal de b sobre F , se verificará que b − w∗ ∈ F ⊥ . Si w es un punto cualquiera de F , se podrá escribir: 2
2
b − w = ( b − w ∗ ) + ( w∗ − w ) =
= b−w
∗ 2
2
∗ 2
+ w∗ − w + 2 b − w∗ , w∗ − w = b − w
2
(43)
∗ 2
+ w∗ − w ≥ b − w
lo que demuestra que w ∗ es la aproximación óptima, que además es única porque para alcanzar la igualdad es necesario que se cumpla w − w∗ = 0 . Recíprocamente, la condición de aproximación óptima en (43) implica que w sea la proyección ortogonal de b sobre F .
Álgebra lineal numérica con Matlab
pág. 10
b
r∈ F ⊥
w
v w*=Pb
F
Sb Figura 9. Teorema de la proyección ortogonal.
1.7.2 Matriz de proyección ortogonal sobre un subespacio Sea F ⊂ E un subespacio vectorial determinado por las columnas linealmente independientes de una matriz A ∈ R m×n , que constituyen una base. Aplicando el teorema de la proyección ortogonal se puede calcular la matriz de proyección sobre dicho subespacio. Sea w ∗ = Ax la mejor aproximación del vector b en el subespacio Im( A), siendo x un vector de coeficientes a determinar. Se define el residuo r en la forma: r ≡ v − Ax
(44)
El teorema de la proyección ortogonal establece que r debe ser ortogonal a Im( A). Esta condición se puede imponer haciendo que r sea ortogonal a una base de Im( A), es decir, a las columnas de A:
A T r = 0 ⇒ A T ( v − Ax ) = 0
⇒ A T Ax = A T v
−1
⇒ x = (A T A ) A T v
(45)
El vector Ax es también la proyección ortogonal de v sobre el subespacio Im( A). Esta condición permite determinar la matriz de proyección P: −1
−1
Ax = Pv ⇒ A ( A T A ) A T v = Pv ⇒ P = A ( A T A ) A T
(46)
Obsérvese la analogía entre esta expresión y la de proyección sobre un subespacio de dimensión 1 vista en la ecuación (33), particularmente si ésta se transforma ligeramente: −1 aaT P = T = a ( aT a ) aT aa
(47)
1.7.3 Simetría ortogonal respecto de un subespacio Como en el caso 2-D, las matrices de proyección y simetría están estrechamente relacionadas. En este no se trata de proyectar un punto sobre un subespacio, sino de hallar su simétrico respecto a dicho subespacio. Se llamará Sb al punto simétrico de b respecto al subespacio Im( A), determinado por las columnas de A. La transformación de simetría puede ponerse en función de la proyección (ver Figura 9): Sb = b − 2r = b − 2 ( b − Pb ) = ( 2P − I ) b
⇒
S = ( 2P − I )
(48)
Introducción
pág. 11
Propiedades de la matriz de simetría S: 1. Es simétrica, pues tanto P como I lo son. 2. Es ortogonal , pues se verifica que su inversa es su transpuesta, es decir, ella misma:
SS = S 2 = ( 2P − I )( 2P − I ) = 4P 2 − 4P + I = I
(49)
A diferencia de la matriz de proyección, la matriz de simetría es invertible, como fácilmente se comprende a partir de su significado geométrico. Si se aplica dos veces se vuelve al origen.
1.7.4 Matriz de Householder La simetría o reflexión más importante se define respecto a un hiperplano (subespacio de dimensión n –1), determinado, no mediante una base (columnas de una matriz A), sino mediante un vector v perpendicular a dicho hiperplano, que genera su complemento ortogonal. Las matrices de simetría de Householder se denotan como H y se definen en la forma: vv T H = I − 2 T , u =1 o bien: H = I − 2uuT (50) v v donde v es un vector cualquiera. La matriz de Householder define una simetría respecto al hiperplano perpendicular al vector v. La Figura 9 justifica la expresión (50): Sea P1 la matriz de proyección sobre el subespacio de dimensión 1 generado por v. La matriz de proyección P sobre el hiperplano ortogonal a v y la correspondiente matriz de simetría de acuerdo con la expresión (48), serán: vvT vv T P = I − P1 = I − T ; H = 2P − I = 2 ( I − P1 ) − I = I − 2P1 = I − 2 T (51) v v v v Es obvio que las matrices de Householder son simétricas, y es fácil demostrar que son ortogonales: HT H = ( I − 2uuT )( I − 2uu T ) = I − 4uu T + 4uu T uu T = I
Como son matrices ortogonales, conservan ángulos y distancias, y por eso son muy estables en su aplicación numérica.
1.7.5 Aplicación de las matrices de Householder Las matrices de Householder se utilizan para, al multiplicarlas por un vector cualquiera x, anular todos los elementos excepto el primero, es decir, para alinear x con el vector e1 de la base natural. Supóngase un vector cualquiera x. Si el vector v se define en la forma: v = x − x e1 ,
e1T ≡ {1 0
0}
(52)
se demostrará que el producto Hx tiene nulos todos los elementos excepto el 1º. Como H es ortogonal conserva las magnitudes, luego se deberá cumplir:
Hx = x e1
(53)
En efecto, con H definida en la ecuación (50) y v definido por la ecuación (52), el producto Hx es:
Álgebra lineal numérica con Matlab
pág. 12 T
2 ( x − x e1 ) x ⎛ vv T ⎞ vv T x Hx = ⎜ I − 2 T ⎟ x = x − 2 T = x − ( x − x e1 ) T T = v v v v ⎝ ⎠ ( x − x e1 ) ( x − x e1 ) 2x T x − 2 x x1 = x − ( x − x e1 ) T 2 = x − ( x − x e1 ) = x e 1 x x − 2 x x1 + x
(54)
En el cálculo del vector v según (52) deben evitarse sustracciones o restas numéricas que puedan disminuir la precisión. Así, cuando x1>0, el cálculo del elemento v1 se puede hacer mediante una expresión alternativa en la que todos los términos son positivos: 2
− ( 2 + ... + xn ) x12 − x v1 = x1 − x = = x1 + x x1 + x 2
2
(55)
Las matrices de Householder se pueden utilizar para hacer ceros debajo de la diagonal en una matriz A, en forma análoga al método de Gauss:
⎡σ 1 ∗ ∗ ∗⎤ ⎢ 0 ∗ ∗ ∗⎥ ⎫ vv T ⎢ ⎥ U1 = H1 = I − 2 T ; v = x − x e1 ⎪ (56) v v ⎬ U1 A = ⎢ 0 ∗ ∗ ∗⎥ ⎢ ⎥ ⎪ σ 1 = x ⎭ ⎢ ⎥ ⎢⎣ 0 ∗ ∗ ∗⎥⎦ La segunda trasformación es similar, pero con una matriz H de tamaño n−1 completada con la primera fila y columna de la matriz identidad: ⎡1 0 0 0 ⎤ ⎡σ 1 ∗ ∗ ∗⎤ ⎢0 ⎥ ⎢ 0 σ ∗ ∗⎥ 2 ⎢ ⎥ ⎢ ⎥ U 2 = ⎢0 U 2U1 A = ⎢ 0 0 ∗ ∗⎥ ⎥ , H2 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢⎣0 ⎥⎦ ⎢⎣ 0 0 ∗ ∗⎥⎦ La Figura 10 y la Figura 11 muestran gráficamente las transformaciones correspondientes.
x
A
Figura 10. Primera transformación de Householder.
x
(57)
A
Figura 11. Segunda transformación de Householder.
Cada transformación conserva los ceros de las anteriores. Esta transformación es más estable que la de Gauss, aunque más costosa en número de operaciones aritméticas. A continuación de incluye una función de Matlab para calcular el vector v y el factor β a partir de un vector cualquiera x. El vector v se normaliza haciendo v(1)=1, de modo que v se pueda almacenar en lugar de los ceros introducidos (el primer elemento vale siempre 1 y no necesita ser almacenado). Se utilizan las expresiones siguientes.
Introducción
pág. 13
v = x − x e1 e1T ≡ {1 0
Si x1 > 0
⎫ ⎫⎪ ⎪ − ( x22 + ... + xn2 ) ⎬ H = I − β vvT ; ⎬ 0}⎪⎭ v1 = ⎪ x1 + x ⎭
β =
2 ; vT v
v = v / v1
(58)
% cál cul o del vect or de Househol der function [v,beta] = householder(x)
n=l engt h( x); % módul o de x ( si n pr i mer a component e) s =x( 2: n) ' * x( 2: n) ; v=[1,x(2:n)']'; if s==0
% l os el ement os x( 2) , . . . x(n) son cer o bet a=0 else
% módul o de x mx=sqr t ( x( 1) ^2+s) ; if x(1)<=0
v( 1) =x( 1) - mx; else
v( 1) =- s/ ( x( 1) +mx) ; end
bet a=2*v( 1) ^2/ ( s+v( 1) ^2) ; v=v/ v(1); end
1.7.6 Almacenamiento de las matrices de Householder: No hace falta almacenar las matrices de Householder H, sino que basta almacenar el vector v. Si se hace v1=1, el vector v puede almacenarse en las posiciones de los ceros introducidos en la matriz A. El producto de una matriz de Householder H por una matriz A se puede hacer, sin tener que formar las matrices H ó U explícitamente y con menos operaciones aritméticas, de la forma siguiente: T vv T H = I − 2 T = I − β vv T , HA = A − β vv T A = A − β v ( A T v ) (59) v v Con la expresión (59) el producto HA no necesita O(n3) operaciones aritméticas como un producto de matrices estándar, sino que basta un nº de operaciones aritméticas de orden O(n2).
1.8
Normas de vectores y matrices
Las normas de vectores y matrices proporcionan una forma sencilla de "medir" el tamaño de un vector y/o una matriz por medio de un número real positivo , en la misma forma en la que se emplea el módulo de un vector en R n.
1.8.1 Normas vectoriales en R n ó Cn Como se ha dicho, el objetivo es disponer de una medida del tamaño de los vectores, con objeto de estudiar convergencias, relaciones de mayor o menor, proximidad a cero, etc. El producto escalar permite definir de modo natural una norma en Cn, que coincide con la longitud de un vector en R 2 y R 3. Esta norma se llama norma euclídea, y responde a la expresión: x2=
x, x =
n
2
∑ x j
(60)
j =1
De un modo más general (hay otros tipos de normas posibles), la norma de un vector x debe satisfacer las tres propiedades siguientes:
Álgebra lineal numérica con Matlab
pág. 14
1. x ≥ 0, y x = 0 si y sólo si x= 0
(61)
2. α x = α ⋅ x , para cualquier α ∈ K ( R o C)
(62)
3. x + y ≤ x + y (desigualdad triangular)
(63)
La norma euclídea no es la única posibilidad. Puede considerarse que es un caso particular de una norma más general, la llamada norma-p, que se define en la forma: 1
⎛ n p ⎞ p x p = ⎜⎜ ∑ x j ⎟⎟ , p ≥ 1 ⎝ j =1 ⎠ Dependiendo del valor de p, existen varios casos particulares muy utilizados de la norma-p: 1. Norma-1: Se toma p=1, n
x1=∑
j
= x1 + x2 + ... + xn
(64)
(65)
j =1
2. Norma-2 ó norma euclídea: Se toma p=2, 1
2 1 2⎞ ⎛ n 2 2 2 2 x 2 = ⎜ ∑ x j ⎟ = ( x1 + x2 + ... + xn ) ⎝ j=1 ⎠
(66)
3. Norma- ó norma máxima: Se toma p=∞,
⎛ n p ⎞ x ∞ = lim x p = lim ⎜⎜ ∑ x j ⎟⎟ p→∞ p→∞ ⎝ j =1 ⎠ x 1 ≤1
1
= max x j
b) p=2
Figura 12. Representación de los puntos
(67)
1≤ j ≤ n
x 2 ≤1
1
a) p=1
1 p
x ∞ ≤ 1
1
c) p=∞
x ≤ 1 con tres normas distintas.
En la Figura 12 se muestra el lugar geométrico de los puntos de R 2 que cumplen la condición x ≤ 1 , con cada una de las tres normas- p consideradas. Es inmediato ver que las normas-p definidas cumplen las condiciones 1 y 2 (ecs. (62) y (63)) de la definición de norma. También se demuestra fácilmente que se cumple la condición 3; por ejemplo, para p=1: n
n
n
j =1
j =1
j =1
x + y 1 = ∑ x j + y j ≤ ∑ x j + ∑ y j = x 1 + y 1 A continuación se enuncian dos teoremas (sin demostración):
(68)
Introducción
−
pág. 15
Teorema 1: Cada norma vectorial x p es una función continua de las componentes del vector x, es decir, de x1, x2, ..., xn,.
−
Teorema 2: Para cada dos normas vectoriales x i y x j , existen dos constantes positivas m y M tales que: m x i ≤ x j ≤ M x i , ∀x ∈ Cn
(69)
Este teorema establece una equivalencia de normas: por ejemplo, si con una norma una secuencia de vectores tiende a cero, también tenderá a cero con cualquier otra norma.
1.8.2 Norma de una matriz Una norma de una matriz es un número real positivo que "mide" el tamaño de una matriz, y permite por ejemplo estudiar cuando una matriz tiende a otra matriz o a la matriz nula. A la norma de una matriz se le exige las mismas tres condiciones que debía satisfacer la norma de un vector, y una condición adicional relacionada con el producto de matrices: 1. A ≥ 0, y A = 0 si y sólo si A= 0
(70)
2. α A = α ⋅ A , para cualquier α ∈ K
(71)
3. A + B ≤ A + B (desigualdad triangular)
(72)
4. AB ≤ A ⋅ B
(73)
Es posible definir normas de matrices directamente. Por ejemplo, en Álgebra II se vio que las matrices A∈R m×n forman un espacio vectorial euclídeo con un producto escalar definido como: m
n
A, B = traza ( A B ) = ∑∑ aij bij T
(74)
i =1 j =1
Este producto escalar induce una norma, la norma de Frobenius , que responde a la expresión: 1
2 ⎛ 2⎞ A F = traza ( A A ) = ⎜ ∑∑ aij ⎟ ⎝ i =1 j =1 ⎠ Sin embargo, las normas matriciales más utilizadas se definen por otro camino. T
m
n
(75)
1.8.3 Norma matricial inducida por una norma vectorial (norma natural): En este caso la norma matricial se va a definir a partir de una norma vectorial en la forma: A ≡ max x ≠0
Ax , o bien: A ≡ max Au = Ay u =1 x
( y = 1)
(76)
donde y es el vector concreto de R n que produce el máximo (y la igualdad). El máximo se alcanza siempre porque el conjunto de vectores u es cerrado y acotado. De esta definición de norma matricial se concluye que:
Ax ≤ A x Se cumplen las cuatro condiciones de la definición de norma matricial:
(77)
Álgebra lineal numérica con Matlab
pág. 16
1. Si x ≠ 0, A ≠ 0, y A = 0 si y sólo si A= 0 2. α A = max
α Ax
= max
x
x ≠0
α
x ≠0
(78)
Ax = α A , ∀α ∈ K ( R o C ) x
(79)
3. A + B = ( A + B ) y = Ay + By ≤ Ay + By ≤ A + B
(80)
4. AB = ( AB ) y = A ( By ) ≤ A ⋅ By ≤ A ⋅ B
(81)
Se van a estudiar a continuación las normas matriciales inducidas por las normas vectoriales x 1 ,
x 2 y x ∞ , que se denotarán con la notación correspondiente. 1.8.4 Norma-1 matricial La expresión de la norma-1 de una matriz A es la siguiente: A 1 = max ∑ j =1 a jk n
(82)
k
La norma-1 es el máximo de la suma de los valores absolutos de los elementos de cada columna. Demostración: Supóngase que la igualdad en la definición de la norma se alcanza para un vector y tal que: A 1 = Ay 1 , siendo y 1 = 1.
(83)
Desarrollando la expresión de la norma-1 del vector Ay: A 1 = Ay 1 = ∑ j =1 ∑ k =1 a jk yk ≤ ∑ j=1 ∑ k =1 a jk yk = ∑ k =1 yk n
n
n
)
n
∑
n
n j= 1
a jk ≤
(84)
≤ ∑ k =1 yk max ∑ j=1 a jm = y 1 max ∑ j=1 a jm = max ∑ j=1 a jm n
(
n
m
n
n
m
m
El desarrollo anterior demuestra que el máximo de las sumas de los valores absolutos de cada columna es un límite la norma-1, pero para que ésta sea verdaderamente la expresión de la norma-1 este límite debe poder ser alcanzado, al menos para un cierto vector x.
superior de
En efecto, supóngase que la columna que da el máximo es la columna K . Tomando como vector x el vector e K de la base natural ( e K 1 = 1) : Ae K 1 = ∑ j =1 ∑ k =1 a jk δ kK = ∑ j=1 a jK = A 1 n
n
(85)
n
1.8.5 Norma- matricial La expresión de la norma-∞ de una matriz A es la siguiente: A ∞ = max ∑ j =1 akj n
(86)
k
La norma-∞ es el máximo de la suma de los valores absolutos de los elementos de cada fila. Demostración: Supóngase que la igualdad en la definición de la norma se alcanza para un vectory tal que: A
∞
= Ay ∞ , siendo y ∞ = 1.
(87)
Desarrollando la expresión de la norma-∞ del vector Ay: A ∞ = Ay
= max ∑ k =1 a jk yk ≤ max ∑ k =1 a jk yk ≤ n
∞
n
j
j
≤ max yk ⋅ max ∑ k =1 a jk = y ∞ max ∑ k =1 a jk = max ∑ k =1 a jk n
k
j
n
j
(88)
n
j
Introducción
pág. 17
Al igual que en el caso anterior, se ha obtenido un límite superior para la norma-∞. Para que la expresión obtenida corresponda a dicha norma, hay que demostrar que dicho límite se puede realmente alcanzar para un cierto vectorx. En efecto, supóngase que la fila que da la máxima suma de valores absolutos es la fila J . Se toma un vector x tal que xk = a Jk a Jk , si aJk ≠ 0; xk = 0, si aJk = 0 . Sustituyendo en (88): Ax
( x ∞ = 1)
= max ∑ k =1 a jk xk = ∑ k =1 aJk xk = ∑ k =1 a Jk2 aJk = ∑ k =1 a Jk , n
∞
n
n
n
j
(89)
1.8.6 Norma espectral La expresión de la norma espectral o norma-2 de una matriz cualquiera A es la siguiente:
A 2 = ρ ( A H A )
(90)
donde ρ ( M ) es el radio espectral , esto es, el valor absoluto del máximo valor propio en valor absoluto. La expresión (90) indica que la norma- 2 es la raíz cuadrada positiva del máximo valor pro pio de la matriz A H A, que es una matriz hermítica y –al menos– semidefinida positiva, cuyos valores propios son por tanto reales y no negativos. Demostración: La matriz A H A tendrá una base ortonormal de vectores propios u1, u2, ..., un: A H Au j = λ ju j ⇒ λ j = u Hj A H Au j
(91)
El vector y R n que produce la igualdad en la definición de la norma inducida (76) se puede expresar como combinación lineal de los n vectores propios u j: y 2 = ∑ k =1α k 2 = 1
y = ∑ k =1α k u k ; n
(92)
n
2
Introduciendo ahora la definición de la norma-2 del vector (Ay) y operando:
( ∑ α u ) A A (∑ α u ) = )(∑ α λ u ) = ∑ ( α λ ) ≤ max λ (∑ α ) = max λ = ρ ( A A ) 2
n
2
A 2 = Ay 2 = y H A H Ay =
= (∑
n
α H j =1 j j
u
n
n
k k
k =1
k
j =1
H j
j
k =1
n
2
k =1
k
n
H
k
i
i
k =1
k
k
2
k
(93)
H
i
i
De nuevo éste es solamente un límite superior, pero dicho límite se puede alcanzar para x=un: Au n 2 = u Hn A H Au n = u nH u nλn = λ n = ρ ( A H A ) 2
(94)
1.8.7 Teoremas relacionados con las normas matriciales (sin demostración) −
Teorema 1: Cualquier norma matricial A es una función continua de los m×n elementos de la matriz A∈Cm×n.
−
Teorema 2: Para cada par de normas matriciales A i y A j existen unas constantes positivas m y M tales que: m A i ≤ A j ≤ M A i
(95)
Este Teorema permite hablar de "equivalencia de normas matriciales". −
Teorema 3: Para cualquier norma natural y cualquier matriz cuadrada A∈Cn×n se verifica: ρ ( A ) ≤
A
pues A = max Ax ≥ Au n = λnu n = λn u n = λn = ρ ( A ) x =1
(96)
Álgebra lineal numérica con Matlab
−
pág. 18
Teorema 4: Para cualquier matriz cuadrada A∈Cn×n y cualquier valor ε arbitrariamente pequeño, existe alguna norma natural A tal que: ρ ( A ) ≤
A ≤ ρ ( A ) + ε
(97)
El radio espectral de una matriz cuadrada no es una norma (salvo que la matriz A sea nor mal , en cuyo caso A 2 = ρ ( A ) ), pero puede ser utilizado como tal dada su "cercanía". −
Corolario: Para cualquier matriz cuadrada A: ρ ( A ) = inf A . −
Sistemas de ecuaciones lineales
pág. 19
2. Sistemas de ecuaciones lineales 2.1
Introducción a los sistemas de ecuaciones lineales
La resolución de sistemas de ecuaciones lineales es uno de los problemas matemáticos más importantes en ingeniería. Hasta la llegada de los computadores digitales (segunda mitad del s. XX) la capacidad de resolver sistemas de ecuaciones lineales estaba muy limitada, no por la dificultad conceptual del problema, sino por el gran número de operaciones aritméticas necesarias. Ahora se puede resolver con un PC un sistema 1000×1000 en menos de 1 seg. Con programas especiales que aprovechan la estructura de la matriz se pueden resolver con PCs, de forma rutinaria, sistemas de decenas ó cientos de miles de ecuaciones lineales. Muchos otros métodos matemáticos (cálculo de valores y vectores propios, integración de ecuaciones diferenciales, optimización, ...) se reducen a la resolución repetida de sistemas de ecuaciones lineales. La resolución de sistemas de ecuaciones lineales tiene además un importante valor didáctico para los métodos numéricos en general y para la programación de ordenadores.
2.2
Interpretaciones del sistema Ax=b
El sistema de ecuaciones lineales Ax=b admite al menos las dos interpretaciones siguientes: 1. Intersección de hiperplanos. La solución es el punto (o conjunto de puntos) que satisface las ecuaciones de todos los hiperplanos, es decir, su intersección. La intersección puede ser un punto, un subespacio de dimensión n−r , o el conjunto vacío. 2. Combinación lineal de vectores. El vector término independiente b es una combinación lineal de las columnas de A, cuyos coeficientes son los valores de x. Con notación de Matlab:
Ax = b;
A(:,1) x1 + A (:,2) x2 + ... + A (:, n ) xn = b
(98)
b A(:,1) x1
A(:,3) x3
A(:,1)
A(:,3) A(:,2) A(:,2) x2
Figura 13. Intersección de hiperplanos.
Figura 14. Combinación lineal de vectores columna.
La expresión (98) indica que para que el sistema de ecuaciones Ax=b tenga solución, es necesario y suficiente que el vector b pertenezca a Im( A), es decir, al subespacio generado por las columnas de A. La solución será única si hay una única forma de expresar b como combinación lineal de las columnas de A.
2.3
Algunos casos posibles del sistema Ax=b en 2-D
Como ejemplo ilustrativo de lo dicho en el apartado anterior, en la Figura 15 se muestran geométricamente algunos casos posibles de sistemas de ecuaciones en R 2, teniendo en cuenta los números de ecuaciones y de incógnitas, el rango de la matriz A y el vector b. La interpretación es inmediata.
Álgebra lineal numérica con Matlab
pág. 20
m=3, n=2, r =2 b∉Im(A)
m=3, n=2, r =2 b∈Im(A)
m=2, n=2,
m=1, n=2, r =1
m=2, n=2, r =1
m=2, n=2, r =1
b∈Im(A)
b∉Im(A)
Figura 15. Algunos casos posibles del sistema Ax=b en R 2.
2.4
Sistema de m ecuaciones con n incógnitas
Se parte de un sistema de ecuaciones lineales expresado en forma matricial del siguiente modo: a11 x1 + a12 x2 + ... + a1 n xn = b1 ⎫ a21 x1 + a22 x2 + ... + a2 n xn = b2 ⎪⎪ ⎬ ⇒ Ax=b ⇒ ... ⎪ am1 x1 + am 2 x2 + ... + amn xn = bm ⎪⎭
⎡ a11 a12 ⎢a ⎢ 21 a22 ⎢ ⎢a a ⎣ m1 m 2
a1n ⎤ ⎧ x1 ⎫ ⎧ b1 ⎫ a2 n ⎥⎥ ⎪⎪ x2 ⎪⎪ ⎪⎪ b2 ⎪⎪ ⎨ ⎬= ⎨ ⎬ ⎥ ⎪ ⎪ ⎪ ⎪ ⎥ amn ⎦ ⎪⎩ xn ⎪⎭ ⎪⎩bm ⎪⎭
(99)
Como se ha dicho anteriormente, para que el sistema tenga solución, el vector b debe ser combinación lineal de las columnas de A (en Matlab, A(:,i) representa la columna i de A):
Ax = b
⇒
A (:,1) x1 + A (:,2) x2 + ... + A (:, n ) xn = b
(100)
El método de eliminación de Gauss está basado en el hecho de que una ecuación cualquiera puede sustituirse por una combinación lineal de esa ecuación y de las demás, sin que varíe la solución del sistema. Así pues, el método de eliminación de Gauss: − Combina ecuaciones (filas de la matriz A y elementos del vector b) de forma que el sistema adopte una forma más sencilla: forma triangular superior o forma de escalera. −
−
Permite entender mejor las características del sistema a resolver: si tiene solución o no, si la solución es única, etc. Es equivalente a una factorización PA=LU, donde P es una matriz de permutación, L es una matriz triangular inferior m×m con unos en la diagonal, que contiene los factores por los que se han multiplicado las filas de los pivots y U es una matriz m×n que resulta de transformar A a la forma de escalera.
Sistemas de ecuaciones lineales
2.5
pág. 21
El método de eliminación de Gauss básico
2.5.1 Operaciones con filas y matrices elementales A continuación se describe el método de eliminación de Gauss con un ejemplo 4×4. En primer lugar se hace un cero en la posición (2,1) multiplicando el sistema por la matriz P21:
⎡ a11 ⎢a A = ⎢ 21 ⎢ a31 ⎢a ⎣ 41
a12 a22 a32 a42
a13 a23 a33 a43
a14 ⎤ a24 ⎥⎥ , a34 ⎥ a44 ⎥⎦
⎡ 1 ⎢ −m P21 = ⎢ 21 ⎢ 0 ⎢ 0 ⎣
0 1 0 0
0 0 1 0
0⎤ 0⎥⎥ , 0⎥ 1⎥⎦
⎡ a11 ⎢0 P21A = ⎢ ⎢a31 ⎢a ⎣ 41
a12 a22′ a32 a42
a13 a23′ a33 a43
a14 ⎤ a24′ ⎥⎥ a34 ⎥ a44 ⎥⎦
(101)
donde: m21 =
a21 a ′ = a21 − m21a11 = a21 − 21 a11 = 0; a2′ j = a2 j − m21 a1 j , j = 2,...,4 ; a21 a11 a11
(102)
Seguidamente se obtienen ceros en los restantes elementos de la 1ª columna pre-multiplicando por las matrices P31 y P41. Se llega a la situación siguiente:
⎡ 1 ⎢ 0 P31 = ⎢ ⎢ −m31 ⎢ 0 ⎣
0 1 0 0
⎡ a11 ⎢0 P41P31P21A = ⎢ ⎢0 ⎢0 ⎣
0 0 1 0
0⎤ 0⎥⎥ , 0⎥ 1 ⎥⎦ a12 ′ a22 ′ a32 ′ a42
⎡ 1 0 0 0⎤ ⎧ m = ai1 ⎢ 0 1 0 0⎥ i1 ⎥ , ⎨⎪ a11 P41 = ⎢ ⎢ 0 0 1 0⎥ ⎪ ⎢ −m 0 0 1⎥ ⎩ aij′ = aij − mi1a1 j ⎣ 41 ⎦ a13 a14 ⎤ ⎡ 1 0 0 0⎤ ⎢ −m 1 0 0⎥ a23′ a24′ ⎥⎥ ⎥ , P1 ≡ P41P31P21 = ⎢ 21 a33′ a34′ ⎥ ⎢ −m31 0 1 0⎥ ⎢ −m 0 0 1 ⎥ a43′ a44′ ⎥⎦ ⎣ 41 ⎦
(103)
Obsérvese la forma peculiar en que se multiplican las matrices elementales Pi1. Si la matriz A es simétrica, la submatriz de P1A en la que hay que seguir haciendo ceros también lo es. De forma análoga se hacen ceros, debajo de la diagonal, en las columnas 2 y 3, pre-multiplicando por unas matrices P2 y P3:
⎡a11 a12 ⎢ 0 a′ 22 P42P32P1A = ⎢ ⎢0 0 ⎢0 0 ⎣
a13 a23′ ′′ a33 ′′ a43
a14 ⎤ a24′ ⎥⎥ , a34′′ ⎥ a44′′ ⎥⎦
⎡1 0 ⎢0 1 P2 ≡ P42 P32 = ⎢ ⎢0 −m32 ⎢0 −m 42 ⎣
⎡ a11 a12 a13 ⎢ 0 a′ a′ 22 23 P3P2P1A = ⎢ ⎢ 0 0 a33′′ ⎢0 0 0 ⎣
a14 ⎤ a24′ ⎥⎥ , a34′′ ⎥ ′′′ ⎥⎦ a44
⎡1 ⎢0 P3 ≡ P43 = ⎢ ⎢0 ⎢ ⎣0
0 0 1 0 0 1 0 −m43
0 0 1 0
0⎤ 0⎥⎥ 0⎥ 1⎥⎦ 0⎤ 0⎥⎥ 0⎥ 1⎦⎥
(104)
(105)
Para llegar a este resultado se ha supuesto que los elementos que aparecen sobre la diagonal −los pivots: a11 , a22′ y a33′′ − son distintos de cero, pues aparecen en los denominadores de los factores m ji. Finalmente se llega a un sistema de ecuaciones lineales equivalente al original, que tiene la forma:
Álgebra lineal numérica con Matlab
pág. 22
( P3P2P1A ) x = P3P2 P1b
(106)
Este sistema es mucho más fácil de resolver que el original, por los ceros introducidos en la matriz. El sistema podría resolverse calculando x4 de la cuarta ecuación, x3 de la tercera, x2 de la segunda, y finalmente x1 de la primera ecuación.
2.5.2 Factorización LU equivalente Las matrices elementales Pi introducidas en las expresiones (103), (104) y (105) tienen propiedades especiales y operan de una forma particular. Por una parte, es fácil comprobar que el producto P3P2 no ofrece ninguna característica especial:
⎡1 ⎢0 P3P2 = ⎢ ⎢0 ⎢0 ⎣
0 0 1 0 0 1 0 − m43
0 ⎤ ⎡1 0 0⎥⎥ ⎢⎢0 1 0⎥ ⎢0 − m32 1⎥⎦ ⎢⎣0 − m42
0 0 1 0
0 ⎤ ⎡1 0 0 0⎥⎥ ⎢⎢0 1 0 = −m32 1 0 ⎥ ⎢0 1 ⎥⎦ ⎢⎣0 m43m32 − m42 −m43
0⎤ 0⎥⎥ 0⎥ 1⎥⎦
(107)
Sin embargo, la inversa de las matrices P j tiene una forma muy sencilla, ya que se obtiene simplemente cambiando el signo de los factores mij, como se comprueba a continuación para P2:
⎡1 0 ⎢0 1 −1 P2P2 = ⎢ ⎢0 −m32 ⎢0 − m ⎣ 42
0 0 1 0
0 ⎤ ⎡1 0 0⎥⎥ ⎢⎢0 1 0⎥ ⎢0 m32 1⎥⎦ ⎢⎣0 m42
0 0 1 0
0 ⎤ ⎡1 0 ⎥⎥ ⎢⎢0 = 0 ⎥ ⎢0 1 ⎥⎦ ⎢⎣0
0 1 0 0
0 0 1 0
0⎤ 0⎥⎥ 0⎥ 1 ⎥⎦
(108)
Ahora se comprobará que el producto de las inversas de P3P2 −en orden inverso− se realiza mediante simple "superposición":
⎡1 0 ⎢0 1 −1 −1 P2 P3 = ⎢ ⎢0 m32 ⎢0 m 42 ⎣
0 0 1 0
0⎤ ⎡1 0⎥⎥ ⎢⎢0 0⎥ ⎢0 1⎥⎦ ⎢⎣0
0 0 1 0 0 1 0 m43
0⎤ ⎡1 0 0 0⎥⎥ ⎢⎢0 1 0 = 0⎥ ⎢0 m32 1 1⎥⎦ ⎢⎣0 m42 m43
0⎤ 0⎥⎥ 0⎥ 1⎥⎦
(109)
Pre-multiplicando la ecuación (105) por las inversas de las matrices Pi y teniendo en cuenta la forma (109) que adopta dicho producto de inversas:
⎡ a11 a12 a13 ⎢ ′ a23′ a22 −1 −1 −1 ⎢ 0 A = P1 P2 P3 ⎢ 0 0 a33′′ ⎢0 0 0 ⎣
a14 ⎤ ⎡ 1 0 0 a24′ ⎥⎥ ⎢⎢ m21 1 0 = a34′′ ⎥ ⎢ m31 m32 1 ′′′ ⎥⎦ ⎢⎣ m41 m42 m43 a44
0⎤ ⎡a11 a12 a13 0⎥⎥ ⎢⎢ 0 a22′ a23′ 0⎥ ⎢ 0 0 a33′′ 1⎥⎦ ⎢⎣ 0 0 0
a14 ⎤ a24′ ⎥⎥ a34′′ ⎥ a44′′′ ⎥⎦
(110)
La conclusión es que el método de eliminación de Gauss equivale a descomponer la matriz A en el producto de una matriz triangular inferior L con "unos" en la diagonal, por una matriz triangular superior U:
A = LU
(111)
La factorización LU tiene una gran importancia en álgebra lineal numérica y se estudiará con detalle en el apartado 4.1.
Sistemas de ecuaciones lineales
k
pág. 23
j
k
k i
k 0
=
k
*
i * k
* =
i 0
Figura 16. Eliminación de Gauss.
Figura 17. Vuelta atrás.
2.5.3 Programa de Gauss básico n× n A continuación se presenta un programa básico para resolver sistemas de ecuaciones lineales con Matlab (los programas en C y Fortran serían muy similares). Se realizan dos tares sucesivas: en primer lugar la triangularización, basada hacer cero el elemento (i,k ) combinando las filas i y k : % t r i angul ar i zaci ón de l a mat r i z A % se hacen cer os en l as n- 1 pr i mer as col umnas f or k=1: n- 1 % se hacen cer os en l a col umna k f or i =k+1: n m=A( i , k)/ A( k, k); % a l a f i l a i s e r e st a l a f i l a k mul t i pl i c ada por m f or j =k+1: n A( i , j ) =A( i , j ) - m* A( k , j ) ; end % se t r ansf or ma del mi smo modo el t ér mi no i ndependi ent e b( i ) =b( i ) - m* b( k) ; end end
En una segunda fase se calculan las incógnitas mediante la vuelta atrás, con el siguiente proceso: % se cal cul a x( n) de l a úl t i ma ecuaci ón x(n)=b( n) / A( n, n) ; % se cal cul a x(k) de l a ecuaci ón k f or k=n- 1: - 1: 1 s=0; f or i =k+1: n s =s +A( k, i ) * x( i ) ; end x( k) =( b( k ) - s ) / A( k, k) ; end
El algoritmo básico puede modificarse si la matriz A es simétrica, haciendo que el número de operaciones se reduzca aproximadamente a la mitad (en negrita las dos sentencias modificadas): % t r i angul ar i zaci ón de l a mat r i z A ( si mét r i ca) % se hacen cer os en l as n- 1 pr i mer as col umnas f or k=1: n- 1 % se hacen cer os en l a col umna k f or i =k+1: n m=A(k,i)/A(k,k);
% se tiene en cuenta la simetría
% a l a f i l a i s e r e st a l a f i l a k mul t i pl i c ada por m % sól o se oper a por enci ma de l a di agonal for j=i:n
A( i , j ) =A( i , j ) - m* A( k , j ) ; end % se t r ansf or ma del mi smo modo el t ér mi no i ndependi ent e b( i ) =b( i ) - m* b( k) ; end end
Álgebra lineal numérica con Matlab
pág. 24
2.5.4 Programa para resolver un sistema en la forma LUx=b Los programas anteriores suponen que las operaciones sobre las filas en el sistema Ax=b se realizan simultáneamente en ambos miembros. En realidad, hay una opción más favorable, que permite se parar las operaciones sobre la matriz A y sobre el vector b. En una primera fase se realiza la triangularización de la matriz A sin operar con el vector b, pero almacenando los factores m que aparecen en la ec. (110) en las posiciones de los nuevos ceros de A, con objeto de disponer al final de las matrices L y U. Una vez hecha la factorización, el sistema LUx=b se transforma definiendo un nuevo vector y≡Ux, con lo que se tiene el sistema Ly=b, del que se puede despejar fácilmente el vector y porque L es triangular inferior. Conocido y, se despeja el vector x del sistema Ux=y teniendo en cuenta que U es triangular superior. A continuación se dan los programas de Matlab que realizan la triangularización y las dos resoluciones con matrices triangulares (ver Figura 18 y Figura 19). Por brevedad, se han reducido al mínimo los comentarios en el código. % f act ori zaci ón A=LU f or k=1: n- 1 f or i =k+1: n m=A( i , k)/ A( k, k); f or j =k+1: n A( i , j ) =A( i , j ) - m* A( k , j ) ; end % se al macena el f act or m en l a posi ci ón del cer o A( i , k) =m; end end % Resol uci ón del si st ema Ly=b y=zer os( n, 1) ; y(1) =b( 1) ; f or k=2: n s=0; f or j =1: k-1 s =s +A( k, j ) * y( j ) ; end y( k) =b( k) - s ; end % Resol uci ón del si st ema Ux=y x=zer os( n, 1) ; x(n)=y(n)/ A( n, n) ; f or k=n- 1: - 1: 1 s=0; f or j =k+1: n s =s +A( k, j ) * x( j ) ; end x( k) =( y( k ) - s ) / A( k, k) ; end
k
k 0
k
1
*
=
* k
k
*
i *k
* =
0 Figura 18. Resolución de Ly=b.
Figura 19. Resolución de Ux=y.
Sistemas de ecuaciones lineales
pág. 25
2.5.5 Programa de Gauss vectorizado Matlab permite realizar una vectorización parcial o total del algoritmo: la velocidad de cálculo con Matlab (sobre todo en versiones anteriores a Matlab 6.5) aumenta considerablemente reemplazando el for más interno por una única instrucción que trabaja con dos filas completas: f or k=1: n- 1 f or i =k+1: n m=A( i , k)/ A( k, k); A(i,k+1:n)=A(i,k+1:n)-m*A(k,k+1:n);
b( i ) =b( i ) - m* b( k) ; end end
La vectorización puede ser doble realizando las operaciones del método de Gauss en otro orden diferente, mostrado en la Figura 20. Obsérvese que en vez de hacer ceros por columnas, se hacen cero los elementos de la fila i que están delante de la diagonal. Los elementos mostrados en gris han alcanzado ya su valor definitivo. Para transformar la fila i (hacer ceros en ella) se crea el vector de factores por los que hay que multiplicar las filas anteriores a la i:
v(1: i − 1) = A (1: i − 1, i )'. / diag(( A (1: i − 1,1: i − 1))'
(112)
La fila i puede sufrir todas las transformaciones a la vez, resultando un método sensiblemente más rápido que los anteriores (con versiones antiguas de Matlab): i −1
A (i, i : n ) = A (i , i : n ) + ∑ A ( j , i : n ) ∗ v ( j ) = A (i , i : n ) + v (1: i − 1) ∗ A (1: i − 1, i : n ) j =1
(113)
A (i,1: i − 1) = v (1: i − 1) i Elementos que ya son 0
i Elementos a hacer 0
Elementos ya transformados Fila a transformar
Elementos con el valor inicial
Figura 20. Vectorización doble del método de Gauss.
2.6
Método de eliminación de Gauss con pivotamiento
2.6.1 Necesidad del pivotamiento En el proceso de eliminación de Gauss básico anteriormente explicado es necesario dividir por los elementos que aparecen en la diagonal, que se denominan " pivots" o "pivotes". En ocasiones puede aparecer un cero (o un elemento de valor muy pequeño) en la posición de un pívot: en ese caso el proceso no puede continuar si se utiliza el algoritmo previo. Si debajo del pívot hay elementos no nulos se pueden permutar las filas correspondientes y proseguir la eliminación (la solución del sistema no varía). Se considerará un sistema 4×5, en el que se ha encontrado un elemento nulo en la posición (2,2):
Álgebra lineal numérica con Matlab
pág. 26
⎡a11 a12 ⎢0 0 P1A = ⎢ ⎢ 0 a32′ ⎢ 0 a′ 42 ⎣
a13 ′ a23 a33′ ′ a43
a14 a24′ a34′ a44′
a15 ⎤ a25′ ⎥⎥ a35′ ⎥ a45′ ⎥⎦
(114)
La permutación de filas puede hacerse pre-multiplicando por una matriz P en la forma:
⎡1 ⎢0 PP1A = ⎢ ⎢0 ⎢0 ⎣
0 0 1 0
0 1 0 0
0⎤ ⎡ a11 a12 0⎥⎥ ⎢⎢ 0 0 ′ 0⎥ ⎢ 0 a32 ′ 1 ⎥⎦ ⎢⎣ 0 a42
a13 ′ a23 a33′ a43′
a14 a24′ a34′ a44′
a15 ⎤ ⎡a11 a12 a25′ ⎥⎥ ⎢⎢ 0 a32′ = a35′ ⎥ ⎢ 0 0 a45′ ⎥⎦ ⎢⎣ 0 a42′
a13 a33′ a23′ a43′
a14 a34′ a24′ a44′
a15 ⎤ a35′ ⎥⎥ a25′ ⎥ a45′ ⎥⎦
(115)
La matriz P puede construirse a partir de un vector p que contiene el orden de los pivots [1,3,2,4]. También puede suceder que algunas columnas carezcan de pívot válido. Por ejemplo, si al hacer ceros en la columna 2 los elementos correspondientes de la columna 3 se han hecho también cero:
⎡ a11 a12 a 13 ⎢ 0 a′ a ′ ¨ 22 23 P2P1A = ⎢ 0 ⎢0 0 ⎢0 0 0 ⎣
a14 a24′ a34′′ ′′ a44
a15 ⎤ a25′ ⎥⎥ , a35′′ ⎥ a45′′ ⎥⎦
⎡1 ⎢0 P44 = ⎢ ⎢0 ⎢0 ⎣
0 0 1 0 0 1 0 −m43
0⎤ 0⎥⎥ ≡P, 0⎥ 4 1 ⎥⎦
m43 =
′′ a44 ′′ a34
(116)
el proceso normal falla y hay dos posibilidades: 1. La eliminación puede proseguir en las columnas posteriores y llegar hasta un último pívot en la última fila ( x3 es una variable libre o independiente, que puede tomar un valor arbitrario).
⎡ a11 a12 a 13 a14 ⎢ 0 a′ a′ a′ ¨ 22 23 24 P4P2P1A = ⎢ ′′ 0 a34 ⎢0 0 ⎢0 0 0 0 ⎣
a15 ⎤ a25′ ⎥⎥ a35′′ ⎥ ′′′ ⎥⎦ a45
(117)
2. La eliminación se detiene porque debajo del último pívot encontrado todas las filas restantes se han hecho cero (una fila dependiente de las demás, y dos variables libres: x3 y x5).
⎡ a11 a12 a 13 a14 a15 ⎤ ⎢ 0 a′ a′ a′ a ′ ⎥ ¨ 22 23 24 25 ⎥ P4P2P1A = ⎢ ′′ a35′′ ⎥ 0 a34 ⎢0 0 ⎢0 0 0 0 0⎥ ⎣ ⎦
(118)
2.6.2 Método de Gauss con pivotamiento por columnas Continuando con lo indicado en la sección anterior, en el caso general el método de Gauss con pivotamiento por columnas procede del siguiente modo, ilustrado con un ejemplo en la Figura 21: 1. Se busca el mayor elemento en valor absoluto de la primera columna. Este elemento servirá como pívot para hacer ceros en dicha columna. 2. El segundo pívot es el máximo elemento en valor absoluto de la 2ª columna, sin contar el elemento de la fila del primer pívot. Con este pívot se hacen ceros en los elementos de la 2ª columna que no pertenecen a una fila de un pívot anterior.
Sistemas de ecuaciones lineales
pág. 27
3. De modo análogo se calculan los pivots en las restantes columnas y se hacen cero los elementos de cada columna, pero sólo en las filas en las que no han aparecido pivots. 4. En el sistema final es fácil hallar las incógnitas xn, ..., x1, partiendo de la fila en la que ha aparecido el último pívot, luego en la del penúltimo pívot, y así hasta la del primero. 5. En la práctica, no es necesario intercambiar filas : basta conocer las filas dónde han ido apareciendo los sucesivos pivots, almacenando en un vector la información correspondiente. ⎡0 ⎢0 ⎢ ⎢0 ⎢⊗ ⎢ ⎢⎣ 0
∗ ∗ ∗ ∗ ∗
∗ ∗ ∗ ∗ ∗
∗ ∗ ∗ ∗ ∗
∗⎤ ⎡ 0 ∗⎥⎥ ⎢⎢ 0 ∗⎥ ⇒ ⎢ 0 ∗⎥⎥ ⎢⎢⊗ ∗⎥⎦ ⎢⎣ 0
0 ∗ ∗ ∗⎤ 0 ∗ ∗ ∗⎥ ⎥
⎡0 ⎢0 ⎢ ⊗ ∗ ∗ ∗⎥ ⇒ ⎢ 0 ∗ ∗ ∗ ∗⎥⎥ ⎢⎢⊗ 0 ∗ ∗ ∗⎥⎦ ⎢⎣ 0
0 ⊗ ∗ ∗⎤ 0 0 ∗ ∗⎥ ⎥
⎡0 ⎢0 ⎢ ⊗ ∗ ∗ ∗⎥ ⇒ ⎢ 0 ∗ ∗ ∗ ∗⎥⎥ ⎢⎢⊗ 0 0 ∗ ∗⎥⎦ ⎢⎣ 0
⊗ 0 ⊗ ∗ ∗ ∗ 0 0 0 0
∗ 0 ∗ ∗ ⊗
∗⎤ ⎡ 0 ∗⎥⎥ ⎢⎢ 0 ∗⎥ ⇒ ⎢ 0 ∗⎥⎥ ⎢⎢⊗ ∗⎥⎦ ⎢⎣ 0
⊗ 0 ⊗ ∗ ∗ ∗ 0 0 0 0
∗ 0 ∗ ∗ ⊗
∗⎤ ⊗⎥⎥ ∗⎥ ∗ ⎥⎥ ∗ ⎦⎥
Figura 21. Método de eliminación de Gauss con pivotamiento por columnas.
Posibles dificultades del método de Gauss con el pivotamiento por columnas descrito: −
Si una ecuación del sistema (o una fila de la matriz) se multiplica por un número muy grande (por ejemplo, 1e10), la solución del sistema no varía, pero se puede alterar el orden en el que se eligen los pivots, y con dicho orden la precisión de la solución calculada.
−
Por ello es conveniente que todas las filas de la matriz sean vectores que tengan aproxima damente la misma norma , o bien que, al elegir los pivots, se comparen los elementos de la columna correspondiente, dividiendo cada uno de ellos por la norma de su fila. A este proceso se le denomina escalamiento o scaling por filas.
A continuación se presenta un programa de Matlab, llamado pivotcol1.m, que selecciona directamente los pivots en cada columna, sin dividir por la norma de la fila. % f i c her o pi vot c ol 1. m % Gauss y pi vot ami ent o por col umnas. Vect or i zaci ón si mpl e % se gener an al eat or i ament e l a mat r i z y el vect or de un t amaño dado n=5; A=r ound( r and( n, n) *20- 10) ; b=r ound( r and( n, 1) *20- 10) ; x=zeros( n, 1) ; % se cal cul a l a sol uci ón exact a con el oper ador ( \ ) xe=AA\ bb; % vect or que cont endr á el orden de l os pi vot s p=[ 1: n] ; for k=1:n-1
% buscar el pí vot en col umna k. No i nt er vi enen l as f i l as de pi vot s ant er i or es [ pk, i k]=max(abs( A( p( k: n) , k)) ) ; % se i nt ercambi an l os el ement os en p para i ndi car el orden de l os pi vot s t emp=p( k) ; p( k)=p( i k+k- 1) ; p( i k+k- 1) =t emp; % hacer cer os en col umna pí vot for i=p(k+1:n)
m=A( i , k) / A( p( k) , k) ; A( i , k: n) =A( i , k: n) m*A( p( k), k: n) ; b( i ) =b( i ) - m* b( p( k) ) ; %A( i , k) =m; end end
% vuel t a at r ás: or den i nver so de apar i ci ón de pi vot s x( n) =b( p( n) ) / A( p( n) , n) ; for i=n-1:-1:1
s =A( p( i ) , i +1: n) * x( i +1: n) ; x( i ) =( b( p( i ) ) - s ) / A( p( i ) , i ) ; end
...
Álgebra lineal numérica con Matlab
pág. 28
di s p( ' Sol uc i ón de Mat l ab' ) ; di s p( xe) ; di s p( ' Sol uci ón cal c ul ada' ) ; di s p( x) ; di sp( ' Ya he t er mi nado' )
2.6.3 Método de Gauss con pivotamiento total En muchas aplicaciones prácticas el pivotamiento por columnas descrito anteriormente es suficiente para obtener una solución correcta y con pequeños errores numéricos. Sin embargo, existen casos en los cuales es necesario realizar pivotamiento total , eligiendo el pívot no entre los elementos de una columna, sino entre todos los elementos de la matriz que no están en una fila y/o columna de los pivots anteriores. Una posible aplicación del pivotamiento total es en los sistemas indeterminados, cuando es necesario elegir las variables libres más adecuadas. El método de Gauss con pivotamiento total se ilustra en la Figura 22 y procede del siguiente modo: 1. En primer lugar se busca el mayor elemento en valor absoluto de toda la matriz . Este elemento será el primer pívot: en el caso de la Figura 22 se supone que es la fila 4, columna 2. 2. Se hacen ceros todos los elementos de la columna 2 (menos el pívot). 3. Se busca el máximo elemento en valor absoluto de toda la matriz, excluyendo la fila 4 y la columna 2. Este elemento será el 2º pívot (fila 3, columna 4). 4. Se hacen ceros en la columna 4 excepto en las filas de pivots anteriores. 5. Se prosigue de la misma forma hasta que todas las filas y/o columnas tienen un pívot. 6. En la vuelta atrás, las incógnitas se calculan en orden decreciente de columnas en las que han aparecido los pivots, y cada una de ellas se calcula en el orden decreciente de filas en las que han aparecido pivots. ⎡∗ ⎢∗ ⎢ ⎢∗ ⎢∗ ⎢ ⎢⎣∗
0 ∗ ∗ ∗⎤ ⎡∗ 0 ∗ 0 ∗⎤ ⎡ 0 0 ∗ 0 ∗⎤ ⎡ 0 0 ∗ 0 ⊗ ⎤ ⎡ 0 0 0 ∗ ∗ ∗⎥⎥ ⎢⎢∗ 0 ∗ 0 ∗⎥⎥ ⎢⎢ ⊗ 0 ∗ 0 ∗⎥⎥ ⎢⎢ ⊗ 0 ∗ 0 ∗ ⎥⎥ ⎢⎢⊗ 0 0 ∗ ∗ ∗⎥ ⇒ ⎢∗ 0 ∗ ⊗ ∗⎥ ⇒ ⎢ ∗ 0 ∗ ⊗ ∗⎥ ⇒ ⎢ ∗ 0 ∗ ⊗ ∗ ⎥ ⇒ ⎢ ∗ 0
⊗ ∗ ∗ ∗⎥⎥ 0 ∗ ∗ ∗⎥⎦
⎢∗ ⊗ ∗ ∗ ∗⎥ ⎢ ⎥ ⎢⎣∗ 0 ∗ 0 ∗⎥⎦
⎢ ∗ ⊗ ∗ ∗ ∗⎥ ⎢ ⎥ ⎢⎣ 0 0 ∗ 0 ∗⎥⎦
⎢ ∗ ⊗ ∗ ∗ ∗⎥ ⎢ ⎥ ⎢⎣ 0 0 ∗ 0 0 ⎥⎦
∗ ∗ ∗ ⎢∗ ⊗ ∗ ⎢ ⎢⎣ 0 0 ⊗
⊗⎤ ∗ ⎥⎥ ⊗ ∗⎥ ∗ ∗⎥⎥ 0 0 ⎦⎥ 0 0
Figura 22. Método de eliminación de Gauss con pivotamiento total.
El pivotamiento total mejora algo la estabilidad del proceso de eliminación de Gauss, pero en general se considera que esta ganancia no merece la pena por el esfuerzo adicional que lleva consigo. Como ya se ha dicho, se utiliza en ocasiones cuando el objetivo es determinar el rango de la matriz o una partición óptima entre variables básicas y variables libres . Con pivotamiento total puede ser conveniente utilizar escalamiento por filas y por columnas, según la siguiente expresión ( D1 y D2 permutaciones de una matriz diagonal):
D1AD 2 y = D1b, siendo x = D 2y
(119)
El escalamiento por columnas implica un cambio de unidades en las incógnitas: en algunos casos puede ser conveniente no mezclar metros con micras, y en otros casos puede ser necesario hacerlo. Hay que señalar que no hay criterios absolutos para hacer un buen escalamiento por columnas. Un buen criterio podría ser hacer que el máximo elemento de cada fila y columna tenga una magnitud parecida, por ejemplo que esté entre ½ y 1. Hay que tener en cuenta que si al hacer el scaling se multiplica y divide siempre por potencias de 2, no se introducen errores de redondeo adicionales.
Sistemas de ecuaciones lineales
2.7
pág. 29
Resumen del método de eliminación de Gauss
2.7.1 Reducción de una matriz a la forma de escalera Mediante el método de eliminación de Gauss con pivotamiento por columnas y permutación de filas un sistema Ax=b, de tamaño m×n, se puede reducir a la forma de "escalera" siguiente:
⎡⊗ ⎢0 ⎢ PA = LU = L ⎢ 0 ⎢0 ⎢ ⎢⎣ 0
∗ ⊗ 0 0 0
∗ ∗ 0 0 0
∗ ∗ ⊗ 0 0
∗ ∗ ∗ 0 0
∗ ∗ ∗ 0 0
∗ ∗ ∗ 0 0
∗ ∗ ∗ 0 0
∗⎤ ∗ ⎥⎥ ∗⎥ ⊗⎥⎥ 0 ⎥⎦
⊗ ∗ P L
→ pivot → elemento ≠ 0 (120) → matriz de permutación m × m → matriz triangular inferior m × m
Esta forma de escalera es muy interesante, pues arroja mucha más luz sobre el sistema de ecuaciones a resolver que el sistema original Ax=b. Aparecen r pivots (elementos distintos de cero), con ceros debajo y a la izquierda. El rango de la matriz A es r . En la práctica las cosas no estarán tan claras, por ejemplo cuando los elementos (3,3), (4,5), (4,6), (4,7) y (4,8) sean muy pequeños, pero no exactamente cero: ¿Cuándo un número pequeño debe pasar a considerarse cero? La mejor solución a este tipo de preguntas no está en la forma de escalera, sino el la DVS (Descomposición de Valores Singulares), que se verá en el apartado 4.2.4. De momento, se supondrá que la forma de escalera no presenta confusión posible entre elementos cero y elementos muy pequeños.
2.7.2 Conclusiones de la forma de escalera PA=LU A partir de la forma de escalera dada por la expresión (120) se pueden extraer las conclusiones siguientes: − La matriz P indica las posibles permutaciones de filas realizadas para que los pivots aparezcan siempre ordenadamente en las filas superiores. La factorización LU realizada corres ponde a la matriz PA, es decir a la matriz A con las filas permutadas. − En el sistema original sólo r de las m ecuaciones eran independientes y por tanto sólo r de los m vectores fila (en R n) eran independientes. Por eso hay m−r filas que se han hecho cero, correspondientes a las filas que eran linealmente dependientes de las demás. − Las columnas (y las incógnitas) aparecen divididas en dos grupos: las que tienen pivots (variables básicas o dependientes) y las que no tienen pivots (variables libres o independien tes). Si hay variables libres, la solución, caso de existir, no será única. Dependiendo de lo que haya sucedido con el vector b (que habrá debido sufrir las mismas permutaciones y combinaciones de filas que la matriz A), se pueden presentar distintas posibilidades. Recuérdese que, matricialmente, el sistema se ha transformado en la forma: PAx = LUx = Pb
⇒ Ux = L−1Pb = ( P4P3P2P1 ) Pb
(121)
Existe solución si las (m – r ) filas que se han hecho cero en el primer miembro Ux se hacen también cero en el segundo miembro L –1Pb, es decir cuando los (m – r ) últimos elementos de ( P4 P3 P2 P1 ) Pb se anulan al igual que las correspondientes filas de U. Además, si el sistema tiene solución y hay variables libres, el sistema tiene infinitas soluciones.
Álgebra lineal numérica con Matlab
pág. 30
2.7.3 La eliminación de Gauss y los cuatro subespacios de la matriz A El sistema de ecuaciones lineales Ax=b transformado a forma de escalera según las ecuaciones (120) y (121), permite relacionar los cuatro subespacios fundamentales de las matrices A y U. 2.7.3.1 SUBESPACIO DE FILAS IM(AT )⊂R N Supóngase que el rango es r ≤m. El subespacio de filas Im( AT ) coincide con Im(UT ) por la forma en la que se ha construido U mediante combinaciones lineales de las filas de A. Las filas no nulas de U constituyen una base de Im( A). Si r =m el sistema siempre tiene solución. Si r
Ax = b ⇒
⎡⊗ ⎢0 −1 PL Ax = Ux = ⎢ ⎢0 ⎢0 ⎣
∗ ⊗ 0 0
∗ ∗ 0 0
∗ ∗ ⊗ 0
∗ ∗ ∗ 0
⎧∗⎫ ∗⎤ ⎪ ⎪ ∗⎥⎥ ⎪∗⎪ −1 x = PL b = ⎨ ⎬ ∗⎥ ⎪∗⎪ 0⎥⎦ ⎩⎪ ? ⎭⎪
(122)
2.7.3.2 SUBESPACIO NULO K ER (A)⊂R N El núcleo de A, Ker(A), está formado por los vectores x∈R n que satisfacen Ax=0. Por tanto, Ker(A) es el complemento ortogonal de Im( AT ) y su rango es n – r . Se verifica que Ker ( U ) = Ker ( A ) . Se puede obtener una base de ambos subespacios resolviendo el sistema de ecuaciones n – r veces, dando alternativamente valor unidad a una de las variables libres y cero a las demás. Los vectores de Ker( A) representan las relaciones de dependencia lineal que existen entre las columnas de A, como se ve de la ecuación Ax=0. 2.7.3.3 SUBESPACIO DE COLUMNAS IM(A)⊂R M La dimensión de Im( A) es r , rango de A. Si r =m el sistema tiene solución para todo b. Si r
2.8
Algunas funciones de Matlab en relación con el método de Gauss
En esta apartado se van a describir brevemente algunas funciones de Matlab en relación con el método de eliminación de Gauss y la resolución de sistemas de ecuaciones lineales. 1. Operador barra invertida (\) Por lo general, la operación A\b produce el mismo resultado que inv(A)*b, con algunas diferencias, entre las que está la forma de realizar los cálculos:
Sistemas de ecuaciones lineales −
pág. 31
Si A es una matriz cuadrada no singular, A\b es la solución de Ax=b calculada por eliminación de Gauss con pivotamiento por columnas.
Si A y b son tales que el sistema tiene infinitas soluciones o ninguna, A\b da la solución de variables libres nulas o la de mínimos cuadrados, respectivamente. − Además, el operador \ es "inteligente", en el sentido de que es capaz de descubrir y aprovechar si A es simétrica y definida-positiva, o si tiene forma triangular. 2. Función [L,U,P]=lu(A) − Calcula con pivotamiento por columnas la factorización LU de la matriz PA, donde P es una matriz de permutación que representa el orden de filas en que han aparecido los pivots. 3. Función [R,c]=rref(A,tol) − Reduce A a forma de escalera con pivots unidad y ceros también encima de cada pívot. El vector c es un vector de enteros cuyo número de elementos indica el rango de la matriz y cuyos elementos indican las columnas de los pivots (variables básicas o dependientes). −
2.9
Errores en la resolución de sistemas de ecuaciones lineales
2.9.1 Número de condición de una matriz A continuación se incluye un recordatorio de ciertos conceptos estudiados en Álgebra, útiles para este apartado y los siguientes. Se llama matriz convergente a una matriz que cumple las condiciones: lim A m = 0 ⇔ ρ ( A ) < 1 ⇔ lim A m = 0
m→∞
(123)
m→∞
donde ρ ( A ) es el radio espectral de A y la serie I + A + A 2 + A 3 + ... converge si A es convergente. En efecto, si A es convergente, la matriz ( I – A) es no singular y se puede expresar como:
S ≡ I + A + A 2 + A 3 + ... = I + A ( I + A + A 2 + A 3 + ... ) = I + AS S (I − A ) = I Proposición: Si
−1
⇒
2
3
( I − A ) = I + A + A + A + ...
(124)
A < 1 , la matriz ( I − A ) es no singular y se verifican las desigualdades: 1 1 −1 ≤a ) ( I − A ) ≤b) 1+ A 1− A
(125)
Demostración de la desigualdad a) I = (I − A ) (I − A ) 1 ≤ I − A ⋅ (I − A)
−1
≤ (1 + A ) ⋅ ( I − A )
−1
−1
1 −1 ⇒ ≤ (I − A) 1+ A
(126)
Demostración de la desigualdad b) −1
−1
−1
−1
−1
(I − A) (I − A ) = I ⇒ (I − A) − (I − A ) A = I ⇒ (I − A ) = I + (I − A ) A (I − A)
−1
≤ 1 + A ⋅ (I − A)
−1
⇒
−1
(I − A )
1 ≤ 1− A
(127)
Álgebra lineal numérica con Matlab
pág. 32
2.9.2 Casos particulares de propagación de errores en sistemas de ecuaciones lineales Se estudiará en primer lugar el error relativo en la solución debido exclusivamente a una perturbación o error en el término independiente : b ≤ A x
Ax = b Ax = b Ax = b ⎫ ⎫ ⎫ ⇒ ⇒ ⎬ ⎬ ⎬ ⇒ A ( x + δ x ) = b + δ b ⎭ Aδ x = δ b ⎭ δ x = A −1δ b ⎭
δx
≤ A −1
⎫⎪ ⎬ δ b ⎪ ⎭
(128)
Multiplicando miembro a miembro las desigualdades anteriores se puede calcular una acotación para el error relativo de la solución, en función del error relativo del término independiente: b δ x ≤ A −1 A x δ b
⇒
δx
x
≤ A −1 A
δ b
b
(129)
A continuación se estudiará el error relativo en la solución debido a una perturbación o error en la matriz del sistema: Ax = b ⎫ ⇒ Aδ x + δ A ( x + δ x ) =0 ⇒ δ x = −A −1δ A (x + δ x ) ⎬ ( A + δ A )( x + δ x ) = b ⎭
(130)
Tomando normas y operando se llega al siguiente resultado: ≤ A −1
δx
δA
x +δx
⇒
δx
x +δx
≤ A −1
δA
⇒
δx
x + δ x
≤ A− 1 A
δ A
A
(131)
Tanto en la expresión (129) como en la (131) el error relativo en los datos se multiplica por un mismo factor: κ (A) ≡
A −1 A ,
κ (A ) ≥ 1
pues I = A −1A, 1 ≤ A −1 A = κ ( A )
(132)
Este factor κ ( ) se llama número de condición o condición numérica y es un número real mayor o igual a uno que controla la transmisión de errores en los sistemas de ecuaciones lineales. Tal y como se ha definido, se exige que la matriz A sea no singular. Otra observación importante es que las cotas para el error relativo (129) y (131) son óptimas, en el sentido de que la igualdad es alcanzable.
2.9.3 Fórmula general de propagación de errores en sistemas de ecuaciones lineales Los desarrollos matemáticos en este caso general son un poco más complicados que en los casos precedentes. Se suponen ahora unas perturbaciones en los datos δ A y δ b que producen un error δ x:
( A + δ A ) ( x + δ x ) = b + δ b
(133)
Se supone también que las perturbaciones en A son pequeñas, de modo que se cumple: δA
<
A A 1 = = ⇒ A −1 A −1 A κ ( A )
δ A
A
≤
1 , κ ( A ) ≡ condición numérica de A κ ( A )
(134)
La matriz ( A −1δ A ) es convergente, pues por la hipótesis (134) se cumple que: A −1δ A ≤ A −1 δ A < 1
⇒
( I + A −1δ A ) no singular
(135)
Aplicando a la matriz ( − A −1δ A ) la desigualdad b) de la proposición (125):
( I + A −1δ A )
−1
≤
1 1 ≤ −1 −1 1 − A δ A 1 − A δ A
(136)
Desarrollando la expresión (133) y simplificando: A −1 ( A + δ A )( x + δ x ) = A −1 ( b + δ b ) ⇒
Despejando δ x de la ecuación (137) y tomando normas:
( I + A−1δ A ) ( x + δ x ) = x + A −1δb
(137)
Sistemas de ecuaciones lineales
δ x = ( I + A −1δ A )
−1
pág. 33
A −1 (δ b − δ Ax ) ⇒
δx
≤ ( I + A −1δ A )
−1
A −1 ( δ b + δ A x )
(138)
Dividiendo por x y teniendo en cuenta la desigualdad (136): δx
x
≤
1 1 − A −1
⎛ δb ⎞ A ⎜ + δA ⎟ δA ⎝ x ⎠ −1
b≤A x
⇒
A −1 ⎛ δb ⎞ A + δ A ⎟ ≤ ⎜ −1 x 1 − A δ A ⎝ b ⎠
δx
(139)
En esta expresión se hace aparecer ahora la condición numérica κ ( ), y se llega a: A −1 A −1 A ⎛ δb ⎞ ⎛ δb δA ⎞ A + δ A ⎟ = ≤ + ⎜ ⎜ ⎟ −1 x A ⎠ 1 − A δ A ⎝ b ⎠ 1 − A −1 A δ A ⎝ b A
δx
(140)
Finalmente: δx
x
≤
⎛ δ b δ A ⎞ + ⎜ ⎟ δ A ⎝ b A ⎠ 1 − κ ( A ) A κ ( A )
(141)
Al igual que en las expresiones (129) y (131), es de nuevo la condición numérica el factor que controla la amplificación en el resultado (siempre es κ ( )≥1) de los errores relativos en los datos. Hay que señalar que la condición numérica no depende de la magnitud de los elementos de una matriz: Si una matriz se multiplica por un escalar α la condición numérica permanece invariable.
2.9.4 Expresión general del número de condición El número de condición κ ( A ) definido en la expresión (132) exige que la matriz A sea invertible. De una forma más general, válida para cualquier tipo de matriz, el número de condición se define en la forma: κ ( A ) =
max Ax x =1
min Ax
(142)
x =1
El resultado (142) es el cociente entre la magnitud de la imagen del vector unitario que más crece con la matriz A, dividido por la magnitud de la imagen del vector unitario que menos crece. Nótese que esta definición se aplica también a una matriz rectangular. Si el rango de A es menor que n, ∃y ∈ R n , no nulo, tal que Ay = 0 y, por tanto, min Ax = 0 . Si A x =1
no es la matriz nula, max Ax > 0 . En esta situación, como extensión a la definición de número de x =1
condición, se dice que κ ( A ) = ∞ .
2.9.5 Número de condición para la norma euclídea El número de condición se ha definido en (132) para una norma cualquiera, aunque suponiendo una matriz invertible. Ahora dicha expresión se va a particularizar para la norma euclídea o espectral. La norma euclídea de una matriz A∈Cn×n es la raíz cuadrada positiva del radio espectral de A H A: A 2 = ρ ( A H A )
(143)
Álgebra lineal numérica con Matlab
pág. 34
La matriz A H A es hermítica y todos sus valores propios son reales y no negativos, pues dicha matriz es –al menos– semidefinida positiva . Además, si la matriz A es normal y sus valores propios son λ j, los valores propios de la matriz A H A son:
U H AU = D, U H A H U = D U H ( AA H ) U = DD = diag ( λ1 ,..., λn 2
2
⇒ DD = U H AUU H A H U = U H AA H U
(144)
)
(145)
2
A 2 = ρ ( A H A ) = λn = λn = ρ ( A )
⇒
Si además la matriz A es regular, los valores propios de A –1 son los inversos de los valores propios de A. Por tanto, se verificará que: 1
A −1 2 = ρ ( A − H A −1 ) =
λ12
=
1 λ 1
= ρ ( A −1 )
(146)
De acuerdo con esto y con (132), el número de condición de una matriz normal y regular será: κ
(2)
( A ) = A −1 2 A 2 = ρ ( A −1 ) ρ ( A ) =
λ n λ 1
(147)
Las matrices ortogonales y unitarias son las únicas que tienen todos sus valores propios de módulo unidad, y por tanto son las únicas que tienen número de condición igual a uno. Estas matrices son óptimas para resolver sistemas de ecuaciones sin errores significativos.
2.9.6 Inversas de matrices perturbadas Ahora se va a estudiar la variación en la inversa de una matriz A con una perturbación δ A. Sea una matriz regular A en la que se introduce una pequeña perturbación δ A, tal que A−1δ A es convergente. Se desea acotar la norma del error: −1 ( A + δ A ) − A −1
(148)
Se parte de la siguiente identidad matricial (se puede comprobar multiplicando porB): B −1 = A −1 − B−1 ( B − A ) A −1
(149)
Sutituyendo B por A+δ A en la expresión (149): −1 −1 ( A + δ A ) = A −1 − ( A + δ A ) (δ A ) A −1
(150)
Pasando A−1 al primer miembro y tomando normas: −1
−1
( A + δ A ) − A −1 ≤ ( A + δ A )
δ A
A −1
(151)
El primer factor del segundo miembro se puede transformar en un producto de inversas, y luego se le puede aplicar la desigualdad (136):
( A + δ A)
−1
= (A (I + A
−1
δA)
)
−1
= (I + A
−1
δ A)
−1
A
−1
≤ A
−1
(I + A
−1
δ A)
A −1 ≤ −1 1 − A δ A
−1
(152)
Sustituyendo este resultado en la expresión (151) de la norma del error se obtiene: −1
( A + δ A ) − A −1 ≤
δ A A −1 1− A δA
κ ( A )
−1
−1 2
⇒
( A + δ A ) − A −1 A
−1
≤
δ A
A 1 − A δ A −1
(153)
Se ve que también en este caso el error relativo en el resultado depende del error en los datos a través del número de condición.
Sistemas de ecuaciones lineales
pág. 35
2.9.7 Conclusiones sobre errores en sistemas de ecuaciones lineales La condición numérica o número de condición κ (A) indica cómo se amplifican los errores δ A y δ b. Como mínimo, los errores relativos correspondientes son del orden de la precisión de los datos (eps=10e–16, en un PC trabajando con doble precisión). Además, en el transcurso de los cálculos se introducen errores de redondeo en las operaciones aritméticas. Existen dos formas de estudiar los efectos de estos errores: El análisis directo considera el error en cada operación y trata de estudiar cómo se propaga y acumula hasta el resultado final. Conduce a resultados muy pesimistas. − El análisis inverso (Wilkinson, 1965) considera los errores de redondeo como efectos equivalentes a perturbaciones en los datos iniciales, y trata de acotar dichas perturbaciones. Los resultados son mucho más precisos. Seguidamente se incluye un resultado del análisis inverso de errores suponiendo la factorización : exacta de una matriz perturbada A + δ A = LU −
aij( k ) (δ A )ij ≤ 2n ⋅ ε ⋅ 1≤max i , j ,k ≤ n
(154)
donde n es el número de ecuaciones, ε es la precisión de la máquina (10e–16 para los PCs, trabajando en doble precisión) y aij( k ) son los elementos que van apareciendo en A a lo largo del proceso de factorización. El pivotamiento por columnas basta para evitar que estos elementos crezcan incontroladamente. Téngase en cuenta que las matrices definidas positivas no necesitan pivotamiento y que por tanto, en general, no suelen dar problemas de errores de redondeo salvo que el número de condición sea muy elevado. Si existen dudas al respecto, siempre es conveniente estimar los errores a través de un número de condición aproximado. Matlab dispone de las funciones cond(A,p) y condest(A). La primera de ellas calcula de modo exacto la condición numérica de A utilizando la norma-p (si se omite p, se supone p=2). Este cálculo puede ser muy laborioso si A es de gran tamaño. La función condest(A) calcula de modo económico un valor aproximado de la condición numérica.
n=r m
=
Figura 23. Sistema de ecuaciones lineales incompatible.
2.10 Sistemas de ecuaciones lineales redundantes ( m> n= r) 2.10.1 Aplicación del teorema de la proyección ortogonal Supóngase que se tiene un sistema de m ecuaciones lineales con n incógnitas (m>n), siendo r =n el rango de la matriz (ver Figura 23): Ax = b (155)
Álgebra lineal numérica con Matlab
pág. 36
Como el subespacio Im( A) (columnas de A) es un subespacio de R m de dimensión r
x∈R
(156)
donde r es el residuo o error con que se satisfacen las ecuaciones. El teorema de la proyección ortogonal (ver apartado 1.7) establece que r debe pertenecer al com plemento ortogonal de Im( A), es decir, que es ortogonal a Im( A). Esta condición se puede imponer haciendo que r sea ortogonal a una base de Im( A), por ejemplo, a las columnas de A:
AT r = 0 ⇒ A T ( b − Ax 0 ) = 0
⇒ A T Ax 0 = A T b
(157)
−1
x 0 = ( A T A ) A T b
(158)
El sistema de ecuaciones obtenido en (157) se conoce como sistema de ecuaciones normales , y es típico del método de los mínimos cuadrados. Es interesante estudiar las propiedades de la matriz AT A que aparece en la solución (158): T T T T T T T − La matriz A A es simétrica: (A A) =A (A ) = A A T − La matriz A A, de tamaño n×n, es invertible pues tiene rango r=n, por ser las n columnas de A linealmente independientes. −
La matriz AT A es definida-positiva ya que: T
2
x T ( A T A ) x = ( Ax ) ( Ax ) = Ax > 0, ∀x ≠ 0 −
(159)
La matriz AT A tiene el mismo núcleo o subespacio nulo que A (también aunque r
En resumen, la solución de mínimo error cuadrático del sistema Ax=b, donde el rango de A coincide con el nº de columnas m, se puede escribir por medio de las ecuaciones normales en la forma: −1
x 0 = ( A T A ) A T b
(160)
Gráficamente, las ecuaciones normales tienen la forma mostrada en la Figura 24:
A
=
AT
A
=
AT
AT A
=
Figura 24. Solución de mínimo error cuadrático de u n sistema incompatible.
2.10.2 Sistema de ecuaciones ampliado Las ecuaciones del método de los mínimos cuadrados se pueden también escribir mediante un sis tema de ecuaciones ampliado. El residuo se ha definido en la forma: r = b − Ax (161)
Sistemas de ecuaciones lineales
pág. 37
Además, el residuo es ortogonal a Im( A) lo que se puede expresar matemáticamente en la forma:
A T r = 0
(162)
Expresando conjuntamente las ecuaciones (161) y (162) se llega a la siguiente ecuación matricial:
⎡ I A ⎤ ⎧ r ⎫ ⎧b ⎫ ⎢ AT 0 ⎥ ⎨x ⎬ = ⎨ 0 ⎬ ⎣ ⎦⎩ 0⎭ ⎩ ⎭
(163)
que es un sistema de m+n ecuaciones con m+n incógnitas. Este sistema de ecuaciones ampliado tiene o puede tener algunas ventajas: −
−
La ventaja más importante es que la condición numérica de la matriz ampliada es menor que la de AT A, lo que permite alcanzar la solución con menor error. Si la matriz A es dispersa o sparse (la mayor parte de sus elementos cero) también los tiem pos de cálculo pueden ser inferiores con la matriz ampliada, al menos para ciertos tipos de matrices.
Q
x
=
b
Figura 25. Ecuaciones incompatibles con matriz de columnas ortogonales.
2.10.3 Problema de mínimos cuadrados con matrices de columnas ortogonales La Figura 25 muestra un sistema de ecuaciones lineales incompatible en el que la matriz del sistema es una matriz de columnas ortogonales Q. Como m>n, Q no es una matriz ortogonal, pues no es cuadrada. Considérense pues los sistemas de ecuaciones m>n incompatibles: Qx = b
b ∉ Im ( Q )
(164)
Planteando las ecuaciones normales, la solución de mínimo error cuadrático es:
QT Qx 0 = QT b
⇒
x 0 = QT b
(165)
pues QT es una inversa por la izquierda de Q (no es una verdadera inversa porque no es cuadrada). La proyección de b sobre el Im( Q) está dada por una matriz P definida en la ecuación (46): −1
P = A ( AT A ) AT ,
para A = Q
⇒
−1
P = Q ( QT Q ) QT = QQT
(166)
Así como QT Q = I porque QT es una inversa por la izquierda de Q, la ecuación (166) proporciona un significado para la matriz QQT , que resulta ser una matriz de proyección ortogonal sobre Im(Q). La Figura 26 y la Figura 27 ilustran los tamaños de ambas matrices y los significados citados:
=
= QT Q = I
Figura 26. Producto QT Q=I.
QQT = P Figura 27. Producto QQT =P.
Álgebra lineal numérica con Matlab
pág. 38
2.11 Sistemas de ecuaciones indeterminados ( r= m< n) Supóngase el siguiente sistema de ecuaciones lineales indeterminado, de rango pleno (r =m
(167)
en el que habrá infinitas soluciones, pues hay n – r variables libres.
=
m=r n
Figura 28. Sistema de ecuaciones indeterminado de rango pleno.
Se trata de hallar la solución de mínima norma euclídea del sistema (167), representado en la Figura 28. Recuérdese que la solución general de un sistema indeterminado se puede construir sumando la solución general del sistema homogéneo a una solución particular del sistema completo. Sea N una matriz n×(n – r) cuyas columnas son una base de Ker( A). Si x P es una solución particular de Ax=b, la solución general de dicho sistema viene dada por:
x = x p + Nu
u ∈ R n − r
(168)
Para hallar la solución de mínima norma se puede tener en cuenta que la ecuación (168) representa todas las soluciones del sistema (167), que dependen de los elementos del vector u, que pueden ser considerados como parámetros independientes. Cualquier solución es la suma de una componente en el núcleo de A y otra componente en su complemento ortogonal, que es Im( AT ). Esta suma es directa, lo que implica que los sumandos están unívocamente determinados. La solución de mínima norma será pues el vector que no tenga componente en Ker( A) y que pertenezca a Im( AT ). Sea v el vector de coeficientes de la solución de mínima norma x∗ expresada como combinación lineal de las filas de A:
x∗ = AT v,
v ∈ R r
(169)
Sustituyendo este resultado en el sistema de ecuaciones a resolver:
Ax∗ = AAT v = b ( AAT invertible )
(170)
Despejando el vector v de la expresión (170) y sustituyendo en la expresión (169) de x∗ : −1
v = ( AAT ) b
⇒
−1
x∗ = A T ( AAT ) b
(171)
que es la expresión de la solución de mínima norma. Para que la matriz AAT sea invertible es necesario que sea de máximo rango (r =m), es decir, que las filas de A sean linealmente independientes.
2.12 Problema general de mínimo error cuadrático y mínima norma En los apartados anteriores se han resuelto los casos de solución de mínimo error cuadrático de un sistema de ecuaciones incompatible (m>n) y de solución de mínima norma de un sistema de ecuaciones indeterminado (mn) o las filas (en el caso m
Sistemas de ecuaciones lineales
pág. 39
Considérese el sistema de ecuaciones Ax= b , en el que la matriz A no es de rango máximo (r
Figura 29. Sistema incompatible de dos ecuaciones con dos incógnitas, sin solución de mínimo error única.
m m
A
r
=
n
C
B
r
Figura 30. Descomposición de A en el producto de dos matrices de rango máximo.
Sea cual sea la matriz A (de rango r ), siempre se puede descomponer en el producto de dos matrices de rango máximo de tamaños (m×r ) y (r ×n), tal como se muestra en la Figura 30:
A =BC
(172)
Las columnas de B son una base de Im( A) y las filas de C lo son de Im( AT ). Cada columna de C contiene los coeficientes con los que hay que combinar las columnas de B para obtener cada columna de A; análogamente, cada fila de B contiene los coeficientes con los que hay que combinar las filas de C para obtener cada fila de A. Utilizando la factorización (172), el sistema Ax=b se puede escribir en la forma:
Ax=b ⇒ BCx = b
(173)
Introduciendo ahora un nuevo vector y definido en la forma:
y ≡ Cx
(174)
el sistema de ecuaciones (173) se puede escribir como:
By =b
(175)
Hay que tener en cuenta que la solución de mínimo error cuadrático de Ax=b se corresponde con la de By=b, pues: minn ( b − Ax ) = minn ( b − BCx ) = minr ( b − By ) x∈R
x∈R
y∈R
(176)
Como el rango de la matriz B es máximo (igual al número de columnas), se puede aplicar la expresión (171) de la solución de mínimo error cuadrático, reemplazando A por B: −1
y = ( BT B ) BT b
(177)
Por otra parte, una vez conocido y, la ecuación (174) constituye un sistema indeterminado con una matriz C de rango máximo. La solución de mínima norma viene dada por la expresión (171):
Álgebra lineal numérica con Matlab
pág. 40 −1
x + = CT ( CCT ) y
(178)
Combinando los resultados de mínimo error (177) y mínima norma (178) se obtiene: −1
−1
x + = CT ( CCT ) y = CT ( CCT )
−1
( BT B ) B T b = A + b
(179)
seudoinversa, que se estudiará en detalle en el siguiente apartado. donde la matriz A + es la matriz seudoinversa
2.13 Matriz seudoinversa A+ La solución óptima (179) obtenida para el caso general del sistema Ax=b se puede escribir como:
x + = A + b,
A + ≡ CT ( CCT )
−1
( BT B )
−1
BT
(180)
La matriz A+ se conoce con el nombre de matriz seudoinversa, por analogía con la solución de un sistema compatible y determinado x=A –1b. Esta matriz juega un papel muy importante en la teoría de sistemas de ecuaciones lineales. Características de la solución dada por la ecuación (180): − Si la matriz A es cuadrada e invertible , también lo serán los factores B y C. En este caso la seudoinversa coincide con la inversa:
x + = A +b = CT ( CCT ) −
−1
( BT B )
−1
BT b = CT C− T C−1B −1B − T B T b = C −1B −1b = A −1b
(181)
Si el sistema es incompatible y la matriz A es de rango máximo, la solución viene dada por las ecuaciones normales. En este caso C es cuadrada e invertible, por lo que se verifica: −1
AT A = C T BT BC ⇒ BT B = C − T A T AC−1 ⇒
(BTB ) = C ( A T A )
−1
CT
(182)
Sustituyendo en la expresión (181) de la solución aproximada x+: −1
−1
−1
x + = CT C−T C −1 ( BT B ) BT b = CT C− T C −1C ( A T A ) C T B T b = ( A T A ) A T b
(183)
que es la solución dada por las ecuaciones normales para mínimos cuadrados. Análogamente se puede demostrar que si el sistema es indeterminado y la matriz A de rango máximo, la solución (180) se reduce a la expresión (171) hallada anteriormente para la solución de mínima norma. La matriz pseudoinversa A+ definida en la expresión (180) se puede expresar también en la forma: −
A + ≡ CT ( CCT )
−1
−1
−1
−1
( BT B ) BT = CT ( BT BCCT ) B T = CT ( B T AC T ) B T
(184)
La expresión x+=A+b representa pues, en cada caso, la solución más conveniente para el sistema general de ecuaciones lineales Ax=b. La matriz inversa A –1 sólo existe cuando A es cuadrada y no semesingular. Sin embargo, la matriz seudoinversa A+ existe siempre y tiene algunas propiedades seme –1 jantes a las de la matriz matriz inversa A . Algunas de estas propiedades son (sin demostración): + − La matriz A es única, a pesar de que la factorización A=BC no está unívocamente definida. + + − La seudoinversa de la seudoinversa es la matriz original: ( A ) =A T + T + + T − Los operadores ( ) y ( ) se pueden permutar: (A ) =(A ) −
Relaciones entre rangos: r (A)=r )=r (A+)=r )=r (AA+)=r )=r (A+A)
Sistemas de ecuaciones lineales
pág. 41
−
Si las columnas de A son independientes, la seudoinversa es una inversa por la izquierda :
−
r (A)=n )=n ⇒ A+A=In Si las filas de A son independientes, la seudoinversa es una inversa por la derecha : r (A)=m )=m ⇒ A A+ =Im
(185) (186)
Sin embargo, existen también diferencias importantes entre la inversa y la seudoinversa. En general: +
( AB ) ≠ B + A + AA + ≠ I, A +A ≠ I, AA + ≠ A + A +
(187)
p
( A p ) ≠ ( A + )
seudoinversa no se hace por medio de la factorización A=BC En la práctica, el cálculo de la matriz seudoinversa singulares (DVS), que se verá en el apartado 4.2.4. Sin sino mediante la descomposición de valores singulares embargo, en ciertas ocasiones (por ejemplo cuando la matriz A es muy grande y dispersa) el método expuesto en este apartado para calcular la seudoinversa A+ puede ser preferible a la DVS.
Propiedades adicionales relacionadas con la matriz seudoinversa: −
La matriz P=AA+ es la matriz de proyección ortogonal sobre sobre Im(A Im(A). En efecto, se cumple que Ax+= Pb, Pb, sien+ + + do x la solución de mínimo error en norma cuadrática. Puesto que también se cumple que x =A b: Ax + = AA + b = Pb
−
P = AA+
⇒
(188)
La matriz Q=A+A es la matriz de proyección ortogonal sobre sobre Im(A Im(AT ). Utilizando AT en lugar de A en la pro piedad anterior se tendrá: Q = AT ( A T ) es la matr matriz iz de proy proyec ecci ción ón sobr sobree Im( A T ); Q = QT = A + A +
−
(189)
Las matrices A y P tienen el mismo espacio de columnas, Im(A Im(A)=Im(P )=Im(P): x ∈ Im( P) ⇒ ∃u; Pu = x ⇒ AA + u = x ⇒ A ( A+ u) = x ⇒ x ∈ Im( A )
(190)
Análogamente se demuestra que: x∈Im(A Im(A) ⇒ x∈Im(P Im(P). −
La solución general del del sistema Ax= Ax=b que minimiza el cuadrado de la norma del residuo (sin minimizar la norma de la solución) viene dada por:
( w ≡ vector arbitrario )
x = A +b + ( I − A + A ) w
(191)
pues (I (I – A+A)=(I )=(I – Q) es la matriz de proyección sobre el complemento ortogonal de Im(A Im(AT ), que es Ker(A Ker(A). −
Las cuatro propiedades siguientes caracterizan a la seudoinversa de A (se puede demostrar que A+ es la única matriz que las cumple): 1. AA+A=A AA + A = BC ⋅ CT ( CCT )
−1
( BT B )
−1
BT ⋅ BC = BC = A
(192)
2. A+AA+=A+ A + AA + = C T ( CC T )
−1
( BT B )
−1
−1
BT ⋅ BC ⋅ CT ( CCT )
−1
( BT B )
BT = A +
(193)
3. (AA+)T =AA+ (simetría del producto) T T T T T ( AA + ) = ( A + ) A T = B ( B T B ) ( CC T ) C ⋅ C T B T = B ( B T B ) B T
−
−1
AA + = BC ⋅ CT ( CCT )
−
( BT B )
−1
−
BT = B ( BT B ) B T −1
(194) (195)
Álgebra lineal numérica con Matlab
pág. 42
4. (A+A)T =A+A (simetría del producto) (demostración análoga a la anterior) −
Conviene finalmente recordar que la seudoinversa obtenida a partir de las expresiones (184) es única, aunque la factorización A=BC no A=BC no lo sea.
2.13.1 Interpretación Interpretación del sistema de ecuaciones general Ax=b En el producto Ax=b la matriz A transforma R n en R m. Los vectores x y b se pueden descomponer respectivamente en sus componentes en los subespacios suplementarios Im( AT ) y Ker(A), e Im(A) y Ker(AT ): x = x f + x n
(x
f
∈ Im ( A T ) , x n ∈ Ker ( A ) )
(196)
b = bc + bn
(b
c
∈ Im ( A ) , b n ∈ Ker ( A T ) )
(197)
Sustituyendo estos valores en el sistema Ax=b y teniendo en cuenta que Axn=0 y que bn nunca puede obtenerse como combinación lineal de las columnas de A:
Ax = Ax f + Ax n = Ax f = b c
(198)
La expresión (198) indica que tanto x como x f se se transforman a Im( A). Por otra parte, cada vector b c∈Im(A) proviene de un y sólo un vector x f ∈Im(AT ), como se puede demostrar por reducción al absurdo. Suponiendo dos vectores distintos en Im( AT ) que se transforman en b c:
Ax f = b c ⎫ Ax′ f = b c ⎬⎭
A ( x f − x′f ) = 0
⇒
( x f − x′f ) ∈ Ker ( A )
¡Absurdo!
(199)
Así pues, cada matriz A transforma el subespacio de filas Im( AT ) en el subespacio de columnas Im(A) y esta parte de la transformación siempre es invertible. Al mismo tiempo, hay una parte del vector b, el vector bn, que es inalcanzable por la transformación Ax=b. ¿Qué hace la matriz seudoinversa en un sistema general de estas características? De la definición de la matriz seudoinversa (180) se concluye concluye que: Im ( A + ) = Im ( CT ) = Im ( A T )
⇒
A +b n = 0
(200)
Si A−1 no existe, la seudoinversa A+ la sustituye e invierte "lo que es invertible" en el sistema Ax=b, obteniendo x f a a partir de bc:
A + Ax = A + A ( x f + x n ) = A +Ax f = x f ⎪⎫ ⎬ A + b = A + ( b c + b n ) = A +b c ⎪⎭
x f = A +b c
(201)
La Figura 31, inspirada en el libro de G. Strang 1, representa gráficamente los 4 subespacios fundamentales de la matriz A e interpreta la transformación Ax=b en función de dichos subespacios, resumiendo gráficamente lo establecido en los párrafos precedentes. Obsérvese que la componente de x en Ker(A) no influye en la imagen Ax, y que la componente de b en Ker( AT ) tampoco interviene porque es inalcanzable. inalcanzable.
1
G. Strang, "Introduction to LinearAlgebra", 3rd edition, Wellesley-Cambridge, ch. 4, p. 188.
Sistemas de ecuaciones lineales
pág. 43
m
r
n
x
x
Im(AT ) n
R
Ax
Im(A) r Ax=bc
Ax
x Ax x
0 Ker(A) n
xn
0
b
bn xn
Axn
0
m
Ker(AT )
n−r m−r Figura 31. Relación entre el sistema Ax=b y los cuatro subespacios fundamentales de la matriz A.
R m
Valores y vectores propios
pág. 45
3. Valores y vectores propios 3.1
Definición del problema de valores y vectores propios
Los vectores propios de una matriz A ∈ Cn×n son los vectores x ∈ Cn , x ≠ 0 que se transforman del modo:
Ax = λ x, λ ∈ C
(202)
El escalar λ , que puede ser real o complejo, es el valor propio asociado al vector propio x. La ecuación (202) también se puede expresar como:
( A − λ I ) x = 0
(203)
El vector nulo no se considera vector propio. Para que exista una solución x distinta de la trivial (x=0), el valor propio λ deberá deberá ser raíz de la ecuación característica, que se obtiene igualando a cero el polinomio característico (de grado n): det ( A − λ I ) = 0
(204)
Características generales del problema de valores y vectores propios: − −
Es un problema no lineal en en los valores propios λ y y lineal en en x. Dado que un polinomio de grado n siempre tiene n raíces reales y/o complejas, siempre existen n valores propios , que pueden ser reales o complejos. Si un valor propio es complejo, su complejo conjugado también es valor propio.
−
Sin embargo, no siempre existen n vectores propios . Los valores propios de multiplicidad algebraica m>1 tienen un subespacio propio asociado de dimensión <=m <= m (multiplicidad geométrica). geométrica). Todos los vectores en este subespacio propio son vectores propios. Las matrices que tienen menos de n vectores propios se llaman matrices defectivas.
−
λI ) y no están unívocamente Los vectores propios pertenecen al subespacio nulo de ( Ax – λ determinados: Si x es un vector propio, ( α⋅ x) también lo es. Si un valor propio es múlti ple y su subespacio subespacio propio tiene dimensión dimensión mayor que uno, hay infinitas posibilidades de elegir una base ortonormal de dicho subespacio propio. Ax 1 = λ 1x 1
x1 Ax 2 = λ 2x 2
x2
Ax x
Figura 32. Interpretación geométrica de los valores y vectores propios.
En lo sucesivo, salvo que se indique lo contrario, se considerarán matrices reales. Para las matrices reales, los vectores propios asociados con valores propios conjugados son también vectores comple jos conjugados entre sí.
Álgebra lineal numérica con Matlab
3.2
pág. 46
Interpretación geométrica de los valores y vectores propios
Por lo general, el vector Ax no tiene la misma dirección que x. Los vectores propios son vectores que se transforman sin cambiar de dirección. El valor propio λ determina determina el cambio de longitud.
3.3
Propiedades de los valores y vectores propios
3.3.1 Subespacios propios Los vectores propios asociados con un mismo valor propio λ forman forman un subespacio vectorial de R n. En efecto, si x1 y x2 son vectores propios asociados a un mismo valor propio λ , se verificará:
A (α x1 + β x 2 ) = α Ax1 + β Ax 2 = λ (α x1 + β x 2 )
(205)
lo que demuestra que una combinación lineal de vectores propios asociados con un mismo valor propio es también vector propio. Estos subespacios vectoriales formados por vectores propios se conocen como subespacios propios. Por otra parte, los subespacios propios correspondientes a valores propios distintos sólo tienen en común el vector nulo. En efecto, si x es un vector propio común, por reducción al absurdo:
Ax = λ 1x ⎫ ⎬ 0 = (λ1 − λ 2 )x ⇒ x = 0 Ax = λ 2x ⎭
(206)
Como consecuencia, los vectores propios correspondientes a valores propios distintos son lineal mente independientes independientes .
3.3.2 Relación de los valores propios con la traza y el determinante Igualando coeficientes en las distintas expresiones del determinante det( A−λ I)=0 de la ecuación (204) se deducen las propiedades siguientes (fórmulas de Vieta): det ( A − λ I ) =
a11 − λ a12 a12 a22 − λ
a1n a2 n
an1
an 2
ann − λ
n
= ∑ ( −1)
j
σ n − j λ j
= ( −1) (λ − λ1 ) (λ − λ2 ) (λ − λn ) n
(207)
j = 0
donde σ n− j es la suma de los menores principales de orden ( n − j ) de la matriz A:
= λ1 + λ2 + ... + λn σ 2 = λ1λ2 + λ1λ3 + ... + λn −1λn σ1
σn
(208)
= λ1λ2λ3 ...λn
En concreto, igualando los coeficientes de λ y de λ n se obtiene que: −
La suma de los valores propios es igual a la traza de la matriz.
−
El producto de los valores propios es igual al determinante de la matriz.
3.3.3 Propiedad de "desplazamiento" de los valores propios Si los valores propios de la matriz A son λ i, los valores propios de la matriz ( A−α I), siendo α un ). Los vectores propios xi de la matriz A siguen siendo vectores propios de la maescalar, son ( λ i−α ). triz (A−α I):
Valores y vectores propios
pág. 47
Axi = λi xi ;
restando α x i
⇒
( A − α I ) x i = ( λi − α ) x i
(209)
Esta propiedad es importante en algunos métodos numéricos de cálculo de valores propios.
3.3.4 Casos particulares del problema de valores propios De la expresión del determinante se deduce que los valores propios de una matriz diagonal , o de una matriz triangular superior o inferior son los elementos de la diagonal. Los vectores propios de las potencias de A son los mismos vectores propios de A; los valores pro pios de las potencias de A son las correspondientes potencias de los valores propios de A: Ax = λ x; A2 x = λ Ax = λ 2 x; ... A n x = λ n x
(210)
Análogamente, para las potencias negativas, suponiendo que A sea invertible:
Ax = λ x; A −1Ax = λ A −1x; A −1x = λ −1x
3.4
(211)
Transformaciones de semejanza
3.4.1 Matrices semejantes Dos matrices A y B son semejantes si existe una matriz no singular M tal que: B = M −1AM (212) Se dice que A y B están relacionadas mediante una transformación de semejanza. Las matrices semejantes tienen los mismos valores propios. Partiendo de la definición de valor pro pio de A: Ax = MBM −1x = λ x
(213)
Pre-multiplicando por la inversa de la matriz M:
B(M −1x ) = M −1λ x = λ (M −1x )
(214)
de donde se deduce que λ es un valor propio de B asociado con el vector propio ( M –1x). Es fácil comprobar que dos matrices semejantes A y B tienen el mismo polinomio característico: det(B − λ I) = det(M −1AM − λ M −1M ) = = det(M −1 ( A − λ I)M ) = det M −1 det( A − λ I) det M = det( A − λ I)
(215)
3.4.2 Diagonalización mediante transformaciones de semejanza Si una matriz A ∈ R n×n tiene n vectores propios linealmente independientes se puede diagonalizar mediante una transformación de semejanza. En efecto, expresando la condición de valor y vector propio para los n valores y vectores propios simultáneamente:
Axi = λ i xi
⇒ A [ x1 x 2
x n ] = [x1 x 2
⎡λ 1 ⎢0 xn ] ⎢ ⎢ ⎢ ⎢⎣ 0
0
λ 2
0
0⎤ 0 ⎥⎥ ⎥ ⎥ λ n ⎦⎥
(216)
Llamando ahora P a la matriz cuyas columnas son los n vectores propios independientes y D a la matriz diagonal formada por los valores propios:
Álgebra lineal numérica con Matlab
pág. 48
P ≡ [x1 x 2 AP = PD
⇒
D ≡ diag ( λ1 , λ2 ,..., λ n )
xn ], −1
P AP = D
y
A = PDP
−1
(217)
Para que una matriz sea diagonalizable, debe tener n vectores propios linealmente independientes. Para que sea invertible, debe tener los n valores propios distintos de cero. Las matrices defectivas tienen menos de n vectores propios independientes. No se pueden diagonalizar mediante transformaciones de semejanza. Se pueden reducir a la forma canónica de Jordan.
3.4.3 Reducción a forma triangular mediante transformaciones de semejanza Toda matriz cuadrada A ∈ Cn×n se puede reducir mediante una transformación de seme janza S –1AS a una matriz triangular superior T, cuyos valores propios aparecen en la diagonal. Teorema:
Demostración: Siempre es posible encontrar para A un vector propio x1 asociado con λ 1 (en general, ambos comple jos). Se puede construir una matriz no singular S1, con el vector x1 en la primera columna, tal que al multiplicar S1 por el primer vector de la base natural e1 se tiene: S1e1 = x1
→
S1−1x1 = e1
(218)
Aplicando a la matriz A la transformación de semejanza definida por la matriz S1: ⎡λ bT ⎤ B1 = S1−1 AS1 = S1−1A [ x1 , X1 ] = S1−1 [ λ1 x1 , AX1 ] = ⎡⎣λ 1e1 , S1−1 AX1 ⎤⎦ = ⎢ 1 ⎥ ⎣ 0 A1 ⎦
(219)
Los valores propios de B1 son los de A, pues son matrices semejantes; los de A1 son los de A con la multiplicidad de λ 1 disminuida en una unidad. Ahora a la matriz A1 se le puede aplicar idéntico razonamiento que a A, pues también tendrá al menos un vector propio asociado con λ 2 . La segunda transformación de semejanza se deberá hacer con una matriz que respete los ceros introducidos por la primera, es decir mediante una matriz que tenga la forma: ⎡1 0T ⎤ P2 = ⎢ ⎥ ⎣0 S 2 ⎦
(220)
La matriz B 2 = P2−1S1−1AS1P2 es semejante a A, y tiene ceros debajo de la diagonal en las dos primeras columnas. Continuando con transformaciones de semejanza de este tipo se llega finalmente a transformar A en una matriz triangular superior T, cuyos valores propios están en la diagonal. La forma canónica de Jordan es un caso particular de matriz triangular superior a la que cualquier matriz es reducible por semejanza.
3.4.4 Transformaciones de semejanza unitarias Las matrices unitarias son una generalización de las matrices ortogonales al campo complejo. Estas matrices U∈Cn×n cumplen la propiedad de que su inversa es la conjugada y transpuesta: U −1 = U H ⇔ U H U = UU H = I (221) Se llama transformación de semejanza unitaria a una transformación de semejanza realizada mediante una matriz unitaria: B = U H AU,
siendo U −1 = U H
(222)
Las transformaciones de semejanza unitarias tienen todas las propiedades de las transformaciones de semejanza (y algunas más). Las matrices A y B relacionadas por la ecuación (222) se dicen uni tariamente semejantes. Se llaman matrices unitariamente diagonalizables (triangularizables) a las matrices que se pueden reducir a forma diagonal (triangular ) mediante una transformación de semejanza unitaria . Hay matrices diagonalizables que no son unitariamente diagonalizables.
Valores y vectores propios
pág. 49
Se cumple que si las matrices A y B son unitariamente semejantes , las matrices A H A y B H B tam bién lo son (con la misma matriz U):
B = U H AU, B H = U H A H U;
⇒
U H A H AU = ( U H A H U )(U H AU ) = B H B
(223)
Las transformaciones de semejanza conservan el polinomio característico y por tanto la traza, luego la suma de los cuadrados de los módulos de los elementos de la matriz se conserva en las transformaciones de semejanza unitarias: n
∑
i , j =1
2
aij = traza ( A A ) = traza ( B B ) = H
n
H
2
∑ bij
(224)
i , j =1
3.4.5 Lema de Schur Enunciado: Para cualquier matriz cuadrada
A Cn×n existe una matriz unitaria U tal que T=U H AU
es triangular superior. Demostración: Se demuestra por inducción. La proposición es cierta para n=1. Supóngase que toda matriz (n –1)×(n –1) es unitariamente semejante a una matriz triangular superior. Para cualquier matriz n×n siempre es posible encontrar un vector propio x1 asociado con un valor propio λ 1 (aunque λ 1 sea múltiple, su subespacio propio tendrá al menos un vector propio). Construyendo una matriz unitaria U1 cuya primera columna es x1 y eligiendo las demás columnas de modo que sean ortonormales, se verificará: ⎡λ1 * ⎢0 * AU1 = U1 ⎢ ⎢ ⎢ ⎣0 *
*⎤ *⎥⎥ ⎥
⎥ *⎦
⎡λ 1 * ⎢0 * ⎢ U H AU = 1 1 ⎢ ⎢ ⎣0 *
*⎤ ⎡λ 1 ⎢ * ⎥⎥ ⎢ 0 = ⎥ ⎢ 0 ⎥ ⎢ * ⎦ ⎢⎣ 0
* * *⎤
B
⎥ ⎥ ⎥ ⎥ ⎥⎦
(225)
Pero por hipótesis la matriz B es una matriz (n –1)×(n –1) unitariamente semejante a una matriz triangular superior luego A es unitariamente semejante a una matriz triangular T. La matriz A se podrá finalmente expresar B = U B TB U H B como: ⎡1 ⎢ 0 U 2 = ⎢⎢ 0 ⎢ ⎢⎣0
0 0 0 ⎤
U B
⎥ ⎥ ⎥ ⎥ ⎥⎦
⎡λ1 ⎢ H H H ⎢ 0 U 2 U1 AU1U 2 = U 2 ⎢ 0 ⎢ ⎢⎣ 0
* * * ⎤
B
⎡λ 1 ⎥ ⎢ ⎥U = ⎢0 ⎥ 2 ⎢0 ⎥ ⎢ ⎥⎦ ⎢⎣ 0
* * * ⎤
TB
⎥ ⎥ ⎥ ⎥ ⎥⎦
(226)
Desde un punto de vista constructivo, la matriz A se puede reducir a forma triangular determinando un vector propio para B y construyendo a partir de él la matriz U B, y continuando sucesivamente hasta que la matriz se haga triangular.
3.5
Matrices normales
3.5.1 Definición, casos particulares y propiedades Una matriz A∈Cn×n se llama normal cuando el producto por su matriz conjugada transpuesta A H es conmutativo:
AA H = A H A Casos particulares de matrices normales: n×n − En C todas las matrices hermíticas, antihermíticas y/o unitarias son normales. −
En R n×n todas las matrices simétricas, antisimétricas u ortogonales son normales.
También existen matrices normales que no pertenecen a las categorías precedentes. Todas las matrices B unitariamente semejantes a una matriz normal A son también normales: −
(227)
Álgebra lineal numérica con Matlab
pág. 50
BB H = ( U H AU )( U H A H U ) = U H AA H U = U H A H AU = U H A H ( UU H ) AU = B H B
(228)
Dada una matriz normal A, los vectores propios de A H son los mismos que los de A, y los valores propios de A H son los conjugados de los valores propios de A:
U H AU = D ⇒ U H A = DU H
⇒ A H U = UD H = UD
(229)
3.5.2 Teorema espectral para matrices normales A∈Cn×n una matriz cuyos valores propios son λ1, λ2, ..., λn. Las siguientes afirmaciones son equivalentes entre sí (cada una de ellas implica todas las demás): a) La matriz A es normal. b) La matriz A es unitariamente diagonalizable.
Enunciado: Sea
c) Se verifica la igualdad:
n
∑
i , j =1
2
n
aij = ∑ λ i
2
i =1
d) La matriz A tiene n vectores propios ortonormales entre sí. Demostración. Se demostrará la equivalencia entre a) y b), entre b) y c) y entre b) y d): a)
b): Una matriz es unitariamente diagonalizable si y sólo si es normal .
Por el Lema de Schur, toda matriz A∈Cn×n es unitariamente semejante a una matriz triangular superior T: U H AU = T
(230)
Si A es normal, también lo será T, pues la condición de normal es conservada por las transformaciones de semejanza unitarias. Se va a demostrar que si una matriz triangular es normal , en realidad es diagonal . La Figura 33 muestra la forma del producto T H T, mientras que la Figura 34 lo hace con el producto TT H . Si T es normal, ambos productos deben dar el mismo resultado.
0
0
H
H
T T
TT
0
0
Figura 33. Forma del producto T H T.
Figura 34. Forma del producto TT H .
De la igualdad de elementos (1,1) se concluye que en la primera fila de T sólo el elemento de la diagonal es distinto de cero: ⎫ ⎪ n n ⎬ ⇒ 2 2 ( TT H )11 = ∑ t1 j t1 j = t11 + ∑ t1 j ⎪ j =1 j =2 ⎭
( T H T )11 = t11t11 = t 11
2
n
∑t j = 2
2
1 j
= 0 ⇔ t1 j = 0, j = 2,3,..., n
(231)
De la igualdad de elementos (2,2) se concluye que en la segunda fila de T sólo el elemento de la diagonal es no nulo: ⎫ ⎪ n n ⎬ ⇒ 2 2 ( TT H )22 = ∑ t2 j t2 j = t22 + ∑ t2 j ⎪ j = 2 j =3 ⎭
( T H T )22 = t12t12 + t22t 22 = t 22
2
n
∑t
2 2 j
=0
⇔ t2 j = 0, j = 3, 4,..., n
(232)
j = 3
Prosiguiendo del mismo modo con los restantes elementos de la diagonal de ambos productos se demuestra que T y T H deben ser diagonales. Por tanto, toda matriz A normal es unitariamente diagonalizable. El recíproco, que toda matriz unitariamente diagonalizable es normal, se cumple porque toda matriz diagonal es nor mal , y por tanto también lo será A que es unitariamente semejante a ella.
Valores y vectores propios
b)
pág. 51
c): En una matriz unitariamente diagonalizable se cumple que:
n
∑a
2
ij
i , j =1
n
2 = ∑ λ i . i=1
La traza se conserva en las transformaciones de semejanza, pues es el 2º coeficiente del polinomio característico, que es un invariante. Si A y D son unitariamente semejantes, aplicando la propiedad vista en la expresión (224): n
∑a
ij
i , j =1
2
n
= traza ( A H A ) = traza ( DH D ) = ∑ λ i i =1
2
(233)
esta propiedad indica que la matriz triangular T = U H AU , unitariamente semejante a A y con los valores propios en la diagonal, debe ser diagonal , pues si se cumplen simultáneamente la condición c) y la ecuación (224), la suma de los cuadrados de los módulos de los elementos de fuera de la diagonal debe ser nula.
Recíprocamente,
b)
d): Toda matriz es unitariamente diagonalizable si y sólo si tiene n vectores propios ortonormales entre sí.
Si A es diagonalizable mediante una matriz unitaria U: U H AU = D
⇔ AU = UD
(234)
Las n columnas de U son vectores propios ortonormales de A y los elementos de D los valores propios. Recíprocamente, si A tiene n vectores propios ortonormales que son las columnas de U, la misma expresión (234) indi-
ca que es unitariamente diagonalizable.
3.5.3 Corolarios del teorema espectral Si A es una matriz normal (y por tanto unitariamente diagonalizable) se cumple que: − La matriz A es hermítica si y sólo si sus valores propios son reales: A = A H ⇔ U H AU = D,
−
U H A H U = D H ⇔ D = D H ⇔ D real
(235) Como caso particular, las matrices reales y simétricas siempre tienen valores propios reales. La matriz A es antihermítica si y sólo si sus valores propios son imaginarios puros:
A = − A H ⇔ U H AU = D, U H A H U = D H = −D ⇔ D : d ii imaginarios puros −
La matriz A es unitaria si y sólo si sus valores propios tienen módulo unidad :
A H A = I ⇔ A H A = ( UD H U H )( UDU H ) = UD H DU H = I ⇔ D H D = I −
(236)
(237)
La matriz A es normal si y sólo si los vectores ( Ax) y (A H x) tienen la misma norma euclí dea (sin demostración):
3.5.4 Descomposición espectral de matrices normales Toda matriz normal –por tanto unitariamente diagonalizable– admite una descomposición espectral mediante la que se expresa como suma de matrices de rango 1, en la forma:
⎤ ⎡u H 1 ⎤ ⎢ ⎥ n ⎥ H H A = UDU = [u1 ⎥ = ∑ λ k u k u k (238) ⎢ ⎥ H k =1 λ n ⎥⎦ ⎢⎣u n ⎥⎦ Esta descomposición es una forma alternativa de considerar la diagonalización unitaria de la matriz A. Se estudiará con algo más de detalle en el apartado 4.3. ⎡λ 1 u n ] ⎢⎢ ⎢⎣
Álgebra lineal numérica con Matlab
3.6
pág. 52
Formas cuadráticas y transformaciones de congruencia
3.6.1 Definición y propiedades de las formas cuadráticas Una forma cuadrática en R n es un polinomio homogéneo de grado dos en n variables, es decir función ϕ : R n R de la forma: ϕ ( x ) = ϕ ( x1 , x2 ,..., xn ) =
n
∑a xx ij i
j
= xT Ax, A ∈ R n×n , A T = A
(239)
i , j =1
Las formas cuadráticas pueden reducirse a forma diagonal mediante cambios de coordenadas en el vector x. Como A es simétrica, una forma particular de reducir la forma cuadrática a forma diagonal es mediante una transformación de semejanza unitaria que la diagonalice. Sin embargo, como no se exige que el cambio de coordenadas sea unitario u ortogonal, una misma forma cuadrática puede transformarse en muchas formas cuadráticas diagonales diferentes, pero sólo en una forma diago nal canónica . La forma diagonal canónica es aquella en la que todos los elementos de la diagonal son "1", "–1" ó "0" (en este orden).
3.6.2 Transformaciones de congruencia Una transformación de congruencia se define en la forma: B = C T AC
( C no singular )
(240)
Estas transformaciones surgen de realizar un cambio de base en una forma cuadrática: ϕ
= xT Ax; x = Cy;
ϕ = y T CT ACy
= y T By ; B = CT AC
(241)
Las transformaciones de semejanza unitaria son también transformaciones de congruencia, pero no cualquier transformación de congruencia es de semejanza o de semejanza unitaria. Permutar unas filas y las correspondientes columnas, multiplicar una fila y la correspondiente columna por un escalar, combinar filas y combinar del mismo modo las correspondientes columnas, son ejemplos de transformaciones de congruencia en una forma cuadrática que pueden llevarla a una forma diagonal o a la forma diagonal canónica.
3.6.3 Ley de Inercia de Sylvester Enunciado: Todas las formas diagonales de una misma forma cuadrática tienen el mismo número de coeficientes positivos, negativos y nulos. Como consecuencia, las transformaciones de congruencia conservan el signo de los valores propios, es decir, el número de valores propios negativos, iguales a cero y positivos. Demostración: Supóngase dos formas diagonales de una misma forma cuadrática (una misma forma cuadrática en dos bases distintas) en un espacio E de dimensión n. Por la igualdad de rangos ambas tienen el mismo número de elementos no nulos. ϕ ( u ) =
y12 + ... + y p2 − y 2p+1 − ... − yq2−1 − yq2 − y q2+1 −... − y r 2 = 2 1
2 p
= z + ... + z + z
2 p +1
+ ... + z
2 q−1
2 q
+z −z
2 q+1
−... − z
(242)
2 r
Hay que demostrar que p=q. Se hará por reducción al absurdo. Supóngase, por ejemplo, que p
zq +1 = z q+ 2 = ... = z n = 0
(243)
Tal vector u existe y es no trivial, pues hay (n – q+p)
Valores y vectores propios
pág. 53 ϕ ( u ) = − y p2 +1 − ... − yr2
< 0,
ϕ ( u ) = z12
+ ... + z q2 > 0
(244)
lo cual es absurdo. De forma análoga se puede demostrar que p>q es imposible, luego p=q.
Una consecuencia práctica importante de la ley de inercia de Sylvester es que para cualquier matriz simétrica con una factorización A=LDLT , los signos de los pivots que aparecen en D coinciden con los signos de los valores propios de A, pues A y D están relacionadas por una transformación de congruencia.
3.6.4 Matrices definidas-positivas Una matriz simétrica A es definida-positiva si cumple: x H Ax > 0, ∀x ≠ 0
(245)
Por otra parte, se dice que A es definida negativa cuando (– A) es definida positiva. La condición (245) es equivalente a cualquiera de las cuatro condiciones siguientes: −
Todos los valores propios de A son estrictamente positivos: λ i > 0, ∀i
(246)
−
Todos los menores principales de A son positivos: det( A k ) > 0, ∀k
(247)
−
Todos los pivots en la descomposición A=LDLT son positivos: dii > 0, ∀i
(248)
−
singular ) La matriz A se puede descomponer en la forma: A = R H R ( R no
(249)
Demostración: 1º)
(245) ⇔ (246). Se pretende demostrar la siguiente equivalencia de propiedades: xT Ax > 0, ∀x ≠ 0
⇔ Ax j = λ j x j , λ j > 0 ∀ j
(250)
En efecto, si xT Ax > 0, ∀x ≠ 0 , se cumplirá: Ax j = λ j x j
⇒
xTj Ax j = λ j x T j x j
⇒
λ j
>0
(251)
Recíprocamente, si todos los valores propios son positivos y los vectores propios ortonormales, la matriz
A es
definida positiva: AQ = QD ⇒ QT AQ = D, d ii = λ i > 0 n
x = Qy, xT Ax = yT QT AQy = y T Dy = ∑ λ i yi2 > 0, ∀x = Qy ≠ 0
(252)
i =1
2º)
(245) ⇔ (247). Ahora se va a demostrar lo siguiente: xT Ax > 0, ∀x ≠ 0 ⇔
det ( A k ) > 0, 1 ≤ k ≤ n
(253)
Si xT Ax > 0, ∀x ≠ 0 , se verifica que det A = λ1λ2 ...λ n > 0 . Hay que demostrar que también se cumple para k 0 ⎣ ∗ ∗⎦ ⎣ 0 ⎦
(254)
Como consecuencia todos los valores propios de Ak serán positivos y también serán positivos sus determinantes det Ak , que son el producto de los k valores propios correspondientes. Recíprocamente, si det Ak >0 "k , la matriz A es definida positiva. Se demostrará por inducción en n. Está claro
que si det A1>0, A1=a11>0
Se supone que An –1 es definida positiva. Entonces existe una transformación de congruencia tal que ST A n – 1S=In –1. La matriz A se puede escribir como:
Álgebra lineal numérica con Matlab
pág. 54
u⎤ ⎡A A = ⎢ nT−1 ⎥ , u ∈ R n−1 , a ∈ R; a⎦ ⎣u
⎡ST 0 ⎤ ⎡ A n−1 u ⎤ ⎡ S 0 ⎤ ⎡I n−1 v ⎤ ⎢ T ⎥ ⎢ uT a ⎥ ⎢0T 1 ⎥ = ⎢ vT a ⎥ ⎦⎣ ⎦ ⎣ ⎦ ⎣0 1⎦ ⎣
(255)
que es una matriz congruente con A y "casi" diagonal. Las últimas fila y columna de esta matriz se pueden anular mediante transformaciones de congruencia (añadiéndoles alguna de las n –1 primeras filas y columnas, multiplicadas por el factor adecuado). Estas transformaciones no alteran las n –1 primeras filas y columnas, por lo que se llega finalmente a una matriz diagonal: 0⎤ ⎡I D = ⎢ nT −1 ⎥ , ⎣ 0 d ⎦
det D = d
(256)
Este determinante (de una matriz congruente con A) debe tener el mismo signo que det A, que es positivo por hipótesis. Así pues la matriz D es definida positiva, y también lo será la matriz A que es congruente con ella. Esta propiedad que caracteriza a las matrices definidas positivas se conoce como criterio de Sylvester. 3º)
(247) ⇔ (248). Ahora se va a demostrar lo siguiente: det A k > 0, 1 ≤ k ≤ n ⇔ pivots eliminación Gauss a (jj j −1) > 0
(257)
Si det Ak > 0 , los pivots a jj( j −1) > 0 . La eliminación de Gauss sin pivotamiento conserva los determinantes de las matrices Ak . Después de triangularizar la matriz A, el pívot k se puede escribir en la forma: akk( k −1) = det A k det A k −1 > 0
(258)
y es positivo si los dos determinantes de Ak y Ak –1 son positivos. Recíprocamente,
si todos los pivots son positivos también los determinantes deberán serlo, pues tanto el primer pívot como el primer determinante son a11.
4º)
(245) ⇔ (249). Se demostrará finalmente la equivalencia siguiente: xT Ax > 0 ∀x ≠ 0 ⇔ A = R T R , rango R = n
(259)
Si A=R T R la matriz A es definida positiva: A = RT R
⇒ xT Ax = xT R T Rx = Rx > 0 ∀x ≠ 0, pues Ker ( R ) =0 2
(260)
si A es simétrica y definida positiva siempre admite la factorización de Choleski, como se puede demostrar a partir de la factorización: Recíprocamente,
1
1
A = LDLT =LD 2 D 2LT = R T R
(261)
pues todos los pivots (elementos de la diagonal de D) son positivos. Es evidente que la matriz R es no singular.
3.6.5 Matrices semi-definidas positivas e indefinidas. Se dice que una matriz simétrica A es semi-definida positiva cuando cumple la condición:
∀x ∈ R n , xT Ax ≥ 0 y ∃y ≠ 0, tal que y T Ay = 0
(262)
Las matrices semi-definidas positivas tienen valores propios positivos y nulos, pero no pueden tener ningún valor propio negativo. Siempre son singulares y tienen determinante nulo. Además, su núcleo Ker(A) tiene al menos dimensión 1, como se deduce directamente de su definición. Las matrices semi-definidas positivas A ∈ R n×n se pueden factorizar en la forma A=R T R , siendo rango(R )
∃y , z ≠ 0 tal que y T Ay > 0 y z T Az < 0
(263)
Valores y vectores propios
pág. 55
Las matrices indefinidas tienen necesariamente valores propios positivos y negativos , y también pueden tener valores propios nulos. En este último caso serán singulares, tendrán determinante nulo y un núcleo no vacío. Estas matrices no admiten factorización real en la forma A=R T R .
3.7
Cociente de Rayleigh
3.7.1 Definición y relación con los valores y vectores propios Para matrices simétricas se define el cociente de Rayleigh en la forma: x T Ax R(x ) = T x x
(264)
Sustituyendo x por el vector propio xi se obtiene el valor propio λ i:
xTi Ax i x iT λ i x i R( xi ) = T = T = λ i xi x i xi xi
(265)
Si los valores propios se ordenan de modo que λ1 < λ2 < ... < λ n , el cociente de Rayleigh es mínimo para x=x1, máximo para x=xn y presenta un valor estacionario para cualquier otro vector propio: 2Ax(xT x) − 2(xT Ax)x ∇ R(x) = =0 (xT x)2
⇒
2Ax(xT x ) − 2(x T Ax )x = 0
(266)
De aquí se concluye que la condición de valor estacionario equivale a que el cociente de Rayleigh cumpla la ecuación de valores propios:
x T Ax Ax − T x = Ax − R(x )x = 0 x x Expresando x en base de vectores propios x=Qy:
(267)
xT Ax y T QT AQy y T Λy λ1 y12 + λ2 y22 + ... + λ n yn2 R (x) = T = = T = = xx yT y y y y12 + y22 + ... + yn2 (268) 2 2 2 2 2 2 2 2 ( λ2 − λ1 ) y2 + ... + ( λn − λ1 ) y n y + y + ... + yn λ1 y1 + λ2 y2 + ... + λ n yn λ = λ1 − λ1 12 22 + = + ≥ λ1 1 y1 + y2 + ... + y2n y12 + y22 + ... + yn2 y12 + y22 + ... + y2n Si los valores propios se suponen ordenados de menor a mayor: λ1 < λ2
< ... < λ n
(269)
todos los términos del numerador de la última fracción son positivos y también la fracción lo será. Por tanto el mínimo de R(x) vale λ 1 y se produce cuando el vector x es x1 ( y1=1, y j=0, j≠1). De forma análoga se puede demostrar que el máximo de R(x) vale λ n y se produce cuando x es el valor propio xn ( yn=1, y j=0, j≠n)
3.7.2 Error en los valores propios estimados mediante el cociente de Rayleigh Si en el cociente de Rayleigh la aproximación en el vector propio es lineal, la aproximación en el valor propio es cuadrática. En efecto, si x=xi+ε v, siendo v el error en xi, se cumple que R(x)=λ i+O(ε 2).
Álgebra lineal numérica con Matlab
pág. 56
Demostración: R ( x i + ε v ) =
( xTi + ε v T ) A(x i + ε v ) x Ti Ax i + 2ε xTi Av + ε 2v T Av = x Ti x i + 2ε x Ti v + ε 2v T v (xTi + ε v T )(x i + ε v )
(270)
Como el error v no tiene componente en xi, se cumplirán las relaciones: v = ∑ α j x j ,
xTi x j = δij ,
Axi = λxi ,
j ≠i
x iT Av = 0,
v T x i = 0,
v T Av = ∑ λ j α 2j
(271)
j ≠i
Sustituyendo en la expresión (264) del cociente de Rayleigh: λi
R ( xi + ε v ) =
+ 0 + ε 2 ∑ λ jα 2j −1 ⎛ j ≠i 2 2 ⎞⎛ 2 2⎞ = ⎜ λi + ε ∑ λ jα j ⎟⎜ 1 + ε ∑ α j ⎟ 1 + ε 2 ∑ α j2 j ≠i j ≠i ⎝ ⎠⎝ ⎠
(272)
j ≠i
(
)
Desarrollando en serie el segundo factor (1 − x )−1 = 1 − x + x 2 − x3 + ... y agrupando términos: ⎛ ⎞ R ( xi + ε v ) = λi + ε 2 ⎜ ∑ λ j α 2j − λi ∑ α 2j ⎟ + ε 4 ( ) + ... j ≠i ⎝ j ≠i ⎠
(273)
3.7.3 Teorema mini-max y maxi-min (Teoremas de Courant-Fischer) Se comenzará explicando el teorema maxi-min para el 2º valor y vector propio. El mínimo valor del cociente de Rayleigh se obtiene para x = x1 y vale λ 1 . Supóngase que ahora se quiere de nuevo calcular el mínimo de R(x) pero con la condición de que x sea ortogonal a un vector dado w. Como hay menos libertad para elegir x que antes, el mínimo que se obtenga será mayor que el mínimo obtenido cuando no había restricciones, que era λ 1 . Supóngase que ahora se toma otro w distinto y se vuelve a repetir el proceso, y así con todos los w posibles. Todos los mínimos obtenidos serán mayores o igual que λ 1 , pero por lo demás cada uno tendrá una magnitud diferente. Para el segundo valor y vector propio, el teorema maxi-min establece que el máximo de los mínimos de R(x), sujeto a la restricción x⊥w, siendo w un vector arbitrario, es λ 2 y se obtiene cuando x=x2. En el caso general, para r ≥2, el teorema maxi-min establece que:
⎛ xT Ax ⎞ λ r = max ⎜ min T ⎟ wi ⎝ x xx ⎠
r = 1, 2,...n
siendo xT w i = 0,
i = 1, 2,...,( r −1)
(274)
Demostración: Expresando x en base de vectores propios y sustituyendo en R(x): ⎛
⎛ α12 λ1 + ... + αr2 λr + αr2+1λr +1 + ... + αn2 λn ⎞ ⎞ ⎟⎟ 2 2 2 2 α1 + ... + α r + αr +1 + ... + αn ⎠ ⎠ ⎝
R = max ⎜ min ⎜ wi
⎝
α j
(275)
Operando de modo que aparezca el valor propio λ r : ⎛
⎛ ⎝
R = max ⎜ min ⎜ λ r − wi
⎝
α j
2
α1 ( λr
− λ1 ) + ... + αr2−1 (λr − λr −1 ) + αr2+1 ( λr − λr +1 ) + ... + αn2 ( λr − λn ) ⎞ ⎞ ⎟⎟ α12 + ... + α r2 + αr2+1 + ... + αn2 ⎠⎠
(276)
Las incógnitas de la minimización son ahora los α j, que están restringidos por la condición de que x sea ortogonal a los vectores wi: w Ti x = w iT ∑ α j x j = 0,
i = 1, 2,..., r − 1
(277)
Los valores de α i deben minimizar el paréntesis interior de la expresión (276), a la vez que satisfacen las restricciones (277). La fracción, que tiene signo negativo (–), deberá ser lo mayor posible. Para ello, α r+ 1=...=α r=0, pues así se anulan los términos negativos (que restan) del numerador de la fracción. Los (r −1) valores de (α 1, α 2, ..., α r− 1) se deben elegir
Valores y vectores propios
pág. 57
de modo que se cumplan las (r −1) restricciones (277). Es evidente que R(x)<λ r pues a λ r se le está restando algo positi vo. Sin embargo, para wi=xi, se cumple que α 1=α 2= ...=α r− 1=0 y se alcanza la igualdad R(x)=λ r. Con esto queda demostrado el teorema maxi-min. De forma análoga se puede demostrar en teorema mini-max. Enunciado conjunto de los teoremas mini-max y maxi-min: El valor propio λ r se puede obtener de
las expresiones siguientes:
⎛ ⎛ ⎛ xT Ax ⎞ ⎞ ⎛ x T Ax ⎞ ⎞ λ r = min ⎜ max ⎜ T ⎟ ⎟ = max ⎜ min ⎜ T ⎟ ⎟ w ⎝ x ⊥w ⎝ x x ⎠ ⎠ w ⎝ x⊥w ⎝ x x ⎠ ⎠ siendo xT w i = 0, r = 2,..., n − 1 i = 1, 2,..., ( r − 1) i
i
i
(278)
i
donde los wi son conjuntos de r vectores arbitrarios de R n.
3.7.4 Interpretación geométrica de los teoremas mini-max y maxi-min La búsqueda de máximos y mínimos del cociente de Rayleigh se puede interpretar como la búsqueda de los puntos de un elipsoide que están a una distancia máxima y mínima de su centro (ver Figura 35). La restricción de ser perpendicular a un vector w obliga a buscar mínimos en la elipse que resulta de cortar el elipsoide con un plano perpendicular a dicho vector w. Se corta el elipsoide con muchos planos diferentes y se calcula el mínimo. El máximo de los mínimos es λ 2 y se obtiene cuando w=x1.
Figura 35. Elipsoide en R 3.
3.7.5 Propiedad de "separación" de los valores propios Si A(m) es la matriz que resulta de suprimir las últimas m filas y columnas de la matriz A, se cumple que los (n−m−1) valores propios de A(m+1) separan –están entre– los ( n−m) valores propios de A(m). Demostración: Esta propiedad se va a demostrar para A y A(1). La matriz A(1) resulta de suprimir la última fila y columna de A. Se trata de demostrar que se cumple la relación: λr
≤ λr(1) ≤ λ r +1
r = 1,2,..., n −1
(279)
Considérense las tres caracterizaciones maxi-min siguientes: ⎛ xT Ax ⎞ ⎟ T ⎝ x ⊥ wi x x ⎠
a) λ r +1 = max ⎜ min wi
xT wi = 0
i = 1,2,..., r
(280)
Álgebra lineal numérica con Matlab
b)
(1) λ r
pág. 58
⎛ xT Ax ⎞ = max ⎜ min ⎟ T wi ⎝ x ⊥ wi x x ⎠
⎧xT wi = 0 i = 1,2,..., r ⎨ w w e arbitrario para i = 1,2,..., r − 1 = r n ⎩ i
⎛ xT Ax ⎞ ⎟ T ⎝ x ⊥ wi x x ⎠
c) λ r = max ⎜ min wi
xT w i = 0
i = 1,2,..., r − 1
(281) (282)
A más libertad para elegir los vectores wi se podrán encontrar mayores máximos de los mínimos. La condición a) deja más libertad para elegir los wi que b), luego λr +1 ≥ λ r (1) . La condición b) deja más libertad para elegir los wi que c), luego λr(1) ≥ λ r .
3.8
Valores y vectores propios generalizados
3.8.1 Introducción a partir del problema estándar El problema generalizado se puede introducir a partir del problema de valores y vectores propios estándar: Supóngase que la matriz A∈R n×n es simétrica. Sea el problema de valores y vectores pro pios estándar, planteado en la forma: Ax = λ x (283) Sea C una matriz relacionada con A mediante la transformación de congruencia:
C = PT AP
(284)
donde la matriz P es una matriz cuadrada no singular. La ley de inercia de Sylvester asegura que los valores propios de C tienen los mismos signos que los de A, aunque en general tendrán diferentes valores numéricos. Sean y y µ los vectores y valores propios de C:
Cy = µ y,
PT APy = µ y
(285)
Haciendo x=Py, la ecuación (285) se puede poner en la forma:
x = Py , y = P −1x
−1
−1
⇒ Ax = µ P −T P −1x = µ ( P T P ) x ⇒ B ≡ ( P T P ) , Ax = µ Bx (286)
que es el llamado problema generalizado de valores y vectores propios . Muchas de las propiedades del problema estándar del que se ha partido se conservan en el problema generalizado. Por ejemplo, si la matriz C es simétrica y semidefinida positiva sus valores propios serán reales y no negativos. Estas mismas características tienen los valores propios generalizados de A y B cuando A es simétrica y semidefinida positiva, y B es simétrica y definida positiva (lo debe ser por la forma en que ha aparecido a partir de P, según (286)).
3.8.2 Planteamiento del problema generalizado de valores y vectores propios Este problema se plantea matricialmente mediante la siguiente ecuación en x y λ : Ax = λ Bx
(287)
donde A y B son ser simétricas y definidas positivas –al menos una de ellas, de ordinario B –. La expresión anterior se puede poner en la forma:
( A − λ B ) x = 0
(288)
Para que exista una solución x distinta de la trivial, el valor propio λ deberá ser raíz del polinomio de grado n que resulta de desarrollar el determinante: det ( A − λ B ) = 0
(289)
Valores y vectores propios
pág. 59
Propiedades del problema de valores y vectores propios generalizado: − −
−
Es un problema no lineal en los valores propios λ y lineal en los vectores x. Siempre existirán n valores propios, que pueden ser reales o complejos según sean las matrices A y B. Los valores propios de multiplicidad m>1 tienen un subespacio propio asociado de dimensión ≤m. Todos los vectores en este subespacio propio son vectores propios.
) y no están unívocamente determinados: Si x Los vectores propios pertenecen a Ker( A – λB es un vector propio, α x también lo es. Este problema se considerará con menos generalidad que el problema estándar: se suele suponer que las matrices son reales y –casi siempre– simétricas y definidas positivas. La interpretación geométrica de este problema generalizado se puede ver en la Figura 36. −
Ax1 = λ 1Bx1
Bx1 x2
x1
Bx 2
Ax 2 = λ 2Bx 2
Bx x
Ax
Figura 36. Interpretación geométrica del problema generalizado de valores y vectores propios.
Se pueden hacer algunos comentarios sobre la Figura 36: − En general, los vectores Ax y Bx no tiene la misma dirección que el vector x. −
Los vectores propios generalizados son vectores que se transforman en vectores de la misma dirección con la matriz A que con la matriz B.
−
El valor propio generalizado λ determina la diferencia en el cambio de longitud con una y otra matriz. La Figura 36 muestra un ejemplo en el que se observa la diferencia en cómo se transforma un vector cualquiera x y cómo se transforman dos vectores propios generalizados x1 y x2.
−
En algunos casos prácticos el problema generalizado surge cuando se alinean fuerzas de distinta naturaleza, por ejemplo en Mecánica cuando se alinean las fuerzas de inercia y las fuerzas elásticas.
3.8.3 Reducción del problema generalizado al problema estándar Supóngase el problema generalizado de valores y vectores propios: Ax = λ Bx
(290)
Se desea transformarlo a la forma estándar. Pre-multiplicando por la inversa de B se obtiene:
B−1Ax = λ x
(291)
Álgebra lineal numérica con Matlab
pág. 60
Este método, aunque válido, no es una forma práctica de resolver el problema porque se pierde la simetría y se hacen más complicados los cálculos. De otra forma, si la matriz B es simétrica y definida positiva se puede factorizar como B=R T R (R no singular) y se tiene:
Ax = λ Bx = λ R T Rx, AR −1y = λ R T y ⇒
y ≡ Rx x = R −1y − T −1 R AR y = λ y
(292)
Consecuencias: − Si A y B son simétricas, la nueva matriz también lo es y los valores propios serán reales. −
−
Los valores propios del problema generalizado tienen los mismos signos que los valores propios de A (por la ley de inercia de Sylvester) Si A y B son simétricas, los vectores propios correspondientes a valores propios distintos son ortogonales respecto a dichas matrices :
Ax1 = λ1Bx1 ⎫ ⎬ Ax 2 = λ2Bx 2 ⎭
−
x2T Ax1 = λ 1x2T Bx1 ⎫ ⎬ x1T Ax 2 = λ 2 x1T Bx2 ⎭
⎧ x2T Bx1 = 0 ( λ1 − λ 2 ) x Bx1 = 0 ⎨ T x Ax = 0 ⎩ 2 1 Los vectores propios se suelen normalizar respecto a la matriz B: T 2
(293)
xT iBxi = 1 −
Como consecuencia, la matriz P cuyas columnas son los vectores propios, diagonaliza por congruencia las matrices A y B simultáneamente:
PT AP = D; −
x ≡ Py
⇒
PT APy = λPT BPy
⇒
Dy = λ y
(296)
⇒
PT BPP −1 = IP −1 ⇒
P −1 = PT B
(297)
Desplazamiento de los valores propios. Si se sustituye la matriz A por A – α B, los valores propios λ pasan a ser ( λ – α) :
AP = BPD; −
(295)
Obsérvese que la matriz de vectores propios P no es una matriz ortogonal puesto que PT P≠I. Para calcular la inversa de la matriz P (que no es PT ) se puede proceder así:
PT BP = I −
PT BP = I
Introduciendo el cambio de variables x=Py y pre-multiplicando por PT la ecuación de valores y vectores propios se reduce a forma diagonal:
Ax = λ Bx, −
(294)
( A − α B ) P = BPD − α BP;
( A − α B ) P = BP (D − α I )
(298)
Los valores propios no cambian ante una transformación de congruencia aplicada a ambas matrices A y B. La ecuación de los valores y vectores propios es:
Ax = λ Bx
(299)
Introduciendo la transformación x=Cy y premultiplicando por CT
CT ACy = λ CT BCy
(300)
de donde se deduce que λ es un valor propio de las matrices CT AC y CT BC asociado con el vector propio y.
Valores y vectores propios −
−
−
pág. 61
Valores propios cuando A tiene rango menor que n. En este caso los vectores de Ker( A) son vectores propios generalizados asociados con λ =0. Esto sucede cuando la matriz A es sólo semidefinida positiva. Si es la matriz B la que tiene rango menor que n, los vectores de Ker( B) son vectores pro pios asociados con valor propio infinito: λ =∞. El papel de las matrices A y B puede intercambiarse escribiendo:
Ax = λ Bx ⇔ Bx = λ −1Ax
(301)
3.8.4 Cociente de Rayleigh para Ax= Bx El cociente de Rayleigh para el problema generalizado se define en la forma: x T Ax R(x) = T x Bx
(302)
Propiedades: −
Derivando respecto a x para imponer la condición de valor estacionario: (xT Bx )Ax − (xT Ax )Bx ∇ R(x) = 2 =0 xT Bx
−
⇒
x T Ax Ax − T Bx = 0 x Bx
(303)
de donde se deduce que los valores estacionarios se obtienen cuando x es un vector propio, en cuyo caso el cociente de Rayleigh es el valor propio correspondiente. Si las matrices A y B son simétricas y definidas-positivas o Todos los valores propios son reales y mayores que cero. o El mínimo se obtiene para el valor propio más pequeño λ 1 y el máximo para el valor
propio más grande λ n. o El caso de las matrices de rango menor que n se deduce fácilmente. −
Cálculo de valores y vectores propios. Si se conoce: o el valor propio, el vector propio se calcula resolviendo ( A – λi B)xi=0 y normalizando. o el vector propio, el valor propio se puede calcular con el cociente de Rayleigh.
−
Existen también versiones de los teoremas mini-max y maxi-min para este problema.
3.8.5 Convergencia a otros valores propios. Técnicas de deflacción matricial Para calcular el k -ésimo valor propio se itera con un vector B-ortogonal a los k –1 vectores propios ya calculados (esta condición se impone en cada paso). z k = z k − ∑ j =1 ( x T j Bz k )x j k −1
(304)
También puede obtenerse la convergencia a otros valores y vectores propios mediante una transformación de congruencia. Supóngase calculado x1 y λ 1 ; sea S una matriz tal que:
S = [ x1 s2 s3 ... s n ]
(305)
donde x1 es como se ha dicho un vector propio ya calculado, y s2, ..., sn son vectores linealmente independientes de x1 y entre sí. Teniendo en cuenta que PT AP=D y PT BP=I:
Álgebra lineal numérica con Matlab
⎡λ 1 0 0 ⎤ ⎢ ⎥ 0 ⎥ ST AS = ⎢ ⎢ ⎥ A1 ⎢ ⎥ 0 ⎣⎢ ⎦⎥
pág. 62
⎡1 0 0 ⎤ ⎢ ⎥ 0 ⎥ ST BS = ⎢ ⎢ ⎥ B1 ⎢ ⎥ 0 ⎣⎢ ⎦⎥
(306)
Resulta que las matrices A1 y B1 ampliadas tienen los mismos valores propios que A y B, pues se relacionan con éstas a través de una transformación de congruencia. Además, tienen una fila y columna menos. Con transformaciones de este tipo se pueden calcular valores propios distintos de los ya calculados.
Factorizaciones de una matriz
pág. 63
4. Factorizaciones de una matriz 4.1
Factorización LU
4.1.1 Ejemplo de factorización LU directa Para una matriz cuadrada de rango r =n la eliminación de Gauss llega al siguiente resultado (ver apartado 2.5.2):
A = LU;
⎡ a11 ⎢a ⎢ 21 ⎢ a31 ⎢a ⎣ 41
a12 a22 a32 a42
a13 a23 a33 a43
a14 ⎤ ⎡ 1 0 0 a24 ⎥⎥ ⎢⎢l21 1 0 = a34 ⎥ ⎢l31 l32 1 a44 ⎥⎦ ⎢⎣l41 l42 l43
0⎤ ⎡ u11 u12 u13 0⎥⎥ ⎢⎢ 0 u22 u23 0⎥ ⎢ 0 0 u33 1 ⎥⎦ ⎢⎣ 0 0 0
u14 ⎤ u24 ⎥⎥ u34 ⎥ u44 ⎥⎦
(307)
Exactamente a este mismo resultado se puede llegar por identificación directa de elementos, con el mismo nº de operaciones aritméticas, pero con una secuencia más favorable. Considérense los pasos siguientes: 1. A partir de la primera fila de A:
u1 j = a1 j 2.
u2 j = a2 j − l21u1 j = a2 j −
u12u1 j u11
j = 2,..., n
(310)
l j2 =
1 (a − l u ) u22 j 2 j1 12
j = 3,..., n
(311)
Para la tercera fila y columna de A: −
Elemento de la diagonal:
a33 = l31u13 + l32 u23 + 1 ⋅ u33 −
−
u33 = a33 − l31 u13 − l32 u23
(312)
Resto de los elementos de la 3ª fila:
a34 = l31u14 + l32 u24 + u34
u34 = a34 − l31 u14 − l32 u24
(313)
Resto de los elementos de la 3ª columna: a43 = l41u13 + l42 u23 + l43 u33
6.
(309)
A partir de la segunda columna de A, los elementos de la segunda columna de L:
a j 2 = l j1u12 + l j 2u22 ; 5.
j = 2,3,...n
l j1 =a j1 / u11
A partir de la segunda fila de A, los elementos de la segunda fila de U: a2 j = l21u1 j + 1 ⋅ u2 j ;
4.
(308)
A partir de la primera columna de A:
a j1 = l j1u11 ; 3.
j = 1,2,..., n
l43 =
1 (a − l u − l u ) u33 43 41 13 42 23
(314)
Para el último elemento de la diagonal: a44 = l41u14 + l42 u24 + l43 u34 + 1 ⋅ u44
u44 = a44 − l41 u14 − l42 u24 − l43 u34
(315)
Álgebra lineal numérica con Matlab
pág. 64
4.1.2 Fórmulas generales para la factorización LU La Figura 37 muestra esquemáticamente la factorización directa A=LU en un estado intermedio. Se supone que, utilizando los elementos de las partes rayadas de A, se han calculado ya las correspondientes partes rayadas de las matrices L y U. En el siguiente paso se van a calcular la columna i de L y la fila i de U, utilizando las correspondientes fila y columna de la matriz A (marcadas con línea gruesa en la Figura). Hay que recordar que los elementos de la diagonal de L son unos.
0
=
0 Figura 37. Estado intermedio de la factorización A=LU.
Para calcular los elementos marcados (columna i de L y fila i de U, a partir de la diagonal) se procede del siguiente modo: 1. A partir de aii se calcula el elemento de la diagonal uii: i −1
aii = ∑ lik uki + uii k =1
⇒
i −1
uii = aii − ∑ lik uki
(316)
k =1
2. A partir de aij ( j=i+1 , …, n) se calculan los elementos de la fila i de U: i −1
aij = ∑ lik ukj + uij
⇒
k =1
i −1
uij = aij − ∑ lik ukj , j = i + 1,..., n
(317)
k =1
3. A partir de a ji ( j=i+1 , …, n) se calculan los elementos de la columna i de L: i −1
a ji = ∑ l jk uki + l jiuii k =1
⇒
i −1 1⎛ ⎞ l ji = ⎜ a ji − ∑ l jk uki ⎟ , j = i + 1,..., n uii ⎝ k =1 ⎠
(318)
Las operaciones indicadas por las ecuaciones (316)-(318) presentan las siguientes ventajas respecto a la eliminación de Gauss: 1. Están fundamentalmente basadas en el producto escalar de vectores (fila de L por columna de U), operación que se realiza de un modo muy eficiente en un computador. 2. Los elementos de las matrices L y U alcanzan su valor definitivo de una vez, sin resultados intermedios. De esta forma se minimiza el trasiego de datos entre el procesador y los distintos niveles de la memoria del ordenador. Téngase en cuenta que en los ordenadores modernos el acceso a memoria es tan costoso o más que las propias operaciones aritméticas. En relación con lo que se acaba de decir, todos los datos que se utilizan a la derecha de la igualdad en las expresiones (316)-(318) tienen ya su valor definitivo, por lo que pueden ser accedidos con operaciones de sólo-lectura y hacer un uso muy eficiente de las memorias intermedias del procesador (memorias cache).
4.1.3 Factorización LU con matrices simétricas Si la matriz A es simétrica las expresiones de l ij y de u ji coinciden, excepto en que los elementos l ij están divididos por los uii. Se puede escribir entonces la factorización en la forma:
Factorizaciones de una matriz
pág. 65
(319) A = LU = LDLT Hay dos posibilidades para almacenar los resultados de la factorización: 1. Almacenar U en la matriz resultado en su posición habitual, incluyendo su diagonal. Se puede recuperar L fácilmente dividiendo las filas de U por el elemento de la diagonal y trasponiendo. 2. Almacenar D en la diagonal de la matriz resultado y LT en la mitad superior. En este caso, las filas de U se guardan ya divididas por el elemento de la diagonal. Los cálculos pueden hacerse por filas o por columnas. A continuai ción se desarrollará con detalle el cálculo por filas, basándose en la representación simbólica de la matriz de la Figura 38. En el cálculo por filas se supone que la parte rayada está ya calculada. Se van a calcular los elementos de la fila i de LT y el elemento de la diagonal de D (o de U). Se parte de las expresiones (316)(318) y se tiene en cuenta la relación existente entre los elementos de U y los de LT y D. En primer lugar se calcula el elemento de la diagonal d ii, distinguiendo según las dos formas de almacenar los resultados: i −1 2 uki u dii = aii − ∑ lik uki = aii − ∑ uki =aii − ∑ ki k =1 k =1 ukk k =1 u kk i −1
i −1
i −1
i −1
(320) Figura 38. Factorización LDLT de una matriz simétrica.
2
dii = aii − ∑ l l d = aii − ∑ ( lkiT ) d kk T T ki ki kk
k =1
(321)
k =1
El cálculo de los elementos de la fila i (con línea gruesa en la Figura 38) recuerda al producto escalar de la columna i por la columna j, excepto que cada sumando está dividido (o multiplicado) por el elemento de la diagonal. Estas divisiones restan eficiencia a los cálculos. Calculando U: uki ukj , u k =1 kk
i −1
i −1
uij = aij − ∑ lik ukj = aij − ∑ k =1
j = i + 1,..., n
(322)
Si se almacena LT , es decir, los elementos de U divididos por el elemento de la diagonal: i −1 i −1 1⎛ ⎞ 1⎛ ⎞ l = ⎜ aij − ∑ lik ukj ⎟ = ⎜ aij − ∑ lkiT lkjT dkk ⎟ , dii ⎝ k =1 k =1 ⎠ d ii ⎝ ⎠ T ij
j = i + 1,..., n
(323)
El algoritmo puede modificarse para resolver la dificultad de tener que dividir (o multiplicar) por el elemento de la diagonal, al mismo tiempo que se vectoriza. La clave está en construir un vector auxiliar v cuyos elementos son los de la columna i de LT multiplicados cada uno de ellos por el correspondiente elemento de la diagonal. Las expresiones (321) y (323) se transforman en: i −1
i −1
d ii = aii − ∑ l l d kk =aii − ∑ lkiT vk k =1
T T ki ki
k =1
i −1 i −1 1⎛ ⎞ 1⎛ ⎞ T T l = ⎜ aij − ∑ lki lkj dkk ⎟ = ⎜ aij − ∑ lkjT vk ⎟ , dii ⎝ k =1 k =1 ⎠ d ii ⎝ ⎠ T ij
(324) j = i + 1,..., n
(325)
La factorización LU, almacenando D y LT sobre la propia matriz A y utilizando las ecuaciones (324) y (325), puede programarse en Matlab como se muestra a continuación:
Álgebra lineal numérica con Matlab
pág. 66
function [D,Lt]=LUsim(A)
n=si ze( A, 1) ; % pr i mer a f i l a de L ' for j=2:n
A( 1, j ) =A( 1, j ) / A( 1, 1) ; end for i=2:n
% se f or ma el vector auxi l i ar v for k=1:i-1;
v( k) =A( k, i ) * A( k, k) ; end
% cál cul o del el ement o di agonal for k=1:i-1
A( i , i ) =A( i , i ) - v( k) * A( k, i ) ; end
% c ál c ul o del r e st o de l a f i l a i for j=i+1:n for k=1:i-1
A( i , j ) =A( i , j ) - v( k) * A( k, j ) ; end
A( i , j ) =A( i , j ) / A( i , i ) ; end end
D=di ag( A) ; Lt =eye( n) +t r i u( A, 1) ;
4.1.4 Factorización LU vectorizada Al igual que en el apartado anterior, el resultado contendrá la matriz D y la matriz LT almacenadas sobre la matriz A original. Las expresiones utilizadas son las mismas que en el programa anterior. El programa de factorización puede ser como sigue (sólo dos bucles for ): function [d,Lt]=LUsimVect1(A)
% se guar dan D y L' ( A=L*D*L' ) n=si ze( A, 1) ; v=zer os( n, 1) ; d=zer os( n, 1) ; d( 1) =A( 1, 1) ; % s e di vi de l a f i l a 1 por A( 1, 1) A( 1, 2: n) =A( 1, 2: n) / A( 1, 1) ; for i=2:n
% vec t or auxi l i ar v(1:i-1)=A(1:i-1,i).*d(1:i-1);
% el ement o di agonal f i l a i di i =A( i , i ) ; dii=dii-v(1:i-1)'*A(1:i-1,i);
d( i ) =di i ; A( i , i ) =di i ; % c ál c ul o del r e st o de l a f i l a i for j=i+1:n A(i,j)=(A(i,j)-v(1:i-1)'*A(1:i-1,j))/dii; end end
Lt=eye( n) +t r i u( A, 1) ;
La vectorización del algoritmo de factorización directa puede llegar aún más lejos: El siguiente programa de Matlab calcula la factorización LU utilizando doble vectorización (con un solo bucle for ): function [d,Lt]=LUsimVect2(A)
% se guar dan D y L' ( A=L*D*L' ) n=si ze( A, 1) ; v=zer os( n, 1) ; d=zer os( n, 1) ; d( 1) =A( 1, 1) ; % pr i mer a f i l a de L ' A( 1, 2: n) =A( 1, 2: n) / A( 1, 1) ;
Factorizaciones de una matriz
pág. 67
% r e st ant es f i l as f or i =2: n % vector auxi l i ar : col umna i di vi di da di agonal v( 1: i - 1) =A( 1: i - 1, i ) . * d( 1: i - 1) ; % el ement o de l a di agonal di i =A( i , i ) ; di i =di i - v( 1: i - 1) ' * A( 1: i - 1, i ) ; A( i , i ) =di i ; d( i ) =di i ; % r es t o de l a f i l a i A( i , i +1: n) =( A( i , i +1: n) - v( 1: i - 1) ' * A( 1: i - 1, i +1: n) ) / di i ; end Lt=eye( n) +t r i u( A, 1) ;
Para ofrecer una idea de la eficiencia de la vectorización y de los propios ficheros *.m de Matlab, se ha construido y ejecutado una función comparativa que llama a las funciones anteriores y a la función lu() de Matlab. Los resultados obtenidos con un procesador Intel Pentium 4 a 3 Ghz han sido los siguientes (tiempo en segundos): n
LUsim 0.2190 1.5940 12.2810
250 500 1000
LUsimVect1 LUsimVect2 0.5470 0.0470 2.4840 0.3280 11.4380 3.4220
lu() 0.0160 0.0620 0.4060
Hay que tener en cuenta que se ha utilizado la versión 6.5 de Matlab, que incorpora un JIT (Just In-Time) Accelerator . Los resultados de esta tabla muestran que el acelerador JIT consigue que el código normal sea más rápido que la vectorización simple, aunque no tan rápido como la vectorización doble. La función lu() de Matlab, escrita en Fortran ó C y compilada, es en cualquier caso mucho más rápida que cualquiera de los ficheros *.m.
4.1.5 Factorización de Choleski Esta factorización se aplica a matrices simétricas y definidas positivas, pues hace falta que esté garantizado el que todos los pivots sean positivos. Está basada en la siguiente expresión:
A = LU = LDLT = LD1/ 2D1/ 2LT = LLT Si se trabaja con L y con la parte inferior de A, las expresiones son las siguientes: − Elemento de la diagonal: i
i −1
k =1
k =1
aii = ∑ lik lkiT = ∑ lik2 + lii2 −
(326)
1
⇒
i −1 2 ⎛ 2 ⎞ lii = ⎜ aii − ∑ lik ⎟ k =1 ⎝ ⎠
(327)
Restantes elementos de la columna i: i
i −1
a ji = ∑ l l = ∑ l jk lik + l ji lii T jk ki
k =1
k =1
⇒
i −1 1⎛ ⎞ l ji = ⎜ a ji − ∑ l jk lik ⎟ l ii ⎝ k =1 ⎠
(328)
La factorización de Choleski requiere un código muy similar al visto anteriormente.
4.2
Factorización QR
4.2.1 Factorización QR por Gram-Schmidt El método de ortogonalización de Gram-Schmidt fue ya considerado en el apartado 1.2.4. Ahora se va a utilizar para pasar de una matriz cuadrada A, de tamaño n×n, cuyas columnas son independientes, a una matriz ortogonal Q:
Álgebra lineal numérica con Matlab
A = [a1
a2
pág. 68
⇒
an ]
Q = [q1
q2
qn ]
(329)
Recuérdese que L [a1 , a 2 ,..., a k ] = L [q1 , q2 ,..., q k ] , es decir, que las k primeras columnas de A general el mismo subespacio que las k primeras columnas de Q. Según esto, cada vector ai sólo tiene componentes en los i primeros vectores ortonormales q, es decir:
a1 = ( q1T a1 ) q1 a2 = ( q1T a2 ) q1 + (q 2T a2 ) q2
(330)
a3 = ( q1T a3 ) q1 + (q2T a3 ) q2 + (q3T a3 ) q3 Las ecuaciones (330) indican que las matrices A y Q se relacionan a través de una matriz triangular superior R , que es cuadrada n×n e invertible:
[a1 a2
an ] = [ q1 q2
⎡q1T a1 q1T a2 ⎢ 0 qT a 2 2 qn ] ⎢ ⎢ ⎢ 0 ⎢⎣ 0
q1T a n ⎤ q T 2 a n ⎥⎥ ⎥ ⎥ qT nan ⎥⎦
A = QR
(331)
La factorización QR suele aplicarse de modo que todos los elementos de la diagonal de R sean positivos (la función qr() de Matlab sin embargo, no cumple esta condición). La fila i de Q representa las componentes de todas las columnas de A según el vector qi de la base ortonormal. La columna j de R representa las componentes de a j en los vectores de la base ortonormal. Como se verá a continuación, esta factorización mantiene pleno sentido aunque el número de columnas de A sea inferior al número de filas. El método de Gram-Schmidt se suele aplicar según el llamado algoritmo modificado, que utiliza una secuencia de operaciones diferente a la convencional. En el algoritmo modificado, tan pronto como se halla un vector de la base ortonormal qi se hacen ortogonales respecto a él todos los vectores a j, j=i+1,...,n. Las expresiones correspondientes son las siguientes:
q1 = a1 a1
a j = a j − (q1T a j ) q1
j = 2,3,..., n
q 2 = a 2 a2
a j = a j − ( qT 2 a j ) q 2
j = 3,4,..., n
a j = a j − ( qT i a j ) qi
j = i + 1,..., n
...
q i = a i ai
(332)
4.2.2 Factorización QR de matrices rectangulares Considérese una matriz rectangular A, de tamaño m×n (m>n). Se puede obtener una matriz Q, tam bién m×n, cuyas columnas son ortonormales y generan el mismo subespacio de R m que las de A. Las fórmulas anteriores (330) y (331) se aplican también en este caso. La matriz Q es m×n, pero la matriz R es cuadrada n×n. A esta factorización se le llama factorización QR incompleta, porque las columnas de Q no constituyen una base completa de R m (en este caso Q no es una matriz ortogo nal , sino lo que se llamó matriz de columnas ortogonales). La factorización QR completa obtiene una matriz Q ortonormal (m×m) aunque la matriz A sea rectangular, añadiendo a la matriz R tantas filas de ceros como columnas se han añadido a Q (la Figura 39 representa gráficamente la factorización QR completa de una matriz rectangular):
Factorizaciones de una matriz
[a1 a2
pág. 69
an ] = [q1 q2
q n q n +1
⎡q1T a1 q1T a2 ⎢ T ⎢ 0 q 2 a2 ⎢ ⎢ qm ] ⎢ 0 0 ⎢ 0 0 ⎢ ⎢ ⎢⎢ 0 0 ⎣
q1T a n ⎤ ⎥ qT 2 a n ⎥ ⎥ ⎥ qT n a n ⎥ 0 ⎥ ⎥ ⎥ 0 ⎥⎥⎦
(333)
R A
Q
=
0 Figura 39. Factorización QR completa de una matriz rectangular.
La factorización QR se puede utilizar para resolver sistemas de ecuaciones lineales. En sí es más cara que la factorización LU, pero es más estable y es ventajosa siempre que por algún motivo se disponga ya de ella. Con la factorización QR la resolución de un sistema de ecuaciones lineales se reduce a un producto por QT y a una vuelta atrás con una matriz triangular superior:
Ax = b
A = QR
QRx = b
Rx = y ⎫ ⇒ y = QT b ⇒ x = R −1y ⎬ Qy = b ⎭
(334)
En el caso de las ecuaciones normales propias del método de los mínimos cuadrados su aplicación es mucho más ventajosa, porque las ecuaciones normales suelen tener números de condición altos. En este caso se llega a:
Ax = b
A T Ax = A T b
A = QR
R T QT QRx = R T Q T b
(335)
y como Q es ortogonal y R triangular superior e invertible, la ecuación final es muy fácil de resolver
RT Rx = R T QT b
⇒
Rx = Q T b
(336)
A continuación se incluye una función de Matlab para la factorización QR incompleta mediante el método de Gram-Schmidt modificado function [Q,R]=QRgsmod(A)
% Fact or i zaci ón QR ( G. - S. modi f i cado) [ m, n] =si ze( Q) ; R=zeros( n, n) ; for k=1:n
% nor ma de l a col umna k R( k, k)=nor m( Q( : , k)) ; % se normal i za l a col umna k Q( : , k) =Q( : , k) / R( k, k) ; % Se el i mi na l a col umna k de l as si gui ent es col umnas for j=k+1:n
R( k , j ) =Q( : , k) ' * Q( : , j ) ; Q( : , j ) =Q( : , j ) - Q( : , k) * R( k , j ) ; end
end
Álgebra lineal numérica con Matlab
pág. 70
4.2.3 Factorización QR mediante matrices de Householder La factorización QR puede hacerse también por medio de transformaciones de Householder, de modo aún más estable que con Gram-Schmidt. Además este método es más adecuado cuando se desea obtener la factorización QR completa. En el apartado 1.7.5 se vio que pre-multiplicando por matrices de Householder se pueden hacer ceros por debajo de la diagonal principal en todas las columnas de una matriz rectangular A: U n −1...U 2U1A = R
(337)
Cada matriz U j tiene las ( j –1) primeras filas y columnas iguales a las de la matriz identidad Im. Los restantes elementos contienen una matriz de simetría de Householder que responde a la expresión:
vvT H = I − 2 T = I − β vvT , v = x − x e1 (338) v v siendo x el vector cuyas componentes 2, 3, ..., n se quiere anular. Pre-multiplicando la expresión (337) por las inversas de las matrices de simetría Ui y teniendo en cuenta que el producto de matrices unitarias es también unitaria: A = U1−1U 2−1...U n−1−1R = U 1U2 ...U n −1R = QR,
Q ≡ U1 U2 ...U n −1
(339)
No es necesario formar explícitamente las matrices Ui y menos aún realizar los productos de matrices que aparecen en las expresiones anteriores. La matriz Hi se puede calcular mediante la expresión (338) a partir del vector v cuando se necesite, pero para hacer un producto HA basta disponer del vector v, como se explicó en el apartado 1.7.6: T
HA = ( I − β vvT ) A = A − β vv T A = A − β v ( A T v )
Por otra parte, el vector v se puede calcular a partir de x mediante la función householder.m, explicada en el apartado 1.7.5. Con este método, aunque la matriz A sea rectangular (m>n), se obtiene la factorización QR completa. A continuación se muestra un programa de Matlab para hacer la factorización QR por Householder. % QR por el mét odo de Househol der function [Q,R]=QRhouse(A)
[ m, n] =si ze( A) ; i f m
% cr ear el vect or de Househol der if j
[ v, bet a( j ) ] =hous ehol der ( A( j : m, j ) ) ; % actual i zar l a mat r i z H A( j : m, j : n) =( A( j : m, j : n) - bet a( j ) * v ( 1: m- j +1) * ( v ( 1: m- j +1) ' * A( j : m, j : n) ) ) ; % al macenar v en l os cer os hechos en l a mat r i z A A( j +1: m, j ) =v( 2: m- j +1) ; end end
% Q y R no se guar dan expl í ci t ament e % Par a obt ener Q se puede mul t i pl i car I ( n, n) por l a secuenci a de t r ansf ormaci ones % Q=H1*H2*. . . *Hn- 1*I % Q=I ; Q=Hj *Q; j =n- 1: - 1: 1
Factorizaciones de una matriz
pág. 71
for j=n:-1:1
% se ext r ae v de H, añadi endo el "1" v( j : m) =[ 1; A( j +1: m, j ) ] ; % se apl i ca l a f ór mul a H*Q=Q- v*w' ; % donde w=bet a*Q' *v; Q( j : m, j : m) =( eye( m- j +1) - bet a( j ) *v(j : m) *v(j : m) ' ) *Q( j : m, j : m) ; end
% R est á en l a part e superi or de A R=[ t r i u( A( 1: n, 1: n) ) ; z er o s ( m- n, n) ] ;
4.2.4 Factorización QR mediante rotaciones de Givens Las matrices de simetría de Householder permiten hacer ceros en una columna, debajo de la diagonal, en un solo producto matricial. Las matrices de rotación de Givens, vistas en el apartado 1.6.3, sólo introducen un cero en cada producto matricial. Las matrices de Givens son también ortogonales y más sencillas que las de Householder, pero en general no resultan rentables para hacer en muchos pasos los ceros que las matrices de Householder hacen en uno solo. Sin embargo, las matrices de Givens se aplican con ventajas respecto a las de Householder en todos aquellos casos en los que se trata de hacer ceros selectivamente, es decir, en sólo unos pocos elementos debajo de la diagonal. Tal sucede por ejemplo en las matrices tridiagonales o con forma de Hessenberg, como las que se muestran en la Figura 40. Para anular los elementos de debajo de la diagonal de estas matrices hay que aplicar ( n –1) rotaciones de Givens, exactamente el mismo número que simetrías de Householder.
⎡* ⎢ ⎢* ⎢0 ⎢ ⎢0 ⎢ ⎢⎣ 0
*
0
*
*
*
*
0
*
0
0
0 0⎤ ⎥ 0 0⎥ * 0 ⎥⎥ , * *⎥ ⎥ * * ⎥⎦
⎡* ⎢ ⎢* ⎢0 ⎢ ⎢0 ⎢ ⎢⎣ 0
*
*
*
*
*
*
0
*
0
0
* *⎤ ⎥ * *⎥ * *⎥⎥ * *⎥ ⎥ * *⎥⎦
Figura 40. Matrices simétrica y con forma de Hessenberg.
La factorización QR se deduce fácilmente a partir de las rotaciones de Givens que convierten A en una matriz triangular superior R . En general: G n,n −1 G n 2 G 32 G n1 G 31G 21 A = R
⇒
T T T T T A = GT21 G31 G n1 G32 Gn2 Gn, n −1 R = QR
T T T G n 2 G n , n −1 Q ≡ G T21G T31 G Tn1G32
4.3
(340)
Descomposición espectral de matrices normales
Sea A ∈ R n×n una matriz normal. Toda matriz normal –por tanto unitariamente diagonalizable– admite una descomposición espectral en la forma:
A = UDU H = [u1
⎡λ 1 u n ] ⎢⎢ ⎢⎣
⎤ ⎡u H 1 ⎤ ⎥ ⎢ ⎥ = n λ u u H k k k ⎥ ⎢ H ⎥ ∑ k =1 λ n ⎥⎦ ⎢⎣u n ⎥⎦
(341)
En el caso de que haya valores propios múltiples, es posible agrupar en una matriz Pi las mi matrices de rango 1 asociadas con el valor propio λ i . En este caso, la descomposición espectral se expresa en la forma:
Álgebra lineal numérica con Matlab n
pág. 72 m1
A = UDU = ∑ λk u u = λ1 ∑ u u + λ2 T
T k k
k =1
T k k
k =1
= λ1P1 + λ2 P2 + ... + λ s Ps
m1 + m2
∑uu
T k k
n
+ ... + λ s
k = m1 +1
∑
u ku Tk =
k = n − m s +1
(342)
( m1 + m2 + ... + ms = n )
Las matrices Pi son matrices de proyección ortogonal sobre los subespacios propios , pues son simétricas e idempotentes. La simetría es evidente, pues cada Pi es una suma de matrices simétricas. La idempotencia se demuestra también fácilmente:
Pi2 = ∑ i = p ui uTi ∑ j = p u j uTj = ∑ i = p ∑ j = p ui uiT u j u Tj = ∑ i = p ∑ j = p ui δ ij u Tj = ∑ i = p ui uiT = Pi q
q
q
q
q
q
q
(343)
Se verifican además las siguientes propiedades:
Pi Pj = ∑ m= p u muTm ∑ k = r u k uTk = ∑ m = p ∑ k = r u mu Tmu ku Tk = q
s
q
= ∑ m= p ∑ k =r u m ( u u ) u = ∑ m = p ∑ k =r u q
(∑
s
T m k
q
T k
s
s
T mδ mk k
)
q
q
(344)
Pj = I
(345)
u = 0 ( m ≠ k siempre )
T T = = P x u u x u u ( ∑ ∑ j i i i i x ) = ∑ i = p u iα i = x j =1 i= p i=p
s
q
⇒
∑
s j =1
La expresión (345) indica que la suma de todas las matrices Pi es la matriz identidad I, pues cualquier vector x se puede expresar como suma directa de sus componentes en la base ui.
4.4
Descomposición de valores singulares (DVS)
Para todas las matrices normales n×n existe una matriz unitaria U tal que A se puede factorizar mediante la descomposición espectral (ver apartado 4.3):
A = UDU H
(346)
donde las columnas de U son los vectores propios de A. Si la matriz A no es cuadrada, no tiene vectores propios ni descomposición espectral. Sin embargo, para toda matriz real m×n siempre existe una factorización análoga a la espectral en la forma:
A = UΣV H
(347)
donde U y V son matrices unitarias (ortogonales, en casi todas las aplicaciones de este curso) de tamaño m×m y n×n, respectivamente, y es una matriz diagonal de tamaño m×n cuyos elementos son tales que σ ii>0 (i=1,2,...,r ) y σ ii=0 (i>r ), siendo r el rango de la matriz A. Es habitual considerar los valores singulares ordenados de mayor a menor: σ1 ≥ σ 2
≥ ... ≥ σ r > σ r +1 = ... = σ p = 0, p = min ( m, n )
(348)
en cuyo caso las columnas de U y V deben estar ordenadas de modo acorde. La forma de estas matrices se puede ver gráficamente en la Figura 41, que distingue los casos en que A tenga más o menos filas que columnas.
A
=
U
V H
A
=
U
V H
Figura 41. Descomposición de valores singulares de una matriz con menos y más filas que columnas.
Factorizaciones de una matriz
pág. 73
A esta factorización se le llama "Descomposición de Valores Singulares o DVS" (Singular Value Decomposition o SVD, en inglés) y tiene aplicaciones muy importantes en ingeniería, tanto teóricas como prácticas. En lo sucesivo se considerarán exclusivamente matrices A reales, por lo que se hará referencia a matrices ortogonales.
4.4.1 Cálculo y existencia de la DVS Sea A ∈ R m×n . La prueba de la existencia de la DVS va a ser de tipo constructivo, es decir, se demostrará que existe definiendo un proceso que permite calcularla. La matriz AT A es simétrica y –al menos– semi-definida positiva. Por tanto tiene n vectores propios ortogonales vi que pueden formar las columnas de la matriz ortogonal V ∈ R n×n :
AT Av i = λ i vi
( AT A ) V = VDn
(349)
Los valores propios λ i son no-negativos. Multiplicando la expresión (349) por vT i , λ i se puede ex presar en la forma (cociente de Rayleigh):
v Ti AT Av i = λi viT v i =λi
λ i = viT A T Av i
2
= Av i
(350)
Si A tiene rango r la matriz AT A tendrá también rango r (pues Ker(A)=Ker(AT A)). Habrá r valores propios λ i mayores que cero y n–r valores propios nulos. Se definen ahora el valor singular σ j y el vector singular u j (u j∈R m) en la forma: σ j
≡+
λj ,
u j ≡ Av j σ j ;
j = 1, 2,..., r
(351)
Los vectores u j así formados constituyen un sistema ortonormal en el espacio R m. En efecto:
uT j u j = vTj AT Av j σ 2j = λ j λ j = 1, j = 1,2,..., r uTi u j = vTi AT Av j σ iσ j = λ j viT v j σ iσ j = 0, i ≠ j
(352)
Los r vectores ortogonales u j se pueden completar por el teorema de la base incompleta (y el método de Gram-Schmidt) con otros (m – r) vectores ortogonales, hasta constituir una base de R m. Este sistema ampliado de m vectores ortogonales u j constituyen las columnas de la matriz ortogonal U ∈ R m×m . Se pretende demostrar la existencia de la DVS ( A=U VT ), pero para ello se estudiará en primer lugar el producto UT AV. El elemento (i, j) de este producto matricial tiene como valor:
( UT AV )ij = uiT Av j = uiT u jσ j =δijσ j =0
si j ≤ r , por la definición de u j si j > r , pues Av j = 0
(353)
En definitiva, el producto UT AV es igual a una matriz , de tamaño m×n, que sólo tiene distintos de cero los r primeros elementos de la diagonal. Estos elementos son los valores singulares σ j. Se ha demostrado pues que existen unas matrices ortogonales U y V que diagonalizan a la matriz A:
UT AV = Σ
(354)
Como las matrices U y V son ortogonales y por tanto invertibles, se puede escribir:
A = UΣVT = ∑ i =1
min ( m , n )
σ i ui vT i
(355)
Álgebra lineal numérica con Matlab
pág. 74
que es la expresión de la descomposición de valores singulares (DVS). Los vectores v j y u j son, respectivamente, los vectores singulares por la derecha y por la izquierda de la matriz A:
Av j = σ j u j , A H u j = σ j v j
j = 1,2,..., r
(356)
Como no se ha presupuesto ninguna condición para la matriz A, la DVS existe para cualquier ma triz rectangular, cualquiera que sea su rango r . La DVS tiene propiedades muy importantes, derivadas de la ortogonalidad de las matrices U y V, y del carácter no negativo de los valores singulares σ j.
4.4.2 Propiedades de las matrices U y V Las matrices ortogonales U y V que aparecen en la descomposición de valores singulares A=U VT tienen –entre otras– las siguientes propiedades: 1. Las columnas de U son vectores propios de AAT
AA T = UΣVT VΣT UT = UΣΣT UT = UD mUT
(357)
(358)
2. Las columnas de V son vectores propios de AT A
A T A = VΣT UT UΣVT = VΣT ΣVT = VD n VT
3. Tanto AAT como AT A tienen los mismos valores propios no nulos , lo que deduce de la relación de Dm y Dn con T y T , respectivamente. 4. Relación de las matrices U y V con los subespacios de A: o Las columnas 1 a r de U son una base ortonormal de Im( A) o Las columnas r +1 a m de U son una base ortonormal de Ker( AT ) o Las columnas 1 a r de V son una base ortonormal de Im( AT ) o Las columnas r +1 a n de V son una base ortonormal de Ker( A)
Las propiedades anteriores se deducen de las relaciones siguientes: o Ker(AT A)=Ker(A), pues si Ax=0 también se verifica AT Ax=0 y si AT Ax=0 se verifica que xT AT Ax=(Ax)T Ax=0, luego Ax=0. o Im(AT A)=Im(AT ), pues ambos son los espacios ortogonales complementarios de Ker(AT A) y Ker(A) en R n. o Análogamente, Ker( AAT )=Ker(AT ), Im(AT A)=Im(A) y Im(AAT )=Im(A)
4.4.3 Cálculo de la DVS Las propiedades 1. y 2. del apartado anterior no son suficientes por sí mismas para determinar las matrices U y V, puesto que los vectores propios de una matriz no están unívocamente determinados (por ejemplo, si u j es vector propio de AAT , el vector – u j también lo es; si hay valores propios múltiples lo que está determinado es el subespacio propio correspondiente, pero no una base ortonormada de dicho subespacio). Las matrices U y V deben calcularse de modo que sean compatibles entre sí , según la expresión (351), que también se puede escribir en la forma: AV = UΣ (359) Para calcular las matrices U, V y se puede proceder del siguiente modo: 1. Se calculan las matrices V y resolviendo el problema de valores y vectores propios:
Factorizaciones de una matriz
pág. 75
A T AV = VDn ,
σ i
= d ii , i = 1,2,..., r
(360)
2. Ahora se calcula una matriz U compatible con la V ya calculada por medio de la ecuación (359). Esta ecuación sólo permite calcular directamente las r primeras columnas de U, pues son las que están multiplicadas por un valor singular distinto de cero. Las m – r columnas restantes se calculan simplemente completando una base ortonormal de R m. De forma análoga, es también posible calcular en primer lugar las matrices U y resolviendo el problema de valores y vectores propios AAT U = UDm , y luego calcular V de modo compatible con la U ya calculada por medio de la ecuación UT A = ΣV T (o bien, trasponiendo: AT U = VΣT ). La DVS de la matriz A no es única. Está claro que los valores singulares, ordenados de mayor a menor, sí están unívocamente determinados. Los vectores singulares por la derecha vi, correspondientes a valores singulares simples, no están unívocamente determinados pues pueden ser multiplicados por un factor escalar de módulo unidad (–1, en el caso real). En el caso de valores singulares múltiples, lo que está unívocamente determinado es un subespacio singular de dimensión igual a la multiplicidad del valor singular. Sin embargo, dentro de este subespacio se puede tomar cualquier base ortonormal de vectores singulares.
4.4.4 Aplicaciones de la DVS Una de las aplicaciones más importantes de la DVS es el cálculo del rango de la matriz A. El producto por matrices unitarias conserva el rango, luego se tendrá: rango ( A ) = rango ( )
(361)
lo que indica que el rango de una matriz es igual al número de valores singulares no nulos. En un próximo apartado se verá cómo tratar el caso de valores singulares distintos de cero pero muy pequeños. Si la matriz A es cuadrada e invertible, los valores singulares serán todos no nulos y positivos. La matriz inversa de A se puede expresar como:
A −1 = V
UT = ∑ i =1 n
1
v iuT i
σ i
(362)
Análogamente, la solución de un sistema de ecuaciones lineales compatible y determinado Ax=b se puede expresar en la forma:
Ax = b
⇒
−1
x=A b=V
U b = ∑ i =1 n
T
1 σi
v u b = ∑ i =1 n
T i i
uT i b σ i
vi
(363)
La DVS conserva las normas espectral y de Frobenius , como se demuestra a continuación: 2
A 2 = ρ ( AT A ) = ρ ( V T UT U V T ) = ρ ( V
T
2
A F = traza ( A H A ) = traza ( V T U H U V H ) = traza ( V
VT ) = ρ ( T
T
2 2
)=
V H ) = traza (
T
= σ 12
)=
2 F
(364)
= ∑ σ i2 (365)
La DVS se puede aplicar también al cálculo de la condición numérica de una matriz, que para una matriz cualquiera se definió como el cociente entre el vector unitario que más "crece" al multiplicarse por la matriz A dividido por el vector unitario que menos crece: κ ( A ) = max u =1
Au min Au u =1
(366)
Álgebra lineal numérica con Matlab
pág. 76
Considerando la norma euclídea, hay que tener en cuenta que el producto por una matriz ortogonal no cambia la norma euclídea, por lo que las normas máxima y mínima de la expresión (366) son respectivamente el máximo y mínimo valor singular. Por ejemplo, para el máximo de la norma: T
max Au 2 = max UΣV u 2 = max u 2 =1
u 2 =1
Σ
u 2 =1
v = VT u
T
Vu
2
= max
v 2 = σ 1e1 2 = σ 1
Σ
v 2 =1
(367)
Por tanto, el número de condición espectral será: κ 2 ( A ) = σ 1 σ p ,
p = min ( m, n )
(368)
La condición numérica de una matriz unitaria u ortogonal es la unidad, pues sus valores singulares son los valores propios de UT U=UUT =I, que son todos la unidad. Operar con estas matrices no am plifica los errores en los datos.
4.4.5 Valores singulares y perturbaciones en una matriz En este apartado se va a considerar la sensibilidad de los valores singulares ante perturbaciones de la matriz. Evaluar el rango de una matriz puede ser muy difícil, pues el rango es un valor entero, mientras que los elementos de la matriz pueden varíar de modo continuo. Para resolver esta dificultad se introduce el concepto de rango numérico de una matriz. Sea ε un número muy pequeño. Se define el rango numérico o rango- de la matriz A como: rango ( A, ε ) ≡ min rango ( A )
(369)
A − A <ε 2
Es importante interpretar correctamente esta definición: Las matrices A son todas las posibles matrices que están a una distancia de A menor que ε , medida en la norma espectral. El rango numérico de A se define como el rango de aquella matriz que menor rango tiene y que está a una distancia de A menor que ε .
Teorema 1: Sea A∈R m×n una matriz cuya DVS es A = U V T = ∑ i =1σ i ui viT . Se define la matriz Ak r
como una matriz de rango k formada por la suma A k = ∑ i =1σ iui vT i . Se cumple que la matriz de rango k que mejor aproxima a A en norma espectral es Ak , es decir: k
min
rango ( A ) = k
A − A 2 = A − A k 2 = σ k +1
(370)
y además el error en la aproximación es el valor singular σ k+ 1, que es el mayor valor singular no incluido en Ak . Demostración: Utilizando la DVS y la definición de Ak se tiene que: U H A k V = diag (σ 1 ,..., σ k , 0,...0 ) , rango ( A k ) = k
U H ( A − A k ) V = diag ( 0,...,0,σ k +1 ,...,σ p ) ,
A − Ak
2
= σ k +1
(371)
(372)
Hay que demostrar que Ak es la matriz de rango k que mejor aproxima a A. Supóngase una matriz B∈R m×n tal que rango(B)=k . El núcleo de B tendrá dimensión (n – k ) y podrá ser generado por (n– k ) vectores independientes: Ker ( B ) = L [ x1 , x 2 ,..., x n − k ]
(373)
Considerando el subespacio F = L [ v1 , v 2 ,..., v k +1 ] , la intersección de F y Ker(B) no podrá ser sólo el vector nulo, sino que tendrá dimensión 1, pues la suma de dimensiones es (n–k )+(k +1)=n+1.
Factorizaciones de una matriz
pág. 77
Sea z un vector de F ∩ Ker ( B ) tal que z 2 = 1 . Por pertenecer a ambos subespacios se tendrá que: Bz = 0 Az =
(∑
r
σ u vT i =1 i i i
)
z = ∑ i =1 σ i ( vT i z ) u i k +1
Estudiando el cuadrado de la norma espectral de (A–B) y recordando que
(374)
∑ (v z) k +1
i =1
T i
2
= 1 por ser z un vector unitario,
se tiene: 2
2
2
A − B 2 ≥ ( A − B ) z 2 ≥ Az 2 − Bz 2 = ∑ i =1 σ i2 ( v Ti z ) uiT ui ≥ σ k2+1 ∑ i =1 ( v iT z ) = σ k2+1 2
2
2
k +1
k +1
(375)
En esta expresión la igualdad se obtiene cuando B=Ak . Con esto queda demostrado el teorema.
Observación: El valor singular más pequeño de A es la distancia, en norma espectral, entre A y el conjunto de matrices de rango inferior al de A.
Teorema 2: La perturbación en los valores singulares de una matriz ante una perturbación E en la matriz A está acotada mediante la expresión: σ k ( A + E ) − σ k ( A )
≤ E 2
(376)
Esto quiere decir que los valores singulares son muy estables ante perturbaciones de la matriz. Demostración: En función de los valores y vectores singulares de la matriz A se puede escribir: A = ∑ j =1σ j u j vTj , A k −1 = ∑ j =1 σ j u j vT j , k −1
r
r ( A k −1 ) = k − 1
(377)
donde r(A) es el rango de la matriz A. En virtud del teorema 1 anteriormente demostrado: min ( A + E ) − B 2 = σ k ( A + E )
(378)
r ( B )≤ k −1
Sustituyendo B=Ak –1 en la expresión (378), se obtiene: σ k ( A + E) ≤
A + E − A k −1 2 ≤ A − A k −1 2 + E 2 = σ k ( A ) + E 2
(379)
A partir de este resultado: σ k ( A + E ) − σ k ( A ) ≤
E 2
(380)
Esto es sólo parte de lo que hay que demostrar según (376). Mediante un razonamiento similar, si σ j , u j y v j son los valores y vectores singulares de la matriz perturbada A+E: A + E = ∑ j =1σ j u j vTj , ( A + E )k −1 = ∑ j =1σ j u j vT j , k −1
r
r ( A + E )k −1 = k −1
(381)
Aplicando ahora el teorema 1 y tomando B = ( A + E )k −1 se obtiene: min A − B 2 = σ k ( A )
(382)
r ( B )≤ k −1
σk (A) ≤
A + E − E − ( A + E )k −1 2 ≤ A + E − ( A + E )k −1 2 + E 2 = σ k ( A + E ) + E 2
(383)
De esta expresión se concluye lo que falta para probar (376), es decir: σ k ( A ) − σ k ( A + E ) ≤
E 2
(384)
4.4.6 Problemas de reducción de la dimensionalidad La descomposición de valores singulares permite expresar una matriz cualquiera A de rango r , como suma de r matrices de rango 1: A = U V T = ∑ i =1σ i ui viT r
(385)
Álgebra lineal numérica con Matlab
pág. 78
En función del valor de los valores singulares, la matriz A puede ser aproximada mediante las primeras k matrices de rango 1 (k
A ≈ A k = ∑ i =1σ iui vT i k
(386)
En virtud del teorema 1, el error de esta aproximación en norma cuadrática es A − A k
2
= σ k +1 .
La aproximación de la matriz A con reducción de la dimensionalidad dada por la expresión (386) tiene importantes ventajas y aplicaciones prácticas, como por ejemplo: − La matriz aproximada requiere mucha menos memoria (2k vectores de m y n elementos) y también es mucho más fácil operar con ella. − Además, la DVS proporciona información precisa y fiable sobre el error cometido en esta aproximación, que viene dado por el primer valor singular no incluido en la aproximación.
4.4.7 Aplicación de la DVS a la solución del sistema general Ax=b Se considerará en primer lugar la solución de mínimos cuadrados del sistema Ax=b (m ≥ n = r ). La solución de mínimo error cuadrático viene dada por las ecuaciones normales (ver apartado 2.10): −1
x0 = ( AT A ) AT b
(387)
Utilizando la DVS ( A=U VT ) esta expresión se puede escribir en la forma:
A = U V T , AT = V T UT
( AT A )
−1
AT = V
−2 T n
−1 n,m
UT = V
UT
⇒
AT A = V
T
VT =V 2n VT
−1 n, m
−1 n, m
−1
⇒ x0 = ( AT A ) A T b = V
U T b,
(388)
∈ R n×m (389)
Desarrollando los sumatorios implicados por esta expresión se llega a:
x0 = ∑ i =1 n
1 σi
v u b = ∑ i =1 r
T i i
uT i b σ i
vi
(390)
A continuación se estudiará la solución de mínima norma para el sistema indeterminado de rango máximo Ax=b (r = m ≤ n), que viene dada por la expresión (171) (ver apartado 2.11): −1
x∗ = AT ( AAT ) b
(391)
Introduciendo la DVS de A, tal como se ha hecho en el caso anterior:
AAT = U V T V T UT = U
T mU ,
−1
A T ( AA T ) b = V T U T U
m
U Tb = V
n, m
U Tb
(392)
Poniendo este resultado en forma de sumatorio:
x = ∑ i =1 ∗
m
1 σi
v u b = ∑ i =1 T i i
r
uT i b σ i
vi
(393)
Obsérvese las expresiones (390) y (393) son idénticas. Esto induce a aplicar la DVS al caso general de un sistema de ecuaciones lineales incompatible y sin solución de mínimo error única, es decir al caso r < min ( m, n ) . Este caso fue estudiado en el apartado 2.12 y condujo a la introducción de la matriz seudoinversa A+. La DVS permite estudiar la pseudoinversa de un modo más sencillo.
Factorizaciones de una matriz
pág. 79
Sea el sistema de ecuaciones lineales Ax=b (A∈R m×n y rango(A)
Ax = b, A = U VT
⇒ U VT x = b
(394)
Haciendo el cambio de variable y=VT x y pre-multiplicando por UT se obtiene:
y = UT b ≡ c,
= diag (σ 1 , σ 2 ,..., σ r ,0,...,0 ) ,
∈ R m×n
(395)
El sistema de ecuaciones y = c es incompatible e indeterminado, pero es también muy fácil de tratar porque la matriz del sistema sólo tiene elementos no nulos en la diagonal. Este sistema se muestra gráficamente en la Figura 42. σ i
=
=
σ i−1
y =c Figura 43. Solución de mínimo error y mínima norma y= +c.
Figura 42. Sistema de ecuaciones y=c.
Las variables libres yi , r + 1 ≤ i ≤ n, están multiplicadas por ceros y no tienen ninguna influencia en la solución. Para que la solución y sea de norma mínima estas variables deben ser nulas: yi = 0, r + 1 ≤ i ≤ n.
(396)
Por otra parte, las r primeras ecuaciones se satisfacen exactamente en la forma: yi = ci σ i , 1 ≤ i ≤ r ,
(397)
mientras que las restantes ecuaciones no se pueden satisfacer en ningún caso, pues son: 0 ⋅ yi = ci , r ≤ i ≤ m.
(398)
Definiendo la matriz pseudoinversa + como una matriz n×m con los inversos de los valores singulares distintos de cero en la diagonal, la solución de mínimo error cuadrático y mínima norma vendrá dada por (ver Figura 43):
y=
c
(399)
Esta solución minimiza también el residuo y la norma en el sistema de ecuaciones original Ax=b, pues las matrices ortogonales no cambian la norma euclídea:
r 2 = c − y 2 = UT b − V T x 2 = b − U V T x 2 = b − Ax 2 ,
y 2 = VTx 2 = x 2
(400)
Se concluye que la matriz que conduce a la solución de mínimo error y mínima en el caso general, es decir, la matriz pseudoinversa de A, es la matriz: y=
c , ( y = V T x, c = U T b )
⇒ VT x =
UT b ⇒ x = V UT b = A + b,
A + = V +UT (401)
Esta matriz seudoinversa es única, coincide con la vista anteriormente y tiene todas sus propiedades.
Álgebra lineal numérica con Matlab
pág. 80
La aplicación convencional del método de los mínimos cuadrados a través de las ecuaciones normales puede dar problemas numéricos, pues la norma espectral de la matriz AT A es el cuadrado de la norma espectral de la matriz A. Como consecuencia la condición numérica de AT A es el cuadrado de la de A, lo que puede crear dificultades numéricas. Anteriormente se ha sugerido la utilización de un sistema de ecuaciones ampliado o de la factorización QR como formas de resolver este problema. Si existen estas dificultades y además es posible que el rango de la matriz A sea menor que min(m,n) el método más estable de resolver el problema es aplicar la DVS según las expresiones (390) (393) vistas anteriormente, que se repiten aquí haciendo intervenir a la matriz seudoinversa:
Ax = b, x = A b = V U b = ∑ i =1 +
+
T
r
1 σi
v u b = ∑ i =1 T i i
r
uT i b σ i
vi
(402)
La DVS permite determinar el rango con precisión numérica, separando la información propia de la matriz de las componentes que se han podido introducir por errores de distinto tipo.