´ APUNTES DE ALGEBRA II INGENIER´IA INDUSTRIAL a t e a’ t i c a s

UNIVERSIDAD CARLOS III DE MADRID ´ DEPARTAMENTO DE MATEMATICAS

Fernando de Ter´an Vergara

´Indice general

1. PRODUCTOS ESCALARES 1.1. Definiciones y propiedades b´asicas . . . . . . . . . . . . . . 1.2. Bases ortogonales y ortonormales. M´etodo de Gram-Schmidt 1.3. Proyecciones ortogonales . . . . . . . . . . . . . . . . . . . 1.4. Factorizaci´on QR . . . . . . . . . . . . . . . . . . . . . . . 1.5. Problemas de m´ınimos cuadrados . . . . . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

3 3 9 13 16 24

2. ESTUDIO ESPECTRAL DE ENDOMORFISMOS 2.1. Matrices unitarias y ortogonales . . . . . . . . . . . . . . . . . 2.2. Matrices de Householder . . . . . . . . . . . . . . . . . . . . . 2.3. El lema de Schur y sus aplicaciones . . . . . . . . . . . . . . . 2.4. Interpretaci´on geom´etrica de las matrices ortogonales en R2 y R3

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

29 29 31 36 44

. . . . .

3. VALORES SINGULARES DE MATRICES 53 3.1. La descomposici´on en valores singulares . . . . . . . . . . . . . . . . . . . . . 53 3.2. Aplicaci´on al problema de m´ınimos cuadrados. La pseudoinversa. . . . . . . . 60 3.3. Aproximaciones de rango bajo. Compresi´on de im´agenes . . . . . . . . . . . . 64 ´ 4. FORMAS CUADRATICAS 4.1. Formas bilineales . . . . . . . . . . . . . . . . . . 4.2. Formas cuadr´aticas. Reducci´on a la forma can´onica 4.3. Formas cuadr´aticas definidas . . . . . . . . . . . . 4.4. La ley de inercia de Sylvester . . . . . . . . . . . . 4.5. Factorizaci´on de Choleski . . . . . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

67 67 71 77 79 80

´ ´ 5. CONICAS Y CUADRICAS 5.1. C´onicas en R2 . Reducci´on a la forma can´onica . . . . . . 5.2. Determinaci´on de los ejes, centro y v´ertice de una c´onica . 5.2.1. Centro y ejes de una elipse o una hip´erbola . . . . 5.2.2. Eje y v´ertice de una par´abola . . . . . . . . . . . . 5.3. Cu´adricas en R3 . Reducci´on a la forma can´onica . . . . . 5.4. Determinaci´on de los ejes, centro y v´ertice de una cu´adrica 5.4.1. Centro de un elipsoide o un hiperboloide . . . . . 5.4.2. V´ertice de un paraboloide . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

83 83 86 87 88 90 96 96 97

III

. . . . .

. . . . .

. . . . .

IV

´ INDICE GENERAL 5.5. Representaci´on gr´afica de las cu´adricas en forma reducida . . . . . . . . . . . 98

6. NORMAS DE VECTORES Y MATRICES. SENSIBILIDAD DE MATRICES 6.1. Normas vectoriales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2. Normas matriciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3. El radio espectral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4. Errores en inversas y en soluciones de sistemas lineales. N´umero de condici´on

. . . .

103 103 104 108 109

´ A LOS METODOS ´ 7. INTRODUCCION ITERATIVOS 117 7.1. M´etodos iterativos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 7.2. Aplicaci´on a los sistemas lineales. M´etodos de Jacobi y Gauss-Seidel . . . . . 119

´ NOTACION En estos apuntes y, salvo que se indique expl´ıcitamente lo contrario, utilizaremos la siguiente notaci´on:

S´ımbolo C Cm×n R Rm×n K dim(V) In 0p×q det(A) rango(A) AT A∗ diag(A1 , . . . , An ) A(i, j) ker (A) ¤ £ col (A) A = a1 . . . ak Gen (S) CRn M (B1 , B2 ) c¯ |c| Re (v) Im (v) [v]B ∀

Significado Cuerpo de los n´umeros complejos Espacio vectorial de las matrices m × n con coeficientes complejos Cuerpo de los n´umeros reales Espacio vectorial de las matrices m × n con coeficientes reales Cuerpo de los n´umeros reales o complejos (indistintamente) Dimensi´on del espacio vectorial V Matriz identidad de tama˜no n × n Matriz de tama˜no p × q con todas sus entradas nulas Determinante de la matriz A Rango de la matriz A Traspuesta de la matriz A Traspuesta conjugada de la matriz A Matriz diagonal por bloques, que son A1 , . . . , An Entrada (i, j) de la matriz A N´ucleo de la matriz A Espacio generado por las columnas de la matriz A Matriz cuyas columnas son a1 , . . . , ak Espacio vectorial (sobre K) generado por el conjunto S Base can´onica de Rn Matriz de cambio de base de B1 a B2 Conjugado del n´umero complejo c M´odulo del n´umero complejo c Parte real del vector (complejo) v Parte imaginaria del vector (complejo) v Vector v expresado en coordenadas en la base B Para todo

Adem´as, seguiremos el convenio de considerar vectores columna, es decir, las coordenadas del vector v se suponen ordenadas en una columna. Para expresar vectores fila, por tanto, utilizaremos la traspuesta, v T . Igualmente, denotaremos por ei al i-´esimo vector can´onico de Rn , es decir ei =

£

i)

¤T 0 ... 1 ...0 .

1 PRODUCTOS ESCALARES

´ 1.1 Definiciones y propiedades basicas A lo largo de estos apuntes, V denotar´a un espacio vectorial sobre el cuerpo, C, de los n´umeros complejos o sobre el cuerpo, R, de los n´umeros reales. En algunos casos el cuerpo base ser´a u´ nicamente C o u´ nicamente R, pero en otros ser´a C o´ R, indistintamente. En este u´ ltimo caso denotaremos al cuerpo gen´ericamente por K. La definici´on fundamental de este cap´ıtulo es la siguiente. Definici´on 1.1.1 (Producto escalar herm´ıtico) Sea V un espacio vectorial sobre C. Una aplicaci´on Φ : V × V −→ C es un producto escalar herm´ıtico en V si cumple los siguientes axiomas: i) Φ(u, u) es un n´umero real no negativo, para todo u ∈ V, y Φ(u, u) = 0 si y s´olo si u=0 (Positiva) ii) Φ(u, v) = Φ(v, u), para cualesquiera u, v ∈ V

(Herm´ıtica)

iii) Φ(αu + βv, w) = αΦ(u, w) + βΦ(v, w), para cualesquiera n´umeros complejos α, β y vectores u, v, w ∈ V (Bilineal) En el caso de espacios vectoriales reales la definici´on de producto escalar se obtiene particularizando la definici´on anterior a un espacio sobre R. Definici´on 1.1.2 (Producto escalar) Sea V un espacio vectorial sobre R. Un aplicaci´on Φ : V × V −→ R es un producto escalar en V si cumple los siguientes axiomas: i) Φ(u, u) es un n´umero real no negativo, para todo u ∈ V, y Φ(u, u) = 0 si y s´olo si u=0 (Positiva) 3

1. PRODUCTOS ESCALARES

4 ii) Φ(u, v) = Φ(v, u), para cualesquiera u, v ∈ V

(Sim´etrica)

iii) Φ(αu + βv, w) = αΦ(u, w) + βΦ(v, w), para cualesquiera n´umeros reales α, β y vectores u, v, w ∈ V (Bilineal) Diremos, por tanto, que un producto escalar herm´ıtico es una “forma bilineal herm´ıtica definida positiva”, mientras que un producto escalar (real) es una “forma bilineal sim´etrica definida positiva”. A partir de este momento, y salvo que sea necesario, nos referiremos simplemente a “productos escalares”, incluyendo con esta denominaci´on tanto a los herm´ıticos como a los reales. Notaci´on 1.1.3 Si Φ es un producto escalar, entonces el producto escalar de u con v es el n´umero Φ(u, v). Habitualmente lo escribiremos < u, v >. Otra notaci´on que puede encontrarse en la literatura es, simplemente, u · v. Obs´ervese que la propiedad iii) del producto escalar s´olo establece la linealidad en la primera componente. En el siguiente ejercicio se pide demostrar que tambi´en es lineal en la segunda y que, por tanto, el calificativo de “bilineal” est´a completamente justificado. Ejercicio 1.1.4 Sea Φ : V × V → K un producto escalar sobre V. Demostrar, usando u´ nicamente la definici´on, que Φ es lineal en la segunda componente, es decir: Φ(u, αv + βw) = αΦ(u, v) + βΦ(u, w) , para cualesquiera α, β n´umeros en K y vectores u, v, w ∈ V. Hay que hacer notar que, en el caso herm´ıtico, la linealidad en la primera componente se da con el conjugado de los coeficientes, mientras que en la segunda componente se da con los mismos coeficientes (sin conjugar). Ejercicio 1.1.5 Demostrar, usando los axiomas i), ii) y iii) de la definici´on, que cualquier producto escalar Φ (real o herm´ıtico) sobre V satsiface Φ(0, u) = Φ(u, 0) = 0 , para todo u ∈ V. Observaci´on 1.1.6 Es habitual encontrar en la literatura la denominaci´on “producto interior” para referirse al producto escalar. A continuaci´on, vamos a ver unos cuantos ejemplos especialmente relevantes de productos escalares. Ejemplo 1.1.7

1. En el espacio vectorial V = Kn definimos el producto escalar < u, v >= u∗ v.

1.1 Definiciones y propiedades b´asicas

5

Es f´acil obtener la expresi´on de dicho producto escalar en la base can´onica de Kn . Si expresamos los vectores u y v en dicha base     u1 v1  u2   v2      u =  .  , v =  . ,  ..   ..  un

vn

entonces < u, v >= u∗ v = u ¯1 v1 + u ¯2 v2 + · · · + u ¯n vn . Si K = R entonces el producto anterior es < u, v >= uT v = u1 v1 + u2 v2 + · · · + un vn . Este producto escalar se conoce como el producto escalar usual. 2. Sea V el espacio vectorial de las sucesiones {xn }n∈N , con xn ∈ K. Definimos el producto escalar ∞ X < {xn } , {yn } >= x ¯ n yn . n=1

3. En el espacio de funciones reales continuas en un intervalo [a, b] definimos el producto escalar Z b

< f, g >=

f (x)g(x)dx. a

4. En el espacio Kn×n de matrices cuadradas n × n con entradas en K definimos < A, B >= traza(A∗ B). Ejercicio 1.1.8 Demostrar que las aplicaciones del ejemplo anterior son, efectivamente, productos escalares. (Indicaci´on: Hay que comprobar que cumplen los axiomas i), ii) y iii) de la definici´on de producto escalar). Un espacio vectorial real V dotado de un producto escalar se suele llamar espacio vectorial eucl´ıdeo, mientras que un espacio vectorial complejo con un producto vectorial herm´ıtico se denomina espacio vectorial unitario. En los espacios vectoriales eucl´ıdeos/unitarios se pueden definir los conceptos de longitud, distancia y a´ ngulo. Es importante destacar que para definir estos conceptos utilizaremos u´ nicamente el producto escalar. En lo sucesivo, denotaremos por (V, < ·, · >) a un espacio vectorial eucl´ıdeo/unitario, donde V es el espacio vectorial y < ·, · > es el producto escalar. Definici´on 1.1.9 Sea (V, < ·, · >) un espacio vetorial eucl´ıdeo/unitario. La norma de un vector u ∈ V es el n´umero real k u k=< u, u >1/2 . Diremos que u ∈ V es un vector unitario si k u k= 1.

1. PRODUCTOS ESCALARES

6

N´otese que la norma de un vector depende del producto escalar asociado. Por tanto, un vector puede ser unitario con una norma y no serlo con otra. Propiedades de la norma: En las siguientes propiedades (V, < ·, · >) es un espacio vectorial eucl´ıdeo/unitario sobre el cuerpo K. 1. k u k> 0, ∀u ∈ V − {0V } , y k 0V k= 0 2. k αu k= |α| k u k, ∀α ∈ K, u ∈ V 3. k u + v k≤k u k + k v k, ∀u, v ∈ V

(Desigualdad triangular)

Ejercicio 1.1.10 Demostrar las propiedades anteriores. Ejemplo 1.1.11 En V = Kn la norma asociada al producto escalar usual est´a dada por k u k2 = (|u1 |2 + |u2 |2 + · · · + |un |2 )1/2 donde u = [u1 . . . un ]T . Esta norma se denomina la norma-2 (o norma usual). Definici´on 1.1.12 Sea (V, < ·, · >) un espacio vectorial eucl´ıdeo/unitario y sean u, v ∈ V dos vectores de V. El a´ ngulo entre u y v es el u´ nico n´umero real ϕ ∈ [0, π) tal que cosϕ =

< u, v > . k u kk v k

La definici´on anterior coincide con el concepto habitual de a´ ngulo en R2 y R3 cuando consideramos el producto escalar usual, como se muestra en el siguiente ejemplo. Ejemplo 1.1.13 Consideremos el plano R2 con el producto escalar usual y sean u, v dos vectores del plano. Sean ϕu y ϕv los a´ ngulos (en el sentido habitual) que forman, respectivamente, u y v con el eje horizontal, y supongamos, sin p´erdida de generalidad, que ϕu ≤ ϕv . y v 

*

u

ϕ ϕv

ϕu

x Llamemos ϕ al a´ ngulo (tambi´en en el sentido habitual) que forman los vectores u, v. Entonces ϕ = ϕv − ϕu y, usando la identidad trigonom´etrica bien conocida del coseno de una diferencia, tenemos cos ϕ = cos(ϕv − ϕu ) = cos ϕv cos ϕu + sen ϕv sen ϕu =

v2 u2 v1 u1 + kv kkuk kv kkuk

1.1 Definiciones y propiedades b´asicas =

7

u1 v1 + u2 v2 < u, v > = , k u kk v k k u kk v k

por tanto ϕ es el a´ ngulo entre u y v, tal y como se ha introducido en la Definici´on 1.1.12. Otro concepto importante es el de perpendicularidad, que definimos a continuaci´on. Definici´on 1.1.14 Sea (V, < ·, · >) un espacio vectorial eucl´ıdeo/unitario. 1) Dos vectores u, v ∈ V son ortogonales (o perpendiculares) si < u, v >= 0. Lo denotaremos por u ⊥ v. 2) Un conjunto de vectores no nulos {u1 , . . . , uk } se dice ortogonal si ui ⊥ uj para cualesquiera i, j con i 6= j. 3) El subconjunto ortogonal de un subconjunto S de V es el conjunto S ⊥ = {v ∈ V : < u, v >= 0 ∀u ∈ S} . Conviene se˜nalar que, en virtud de la simetr´ıa/hermiticidad del producto escalar/herm´ıtico, < u, v >= 0 si y s´olo si < v, u >= 0 y, por tanto, podemos decir que u, v son perpendiculares sin tener en cuenta el orden en que se indican ambos vectores. Ejercicio 1.1.15 Demostrar que si S es un subconjunto cualquiera de un espacio vectorial eucl´ıdeo/unitario entonces S ∩ S ⊥ = {0}. Proposici´on 1.1.16 Sea (V, < ·, · >) un espacio vectorial eucl´ıdeo/unitario sobre K y S un subconjunto cualquiera de V. Entonces el subconjunto ortogonal de S, S ⊥ , es un subespacio vectorial de V. Demostraci´on. Dados dos vectores v1 , v2 ∈ S ⊥ , cualquier combinaci´on lineal de ellos, αv1 + βv2 , con α, β ∈ K, satisface < u, αv1 + βv2 >= α < u, v1 > +β < u, v2 >= α0 + β0 = 0, para todo u ∈ S. Por tanto αv1 +βv2 ∈ S ⊥ , lo que significa que S ⊥ es un subespacio vectorial. N´otese que S ⊥ siempre es un subespacio vectorial, aunque S no lo sea. En el caso de que el conjunto S sea un subespacio vectorial, la dimensi´on del subespacio ortogonal es complementaria a la de S. Lema 1.1.17 Sea (V, < ·, · >) un espacio vectorial eucl´ıdeo/unitario de dimensi´on n y W un subespacio vectorial de V. Entonces dim W + dim W ⊥ = n. En el siguiente ejercicio se pide demostrar que, en general, para comprobar si un vector est´a en un subespacio ortogonal no es necesario comprobar la ortogonalidad con todos los vectores del conjunto (que pueden ser infinitos). Concretamente, si partimos de un subespacio vectorial W , para saber si un vector cualquiera est´a en el subespacio ortogonal W ⊥ es suficiente con que dicho vector sea ortogonal a un conjunto de generadores de W .

1. PRODUCTOS ESCALARES

8

Ejercicio 1.1.18 Demostrar que si W es un subespacio vectorial entonces v ∈ W ⊥ si y s´olo si u es ortogonal a un conjunto de generadores de W . Por otra parte, los subconjuntos ortogonales tienen la siguiente propiedad Lema 1.1.19 Sea (V, < ·, · >) un espacio vectorial eucl´ıdeo/unitario sobre K y {u1 , . . . , uk } un subconjunto ortogonal. Entonces {u1 , . . . , uk } es linealmente independiente. Demostraci´on. Procederemos por reducci´on al absurdo. Supongamos, por tanto, que hay una relaci´on de dependencia entre los vectores del conjunto ortogonal: α1 u1 + α2 u2 + · · · + αk uk = 0, con αi 6= 0, para alg´un 1 ≤ i ≤ k. Multiplicando por ui la anterior expresi´on obtenemos < α1 u1 + α2 u2 + . . . + αk uk , ui >= 0 y, usando la linealidad en la primera componente, α1 < u1 , ui > +α2 < u2 , ui > + . . . + αi < ui , ui > + . . . + αk < uk , ui >= 0. Como {u1 , . . . , uk } es ortogonal, tenemos que < uj , ui >= 0 para j 6= i, luego la anterior identidad es equivalente a αi < ui , ui >= 0. Finalmente, como ui 6= 0 (los vectores de un conjunto ortogonal son, por definici´on, no nulos), se tiene que < ui , ui >6= 0 (por el axioma i) del producto escalar), luego tiene que ser αi = 0, en contradicci´on con la hip´otesis. La u´ ltima definici´on de esta secci´on es el concepto de distancia entre dos vectores. Definici´on 1.1.20 Sea (V, < ·, · >) un espacio vectorial eucl´ıdeo/unitario y sean u, v dos vectores de V. La distancia entre u y v es el n´umero real no negativo d(u, v) =k u − v k . N´otese que d(u, v) = d(v, u) y que d(u, v) = 0 ⇔ u = v. ´ matricial en una base Expresion Sea (V, < ·, · >) un espacio vectorial eucl´ıdeo/unitario de dimensi´on n y B = {u1 , . . . , un } una base de V. Vamos a encontrar una expresi´on matricial para el producto escalar de dos vectores cualesquiera expresados en esa base. Sean, por tanto, v, w dos vectores de V cuyas coordenadas en la base B son     v1 w1     v =  ...  , w =  ...  , vn wn es decir v = v1 u1 + . . . + vn un , w = w1 u1 + . . . + wn un .

1.2 Bases ortogonales y ortonormales. M´etodo de Gram-Schmidt

9

(N´otese que vi , wi , para i = 1, . . . , n son n´umeros, mientras que ui , para i = 1, . . . , n, son vectores). Usando la bilinealidad   n n n n n n X X X X X X < v, w >=< vi ui , wj uj >= v¯i < ui , wj uj >= v¯i  wj < ui , uj > i=1

j=1

i=1

=

j=1

n X

i=1

j=1

v¯i wj < ui , uj >,

i,j=1

que, matricialmente, se puede expresar as´ı  < u1 , u1 > < u1 , u2 > . . . < u1 , un > £ ¤  < u2 , u1 > < u2 , u2 > . . . < u2 , un > < v, w >= v¯1 . . . v¯n  .. .. ..  . . . < un , u1 > < un , u2 > . . . < un , un >

    

w1 w2 .. .

    

wn

= v ∗ GB w, donde

   GB =  

< u1 , u1 > < u1 , u2 > . . . < u1 , un > < u2 , u1 > < u2 , u2 > . . . < u2 , un > .. .. .. . . . < un , u1 > < un , u2 > . . . < un , un >

    

se conoce como la matriz de Gram del producto escalar < ·, · > en la base B.

´ 1.2 Bases ortogonales y ortonormales. Metodo de Gram-Schmidt A largo de este apartado (V, < ·, · >) es un espacio vectorial eucl´ıdeo/unitario. Una base ortogonal de V es una base de V formada por vectores ortogonales. Si, adem´as, los vectores son unitarios, diremos que forman una base ortonormal. Lo que vamos a ver a continuaci´on es que todo espacio vectorial eucl´ıdeo/unitario (o cualquier subespacio vectorial de ese espacio) tiene, al menos, una base ortonormal. M´as a´un, veremos un procedimiento para obtener una base ortogonal/ortonormal a partir de una base cualquiera del espacio (o subespacio). El procedimiento que estudiaremos es el que se conoce con el nombre de ortogonalizaci´on de Gram-Schmidt. Algoritmo 1: Ortogonalizaci´on de Gram-Schmidt Entrada: {v1 , v2 , . . . , vn } conjunto de n vectores. Salida: {u1 , u2 , . . . , un } conjunto de n vectores ortogonales dos a dos. u 1 = v1 u2 = v2 − .. . ³ uk = vk −

u1

+ ... +

uk−1

´

1. PRODUCTOS ESCALARES

10

Las propiedades fundamentales del conjunto de vectores que produce el m´etodo anterior son las descritas en el siguiente teorema. Teorema 1.2.1 El Algoritmo 1 produce un conjunto de vectores ortogonales dos a dos y que satisfacen Gen(v1 , . . . , vk ) = Gen(u1 , . . . , uk ), k = 1, 2, . . . , n.

(1.1)

M´as a´un, el n´umero de vectores linealmente independientes en el conjunto de partida {v1 , v2 , . . . , vn } coincide con el n´umero de vectores no nulos del conjunto de salida {u1 , u2 , . . . , un } (que, adem´as, son linealmente independientes). En particular, si el conjunto de partida es linealmente independiente, entonces el Algoritmo 1 produce un conjunto ortogonal de n vectores (no nulos). Demostraci´on. Vamos a demostrar que los vectores que produce el Algoritmo 1 son ortogonales dos a dos y que satisfacen la propiedad (1.1). Pospondremos para la Secci´on 1.3 la demostraci´on acerca de la relaci´on entre la independencia lineal del conjunto de partida y el n´umero de vectores no nulos del conjunto de salida. Para demostrar (1.1) n´otese, en primer lugar, que Gen(v1 , . . . , vk ) ⊆ Gen(u1 , . . . , uk ), para k = 1, . . . , n, pues vk ∈ Gen(u1 , . . . , uk ) (basta despejar vk en la expresi´on correspondiente a uk en el Algoritmo 1). La inclusi´on opuesta, Gen(v1 , . . . , vk ) ⊇ Gen(u1 , . . . , uk ), se puede obtener f´acilmente por inducci´on: para k = 1 es obvia, pues u1 = v1 . Si la suponemos cierta para k − 1, tenemos que uk ∈ Gen{vk , u1 , . . . , uk−1 } ⊆ Gen{vk , v1 , . . . , vk−1 } (donde la hip´otesis de inducci´on se ha usado en esta u´ ltima inclusi´on). Veamos ahora que los vectores del conjunto de salida son ortogonales dos a dos. Lo que vamos a ver, por inducci´on sobre k, es que el vector uk producido en el paso k-´esimo es ortogonal a todos los vectores anteriores, u1 , . . . , uk−1 . El caso inicial es para k = 2. En este caso tenemos que < u2 , u1 >=< v2 −

< u 1 , v2 > < v2 , u1 > u1 , u1 >=< v2 , u1 > − < u1 , u1 >= 0. < u 1 , u1 > < u1 , u1 >

(Para obtener la u´ ltima igualdad hemos usado los axiomas del producto escalar). Supongamos ahora que ui ⊥ uj , para i, j = 1, . . . , k − 1 con i 6= j (hip´otesis de inducci´on). Entonces, para j ≤ k − 1 tenemos que < uk , uj >=< vk −

k−1 k−1 X X < vk , ui > < ui , vk > ui , uj >=< vk , uj > − < ui , uj > < ui , ui > < ui , ui > i=1

i=1

=< vk , uj > − < vk , uj >= 0, donde, en la u´ ltima igualdad hemos usado la hip´otesis de inducci´on < ui , uj >= 0 para i = 1, . . . , k − 1 , i 6= j. Esto demuestra que el conjunto {u1 , . . . , un } es ortogonal. Observaci´on 1.2.2 Si eliminamos los vectores nulos producidos en el proceso de Gram-Schmidt, obtendremos una base ortogonal del espacio generado por los vectores de partida {v1 , v2 , . . . , vn }.

1.2 Bases ortogonales y ortonormales. M´etodo de Gram-Schmidt

11

Para obtener una base ortonormal tenemos el siguiente algoritmo, donde en cada paso se normalizan los vectores ui que se obtienen en el Algoritmo 1 con el fin de obtener unos nuevos vectores unitarios, que llamaremos qi . Algoritmo 2: Ortonormalizaci´on de Gram-Schmidt Entrada: {v1 , v2 , . . . , vn } conjunto de n vectores. Salida: {q1 , q2 , . . . , qn } conjunto de n vectores unitarios (o nulos) ortogonales dos a dos. u1 = v1 , q1 =

u1 ku1 k

u2 = v2 − < q1 , v2 > q1 , q2 = kuu22 k .. . uk = vk − (< q1 , vk > q1 + . . . + < qk−1 , vk > qk−1 ) , qk =

uk kuk k

Por supuesto, tenemos el resultado an´alogo al Teorema 1.2.1, cuya demostraci´on es exactamente la misma, pues el hecho de que los vectores qi son unitarios es evidente. Teorema 1.2.3 El Algoritmo 2 produce un conjunto de vectores unitarios que son ortogonales dos a dos y que satisfacen Gen(v1 , . . . , vk ) = Gen(q1 , . . . , qk ), k = 1, 2, . . . , n.

(1.2)

M´as a´un, el n´umero de vectores linealmente independientes en el conjunto de partida {v1 , v2 , . . . , vn } coincide con el n´umero de vectores no nulos del conjunto de salida {q1 , q2 , . . . , qn } (que, adem´as, son linealmente independientes). En particular, si el conjunto de partida es linealmente independiente, entonces el Algoritmo 2 produce un conjunto ortonormal de n vectores (no nulos). Es importante observar que a la hora de normalizar suelen aparecer ra´ıces cuadradas que no estaban presentes en los vectores de partida. Por tanto, si se quiere calcular una base ortonormal y se van a hacer los c´alculos a mano, lo m´as pr´actico para evitar c´alculos engorrosos con ra´ıces cuadradas es obtener una base ortogonal con el Algoritmo 1 y despu´es normalizarla (es decir, normalizar cada vector). La u´ nica diferencia entre el Algoritmo 2 y este procedimiento es que en el primer caso los vectores se normalizan en cada paso, mientras que en el segundo se normalizan todos al final. Esto significa, en particular, que los vectores ui que se obtienen en ambos casos son exactamente los mismos. Ejercicio 1.2.4 Demostrar que los vectores uj que se obtienen al aplicar el Algoritmo 1 y los que se obtienen al aplicar el Algoritmo 2 son exactamente los mismos si en ambos casos se parte del mismo conjunto de vectores {v1 , . . . , vn }. Problema Resuelto 1.2.5 Calcular una base ortonormal, con el producto escalar usual, del subespacio de R4 generado por los vectores       2 −1 1  4   3   2       v1 =   0  , v2 =  1  , v3 =  6  . 2 1 1

1. PRODUCTOS ESCALARES

12

Soluci´on. Obtendremos la base ortonormal de las dos maneras descritas anteriormente: la primera aplicando el Algoritmo 1 y normalizando los vectores al final y la segunda aplicando el Algoritmo 2. Primera forma: Primero aplicamos el Algoritmo 1 al conjunto {v1 , v2 , v3 }. 

 1  2   u1 = v1 =   0  1

     −2 −1 1  3  6 2   1        u = u2 = v2 − 1 0   1  2 1      4  − 12  2  − u3 = v3 − 6    1 1 2 2 6 0  2 1 



  −2 2    6 1   −1 6 1 = 5 0 0

  . 

Ahora normalizamos los vectores u1 , u2 y u3 para obtener la base ortonormal        1 −2 2      1   1   1  −1  2 1       B = √  , √  . ,√  0 1  5   6 30    6  1 0 0 Segunda forma: Aplicamos el Algoritmo 2 al conjunto {v1 , v2 , v3 }.     1 1  2     , q1 = u1 = √1  2  u1 = v1 =  ku1 k  0  6 0  1 1         1 −2 −1 −2      1   3   u2   √6  √1  2   1  √1  u2 = v2 − < q1 , v2 >=   1  − 6  6  0  =  1  , q2 = ku2 k = 6  1  0 1  0    1 1 2  4  12  1  2   √  √   u3 = v3 − < q1 , v3 > q1 − < q2 , v3 > q2 =   6  − 6  6  0  1 2        2 −2 2         −1 1 −1 u  , q3 = 3 = √1   =   . √1  − √66  ku3 k  6  1   5  30  5  0 0 0 La base que obtenemos, B = {q1 , q2 , q3 } , es la misma que en procedimiento anterior. M´as a´un, como ya hemos comentado, los vectores ui que se obtienen en cada paso son los mismos en ambos procedimientos.

1.3 Proyecciones ortogonales

13

1.3 Proyecciones ortogonales El concepto de proyecci´on ortogonal ya es conocido. Por ejemplo, en el Bachillerato se estudian las proyecciones de un punto o un vector sobre una recta o un plano. La generalizaci´on al concepto de proyecci´on ortogonal de un vector sobre un subespacio vectorial cualquiera es el objeto del siguiente Teorema. Teorema 1.3.1 (de la proyecci´on ortogonal) Sean (V, < ·, · >) un espacio vectorial eucl´ıdeo/ unitario, W un subespacio vectorial de V y v un vector de V. Entonces existe un u´ nico vector, llamado proyecci´on de v sobre W , que denotaremos por proyW (v), tal que i) proyW (v) ∈ W . ii) v − proyW (v) ∈ W ⊥ . El Teorema 1.3.1 permite descomponer el vector v, de manera u´ nica, en la forma v = proyW (v) + v ⊥ ,

(1.3)

donde proyW (v) es un vector de W y v ⊥ es un vector perpendicular a W . 6

v⊥

v -

W

proyW (v)

La igualdad (1.3) se conoce como la descomposici´on ortogonal de v en W . N´otese que, en el caso particular en que v ∈ W , se tiene que proyW (v) = v y v ⊥ = 0, mientras que si v ∈ W ⊥ entonces proyW (v) = 0 y v ⊥ = v. Observaci´on 1.3.2 No hay que confundir v ⊥ con {v}⊥ o´ Gen(v)⊥ . En el primer caso se trata del vector v − proyW (v), que depende del subespacio W , mientras que en el segundo caso se trata del subespacio ortogonal al (subespacio generado por el) vector v. Demostraci´on. (del Teorema 1.3.1) La siguiente demostraci´on es constructiva, es decir, daremos una f´ormula para construir el vector proyW (v). El procedimiento es el siguiente. Sea B = {u1 , . . . , uk } una base ortogonal de W . Entonces, el vector proyW (v) =

< u1 , v > < u2 , v > < uk , v > u1 + u2 + . . . + uk < u 1 , u1 > < u2 , u2 > < u k , uk >

(1.4)

satisface las propiedades i) y ii) del enunciado. Ve´amoslo. La condici´on i) es inmediata, porque el vector definido en (1.4) es una combinaci´on lineal de vectores de W .

1. PRODUCTOS ESCALARES

14

Para demostrar la propiedad ii) ser´a suficiente, en virtud del Ejercicio 1.1.18, demostrar que el vector v − proyW (v) es ortogonal a los vectores de la base B. Fij´emonos en un ´ındice cualquiera j ∈ {1, . . . , k}. Para este ´ındice se tiene k k X X < ui , v > < v, ui > < v−proyW (v), uj >=< v− ui , uj >=< v, uj > − < ui , uj > . < u i , ui > < u i , ui > i=1

i=1

Como B es un conjunto ortogonal, lo anterior es igual a < v, uj > −

< v, uj > < uj , uj >=< v, uj > − < v, uj >= 0, < uj , uj >

lo que significa que v − proyW (v) es ortogonal a cualquier vector de B. ´ Finalmente, vamos a demostrar que la descomposici´on (1.3) es unica. Por supuesto, la f´ormula (1.4) est´a u´ nivocamente definida, pero podr´ıa existir otra “proyecci´on” dada por otra f´ormula o, incluso, bases distintas podr´ıan dar lugar a “proyecciones” distintas. Supongamos que hay dos descomposiciones distintas v = v1 + v1⊥ , v = v2 + v2⊥ , donde v1 , v2 son vectores de W y v1⊥ , v2⊥ son vectores de W ⊥ . Entonces, de la igualdad v1 + v1⊥ = v2 + v2⊥ se obtiene que v1 − v2 = v2⊥ − v1⊥ . El vector v1 − v2 es un vector de W , mientras que v2⊥ − v1⊥ lo es de W ⊥ . Como ambos son el mismo, se trata de un vector en W ∩ W ⊥ . En virtud del Ejercicio 1.1.15, dicho vector ha de ser el vector nulo. Por tanto v1 = v2 y v1⊥ = v2⊥ , lo que significa que la descomposici´on es u´ nica. En el siguiente cuadro conclusivo se resume el procedimiento que acabamos de ver para calcular la proyecci´on ortogonal. C´ alculo de proyW (v) 1. Calcular una base ortogonal de W : B = {u1 , . . . , uk } 2. Aplicar la f´ormula: proyW (v) =

< u1 , v > < uk , v > u1 + . . . + uk < u 1 , u1 > < uk , uk >

Es importante destacar que la base B ha de ser ortogonal, pues la f´ormula anterior no es v´alida si B no es ortogonal (v´ease la demostraci´on del Teorema 1.3.1, ¿en qu´e momento se usa que la base es ortogonal?). Ejercicio 1.3.3 Comprobar, con un ejemplo, que la f´ormula (1.4) no es v´alida si la base no es ortogonal (Indicaci´on: Por ejemplo, en R3 consid´erese el plano horizontal z = 0 y una base de dicho plano que no sea ortogonal. Apl´ıquese la f´ormula (1.4) para esa base y un vector v cualquiera de R3 y compru´ebese que el vector que resulta no es la proyecci´on ortogonal de v sobre el plano).

1.3 Proyecciones ortogonales

15

Las proyecciones ortogonales permiten obtener una interpretaci´on geom´etrica del proceso de ortogonalizaci´on de Gram-Schmidt, lo que ayudar´a a recordarlo m´as f´acilmente. Para obtener esta interpretaci´on, denotaremos por Vj−1 = Gen(v1 , . . . , vj−1 ) al subespacio vectorial generado por los primeros j − 1 vectores del conjunto de partida (v´ease el Algoritmo 1). Los vectores uj del nuevo conjunto que nos porporciona el Algoritmo 1 est´an dados por uj = vj − proyVj−1 (vj ). Es decir, el vector uj que proporciona el m´etodo de Gram-Schmidt en el paso j-´esimo es el vector que anteriormente hemos denotado por vj⊥ , y que, como hemos visto, cumple ⊥ vj⊥ ∈ Vj−1 = Gen(u1 , . . . , uj−1 )⊥ ,

es decir, uj es ortogonal a u1 , . . . , uj−1 . Problema Resuelto 1.3.4 Calcular la proyecci´on ortogonal del vector v = sobre el subespacio vectorial de R4 del Problema Resuelto 1.2.5.

£

1 1 1 −3

¤T

Soluci´on. Llamemos W al subespacio vectorial del Problema Resuelto 1.2.5. Como ya hemos calculado en dicho problema una base ortogonal de W , podemos aplicar directamente la f´ormula (1.4)     1 −2   < u2 , v > < u3 , v > 0 2  < u1 , v > + 0 1  u1 + u2 + u3 =  proyW (v) = < u1 , u1 > < u 2 , u2 > < u3 , u3 > 6 0  6 1  1 0     2 2/5   −1/5  6  −1 = . +  30  5   1  0 0 Para terminar este apartado, vamos a demostrar la parte que ha quedado pendiente del enunciado de los Teoremas 1.2.1 y 1.2.3, es decir, que el n´umero de vectores linealmente independientes del conjunto {v1 , . . . , vn } coincide con el n´umero de vectores no nulos del conjunto {u1 , . . . , un } que produce el Algoritmo 1 y del conjunto {q1 , . . . , qn } que produce el Algoritmo 2. Supongamos que, para alg´un k ∈ {1, . . . , n}, tenemos que vk ∈ Gen(v1 , . . . , vk−1 ) = Gen(u1 , . . . , uk−1 ) = Gen(q1 , . . . , qk−1 ) Entonces (siguiendo la notaci´on anterior) uk = vk − proyVk−1 (vk ) = vk⊥ , y, como vk ∈ Gen(v1 , . . . , vk−1 ), se tiene que uk = 0 y, por tanto, qk = 0. A la inversa, si qk = 0, para alg´un sub´ındice k ∈ {1, . . . , n}, entonces tambi´en uk = 0 y eso significa que vk⊥ = 0, es decir, vk ∈ Gen(v1 , . . . , vk−1 ). As´ı pues uk = 0 ⇔ qk = 0 ⇔ vk ∈ Gen(v1 , . . . , vk−1 )

1. PRODUCTOS ESCALARES

16 ´ matricial de la proyeccion ´ ortogonal Expresion

£ ¤ En el espacio vectorial Kn con el producto escalar usual, si Q = q1 . . . qk es una matriz cuyas columnas son una base ortonormal de W , entonces la f´ormula (1.4) es equivalente a la expresi´on matricial proyW (v) = QQ∗ v Por este motivo, la matriz QQ∗ se conoce como matriz de proyecci´on (sobre el espacio generado por las columnas de Q).

´ QR 1.4 Factorizacion En este apartado vamos a aplicar el m´etodo de ortonormalizaci´on de Gram-Schmidt al espacio de columnas de una matriz dada, A. En primer lugar, trataremos el caso en que todas las columnas de A son linealmente independientes, lo que nos permite tomar dichas columnas como base del espacio col(A). Despu´es estudiaremos el caso en que las columnas no son linealmente independientes. Antes de comenzar, fijaremos la siguiente definici´on. Definici´on 1.4.1 Diremos que una matriz U ∈ Cn×m (respectivamente, P ∈ Rn×m ) tiene columnas ortonormales si U ∗ U = Im (resp. P T P = Im ). N´otese que la definici´on anterior equivale a decir que las columnas de U (resp. P ) forman un conjunto ortonormal con respecto al producto escalar usual. Es importante destacar que una matriz de tama˜no n × m con m > n no puede tener columnas ortonormales, pues eso supondr´ıa que existe un conjunto ortogonal (y, por lo tanto, linealmente independiente) de m vectores en Cn (resp. Rn ), lo cual es imposible si m > n. Por tanto, para que la Definici´on 1.4.1 tenga sentido es necesario que n ≥ m. Por otra parte, tambi´en es importante tener en cuenta que la igualdad U ∗ U = Im no implica necesariamente la igualdad U U ∗ = In . De hecho, como hemos visto en el p´arrafo anterior, ´ esto u´ ltimo no puede ocurrir si n > m. Unicamente en el caso cuadrado, n = m, la igualdad ∗ ∗ U U = In s´ı que implica U U = In y, adem´as, en este caso, U ∗ = U −1 . Este caso ser´a tratado en el Tema 2. Recordamos la siguiente definici´on, que ser´a utilizada en lo sucesivo. Definici´on 1.4.2 Una matriz T ∈ Cn×m se dice triangular superior si T (i, j) = 0 para todo i > j, y se dice triangular inferior si T (i, j) = 0 para todo i < j. Una matriz D ∈ Cn×m es diagonal si D(i, j) = 0 para todo i 6= j. Obs´ervese que una matriz diagonal no es necesariamente cuadrada y que, por otro lado, una matriz es diagonal si y s´olo si es triangular inferior y superior al mismo tiempo. Ahora ya podemos enunciar el resultado fundamental de este apartado. Teorema 1.4.3 (Factorizaci´on QR) Dada una matriz A ∈ Cn×m , con n ≥ m, existen dos matrices Q y R (que dependen de A) de manera que A = QR y

1.4 Factorizaci´on QR

17

i) Q ∈ Cn×m tiene columnas ortonormales (Q∗ Q = Im ). ii) R ∈ Cm×m es triangular superior con entradas diagonales no negativas. Si, adem´as, las columnas de A son linealmente independientes entonces R puede tomarse con todos sus elementos diagonales reales positivos (es decir, R(i, i) > 0, para todo i = 1, . . . , m). Si A ∈ Rn×m entonces tanto Q como R pueden tomarse con entradas reales y QT Q = Im . £ ¤ Demostraci´on. Expresemos A = a1 . . . am , con ai ∈ Cn . Si aplicamos el Algoritmo 2 a las columnas de A, obtenemos u1 = a1 , q1 = kuu11 k uk = ak − (< q1 , ak > q1 + . . . + < qk−1 , ak > qk−1 ) , qk =

uk kuk k

,

k = 2, . . . , m.

De aqu´ı despejamos el vector ak para obtener ak =< q1 , ak > q1 + . . . + < qk−1 , ak > qk−1 + k uj k qk . As´ı pues, si definimos la matriz R cuya entrada (i, j) es ½ < qi , aj > , i 6= j R(i, j) = , k ui k , i=j tenemos que A = QR,

£ ¤ donde, si A tiene columnas linealmente independientes, entonces la matriz Q = q1 . . . qm tiene columnas ortonormales, la matriz R es triangular superior (y cuadrada de tama˜no m × m) y R(i, i) =k ui k> 0. Analicemos ahora el caso en que las columnas de A no son linealmente independientes. Como ya hemos visto al final de la Secci´on 1.3, ak ∈ Gen(a1 , . . . , ak−1 ) si y s´olo si qk = 0. Supongamos que el rango de A es r, lo que significa que A tiene un total de m − r columnas linealmente dependientes. Por tanto, hay un total de m − r vectores qi que son nulos. El procedimiento para hallar la factorizaci´on QR de A en este caso es el siguiente. Una vez conclu´ıdo el Algoritmo 2 , los vectores ortonormales no nulos que hemos obtenido, qi1 , . . . , qir , se completan hasta obtener un conjunto ortonormal de Cn de m vectores, qi1 , . . . , qir , qej1 , . . . , qejm−r . Ahora, los vectores nulos se van sustituyendo por los vectores nuevos (los denotados por qej ). Por ejemplo, si el primer vector nulo es qs = 0 entonces tendremos la sucesi´on q1 , q2 ,£. . . , qs−1 , qej1 , . ¤. . . Procedemos de esta manera hasta completar el total de m vectores. Si Q = q1 . . . qm es la matriz obtenida por el Algoritmo 2 (incluyendo las e es la matriz que se obtiene al sustituir las columnas nulas por los vectores columnas nulas) y Q qej , entonces e QR = QR. La justificaci´on de esta identidad es la siguiente. Si qs = 0 es uno de los vectores nulos, entonces < qs , ai >= 0, para todo i = 1, . . . , m, lo que significa que la s-´esima fila de R en nula y, por lo tanto, al cambiar qs por otro vector (en particular, un vector qej ) el producto por la matriz R no var´ıa (n´otese que, en el producto QR las entradas de la s-´esima fila de R multiplican a las de la columna s-´esima de Q).

1. PRODUCTOS ESCALARES

18

Finalmente, si la matriz A tiene entradas reales, entonces tanto las coordenadas de los vectores qi , como los coficientes que produce el Algoritmo 2 son n´umeros reales. La demostraci´on del Teorema 1.4.3 nos muestra que la factorizaci´on QR de A se puede obtener mediante el Algoritmo 2 aplicado a los vectores columna de A. Este procedimiento nos da, en realidad, UNA factorizaci´on QR, porque la factorizaci´on QR de una matriz dada ´ NO ES UNICA en general. No obstante, hay un caso particular en el que s´ı podemos afirmar que es u´ nica. Este caso es el objeto del siguiente ejercicio. Ejercicio 1.4.4 Demostrar que si n = m y A es invertible, entonces la factorizaci´on QR de A es u´ nica. (Indicaci´on: suponer, por reducci´on al absurdo, que hay dos factorizaciones distintas de A = QR = Q0 R0 y deducir, usando que A es invertible, que Q = Q0 y R = R0 ). A continuaci´on, vamos a ver dos ejercicios resueltos en los que se pide calcular la factorizaci´on QR de una matriz A. En el primero de ellos la matriz A tiene columnas linealmente independientes, mientras que en el segundo no. Problema Resuelto 1.4.5 Hallar la factorizaci´on QR de   1 1 2 A =  1 −1 1  2 1 2 Soluci´on. Aplicamos el Algoritmo 2 a las columnas de la matriz A (que denotaremos por a1 , a2 , a3 ).   1 a1 1  √ u1 = a1 , q1 = ka1 k = 6 1  2        2 1 1 2 u2 = a2 − < a2 , q1 > q1 =  −1  − 13  1  = 13  −4  , q2 = kuu22 k = √121  −4  1 2  1    1    2 1 2 9 2  1  −4  = 14 3 , u3 = a3 − < a3 , q1 > q1 − < a3 , q2 > q2 =  1  − 76  1  − 21 2 2 1 −6   3 q3 = kuu33 k = √114  1  . −2 Ahora despejamos los vectores ai en las expresiones anteriores: a1 =k a1 k q1 a2 =< a2 , q1 > q1 + k u2 k q2 a3 =< a3 , q1 > q1 + < a3 , q2 > q2 + k u3 k q3 , lo que, expresado matricialmente, nos da la factorizaci´on QR de A √ √  1  √  7 6 √ √2 √3 6 √36 6 21 14 √6   21 2 21  . √1  A =  √16 √−4 21 14   0 3 21  √ 1 −2 2 3 14 √ √ √ 0 0 6 21 14 14

1.4 Factorizaci´on QR

19

Problema Resuelto 1.4.6 Hallar la factorizaci´on QR de la matriz   1 1 0 0  2 1 −1 1     A=  −2 1 3 0  .  0 1 1 0  0 0 0 1 Soluci´on. Aplicamos el Algoritmo 2 al conjunto de a1 , a2 , a3 , a4 ).   1  2    a1 1 u1 = a1 , q1 = ka1 k = 3  −2    0  0      1 1  1       1 2  1     u2 = a2 − < a2 , q1 > q1 =   1  − 9  −2  = 9   1   0   0 0   0  −1     u3 = a3 − < a3 , q1 > q1 − < a3 , q2 > q2 =   3 +  1  0   0  0    u3  q3 = ku3 k =   0   0  0

las columnas de A (que denotamos

 8 7   11  , 9  0  1  2  8 9  −2  0 0



 8  7    u2 1  q2 = ku2 k = 3√35   11   9  0      8 0       1 7   0   −  11  =  0  ,  9     9   0   0 0



 1       2 2     u4 = a4 − < a4 , q1 > q1 − < a4 , q2 > q2 − < a4 , q3 > q3 =   − 9  −2    0    0           0 8 1 8 −2    7   1    7      2 2   1 2    1 1 7          − 3√35 3√35  11  =  0  − 9  −2  − 45  11  = 5  1  ,  9   0   0   9   −1  1 0 0 0 5   −2  2    u4 1  q4 = ku4 k = √35   1 .  −1  5 Como q3 = 0 debemos substituir el vector q3 por otro vector, qe3 , de manera que el conjunto {q1 , q2 , qe3 , q4 } sea ortonormal. Para hallar un vector unitario que sea ortogonal a 

0 1 0 0 1



1. PRODUCTOS ESCALARES

20

eneo U T x = 0, donde U = £Gen(q1 , q2 , q4 ) ¤= Gen{u1 , u2 , u4 } resolvemos el sistema homog´ u1 u2 u4 . Para ello pasamos la matriz del sistema, U T , a su forma escalonada reducida     1 2 −2 0 0 1 0 0 2/3 −4/3 1 . U T =  8 7 11 9 0  ∼  0 1 0 0 −2 2 1 −1 5 0 0 1 1/3 1/3 A partir de aqu´ı se obtiene de forma inmediata una soluci´on    4   −3    u e3 1    u e3 =  −1  , qe3 = =√  k u e k 35  3   0  3

4 −3 −1 0 3

   .  

Ahora despejamos los vectores ai en las expresiones de la p´agina anterior: a1 a2 a3 a4

=k a1 k q1 =< a2 , q1 > q1 + k u2 k q2 =< a3 , q1 > q1 + < a3 , q2 > q2 + k u3 k q3 =< a3 , q1 > q1 + < a3 , q2 > q2 + 0 · qe3 =< a4 , q1 > q1 + < a4 , q2 > q2 + < a4 , q3 > q3 + k u4 k q4 =< a4 , q1 > q1 + < a4 , q2 > q2 + 0 · qe3 + k u4 k q4 .

Sustituyendo los coeficientes y expresando matricialmente las igualdades anteriores llegamos a la siguiente factorizaci´on QR de A  1  −2 √8 √4 √   3 35 35 −8 2 3 √13  2 3√735 √−3 2  3 3 √ √ √  3 3 35  35 35 35 35 35   −2 0 11 −1 1  3 3 15    √ √ √ A =  3 3 35 .  35 35  0 0 0  −1   0  0 √ √3 √ 0  35 35 35  0 0 0 5 √3 √5 0 0 35 35

´ QR completa y factorizacion ´ QR reducida Factorizacion La factorizaci´on QR completa de A ∈ Kn×m (con n ≥ m) consiste en a˜nadir a la matriz Q de la factorizaci´on QR (que, a partir de ahora, llamaremos factorizaci´on QR reducida) un total de n − m columnas de modo que la matriz resultante, Qc , tenga un total de n columnas ortonormales. La matriz Qc es, por tanto, cuadrada n × n, invertible y con todas sus columnas ortonormales. Estas matrices se llaman unitarias. Tienen una gran importancia y ser´an objeto de estudio en el Cap´ıtulo 2. Pero la factorizaci´on QR completa no consiste s´olo en alterar la matriz Q pues, en tal caso, no se podr´ıa hacer el producto de la nueva matriz Qc con la antigua matriz R. Para que las dimensiones permitan hacer el producto, la matriz R tambi´en se modifica a˜nadi´endole un total de n − m filas nulas en las mismas posiciones en que se han a˜nadido las columnas a la matriz Q (lo natural es hacerlo al final en ambas matrices). De esta manera, se obtiene una nueva matriz Rc de modo que el producto de las dos nuevas matrices es el mismo que el de las dos matrices antiguas, es decir A = QR = Qc Rc .

1.4 Factorizaci´on QR

21

La relaci´on entre la QR completa y la QR reducida de A est´a descrita gr´aficamente en la siguiente figura. R =

A

Q

Qc

Rc

´ de la factorizacion ´ QR a la resolucion ´ de sistemas lineales Aplicacion La factorizaci´on QR puede aplicarse a la resoluci´on del sistema lineal Ax = b ,

(1.5)

donde A ∈ Kn×m y b ∈ Kn , de la siguiente manera. Dada la factorizaci´on QR de A, A = QR, el sistema (1.5) es equivalente al sistema (QR)x = b ⇔ Rx = Q∗ b. El u´ ltimo sistema Rx = Q∗ b es un sistema escalonado (la matriz R es triangular superior) y, por lo tanto, f´acil y r´apido de resolver. As´ı pues, tenemos el siguiente m´etodo para resolver el sistema (1.5). Resoluci´on de Ax = b usando QR 1. Calcular una factorizaci´on QR de A: A = QR. 2. Calcular y = Q∗ b. 3. Resolver Rx = y (en x). ´ ´ QR con Matlab Calculo de la factorizacion El comando Matlab que calcula la factorizaci´on QR completa de una matriz A es qr(A) mientras que el comando que da la factorizaci´on QR reducida es qr(A,0) . Para obtener las dos matrices Q y R escribiremos [q r]=qr(A) o´ [q r]=qr(A,0) (por supuesto, habiendo introducido previamente la matriz A) y el programa nos devolver´a dos matrices: una matriz q , que es la matriz Q, y una segunda matriz, r , que es la matriz R. Para m´as informaci´on se puede introducir la orden help qr. Ejemplo 1.4.7 La factorizaci´on QR de la matriz A del Problema Resuelto 1.4.5 se obtiene en Matlab introduciendo la siguiente secuencia: primero introducimos la matriz A y despu´es el comando qr(A) . >>A=[1 1 2; 1 -1 1;2 1 2]; >>[q r]=qr(A)

1. PRODUCTOS ESCALARES

22 q= -0.4082 -0.4082 -0.8165

0.4364 -0.8729 0.2182

-0.8018 -0.2673 0.5345

r= -2.4495 0 0

-0.8165 1.5275 0

-2.8577 0.4364 -0.8018

En este caso, al ser A una matriz cuadrada, la factorizaci´on QR reducida coincide con la factorizaci´on QR completa. Para la matriz del Problema Resuelto 1.4.6 tenemos >>A=[1 1 0 0; 2 1 -1 1;-2 1 3 0; 0 1 1 0;0 0 0 1]; >>[q r]=qr(A) q= -0.3333 -0.4507 -0.7778 0.0604 -0.2777 -0.6667 -0.3944 0.5938 0.0463 -0.2127 0.6667 -0.6198 0.2049 0.0765 -0.3515 0 -0.5071 -0.0209 -0.1833 0.8419 0 0 0 0.9771 0.2127 r= -3.0000 -0.3333 2.6667 -0.6667 0 -1.9720 -1.9720 -0.3944 0 0 0.0000 0.5938 0 0 0 1.0234 0 0 0 0 La factorizaci´on QR reducida es >>[q r]=qr(A,0) q= -0.3333 -0.4507 -0.6667 -0.3944 0.6667 -0.6198 0 -0.5071 0 0 r= -3.0000 -0.3333 0 -1.9720 0 0 0 0

-0.7778 0.5938 0.2049 -0.0209 0

0.0604 0.0463 0.0765 -0.1833 0.9771

2.6667 -1.9720 0.0000 0

-0.6667 -0.3944 0.5938 1.0234

Es importante observar que las matrices que nos da el programa Matlab son las mismas que hemos obtenido anteriormente en el Problema Resuelto 1.4.5 excepto por algunos cambios de signo. Esto se debe a que el programa Matlab no calcula necesariamente la matriz R que

1.4 Factorizaci´on QR

23

tiene todas sus entradas diagonales no negativas (obs´ervese que hay, efectivamente, entradas diagonales negativas en la matriz R). Para obtener la factorizaci´on QR que hemos estudiado (es decir, con las entradas diagonales de R no negativas) a partir de la que nos da el programa Matlab procederemos de la siguiente manera. bR b a la factoriDenotaremos por QR a la factorizaci´on que da el programa Matlab y por Q b zaci´on que queremos encontrar (en la que las entradas diagonales de R son nonegativas). Como es habitual, qi denota a la columna i-´esima de Q, y R(i, j) a la entrada (i, j) de R. Supongamos que R(i, i) < 0. Entonces, tendremos que cambiar de signo esta entrada, es decir: b i) = −R(i, i) . R(i, Este cambio implica un cambio de signo en la columna qi , es decir: qbi = −qi . Pero, ahora, este nuevo cambio tiene consecuencia en los siguientes pasos (recu´erdese el m´etodo de Gram-Schmidt) en los que aparece la columna qi . Concretamente, e´ sta aparece multiplicada por las entradas de la i-´esima fila de R. Por tanto, hay que cambiar de signo esta fila, es decir: b j) = −R(i, j) , j = i + 1, ..., m . R(i, En resumen: Por cada entrada R(i, i) < 0 de R hay que cambiar de signo: La i-´esima fila de R. La i-´esima columna de Q. Interpretaci´on matricial: Lo que estamos haciendo en el procedimiento anterior es multiplicar la Q y la R por la matriz diagonal   ±1   .. S= ,  . ±1 m×m donde el signo de la i-´esima entrada diagonal es (+) si R(i, i) ≥ 0 y (−) si R(i, i) < 0. Claramente S 2 = Im , por lo que A = QR = (QS)(SR) b = QS y R b = SR. Recu´erdese que multiplicar a de manera que, con la notaci´on anterior, Q la derecha es hacer operaciones con las columnas de una matriz, mientras que multiplicar a la izquierda es hacerlo con las filas. Eso significa que QS cambia de signo las columnas de Q y SR cambia de signo las filas de R. Esta interpretaci´on matricial nos permite generalizar el procedimiento anterior al caso complejo, en el que algunas de las entradas diagonales de R son n´umeros complejos no reales. En este caso, har´ıamos lo siguiente: A = QR = (QT1 )(T2 R) ,

1. PRODUCTOS ESCALARES

24 donde   T1 =  

R(1,1) |R(1,1)|

 ..

. R(m,m) |R(m,m)|

  

 ,

 T2 =  

|R(1,1)| R(1,1)

 ..

.

m×m

|R(m,m)| R(m,m)

  

,

m×m

b = QT1 , R b = T2 R. Por supuesto, este u´ ltimo procede modo que T1 T2 = Im , y tomar´ıamos Q dimiento incluye el anterior de las entradas negativas ya que, en ese caso R(i, i)/|R(i, i)| = ±1 = |R(i, i)|/R(i, i). Finalmente, observemos que la factorizaci´on QR que da Matlab en el caso del Problema Resuelto 1.4.6 es distinta de la que nosotros estudiamos (aparte de los signos). Este ejemplo demuestra que, en el caso de que las columnas de A no sean linealmente independientes, la factorizaci´on QR de A no es u´ nica.

1.5 Problemas de m´ınimos cuadrados Sean A ∈ Cn×m , b ∈ Cn y x ∈ Cm un vector de inc´ognitas. Supongamos que el sistema lineal Ax = b, (1.6) no tiene soluci´on. En esta situaci´on es natural buscar un vector, x b, que sea “lo m´as parecido” a una soluci´on. Naturalmente que lo primero que hay que concretar es qu´e entendemos por “lo m´as parecido” . Para ello usamos el concepto de distancia que hemos introducido en la Definici´on 1.1.20. Concretamente: dos vectores ser´an parecidos si su distancia es peque˜na (diremos, en este caso, que ambos vectores son pr´oximos). As´ı pues, para el sistema lineal (1.6) razonamos de la siguiente manera: si x b es un vector pr´oximo a una soluci´on de (1.6) entonces Ab x es pr´oximo al vector b o, lo que es igual, Ab x − b es pr´oximo al vector 0Cn . Eso significa, en virtud de la Definici´on 1.1.20, que la norma de Ab x − b es peque˜na. As´ı pues, nuestro problema es el siguiente: Minimizar : k Ab x − b k2 , para x b ∈ Cm . Gr´aficamente, la situaci´on es: b 

col(A)

El vector b no est´a en el espacio col(A), y por eso el sistema (1.6) no tiene soluci´on. Si llamamos bb = Ab x, entonces bb ∈ col(A), y el problema se puede plantear as´ı:

1.5 Problemas de m´ınimos cuadrados

25

Encontrar el vector de col(A), bb, m´as pr´oximo a b. En general, si en lugar del espacio col(A) consideramos un subespacio vectorial cualquiera, llegamos a la siguiente definici´on. Definici´on 1.5.1 Sea Kn con el producto escalar usual y k · k2 la norma asociada. Dado un subespacio vectorial W y un vector cualquiera v ∈ Kn , la mejor aproximaci´on de W a v en sentido de m´ınimos cuadrados es el vector vb tal que k v − vb k2 = ´ınf k v − w k2 . w∈W

La distancia de v a W , que denotaremos por d(v, W ), es el n´umero real d(v, W ) = d(v, vb). El siguiente Teorema nos dice que la mejor aproximaci´on de W a v est´a dada por la proyecci´on ortogonal de v sobre W . Teorema 1.5.2 (Soluci´on al problema de m´ınimos cuadrados) La mejor aproximaci´on de W a v en el sentido de m´ınimos cuadrados es el vector vb = proyW (v). El resultado del Teorema 1.5.2 no debe resultar sorprendente. En realidad, ya conocemos un caso particular de dicho teorema: la distancia de un punto a una recta. Como ya sabemos, la distancia de un punto P a una recta r se define como la menor distancia de P a cualquier punto de r, y est´a dada por la distancia de P a Pb, donde Pb es la proyecci´on ortogonal de P sobre r. En virtud del Teorema 1.5.2, el procedimiento est´andar para hallar la soluci´on al problema de m´ınimos cuadrados consiste en calcular la proyecci´on de v sobre W . Recordemos que para ello necesitamos calcular una base ortogonal de W , lo que en general es un procedimiento muy costoso si se calcula a mano. A continuaci´on, vamos a estudiar un procedimiento menos costoso, que no requiere el c´alculo de una base ortogonal para obtener la proyecci´on. Las ecuaciones normales En este apartado W es un subespacio vectorial de V, del que conocemos una base B = {v1 , . . . , vk } (no necesariamente ortogonal) y v es un vector cualquiera de V. Recordemos que el vector v ⊥ = v − proyW (v) es ortogonal a W . Si expresamos la proyecci´on ortogonal en la base B proyW (v) = x1 v1 + . . . + xk vk entonces, para cada j = 1, . . . , k, tenemos que < vj , v −

k X

xi vi >= 0,

i=1

lo que equivale a < vj , v >=

k X i=1

xi < vj , vi > .

1. PRODUCTOS ESCALARES

26

La igualdad anterior se expresa matricialmente de la siguiente manera z    

G

}| {    < v1 , v1 > . . . < v1 , vk > x1 < v1 , v >     < v2 , v1 > . . . < v2 , vk >   x2   < v2 , v >   ..  =  .. .. ..  .   . . . < vk , v1 > . . . < vk , vk > xk < vk , v >

   . 

Si V = Kn y el producto escalar < ·, · > es el producto escalar usual, la igualdad matricial anterior se puede expresar en la forma A∗ Ax = A∗ v

(1.7)

o bien Gx = A∗ v, £ ¤ donde A = v1 . . . vk es la matriz cuyas columnas son los vectores de la base B y G es la matriz de Gram del producto escalar < ·, · > en la base B, que hemos definido al final de la Secci´on 1.1. Las ecuaciones (1.7) se conocen con el nombre de ecuaciones normales para las coordenadas de la proyecci´on y s´olo involucran a los vectores de la base B (ordenados por columnas en la matriz A) y al vector v. Si G = A∗ A es invertible, podemos despejar x en (1.7) para obtener la siguiente f´ormula de la proyecci´on ortogonal en funci´on de una base cualquiera (no necesariamente ortogonal) de W proyW (v) = A(A∗ A)−1 A∗ v (1.8) La pregunta que queda por resolver es: ¿cu´ando la matriz A∗ A es invertible? El siguiente resultado nos proporciona una condici´on suficiente. Lema 1.5.3 Si las columnas de A son linealmente independientes entonces la matriz A∗ A es invertible. Por lo tanto, si partimos de una base B = {v1 , . . . , vk } cualquiera de W , la matriz A∗ A es invertible, porque las columnas de A son linealmente independientes. En este caso, la f´ormula (1.8) nos permite calcular la proyecci´on de v sobre W sin necesidad de ortogonalizar B. Resumimos este procedimiento alternativo para el c´alculo de la proyecci´on ortogonal (con el producto escalar usual) en el siguiente cuadro C´ alculo de proyW (v) (con el producto escalar usual) 1. Calcular una base de W : B = {v1 , . . . , vk } £ ¤ 2. Construir la matriz A = v1 . . . vk 3. Aplicar la f´ormula: proyW (v) = A(A∗ A)−1 A∗ v

1.5 Problemas de m´ınimos cuadrados

27

´ QR M´ınimos cuadrados usando la factorizacion £ ¤ En el caso en que las columnas de la matriz A = v1 . . . vk sean linealmente independientes (en otras palabras, cuando los vectores v1 , . . . , vk formen una base de W ), podemos usar la factorizaci´on QR de A para hallar una f´ormula m´as elemental de la proyecci´on ortogonal sobre W . Por la propiedad iv) del Lema 2.1.1, que veremos m´as adelante, tenemos que A∗ = R∗ Q∗ , luego A(A∗ A)−1 A∗ = QR(R∗ Q∗ QR)−1 R∗ Q∗ = QR(R∗ R)−1 R∗ Q∗ , y, dado que R es tambi´en invertible (pues lo es A), la anterior expresi´on se simplifica a QQ∗ . Por lo tanto, nos queda la f´ormula elemental proyW (v) = QQ∗ v Esta f´ormula es la misma que hemos obtenido al final de la Secci´on 1.3, pues las columnas de Q son una base ortonormal del espacio de columnas de A. Hay que destacar que la simplicidad de la f´ormula anterior se basa en las propiedades de la matriz Q, cuyo c´alculo equivale al c´alculo de una base ortogonal de W . Dicho c´alculo es costoso si se hace a mano pero, por el contrario, la matriz Q (y la factorizaci´on QR) se pueden obtener de una manera r´apida y precisa con el programa Matlab. La soluci´on que ofrece el programa mediante este procedimiento es mucho m´as precisa que si se resuelve el problema (con el mismo programa) usando las ecuaciones normales. ´ de Ax = b por m´ınimos cuadrados Solucion El Teorema 1.5.2 nos dice que bb = proy col(A) (b) es la mejor aproximaci´on de col(A) al vector b. Por tanto, para resolver (1.6) por m´ınimos cuadrados debemos resolver Ab x = bb . (1.9) ´ Este sistema tiene soluci´on porque bb ∈ col(A), pero la soluci´on puede no ser unica. En este apartado nos interesa el caso en el que hay soluci´on u´ nica. En la Secci´on 3.2 estudiaremos un m´etodo para el c´alculo de la soluci´on o´ ptima cuando hay m´as de una soluci´on. ´ Recordamos el siguiente resultado de Algebra elemental que es clave para nuestro problema. Teorema 1.5.4 La soluci´on del sistema lineal (1.9) es u´ nica si y s´olo si la matriz de coeficientes A tiene columnas linealmente independientes. En el caso de que las columnas de A sean linealmente independientes, en virtud del Lema 1.5.3, la matriz A∗ A es invertible y, por tanto, la soluci´on por m´ınimos cuadrados del sistema (1.6) es la soluci´on del sistema Ab x = A(A∗ A)−1 A∗ b.

1. PRODUCTOS ESCALARES

28

Si multiplicamos en ambos t´erminos a la izquierda por A∗ obtenemos A∗ Ab x = A∗ b. As´ı pues, hemos llegado al siguiente resultado. Teorema 1.5.5 Si la matriz A tiene columnas linealmente independientes entonces el sistema lineal (1.6) tiene soluci´on u´ nica por m´ınimos cuadrados y dicha soluci´on, x b, es la soluci´on del sistema lineal (A∗ A)b x = A∗ b. £ ¤T Problema Resuelto 1.5.6 Calcular la proyecci´on ortogonal del vector v = 1 −1 1 sobre el plano Π ≡ x − y + 2z = 0 y usar dicha proyecci´on para resolver, por m´ınimos cuadrados, el sistema lineal   x + 3y = 1 x + y = −1 .  −y = 1 Soluci´on. Para calcular la proyecci´on ortogonal usaremos la f´ormula (1.8), para lo que nece£ ¤T sitamos una base del plano Π. El vector normal al plano es n = 1 −1 2 , por lo que una base del plano estar´a formada por dos vectores perpediculares a n. Podemos tomar, por ejemplo, los vectores que determinan los coeficientes del sistema del enunciado, lo que nos da la matriz   1 3 A =  1 1 . 0 −1 Ahora, aplicando la f´ormula (1.8) 

   5 1 −2 1/3 1 proyΠ (v) = A(A∗ A)−1 A∗ v =  1 5 2  v =  −1/3  . 6 −2 2 2 −1/3 El vector v no pertenece al plano Π, por lo que el sistema lineal no tiene soluci´on. Como las columnas de A son linealmente independientes, la soluci´on por m´ınimos cuadrados es u´ nica y coincide con la soluci´on, x b, del sistema   1/3 Ab x = proyΠ (v) =  −1/3  . −1/3 Dicha soluci´on se obtiene f´acilmente escalonando la matriz A y es · ¸ −2/3 x b= . 1/3

2 ESTUDIO ESPECTRAL DE ENDOMORFISMOS

2.1 Matrices unitarias y ortogonales En este tema estudiaremos dos clases fundamentales de matrices: las matrices unitarias y las matrices ortogonales. La relevancia de estas clases de matrices radica en que preservan el producto escalar usual y, en consecuencia, la norma asociada. Las matrices unitarias preservan el producto escalar herm´ıtico y las ortogonales preservan el producto escalar real. Este hecho tiene implicaciones tanto desde el punto de vista de la geometr´ıa como desde el punto de vista del c´alculo num´erico. El primer punto ser´a estudiado en este tema, mientras que el segundo se tratar´a en el Tema 6. Comenzamos este tema con unas propiedades elementales de la traspuesta conjugada de matrices. Su demostraci´on es inmediata. Lema 2.1.1 Sean A, B dos matrices cualesquiera (de las dimensiones adecuadas) y α ∈ C. Entonces i) (A∗ )∗ = A. ii) (A + B)∗ = A∗ + B ∗ . iii) (αA)∗ = αA∗ . iv) (AB)∗ = B ∗ A∗ . Recordemos que ya hemos usado la propiedad iv) en el Tema 1 para obtener la f´ormula de la proyecci´on en t´erminos de la matriz Q de la factorizaci´on QR. En la siguiente definici´on introducimos las clases fundamentales de este tema y tambi´en recordamos otras dos definiciones ya conocidas que tienen, como veremos, una interpretaci´on relevante desde el punto de vista de los productos escalares. 29

2. ESTUDIO ESPECTRAL DE ENDOMORFISMOS

30

Definici´on 2.1.2 i) Una matriz U ∈ Cn×n se dice unitaria si U ∗ U = In . En este caso se tiene tambi´en U U ∗ = In y, por tanto, U es invertible y U ∗ = U −1 . ii) Una matriz P ∈ Rn×n se dice ortogonal si P T P = In . En este caso se tiene tambi´en P P T = In y, por tanto, P es invertible y P ∗ = P −1 . iii) Una matriz A ∈ Cn×n se dice herm´ıtica (o autoadjunta) si A∗ = A. iv) Una matriz A ∈ Cn×n se dice sim´etrica si AT = A. Observaci´on 2.1.3 N´otese que, seg´un la definici´on que hemos introducido, una matriz ortogonal es real. No obstante, tambi´en se puede definir el concepto de matriz ortogonal para matrices complejas (la definici´on ser´ıa exactamente la misma, pero permitiendo que la matriz tenga entradas complejas). Sobre este particular no hay un criterio universal, pero s´ı mayoritario. La mayor´ıa de los autores entienden por matriz ortogonal (a secas) una matriz real, mientras que se utiliza el t´ermino matriz compleja ortogonal para las matrices con entradas complejas. En este curso las matrices ortogonales complejas no ser´an necesarias, por lo que hemos preferido definir las matrices ortogonales como matrices reales. En este apartado nos centraremos en las matrices ortogonales y unitarias. S´olo en un ap´endice final daremos la interpretaci´on de las matrices sim´etricas y herm´ıticas a trav´es de los productos escalares/herm´ıticos. En el siguiente lema se enuncian algunas de las propiedades b´asicas de las matrices unitarias y ortogonales. Lema 2.1.4 (Propiedades de las matrices unitarias y ortogonales) En Cn consideramos el producto herm´ıtico usual < u, v >= u∗ v. Dada una matriz U ∈ Cn×n i) U es unitaria si y s´olo si las columnas de U son ortonormales, ii) U es unitaria si y s´olo si las filas de U son ortonormales, iii) Si U es unitaria entonces |det(U )| = 1, iv) Si U es unitaria entonces < U u, U w >=< u, v >, para cualquier par de vectores u, v ∈ Cn . Es decir, las matrices unitarias conservan el producto herm´ıtico usual, en particular v) Si U es unitaria entonces k U v k2 =k v k2 , para todo v ∈ Cn . En Rn consideramos el producto escalar usual < u, v >= uT v. Dada una matriz P ∈ Rn×n vi) P es ortogonal si y s´olo si las columnas de P son ortonormales, vii) P es ortogonal si y s´olo si las filas de P son ortonormales, viii) Si P es ortogonal entonces det(P ) = ±1,

2.2 Matrices de Householder

31

ix) Si P es unitaria entonces < P u, P w >=< u, v >, para cualquier par de vectores u, v ∈ Rn . Es decir, las matrices ortogonales conservan el producto escalar usual, en particular x) Si P es ortogonal entonces k P v k2 =k v k2 , para todo v ∈ Rn . xi) Si U, V ∈ Cn×n son matrices unitarias entonces U V es unitaria. xii) Si P, Q ∈ Rn×n son matrices ortogonales entonces P Q es ortogonal. Observaci´on 2.1.5 La matriz Q de la factorizaci´on QR de una matriz A real cuadrada es una matriz ortogonal (es cuadrada y tiene columnas ortonormales), mientras que si la matriz A es compleja cuadrada entonces la matriz Q es unitaria. En las propiedades iv) y ix) radica gran parte de la importancia que tienen las matrices unitarias y ortogonales tanto desde el punto de vista geom´etrico como desde el punto de vista num´erico. Seg´un estas propiedades, las matrices unitarias/ortogonales, o sus transfromaciones asociadas, conservan tanto el tama˜no (norma) de los vectores como su “posici´on relativa”, dada por sus a´ ngulos. Este hecho no s´olo tiene importancia geom´etrica (que explotaremos en el Tema 5), sino tambi´en num´erica. Cuando operamos con matrices que no son unitarias/ortogonales corremos el riesgo de introducir grandes errores en el c´alculo de vectores. Por ejemplo, al multiplicar dos vectores perpendiculares (o casi perpendiculares) y, por lo tanto, lo m´as alejados posible de la dependencia lineal, por una matriz cualquiera invertible no podemos obtener dos vectores paralelos (debido a la invertibilidad de la matriz) pero s´ı dos vectores casi paralelos, es decir, casi dependientes. Es decir, dos vectores de direcciones muy diferentes pueden transformarse en vectores de direcciones parecidas al multiplicarlos por una matriz cualquiera. Esto, en cambio, no puede ocurrir al multiplicarlos por una matriz unitaria/ortogonal, pues los a´ ngulos se conservan. Por este motivo las matrices unitarias/ortogonales son preferibles desde el punto de vista num´erico.

2.2 Matrices de Householder En esta secci´on vamos a estudiar una clase importante de matrices unitarias, conocidas como matrices de Householder, en honor al matem´atico norteamericano Alston Householder (1904-1993). Estas matrices se utilizan en el c´alculo de la factorizaci´on QR y, adem´as, tienen una interpretaci´on geom´etrica elemental, como veremos a continuaci´on. Definici´on 2.2.1 (Matriz de Householder) Dado un vector v ∈ Rn , la matriz de Householder asociada a v es la matriz 2 Hv = In − T vv T . v v Conviene observar que Hv ∈ Rn×n . As´ı, la matriz de Householder Hv puede verse como una transformaci´on lineal de Rn , dada por Hv : Rn −→ Rn v 7→ Hv (w) = Hv · w.

2. ESTUDIO ESPECTRAL DE ENDOMORFISMOS

32

A estas transformaciones las llamaremos transformaciones de Householder. A continuaci´on, vamos a dar una interpretaci´on geom´etrica de dichas aplicaciones. Para entenderla, recordamos que un hiperplano en el espacio Rn es un subespacio vectorial de dimensi´on n − 1. Su nombre se debe a que un hiperplano es la generalizaci´on a Rn del concepto de plano en R3 . Por ejemplo, en el plano R2 los hiperplanos son las rectas. Tambi´en conviene tener en cuenta que, en virtud del Lema 1.1.17, el subespacio ortogonal de una recta (o un vector) en Rn es un hiperplano. De hecho, todo hiperplano es el subespacio ortogonal de alg´un vector. La matriz de Householder Hv representa una simetr´ıa en Rn respecto al hiperplano ortogonal al vector v. La comprobaci´on de este hecho es muy sencilla: basta darse cuenta de que la matriz que aparece en el segundo t´ermino de Hv (sin el factor 2), Pv =

1 vv T , vT v

es la matriz asociada a la proyecci´on sobre el espacio generado por v, Gen(v), pues, para cualquier vector w ∈ V vT w < v, w > 1 v, Pv (w) = T vv T w = T v = v v v v < v, v > donde < ·, · > denota el producto escalar usual. Esta es, precisamente, la f´ormula de la proyecci´on de w sobre Gen(v) (comp´arese con la f´ormula (1.4) para un u´ nico vector u1 ). Por tanto Hv (w) = (In − 2Pv )(w) = w − 2 · proyGen(v) (w) es el vector sim´etrico de w con respecto al hiperplano ortogonal a Gen(v). Gr´aficamente proyGen(v) (w) 6

w 

Gen(v)⊥ Gen(v)

^ Hv (w)

Lema 2.2.2 (Propiedades de las matrices de Householder) Sea v ∈ Rn y Hv la matriz de Householder asociada a v. Entonces i) Hv es sim´etrica. ii) Hv es ortogonal. iii) Hαv = Hv , para todo α ∈ R − {0}.

2.2 Matrices de Householder

33

iv) det(Hv ) = −1. v) La transformaci´on lineal asociada a Hv deja fijos • La recta Gen(v). • Todos los vectores del hiperplano Gen(v)⊥ . vi) Dados dos vectores x, y ∈ Rn con k x k2 =k y k2 , Hx−y (x) = y. La interpretaci´on geom´etrica del apartado vi) del Lema 2.2.2 es la siguiente: x−y Los vectores x, y son sim´etricos respecto del espacio ortogonal al vector x − y.

o

x

Gen(x − y)⊥

 - y

Problema Resuelto 2.2.3 Calcular la matriz de Householder que transforma el vector £ ¤T x = 3 0 −4 en un vector paralelo al primer vector can´onico e1 . Soluci´on. Seg´un la propiedad vi) del Lema 2.2.2, la transformaci´on de Householder llevar´a el vector x a un vector de la misma norma que x. El vector paralelo a e1 de la misma norma que x es k x k2 e1 = 5e1 . Por tanto, la transformaci´on de Householder est´a asociada al vector £ ¤T v = x − 5e1 = −2 0 −4 . Se trata, por lo tanto, de la transformaci´on asociada a la matriz   3/5 0 −4/5 2 T  0 1 0 . Hv = Iv − vv = 20 −4/5 0 −3/5 El Problema Resuelto 2.2.3 es la base del procedimiento que vamos a describir a continuaci´on para el c´alculo de la factorizaci´on QR de una matriz, y que explota la propiedad ii) del Lema 2.2.2. Teorema 2.2.4 Dada una matriz A ∈ Rn×m con n ≥ m, existen m matrices de Householder, H1 , . . . , Hm de modo que Hm · · · H1 A = R, donde R es triangular superior con todas sus entradas diagonales no negativas. Por tanto A = QR, con Q una matriz unitaria que es el producto de m matrices de Householder.

34

2. ESTUDIO ESPECTRAL DE ENDOMORFISMOS

Demostraci´on. Explicaremos esquem´aticamente el procedimiento, que consta de m pasos, en cada uno de los cuales se aplica una transformaci´on de Householder a la matriz obtenida en el paso anterior. El primer paso consiste en hacer ceros en la primera columna de A por debajo de la entrada (1, 1). Si llamamos, como es habitual, a1 a la primera columna de A, sabemos que hay una matriz de Householder, H1 , tal que H1 a1 =k a1 k2 e1 . Por lo tanto · ¸ k a1 k2 ∗ H1 A = e2 , 0(n−1)×1 A donde ∗ es una fila de m − 1 entradas que no nos interesan, el 0 que hay bajo la entrada (1, 1) e2 es una matriz de tama˜no (n − 1) × (m − 1). N´otese es en realidad una columna de ceros y A que la entrada en la posici´on diagonal (1, 1) es no negativa. Esta entrada no se ver´a alterada en el resto del proceso. e2 , es El segundo paso consiste en hacer lo mismo que en el primero pero sobre la matriz A decir, aplicarle una transformaci´on de Householder que convierta en ceros todas las entradas e 2 , es de tama˜no (n − por debajo de la entrada (1, 1). La matriz de Householder asociada, H 1) × (n − 1). La matriz de Householder del segundo paso, que actuar´a sobre la matriz completa H1 A y es de tama˜no n × n, ser´a · ¸ 1 01×(n−1) H2 = . e2 0(n−1)×1 H El proceso contin´ua de la misma manera: en el paso j-´esimo partimos de una matriz de la forma ¸ · Rj−1 ∗ ej , 0(n−j+1)×(j−1) A donde Rj−1 es triangular superior de tama˜no (j − 1) × (j − 1) y ∗ representa una submaej la transformaci´on de triz que carece de relevancia en el proceso. Aplicaremos a la matriz A e j , que lleva la primera columna a un m´ultiplo de e1 . Esta matriz es de tama˜no Householder, H (n − j + 1) × (n − j + 1). La matriz Hj de tama˜no n × n que buscamos en este paso ser´a · ¸ Ij−1 0(j−1)×(n−j+1) Hj = . ej 0(n−j+1)×(j−1) H N´otese que esta matriz no altera la matriz Rj−1 obtenida en el paso anterior. Tras aplicar m−1 transformaciones de Householder, llegamos a una matriz de la forma ¸ · Rm−1 ∗ , 0(n−m+1)×1 e am donde e am es una u´ nica columna (si n = m es simplemente una entrada). Al aplicar el u´ ltimo paso llegamos a una matriz R, que es triangular superior y con entradas diagonales no negativas. e = Hm · · · H1 tal que QA e = R. Como Hi , para As´ı pues, hemos obtenido una matriz Q i = 1, . . . , m, es sim´etrica y ortogonal (apartados i) y ii) del Lema 2.2.2) tenemos que Hi−1 = e −1 = H −1 · · · H −1 = H1 · · · Hm . Si llamamos Q = H1 · · · Hm entonces HiT = Hi , luego Q m 1 Q es unitaria y A = QR, lo que demuestra la u´ ltima afirmaci´on del enunciado.

2.2 Matrices de Householder

35

Observaci´on 2.2.5 En el Teorema 2.2.4 consideramos matrices de Householder en sentido generalizado, incluyendo la matriz identidad, que ser´ıa la matriz de Householder asociada al vector nulo. Esto quiere decir que, si exclu´ımos esta matriz, el n´umero de matrices Hj del enunciado puede ser menor que m. N´otese que la matriz Q en el Teorema 2.2.4 es un producto de matrices sim´etricas (las matrices de Householder lo son, en virtud de la propiedad i) del Lema 2.2.2). En cambio, Q no es siempre sim´etrica. Este hecho nos indica que el producto de matrices sim´etricas no siempre da como resultado una matriz sim´etrica. Esto es precisamente lo que se pide demostrar en el siguiente ejercicio. Ejercicio 2.2.6 Demostrar, con un contraejemplo, que el producto de matrices sim´etricas no es siempre una matriz sim´etrica. El Teorema 2.2.4 nos proporciona un m´etodo para obtener la matriz R de la factorizaci´on QR de una matriz A mediante una serie de transformaciones de Householder aplicadas a dicha matriz A. La ventajas, desde el punto de vista num´erico, de operar con matrices de Householder (que son ortogonales) se ha comentado ya tras el Lema 2.1.4. Por las razones all´ı mencionadas el procedimiento descrito en el Teorema 2.2.4 para calcular la factorizaci´on QR es preferible al descrito en el Teorema 1.4.3 basado en el m´etodo de Gram-Schmidt. Precisamente el m´etodo del Teorema 2.2.4 es el que utiliza el programa Matlab para el c´alculo de la factorizaci´on QR. Recordemos que dicho programa no proporciona, en general, una matriz R con entradas diagonales no negativas. Ahora podemos dar una explicaci´on de este hecho. Como hemos visto en la demostraci´on del Teorema 2.2.4, en cada paso del proceso para llegar a la matriz R aplicamos una transformaci´on de Householder que lleva una columna, llam´emosla a, al vector k a k2 e1 . No obstante, lo importante del procedimiento es llegar a un vector que tenga todas las entradas nulas salvo la primera y, como queremos llegar a e´ l mediante matrices ortogonales, su norma ha de coincidir con la norma de a. El vector k a k2 e1 no es el u´ nico que cumple estas dos condiciones. De hecho, hay otro vector con estas caracter´ısticas, que es precisamente el vector opuesto al anterior, − k a k2 e1 . Por cuestiones num´ericas, hay ocasiones en que no es aconsejable elegir el primero, sino el segundo. Es en estos casos en los que aparece una entrada negativa en la diagonal de la matriz R. ´ de las matrices simetricas ´ Interpretacion y herm´ıticas Debido a la simetr´ıa del producto escalar real < ·, · >, la matriz de Gram de un producto escalar en una base cualquiera (v´ease la definici´on al final de la Secci´on 1.1) es una matriz sim´etrica con entradas reales (diremos sim´etrica real). A lainversa, aunque cualquier producto escalar queda determinado por los valores que asigna a los vectores de una base, no es verdad que cualquier matriz sim´etrica real es la matriz de Gram de un determinado producto escalar real en una base prefijada. Para ello la matriz ha de ser definida positiva (v´ease la Secci´on 4.3). En concreto, si B = {v1 , . . . , vn } es una base de V y A = [aij ]ni,j=1 es una matriz sim´etrica real definida positiva, definimos el producto escalar de los vectores vi , vj de la base como < vi , vj >= aij . De este modo, el producto escalar de dos vectores u, v ∈ V cualesquiera, expresados en la base B, est´a definido por < u, v >= [u]TB A[v]B .

2. ESTUDIO ESPECTRAL DE ENDOMORFISMOS

36

De forma an´aloga, la matriz de Gram de un producto herm´ıtico es una matriz herm´ıtica mientras que rec´ıprocamente, toda matriz herm´ıtica definida positiva H define un producto herm´ıtico en una base prefijada B, concretamente < u, v >= [u]∗B H[v]B . Resumiendo: Teorema 2.2.7 a) Sea V un espacio vectorial real y B una base de V. Cualquier matriz sim´etrica definida positiva A ∈ Rn×n es la matriz de Gram de un producto escalar sobre V en la base B, que est´a definido por < u, v >= [u]TB A[v]B . b) Sea V un espacio vectorial complejo y B una base de V. Cualquier matriz herm´ıtica definida positiva H ∈ Cn×n es la matriz de Gram de un producto herm´ıtico sobre V en la base B, que est´a definido por < u, v >= [u]∗B H[v]B .

2.3 El lema de Schur y sus aplicaciones ˜ n×n. Cuando En este apartado todas las matrices que aparecen son cuadradas de tamano hablemos de transformaci´on unitaria (respectivamente, ortogonal) nos referiremos a la transformaci´on lineal asociada a una matriz unitaria (resp. ortogonal). Comenzamos con la siguiente definici´on fundamental. Definici´on 2.3.1 Una matriz A ∈ Kn×n se dice que es semejante a otra matriz B ∈ Kn×n si existe una matriz invertible P ∈ Cn×n tal que P −1 AP = B. La matriz A ∈ Cn×n es unitariamente semejante a la matriz B ∈ Cn×n si existe una matriz unitaria U ∈ Cn×n tal que U ∗ AU = B. Finalmente, una matriz A ∈ Rn×n es ortogonalmente semejante a una matriz B ∈ Rn×n si existe una matriz ortogonal P ∈ Cn×n tal que P T AP = B. Observaci´on 2.3.2 Los tres conceptos introducidos en la definici´on anterior son sim´etricos. Esto significa que, por ejemplo, si una matriz A es semejante a otra matriz B entonces B es a su vez semejante a la matriz A. La justificaci´on de este hecho es la siguiente. Si P¢−1 AP = ¡ −1 B, para alguna matriz invertible P , entonces P BP −1 = A. Como P = P −1 , B es semejante a A, donde la matriz que da la semejanza, seg´un la definici´on, es P −1 . En realidad, obtendr´ıamos la misma definici´on de semejanza si escribi´eramos P AP −1 = B en lugar de P −1 AP = B. Geom´etricamente, dos matrices semejantes representan la misma transformaci´on lineal pero en dos bases distintas. Es decir, si A es semejante a B entonces existe una transformaci´on lineal T de Kn y dos bases de Kn , B y B 0 de manera que A = M (T, B, B)

y

B = M (T, B0 , B 0 )

2.3 El lema de Schur y sus aplicaciones

37

(donde M (T, B1 , B2 ) representa la matriz de la transformaci´on lineal T en las bases B1 del espacio de partida y B2 del espacio de llegada). En este caso, la matriz P que da la semejanza es la matriz de cambio de base de B0 a B. No obstante, aunque dos matrices semejantes representen la misma transformaci´on puede ocurrir que al cambiar de base el producto escalar de dos vectores se vea alterado. Esto, en cambio, no ocurre con la semejanza unitaria/ortogonal, que es una semejanza mediante matrices unitarias/ortogonales. Por tanto, las matrices de cambio de base en la semejanza unitaria/ortogonal conservan el producto escalar. Puesto que, como acabamos de ver, una transformaci´on lineal puede expresarse en diferentes (en realidad infinitas) bases a trav´es de matrices semejantes, es natural preguntarse por la expresi´on m´as sencilla posible de dicha transformaci´on. Las matrices m´as elementales, desde el punto de vista operativo, son las matrices diagonales. Esto nos lleva a la siguiente definici´on. Definici´on 2.3.3 Una matriz es diagonalizable si es semejante a una matriz diagonal. Una matriz es unitariamente diagonalizable (resp. ortogonalmente diagonalizable) si es unitariamente (resp. ortogonalmente) semejante a una matriz diagonal. En otras palabras, las matrices diagonalizables son aquellas que, mediante un cambio de base (el mismo en el espacio de partida y en el de llegada) pueden transformarse en una matriz diagonal. La cuesti´on que surge de forma inmediata es: ¿Todas las matrices son diagonalizables? Y la respuesta, ya conocida, es: NO. El siguiente resultado, que es un mero recordatorio del curso ´ de Algebra elemental, nos da una caracterizaci´on de las matrices diagonalizables. Teorema 2.3.4 Sea A ∈ Cn×n . Los siguientes enunciados son equivalentes. i) A es diagonalizable. ii) Existe una base de Cn formada por vectores propios de A. iii) La multiplicidad geom´etrica de cada autovalor de A coincide con su multiplicidad algebraica. Nuestro inter´es se centra ahora en resolver la misma pregunta para la diagonalizaci´on unitaria y ortogonal, es decir: ¿Todas las matrices son unitariamente/ortogonalmente diagonalizables? En caso negativo, ¿qu´e matrices lo son y cu´ales no? La respuesta a esta pregunta es una de las aplicaciones del lema de Schur, que enunciamos a continuaci´on. Teorema 2.3.5 (Lema de Schur) Dada una matriz A ∈ Cn×n , existen una matriz U unitaria y una matriz S triangular superior de manera que U ∗ AU = S.

(2.1)

Dicho de otra forma: toda matriz es unitariamente semejante a una matriz triangular superior. M´as a´un, si la matriz de partida A es real con todos sus autovalores reales, entonces la matriz U puede tomarse ortogonal.

2. ESTUDIO ESPECTRAL DE ENDOMORFISMOS

38

La matriz S en (2.1) se conoce como forma de Schur de A, y, puesto que es triangular y semejante a la matriz A, sus entradas diagonales son los autovalores de A. En cambio, las columnas de la matriz U en (2.1) no son necesariamente los autovectores de A (s´olo podemos asegurar que la primera lo es). Como hemos visto en el Lema de Schur, si la matriz de partida es real con autovalores reales, entonces la matriz de paso, U , puede tomarse ortogonal. La condici´on de que los autovalores sean reales es necesaria para que esto ocurra. En el caso, m´as general, en que los autovalores son complejos, tenemos el siguiente resultado. Teorema 2.3.6 (Lema de Schur real) Dada una matriz A ∈ Rn×n existe una matriz ortogonal tal que   A1 ∗ . . . ∗  ..   0 A2 .  T ,  (2.2) Q AQ =  .  ..  .. . ∗  0 . . . 0 Ak donde Ai , para i = 1, . . . , k, es una matriz 1 × 1 o´ 2 × 2. En este u´ ltimo caso, Ai tiene dos valores propios conjugados. La matriz derecha de la f´ormula (2.2) es una matriz triangular superior por bloques y se conoce como la forma de Schur real deA. Los s´ımbolos ∗ en dicha matriz representan bloques que no son relevantes. La siguiente definici´on es clave para responder a la pregunta sobre qu´e matrices son unitariamente diagonalizables. Definici´on 2.3.7 Una matriz A ∈ Cn×n es normal si AA∗ = A∗ A. Ejercicio 2.3.8 Demostrar que si una matriz es normal y triangular superior entonces es diagonal. El conjunto de las matrices normales incluye a (casi) todos los tipos de matrices de la definici´on 2.1.2. En el siguiente cuadro mostramos algunos de estos tipos de matrices.

Matrices normales:

Sim´etricas reales Herm´ıticas Unitarias Ortogonales (reales)

N´otese que, en cambio, las matrices sim´etricas complejas no son necesariamente normales. En el siguiente ejercicio se pide encontrar un contraejemplo. Ejercicio 2.3.9 Encontrar un ejemplo de matriz sim´etrica (compleja) que no sea normal. Ahora ya podemos enunciar la caracterizaci´on del conjunto de matrices que son unitariamente diagonalizables. Se trata de una caracterizaci´on, como se ver´a, muy elemental en su enunciado.

2.3 El lema de Schur y sus aplicaciones

39

Teorema 2.3.10 Una matriz es unitariamente diagonalizable si y s´olo si es normal. Demostraci´on. Demostraremos primero que una matriz unitariamente diagonalizable es normal. Sea A, por tanto, una matriz unitariamente diagonalizable y U ∗ AU = D una diagonalizaci´on unitaria de A (U es, por tanto, una matriz unitaria y D una matriz diagonal). Entonces U ∗ A∗ U = D∗ y, multiplicando ambas identidades tenemos (U ∗ AU )(U ∗ A∗ U ) = DD∗ = D∗ D = (U ∗ A∗ U )(U ∗ AU ), donde hemos usado, para obtener la segunda igualdad, que dos matrices diagonales siempre conmutan. Como U ∗ U = U U ∗ = I (pues U es unitaria), la anterior igualdad es equivalente a U ∗ AA∗ U = U ∗ A∗ AU. Finalmente, multiplicando ambos t´erminos a la izquierda por U y a la derecha por U ∗ llegamos a la igualdad AA∗ = A∗ A, es decir, A es normal. Supongamos ahora que A es una matriz normal. Por el lema de Schur, existen una matriz U unitaria y una matriz S triangular superior tales que U ∗ AU = S, luego U ∗ A∗ U = S ∗ y, multiplicando ambas expresiones tenemos que U ∗ AA∗ U = SS ∗ . Como AA∗ = A∗ A, lo anterior es igual a U ∗ A∗ AU = U ∗ A∗ U U ∗ AU = S ∗ S, de modo que SS ∗ = S ∗ S, es decir, S es normal. Como S es triangular superior, el Ejercicio 2.3.8 nos permite concluir que S es diagonal, y, por tanto, A es unitariamente diagonalizable. Ejercicio 2.3.11 Demostrar que los autovalores de una matriz real sim´etrica son n´umeros reales. El teorema 2.3.10 caracteriza a las matrices unitarias como aquellas que tienen n autovectores ortonormales con el producto escalar usual herm´ıtico. Cabe preguntarse qu´e ocurre en el caso real, es decir, ¿cuales son las matrices reales n × n que tienen un total de n autovectores ortonormales con el producto escalar usual? La respuesta la daremos en el Teorema 2.3.15 y es muy sencilla de formular: se trata de las matrices real sim´etricas. Este resultado ser´ıa muy sor´ prendente si no lo conoci´eramos ya del curso elemental de Algebra lineal. El Teorema 2.3.15 nos da, por tanto, una caracterizaci´on de una de las clases importantes de matrices normales que estamos tratando en este curso (las matrices reales sim´etricas). Este resultado forma parte de una serie de teoremas que veremos a continuaci´on, en los que obtendremos caracterizaciones

40

2. ESTUDIO ESPECTRAL DE ENDOMORFISMOS

de otras clases importantes de matrices normales, concretamente de las matrices herm´ıticas, unitarias y ortogonales. En estas caracterizaciones, a la propiedad fundamental de las matrices normales (ser unitariamente diagonalizables), se a˜naden ciertas propiedades espectrales (es decir, relativas a los autovalores y autovectores). Vamos a considerar por separado el caso complejo y el caso real. Es importante destacar que todas estas caracterizaciones se derivar´an del Lema de Schur, lo que da una idea de la importancia de este resultado. Matrices complejas Los dos resultados que enunciamos a continuaci´on caracterizan a los conjuntos de matrices herm´ıticas y unitarias. Concretamente, el primero de ellos establece que las matrices herm´ıticas son aquellas que se pueden diagonalizar unitariamente y tienen todos sus autovalores reales. El segundo establece que las matrices unitarias son aquellas que pueden diagonalizarse unitariamente y tiene todos sus autovalores de norma 1. Teorema 2.3.12 A ∈ Cn×n tiene n autovectores ortogonales y todos sus autovalores reales si y s´olo si A es herm´ıtica. Demostraci´on. En primer lugar, observamos que una matriz A ∈ Cn×n tiene n autovectores ortogonales si y solo si tiene n autovectores ortonormales, pues los autovectores pueden normalizarse y siguen siendo autovectores (recu´erdese que cualquier m´ultiplo de un autovector sigue siendo autovector). Por tanto, la condici´on “A tiene n autovectores ortogonales” es equivalente a la condici´on “A es unitariamente diagonalizable”. Para demostrar el Teorema, supongamos primero que A es herm´ıtica, es decir, A = A∗ . Por el Lema de Schur, U ∗ AU = S, donde U es unitaria y S es triangular superior. Como A∗ = A, aplicando la traspuesta conjugada a la igualdad anterior tenemos que U ∗ AU = S ∗ , de donde se concluye que S = S ∗ . Pero como S es triangular superior, S ∗ es triangular inferior, luego S es, al mismo tiempo, triangular superior e inferior, es decir, es diagonal. Adem´as, la igualdad S = S ∗ implica que los autovalores de S (entradas diagonales) son n´umeros reales. Como los autovalores de S y los de A coinciden, se concluye que los autovalores de A son reales. Para demostrar el rec´ıproco, supongamos que A es unitariamente diagonalizable y con todos sus autovalores reales. Entonces, existen una matriz unitaria, U , y una matriz diagonal con entradas reales, D, tales que U ∗ AU = D. Tomando la traspuesta conjugada en esta igualdad llegamos a que U ∗ A∗ U = D, de donde, igualando ambas expresiones, U ∗ AU = U ∗ A∗ U . Ahora, basta multiplicar a la izquierda por U y a la derecha por U ∗ para concluir que A = A∗ , es decir, A es herm´ıtica. Teorema 2.3.13 A ∈ Cn×n tiene n autovectores ortogonales y todos sus autovalores de m´odulo 1 si y s´olo si A es unitaria. Demostraci´on. Supongamos primero que A es unitaria, es decir, A∗ A = In = AA∗ . Por el Lema de Schur, U ∗ AU = S, donde U es unitaria y S es triangular superior. Tomando la traspuesta conjugada en esta igualdad llegamos a que U ∗ A∗ U = S ∗ y, multiplicando ambas expresiones, nos queda In = SS ∗ . Por tanto, S es tambi´en unitaria y, en particular, es normal. Por el Ejercicio 2.3.8 se concluye que S es diagonal. Adem´as, la igualdad S ∗ S = In implica, igualando las entradas diagonales, que los autovalores de S tienen norma (m´odulo) 1. Como los autovalores de A y S coinciden, tambi´en los autovalores de A tienen m´odulo 1.

2.3 El lema de Schur y sus aplicaciones

41

Rec´ıprocamente, supongamos que A es unitariamente diagonalizable y todos sus autovalores tiene m´odulo 1. Entonces, existen una matriz U unitaria y una matriz diagonal D con entradas diagonales de m´odulo 1 tales que U ∗ AU = D. Tomando la traspuesta conjugada tenemos que U ∗ A∗ U = D∗ y, multiplicando ambas expresiones, llegamos a que U ∗ AA∗ U = DD∗ . Como los autovalores de D (entradas diagonales) tienen m´odulo 1, DD∗ = In , de modo que U ∗ AA∗ U = In y, multplicando a la derecha por U y a la izquierda por U ∗ llegamos a que AA∗ = In , es decir, A es unitaria. Problema Resuelto 2.3.14 Estudiar si es posible diagonalizar unitariamente la matriz 

 0 0 1 A= 1 0 0  0 1 0 y, en caso afirmativo, hallar una diagonalizaci´on unitaria. Soluci´on. Puede comprobarse que AAT = AT A, lo que significa (A es real) que A es normal y, por tanto, es unitariamente diagonalizable. Los autovalores de A se calculan a partir del polinomio caracter´ıstico, y son λ1 = 1,

λ2 = e2πi/3 ,

λ3 = e4πi/3 .

Los autovectores asociados (normalizados) son √  1/√3 v1 =  1/√3  , 1/ 3 

√  1/ 3√ v2 =  e−2πi/3 /√ 3  , e2πi/3 / 3 

√  1/ 3√ v3 =  e−4πi/3 /√ 3  , e4πi/3 / 3 

£ ¤ que, como puede comprobarse, son, efectivamente, ortogonales. As´ı pues, si U = v1 v2 v3 y D = diag(λ1 , λ2 , λ3 ), entonces U ∗ AU = D es una diagonalizaci´on unitaria de A. La matriz A es ortogonal, pues AAT = AT A = I3 . Seg´un el Teorema 2.3.13 los autovalores de A tienen m´odulo 1, que es algo que podemos comprobar inmediatamente. La matriz A del Problema Resuelto anterior es una matriz real. En cambio, para diagonalizarla hemos utilizado n´umeros complejos, que son inevitables porque dos de los autovalores de A son complejos (no reales). Podr´ıamos plantearnos la posibilidad de hallar una forma can´onica de A que, aunque no fuera diagonal (ya que esto es imposible porque, como hemos dicho, dos de los autovalores no son reales) fuera “lo m´as parecido” a una matriz diagonal y que no involucrase valores complejos. Es posible hallar esta matriz y eso lo veremos en el Teorema 2.3.17. Matrices reales El primer resultado nos proporciona la caracterizaci´on que ya hemos aludido de las matrices reales ortogonalmente diagonalizables. Teorema 2.3.15 Una matriz real es ortogonalmente diagonalizable si y s´olo si es sim´etrica.

42

2. ESTUDIO ESPECTRAL DE ENDOMORFISMOS

Demostraci´on. Si A es una matriz real que es ortogonalmente diagonalizable, entonces P T AP = D, donde P es ortogonal y D es diagonal. Trasponiendo la igualdad anterior tenemos que P T AT P = DT = D, lo que implica que P T AP = P T AT P y, multiplicando a la izquierda por P y a la derecha por P T llegamos a que A = AT , es decir, A es sim´etrica. A la inversa, si A es una matriz real sim´etrica, por el ejercicio 2.3.11 A tiene autovalores reales y, por el Teorema 2.3.5 la matriz U en (2.1) puede tomarse real ortogonal. Tenemos, por tanto, U T AU = S, donde S es tambi´en real porque es el producto de matrices reales. Trasponiendo esta identidad llegamos a que U T AU = S T , e igualando, se tiene que S es sim´etrica. En consecuencia, S es normal y, por el Ejercicio 2.3.8, conclu´ımos que S es diagonal. Como la matriz U en (2.1) es ortogonal, A es ortogonalmente diagonalizable. En el segundo resultado obtenemos una forma can´onica real para las matrices reales normales. Teorema 2.3.16 Una matriz A ∈ Rn×n es normal si y s´olo si A es ortogonalmente semejante a una matriz diagonal por bloques de tama˜no 1 × 1 y 2 × 2   λ1   ..   .     λ p  , (2.3)   A 1     ..   . As donde los bloques diagonales de tama˜no 2 × 2 son de la forma · ¸ aj bj Aj = , j = 1, . . . , s, −bj aj y corresponden a un par de autovalores complejos conjugados de A, aj ± ibj . Finalmente, el u´ ltimo resultado de esta secci´on, nos da una forma can´onica real para las matrices ortogonales. Teorema 2.3.17 Una matriz A ∈ Rn×n es ortogonal si y s´olo si es ortogonalmente semejante a una matriz diagonal por bloques de tama˜no 1 × 1 y 2 × 2   λ1   ..   .     λp ,  (2.4)   A1     ..   . As

2.3 El lema de Schur y sus aplicaciones

43

con λi = ±1, para i = 1 . . . , p, y los bloques diagonales de tama˜no 2 × 2 son de la forma · ¸ cos θj sen θj Aj = , j = 1, . . . , s, − sen θj cos θj y corresponden a un par de autovalores complejos conjugados de A (de norma 1), cos θj ± i sen θj . La matriz (2.3) o (2.4) se conoce como la diagonalizaci´on ortogonal de A. Hay que tener en cuenta que no se trata de una matriz diagonal en sentido estricto, sino de una matriz diagonal por bloques. Problema Resuelto 2.3.18 Hallar una diagonalizaci´on ortogonal de la matriz A del Problema Resuelto 2.3.14. Soluci´on. Para hallar la diagonalizaci´on ortogonal (real) de A, conviene recordar c´omo se obtiene la diagonalizaci´on real de una matriz real a partir de la diagonalizaci´on compleja: la matriz diagonal por bloques, Dr , contiene los autovalores reales (bloques 1×1) y unos bloques diagonales 2 × 2 de la forma · ¸ a b , −b a donde a, b son, respectivamente, las partes real e imaginaria de cada uno de los pares de autovalores conjugados (es decir, λ = a + bi, λ = a − bi son dos autovalores conjugados de A). Por otro lado, la matriz de paso, P , se obtiene tomando las partes real e imaginaria de los autovectores (y orden´andolas de acuerdo a la ordenaci´on de los autovalores en Dr ). Como A tiene un autovalor real, λ1 = 1, y dos autovalores complejos conjugados, λ2 = 2πi/3 e = cos(2π/3)+i sen(2π/3), λ3 = e4πi/3 = cos(2π/3)−i sen(2π/3), la matriz diagonal ser´a     1 0 1 0 0 √0 3/2  . Dr =  0 cos(2π/3) sen(2π/3)  =  0 −1/2 √ 0 − sen(2π/3) cos(2π/3) 0 − 3/2 −1/2 Por otro lado, la matriz de paso, P , se obtiene a partir de los autovectores v1 , v2 , v3 de la siguiente manera (n´otese que v3 = v 2 y no lo usaremos). Primero calculamos los vectores √  √    3 1/√3 1/ √ 1 u1 = v1 =  1/√3  , u2 = (v2 + v 2 ) = Re(v2 ) =  −1/(2√3)  , 2 1/ 3 −1/(2 3)   0 i u3 = − (v2 − v 2 ) = Im(v2 ) =  −1/2  , 2 1/2 que son ortogonales pero no unitarios, as´ı que los normalizamos para obtener √   √ 1/√3 2/ √6 0√ P =  1/√3 −1/√6 −1/√ 2  . 1/ 3 −1/ 6 1/ 2

2. ESTUDIO ESPECTRAL DE ENDOMORFISMOS

44

Puede comprobarse que P T AP = Dr . El siguiente cuadro-resumen contiene los resultados de los teoremas 2.3.10, 2.3.12, 2.3.13 y 2.3.15. Conjunto de matrices Normales Herm´ıticas Unitarias Reales sim´etricas

Caracterizaci´on Unitariamente diagonalizables Unitariamente diagonalizables y autovalores reales Unitariamente diagonalizables y autovalores de norma 1 Ortogonalmente diagonalizables

´ geometrica ´ 2.4 Interpretacion de las matrices ortogonales en R2 y R3 Hemos visto en el Teorema 2.3.17 que toda matriz ortogonal es semejante (ortogonalmente) a una matriz diagonal por bloques de tama˜no 1 × 1 y 2 × 2, y que estos bloques son de la forma · ¸ cos θ − sen θ 1, −1, (2.5) sen θ cos θ asociados, respectivamente, con los autovalores 1, −1 y el par conjugado cos θ+i sen θ, cos θ− i sen θ. Observaci´on 2.4.1 En la secci´on anterior hemos escrito los bloques asociados al par de au£ ¤ θ sen θ tovalores conjugados cos θ + i sen θ, cos θ − i sen θ en la forma −cos para obtener los sen θ cos θ autovectores de la forma real tomando directamente la parte real y la parte imaginaria de los autovectores de la forma compleja. Ahora, para favorecer la interpretaci´ geom´ £ cos θ o−nsen ¤ etrica de las θ transformaciones asociadas, nos interesa expresarlos en la forma sen θ cos θ , con el signo (−) en la entrada (1, 2) en lugar de la entrada (2, 1). N´otese que una matriz se obtiene de la otra cambiando θ por −θ y viceversa. Estamos interesados en aquellas transformaciones lineales del plano R2 o del espacio R3 que conservan la distancia entre vectores. Si extendemos este concepto a un espacio vectorial eucl´ıdeo/unitario cualquiera, tenemos la siguiente definici´on. Definici´on 2.4.2 Una isometr´ıa de un espacio vectorial euclideo/unitario V es una transformaci´on lineal, T , de V tal que d(T (u), T (v)) = d(u, v), para cualesquiera dos vectores u, v ∈ V. Es decir, una isometr´ıa es un transformaci´on lineal que preserva la distancia entre vectores. Como ya sabemos, las transformaciones lineales asociadas a las matrices ortogonales preservan la distancia asociada al producto escalar usual. Son, por tanto, isometr´ıas. En realidad, no es dif´ıcil demostrar que estas son las u´ nicas isometr´ıas de Rn .

2.4 Interpretaci´on geom´etrica de las matrices ortogonales en R2 y R3

45

Teorema 2.4.3 Toda isometr´ıa en Rn , con el producto escalar usual, es una transfomaci´on lineal asociada a una matriz ortogonal. En este apartado nos centramos en las isometr´ıas de los espacios vectoriales R2 y R3 con el producto escalar usual. Vamos a ver que todas las isometr´ıas del plano R2 pueden clasificarse en tres grupos, mientras que las del espacio R3 pueden clasificarse en cuatro grupos, que son los indicados en la siguiente tabla.

Isometr´ıas de R2

  I Identidad I Giro de a´ ngulo θ  I Simetr´ıa respecto a un eje

Isometr´ıas de R3

 I    I I    I

Identidad Giro de eje r y a´ ngulo θ Simetr´ıa respecto a un plano Giro-reflexi´on respecto a un eje r y el plano r⊥

Tabla 2.4: Isometr´ıas del plano y el espacio

Antes de comenzar con la clasificaci´on, vamos a explicar brevemente aquellas isometr´ıas de la tabla anterior que requieren alguna explicaci´on. En concreto, el giro en R3 y el giroreflexi´on (tambi´en en R3 ). En lo sucesivo, entenderemos que un giro en un plano se efect´ua en sentido positivo si lo hace en sentido opuesto al movimiento de las agujas del reloj. El giro de eje r y a´ ngulo θ en R3 es un giro en el plano perpedicular al eje r. Es decir, para conocer la imagen de un vector v cualquiera, consideramos el extremo de dicho vector como un punto P y el plano perpendicular a r por dicho punto. Tomando como centro de giro el punto de corte de r con el plano perpendicular, giramos el punto P un a´ ngulo de θ radianes en sentido positivo y el punto que obtenemos es el extremo de la imagen del vector v. El giro-reflexi´on de eje r y plano perpendicular a r es la composici´on de un giro de eje r (con un a´ ngulo θ, como se ha descrito en el p´arrafo anterior) y una simetr´ıa (reflexi´on) respecto al plano perpendicular a r. La composici´on puede hacerse en los dos sentidos, es decir, primero el giro y luego la simetr´ıa o a la inversa, porque ambos movimientos conmutan. El procedimiento que seguiremos para clasificar las isometr´ıas de R2 y R3 se basa en los siguientes dos argumentos: (1) Cada uno de los tipos de isometr´ıas descritos en la Tabla 2.4 es invariante por semejanza ortogonal. (2) Todas las posibles formas can´onicas de matrices ortogonales descritas en el Teorema 2.3.17 corresponden a alguna de las isometr´ıas de la Tabla 2.4. La justificaci´on del punto (1) es la siguiente: recordemos que una semejanza consiste simplemente en un cambio de base. Si, adem´as, la semejanza es ortogonal, dicho cambio de base preserva el producto escalar y, en particular, las distancias. Por tanto, una semejanza ortogonal representa la misma isometr´ıa pero expresada en bases distintas.

2. ESTUDIO ESPECTRAL DE ENDOMORFISMOS

46

Para demostrar (2) analizaremos todas las posibles formas can´onicas que se derivan del Teorema 2.3.17. Trataremos por separado las del plano R2 y las del espacio R3 . Isometr´ıas b´asicas de R2 Como ya hemos visto en el Teorema 2.3.17, la forma can´onica de una isometr´ıa de R2 consiste en una matriz 2 × 2 que es diagonal por bloques, y los bloques diagonales son de la forma (2.5). Las posibilidades se reducen a · ¸ 1 0 (i) Autovalor 1 (doble) −→ I2 = · ¸ 0· 1 ¸ 1 0 −1 0 (ii) Autovalor 1 y autovalor −1 −→ , 0 −1 · 0 1¸ −1 0 (iii) Autovalor −1 (doble) −→ −I2 = 0 −1 ¸ · cos θ − sen θ (iv) Dos autovalores complejos conjugados −→ sen θ cos θ La matriz de (i) corresponde a la transformaci´on identidad en cualquier base (por tanto, este caso se dar´a s´olo cuando la matriz de partida sea la propia identidad). Las matrices del caso (ii) corresponden, respectivamente, a una simetr´ıa respecto del eje Y (vertical) y una simetr´ıa respecto del eje X (horizontal). Si llegamos a alguna de estas dos formas can´onicas es porque la transformaci´on de partida corresponde a una simetr´ıa. La matriz del caso (iii) corresponde a un giro de 180o (π radianes) y, como en el caso (i), s´olo se dar´a este caso si la propia matriz de partida es −I2 . La matriz del caso (iv) corresponde a un giro de a´ ngulo θ (en sentido positivo). Por lo tanto, el caso (iii), considerado como isometr´ıa, es un caso particular del caso (iv) (con θ = π). Como se ve, una isometr´ıa puede identificarse a trav´es de los autovalores de la matriz ortogonal de partida. Pero los autovalores por s´ı solos no determinan completamente la isometr´ıa, pues, por ejemplo, en el caso de una simetr´ıa axial, no dicen cu´al es el eje. No obstante, a partir de ellos puede obtenerse el eje de simetr´ıa: se trata de la recta generada por el autovector asociado al autovalor 1 (que es el autovector que queda fijo por la transformaci´on). En la siguiente tabla conclusiva resumimos todo lo que acabamos de estudiar. Esta tabla contiene la informaci´on fundamental que se necesita para poder describir la isometr´ıa que representa una matriz ortogonal dada A ∈ R2×2 . Autovalores

Tipo de transformaci´on

1, 1

Identidad

1, −1

Simetr´ıa de eje r

cos θ + i sen θ, cos θ − i sen θ

Giro

Forma · can´o¸nica 1 0 · 0 1 ¸ 1 0 0 −1 · ¸ cos θ − sen θ sen θ cos θ

Elementos − r ≡ ker(A − I) ´ Angulo: θ

2.4 Interpretaci´on geom´etrica de las matrices ortogonales en R2 y R3

47

Isometr´ıas b´asicas de R3 Razonando como en el caso de R2 , las posibles formas can´onicas ahora son: 

(i)

Autovalor 1 (triple)

(ii)

Autovalor 1 (doble) y autovalor −1 (simple) −→

(ii)

Autovalor 1 (simple) y autovalor −1 (doble) −→

(iv) Autovalor −1 (triple)

−→

−→

(v)

Autovalor 1 (simple) y dos autovalores complejos conjugados

−→

(vi)

Autovalor −1 (simple) y dos autovalores complejos conjugados

−→

1 I3 =  0 0  1 0  0 1  0 0 1 0  0 −1 0 0  −1 −I3 =  0 0  1 0  0 cos θ  0 sen θ −1 0  0 cos θ 0 sen θ

 0 0 1 0  0 1 0 0  −1  0 0  −1  0 0 −1 0  0 −1 0 − sen θ  cos θ  0 − sen θ  cos θ

La matriz del caso (i) corresponde a la transformaci´on identidad y, como en el caso de R2 , s´olo se dar´a este caso si la propia matriz de partida es la identidad. La matriz del caso (ii) corresponde a una simetr´ıa respecto del plano z = 0 y es una de las tres posibles matrices can´onicas con autovalor 1 doble y autovalor −1 simple. Las dos restantes se obtienen colocando el autovalor −1, respectivamente, en la entrada (1, 1) y en la entrada (2, 2), y correponden, respectivamente, a una simetr´ıa respecto al plano x = 0 y una simetr´ıa respecto al plano y = 0. Cualquier simetr´ıa plana se reduce a una de estas tres formas can´onicas. La matriz del caso (iii) corresponde a una simetr´ıa respecto de la recta y = z = 0 (eje X), pero, como isometr´ıa, es un caso particular del caso (v), tomando θ = π. Dicha matriz es, como en el caso anterior, una de las tres posibles matrices can´onicas con autovalor 1 simple y autovalor −1 doble. Las dos restantes se obtienen situando el autovalor 1, respectivamente, en la entrada (2, 2) y en la entrada (3, 3), y correponden, respectivamente, a una simetr´ıa respecto al eje x = z = 0 (eje Y ) y una simetr´ıa respecto al eje x = y = 0 (eje Z). La matriz del caso (iv) corresponde a una simetr´ıa respecto del origen de coordenadas. Es un caso particular del tipo (vi), tomando θ = π. La matriz del tipo (v) corresponde a un giro de eje y = z = 0 (eje X) y a´ ngulo θ o, lo que es igual, un giro de a´ ngulo θ en el plano x = 0. La matriz de cualquier giro respecto a un eje en R3 se reduce por semejanza ortogonal a una matriz de este tipo. Finalmente, la matriz del caso (vi) corresponde a un giro-reflexi´on de eje y = z = 0 y plano perpendicular x = 0. Igual que en el caso del plano R2 , para clasificar las isometr´ıas basta conocer los autovalores de la matriz de partida. Esta informaci´on, al igual que en el caso anterior, no es suficiente para deteminar por completo la isometr´ıa. Para eso necesitamos conocer, en algunos casos, los autovectores asociados. Por ejemplo, en el caso (ii) el plano de simetr´ıa es el autoespacio

2. ESTUDIO ESPECTRAL DE ENDOMORFISMOS

48

asociado al autovalor 1, que corresponde al subespacio fijo por la isometr´ıa. En el caso (v) el eje de giro corresponde igualmente al autoespacio asociado al autovalor 1 (que ahora es una recta), porque es, igualmente, el subespacio que queda fijo por la isometr´ıa. Finalmente, en el caso (vi), en el que no hay vectores fijos, el eje de giro es el autoespacio asociado al autovalor −1 y el plano de reflexi´on es el plano ortogonal a este eje. El cuadro siguiente resume todo lo anterior. En e´ l se encuentra la informaci´on fundamental para describir la isometr´ıa que corresponde a una matriz ortogonal A ∈ R3×3 . Autovalores

Tipo de transformaci´on

1, 1, 1

Forma can´onica " #

Identidad "

1, 1, −1 1, cos θ + i sen θ, cos θ − i sen θ −1, cos θ + i sen θ, cos θ − i sen θ

Simetr´ıa plana " Giro respecto al eje r " Giro-reflexi´on respecto al eje r

1 0 0

1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 −1

0 cos θ sen θ

−1 0 0

0 cos θ sen θ

Elementos −

#

0 − sen θ cos θ

Plano : Π ≡ ker(A − I) #

0 − sen θ cos θ

Eje: r ≡ Gen(A − I) ´ Angulo: θ # Eje: r ≡ Gen(A + I) ´ Angulo: θ

Problema Resuelto 2.4.4 En R3 con el producto escalar usual, se consideran las isometr´ıas: S: simetr´ıa respecto del plano de ecuaci´on x = y. G: giro de eje Y y a´ ngulo π/4 (orientaci´on positiva). R = S ◦ G. Se pide: 1. Hallar la matriz de R respecto de la base can´onica de R3 . Estudiar si dicha matriz es ortogonal. 2. Comprobar que R tiene un u´ nico subespacio propio (autoespacio) de dimensi´on 1 y determinar una base, B1 , de dicho subespacio, L. 3. Hallar una base ortonormal, B2 , del subespacio ortogonal de L. 4. Determinar la matriz de R respecto de la base B = (B1 ; B2 ) de R3 . Interpretar geom´etricamente la aplicaci´on R. Soluci´on. 1. Para hallar la matriz de R respecto a la base can´onica necesitamos conocer la imagen por R de la base can´onica de R3 . Hay varias formas de hacerlo. El procedimiento que seguiremos aqu´ı es calcular primero la matriz de G y la matriz de S respecto de la base can´onica y despu´es multiplicarlas (recu´erdese que la matriz de una composici´on de aplicaciones es el producto de las matrices de esas aplicaciones).

2.4 Interpretaci´on geom´etrica de las matrices ortogonales en R2 y R3

49

• Matriz de G en la base can´onica: El vector e2 = [0 1 0]T es el vector director del eje de giro, luego es invariante por G. As´ı pues, G(e2 ) = e2 . Por otra parte, los otros dos vectores can´onicos e1 = [1 0 0]T y e3 = [0 0 1]T est´an en el plano perpendicular al eje de giro y, por tanto, su imagen se obtiene girando los dos vectores can´onicos de ese plano un a´ ngulo de π/4 radianes en sentido positivo. As´ı   √     √  2/2 cos(π/4) − sen(π/4) − 2/2  =  0  , G(e3 ) =  = . 0 0 G(e1 ) =  √ √0 sen(π/4) cos(π/4) 2/2 2/2 

Por tanto



√

2 2

 M (G, CR3 , CR3 ) =  √0

2 2

√  0 − 22  1 √0  . 2 0 2

• Matriz de S en la base can´onica: Recordamos que la matriz de Householder asociada al vector v = [1 − 1 0]T , ortogonal al plano x = y, corresponde, precisamente, a la simetr´ıa respecto a dicho plano. Por lo tanto, 

     1 0 0 1 −1 0 0 1 0 M (S, CR3 , CR3 ) = Hv = I3 − T vv T =  0 1 0 − −1 1 0  =  1 0 0  . v v 0 0 1 0 0 0 0 0 1 2

• Matriz de R en la base can´onica   M (R, CR3 , CR3 ) = M (S, CR3 , CR3 )M (G, CR3 , CR3 ) = 

0

√ 2 √2 2 2

 1 0√  0 −√ 22  . 2 0 2

El apartado xii) del Lema 2.1.4 nos asegura que M (R, CR3 , CR3 ) es ortogonal. Vamos a comprobarlo, no obstante. Un c´alculo elemental permite comprobar que M (R, CR3 , CR3 )T M (R, CR3 , CR3 ) = I3 , lo que significa que M (R, CR3 , CR3 ) es, efectivamente, ortogonal. En otras palabras, R es una isometr´ıa. 2. Para hallar los subespacios propios de R primero calculamos su polinomio caracter´ıstico: Ã ! Ã ! √ √ √ 2 2 2 pT (λ) = − λ3 − λ2 − λ + 1 = −(λ + 1) λ2 − (1 + )λ + 1 . 2 2 2 √

La u´ nica raiz real es λ1 = −1, pues el polinomio λ2 −(1+ 22 )λ+1 no tiene ra´ıces reales. Por tanto, M (R, CR3 , CR3 ) tiene un u´ nico autovalor real, con multiplicidad algebraica 1. El autoespacio asociado a λ1 = −1 est´a definido por   L = ker(M (R, CR3 , CR3 ) + I) = ker 

1 √

2 √2 2 2

 1 0√  1 − 2√2  . 0 1 + 22

2. ESTUDIO ESPECTRAL DE ENDOMORFISMOS

50

Resolviendo el sistema homog´eneo correspondiente llegamos a que   1 L = Gen  −1√  . 1− 2 Si normalizamos el vector anterior llegamos a    1 B1 = p √   5+2 2

 1  −1√  .  1− 2

3. Primero calculamos una base del subespacio ortogonal L⊥ . Por inspecci´on, es f´acil encontrar dos vectores linealmente independientes perpendiculares al vector generador de L, concretamente    √  2 √ 1    2+ 2  ⊥ L = Gen  1 ,  0  , 0 1 Ahora aplicamos el m´etodo de Gram-Schmidt a la anterior base de L⊥ :  1   q1 = 

 

0



 1 −1√  , u2 =  2+2 2 Por tanto,

√ 2 √1 2



 1 1√  −1√  . q2 = 14+8 2 2+2 2

    1   1  √2  1  √1    p −1 B2 =  2  , . √ √   0  14 + 8 2 2 + 2 2 

4. La matriz de R respecto a la base B se puede calcular usando la f´ormula    = 

1 √ 5+2 2 −√ 1 √ 5+2 √ 2 2 √1− √ 5+2 2

√

M (R, B, B) = M (B, CR3 )−1 M (R, CR3 , CR3 )M (B, CR3 ) −1    √1 √1 √ 1 √ √ 1√ 2 2 0 1 0 14+8 2  5+2 2  √ √ 1 1 1     √1 √ √ √ − √  √  √22 0 −√ 22   − 5+2 2 2 14+8 √ 2 √ 2   2 2 1− 2 2 2+2 0 √ √ √ 0 0 √ 2 2 14+8 2



−1  =  0 0

0

√ 2+ 2 √ 4 √ 10−4 2 4

√0 −

√ 10−4 2 4 √ 2+ 2 4

  . 

5+2 2

1 √ 14+8 2 −√ 1 √ 14+8 √ 2 2 √2+2 √ 14+8 2

√

(2.6)

Naturalmente que ni la matriz inversa M (B, CR3 )−1 ni el producto anterior se han calculado a mano (eso requerir´ıa unas grandes dosis de valor y desesperaci´on a partes iguales).

    

2.4 Interpretaci´on geom´etrica de las matrices ortogonales en R2 y R3

51

De hecho, ni siquiera se han calculado porque no es necesario. Para obtener la matriz de la transformaci´on R en la base B s´olo necesitamos conocer los autovalores conjugados de la matriz de R en la base can´onica. Estos autovalores son p p √ √ √ √ 2+ 2 10 − 4 2 2+ 2 10 − 4 2 λ2 = +i y λ3 = λ2 = −i . 4 4 4 4 La teor´ıa que acabamos de estudiar nos asegura que R es un giro-reflexi´on con respecto √ al eje L y de a´ ngulo θ = arc cos( 2+4 2 ), y que la forma can´onica de R es la matriz (2.6), que es, precisamente, la matriz de R en la base B.

52

2. ESTUDIO ESPECTRAL DE ENDOMORFISMOS

3 VALORES SINGULARES DE MATRICES

´ en valores singulares 3.1 La descomposicion El resultado fundamental de este tema es el siguiente. Teorema 3.1.1 (Descomposici´on en valores singulares) Sea A ∈ Cn×m . Entonces existen: Una matriz unitaria U ∈ Cn×n (ortogonal si A es real), una matriz unitaria V ∈ Cn×n (ortogonal si A es real), y una matriz diagonal Σ ∈ Rn×m cuyas entradas diagonales son n´umeros reales no negativos: Σ(i, i) = σi , que supondremos ordenados en orden no creciente σ1 ≥ σ2 ≥ . . . ≥ σs ≥ 0 (donde s = m´ın{m, n}) de manera que A = U ΣV ∗

(A = U ΣV T si A es real).

(3.1)

Los n´umeros σ1 ≥ σ2 ≥ . . . ≥ σs ≥ 0 est´an un´ıvocamente determinados por A (si se dan en ese orden) y se llaman los valores singulares de A. La descomposici´on (3.1) se conoce como la descomposici´on en valores singulares (SVD) de A. ´ Observaci´on 3.1.2 i) La matriz diagonal Σ de la SVD de una matriz A es unica, pero las ´ matrices U , V no son unicas. ii) Si la matriz A es cuadrada, entonces tenemos U ∗ AV = Σ, siendo Σ una matriz diagonal cuadrada. Esto no es una diagonalizaci´on de A, porque, en general U 6= V . A continuaci´on vamos a estudiar un procedimiento para calcular la SVD de una matriz dada, A. Para ello, necesitaremos el siguiente lema. 53

3. VALORES SINGULARES DE MATRICES

54

Lema 3.1.3 Dada una matriz A ∈ Cn×m cualquiera, los autovalores de A∗ A y A∗ A son n´umeros reales no negativos. Demostraci´on. Las matrices A∗ A y AA∗ son ambas herm´ıticas, y ya sabemos (Teorema 2.3.12) que las matrices herm´ıticas tienen todos sus autovalores reales. Si λ es un autovalor de A∗ A con autovector asociado v, entonces 0 ≤k Av k22 = (Av)∗ Av = v ∗ (A∗ Av) = v ∗ (λv) = λ(v ∗ v) = λ k v k22 . Como k v k22 ≥ 0, necesariamente λ ≥ 0. As´ı pues, los autovalores de A∗ A son n´umeros reales no negativos. De forma an´aloga se procede con los autovalores de la matriz AA∗ . Antes de describir el procedimiento para el c´alculo de la SVD, sintetizaremos las dos propiedades fundamentales en las que descansa dicho procedimiento. La propiedad 1 est´a relacionada con el lema anterior. Seg´un este lema, los autovalores de A∗ A son los cuadrados de ciertos n´umeros reales positivos. Estos n´umeros son, precisamente, los valores singulares. 1 Los cuadrados de los valores singulares no nulos, σ12 ≥ . . . ≥ σr2 (siendo r el mayor ´ındice k tal que σk > 0), son los autovalores de A∗ A. 2 Los autovectores asociados a los autovalores de A∗ A son las columnas de V (ordenadas de acuerdo al orden que tienen los valores singulares en Σ). Se tienen las propiedades an´alogas para la matriz AA∗ . 1’ Los cuadrados de los valores singulares no nulos, σ12 ≥ . . . ≥ σr2 (siendo r el mayor ´ındice k tal que σk > 0), son los autovalores de AA∗ . 2’ Los autovectores asociados a los autovalores de AA∗ son las columnas de U (ordenadas de acuerdo al orden que tienen los valores singulares en Σ). Las columnas de V y U se llaman, respectivamente, vectores singulares derechos y vectores singulares izquierdos de A. Est´an relacionadas por la igualdad Avi = σi ui , 1 ≤ i ≤ s, o por la igualdad equivalente vi =

1 ∗ σi A ui

,

1 ≤ i ≤ s.

Ahora ya estamos en disposici´on de explicar el procedimiento para el c´alculo de la SVD.

3.1 La descomposici´on en valores singulares

55

C´alculo de la SVD (v´ıa derecha) 1. Diagonalizar unitariamente el producto: m−r

z }| { V (A A)V = diag(λ1 , . . . , λr , 0, . . . , 0) ∗

∗

(ordenamos los autovalores: λ1 ≥ . . . ≥ λr > λr+1 = . . . = λm = 0) 2. Los valores singulares de A son σk =    Σ= 

σ1

√ λk , k = 1, . . . , m. 

..

   

. σr 0(n−r)×(m−r)

n×m

3. Para cada uno de los valores singulares no nulos, σi , calcular ui =

1 σi Avi

,

i = 1, . . . , r

4. Para calcular las restantes columnas de U , ur+1 , . . . , um , como sabemos que A∗ U = V ΣT , y las u´ ltimas m − r columnas de ΣT son nulas, se tiene que A∗ ur+1 = . . . = A∗ um = 0, luego ur+1 , . . . , um son una base ortonormal de ker(A∗ ) Por tanto, calculamos una base de ker(A∗ ) y luego la ortonormalizamos.

La condici´on del apartado 4 referente a los vectores ur+1 , . . . , um (es decir, que son una base de ker(A∗ )) se ha inclu´ıdo con el fin de facilitar la b´usqueda de estos vectores, porque el ´ c´alculo de una base de ker(A∗ ) se ha estudiado ya en el curso previo de Algebra Lineal. No obstante, los vectores ur+1 , . . . , um pueden elegirse arbitrariamente con la condici´on de que el conjunto resultante u1 , . . . , um sea una base ortonormal de Cm . Es decir: otra alternativa al paso 4 consiste en completar el conjunto u1 , . . . , ur , obtenido en el paso 3, a una base de Cm y despu´es ortogonalizarla mediante el m´etodo de Gram-Scmidt (s´olo ser´a necesario aplicar el m´etodo de ortonormalizaci´on a los u´ ltimos vectores, porque los obtenidos en el paso 3 ya son ortonormales). El procedimiento an´alogo al anterior, pero comenzando con la matriz AA∗ , es el siguiente.

3. VALORES SINGULARES DE MATRICES

56

C´alculo de la SVD (v´ıa izquierda) 1. Diagonalizar unitariamente el producto:

m−r

z }| { U (AA )U = diag(λ1 , . . . , λr , 0, . . . , 0) ∗

∗

(ordenamos los autovalores: λ1 ≥ . . . ≥ λr > λr+1 = . . . = λm = 0) 2. Los valores singulares de A son σk =    Σ= 

σ1

√ λk , k = 1, . . . , m. 

..

   

. σr 0(n−r)×(m−r)

n×m

3. Para cada uno de los valores singulares no nulos, σi , calcular vi =

1 ∗ σi A u i

,

i = 1, . . . , r

4. Para calcular las restantes columnas de V , vr+1 , . . . , vn , como sabemos que AV = U Σ, y las u´ ltimas n − r columnas de Σ son nulas, se tiene que Avr+1 = . . . = Avn = 0, luego vr+1 , . . . , vn son una base ortonormal de ker(A) Por tanto, calculamos una base de ker(A) y luego la ortonormalizamos.

Naturalmente, cualquiera de los dos procedimientos anteriores (v´ıa derecha y v´ıa izquierda) pueden aplicarse indistintamente a una matriz dada A. No obstante, en algunas ocasiones es m´as ventajoso, a efectos de c´alculo, un procedimiento que otro. El criterio recomendable para elegir la v´ıa derecha o la v´ıa izquierda es el tama˜no de la matriz sobre la que se aplica el primer paso. Es natural, a efectos de c´alculo, elegir aquel procedimiento para el cual dicho tama˜no sea menor. Si elegimos la v´ıa derecha, la matriz a la que aplicamos el paso 1 es A∗ A, cuyo tama˜no es m × m, mientras que si elegimos la v´ıa izquierda trabajaremos sobre la matriz AA∗ , que es de tama˜no n × n. Por tanto, tenemos el siguiente criterio: Si n ≥ m: elegir la v´ıa derecha. Si n ≤ m: elegir la v´ıa izquierda.

3.1 La descomposici´on en valores singulares

57

Problema Resuelto 3.1.4 Calcular la descomposici´on en valores singulares de la matriz   1 −1 A =  −1 1  . 2 −2 Soluci´on. Siguiendo el criterio que acabamos de enunciar, como n = 3 ≥ m = 2, elegiremos la v´ıa derecha. Aplicamos el procedimiento que hemos descrito en la p´agina anterior: 1. Para diagonalizar unitariamente la matriz · ¸ 6 −6 A∗ A = . −6 6 primero calculamos sus autovalores. El polinomio caracter´ıstico es p(λ) = λ(λ − 12), luego los autovalores son λ1 = 12, λ2 = 0. El autoespacio asociado a λ1 = 12 es µ· ∗

ker(A A − 12I2 ) = Gen

1 −1

¸¶ ,

mientras que el autoespacio asociado a λ2 = 0 es µ· ¸¶ 1 ∗ ker(A A) = Gen . 1 Ahora normalizamos la base formada por los dos autovectores para obtener la matriz " 1 # 1 V =

√ 2 − √12

√ 2 √1 2

2. En virtud del paso anterior, la matriz Σ es  √  12 0 Σ= 0 0 . 0 0 3. La primera columna de U (vector singular izquierdo) es el vector  1    √ 1 −1 " √1 # 6 1 1   2 √1  − −1 1  u1 = Av1 = √ = .  6  − √12 σ1 12 2 √ 2 −2 6

3. VALORES SINGULARES DE MATRICES

58

4. La segunda y tercera columnas de U son una base ortonormal de ker(A∗ ). Para hallarla resolvemos el sistema A∗ x = 0, cuya matriz de coeficientes escalonamos f´acilmente · ¸ · ¸ 1 −1 2 1 −1 2 ∼ −1 1 −2 0 0 0 lo que significa que 

   2 1 ∗     0 ker(A ) = Gen , 5  . −1 2 Ahora ortonormalizamos la base anterior de ker(A∗ ). Como los vectores ya son ortogonales, basta con normalizarlos, y el resultado son los dos vectores ortonormales  2   1   u2 =  Finalmente

√ 5

0 − √15 

 U =

   , u3 = 

√1 6 − √16 √2 6

√2 5

0 − √15

√ 30 √5 30 √2 30

√1 30 √5 30 √2 30

 .

  

Puede comprobarse que, efectivamente, se tiene A = U ΣV ∗ . El siguiente resultado caracteriza las columnas de U y V en t´erminos de ciertos espacios invariantes de la matriz de partida A. Teorema 3.1.5 Sea A ∈ Cn×m y A = U ΣV ∗ la SVD de A. Si Σ tiene exactamente r entradas diagonales (valores singulares de A) no nulas, entonces (i) r = rango(A). (ii) Las primeras r columnas de U son una base ortonormal de col(A). (iii) Las primeras r columnas de V son una base ortonormal de col(A∗ ). (iv) Las u´ ltimas n − r columnas de U son una base ortonormal de ker(A∗ ). (v) Las u´ ltimas n − r columnas de V son una base ortonormal de ker(A). Los valores singulares, a diferencia de los autovalores, no son invariantes por semejanza. En cambio, s´ı son invariantes por multiplicaci´on, tanto a derecha como a izquierda, por matrices unitarias. Lema 3.1.6 Si P ∈ Cn×n y Q ∈ Cm×m son matrices unitarias y A ∈ Cn×m entonces los valores singulares de A coinciden con los valores singulares de U AV .

3.1 La descomposici´on en valores singulares

59

Demostraci´on. Si A = U ΣV ∗ es la SVD de A, entonces P AQ = (P U )Σ(Q∗ V )∗ es la SVD de P AQ, porque P U y Q∗ V son unitarias. En ambos casos, la matriz diagonal Σ es la misma. Este lema tiene una aplicaci´on inmediata al c´alculo de autovalores de matrices normales. Dicha aplicaci´on est´a enunciada en el siguiente corolario. Corolario 3.1.7 Si A es una matriz normal entonces los valores singulares de A son los m´odulos de los autovalores de A. El Corolario 3.1.7 nos dice que los valores singulares de una matriz normal A se pueden obtener calculando el m´odulo de los autovalores de A. Este procedimiento es, en general, menos costoso, si se hacen los c´alculos a mano, que el descrito anteriormente para obtener la SVD de una matriz arbitaria. ´ La descomposici´on en valores singulares es una herramienta fundamental en el Algebra Lineal Num´erica, que tiene, adem´as, relevantes aplicaciones en otras a´ reas. En este curso vamos a estudiar algunos de los problemas o conceptos que pueden reinterpretarse mediante el uso de la SVD. Las aplicaciones que aparecer´an en este curso (y que son s´olo una parte del conjunto de aplicaciones de la SVD) est´an resumidas en la siguiente tabla.  I Normas matriciales (Secci´on 6.2)     I C´alculo de la Pseudoinversa (Secci´on 3.2)    I Problemas de m´ınimos cuadrados (Secci´on 3.2) I Perturbaci´on de matrices y sistemas lineales (Secci´on 6.4)      Aproximaciones de rango bajo.   I (Secci´on 3.3) Compresi´on de im´agenes

Aplicaciones de la SVD

Tabla 3.1: Aplicaciones de la SVD

´ en valores singulares reducida La descomposicion La SVD de A (3.1) contiene informaci´on innecesaria cuando n 6= m. Supongamos que n > m (el caso n < m es an´alogo). En este caso, las u´ ltimas n − m filas de la matriz Σ son nulas. Si eliminamos en (3.1) las u´ ltimas n − m columnas de U y las u´ ltimas n − m filas de Σ la matriz resultante al hacer el producto sigue siendo la misma (porque el producto de entradas Σ nulas no aporta nada): z }|  { σ1 V∗   z }| { . ..   ∗  U v   1 }| {¤  £z σm  ..   A = u1 . . . um um+1 . . . un   .   0    ∗   vm ..   . 0

=

£z

b U

u1

z

}| {¤ σ1  . . . um 

b Σ

V

}| ..

∗

{ z }|∗ { v1   ..   . .

. σm

∗ vm

(3.2)

3. VALORES SINGULARES DE MATRICES

60

La expresi´on (3.2) se conoce como la descomposici´on en valores singulares reducida de A, y contiene toda la informaci´on relevante de la SVD (3.1) (que, para distinguirla de la anterior, llamaremos SVD completa). La relaci´on entre las SVD completa y reducida se representa en el siguiente dibujo (para n > m). b Σ b U

=

A

U

Σ

V∗

En el caso n < m la SVD reducida se obtiene de la SVD completa eliminando la u´ ltimas m − n columnas de Σ y las u´ ltimas m − n filas de V ∗ . Ejemplo 3.1.8 La SVD reducida de la matriz A del Problema Resuelto 3.1.4 es z  A=

b U

√1 6 − √16 √2 6

}|

√2 5

0 − √15

{

V∗

b Σ

z z· √ }| ¸{ " 12 0   0 0

√1 2 √1 2

}| #{ 1 √ − 2 √1 2

.

´ al problema de m´ınimos cuadrados. La pseudoin3.2 Aplicacion versa. Volvamos al problema de resolver por m´ınimos cuadrados el sistema lineal (1.5), es decir: Minimizar k Ax − b k2 con x ∈ Cm

(3.3)

donde A ∈ Cn×m , b ∈ Cn y n ≥ m. En lo sucesivo, supondremos que rango(A) = r. Utilizando la SVD de A (3.1) podemos escribir k Ax − b k2 =k (U ΣV ∗ )x − b k2 =k U Σ (V ∗ x) −b k2 | {z } y

=k (U Σ)y − b k2 =k U ∗ (U Σy − b) k2 =k Σy − U ∗ b k2 donde y = V ∗ x y, para obtener la pen´ultima desigualdad, hemos usado la propiedad v) del Lema 2.1.4 aplicada a la matriz unitaria U ∗ . Por tanto, el problema (3.3) es equivalente, con el cambio de variable y = V ∗ x, al problema Minimizar k Σy − U ∗ b k2 con y ∈ Cm . Por simplicidad, denotaremos eb =

h

eb1 . . . ebn

iT

= U ∗ b. Como

k Σy − eb k2 = |σ1 y1 − eb1 |2 + · · · + |σr yr − ebr |2 + |ebr+1 |2 + · · · + |ebn |2

(3.4)

3.2 Aplicaci´on al problema de m´ınimos cuadrados. La pseudoinversa.

61

es una suma de t´erminos no negativos, la norma de Σy − eb se minimiza para el mayor n´umero posible de t´erminos nulos. Es decir, para y1 =

e b1 σ1 , . . . , yr

=

e br σr

.

Deshaciendo el cambio de variable, x = V y, obtenemos la soluci´on al problema (3.3) x1 = V y1 , . . . , xr = V yr . en funci´on de la soluci´on al problema (3.4). Observemos que sobre las entradas xr+1 , . . . , xm no hay ninguna restricci´on, lo que significa que si m > r entonces el problema (3.3) tiene infinitas soluciones. El caso m > r corresponde al caso en que las columnas de la matriz A de coeficientes del sistema (1.5) son linealmente dependientes. Como recordaremos (v´ease la Secci´on 1.5) si las columnas de A son linealmente independientes el problema de m´ınimos cuadrados (3.3) tiene soluci´on u´ nica. Ahora sabemos que el rec´ıproco tambi´en es cierto. En caso de que las columnas de la matriz A no sean linealmente independientes y, por tanto, el problema (3.3) tenga infinitas soluciones, nos interesa encontrar la soluci´on de norma m´ınima. Como y = V ∗ x y la matriz V ∗ es unitaria, la norma de x coincide con la norma de y, luego x es la soluci´on de norma m´ınima de (3.3) si y s´olo si y es la soluci´on de norma m´ınima de (3.4). Claramente dicha soluci´on se obtiene tomando yr+1 = . . . = ym = 0. En el siguiente cuadro se resume el procedimiento descrito en esta secci´on para calcular la soluci´on por m´ınimos cuadrados del sistema (1.5) a trav´es de la SVD de A. Minimizar k Ax − b k2 usando la SV D de A A ∈ Cn×m con n ≥ m r = rango(A) 1. Calcular la SVD de A : A = U ΣV ∗ 2. Calcular : eb = U ∗ b e

3. Calcular : y1 = σb11 , . . . , yr = £ ¤T y0 = y1 . . . ym

e br σr ,

yr+1 = 0, . . . , ym = 0

4. Calcular : x0 = V y0 x0 es la soluci´on de menor norma entre todas las soluciones que minimizan k Ax − b k2 Todas las soluciones que minimizan k Ax − b k2 son de la forma x = x0 + αr+1 vr+1 + · · · + αm vm con αr+1 , . . . , αm ∈ C.

3. VALORES SINGULARES DE MATRICES

62 La pseudoinversa

En este breve apartado continuamos con la notaci´on del apartado anterior. En particular, A ∈ Cn×m y (3.1) es la SVD de A. Si definimos la matriz  1    † Σ = 

σ1

..

   

. 1 σr

0(m−r)×(n−r)

m×n

entonces la soluci´on de norma m´ınima del problema (3.4) se puede expresar como y0 = Σ† U ∗ b , y la soluci´on de norma m´ınima de (3.3) es x0 = V Σ† U ∗ b . Esto nos conduce a la siguiente definici´on. Definici´on 3.2.1 Si A = U ΣV ∗ es la SVD de de A, entonces la matriz A† = V Σ† U ∗ es la pseudoinversa de A. La pseudoinversa se conoce tambi´en con el nombre de inversa generalizada. Su nombre se debe, entre otras razones, a que, si A es invertible, entonces la pseudoinversa A coincide con la inversa de A. Adem´as, la pseudoinversa de A satisface las siguientes igualdades (para cualquier matriz A, no necesariamente invertible) que nos recuerdan a la inversa de A (en particular las dos primeras) AA† A = A . A† AA† = A† . AA† y A† A son herm´ıticas. Las tres propiedades anteriores se conocen como las condiciones de Moore-Penrose y son una definici´on alternativa de la pseudoinversa de A. Es decir: la u´ nica matriz A† que satisface esas tres condiciones es la pseudoinversa de A. Como hemos visto, la pseudoinversa de A nos proporciona la soluci´on por m´ınimos cuadrados del sistema lineal Ax = b. Teorema 3.2.2 Sean A ∈ Cn×m una matriz de rango r, A = U ΣV ∗ la SVD de A y A† la pseudoinversa de A. Entonces a) x0 = A† b minimiza k Ax − b k2 para cualquier x ∈ Cm . b) Entre todos los vectores que minimizan k Ax − b k2 , x0 = A† b es el que tiene menor norma-2.

3.2 Aplicaci´on al problema de m´ınimos cuadrados. La pseudoinversa.

63

c) Cualquier otro vector x e que minimiza k Ax − b k2 es de la forma x e = x0 + v, donde v es una combinaci´on lineal de las u´ ltimas m − r columnas de V . Problema Resuelto 3.2.3 Resolver, por m´ınimos cuadrados, el sistema Ax = donde A es la matriz del Problema Resuelto 3.1.4.

£

1 1 1

¤T

,

Soluci´on. Aplicamos el procedimiento descrito en la p´agina anterior. La SVD de A ya ha sido calculada en el Problema Resuelto 3.1.4. Continuamos por el segundo paso del procedimiento descrito anteriormente:  eb = U ∗ b =  

√2 5

√1 6 − √16 √2 6

0 − √15

de modo que

"

" x0 = V y0 =

√1 2 − √12

2 6

  ,

,

0 √1 2 √1 2

√2 6 √1 5 √8 30

#

√

y0 = y

T    1     1 =  1

√1 30 √5 30 √2 30

#"

#

√

2 6

·

1 6

=

0

− 61

¸ .

x0 es la soluci´on de norma m´ınima. Las restantes soluciones se obtienen a˜nadiendo a e´ sta un m´ultiplo de la segunda columna de V : " 1 # x = x0 + α

√

2 √1 2

.

Problema Resuelto 3.2.4 Resolver el Problema Resuelto 3.2.3 usando la pseudoinversa de A. Soluci´on. Como ya hemos calculado la SVD de A, su pseudoinversa se obtiene de forma inmediata haciendo el c´alculo  1  " 1 #" # √ √1 √2 − 1 6 6 √ √ √1 0 0  √26 2 2 12 √1  0 − A† = V Σ† U ∗ =  1 1 5 5  − √2 √2 0 0 0 √1 √5 √2 30

=

1 12

·

1 −1 2 −1 1 −2

Por tanto x0 =

1 T A b= 12

·

¸ = 1 6

− 61

30

30

1 T A . 12 ¸ .

Como en el Problema Resuelto 3.2.3, todas las soluciones se obtienen a˜nadiendo a x0 un m´ultiplo de la segunda columna de V .

3. VALORES SINGULARES DE MATRICES

64 ´ Calculo de la SVD con Matlab

El comando Matlab para calcular la SVD (completa) de una matriz A es svd(A), mientras que el comando que devuelve la SVD reducida de A es svd(A,0). Para obtener las tres matrices de la descomposici´on, escribiremos [u s v]=svd(A) o [u s v]=svd(A,0) , y el programa nos devolver´a tres matrices: la primera (u) es la matriz U , la segunda (s) es la matriz diagonal Σ y la tercera (v) es la matriz V (en ese orden). Para m´as informaci´on teclear el comando help svd. Ejemplo 3.2.5 Para obtener con Matlab la SVD de la matriz A del Problema Resuelto 3.1.4 primero introducimos la matriz A y despu´es el comando svd(A). )) A=[1 -1;-1 1;2 -2]; )) [u s v]=svd(A) u = 0.4082 0.4082 -0.4082 0.8816 0.8165 0.2367 s = 3.4641 0 0 0 0 0 v = 0.7071 -0.7071 -0.7071 -0.7071

-0.8165 0.2367 0.5266

Obs´ervese que las columnas segunda y tercera de la matriz U que nos devuelve el programa Matlab son distintas de las que hemos obtenido en el Problema Resuelto 3.1.4, lo que demuestra que la matriz U no es u´ nica.

´ de imagenes ´ 3.3 Aproximaciones de rango bajo. Compresion La SVD de A (3.1) puede escribirse como suma de matrices de rango 1 ∗

A = U ΣV =

r X

σi ui vi∗ .

(3.5)

i=1

Cada uno de los sumandos que aparecen en la expresi´on derecha de (3.5), ui vi∗ , es una matriz de rango 1 que se obtiene multiplicando un vector columna por un vector fila (todas las matrices de rango 1 se obtienen de esta manera). En la expresi´on derecha de (3.5) se prescinde de los valores singulares nulos de A y tambi´en de las correspondientes columnas de U y filas de V ∗ , con lo que se reduce al m´ınimo la informaci´on necesaria para “reconstruir” la matriz A. Por otro lado, la descomposici´on (3.5) de A como suma de matrices de rango 1 nos permite resolver el siguiente problema: Problema 3.3: Minimizar k A−B k2 entre todas las matrices B de rango k, con k < r.

3.3 Aproximaciones de rango bajo. Compresi´on de im´agenes

65

Para entender el planteamento exacto del problema es necesario conocer el concepto de norma-2 de una matriz, que se introducir´a en la Secci´on 6.2. Por el momento, es suficiente con saber que buscamos una matriz B de rango inferior al de A que es la ”m´as parecida a la matriz A” en un cierto sentido. De hecho, podemos enunciar el problema como el de “encontrar la matriz B de rango k, con k < r, que mejor se aproxima a la matriz A en sentido de m´ınimos cuadrados”. Este problema tiene diversas aplicaciones, entre ellas la reconstrucci´on de im´agenes, que trataremos a continuaci´on. Pero antes, enunciaremos la soluci´on al problema. Teorema 3.3.1 Sea A ∈ Cn×m de rango r, A = U ΣV ∗ la SVD de A y σ1 ≥ σ2 ≥ . . . ≥ σr > 0 los valores singulares no nulos de A. Sea B otra matriz del mismo tama˜no y de rango k < r. Entonces k A − B k2 ≥ σk+1 . Por tanto, la matriz B de rango k < r que m´as se aproxima a la matriz A en sentido de m´ınimos cuadrados (es decir, la soluci´on al Problema 3.3) es B = U diag(σ1 , . . . , σr , 0, . . . , 0)V ∗ . A continuaci´on vamos a describir qu´e relaci´on tiene lo anterior con la compresi´on y reconstrucci´on de im´agenes. La idea b´asica consiste en identificar una imagen con una matriz. Por simplicidad, s´olo vamos a considerar im´agenes dibujadas en la escala del gris, que supondremos “pixeladas”, o cuadriculadas en una cuadr´ıcula n × m. Estas im´agenes se identifican de manera natural con matrices cuyas entradas son n´umeros entre 0 y 1, donde 0 es el tono m´as oscuro (negro) y 1 es el tono m´as claro (blanco). As´ı pues, la informaci´on de la imagen se cifra en una matriz A de tama˜no n × m. En muchas ocasiones, no es posible, o no es aconsejable almacenar o transmitir toda las entradas de la matriz A, sino comprimirla en una cantidad inferior de entradas que nos permita reconstruir de un modo fiable la imagen original. Esto puede hacerse mediante aproximaciones de A de rango menor, tal y como se describe en el Teorema 3.3.1 y es lo que vamos a explicar a continuaci´on. Sea A = U ΣV T la SVD de A (las matrices U , V son reales porque A tiene entradas P reales). El Teorema 3.3.1 nos dice que Bk = ki=1 σi ui viT es la mejor aproximaci´on de rango k a la matriz A. Para almacenar la matriz Bk necesitamos un total de m · k + n · k = (n + m) · k posiciones de memoria, que corresponden a las n entradas de cada uno de los k vectores ui y a las m entradas de cada uno de los k vectores σi vi . Si k es un valor relativamente peque˜no y m, n son valores grandes, este n´umero es sensiblemente inferior al n´umero total de posiciones de memoria que ocupa la matriz completa, que es n · m. En el siguiente ejemplo comparamos los resultados de almacenar la matriz (imagen) completa con varias aproximaciones de rango bajo. Ejemplo 3.3.2 Consideremos la imagen de la figura 3.3.2 (a). Se trata de una imagen de 320× 200 p´ıxeles en blanco (1) y negro (0). La matriz asociada, A, es, por tanto, una matriz 320×20 cuyas entradas son ceros y unos. El rango de dicha matriz es 200. Las im´agenes (b), (c) y (d) corresponden a las aproximaciones de A por matrices de rangos respectivos 40, 15 y 5, seg´un el enunciado del Teorema 3.3.1.

3. VALORES SINGULARES DE MATRICES

66

(a)

(b)

(c)

(d)

Figura 3.3.2:(a) Imagen original. (b) Aproximaci´on de rango k = 40. (c) Aproximaci´on de rango k = 15 (d) Aproximaci´on de rango k = 5.

En la siguiente tabla comparamos el ahorro en memoria que supone cada una de las anteriores aproximaciones. Este ahorro se indica en la tercera columna, que muestra el cociente (ratio) entre el n´umero de posiciones de memoria de la aproximaci´on y el n´umero de posiciones de la figura original, seg´un la f´ormula (m + n) · k/(m · n) = 520k/6400 ≈ k/123. En la segunda columna aparece el error relativo que se comete al aproximar la matriz de la figura original mediante cada una de las matrices de las restantes figuras. Dicho error es el cociente σk+1 /σ1 , y se entender´a mejor cuando lleguemos a la Secci´on 6.2. k 5 15 40

Error relativo= σk+1 /σ1 0.131 0.051 0.022

Ratio de compresi´on= 520k/64000 0.041 0.122 0.325

Como es natural, el error relativo y el ratio de compresi´on est´an en relaci´on inversa. N´otese que para k = 15 se obtiene una imagen bastante fiable, con un error de un 5 %, mientras que el ahorro en memoria es de aproximadamente el 88 %. La imagen original est´a disponible en Matlab en el archivo de demostraci´on clown.mat. Por otra parte, las cuatro figuras de la portada corresponden a una figura original (la palabra “LUZ”) de 12×26 p´ıxeles y tres aproximaciones de rango bajo. La figura original est´a asociada a una matriz de rango 10, mientras que las aproximaciones (de izquierda a derecha y de arriba abajo) corresponden a matrices de rangos 7, 5 y 3.

4 ´ FORMAS CUADRATICAS

4.1 Formas bilineales En este tema y en el siguiente abandonamos el cuerpo de los n´umeros complejos y nos centramos en el de los n´umeros reales. Hemos visto en la Secci´on 1.1 que un producto escalar (real) es una forma bilineal sim´etrica definida positiva. En la presente secci´on estudiaremos el concepto de forma bilineal, que ser´a fundamental para introducir las formas cuadr´aticas en la Secci´on 4.2 y, como consecuencia, las c´onicas y cu´adricas en el Tema 5. Definici´on 4.1.1 Sea V un espacio sobre R. Una aplicaci´on f : V × V −→ R es una forma bilineal sobre V si f es lineal en cada componente, es decir, si satisface las dos condiciones siguientes i) f (αu + βv, w) = αf (u, w) + βf (v, w) ii) f (u, αv + βw) = αf (u, v) + βf (u, w) para cualesquiera n´umeros reales α, β ∈ R y vectores u, v, w ∈ V. ´ El ejemplo m´as inmediato de forma bilineal es un producto escalar real. Este es un caso particular de forma bilineal y tambi´en un caso particular de las formas bilineales del segundo apartado en el siguiente ejemplo. Ejemplo 4.1.2 Sea V un espacio vectorial real de dimensi´on n. Son formas bilineales (i) Cualquier producto escalar en V. 67

´ 4. FORMAS CUADRATICAS

68

(ii) Dada una matriz A ∈ Rn×n y una base B de V, la aplicaci´on f : V × V −→ R definida por f (x, y) = [x]TB A[y]B , donde x, y son vectores cualesquiera de V. A partir de ahora emplearemos las letras x, y para los vectores de V, en lugar de las letras que venimos empleando u, v, w. Este cambio de notaci´on es m´as acorde con la notaci´on que se encuentra habitualmente en la literatura de formas cuadr´aticas, que se interpretan como aplicaciones, y de c´onicas y cu´adricas (que son ecuaciones). Ejercicio 4.1.3 Demostrar que la aplicaci´on (ii) del Ejercicio 4.1.2 es, efectivamente, una forma bilineal. ´ matricial en una base Expresion A continuaci´on vamos a ver que el apartado (ii) del Ejemplo 4.1.2 es algo m´as que un ejemplo. Concretamente, veremos que cualquier forma bilineal sobre V puede expresarse con respecto a una base B de la forma descrita en dicho ejemplo, para alguna matriz A ∈ Rn×n . Sean, por tanto, B = {v1 , . . . , vn } una base de V, f : V × V −→ R una forma bilineal sobre V y x, y ∈ V dos vectores cualesquiera de V, cuyas coordenadas en la base B son, £ ¤T £ ¤T respectivamente, [x]B = x1 . . . xn e [y]B = y1 . . . yn . Esto significa que u = x1 v1 + . . . + xn vv ,

y = y1 v1 + . . . + yn vv

Aplicando la bilinealidad, f (x, y) = f (

n X

xi vi ,

i=1

n X j=1

yj vj ) =

n X

xi yj f (vi , vj )

i,j=1

que, matricialmente, se expresa    f (v1 , v1 ) . . . f (v1 , vn ) y1 £ ¤T    ..  .. .. T f (x, y) = x1 . . . xn    .  = [x]B M (f, B)[y]B , . . f (vn , v1 ) . . . f (vn , vn ) yn donde



 f (v1 , v1 ) . . . f (v1 , vn )   .. .. M (f, B) =   . . f (vn , v1 ) . . . f (vn , vn )

es la expresi´on matricial de f en la base B. Es importante notar que no existe ninguna restricci´on en la matriz anterior. Por tanto, tenemos el siguiente resultado: Sea V un espacio vectorial real de dimensi´on n y B una base prefijada de V. Toda matriz A ∈ Rn×n es la expresi´on matricial de una forma bilineal sobre V expresada en la base B. Dicha forma bilineal est´a definida por f (x, y) = [x]TB A[y]B .

4.1 Formas bilineales

69

Si desarrollamos el producto indicado en la expresi´on matricial de una forma bilineal en una base obtendremos un polinomio homog´eneo (es decir, con todos los monomios del mismo grado) de segundo grado en las variables x1 , . . . , xn , y1 , . . . , yn , de tal modo que cada variable xi aparece multiplicada por una variable yj . Esta expresi´on polin´omica es equivalente a la expresi´on matricial, de tal manera que toda forma bilineal tiene una expresi´on polin´omica asociada (en una base prefijada). En lo sucesivo, tanto para la expresi´on matricial como para la expresi´on polin´omica, si no se indica ninguna base expl´ıcitamente, se supondr´a, por defecto, que la base es la can´onica. Ejemplo 4.1.4 La forma bilineal sobre R3 definida, en la base can´onica, por    1 2 −1 y1 £ ¤ f (x, y) = x1 x2 x3  0 1 3   y2  −2 0 4 y3 tiene la siguiente expresi´on polin´omica: f (x, y) = x1 y1 + 2x1 y2 − x1 y3 + x2 y2 + 3x2 y3 − 2x3 y1 + 4x3 y3 . Cambio de base Sean B = {v1 , . . . , vn } y B0 = {w1 , . . . , wn } dos bases del espacio vectorial real V y denotemos por M (B, B 0 ) a la matriz de cambio de base de B a B0 . Eso significa que, dados dos vectores x, y ∈ V expresados en la base B como     x1 y1     [x]B =  ...  , [y]B =  ...  xn yn y en la base B0 como



[x]B0

 x01   =  ...  , x0n



[y]B0

 y10   =  ...  yn0

se tiene que 

  x01  ..  0   .  = M (B, B )  x0n

 x1 ..  , .  xn



  y10  ..  0   .  = M (B, B )  yn0

 y1 ..  . .  yn

Comparando las dos expresiones para la forma bilineal f en las bases B y B0 f (x, y) = [x]TB M (f, B)[y]B , f (x, y) = [x]TB0 M (f, B 0 )[y]B0 = [x]TB M (B, B 0 )T M (f, B0 )M (B, B 0 )[y]B conclu´ımos que M (f, B) = M (B, B0 )T M (f, B0 )M (B, B0 ) .

(4.1)

Recordemos que la matriz de cambio de base M (B, B0 ) es invertible. Esto nos conduce a la siguiente definici´on.

´ 4. FORMAS CUADRATICAS

70

Definici´on 4.1.5 Dos matrices A, B ∈ Rn×n son congruentes si existe una matriz invertible P ∈ Rn×n tal que B = P T AP . Al igual que en el caso de la relaci´on de semejanza, la relaci´on de congruencia es sim´etrica (es decir, A, B son congruentes si y s´olo si B, A son congruentes). As´ı pues, una forma bilineal puede ser representada por infinitas matrices (cada una de ellas representa la forma bilineal en una base), pero todas ellas son congruentes. El rec´ıproco es tambi´en cierto: dos matrices congruentes representan la misma forma bilineal en dos bases distintas. As´ı pues: Dos matrices son congruentes si y s´olo si representan la misma forma bilineal en dos bases distintas. Como ya sabemos, el rango de una matriz no var´ıa al multiplicarla por una matriz invertible. En consecuencia, dos matrices congruentes tienen el mismo rango. Por tanto, tiene sentido la siguiente definici´on. Definici´on 4.1.6 El rango de una forma bilineal f sobre V es el rango de la matriz M (f, B), para cualquier base B de V. Problema Resuelto 4.1.7 Dada la forma bilineal sobre R3 f (x, y) = x1 y2 − 2x1 y3 + 7x2 y1 − 9x2 y3 , definida en la base can´onica, hallar su rango y la expresi´on de f en la base       1 0   1      0 , −1 , 2  . B=   1 0 1 Soluci´on. La matriz de f en la base can´onica es   1 −2 0 M (f, CR3 ) =  7 0 −9  , 0 0 0 que tiene rango 2, por lo que el rango de f es 2. Como 

 1 1 0 M (B, CR3 ) =  0 −1 2  1 0 1

tenemos M (f, B) = M (B, CR3 )T M (f, CR3 )M (B, CR3 )     T   1 3 −4 1 1 0 1 −2 0 1 1 0 5 . =  0 −1 2   7 0 −9   0 −1 2  =  3 −4 −4 14 −18 1 0 1 0 0 0 1 0 1

4.2 Formas cuadr´aticas. Reducci´on a la forma can´onica

71

Por tanto, la expresi´on de f en la base B es f (x, y) = x1 y1 + 3x1 y2 − 4x1 y3 + 3x2 y1 − 4x2 y2 + 5x2 y3 − 4x3 y1 + 14x3 y2 − 18x3 y3 . Dado que cualquier matriz cuadrada es la matriz de una forma bilineal, cabe preguntarse si se obtiene alg´un tipo particular de forma bilineal cuando la matriz asociada tiene alguna de las estructuras que hemos visto a lo largo del curso (sim´etrica, unitaria, ortogonal). La que nos interesa, en particular, es la simetr´ıa. Definici´on 4.1.8 Una forma bilineal f sobre V es sim´etrica si f (x, y) = f (y, x) ,

∀ x, y ∈ V.

La matriz asociada a una forma bilineal sim´etrica es, precisamente, una matriz sim´etrica (y viceversa). f es una forma bilineal sim´etrica sobre V si y s´olo si M (f, B) es una matriz sim´etrica para cualquier base B de V.

´ ´ a la forma canonica ´ 4.2 Formas cuadraticas. Reduccion En esta secci´on introduciremos y estudiaremos las formas cuadr´aticas asociadas a las formas bilineales sim´etricas, que ser´an fundamentales en el Tema 5. Definici´on 4.2.1 Dado un espacio vectorial real V y una forma bilineal sim´etrica f : V × V −→ R sobre V, la forma cuadr´atica asociada a f es la aplicaci´on q : V −→ R definida por q(x) = f (x, x). Si la forma bilineal f est´a asociada a la matriz A (en una determinada base) entonces la expresi´on de la forma cuadr´atica asociada es q(x) = xT Ax. Por lo tanto, una forma cuadr´atica puede expresarse mediante un polinomio homog´eneo de segundo grado en n variables x1 , . . . , xn , donde n es la dimensi´on del espacio vectorial V. Ejemplo 4.2.2 Sea f : R2 × R2 −→ R la forma bilineal sim´etrica definida, en la base can´onica, por · ¸· ¸ £ ¤ 1 2 y1 f (x, y) = x1 x2 . 2 −3 y2 La forma cuadr´atica asociada a f , en la base can´onica de R2 , est´a dada por · ¸· ¸ £ ¤ 1 2 x1 q(x) = x1 x2 = x21 + 4x1 x2 − 3x22 . 2 −3 x2

´ 4. FORMAS CUADRATICAS

72

Rec´ıprocamente, todo polinomio homog´eneo de segundo grado en n variables es la expresi´on polin´omica de una forma cuadr´atica sobre el espacio Rn , cuya matriz asociada puede obtenerse f´acilmente, como se muestra en el siguiente ejemplo. Ejemplo 4.2.3 El polinomio homog´eneo de segundo grado en tres variables q(x1 , x2 , x3 ) = −3x21 + 4x1 x2 − 6x1 x3 + 5x22 − x23 est´a asociado a la forma cuadr´atica sobre R3 cuya matriz, en la base can´onica, es 

 −3 2 −3 A= 2 5 0  −3 0 −1 (n´otese que para obtener las entradas no diagonales se dividen entre 2 los coeficientes correspondientes). ´ ´ Forma canonica de una forma cuadratica Del mismo modo que una forma bilineal puede tener distintas expresiones matriciales, tambi´en una forma cuadr´atica, que se obtiene de una forma bilineal haciendo coincidir el primer vector y el segundo, se puede expresar en diferentes bases mediante diferentes matrices. E, igualmente, todas las matrices que representan a una misma forma cuadr´aticas son semejantes. Nuestro inter´es ahora se centra en encontrar la expresi´on (matriz) m´as elemental de una forma cuadr´atica. Puesto que las matrices m´as elementales que existen son las matrices diagonales, diremos que una forma cuadr´atica q est´a expresada en forma can´onica si su matriz asociada es diagonal, es decir, si la expresi´on polin´omica es del tipo q(x) = xT Dx = λ1 x21 + . . . + λn x2n , donde D = diag(λ1 , . . . , λn ). La pregunta que surge inmediatamente es: ¿toda forma cuadr´atica se puede expresar en forma can´onica? La respuesta es SI. Teorema 4.2.4 Toda forma cuadr´atica tiene, al menos, una forma can´onica. Demostraci´on. Sea q una forma cuadr´atica sobre V cuya matriz asociada, en alguna base de V, es A. Como A es sim´etrica, por el Teorema 2.3.15, A es ortogonalmente diagonalizable, es decir, existe una matriz P invertible tal que A = P DP T , donde D es una matriz diagonal. Por tanto: q(x) = xT Ax = xT (P DP T )x = (P T x)T DP T x = y T Dy = q(y), donde y = P T x. La u´ ltima expresi´on est´a en forma can´onica. La segunda cuesti´on, una vez resuelta afirmativamente la primera, consiste en preguntarse si la forma can´onica es u´ nica. En este caso, la respuesta es NO, como veremos en el siguiente contraejemplo.

4.2 Formas cuadr´aticas. Reducci´on a la forma can´onica

73

Ejemplo 4.2.5 Sea q la forma cuadr´atica cuya matriz en la base can´onica es · ¸ 1 0 A= , 0 2 que es diagonal y, por tanto, q(x) = x21 + 2x22 est´a dada en forma can´onica. La matriz · ¸· ¸· ¸ · ¸ −1 1 1 0 −1 4 3 0 B= = 4 2 0 2 1 2 0 24 | {z } | {z } | {z } PT

A

P

representa la misma forma cuadr´atica y es diagonal, luego q(y) = 3y12 + 24y22 es otra forma can´onica de q. Las dos formas can´onicas del ejemplo anterior est´an relacionadas por el cambio de base · ¸ · ¸· ¸ y1 −1 4 x1 = . y2 1 2 x2 El ejemplo anterior pone de manifiesto la importancia de indicar el cambio de base, es decir, la base en la que se est´a expresando la forma cuadr´atica. Esto tiene especial relevancia en la forma can´onica que, como acabamos de ver, no es u´ nica. Es importante destacar que los coeficientes que aparecen en la forma can´onica de q no son (necesariamente) los autovalores de la matriz A que representa la forma cuadr´atica. Este hecho es consecuencia de que los autovalores no son invariantes por congruencia. La demostraci´on del Teorema 4.2.4 nos proporciona un m´etodo para reducir q a su forma can´onica: diagonalizarla por congruencia. Pero este no es el u´ nico m´etodo posible. A continuaci´on, vamos a estudiar por separado los tres m´etodos que recogemos en la siguiente tabla: M´etodos para obtener una forma can´onica de q A = matriz asociada a q 1 Diagonalizar A ortogonalmente: P T AP = D (P −1 = P T ) 2 Diagonalizar A por congruencia: P T AP = D (P −1 6= P T ) 3 Completar cuadrados (m´etodo de Gauss) Hay que destacar que, en el m´etodo 1 , los coeficientes de la forma can´onica (las entradas diagonales de la matriz D) son los autovalores de A, porque la congruencia es, al mismo tiempo, una equivalencia, mientras que que en los m´etodos 2 y 3 no son los autovalores.

´ 4. FORMAS CUADRATICAS

74

El m´etodo 1 es el m´etodo usual de diagonalizaci´on que se conoce del curso previo de ´ Algebra Lineal y que ya se ha tratado en la Secci´on 2.3. A continuaci´on vamos a explicar brevemente los otros dos m´etodos. 2 Diagonalizaci´on por congruencia: Este procedimiento se sintetiza en los siguientes dos pasos: 1) Escalonar la matriz A por filas: P A = T , donde T es triangular superior. 2) Aplicar sobre las columnas de A las mismas operaciones que se han aplicado sobre las filas en el paso 1): P AP T = T P T = D

3 Diagonalizar completando cuadrados: La idea b´asica de este procedimiento es completar cuadrados en cada una de las variables para “absorber” los t´erminos xi xj , con i 6= j. El proceso se completa en n pasos, donde en el paso i-´esimo se completan cuadrados en todos los t´erminos xi xj , para j > i (los otros ya se han completado en pasos anteriores). En el siguiente ejemplo ilustramos los tres m´etodos. Problema Resuelto 4.2.6 Sea q la forma cuadr´atica sobre R3 cuya representaci´on en la base can´onica es q(x) = x21 + x22 + x23 − 4(x1 x2 + x1 x3 + x2 x3 ). Reducir q a una forma can´onica de tres formas diferentes empleando cada uno de los tres m´etodos anteriores e indicar las bases en las que est´a dada cada una de las formas can´onicas. Soluci´on. 1 La matriz asociada a q en la base can´onica es   1 −2 −2 A =  −2 1 −2  . −2 −2 1 Sus autovalores se obtienen a partir del polinomio caracter´ıstico p(λ) = λ3 − 3λ2 − 9λ + 27 −→ λ1 = 3 (doble), λ2 = −3. En este punto ya sabemos que una forma can´onica de q es q(x) = 3x21 + 3x22 − 3x23 . Para hallar la base necesitamos conocer una base ortonormal de autovectores de A. El espacio propio asociado a λ1 = 3 es     1 1 V{λ1 =3} = ker(A − 3I3 ) = Gen  −1  ,  1  . −2 0

4.2 Formas cuadr´aticas. Reducci´on a la forma can´onica

75

La base anterior es ortogonal, as´ı que es inmediato obtener una base ortonormal de autovectores asociados a λ1 = 3:  1   1  √

√

2   u1 =  − √12  , 0

 u2 = 

6 √1 6 − √26

 .

Por otro lado, el espacio propio asociado a λ2 = −3 es 

V{λ2 =−3}

 1 = ker(A + 3I3 ) = Gen  1  , 1

de manera que una base ortonormal de dicho espacio es  1   u3 = 

√ 3 √1 3 √1 3

 .

Por tanto, la base en la que est´a expresada la forma can´onica es B1 = {u1 , u2 , u3 }. £ ¤ Puede comprobarse que, si P1 = u1 u2 u3 entonces P1T AP1 = diag(3, 3, −3). 2 Primero pasamos A a forma escalonada mediante operaciones por filas:       1 −2 −2 1 −2 −2 1 −2 −2 −→ A =  −2 1 −2  F2 + 2F1  0 −3 −6  F3 −−→2F2  0 −3 −6  = T. F3 + 2F1 −2 −2 1 0 −6 −3 0 0 9 Aplicamos a las columnas de T las mismas operaciones que hemos aplicado a las filas de A (en el mismo orden) o, lo que es igual, multiplicamos T a la derecha por P2T , donde   1 0 0 1 0 , P2 =  2 −2 −2 1 para obtener la matriz diagonal 

P2 AP2T = T P2T

 1 0 0 =  0 −3 0  . 0 0 9

As´ı, obtenemos la forma can´onica q(x) = x21 − 3x22 + 9x23 , que, seg´un (4.1), est´a dada en la base que determinan las columnas de P2T , es decir       −2  2  1 B2 =  0  ,  1  ,  −2  .   1 0 0

´ 4. FORMAS CUADRATICAS

76 3 Primero completamos los cuadrados de la variable x1 :

q(x) = (x21 − 4x1 x2 − 4x1 x3 ) + x22 + x23 − 4x2 x3 = (x1 − 2x2 − 2x3 )2 − 4x22 − 4x33 − 8x2 x3 + x22 + x23 − 4x2 x3 = (x1 − 2x2 − 2x3 )2 − 3x22 − 3x33 − 12x2 x3 . Con el cambio de variable   y1 = x1 − 2x2 − 2x3 y2 = x 2  y3 = x 3 obtenemos q(y) = y12 − 3y22 − 3y32 − 12y2 y3 , que a´un no est´a en forma can´onica. El siguiente paso consiste en completar cuadrados en y2 : q(y) = y12 −3(y22 +4y2 y3 )−3y32 = y12 −3(y2 +2y3 )2 +12y32 −3y32 = y12 −3(y2 +2y3 )2 +9y32 . Con el cambio de variable

  z1 = y1 z2 = y2 + 2y3  z3 = y3

obtenemos q(z) = z12 − 3z22 + 9z32 , que ya est´a en forma can´onica. La relaci´on entre las forma    z1 1  z2  =  0 z3 0

nuevas variables y las antiguas se expresa matricialmente en la       0 0 y1 1 0 0 1 −2 −2 x1 1 2   y2  =  0 1 2   0 1 0   x2  0 1 y3 0 0 1 0 0 1 x3   x1 1 −2 −2 2   x2  . = 0 1 x3 0 0 1 {z } | 

P3

Si D = diag(1, −3, 9) tenemos que q(z) = z T Dz = (P3 x)T DP3 x = xT (P3T DP3 )x = xT Ax, de modo que P3T DP3 = A y, seg´un la f´ormula (4.1), P3 = M (CR3 , B3 ), lo que significa que   1 2 −2 M (B, CR3 ) = M (CR3 , B3 )−1 =  0 1 −2  0 0 1

4.3 Formas cuadr´aticas definidas y, por tanto,

77

      2 −2   1 B3 =  0  ,  1  ,  −2    0 0 1

es la base en la que est´a expresada la forma can´onica de q. N´otese que obtenemos la misma forma can´onica y la misma base que usando el procedimiento anterior.

´ 4.3 Formas cuadraticas definidas Como se ha mencionado al principio de la Secci´on 1.1, un producto escalar es una forma bilineal sim´etrica definida positiva. Dicha positividad se refiera en realidad a la forma cuadr´atica asociada a la forma bilineal. A continuaci´on enunciamos la definici´on precisa de esta clase de formas cuadr´aticas, junto con las restantes clases an´alogas. Definici´on 4.3.1 Sea q una forma cuadr´atica sobre V. Diremos que q es (i) Definida positiva si q(x) > 0 para todo x ∈ V − {0V }. (ii) Definida negativa si q(x) < 0 para todo x ∈ V − {0V }. (iii) Semidefinida positiva si q(x) ≥ 0 para todo x ∈ V. (iv) Semidefinida negativa si q(x) ≤ 0 para todo x ∈ V. (v) Indefinida si no cumple ninguna de las anteriores, es decir, si existen x1 , x2 ∈ V tales que q(x1 ) > 0 y q(x2 ) < 0. La definici´on anterior establece una clasificaci´on completa de las formas cuadr´aticas, pues toda forma cuadr´atica est´a en, al menos, uno de los casos (i) − (v). Algunos de ellos incluyen otros. En concreto, si q es definida positiva entonces es semidefinida positiva, y lo mismo ocurre con las negativas. En cambio, otros son disjuntos. Por ejemplo, si q es indefinida no puede ser de ning´un otro tipo. El problema que ahora nos concierne es el de clasificar una forma cuadr´atica dada, q, seg´un la clasificaci´on de la Definici´on 4.3.1. El siguiente Teorema nos da un criterio completo si conocemos una forma can´onica de q. Teorema 4.3.2 Sea q(x) = λ1 x21 + . . . + λn x2n una forma cuadr´atica en forma can´onica. Entonces q es (i) definida positiva si y s´olo si λi > 0 para todo i = 1, . . . , n, (ii) definida negativa si y slo si λi < 0 para todo i = 1, . . . , n, (iii) semidefinida positiva si y s´olo si λi ≥ 0 para todo i = 1, . . . , n y λi0 = 0 para alg´un 1 ≤ i0 ≤ n, (iv) semidefinida negativa si y s´olo si λi ≤ 0 para todo i = 1, . . . , n y λi0 = 0 para alg´un 1 ≤ i0 ≤ n,

´ 4. FORMAS CUADRATICAS

78

(v) indefinida si y s´olo si hay dos sub´ındices 1 ≤ i1 , i2 ≤ n de modo que λi1 < 0 y λi2 > 0. Ejercicio 4.3.3 Demostrar el Teorema 4.3.2. El Teorema 4.3.2 es muy u´ til una vez que hemos reducido q a una forma can´onica, pero podemos plantearnos si existe alg´un criterio para clasificar la forma can´onica a partir de cualquier expresi´on equivalente. La respuesta es afirmativa y para enunciarla necesitamos introducir el siguiente concepto. Definici´on 4.3.4 El menor principal k-´esimo de una matriz cuadrada A = [aij ]ni,j=1 , que denotaremos por Mk (A), es el determinante de la submatriz formada por las primeras k filas y k columnas de A, es decir: ³ ´ Mk (A) = det [aij ]ki,j=1 . El criterio para clasificar una forma cuadr´atica a partir de la matriz en cualquier base es e´ ste. Teorema 4.3.5 (Criterio de los menores principales) Sea q una forma cuadr´atica sobre un espacio vectorial de dimensi´on n cuya expresi´on matricial en una base cualquiera es q(x) = xT Ax. Entonces q es (a) definida positiva si y s´olo si Mk (A) > 0, para todo k = 1, . . . , n, (b) definida negativa si y s´olo si M1 (A) < 0, M2 (A) > 0, M3 (A) < 0, . . . , (−1)n Mn (A) > 0 (es decir, los menores principales tienen signos alternos y el primero es negativo). El criterio del Teorema 4.3.5 no es un criterio completo, pues s´olo permite reconocer dos casos de los cinco posibles, pero es decisivo si la forma cuadr´atica est´a en uno de estos dos casos. Ejemplo 4.3.6 La forma cuadr´atica sobre R3 cuya expresi´on en la base can´onica es q(x) = x21 + 6x1 x2 + 2x1 x3 + 10x22 + 4x2 x3 + 3x23 es definida positiva, porque la matriz asociada   1 3 1 A =  3 10 2  1 2 3 tiene sus tres menores principales positivos: M1 (A) = 1, M2 (A) = 1, M3 (A) = 1. En cambio, la forma cuadr´atica s(x) = x21 + 6x1 x2 + 2x1 x3 + 10x22 + 4x2 x3 + x23 no es ni definida positiva ni definida negativa, pues los menores principales de la matriz asociada   1 3 1 B =  3 10 2  1 2 1 son M1 (B) = 1 > 0, M2 (B) = 1 > 0, M3 (B) = −1 < 0.

4.4 La ley de inercia de Sylvester

79

4.4 La ley de inercia de Sylvester Introducimos, en primer lugar, el concepto fundamental de esta secci´on. Definici´on 4.4.1 La inercia de una matriz sim´etrica A es la terna i(A) = (i+ (A), i0 (A), i− (A)), donde i+ (A) es el n´umero de autovalores positivos de A, i0 (A) es el n´umero de autovalores nulos de A e i− (A) es el n´umero de autovalores negativos de A. La signatura de A es el n´umero σ(A) = i+ (A) − i− (A). En la secci´on anterior hemos visto que la forma can´onica de una forma cuadr´atica q no es u´ nica pero, por otro lado, el Teorema 4.3.2 nos dice que hay ciertos invariantes en todas las formas can´onicas. En concreto, si todos los coeficientes son positivos en una de ellas tambi´en lo ser´an en cualquier otra, y lo mismo ocurre si todos son negativos, o todos positivos (o negativos) y alguno nulo o si, finalmente, hay al menos uno positivo y uno negativo. El siguiente teorema nos dice m´as a´un. Teorema 4.4.2 (Ley de inercia de Sylvester) Dos matrices congruentes tienen la misma inercia. Como ya hemos mencionado anteriormente, los autovalores no se conservan por congruencia. El teorema anterior nos dice que los signos, en cambio, s´ı se conservan. Traducido al lenguaje de las formas cuadr´aticas, lo que nos dice la ley de inercia de Sylvester es que el n´umero de coeficientes positivos, negativos y nulos de una forma can´onica cualquiera de una forma cuadr´atica dada q, no depende de dicha forma can´onica, sino que es un invariante de la forma cuadr´atica. Por tanto, llamaremos inercia de una forma cuadr´atica q a la terna i(q) = (i+ (q), i0 (q), i− (q)), donde i+ (q) es el n´umero de coeficientes positivos de cualquier forma can´onica de q, i0 (q) es el n´umero de coeficientes nulos de cualquier forma can´onica de q e i− (q) es el n´umero de coeficientes negativos de cualquier forma can´onica de q. Ejemplo 4.4.3 La inercia de la forma cuadr´atica del Ejemplo 4.2.6 es i(q) = (2, 0, 1). N´otese que en dos de las formas can´onicas que se han obtenido en dicho problema los coeficientes son distintos, pero no as´ı la inercia. ´ Forma polar asociada a una forma cuadratica Una forma cuadr´atica q sobre un espacio vectorial V est´a asociada, por definici´on, a una forma bilineal sim´etrica f sobre V mediante la relaci´on q(x) = f (x, x). La forma bilineal f se conoce como la forma polar asociada a q y es u´ nica, es decir, que si f1 y f2 son dos formas bilineales sim´etricas de manera que q(x) = f1 (x, x) = f2 (x, x) ∀x ∈ V entonces f1 (x, y) = f2 (x, y),

∀x, y ∈ V.

Obtener q a partir de f es inmediato. Para recuperar f a partir de q lo que haremos ser´a:

´ 4. FORMAS CUADRATICAS

80

1) Calcular la matriz A asociada a q en una base determinada: q(x) = xT Ax. 2) Entonces: f (x, y) = xT Ay (en la misma base).

´ de Choleski 4.5 Factorizacion Recordemos que una forma cuadr´atica q sobre V es definida positiva si q(x, x) > 0 para todo x ∈ V. Si A es la matriz asociada a q en una base cualquiera de V, esto equivale a decir que xT Ax > 0 para todo x ∈ Rn (donde n es la dimensi´on del espacio vectorial V). Por tanto, diremos que una matriz A ∈ Rn×n sim´etrica es definida positiva si xT Ax > 0 para todo x ∈ Rn . En el siguiente ejercicio se pide demostrar que la propiedad de ser definida positiva es invariante por congruencia. Ejercicio 4.5.1 Demostrar que una matriz A ∈ Rn×n sim´etrica es definida positiva si y s´olo si P T AP es sim´etrica definida positiva, para cualquier matriz P ∈ Rn×n invertible. Ya hemos visto que toda matriz sim´etrica es la matriz asociada a una forma bilineal sim´etrica en una base determinada. De la misma forma, toda matriz sim´etrica definida positiva A ∈ Rn×n es la matriz de Gram en alguna base del producto escalar usual en Rn . Como la matriz de Gram del producto escalar usual en Rn respecto a la base B = {v1 , . . . , vn } se obtiene como el producto GB = V T V ¤ £ donde V = v1 . . . vn , llegamos a A ∈ Rn×n es una matriz sim´etrica definida positiva si y s´olo si existe una matriz G ∈ Rn×n invertible tal que A = GT G. La matriz G anterior es una matriz real cualquiera (invertible) de tama˜no n×n. El resultado fundamental de esta secci´on proporciona una descomposici´on similar a la anterior pero usando una matriz triangular. Teorema 4.5.2 (Factorizaci´on de Choleski) Toda matriz A sim´etrica definida positiva se puede factorizar como A = LLT , donde L es triangular inferior con todas sus entradas diagonales positivas. La factorizaci´on anterior se conoce como la factorizaci´on de Choleski de A. Demostraci´on. La siguiente demostraci´on es constructiva. Hemos visto en la Secci´on 4.2 que la matriz sim´etrica A puede diagonalizarse de la forma P AP T = D = diag(λ1 , . . . , λn ) donde las entradas diagonales de D son positivas (λi > 0, i = 1, . . . , n), porque A es definida positiva y esta propiedad es invariante por congruencia (Ejercicio 4.5.1). La matriz P comprende las transformaciones elementales de filas que llevan la matriz A a su forma escalonada.

4.5 Factorizaci´on de Choleski

81

En estas operaciones a cada fila se le a˜nade una combinaci´on lineal de las filas anteriores. La matriz P es, por tanto, triangular inferior con 1’s en las entradas diagonales. As´ı pues √ √ A = P −1 D(P −1 )T = P −1 D(P −1 D)T = LLT , √ √ donde L = P −1 D, y D es la matriz diagonal cuyas entradas ra´ıces √ diagonales √ son las √ cuadradas (positivas) de las entradas diagonales de D, es decir, D = diag( λ1 , . . . , λn ). Problema Resuelto 4.5.3 Encontrar la factorizaci´on de Choleski de la matriz   1 3 1 A =  3 10 2  . 1 2 3 Soluci´on. Como hemos visto en el Ejemplo 4.3.6, la matriz A es definida positiva. Por tanto, existe la factorizaci´on de Choleski de A. Para obtenerla seguimos el procedimiento descrito en la demostraci´on del Teorema 4.5.2. Primero reducimos la matriz A a forma escalonada       1 3 1 1 3 1 1 3 1 −→ A =  3 10 2  F2 − 3F1  0 1 −1  F3−+→F2  0 1 −1  = T. F3 − F1 1 2 3 0 −1 2 0 0 1 Entonces



P AP T = T P T donde

 1 0 0 = D =  0 1 0 , 0 0 1



 1 0 0 P =  −3 1 0  . −4 1 1

Si ahora llamamos



    1 0 0 1 0 0 1 0 0 √ L = P −1 D =  3 1 0   0 1 0  =  3 1 0  , 1 −1 1 0 0 1 1 −1 1

entonces A = LLT es la factorizaci´on de Choleski de A. En el siguiente cuadro resumimos algunos resultados relevantes de este tema, en los que se establece la relaci´on (biun´ıvoca) entre formas bilineales (q) con las propiedades que hemos estudiado (simetr´ıa y positividad) y las estructuras correspondientes de las matrices asociadas (A). q Forma bilineal Forma bilineal sim´etrica Forma bilineal sim´etrica definida positiva

A Cualquiera Sim´etrica Sim´etrica definida positiva

82

´ 4. FORMAS CUADRATICAS

5 ´ ´ CONICAS Y CUADRICAS

´ ´ a la forma canonica ´ 5.1 Conicas en R2 . Reduccion Una c´onica en el plano R2 es el conjunto de soluciones de una ecuaci´on cuadr´atica. Es £ ¤T decir, es el conjunto de puntos del plano x = x1 x2 tales que xT Ax + aT x + b = 0 ,

(5.1)

donde A es una matriz sim´etrica, · A=

a11 a12 a12 a22

¸ ,

£ ¤T a = a1 a2 es un vector (de coeficientes) de R2 y b ∈ R es un n´umero real. Asumiremos, desde un principio, que la matriz A es no nula, pues si A = 0 entonces (5.1) es una ecuaci´on de primer grado. El objetivo de esta secci´on es clasificar las c´onicas, describirlas geom´etricamente y dar un procedimiento para poder identificar una c´onica a partir de su ecuaci´on (5.1). La herramienta fundamental para hacer todo esto es la forma can´onica de la forma cuadr´atica asociada a la ecuaci´on (5.1), es decir, q(x) = xT Ax. Recordemos que la forma can´onica se obtiene a trav´es de una congruencia de la matriz original A. Esta congruencia est´a asociada a un cambio de base (o “cambio de coordenadas” ) en el plano (o en el espacio R3 en el caso que trataremos en la pr´oxima secci´on). Como hemos comentado en el Tema 4, los cambios de base arbitrarios no respetan las distancias. Para eso necesitamos realizar un cambio de base ortonormal, es decir, en t´erminos matriciales, realizar una semejanza ortogonal. Como las c´onicas son objetos definidos en t´erminos de distancias (v´ease 83

´ ´ 5. CONICAS Y CUADRICAS

84

la Secci´on 5.2), para mantener invariante el tipo de c´onica debemos, en principio, restringirnos a las transformaciones ortogonales. Lo que haremos, por tanto, ser´a diagonalizar ortogonalmente la matriz A. Este ser´a el primer paso del procedimiento. El segundo (y u´ ltimo) paso consistir´a en una traslaci´on, que es una transformaci´on que no altera el tipo de c´onica, y nos dar´a su ecuaci´on reducida. Clasificaci´on de las c´onicas (ecuaci´on reducida) I Partimos de la c´onica (5.1). I Paso 1: Diagonalizar ortogonalmente la matriz A : · T

A = P DP ,

D=

λ1 0 0 λ2

¸ .

Con el cambio de variable y = P T x la ecuaci´on (5.1) se transforma en y T Dy + aT P y + b = 0 ⇔ λ1 y12 + λ2 y22 + e a1 y1 + e a 2 y2 + b = 0 , donde e a=

£

e a1 e a2

¤

= aT P . Ahora distinguimos los tres casos siguientes:

CASO I (tipo el´ıptico): λ1 · λ2 > 0 I Paso 2: Completar cuadrados en (5.2):

c

zµ }| ¶ µ ¶ ¶{ µ e a21 e e a2 2 e a22 a1 2 + λ 2 y2 + − + − b = 0. λ1 y1 + 2λ1 2λ2 4λ1 4λ2 Con el cambio de variable (que corresponde a una traslaci´on) (

z1 = y1 +

e a1 2λ1

z2 = y2 +

e a2 2λ2

,

llegamos a la ecuaci´on reducida de la c´onica λ1 z12 + λ2 z22 − c = 0 Entonces

Si c > 0 Si c = 0 Si c < 0

−→ −→ −→

Elipse Un punto Conjunto vac´ıo

(5.2)

5.1 C´onicas en R2 . Reducci´on a la forma can´onica

85

CASO II (tipo hiperb´olico): λ1 · λ2 < 0 I Paso 2: Completando cuadrados en (5.2) y haciendo el mismo cambio de variable que en el CASO I, llegamos a la misma ecuaci´on reducida λ1 z12 + λ2 z22 − c = 0 Ahora Si c 6= 0

−→

Hip´erbola

Si c = 0

−→

Dos rectas

± ± z12 (c/λ1 )³+ z22 (c/λ2´) = 1 p z2 = ± −λ1 /λ2 z1

CASO III (tipo parab´olico): λ1 · λ2 = 0 Sin p´erdida de generalidad, supondremos que λ1 = 0. Distinguimos los casos: i) e a1 6= 0 : I Paso 2: Completamos cuadrados en (5.2): ¶ ¶ µ µ e a22 e a2 2 b − = 0, λ2 y2 + +e a1 y1 + 2λ2 e a1 4e a1 λ2 y con el cambio de variable (que corresponde a una traslaci´on) ( e a2 z1 = y2 + 2λ 2 z2 = y1 +

b e a1

−

e a22 4e a1 λ2

llegamos a la ecuaci´on reducida λ2 z12 + e a1 z2 = 0 La c´onica es, por tanto, una par´abola. ii) e a1 = 0 : I Paso 2: Completamos cuadrados en (5.2):

c

zµ }| ¶{ ¶ µ e a22 e a2 2 − − b = 0, λ2 y2 + 2λ2 4λ2 y con el cambio de variable (que corresponde a una traslaci´on) ½ z1 = y1 e a2 z2 = y2 + 2λ 2 llegamos a la ecuaci´on reducida p λ2 z22 − c = 0 ⇔ z2 = ± c/λ2 Por tanto

´ ´ 5. CONICAS Y CUADRICAS

86 Si (c/λ2 ) > 0 Si c = 0 Si (c/λ2 ) < 0

−→ −→ −→

Dos rectas Una recta Conjunto vac´ıo

N´otese que el caso λ1 = λ2 = 0 ha sido exclu´ıdo de nuestro an´alisis. Este caso corresponde al caso A = 0, que ha sido exclu´ıdo desde el principio. Observaci´on 5.1.1 Para saber el tipo de c´onica que corresponde a la ecuaci´on (5.1) (el´ıptico, hiperb´olico o parab´olico) es suficiente con conocer la inercia de A. Como ya sabemos por la Ley de inercia de Sylvester (Teorema 4.4.2) la inercia es invariante por congruencia. Esto significa que para clasificar la c´onica al primer nivel de clasificaci´on podemos diagonalizar A usando cualquiera de los tres procedimientos que se han descrito en la Secci´on 4.2, es decir, no es necesario emplear transformaciones ortogonales. No obstante, para conocer algunos de los elementos que describiremos en la pr´oxima secci´on (concretamente los ejes) debemos hacer semejanza ortogonal.

´ de los ejes, centro y vertice ´ ´ 5.2 Determinacion de una conica Un sistema de referencia en R2 es una terna R = {O; v1 , v2 } , donde O es un punto (llamado centro u origen del sistema) y v1 , v2 son dos vectores, llamados ejes del sistema. Los sistemas de referencia permiten localizar puntos en el plano af´ın, que como conjunto es el mismo plano R2 pero formado por puntos y vectores (no necesariamente apoyados en el origen de coordenadas). Puesto que las c´onicas no son variedades lineales, hay que estudiarlas en el plano af´ın, y no en el plano vectorial. En este apartado vamos a presentar los elementos singulares de una c´onica y vamos a aprender a localizarlos en el plano af´ın a partir de la ecuaci´on de la c´onica. Dichos elementos singulares son los siguientes. Si la c´onica (5.1) es una elipse o una hip´erbola, entonces son el centro y los ejes, mientras que si es una par´abola, entonces se trata del eje y del v´ertice. Recordemos que, para pasar a la ecuaci´on reducida, primero diagonalizamos A por congruencia ortogonal, de modo que A = P DP T ,

o bien: P −1 AP = D ,

donde P es una matriz ortogonal y D es diagonal: · ¸ λ1 D= , λ1 , λ2 son los autovalores de A. λ2 Recordemos tambi´en que las columnas de P son los autovectores de A. La ecuaci´on reducida que obtenemos al hacer esta transformaci´on es λ1 y12 + λ2 y22 + b1 y1 + b2 y2 + b = 0 , donde, llamando y = [y1 y2 ]T , hemos efectuado el cambio de variable x = P y (⇔ y = P T x) .

(5.3)

5.2 Determinaci´on de los ejes, centro y v´ertice de una c´onica

87

Esto significa que la nueva ecuaci´on est´a dada en una nueva base B, y que P es la matriz de cambio de base de B a Bc (donde Bc es la base can´onica). Es decir: las columnas de P son los vectores de B expresados en la base can´onica. Por otra parte, la traslaci´on que efectuamos en el Paso 2 no altera la direcci´on de los vectores. En resumen, al llevar a cabo los pasos 1 y 2 que transforman la ecuaci´on (5.1) en la ecuaci´on reducida, lo que hacemos es expresar la c´onica en un nuevo sistema de referencia {O; B} que sustituye al sistema de referencia can´onico, en el que est´a expresada la c´onica inicialmente, y que consta del origen de coordenadas y la base can´onica. Los ejes del nuevo sistema de referencia son, por tanto, las columnas de P . Queda por determinar el origen del sistema. Veremos que, en el caso de una elipse o hip´erbola, el origen O coincide con el centro de la c´onica.

´ 5.2.1 Centro y ejes de una elipse o una hiperbola En este apartado escribiremos la ecuaci´on de la c´onica en la forma f (x1 , x2 ) = 0 (es decir, f (x) = xT Ax + aT x + b). Recordemos la definici´on m´etrica de la elipse y la hip´erbola: Elipse: Conjunto de puntos del plano cuya SUMA de distancias a dos puntos fijos (llamados focos) es constante. Hip´erbola: Conjunto de puntos del plano cuya DIFERENCIA de distancias a dos puntos fijos (llamados focos) es constante. En ambos casos, el centro de la c´onica es el punto medio entre los focos (en la recta que los une), uno de los ejes es la recta que contiene los focos y el otro es la recta perpendicular a e´ ste por el centro de la c´onica. A continuaci´on explicamos c´omo hallar el centro y los ejes. Centro: Llamemos C = (c1 , c2 ) al centro de la elipse/hip´erbola. Cuando hacemos una traslaci´on cualquiera de la forma · ¸ · ¸ · ¸ x1 y1 c1 = + , x2 y2 c2 en la ecuaci´on (5.1), el centro de la elipse/hip´erbola va al origen de coordenadas, lo que significa que la nueva ecuaci´on que obtenemos, en las coordenadas y1 , y2 , no tiene t´ermino en y1 ni en y2 (es decir: los t´erminos de primer grado son cero). Por otro lado, puede comprobarse que al hacer el cambio anterior los coeficientes de los t´erminos de primer grado que se obtienen son: Coeficiente de y1 : 2a11 c1 + 2a12 c2 + a1 =

∂f ∂x1 (c1 , c2 )

Coeficiente de y2 : 2a22 c2 + 2a12 c1 + a2 =

∂f ∂x2 (c1 , c2 )

y, como estos dos coeficientes han de ser nulos, llegamos a que (c1 , c2 ) es la soluci´on del sistema lineal: Centro:

∂f ∂x1

=

∂f ∂x2

=0

´ ´ 5. CONICAS Y CUADRICAS

88

Ejes: La traslaci´on no altera la direcci´on de los ejes, por lo que e´ stos tendr´an las direcciones y1 = 0 e y2 = 0 de la ecuaci´on (5.3), y que corresponden a las direcciones de los autovectores de A. Como adem´as los ejes pasan por el centro, tendremos que las ecuaciones param´etricas de los ejes son: · Eje 1: · Eje 2:

x1 x2 x1 x2

¸

· =

¸

· =

c1 c2 c1 c2

¸ + t u1 ¸ + t u2

donde u1 , u2 son los autovectores de A (es decir: las columnas de P ). Si conocemos los autovalores de A y el centro de la elipse/hip´erbola, conocemos su ecuaci´on reducida, que viene dada por λ1 x21 + λ2 x22 + f (c1 , c2 ) = 0. No obstante, para conocer el sistema de referencia en el que est´a expresada la ecuaci´on anterior, es necesario conocer los autovectores de A.

´ ´ 5.2.2 Eje y vertice de una parabola Recordemos la definici´on m´etrica de par´abola: Par´abola: Conjunto de puntos del plano que equidistan de un punto (llamado foco) y una recta (llamada eje). El v´ertice de la par´abola es el punto medio entre el foco y el punto de corte con el eje de la recta perpendicular al eje por el foco. Es un punto de la c´onica. A continuaci´on explicamos c´omo hallar el eje y el v´ertice. Supondremos, como antes, que λ1 = 0. V´ertice: En el v´ertice, la pendiente de la c´onica coincide con la pendiente del eje, que es la pendiente del autovector asociado al autovalor λ2 (el que no es nulo). Llamaremos a esta pendiente m2 y recordemos que la pendiente es la derivada de x2 con respecto a x1 : ∂x2 = m2 ∂x1 Por otro lado, si derivamos (impl´ıcitamente) la ecuaci´on f (x1 , x2 ) = 0 obtenemos ∂f ∂f ∂x2 + = 0. ∂x1 ∂x2 ∂x1 Como, adem´as, el v´ertice est´a en la par´abola, conclu´ımos que el v´ertice, (v1 , v2 ), es la soluci´on del sistema (no lineal)

5.2 Determinaci´on de los ejes, centro y v´ertice de una c´onica

89

 ∂f ∂f   ∂x + m2 ∂x = 0 1

V´ertice:

 

2

f (x1 , x2 ) = 0

Eje: El eje (que corresponde al eje de simetr´ıa de la par´abola), est´a dado por la ecuaci´on y2 = 0, que corresponde al autovector asociado al autovalor cero. Como la traslaci´on no altera la direcci´on del eje, se concluye que el eje tiene la direcci´on del autovector asociado al autovalor cero (que llamaremos u1 ). Como, adem´as, pasa por el v´ertice, su ecuaci´on param´etrica es: · Eje:

x1 x2

¸

· =

v1 v2

¸ + t u1

(donde (v1 , v2 ) son las coordenadas del v´etice). Problema Resuelto 5.2.1 Hallar la ecuaci´on reducida de la c´onica 4x21 − 4x1 x2 + x22 − 6x1 − 2x2 + 18 = 0 y clasificarla. Hallar sus elementos notables (centro/v´ertice y eje(s)). Soluci´on. Matricialmente la c´onica se expresa en la forma · ¸· ¸ · ¸ £ ¤ £ ¤ x1 4 −2 x1 C ≡ x1 x2 − 6 2 + 18 = 0. −2 1 x2 x2 | {z } A

Diagonalizamos la matriz A ortogonalmente como ya se ha descrito en la Secci´on 4.2: √ √ ¸· √ ¸ ¸· √ · 2/ √5 1/√5 5 0 2/√5 −1/√ 5 . A= 0 0 −1/ 5 2/ 5 1/ 5 2/ 5 {z } {z } | | P

Con el cambio de variable

PT

·

y1 y2

¸

· =P

T

x1 x2

¸

obtenemos la ecuaci´on de la c´onica · ¸· ¸ · ¸ √ £ ¤ 5 0 £ ¤ y1 y1 −2 5 1 1 C ≡ y1 y2 + 18 = 0. 0 0 y2 y2 Se trata, por tanto, de una c´onica de tipo parab´olico. Completando cuadrados en la ecuaci´on anterior: Ã √ !2 µ ¶ √ 5 17 C ≡ 5 y1 − − 2 5 y2 − √ 5 2 5

´ ´ 5. CONICAS Y CUADRICAS

90 y efectuando el cambio de variable   z1 = y1 −  z2 = y2 − llegamos a la ecuaci´on reducida

√

5 5 17 √ 2 5

√ 5z12 − 2 5z2 = 0.

El v´ertice de C es el punto ·

z1 z2

¸

· =

0 0

¸

· ⇒

y1 y2

¸

 =

√

5 5 17 √ 2 5

 ⇒

·

x1 x2



¸

= PT 

√

5 5 17 √ 2 5

 =

"

21 10 16 5

#

mientras que el eje es la recta que pasa por el v´ertice y cuya direcci´on es el autovector asociado al autovalor 0, es decir, " 21 # · ¸ 1 10 +t r≡ 16 2 5

Tambi´en pod´ıamos haber calculado el v´ertice mediante el procedimiento que acabamos de estudiar, es decir, como la soluci´on del sistema ∂f ∂f + m2 = 0 = f (x1 , x2 ). ∂x1 ∂x2 En este caso, dicho sistema (no lineal) es ½ 8x1 − 4x2 − 6 − 12 (−4x1 + 2x2 − 2) = 0 4x21 − 4x1 x2 + x22 − 6x1 − 2x2 + 18 = 0 Despejando x2 = 2x1 − 1 en la primera ecuaci´on y sustituyendo en la segunda obtenemos 4x21 −4x1 (2x1 −1)+(2x1 −1)2 −6x1 −2(2x1 −1)+18 = 0 ⇔ −10x1 +21 = 0 ⇔ x1 =

21 , 10

y sustituyendo el valor de x1 tenemos que x2 = 16/5.

´ ´ a la forma canonica ´ 5.3 Cuadricas en R3 . Reduccion Una cu´adrica en el espacio R3 es el conjunto de soluciones de una ecuaci´on cuadr´atica. Es £ ¤T decir, es el conjunto de puntos del espacio x = x1 x2 x3 tales que xT Ax + aT x + b = 0 , donde A es una matriz sim´etrica,  a11 a12 a13 A =  a12 a22 a23  , a13 a23 a33 

(5.4)

5.3 Cu´adricas en R3 . Reducci´on a la forma can´onica

91

£ ¤T a = a1 a2 a3 es un vector (de coeficientes) de R3 y b ∈ R es un n´umero real. Como en el caso de las c´onicas, supondremos que la matriz A es no nula, pues si A = 0 entonces (5.4) es una ecuaci´on de primer grado. El objetivo de esta secci´on es clasificar las cu´adricas, describirlas geom´etricamente y dar un procedimiento para poder identificar una cu´adrica a partir de su ecuaci´on (5.4). La herramienta fundamental para hacer todo esto es, como en el caso de las c´onicas, la forma can´onica de la forma cuadr´atica q(x) = xT Ax. El primer paso consistir´a, de nuevo, en una semejanza ortogonal para diagonalizar la matriz A y el segundo en una traslaci´on despu´es de completar cuadrados. Clasificaci´on de las cu´adricas (ecuaci´on reducida) I Partimos de la cu´adrica (5.4). I Paso 1: Diagonalizar ortogonalmente la matriz A :   λ1 0 0 A = P DP T , D =  0 λ2 0  . 0 0 λ3 Con el cambio de variable y = P T x la ecuaci´on (5.1) se transforma en y T Dy + aT P y + b = 0 ⇔ λ1 y12 + λ2 y22 + λ3 y32 + e a1 y1 + e a 2 y2 + e a3 y3 + b = 0 , donde e a=

£

e a1 e a2 e a3

¤

(5.5)

= aT P . Ahora distinguimos los tres casos siguientes:

CASO I : λ1 · λ2 · λ3 6= 0 Tenemos los casos (a) Los tres autovalores del mismo signo. Supondremos que λ1 > 0, λ2 > 0, λ3 > 0. I Paso 2: Completamos cuadrados en (5.5) µ ¶ µ ¶ µ ¶ µ 2 ¶ e a1 2 e a2 2 e a3 2 e a1 e a22 e a23 λ1 x1 + +λ2 x2 + +λ3 x3 + − + + −b =0 2λ1 2λ2 2λ3 4λ1 4λ2 4λ3 | {z } c

y hacemos el cambio de variable     y1 = x 1 + y2 = x 2 +    y =x + 3 3

e a1 2λ1 e a2 2λ2 e a3 2λ3

para llegar a la ecuaci´on reducida λ1 y12 + λ2 y22 + λ3 y32 = c En funci´on del signo de c tenemos

´ ´ 5. CONICAS Y CUADRICAS

92 Si c > 0 Si c = 0 Si c < 0

−→ −→ −→

Elipsoide Un punto Conjunto vac´ıo

(b) Dos autovalores del mismo signo y el tercero de signo distinto. Supondremos que λ1 > 0, λ2 > 0, λ3 < 0. I Paso 2: Completamos cuadrados en (5.5) y, con el mismo cambio de variable que en el caso anterior, llegamos a la misma ecuaci´on reducida. Ahora tenemos los casos siguientes en funci´on del valor de c: Si c 6= 0 Si c = 0

−→ −→

Hiperboloide Cono El´ıptico

CASO II : λ1 · λ2 6= 0 , λ3 = 0 (Un autovalor nulo y los restantes no nulos) Tenemos los dos casos siguientes (c) Los dos autovalores del mismo signo. Supondremos que λ1 > 0, λ2 > 0, λ3 = 0. I Paso 2: Completamos cuadrados en (5.5) µ ¶ µ ¶ µ 2 ¶ e a1 2 e a2 2 e a1 e a22 λ1 y1 + + λ2 y2 + +e a3 y3 − + −b =0 2λ1 2λ2 4λ1 4λ2 {z } | c

y hacemos el cambio de variable    z1 = y1 + z2 = y2 +   y3 fijo

e a1 2λ1 e a2 2λ2

para llegar a la ecuaci´on reducida λ1 z12 + λ2 z22 + e a3 y3 = c Ahora distinguimos los casos (c1) e a3 6= 0: Haciendo el cambio de variable (traslaci´on)   z1 fijo z2 fijo  z3 = y3 − eac3 llegamos a λ1 z12 + λ2 z22 + e a3 z3 = 0 La cu´adrica es un paraboloide el´ıptico.

5.3 Cu´adricas en R3 . Reducci´on a la forma can´onica

93

(c2) e a3 = 0: La ecuaci´on reducida es λ1 z12 + λ2 z22 = c, y en funci´on del signo de c Si c > 0 Si c = 0 Si c < 0

−→ −→ −→

Cilindro el´ıptico Recta Conjunto vac´ıo

(d) Los dos autovalores de distinto signo. Supondremos que λ1 > 0, λ2 < 0, λ3 = 0. I Paso 2: Completamos cuadrados en (5.5) y haciendo el mismo cambio de variable que en el apartado (c) llegamos a la misma ecuaci´on reducida. Ahora distinguimos los casos (d1) e a3 6= 0: Haciendo el cambio de variable (traslaci´on)   z1 fijo z2 fijo  z3 = y3 − eac3 llegamos a λ1 z12 + λ2 z22 + e a3 z3 = 0 La cu´adrica es un paraboloide hiperb´olico. (d2) e a3 = 0: La ecuaci´on reducida es λ1 z12 + λ2 z22 = c, y en funci´on del signo de c Si c 6= 0 Si c = 0

−→ −→

Cilindro hiperb´olico Dos planos no paralelos

CASO III : λ1 6= 0, λ2 = λ3 = 0 (Dos autovalores nulos y el tercero no nulo) I Paso 2: Completamos cuadrados en (5.5) µ ¶ µ 2 ¶ e a1 2 e a1 λ1 y1 + +e a2 y2 + e a3 y3 − −b =0 2λ1 4λ1 | {z } c

y hacemos el cambio de variable   z1 = y1 + y fijo  2 y3 fijo

e a1 2λ1

para llegar a la ecuaci´on reducida λ1 z12 + e a2 y2 + e a3 y3 = c Supondremos que

´ ´ 5. CONICAS Y CUADRICAS

94 (e) λ1 > 0, λ2 = λ3 = 0 y distinguimos los casos

(e1) e a2 6= 0 o´ e a3 6= 0: Podemos encontrar una transformaci´on unitaria y una traslaci´on de manera que al aplicar estos cambios de variable llegamos a la ecuaci´on λ1 z12 + nz2 = 0, donde n ∈ R. La c´onica es un cilindro parab´olico. Vamos a describir el cambio de variable. La variable z1 quedar´a fija y s´olo cambiaremos las otras dos. El primer cambio de variable corresponde a una transformaci´on ortogonal y es   z1 fijo w2 = ean2 y2 + ean3 y3 ,  e a3 e a2 w3 = − n y2 + n y3 √ donde n = e a2 + e a3 . Con este cambio la c´onica se reduce a λ1 z12 + nw2 = c. Finalmente, con la traslaci´on   z1 fijo z2 = w2 −  z3 = w3

c n

,

obtenemos la ecuaci´on indicada arriba. (e2) e a2 = e a3 = 0: La ecuaci´on reducida es λ1 z12 = c y en funci´on del signo de c tenemos Si c > 0 Si c = 0 Si c < 0

−→ −→ −→

Dos planos paralelos Plano (doble) Conjunto vac´ıo

Problema Resuelto 5.3.1 Clasificar la cu´adrica C ≡ x21 + 2x22 + 4x23 − 2x1 + 8x2 + 5 = 0. Soluci´on. La expresi´on matricial de la cu´adrica es      x 1 0 0 x 1 1 £ ¤ £ ¤ C ≡ x1 x2 x3  0 2 0   x2  + −2 8 0  x2  + 5 = 0 . x3 0 0 4 x3 | {z } A

I Paso 1: La matriz A ya est´a en forma diagonal. Como los tres autovalores son positivos, C es de tipo el´ıptico. I Paso 2: Completamos cuadrados en la ecuaci´on de partida (x1 − 1)2 + 2(x2 + 2)2 + 4x23 − 4 = 0.

5.3 Cu´adricas en R3 . Reducci´on a la forma can´onica

95

Con el cambio de variable (que corresponde a una traslaci´on)   y1 = x1 − 1 y2 = x2 + 2  y3 = x3 la cu´adrica se convierte en C ≡ y12 + 2y22 + 4y32 − 4 = 0. Se trata, por lo tanto, de un elipsoide. Problema Resuelto 5.3.2 Clasificar la cu´adrica C ≡ 2x1 x2 − 6x1 + 10x2 − x3 + 6 = 0. Soluci´on. La expresi´on matricial de la ecuaci´on de C en la base can´onica es      0 1 0 x x 1 1 £ ¤ £ ¤ C ≡ x1 x2 x3  1 0 0   x2  + −6 10 1  x2  + 6 = 0 . x3 0 0 0 x3 | {z } A

I Paso 1: Los autovalores de A se obtienen f´acilmente y son λ1 = 1, λ2 = −1 y λ3 = 0. Por tanto la cu´adrica es del tipo (d) en la clasificaci´on anterior. Los autoespacios correspondientes son     1 1 1 1 V{λ1 =1} = ker(A − I) = √  1  , V{λ2 =−1} = ker(A + I) = √  −1  , 2 0 2 0 

V{λ3 =0}

 0 = ker(A) =  0  , 1

donde ya hemos calculado en cada caso un vector generador unitario. Como puede verse, los autovectores son ortogonales dos a dos. Por tanto, si   1 √ √1 − 0 2   2 √1 0  P =  √12 2 0 0 1 entonces P T AP = D = diag(1, −1, 0), o bien, A = P DP T , lo que significa que con el cambio de variable     x1 y1  y2  = P T  x2  x3 y3 la cu´adrica se convierte en  x1 4 1 P  x2  +6 = 0 = 0 ⇔ y12 −y22 + √ y1 − √ y2 +y3 +6 = 0. 2 2 x3 

£

C ≡ y12 −y22 + −6 10 1

¤

´ ´ 5. CONICAS Y CUADRICAS

96

I Paso 2: Completamos cuadrados en la ecuaci´on anterior µ ¶ µ ¶ 2 2 8 2 y1 + √ − y2 − √ − y3 + 36 = 0, 2 2 y con el cambio de variable (traslaci´on)  2   z 1 = y1 + √ 2 z2 = y2 − √82   z = y − 36 3 3 la ecuaci´on reducida de la c´onica es C ≡ z12 − z22 − z3 = 0. Se trata, por lo tanto, de un paraboloide hiperb´olico.

´ de los ejes, centro y vertice ´ ´ 5.4 Determinacion de una cuadrica De igual forma que en la secci´on anterior hemos introducido el concepto de sistema de referencia en el plano, podemos definir un sistema de referencia en R3 como una cuaterna R = {O; v1 , v2 , v3 }, donde O es un punto de R3 (llamado origen del sistema) y v1 , v2 , v3 son tres vectores, llamados ejes del sistema. Diremos que el sistema es ortogonal si {v1 , v2 , v3 } es un conjunto ortogonal, y que es ortonormal si los vectores son, adem´as, unitarios. Un cambio de coordenadas y = P x, donde P ∈ R3×3 , est´a asociado a un cambio de referencia en R3 donde el centro queda fijo. Si la matriz P es ortogonal y el sistema de referencia de partida es ortonormal, entonces el nuevo sistema de referencia es tambi´en ortonormal. Por otro lado, un cambio de base z = y + b, donde b ∈ R3 , debido a una traslaci´on, tambi´en corresponde a un cambio de sistema de referencia, pero en este caso los ejes quedan fijos y cambia u´ nicamente el centro. En consecuencia, los cambios de variable que llevamos a cabo para alcanzar la forma reducida de una cu´adrica corresponden a cambios de referencia en los que la nueva referencia es ortonormal, de tal modo que, en esta nueva referencia, la cu´adrica tiene una expresi´on elemental y, tambi´en, una representaci´on gr´afica sencilla (que veremos en la Secci´on 5.5). El origen de coordenadas del nuevo sistema de referencia es un punto distinguido. Si la cu´adrica es un elipsoide o un hiperboloide (CASO I de la clasificaci´on anterior) entonces dicho punto se conoce como el centro de la cu´adrica (por razones geom´etricas que son evidentes a partir de la representaci´on gr´afica) y no es un punto de la cu´adrica. En cambio, si la cu´adrica es un paraboloide (Casos II (c1) y II (d1) en la clasificaci´on anterior), el origen de coordenadas s´ı forma parte de la cu´adrica y se le llama v´ertice del paraboloide. A continuaci´on, vamos a ver c´omo se pueden calcular estos dos puntos.

5.4.1 Centro de un elipsoide o un hiperboloide Si la cu´adrica C ≡ f (x1 , x2 , x3 ) = 0 es un elipsoide o un hiperboloide el centro de la cu´adrica es el origen de coordenadas del sistema de referencia en el que est´a expresada la ecuaci´on reducida, es decir, z = 0R3 . Deshaciendo los cambios para obtener las coordenadas en el sistema de referencia can´onico tenemos y = z − u = −u ⇒ x = P y = −P u,

5.4 Determinaci´on de los ejes, centro y v´ertice de una cu´adrica

97

donde u es el vector de traslaci´on y P es la matriz ortogonal que diagonaliza la matriz A de la ecuaci´on de partida. No obstante, hay otra forma de calcular el centro sin necesidad de conocer las transformaciones que llevan la cu´adrica a su ecuaci´on reducida. Igual que en el caso de la elipse y la hip´erbola, estudiado en la Secci´on 5.2.1, el centro es el punto para el cual el gradiente de f se anula, es decir, es la soluci´on del sistema lineal de 3 ecuaciones y 3 inc´ognitas Centro:

∂f ∂x1

=

∂f ∂x2

=

∂f ∂x3

=0

A partir de los autovalores de A y del centro de la cu´adrica podemos obtener la ecuaci´on reducida. Dicha ecuaci´on es C ≡ λ1 z12 + λ2 z22 + λ3 z32 + f (c1 , c2 , c3 ) = 0, donde (c1 , c2 , c3 ) son las coordenadas del centro (en la base can´onica).

´ 5.4.2 Vertice de un paraboloide Si la cu´adrica C ≡ f (x1 , x2 , x3 ) = 0 es un paraboloide el v´ertice de la cu´adrica, igual que en el caso anterior, es el origen de coordenadas del sistema de referencia en el que est´a expresada la ecuaci´on reducida, es decir, z = 0R3 . Deshaciendo los cambios para obtener las coordenadas en el sistema de referencia can´onico tenemos y = z − u = −u ⇒ x = P y = −P u, donde u es el vector de traslaci´on y P es la matriz ortogonal que diagonaliza la matriz A de la ecuaci´on de partida. Tambi´en lo podemos obtener teniendo en cuenta que el gradiente de f evaluado en el v´ertice es paralelo al autovector asociado al autovalor λ = 0 y que, adem´as, es un punto de la cu´adrica. Por tanto, el v´ertice es la soluci´on del sistema (no lineal) ½ ∇f (x1 , x2 , x3 ) = tu0 V´ertice: f (x1 , x2 , x3 ) = 0 donde u0 es un autovector de A asociado al autovalor λ = 0 y t es un par´ametro que hay que calcular. Problema Resuelto 5.4.1 Hallar, de dos formas distintas, el centro del elipsoide del Problema Resuelto 5.3.1. Soluci´on. La primera forma de hallar el centro es deshaciendo el cambio de variable:         x1 0 1 y1  y2  =  0  ⇒  x2  =  −2  x3 0 0 y3 La segunda forma es resolviendo el sistema lineal ∇f = 0, es decir:      1 x1  2x1 − 2 = 0 4x2 + 8 = 0 ⇒  x2  =  −2   0 8x3 = 0 x3

´ ´ 5. CONICAS Y CUADRICAS

98

Problema Resuelto 5.4.2 Hallar, de dos formas distintas, el v´ertice paraboloide del Problema Resuelto 5.3.2. Soluci´on. La primera forma de hallar el centro es deshaciendo los cambios de variable:   √2             − 2 0 y1 z1 −5 x1 y1   z2  =  0  ⇒  y2  =   √82  ⇒  x2  = P  y2  =  3  0 y3 z3 x3 y3 36 36 La segunda forma consiste en resolver el sistema no lineal   2x2 − 6 = t · 0 x2 = 3       2x1 + 10 = t · 0 x1 = 5 ⇔ −1 = t · 1 t = −1       2x1 y1 − 6x1 + 10x2 − x3 + 6 = 0 x3 = 36

´ grafica ´ ´ 5.5 Representacion de las cuadricas en forma reducida A continuaci´on mostramos las gr´aficas de las ecuaciones reducidas de las cu´adricas estudiadas en la secci´on anterior. CASO I: λ1 · λ2 · λ3 6= 0 (Cu´adricas con centro) Elipsoide

λ1 x2 + λ2 y 2 + λ3 z 2 = c (λ1 > 0, λ2 > 0, λ3 > 0, c > 0) Hiperboloide de una hoja

λ1 x2 + λ2 y 2 + λ3 z 2 = c (λ1 > 0, λ2 > 0, λ3 < 0, c > 0)

5.5 Representaci´on gr´afica de las cu´adricas en forma reducida Hiperboloide de dos hojas

λ1 x2 + λ2 y 2 + λ3 z 2 = c (λ1 < 0, λ2 < 0, λ3 > 0, c > 0)

Cono el´ıptico

λ1 x2 + λ2 y 2 + λ3 z 2 = 0 (λ1 > 0, λ2 > 0, λ3 < 0)

CASO II: λ1 · λ2 6= 0, λ3 = 0 (Cu´adricas con v´ertice) Paraboloide hiperb´olico

λ1 x2 + λ2 y 2 = bz (λ1 · λ2 < 0, b 6= 0)

99

100

´ ´ 5. CONICAS Y CUADRICAS

Paraboloide el´ıptico

λ1 x2 + λ2 y 2 = bz (λ1 · λ2 > 0, b 6= 0)

CASO III: λ1 6= 0, λ2 = λ3 = 0

Cilindro el´ıptico

λ1 x2 + λ2 y 2 = c (λ1 > 0, λ2 > 0, b > 0)

Cilindro hiperb´olico

λ1 x2 + λ2 y 2 = c (λ1 · λ2 < 0, c 6= 0)

5.5 Representaci´on gr´afica de las cu´adricas en forma reducida

101

Cilindro parab´olico

λ1 x2 + by = 0 (λ1 6= 0, b 6= 0)

Ejercicio 5.5.1 (a) La gr´afica del paraboloide hiperb´olico que aparece entre las figuras anteriores corresponde al caso en que el signo del coeficiente b coincide con el signo de uno de los dos autovalores. ¿Cu´al de ellos? ¿C´omo ser´ıa la gr´afica en el otro caso? Justifica la respuesta. (b) La gr´afica del paraboloide el´ıptico corresponde al caso en que el signo del coeficiente b es ¿igual o distinto al signo de los autovalores? ¿C´omo ser´ıa la gr´afica en el otro caso? Justifica la respuesta. (c) La gr´afica del cilindro hiperb´olico corresponde al caso en que el signo del t´ermino independiente c es ¿igual o distinto al signo de los autovalores? ¿C´omo ser´ıa la gr´afica en el otro caso? Justifica la respuesta. (d) La gr´afica del cilindro parab´olico corresponde al caso en que el signo del coeficiente b es ¿igual o distinto al signo del autovalor λ1 ? ¿C´omo ser´ıa la gr´afica en el otro caso? Justifica la respuesta.

102

´ ´ 5. CONICAS Y CUADRICAS

6 NORMAS DE VECTORES Y MATRICES. SENSIBILIDAD DE MATRICES

6.1 Normas vectoriales Una norma sobre un espacio vectorial es una herramienta que permite “medir”. Por ejemplo, medir el tama˜no de un vector o la distancia entre dos vectores. Tambi´en, como hemos visto en la Secci´on 1.5, a trav´es de las normas podemos medir la distancia de un vector a un subespacio vectorial. En esta secci´on nos vamos a centrar en las normas sobre el espacio Rn . Las llamaremos normas vectoriales para distinguirlas de las normas matriciales que estudiaremos en siguiente la secci´on. Definici´on 6.1.1 Una norma en Rn es una aplicaci´on k · k: Rn −→ [0, ∞) que cumple las tres propiedades siguientes i) k v k> 0 si v 6= 0Rn y k 0Rn k= 0. ii) k λv k= |λ| k v k, para todo v ∈ Rn y λ ∈ R. iii) k u + v k≤k u k + k v k, para cualesquiera u, v ∈ Rn .

(desigualdad triangular)

El concepto de norma asociada a un producto escalar ha sido ya introducido en la Definici´on 1.1.9 (para espacios vectoriales cualesquiera). Estas son un tipo particular de normas, pero no son las u´ nicas. Ejemplo 6.1.2 En este ejemplo consideramos el espacio Rn y un vector cualquiera v = ¤T £ v1 . . . vn . Las siguientes aplicaciones son ejemplos de normas sobre Rn : 103

104

6. NORMAS DE VECTORES Y MATRICES. SENSIBILIDAD DE MATRICES

1) Las normas inducidas por productos escalares k v k=< v, v >1/2 . 2) La norma-1: k v k1 = |v1 | + . . . + |vn |. 3) La norma-p, donde 1 < p < ∞: k v kp = (|v1 |p + . . . + |vn |p )1/p . 4) La norma-∞: k v k∞ = m´ax{|v1 | + . . . + |vn |}. Las normas k · k1 y k · k∞ son dos ejemplos de normas vectoriales que no est´an ´ producto escalar. inducidas por ningun

6.2 Normas matriciales El conjunto Rn×m de matrices reales de tama˜no n × n es un espacio vectorial sobre R de dimensi´on n · m. Por tanto, se puede identificar con Rn·m y utilizar cualquier norma vectorial para ”medir” el tama˜no de una matriz en Rn×m . No obstante, la identificaci´on anterior de Rn×m con Rn·m deja fuera una operaci´on importante en el conjunto de matrices, que es el producto de matrices. Esta operaci´on es fundamental porque, recordemos, se corresponde con la composici´on de aplicaciones lineales. La definici´on formal de norma matricial tiene en cuenta esta operaci´on. Definici´on 6.2.1 Una funci´on k · k Rn×n −→ R es una norma matricial si para dos matrices cualesquiera A, B ∈ Rn×n se cumplen los siguientes cuatro axiomas: i) k A k> 0 si A 6= 0Rn×m y k 0Rn×m k= 0 ii) k λA k= |λ| k A k, para cualquier λ ∈ R. iii) k A + B k≤k A k + k B k.

(desigualdad triangular)

iv) k AB k≤k A k · k B k Los axiomas i) − iii) de la definici´on anterior son los mismos axiomas que definen a las normas vectoriales. El axioma iv) permite comparar el tama˜no del producto AB con los tama˜nos de A y B (a trav´es de su producto). Ejercicio 6.2.2 Demostrar que para cualquier norma matricial k · k en Rn×n se cumple que k In k≥ 1. Nuestro primer ejemplo de norma matricial es un ejemplo de norma vectorial sobre Rn·m que es norma matricial porque satisface la propiedad iv) de la definici´on. Definici´on 6.2.3 La norma de Frobenius de una matriz A = [aij ] ∈ Rn×n es el n´umero real  k A kF = 

n X

i,j=1

1/2 a2ij 

6.2 Normas matriciales

105

La norma de Frobenius es la norma-2 en el espacio Rn·m . Consideremos ahora la matriz A ∈ Rn×n como una aplicaci´on de Rn A : Rn −→ Rn , v 7→ Av y sea k · k una norma vectorial en Rn . Queremos “medir” c´omo var´ıan de tama˜no los vectores de Rn al aplicarles la transfomaci´on A. Es decir, queremos estimar el cociente k Av k . kvk Para acotar la variaci´on de tama˜no que puede experimentar un vector cualquiera calculamos el supremo del anterior cociente sobre todos los vectores de Rn . Definici´on 6.2.4 La norma matricial inducida por la norma vectorial k · k es la aplicaci´on k · k: Rn×n −→ R definida por k A k= m´axn x∈R x 6= 0

k Ax k . kxk

Hemos empleado la misma notaci´on para la norma vectorial y para la norma matricial inducida. Es importante tener presente, no obstante, que representan normas distintas, pues una es una norma matricial (aunque eso a´un no lo hemos demostrado) y la otra es una norma vectorial. La norma matricial inducida de A es una cota superior (ajustada) de la variaci´on de tama˜no que puede experimentar un vector cualquiera de Rn cuando sobre dicho espacio act´ua la aplicaci´on lineal asociada a la matriz A. El siguiente ejercicio nos dice que dicha cota se alcanza en el conjunto de vectores unitarios y, por tanto, para calcular la norma inducida es suficiente restringirse a ese conjunto. Ejercicio 6.2.5 Demostrar que la definici´on de norma inducida es equivalente a k A k=

m´ax

x ∈ Rn k x k= 1

k Ax k .

(Pista: Usar el axioma ii) de las normas vectoriales). Entre las siguientes propiedades de las normas inducidas se incluyen los cuatro axiomas que muestran que, efectivamente, son normas matriciales. Propiedades de las normas inducidas: Dadas dos matrices cualesquiera A, B ∈ Rn×n , un vector v ∈ Rn y k · k una norma matricial inducida en Rn×n se tiene 1. k Av k≤k A kk v k 2. k In k= 1

106

6. NORMAS DE VECTORES Y MATRICES. SENSIBILIDAD DE MATRICES

3. k A k≥ 0 y k A k= 0 si y s´olo si A = 0 4. k λA k= |λ| k A k, para todo λ ∈ R 5. k A + B k≤k A k + k B k 6. k AB k≤k A k · k B k Demostraci´on. (de las propiedades anteriores) 1. Es inmediato a partir de la definici´on de norma inducida. 2. k In k= m´axkxk=1 k In x k= m´axkxk=1 k x k= 1. 3. k A k= 0 ⇔k Ax k= 0, ∀x 6= 0 ⇔ Ax = 0, ∀x 6= 0 ⇔ A = 0. 4. Usando el Ejercicio 6.2.5: k λA k= m´ax k λAx k= m´ax |λ| k Ax k= |λ| m´ax k Ax k= |λ| k A k . kxk=1

kxk=1

kxk=1

5. Es consecuencia inmediata de la desigualdad |x + y| ≤ |x| + |y|, que implica m´ ax (|x + y|) ≤ m´ ax (|x| + |y|)

x,y∈S

x,y∈S

donde S es un subconjunto cualquiera de Rn . 6. Tenemos k AB k= m´ax k (AB)x k= m´ax k A(Bx) k≤ m´ ax k A kk Bx k=k A kk B k, kxk=1

kxk=1

kxk=1

donde para la segunda desigualdad hemos usado que k A(Bx) k≤k A kk B k, que es consecuencia de la propiedad 1.

Obs´ervese que, en particular, la propiedad 2 parece indicar que la igualdad k In k= 1 no es cierta para cualquier norma matricial, y as´ı es. Un ejemplo en el que esto no ocurre es la norma de Frobenius (para n > 1). Algunas de las normas inducidas tienen una expresi´on sencilla, como veremos en el siguiente resultado. Teorema 6.2.6 (1) La norma matricial en Rn×n inducida por la norma vectorial k · k1 est´a dada por n X k A k1 = m´ax |aij | = m´ax k ai k1 , 1≤j≤n

donde ai es la i-´esima columna de A.

i=1

1≤i≤n

6.2 Normas matriciales

107

(2) La norma matricial en Rn×n inducida por la norma vectorial k · k∞ est´a dada por k A k∞ = m´ax

n X

1≤i≤n

|aij | = m´ax k fiT k1

j=1

1≤i≤n

donde fi es la i-´esima fila de A. (3) La norma matricial en Rn×n inducida por la norma vectorial k · k2 est´a dada por k A k1 = mayor valor singular de A = (mayor autovalor de A∗ A)1/2 . Demostraci´on. (1) Para cualquier x ∈ Rn se tiene k Ax k1 =k

n X

xj aj k1 ≤

j=1

n X

n X |xj | k aj k≤ ( |xj |) m´ax k aj k1 .

j=1

j=1

1≤j≤n

Pn on Si x es unitario (es decir, k x k1 = j=1 |xj | = 1) entonces la anterior expresi´ est´a acotada por m´ ax k aj k1 , (6.1) 1≤j≤n

lo que significa que m´ax k Ax k1 ≤ m´ax k aj k1 ,

kxk1 =1

1≤j≤n

(6.2)

es decir, la norma inducida est´a acotada por (6.1). Por otro lado, si j0 es el ´ındice de la columna de A en el que la norma-1 alcanza el m´aximo (es decir, k aj0 k1 = = m´ax1≤j≤n k aj k1 ), entonces k Aej0 k1 =k aj0 k1 = m´ax k aj k1 , 1≤j≤n

es decir, para elj0 -´esimo vector can´onico se alcanza la igualdad en (6.2). Por lo tanto, hemos llegado a que k A k1 = m´ax k Ax k1 = m´ax k aj k1 . kxk1 =1

1≤j≤n

(2) Para cualquier x ∈ Rn se tiene k Ax k∞ = m´ax {|(Ax)1 |, . . . , |(Ax)n |} = m´ax {|f1 x|, . . . , |fn x|} . Si x es un vector unitario, k x k∞ = m´ax{|x1 |, . . . , |xn |} = 1, es f´acil comprobar que |fi x| ≤k fiT k1 ,

∀i = 1, . . . , n.

De aqu´ı se sigue que m´ ax k Ax k∞ ≤ m´ax k fiT k1 .

kxk∞ =1

1≤i≤n

(6.3)

108

6. NORMAS DE VECTORES Y MATRICES. SENSIBILIDAD DE MATRICES Por otra parte, si i0 es el ´ındice de la columna de A en el que la norma-1 alcanza el m´aximo (es decir, k fiT0 k1 = m´ax1≤i≤n k fiT k1 ), y definimos el vector x0 cuya coordenada j-´esima es ½ 1 , si ai0 j > 0 (x0 )j = , −1 , si ai0 j < 0 entonces k x0 k∞ = 1 y k Ax0 k∞ = m´ax {|f1 x0 |, . . . , |fn x0 |} =k fiT0 k1 = m´ax k fiT k1 . 1≤i≤n

Por lo tanto, en (6.3) se alcanza la igualdad para x0 , lo que significa que k A k∞ = m´ax k Ax k∞ = m´ax k fiT k1 . 1≤i≤n

kxk∞ =1

(3) Sea A = U ΣV ∗ la SVD de A. Entonces k A k2 = m´ax k Ax k2 = m´ax k (U ΣV ∗ )x k2 = m´ax k Σ (V ∗ x) k2 | {z } kxk2 =1 kxk2 =1 kxk2 =1 y

  = m´ax k Σy k2 = m´ax k  kyk2 =1

kyk2 =1



σ1 y1 ..  k = m´ax p(σ y )2 + . . . + (σ y )2 , 1 1 n n .  2 kyk =1 2 σn yn

donde en la tercera y cuarta igualdades hemos usado, respectivamente, que U y V son unitarias. Como σ1 ≥ . . . ≥ σn ≥ 0, la u´ ltima expresi´on est´a acotada superiormente por σ1 m´ax k y k2 = σ1 . kyk2 =1

Por otro lado, la igualdad se da para x = v1 , donde v1 es el primer vector singular derecho, pues k Av1 k2 =k (U ΣV ∗ )v1 k2 =k U (Σe1 ) k2 =k Σe1 k2 = σ1 , con lo que el resultado queda demostrado. La norma matricial inducida por k · k1 es, por tanto, la mayor norma-1 de cualquier columna de A, mientras que la norma matricial inducida por k · k∞ es la mayor norma-1 de cualquier fila de A.

6.3 El radio espectral El radio espectral de una matriz A ∈ Rn×n , que denotaremos por ρ(A), es el mayor de los m´odulos de los autovalores de A, es decir ρ(A) = m´ax{|λ| : λ es autovalor de A}. La relaci´on del radio espectral con las normas matriciales se pone de manifiesto en el siguiente resultado.

6.4 Errores en inversas y en soluciones de sistemas lineales. N´umero de condici´on

109

Teorema 6.3.1 Sea A ∈ Rn×n . Entonces a) ρ(A) ≤k A k, para cualquier norma matricial k · k. b) Dado un valor ε > 0, existe una norma matricial k · k tal que ρ(A) ≤k A k≤ ρ(A) + ε. Lo que nos dice el teorema anterior es que el radio espectral est´a acotado por cualquier norma de A, pero hay normas que permiten aproximar el radio espectral de A tanto como se desee. Desgraciadamente, estas normas no son f´acilmente constru´ıbles, por lo que el Teorema 6.3.1 no es un resultado pr´actico para aproximar el radio espectral. La utilidad pr´actica del Teorema 6.3.1 se encuentra en el siguiente corolario, que usaremos en el Tema 7. Corolario 6.3.2 Sea A ∈ Rn×n . Entonces l´ımk→∞ Ak = 0 si y s´olo si ρ(A) < 1. La condici´on l´ımk→∞ Ak = 0 equivale a decir que todas las entradas de la k-´esima potencia de A tienden a cero.

6.4 Errores en inversas y en soluciones de sistemas lineales. ´ Numero de condicion ´ En c´alculos num´ericos es casi inevitable la presencia de errores. Por ejemplo, los c´alculos con un ordenador, que tiene una memoria finita y, por tanto, una “aritm´etica” finita, est´an sujetos a errores de redondeo y truncamiento. Por otra parte, en muchas ocasiones los datos de entrada sobre los que hay que efectuar los c´alculos no son los datos exactos, sino meras aproximaciones de e´ stos. Esta inexactitud puede obedecer a diversas causas, como la imprecisi´on en las herramientas de medici´on o errores en c´alculos previos. En estos casos es importante hacer un an´alisis de c´omo pueden afectar estos errores iniciales al resultado final. En los algoritmos m´as comunes los errores de redondeo se pueden tratar de la misma forma que los errores en los datos. Esta es la situaci´on que nosotros asumiremos, despreocup´andonos de los errores de redondeo y centr´andonos en los errores iniciales. Por tanto, podemos esquematizar la situaci´on de la siguiente manera: D + δD −→ S(D + δD) = S(D) + ? Con D representamos los datos del problema, que han sido alterados por una ligera perturbaci´on δD. El algoritmo num´erico que se emplea para resolver el problema nos devuelve la soluci´on (supuesto que no hay error en los c´alculos) del problema perturbado, S(D + δD). Lo que nos interesa estudiar es la desviaci´on que existe entre esa soluci´on y la soluci´on del problema con los datos exactos, S(D). Lo que haremos ser´a estimar esa diferencia mediante el error relativo k S(D + δD) − S(D) k e= . k S(D) k En este curso vamos a tratar el problema del c´alculo de la matriz inversa de una matriz A y el del c´alculo de la soluci´on de un sistema lineal (1.5). En este u´ ltimo problema supondremos que la matriz de coeficientes es invertible y, por tanto, hay soluci´on u´ nica. Ambos problemas

110

6. NORMAS DE VECTORES Y MATRICES. SENSIBILIDAD DE MATRICES

est´an ´ıntimamente relacionados, pues la soluci´on del sistema (1.5), supuesto que esta soluci´on es u´ nica, est´a dada por x = A−1 b. (6.4) Comenzaremos estudiando el problema de la matriz inversa y acto seguido aplicaremos los resultados obtenidos a la soluci´on del sistema lineal a partir de la relaci´on (6.4). Nuestro primer resultado es el siguiente, que proporciona una cota del error relativo. Teorema 6.4.1 (Cota para el error relativo en la inversa) Sean A, δA ∈ Rn×n , con A invertible, e = A + δA es invertible y se y k · k una norma matricial. Si α =k A−1 δA k< 1, entonces A tiene la siguiente cota del error relativo para la inversa e−1 k k A−1 − A α ≤ . −1 kA k 1−α

(6.5)

La condici´on α < 1 del teorema anterior no es f´acil de verificar si desconocemos la perturbaci´on δA. Desgraciadamente, dicha perturbaci´on no suele conocerse en la pr´actica. No obstante, podemos obtener una condici´on suficiente que garantiza la condici´on α < 1 y que est´a dada en t´erminos de cantidades que s´ı pueden conocerse en la pr´actica. Aunque no conozcamos expl´ıcitamente la matriz δA, s´ı es posible, en muchas ocasiones, conocer su norma o, al menos, una cota superior de su norma. Podemos sacar provecho de esta informaci´on razonando de la siguiente manera. Recordemos que, por el axioma iv) de las normas matriciales, k A−1 δA k≤k A−1 k · k δA k . Esta desigualdad tiene, al menos, las siguientes tres consecuencias 1) Si k A−1 k · k δA k< 1 entonces k A−1 δA k< 1. Por lo tanto, la condici´on α < 1 del enunciado del Teorema 6.4.1 se satisface si la norma de la perturbaci´on δA satisface la acotaci´on 1 k δA k< . (6.6) k A−1 k 2) La condici´on (6.6) junto con el Teorema 6.4.1 nos dicen que si δA es suficientemente peque˜na (al menos tanto como para verificar (6.6)) entonces la matriz perturbada e = A + δA sigue siendo invertible. En consecuencia, para que una matriz deje de A ser invertible ha de ser modificada por una perturbaci´on suficientemente grande. 3) La cantidad α/(1 − α) en (6.5) puede acotarse por k A−1 k · k δA k α ≤ , 1−α 1− k A−1 k · k δA k donde para calcular el cociente de la derecha no es necesario conocer expl´ıcitamente δA, sino u´ nicamente su norma (o una cota superior de ella). La cota superior que aparece en 3) permite obtener, de manera inmediata a partir de (6.5), una cota superior para el error relativo de la inversa que no requiere el conocimiento expl´ıcito de δA (pero s´ı de su norma). No obstante, enunciaremos este resultado utilizando el concepto de n´umero de condici´on de A, lo que nos dar´a una cota que permite relacionar el error en las soluciones (inversas) con el error en los datos (matrices A, δA).

6.4 Errores en inversas y en soluciones de sistemas lineales. N´umero de condici´on

111

Definici´on 6.4.2 Dada una matriz invertible A ∈ Rn×n y k · k una norma matricial, el n´umero de condici´on de A es el n´umero real no negativo c(A) =k A k · k A−1 k . El n´umero de condici´on depende de la norma. As´ı, tendremos, por ejemplo, el n´umero de condici´on asociado a la norma-2, c2 (A), o a la norma infinito, c∞ (A), etc. El caso de la norma-2 es particularmente relevante. Lo podemos expresar como σ

(A)

c2 (A) = σm´ax(A)

(6.7)

m´ın

donde σm´ax (A) y σm´ın (A) son, respectivamente, el mayor y el menor valor singular de A (ambos son no nulos porque A es invertible). La demostraci´on de dicha f´ormula pasa por resolver el siguiente ejercicio. Ejercicio 6.4.3 Demostrar que, si A es invertible, entonces k A−1 k2 = 1/σm´ın (A). Si A es una matriz normal entonces, como sabemos del Corolario 3.1.7, los valores singulares coinciden con los m´odulos de los autovalores. En este caso Si A es normal:

|λ

(A)|

c2 (A) = |λm´ax(A)| m´ın

Como hemos visto en el Ejercicio 6.4.3 la norma-2 de la inversa coincide con el inverso del menor valor singular de A. Eso quiere decir que si dicho valor singular es peque˜no (pr´oximo a 0) entonces el n´umero de condici´on de A es elevado. Una matriz con un n´umero de condici´on elevado es una matriz que se “comporta mal” respecto a la inversi´on, es decir, que peque˜nas alteraciones en las entradas de la matriz pueden provocar grandes alteraciones en su inversa o, incluso, hacerla no invertible. Ahora podemos enunciar el siguiente resultado, que es consecuencia inmediata del Teorema 6.4.1, pero tiene las dos ventajas que se han mencionado arriba: por un lado, la cota del error relativo para las inversas est´a dada en t´erminos que no requieren el conocimiento expl´ıcito de la perturbaci´on y, por otro, establece una relaci´on entre los errores en los datos y los errores en las inversas. Dicha relaci´on involucra al n´umero de condici´on en los t´erminos que explicaremos tras enunciar el resultado. Corolario 6.4.4 Sean A, δA ∈ Rn×n , k · k una norma matricial y c(A) el n´umero de cone = A + δA es dici´on de A en esa norma. Si k δA k< 1/ k A−1 k la matriz perturbada A invertible y e−1 k k A−1 − A c(A) k δA k ≤ . (6.8) −1 kδAk kA k kAk 1 − c(A) kAk

Si k δA k es peque˜na, el denominador de la expresi´on derecha en (6.8) es aproximadamente 1, y el t´ermino completo es aproximadamente igual a c(A) ·

k δA k . kAk

6. NORMAS DE VECTORES Y MATRICES. SENSIBILIDAD DE MATRICES

112

Esto significa que el error relativo en las inversas est´a acotado por el error en los datos a trav´es de un factor igual al n´umero de condici´on de la matriz de partida A. Por lo tanto, si el n´umero de condici´on de A es peque˜no, el error en la soluci´on (inversa) es del mismo orden que el error en los datos. En cambio, si el n´umero de condici´on es elevado, peque˜nos errores en los datos pueden provocar grandes errores en la inversa. Esto nos lleva a decir que una matriz est´a bien condicionada si su n´umero de condici´on es peque˜no (pr´oximo a 1), y que est´a mal condicionada si el n´umero de condici´on es elevado. Si c(A) = 1 se acostumbra a decir que A est´a perfectamente condicionada (con respecto a la norma matricial asociada). Ejercicio 6.4.5 Demostrar que c(A) ≥ 1 para cualquier A ∈ Rn×n y cualquier norma matricial en Rn×n . Problema Resuelto 6.4.6 Sea la matriz · A² =

1 −1 −1 1 + ²

¸

donde ² > 0 es un n´umero real arbitrario. a) Calcular c2 (A² ). b) Estudiar el condicionamiento de A² para valores peque˜nos de ². Soluci´on. a) La matriz A² es normal, de modo que c2 (A² ) =

|λm´ax (A² )| . |λm´ın (A² )|

Los autovalores de A² se obtienen f´acilmente a partir del polinomio caracter´ıstico √ √ 2 + ² + ²2 + 4 2 + ² − ²2 + 4 λ1 = , λ2 = . 2 2 Entonces

à !2 √ √ 2 + ² + ²2 + 4 1 2 + ² + ²2 + 4 √ c2 (A² ) = = . ² 2 2 + ² − ²2 + 4

b) El factor de la derecha en la expresi´on final de c2 (A² ) est´a acotado para valores suficientemente peque˜nos de ². Por lo tanto c2 (A² ) es de orden 1/². Eso significa que para valores peque˜nos de ² el n´umero de condici´on de A² es grande (del orden de 1/²) y, por tanto, la matriz A² est´a mal condicionada. De hecho, la perturbaci´on · ¸ 0 0 δA = , 0 −² cuya norma-2 es igual a ², convierte a la matriz A en una matriz singular A + δA.

6.4 Errores en inversas y en soluciones de sistemas lineales. N´umero de condici´on

113

Si calculamos el determinante de la matriz A² en el Problema Resuelto 6.4.6 obtendremos inmediatamente det(A² ) = ². Por tanto, para valores peque˜nos de ² el determinante es peque˜no. Esto nos puede llevar a pensar que el determinante es una medida del condicionamiento de una matriz. Esta idea es equivocada. El condicionamiento es una estimaci´on de la proximidad de una matriz a la singularidad en t´erminos relativos. El determinante, en cambio, mide la proximidad a la singularidad en t´erminos absolutos. En otras palabras: el determinante puede ser peque˜no porque las entradas de la matriz son peque˜nas. En este caso, la matriz no est´a necesariamente mal condicionada. Por otro lado, el hecho de que una matriz invertible est´e mal condicionada nos indica que hay perturbaciones “malas” , pero no que cualquier perturbaci´on sea “mala” . Podemos detectar de antemano algunas perturbaciones malas si conocemos la SVD de la matriz. Para comprender esto mejor, comencemos con el caso elemental de una matriz que es diagonal (invertible) y, por simplicidad, con autovalores positivos. En este caso, los valores singulares de A son las entradas diagonales, es decir A = diag(σ1 , . . . , σn ), con σ1 ≥ . . . ≥ σn . Si el menor valor singular σn es muy peque˜no, una perturbaci´on mala es δA = diag(0, . . . , 0, −σn ), cuya norma-2 es igual a σn (peque˜na, por tanto) y A + δA es singular. En otras palabras: la matriz A es muy sensible a modificaciones en la u´ ltima entrada diagonal. Este razonamiento puede generalizarse a matrices no necesariamente diagonales a partir de la SVD. Concretamente, si A = U ΣV ∗ es la SVD de A, y el valor singular σi es peque˜no, una perturbaci´on “mala” ser´a δA = U diag(0, . . . , 0, −σi , 0, . . . , 0)V ∗ , donde el valor −σi est´a en la entrada (i, i). La matriz perturbada A + δA es singular (hemos “eliminado” el i-´esimo valor singular) y la norma-2 de la perturbaci´on es igual a σi (peque˜na, por tanto). Este procedimiento, que consiste en localizar perturbaciones “malas” de rango bajo de A, est´a ´ıntimamente relacionado con lo que vimos en la Secci´on 3.3 acerca de las aproximaciones de rango bajo de una matriz A. Consideremos ahora la soluci´on del sistema lineal (1.5). Si la matriz A de coeficientes del e = A + δA, el sistema perturbado es sistema se perturba por una matriz δA y denotamos por A ex = b, Ae

(6.9)

donde hemos usado una notaci´on distinta para la soluci´on de (6.9) para evitar confundirla con la soluci´on de (1.5). Como consecuencia inmediata del Teorema 6.4.1 y del Corolario 6.4.4, tenemos Teorema 6.4.7 Sean A ∈ Rn×n una matriz invertible y k · k una norma vectorial en Rn y su norma matricial inducida en Rn×n . Sea δA ∈ Rn×n una matriz tal que δA < 1/ k A−1 k y e = A + δA. Sean x, x ex = b. denotemos A e las soluciones de los sistemas lineales Ax = b y Ae Entonces kx−x ek c(A) k δA k ≤ . (6.10) kδAk kxk kAk 1 − c(A) kAk

Igual que en el caso de la f´ormula (6.8), la f´ormula (6.10) nos dice que si A est´a bien condicionada, los errores en la soluci´on del sistema (1.5) no difieren mucho de los errores en los datos (matriz de coeficientes). Si se perturba tambi´en el t´ermino independiente, tenemos el siguiente resultado, que proporciona una cota para el error relativo de las soluciones.

6. NORMAS DE VECTORES Y MATRICES. SENSIBILIDAD DE MATRICES

114

Teorema 6.4.8 Sean A ∈ Rn×n una matriz invertible y k · k una norma vectorial en Rn y la norma matricial inducida en Rn×n . Sea δA ∈ Rn×n una matriz tal que δA < 1/ k A−1 k y e = A + δA, eb = b + δb. Si x, x b, δb dos vectores de Rn (con b 6= 0) y denotemos A e son las ex = eb, entonces soluciones de los sistemas lineales Ax = b y Ae µ ¶ kx−x ek c(A) k δA k k δb k ≤ + . (6.11) kxk kbk 1 − c(A) kδAk k A k kAk

Como se ve, la cota del error relativo en (6.11) consta de dos t´erminos. El primero de ellos corresponde al error en la matriz de coeficientes A y el segundo de ellos al error en el t´ermino independiente b. Problema Resuelto 6.4.9 Fijemos ² = 10−4 en la matriz A² del Problema Resuelto 6.4.6, llamemos A = A10−4 y consideremos el sistema lienal Ax = [1 1]T . a) Si δA es una perturbaci´on de A con k δA k2 = 10−4 , obtener una cota te´orica del error cometido en las soluciones del sistema. b) Obtener una perturbaci´on particular para la cual el error en las soluciones sea aproximadamente igual a la cota del apartado anterior. Soluci´on. a) Vamos a calcular la cota te´orica (6.10) establecida en el Teorema 6.4.7. Los valores que aparecen en dicha cota son los siguientes c2 (A) = 4002 , k A k2 = |λm´ax (A)| = 2.0005 , por lo que

k δA k2 = 10−4 ,

kx−x e k2 4002 10−4 ≤ · = 0.2500781494 . k x k2 0.79995 2.0005

Esto significa que una perturbaci´on de una diezmil´esima en la matriz de coeficientes del sistema puede provocar un error de m´as de dos d´ecimas en la soluci´on. Naturalmente que esto es s´olo una cota, pero en el siguiente apartado veremos que, de hecho, existe una perturbaci´on para la que esta cota se alcanza. b) La soluci´on del sistema no perturbado Ax = [1 1]T es · ¸ · ¸ x1 2001 = . x2 2000 Si

· δA =

0.0001 0 0 0.0001

¸

entonces k δA k2 = 10−4 y la soluci´on del sistema perturbado (A + δA)e x = [1 1]T es · ¸ · ¸ x e1 1667.430 = . x e2 1666.597

6.4 Errores en inversas y en soluciones de sistemas lineales. N´umero de condici´on

115

El error relativo cometido es kx−x e k2 = 0.16670 . k x k2 No obstante, podemos encontrar una perturbaci´on para la cual el error cometido es pr´acticamente igual a la cota obtenida en el apartado anterior. M´as a´un, podemos obtener una perturbaci´on de rango 1. Para ello, seguimos el procedimiento que hemos descrito en la p´agina 113. Para ello necesitamos en primer lugar conocer la SVD de A. El programa Matlab nos la ofrece con poco esfuerzo (el de teclear unas cuantas veces). A = U ΣV ∗ , donde · ¸ · ¸ −0.7069299824 0.7072835358 2.00050 0 U= , Σ= , 0.7072835358 0.7069299824 0 0.000499875 · ¸ −0.7069299824 0.7072835358 V = . 0.7072835358 0.7069299824 Tomemos la siguiente perturbaci´on de norma 10−4 y rango 1: · ¸ · ¸ 0 0 −0.5002500000 −0.4999999375 ∗ −4 δA = U V = 10 · . 0 −10−4 −0.4999999375 −0.4997500000 Para esta perturbaci´on la soluci´on del sistema (A + δA)e x = b es · ¸ · ¸ x e1 2501.4064082275 = , x e2 2500.1562675742 y el error relativo cometido es kx−x e k2 = 0.2502188712 . k x k2 Aunque la perturbaci´on anterior es muy espec´ıfica, la mayor´ıa de las perturbaciones implican un error relativo en la soluci´on de aproximadamente el orden del n´umero de condici´on de A. En la siguiente tabla mostramos los errores relativos en las soluciones de los sistemas perturbados para 10 perturbaciones aleatorias de norma 10−4 . 1 2 3 4 5 6 7 8 9 10

Error relativo 0.1755514003 0.2109933108 0.1831942131 0.2399546841 0.1098550106 0.1026737324 0.1097760122 0.1844815833 0.1834607827 0.1610196936

Como se puede ver, en todos los casos el error es del orden de entre 1000 y 2000 veces m´as grande que el tama˜no de la perturbaci´on.

116

6. NORMAS DE VECTORES Y MATRICES. SENSIBILIDAD DE MATRICES

7 ´ A LOS METODOS ´ INTRODUCCION ITERATIVOS

Vamos a estudiar dos m´etodos iterativos que se aplican a la resoluci´on de sistemas de ecuaciones lineales: el m´etodo de Jacobi y el m´etodo de Gauss-Seidel. En la Secci´on 7.1, veremos por encima en qu´e consisten los m´etodos iterativos en general y, posteriormente, en la Secci´on 7.2 veremos c´omo se aplican a la resoluci´on de sistemas lineales y estudiaremos en particular los dos m´etodos que acabamos de mencionar.

´ 7.1 Metodos iterativos Dada una matriz B ∈ Cn×n y un vector c ∈ Cn , generamos la sucesi´on x0 ∈ Cn ,

xk+1 = Bxk + c .

(7.1)

Nos referiremos a dicha sucesi´on como m´etodo iterativo. Diremos que el m´etodo iterativo (7.1) es convergente si existe l´ımk→∞ xk y el l´ımite, que denotaremos por x, no depende del valor inicial x0 . Nuestra intenci´on es aplicar un m´etodo iterativo a la resoluci´on del sistema lineal Ax = b ,

(7.2)

donde A ∈ Cn×n y b ∈ Cn . En lo sucesivo, supondremos que A es invertible, con lo cual el sistema tiene soluci´on u´ nica. Para aplicar el m´etodo (7.1) al sistema debemos asegurarnos de que el l´ımite, caso de existir, ha de ser la soluci´on de dicho sistema. Para que esto ocurra, si llamamos x al l´ımite de la sucesi´on (7.1) y a la soluci´on del sistema (7.2), entonces, por un lado, x = A−1 b y, por otro, tomando lmites en (7.1), ha de ser x = Bx + c. Igualando, llegamos a la siguiente definici´on. 117

118

´ A LOS METODOS ´ 7. INTRODUCCION ITERATIVOS

Definici´on 7.1.1 Diremos que el m´etodo (7.1) es consistente con el sistema (7.2) si c = (I − B)A−1 b. El siguiente resultado nos da una caracterizaci´on de los m´etodos iterativos convergentes que son consistentes con el sistema (7.2). Siguiendo con la notaci´on del tema anterior, ρ(B) denotar´a el radio espectral de B. Teorema 7.1.2 Si (7.1) es un m´etodo consistente con el sistema (7.2) entonces los siguientes enunciados son equivalentes. i) El m´etodo iterativo (7.1) es convergente. ii) ρ(B) < 1. iii) Existe una norma matricial k · k tal que k B k< 1. Demostraci´on. La equivalencia entre los apartados ii) y iii) es consecuencia inmediata del Teorema 6.3.1. Vamos a demostrar la equivalencia entre i) y ii). Sea x la soluci´on del sistema (7.2). Por la hip´otesis de consistencia tenemos que B(xk − x) = Bxk − Bx = xk+1 − c − BA−1 b = xk+1 − (I − B)A−1 b − BA−1 b = xk+1 − x. Iterando la identidad anterior llegamos a xk+1 − x = B k (x0 − x).

(7.3)

En primer lugar, supongamos que ρ(B) < 1. Como ρ(B) < 1, el Corolario 6.3.2 nos dice que l´ımk→∞ B k = 0, luego l´ım xk+1 = x, es decir, el m´etodo es convergente a x. Rec´ıprocamente, supongamos que ρ(B) > 1. Entonces, existe al menos un autovalor de B, λ0 , de m´odulo mayor que 1. Si v0 es el autovector asociado a λ0 , tenemos que Bv0 = λ0 v0 . Si tomamos x0 = v0 + x como vector inicial del m´etodo iterativo, se tiene que, usando (7.3) xk+1 − x = B k (x0 − x) = B k v0 = λk v0 . Como |λ| > 1, λk no converge a 0, luego el vector λk v0 tampoco converge al vector 0. Eso significa, por la u´ ltima igualdad, que el m´etodo (7.1) no converge a x. En caso de que el m´etodo iterativo (7.1) sea convergente, el k-´esimo t´ermino de la sucesi´on, xk , es una aproximaci´on al l´ımite, x. El siguiente resultado nos da una cota del error que cometemos al aproximar x por xk+1 . Proposici´on 7.1.3 Sea (7.1) un m´etodo consistente con el sistema (7.2) y k · k una norma vectorial y su correspondiente norma matricial inducida. Entonces, si k B k< 1, el m´etodo iterativo (7.1) es convergente a un vector x con la siguiente cota de error k xk+1 − x k≤

k B kk+1 k x1 − x0 k . 1− k B k

(7.4)

7.2 Aplicaci´on a los sistemas lineales. M´etodos de Jacobi y Gauss-Seidel

119

´ a los sistemas lineales. Metodos ´ 7.2 Aplicacion de Jacobi y GaussSeidel Ahora vamos a ver c´omo podemos aplicar los m´etodos iterativos, definidos en la secci´on anterior, a la resoluci´on de sistemas lineales. Para aplicar los m´etodos anteriores a la resoluci´on del sistema (7.2) haremos lo siguiente: 1. Descomponer A = M − N , donde M es una matriz invertible y cuya inversa sea f´acil de obtener (por ejemplo, una matriz diagonal, triangular...). 2. El sistema (7.2) es equivalente a (M − N )x = b ⇔ M x = N x + b ⇔ x = M −1 N x + M −1 b = Bx + c , donde B = M −1 N , c = M −1 b. N´otese que, si B es la matriz reci´en definida, entonces I − B = I − M −1 N = M −1 (M − N ) = M −1 A, que es invertible por serlo A. Adem´as (I − B)A−1 b = M −1 b = c, luego el m´etodo iterativo (7.1) es consistente con el sistema (7.2). Tambi´en es importante notar que, para una matriz concreta, A, hay diferentes elecciones de la matriz B. Los m´etodos de Jacobi y Gauss-Seidel se obtienen precisamente a partir de dos elecciones distintas de la matriz B. Ambos son v´alidos para una matriz A con entradas diagonales no nulas y se basan en la siguiente descomposici´on de la matriz A: A=D−L−U, donde D es una matriz diagonal con entradas diagonales no nulas (por tanto, invertible), L es una matriz triangular inferior con entradas diagonal nulas y U es una matriz triangular superior con entradas diagonales nulas. Naturalmente que una descomposici´on as´ı es siempre posible (separando la parte diagonal de la parte superior y de la parte inferior de A) si que A tiene entradas diagonales no nulas. M´etodo de Jacobi: Consiste en tomar M = D, N = L + U , es decir: BJ = D−1 (L + U ) = In − D−1 A. M´etodo de Gauss-Seidel: Consiste en tomar M = D − L, N = U , es decir: BGS = (D − L)−1 U = In − (D − L)−1 A. Los m´etodos anteriores son utilizados frecuentemente en la resoluci´on de sistemas lineales con muchas ecuaciones. Recu´erdese que la matriz A debe tener entradas diagonales no nulas, pues de lo contrario la matriz D no ser´ıa invertible. La pregunta que debemos formularnos ahora es: ¿Para qu´e matrices son convergentes los m´etodos de Jacobi y Gauss-Seidel? Desgraciadamente, no existe una respuesta que abarque a cualquier matriz A (con entradas diagonales no nulas). Por supuesto que, para una matriz dada, siempre se puede usar la Proposici´on 7.1.3. Para ello hay que calcular el radio espectral de la matriz B asociada al m´etodo. No obstante, hay un tipo de matrices para las cuales esto no es

´ A LOS METODOS ´ 7. INTRODUCCION ITERATIVOS

120

necesario porque podemos asegurar que los m´etodos son convergentes. Estas matrices son un subconjunto particular de las matrices tridiagonales. Diremos que una matriz A es tridiagonal si A(i, j) = 0 para |i − j| > 1. En otras palabras, las matrices tridiagonales son aquellas que todas las entradas que no est´an ni en la diagonal principal ni en las dos diagonales adyacentes (la superior y la inferior) son nulas. Proposici´on 7.2.1 Si A es una matriz tridiagonal herm´ıtica definida positiva, entonces los m´etodos de Jacobi y Gauss-Seidel para A son convergentes. Diremos que una matriz A ∈ Cn×n es diagonalmente dominante por filas si n X

|aii | ≥

|aij |,

∀i = 1, . . . , n

j =1 j 6= i

y que A es estrictamente diagonalmente dominante por filas si |aii | >

n X

|aij |,

∀i = 1, . . . , n.

j =1 j 6= i

En otras palaras, A es diagonalmente dominante si el m´odulo de la entrada diagonal de cada fila es mayor o igual que la suma de los m´odulos de las restantes entradas de la fila. Proposici´on 7.2.2 Si A es estrictamente diagonalmente dominante por filas entonces los m´etodos de Jacobi y Gauss-Seidel para A son convergentes. Problema Resuelto 7.2.3 1) Para cada una de las matrices Ai siguientes estudiar si los m´etodos de Jacobi y Gauss-Seidel convergen a la soluci´on del sistema lineal Ai x = b, donde b ∈ R3 es un vector cualquiera. 

 4 1 2 i) A1 =  −1 2 0  1 1 3   5 2 0 ii) A2 =  2 1 1  0 1 6   1 −1 −1 1 −2  iii) A3 =  2 −1 2 1 £ ¤T 2) Tomemos b = 4 2 3 . Para cada uno de los m´etodos convergentes, efect´ua dos iteraciones en cada uno de los dos m´etodos y estima el n´umero de iteraciones necesario para aproximar la soluci´on con un error menor que una mil´esima partiendo del vector x0 = 0R3 . Soluci´on.

7.2 Aplicaci´on a los sistemas lineales. M´etodos de Jacobi y Gauss-Seidel 1)

121

i) La matriz A1 es estrictamente diagonalmente dominante por filas, pues 4 > 1 + 2, 2 > | − 1| + 0, 3 > 1 + 1, luego, por la Proposici´on 7.2.1, los m´etodos de Jacobi y Gauss-Seidel asociados a la matriz A1 son convergentes. ii) La matriz A2 es tridiagonal y herm´ıtica (sim´etrica). Como los menores principales son positivos ¯ ¯ ¯ ¯ ¯ 5 2 0 ¯ ¯ 5 2 ¯ ¯ ¯ ¯ = 1 > 0, ¯ 2 1 1 ¯ = 1 > 0 5 > 0, ¯¯ ¯ ¯ ¯ 2 1 ¯ 0 1 6 ¯ sabemos (Teorema 4.3.5) que tambi´en es definida positiva. Por la Proposici´on 7.2.2 los m´etodos de Jacobi y Gauss-Seidel asociados a la matriz A2 son convergentes. iii) Como

     0 0 0 0 1 1 1 0 0 A =  0 1 0  −  −2 0 0  −  0 0 2 , 1 −2 0 0 0 0 0 0 1 {z } | {z } | {z } | 

D

L

U

la matriz del m´etodo de Jacobi es 

 0 1 1 BJ = D−1 (L + U ) = L + U =  −2 0 2  . 1 −2 0 Esta matriz tiene un autovalor igual a 1, luego ρ(BJ ) ≥ 1 y, por el Teorema 7.1.2 el m´etodo no es convergente. La matriz de Gauss-Seidel es   0 1 1 BGS = (D − L)−1 U =  0 −2 0  , 0 5 1 cuyos autovalores son λ1 = 0, λ2 = 1, λ3 = −2, por lo que ρ(BGS ) = 2 y el m´etodo no es convergente (Teorema 7.1.2). 2) La matriz de Jacobi asociada a la matriz A1 es 

 0 −1/4 −1/2 0 0  BJ = D−1 (L + U ) =  1/2 −1/3 −1/3 0

Por otro lado,

 1 c = D−1 b =  1  . 1 

´ A LOS METODOS ´ 7. INTRODUCCION ITERATIVOS

122 Tomando x0 = 0 tenemos



 1 x1 = BJ x0 + c = c =  1  , 1



 1/4 x2 = BJ x1 + c =  3/2  . 1/3

La matriz de Gauss-Seidel para la matriz A1 es  BGS

 0 −1/4 −1/2 = (D − L)−1 U =  0 −1/8 −1/4  0 1/8 1/4

y



 1 c = (D − L)−1 b =  3/2  1/6

Entonces, para x0 = 0, 

 1 x1 = BGS x0 + c = c =  1  , 1



 13/24 x2 = BGS x1 + c =  61/48  . 19/48

Para estimar el error usaremos la cota (7.4). Como, por el Teorema 6.3.1 hay una norma tan pr´oxima como se desee al radio espectral, podemos tomar ρ(B) en lugar de k B k en dicha f´ormula. Para el m´etodo de Jacobi, primero calculamos con Matlab la norma-2 de la matriz BJ usando el comando norm . El resultado es k BJ k2 = 0.65758 Sustituyendo en la cota (7.4) obtenemos k xk+1 − x k2 ≤

0.65758k+1 √ 3 = 5.05829 · 0.65758k+1 . 1 − 0.65758

Para encontrar el primer valor de k para el que la cota anterior es menor que una mil´esima resolvemos la inecuaci´on en k : 5.05829·0.65758k+1 ≤ 0.001 ⇒ k+1 ≥ log (0.001/5.05829) / log(0.65758) = 20.34604. As´ı pues, podemos asegurar que x20 se aproxima a la soluci´on del sistema con un error menor que una mil´esima. Para el m´etodo de Gauss-Seidel procedemos como antes, calculando primero con Matlab la norma-2 de BGS . El resultado es k BGS k2 = 0.68465. Ahora k xk+1 − x k2 ≤

0.68465k+1 · 1.81046 = 5.74118 · 0.68465k+1 . 1 − 0.68465

De la misma forma, resolvemos 5.74118·0.68465k+1 ≤ 0.001 ⇒ k+1 ≥ log (0.001/5.74118) / log(0.68465) = 22.84671,

7.2 Aplicaci´on a los sistemas lineales. M´etodos de Jacobi y Gauss-Seidel

123

lo que nos asegura que x22 es una aproximaci´on a la soluci´on con un error menor que una mil´esima. En realidad, estas estimaciones no son ajustadas. De hecho, para el m´etodo de Jacobi x9 ya aproxima la soluci´on con un error menor que una mil´esima, mientras que en el m´etodo de Gauss-Seidel esto ocurre con x5 . El apartado ii) se resuelve de manera an´aloga.