Cap´ıtulo 6 Proyecciones y bases ortonormales
6.1.
Introducci´ on
Este es el primero de los dos cap´ıtulos dedicados a la factorizaci´on QR y a la soluci´on del problema de m´ınimos cuadrados. La factorizaci´on QR es una manera matricial de presentar el ya conocido procedimiento de ortonormalizaci´on de Gram-Schmidt. No perseguimos volver a presentar este procedimiento sino mostrar que hay varias formas de conseguir bases ortonormales que producen algoritmos con distintas propiedades. En este cap´ıtulo nos centraremos en algoritmos ligados al m´etodo de Gram-Schmidt y en el siguiente estudiaremos las reflexiones de Householder. Veremos que los principios que fundamentan ambos tipos de algoritmos son completamente diferentes. Debemos recordar que, b´asicamente, el m´etodo de ortonormalizaci´on de GramSchmidt consiste en lo siguiente: Dado un subespacio del que se conoce una base ortonormal y un vector que no est´a en el subespacio, se proyecta ´este sobre el ortogonal de aqu´el a fin de conseguir una base ortonormal para un subespacio de dimensi´on una unidad mayor. Procediendo recursivamente desde un vector dado se puede conseguir una base ortonormal de todo el espacio. Todo este procedimiento 121
122
Proyecciones y bases ortonormales
est´a basado en el concepto de proyecci´on ortogonal que repasamos brevemente a continuaci´on, para volver luego sobre dos algoritmos diferentes ligados, ambos, al m´etodo de Gram-Schmidt.
6.2.
Proyecciones
Comenzamos definiendo lo que entendemos por proyecci´on: Definici´ on 6.1 Una matriz P P Cnn se dice que es una proyecci´ on si es una 2 matriz idempotente. Es decir, si P P . Algunas consecuencias inmediatas de esta definici´on son las siguientes Proposici´ on 6.2 Si P P Cnn es una proyecci´on tambi´en lo es In P . A ´esta se le llama proyecci´ on complemetaria de P y cumple que Ker P ImpIn P q y KerpIn P q Im P . Demostraci´ on.- En primer lugar pIn P q2 P P . As´ı pues, In P es una proyecci´on. 2
In 2P
P2
In P
porque
Adem´as, y P ImpIn P q ô pIn P qx y para alg´ un x. Y de aqu´ı, P y P pIn P qx pP P 2 qx pP P qx 0, lo que indica que y P Ker P . Rec´ıprocamente, y P Ker P ñ P y 0 ñ pIn P qy y ñ y P ImpIn P q. La identidad KerpIn P q Im P se demuestra cambiando los papeles de In P y P.
Proposici´ on 6.3 Si P
P Cnn es una proyecci´on entonces Im P ` Ker P Cn
Demostraci´ on.- En primer lugar, la suma Im P Ker P es directa porque, de acuerdo con la Proposici´on anterior, Im P KerpIn P q y si x P KerpIn P qX Ker P entonces pIn P qx 0 y P x 0. De aqu´ı, x 0.
6.3 Proyecciones ortogonales
123
Por otra parte Im P Ker P Cn . Y si x P Cn es un vector cualquiera entonces x P x pIn P qx. Pero P x P Im P y pIn P qx P ImpIn P q Ker P . Esta simple Proposici´on tiene una consecuencia importante: si P es una proyecci´on, ´esta determina dos subespacios mutuamente suplementarios: Im P y Ker P . El rec´ıproco tambi´en es cierto. Si S1 y S2 son subespacios de Cn tales que S1 ` S2 Cn entonces hay una proyecci´on P tal que S1 Im P y S2 Ker P . En efecto, siendo S1 ` S2 Cn , cualquier matriz queda completamente determinada por su acci´on sobre los vectores de S1 y S2 . En concreto, si definimos P P Cnn como la matriz que cumple " x si x P S1 Px 0 si x P S2 resulta que en cualquier base de Cn obtenida como reuni´on de bases de S1 y S2 , la matriz de P , como aplicaci´on lineal, es P
Ip 0 0 0
siendo p dim S1 . En cualquier otra base de Cn , la matriz ser´ıa semejante a ´esta. Es claro que P es una proyecci´on y que Im P S1 y Ker P S2 . Es la proyecci´on sobre S1 a lo largo de (o paralelamente a) S2 . As´ı pues Conclusi´ on: Si P es una proyecci´on, proyecta Cn sobre Im P a lo largo de Ker P .
6.3.
Proyecciones ortogonales
Como su nombre indica las proyecciones ortogonales son las que proyectan a lo largo del ortogonal. Definici´ on 6.4 Una proyecci´on P
P Cnn es ortogonal si Ker P pIm P qK.
Es f´acil demostrar que para cualquier matriz A P Cmn , pIm AqK Ker A . En efecto, x P Ker A ô A x 0 ô x A 0 ô x Ay 0 para todo y P Cn ô x K z para todo z P Im A ô x P pIm AqK .
124
Proyecciones y bases ortonormales
Proposici´ on 6.5 Una proyecci´on P es ortogonal si y s´olo si P herm´ıtica).
P (i.e., es
Demostraci´ on.- En primer lugar, si P es proyecci´on tambi´en P lo es. Ahora, si P es una proyecci´on ortogonal entonces Ker P pIm P qK KerpP q. Tambi´en se cumple que Im P ImpP q. As´ı es,
pIm P qK KerpP q Ker P pIm P qK,
por definici´on de proyecci´on ortogonal. Como para cualquier subespacio de Cn , pS KqK S, conlcu´ımos que Im P Im P . Ahora bien, dos aplicaciones lineales que coinciden sobre dos subespacios suplementarios son la misma. Como P x P x x para x P Im P ImpP q y P x P x 0 para x P Ker P KerpP q, tenemos que P P . Resultar´a m´as u ´til la siguiente caracterizaci´on de las proyecciones ortogonales. Proposici´ on 6.6 Una matriz P P Cnn de rango r es una proyecci´on ortogonal si y s´olo si existe una matriz Q P Cnr , con columnas ortonormales, tal que P QQ . Demostraci´ on.- Supongamos primero que P nr C con columnas ortonormales. Entonces P2
QQ para alguna matriz Q P
QQQQ QpQQqQ QQ,
porque Q Q Ir por tener Q las columnas ortonormales. Por lo tanto P 2 decir, P es una proyecci´on. Adem´as P
P ; es
pQQq pQqQ QQ P.
As´ı que, por la Porposici´on anterior P es una proyecci´on ortogonal. Supongamos ahora que P es una proyecci´ on ortogonal y sea tq1 , . . . , qr u una base ortonormal de Im P . Pongamos Q q1 q2 qr . Veamos que P QQ . Para ello probaremos que P x QQ x para todo x P Cnn . Sea x un vector cualquiera de Cn . Como P es una proyecci´on Cn Ker P ` Im P . Por lo tanto existen unos u ´nicos vectores x1 P Im P y x2 P Ker P tales que x x1 x2 . Entonces P x P x1 x1 y P x2 0. Es decir, P x x1 . Calculemos QQ x.
6.3 Proyecciones ortogonales
Por una parte
125
q1 x r q x ¸ 2 QQ x Q .. pqi xqqi . . qr x
Ahora bien, como qi
P Im P
y x2
P Ker P pIm P qK
qi x qi x1 Pero, Im P q1 , . . . , qr yendo en (6.1)
(6.1)
i 1
qi x2
qix1.
¡, por lo que x1 a1q1
ar qr y qi x1
ai. Sustitu-
r ¸ QQ x ai q i x 1 . j 1
En conclusi´on P x QQ x para todo x P Cn y P
QQ.
Observaciones 6.7 1. Si Q tiene columnas ortonormales entonces P QQ es la proyecci´on ortogonal sobre Im P a lo largo de Ker P pIm P qK Ker P . Ahora bien, como P QQ y Q es de rango completo, Im P Im Q y Ker P Ker Q . En conclusi´on QQ es la proyecci´on ortogonal sobre Im Q a lo largo de Ker Q pIm QqK . 2. In QQ tambi´en es una proyecci´on, la complementaria de QQ . Y tambi´en es ortogonal porque pIn QQ q In QQ . Por la Porposici´on 6.2, ImpIn QQ q KerpQQ q Ker Q pIm QqK y KerppIn QQ q ImpQQ q Im Q. Por lo tanto, es la proyecci´on ortogonal sobre pIm QqK a lo largo de Im Q. 3. Nos interesan muy especialmente las proyecciones ortogonales de rango 1. Son las de la forma Pq qq siendo q un vector unitario. Por lo dicho anteriormente, qq es la proyecci´on ortogonal sobre Im q q ¡ a lo largo de pIm q qK q ¡K .
De la misma forma In qq es una proyecci´on de rango n 1: es la proyecci´on ortogonal sobre q ¡K a lo largo de q ¡.
126
Proyecciones y bases ortonormales
6.4.
El algoritmo cl´ asico de Gram-Schmidt
Recordemos que el m´etodo de ortonormalizaci´on de Gram-Schmidt consiste, a grandes rasgos, en lo siguiente: Dado un sistema de vectores a1 , . . . , an P Cm , que supondremos linealmente independientes, se trata de conseguir una base ortonormal tq1, . . . , qnu del subespacio a1, . . . , an ¡ de modo que
q1, . . . , qj ¡ a1, . . . , aj ¡,
j
1, . . . , n.
Para ello se procede de la siguiente forma: 1. q1
}aa1}
1 2
2. Para j 2, 3, . . ., qj es la proyecci´on ortogonal de aj sobre el subespacio a1, . . . , aj1 ¡K q1, . . . , qj1 ¡K dividido por su norma. Ahora bien, hemos visto en la Observaci´on 6.7 que si Qj 1
q1
q j 1
entonces la proyecci´on ortogonal sobre q1 , . . . , qj 1 ¡K es Im Qj 1 Qj 1 . Por lo tanto, vj pIm Qj 1 Qj 1 qaj
q1, . . . , qj1 ¡K, y as´ı pIm Qj1Qj1qaj vj qj }vj }2 }pIm Qj1Qj1qaj }2 .
es la proyecci´on ortogonal de aj sobre
Dado que esta construcci´on es inductiva puede implementarse algor´ıtmicamente. En primer lugar
vj pIm Qj 1 Qj 1 qaj
aj Qj1Qj1aj
y qj
aj
}vvj}
j 2
.
Qj 1
q1 aj q2 aj .. .
qj1 aj
j¸1
aj pqiaj qqi.
i 1
6.4 El algoritmo cl´asico de Gram-Schmidt
127
Por lo tanto, si ponemos v1 a1 rjj }vj }2 rij qi aj tenemos que para j
1, 2, . . . , n aj
vj
j°1
i 1
}vj }2qj
j °
rij qi
i 1
q1 q2
rij qi
j°1
i 1
rij qi
r1j r2j qj .. . rjj
Si constru´ımos una matriz con los vectores ai , A entonces A tiene rango completo y
A a1 a2
an
q1 q2
a1 a2
r11 r12 0 r22 qj .. .. . . 0 0
...
an
P
Cmn ,
r1n r2n .. . . rnn
Es decir, A QR con Q P Cmn una matriz con columnas ortonormales y R rrij s P Cnn una matriz triangular superior con rii ¡ 0. A esta forma de expresar A se le llama factorizaci´ on QR reducida de A. Y, tal y como hemos visto, el m´etodo de ortonormalizaci´on de Gram-Scmidt nos proporciona un algoritmo par obtener esta factorizaci´on:
128
Proyecciones y bases ortonormales
Algoritmo cl´ asico de Gram-Schmidt
Dada A a1 a2
an
P Fmn, m ¥ n, rangpAq n, (F R o C)
R=zeros(n,n) Q=A for j 1:n for i 1:j 1 rij qi aj qj qj rij qi end for rjj }qj }2 qj qj rjj end for La salida de este algoritmo debe ser una matriz triangular superior R con rii y una matriz Q con columnas ortonormales.
6.5.
¡0
El algoritmo modificado de Gram-Schmidt
Hay otro procedimiento te´orico que produce tambi´en una base ortonormal como la del m´etodo de Gram-Schmidt. Partimos de la siguiente observaci´on: Si Q P Cmn tiene columnas ortonormales la proyecci´on ortogonal QQ es suma de n proyecciones ortogonales de rango 1. En efecto
QQ
q 1 q2
q1
q 2 qn .. .
q1q1
q2 q2
qn qna st.
qn
Ahora recordemos que en el m´etodo cl´asico de Gram-Schmidt cada qj se obtiene proyectando aj sobre q1 , . . . , qj 1 ¡K , y despu´es normalizando. Es decir, qj
pI Q Q qa }pI m Q j1Qj1qa j} m
j 1
j 1
j 2
.
6.5 El algoritmo modificado de Gram-Schmidt
129
De acuerdo con lo que acabamos de ver j¸ 1 qi qi . Im Qj 1 Qj 1 Im i1
Pero si hacemos el producto
pIm qj1qj1qpIm qj2qj2q . . . pIm q1q1q y tenemos en cuenta que qi qj 0 si i j, resulta que este producto es Im qj 1 qj1 qj 2 qj2 . . . q1 q1 I Qj 1 Qj 1 . Im qiqi entonces Im Qj 1 Qj 1 Pj 1 Pj 2 . . . P1 .
En definitiva, si ponemos Pi
(6.2)
¿Cu´al es la interpretaci´on geom´etrica de esta igualdad?. En primer lugar, Im Qj 1 Qj 1 es la proyecci´on ortogonal sobre pIm Qj 1 qK q1 , . . . , qj 1 ¡K . En segundo lugar, Pi Im qi qi es la proyecci´on sobre qi ¡K . Por lo tanto, para un vector x P Cn cualquiera, pIm Qj 1 Qj 1 qx es la proyecci´on ortogonal de x sobre q1, . . . , qj1 ¡K y Pix es la proyecci´on ortogonal de x sobre qi ¡K. As´ı pues, Pj 1 Pj 2 . . . P1 x es el vector resultante de hacer lo siguiente: 1. Proyectamos x ortogonalmente sobre
q1 ¡K. Esto es P1x
2. El vector obtenido, P1 x, lo proyectamos ortogonalemnte sobre es P2 P1 x.
q2 ¡K. Esto
3. El vector obtenido, P2 P1 x lo proyectamos ortogonalmente sobre q3 es P3 P2 P1 x. 4. As´ı sucesivamente hasta proyectar Pj 2 . . . P2 P1 x sobre
¡K. Esto
qj1 ¡K.
La identidad (6.2) nos dice que es lo mismo hacer este proceso de proyecciones sucesivas que hacer directamente la proyecci´on de x sobre q1 , . . . , qj 1 ¡K . Esto segundo es lo que hace el m´etodo cl´asico de Gram-Schmidt. Pero el resultado te´orico es el mismo si se hacen las proyecciones sucesivas de rango n 1. La implementaci´on
130
Proyecciones y bases ortonormales
de este procedimiento es lo que denominaremos algoritmo modificado de GramSchmidt. La implementaci´on pr´actica del algoritmo difiere muy poco de la del algoritmo cl´asico. De hecho s´olo hay que cambiar una l´ınea:
Algoritmo modificado de Gram-Schmidt
Dada A a1 a2
an
P Fmn, m ¥ n, rangpAq n, (F R o C)
R=zeros(n,n) Q=A for j 1:n for i 1:j 1 rij qi qj (M´etodo Cl´asico: rij qj qj rij qi end for rjj }qj }2 qj qj rjj end for
qiaj )
Puede no ser evidente por qu´e la simple sustituci´on de rij qi aj por rij qi qj produce el mismo resultado en aritm´etica exacta, y puede que tampoco sea evidente por qu´e este algoritmo produce qj
pIn qj1qj1q pIn q1q1qaj .
Para comprenderlo analicemos con detalle el progreso de las operaciones rij qi qj seguida de qj qj rij qi para un j fijo y para i 1, . . . , j 1. La notaci´on que usaremos es la del ordenador; es decir, las variables no lo son en el sentido matem´atico sino en el del ordenador (lugares f´ısicos en la memoria del mismo) y “iguladad”significa “asignaci´on”. As´ı, para la expresi´on: qj
qj r1j q1
se debe entender que a la variable qj se le asigna el valor de restar al valor previo de qj el producto r1j q1 . Supongamos entonces fijado j
P t1, . . . , nu y procedamos con el bucle:
6.6 N´ umero de operaciones
131
for i 1:j 1 rij qi qj qj qj rij qi end for recordando que el valor de qj al comienzo del mismo es aj . Para i 1 r1j q1 qj q1 aj qj qj r1j q1 aj
Para i 2
pq1aj qq1 aj pq1q1qaj pIn q1q1qaj .
q2qj q2pIn q1q1qaj qj r2j q2 pIn q1q1qaj q2pIn q1q1qaj q2 pIn q1q1qaj q2q2pIn q1q1qaj pI q2q2qpIn q1q1qaj . Siguiendo as´ı sucesivamente obtendr´ıamos para i j 1 qj pI qj 1 qj1 q pI q2 q2 qpIn q1 q1 qaj . r2j qj
Por (6.2) resulta que qj
pIn Qj1Qj1qaj
de modo que el valor de qj , en aritm´etica exacta, obtenido por los procedimientos cl´asico y modificado de Gram-Schmidt es el mismo. Debe notarse que con ambos m´etodos, los elementos diagonales de R son n´ umeros reales positivos independientemente de si la matriz original es de n´ umeros reales o complejos. Observaciones 6.8 Cuando el tama˜ no de la matriz A es muy grande y puede haber problemas para alojar en memoria las matrices Q y R, en vez de asignar Q=A, se pueden ir sustituyendo las columnas de A por las de Q a medida que ´esta se va construyendo. Con los ordenadores actuales ´este es un problema menor, pero todav´ıa hay algoritmos que contemplan esta posibilidad.
6.6.
N´ umero de operaciones
Al igual que con el algoritmo LU , para los algoritmos que estudiemos en este curso iremos calculando su coste operacional; es decir, una estimaci´on del n´ umero
132
Proyecciones y bases ortonormales
de flops (“floating point operations”) que conlleva la implementaci´on pr´actica del mismo. Para ambos algoritmos de Gram–Schmidt tenemos el siguiente resultado: Teorema 6.9 Los algoritmos cl´asico y modificado de Gram–Schmidt requieren 2mn2 flops para calcular la factorizaci´on QR de una matriz m n.
Recordemos que el n´ umero de flops es un polinomio en el tama˜ no de la matriz y que de este polinomio s´olo se suele destacar el t´ermino de grado mayor dado que el resto son poco significativos salvo que n y m sean peque˜ nos. As´ı, el Teorema 7.5 nos dice que el coste de los algoritmos de Gram–Schmidt es del orden 2mn2 o Op2mn2 q; teniendo estos s´ımbolos el habitual significado asint´otico: n´ umero de flops m,nÑ8 2mn2 l´ım
1.
Demstraci´ on. El Teorema puede demostrarse de la siguiente forma: Cuando m y n son grandes, el mayor n´ umero de operaciones se encuentran en las dos siguiente sentencias: rij qi qj (M´etodo Cl´asico: rij qi aj ) qj qj rij qi La imera l´ınea calcula el product escalar qi qj y requiere m multiplicaciones y m 1 sumas. Y la segunda fila calcula qj qj rij qi precisando m multiplicaciones y n restas. As’ı pues, el trabajo que se hace en cada iteraci´on es del orden de 4m. Como hay que hacer i 1 : j 1 iteraciones para j 1 : n, el n´ umero total de flops necesarios es asint´otico a
1 n j¸ ¸
j 1i 1
6.7.
4m 4m
n ¸
pj 1q 4m
j 1
n¸1
j 1
j
4m pn 2 1qn 2mn2
Opmnq.
Existencia y unicidad de la factorizaci´ on QR
Hemos exigido que los vectores a1 , . . . , an sean linealmente independientes. En estas circunstancias la descomposici´on QR de A a1 an es u ´nica. En efecto, ambos m´etodos utilizan proyecciones ortogonales, y ´estas son u ´nicas. Por ejemplo, a1 el m´etodo cl´asico, empieza construyendo q1 }a1}2 , de modo que est´a definido
6.7 Existencia y unicidad de la factorizaci´on QR
133
de forma u ´nica; y para j 2, . . . , n , qj es la proyecci´on ortogonal de aj sobre el subespacio a1 , . . . , aj 1 ¡K y luego normalizado. La proyecci´on ortogonal de aj sobre un subespacio fijo es u ´nica; y, por supuesto, la norma del vector obtenido tambi´en. As´ı pues, la matriz Q obtenida es u ´nica. Por su parte, los elementos de R son las componentes de cada aj en la base
q1, . . . , qn ¡, que son u´nicos. Estos argumentos prueban Teorema 6.10 Si A factorizaci´on QR.
P
Cmn (m
¥
n) tiene rango completo admite una u ´nica
No obstante, debe notarse que hay varias formas de escribir A QR con Q una matriz con columnas ortonormales y R triangular superior. Por ejemplo, si A QR ˜ peiθ Qq tiene columnas ortonormales y tambi´en A peiθ Qqpeiθ Rq. cumple que Q iθ ˜ e R es triangular superior. Pero la condici´on rii P R exige que θ 0. R Analizamos ahora el caso en que rangpAq no es completo. En este caso, habr´ıa una columna, digamos aj , que ser´ıa combinaci´on lineal de a1 , . . . , aj 1 . Es decir, aj
P a1, . . . , aj1 ¡ q1, . . . , qj1 ¡ .
Como qj es la proyecci´on ortogonal de aj sobre el subespacio a1 , . . . , aj 1 ¡K una vez normalizado, y aj P q1 , . . . , qj 1 ¡ tenemos que qj 0. Mejor dicho, la proyecci´on ortogonal, vj , de aj sobre a1 , . . . , aj 1 ¡K es el vector cero y as´ı rjj }vj }2 0. Esto imposibilita construir el vector unitario y continuar el algoritmo en el punto en que se hace la divisi´on qj
rqj . jj
Este hecho parecer´ıa indicar que no existe factorizaci´on QR. Sin embargo, te´oricamente, podr´ıamos escoger para qj un vector unitario en q1 , . . . .qj 1 ¡K y proseguir con el proceso de Gram-Schmidt; y hacerlo tantas veces cuantas sea necesario. En la pr´actica no hay posibilidad de escoger de manera adecuada el vector qj ortonormal a q1 , . . . , qj 1 ¡. Esto s´olo es una posibilidad te´orica que nos permite asegurar que, incluso cuando A no es de rango completo, siempre admite una factorizaci´on QR.
134
Proyecciones y bases ortonormales
Teorema 6.11 Toda matriz A P Cmn (m ¥ n) admite una factorizaci´on QR. Debe observarse que si A no tiene rango completo entonces habr´ıa factorizaci´on QR de A, pero no se cumplir´ıa que Im Q Im A. En la pr´actica, dif´ıcilmente el valor calculado por el ordenador ser´a exactamente rjj 0, debido al redondeo. Muchos programas, como MATLAB, si el n´ umero por el que se va a dividir es muy peque˜ no, emiten un aviso de divisi´on por cero. Pero los errores de redondeo pueden ser de tal magnitud que MATLAB no detecte que tal divisi´on por cero se ha producido y nos devuelva un resultado. Por ejemplo >> A=ones(3); >> [Q,R]=clgs(A) Q = 5.7735e-01 -5.7735e-01 5.7735e-01 -5.7735e-01 5.7735e-01 -5.7735e-01 R = 1.7321e+00 1.7321e+00 0 3.8459e-16 0 0
5.7735e-01 5.7735e-01 5.7735e-01 1.7321e+00 3.8459e-16 8.5397e-32
Claramente rangpAq 1 pero MATLAB lleva a cabo el algoritmo sin ning´ un tipo de aviso. Vemos, enseguida, que las dos u ´ltimas filas de R son casi cero. Esto nos indica que algo raro sucede. Si hacemos >> Q’*Q ans = 1.0000e+00 -1.0000e+00 1.0000e+00 >> Q*R ans = 1 1 1 1 1 1
-1.0000e+00 1.0000e+00 -1.0000e+00
1 1 1
1.0000e+00 -1.0000e+00 1.0000e+00
6.8 Factorizaci´on QR reducida y completa
135
observamos que Q no es unitaria. El algoritmo nos devuelve una factorizaci´on de A en dos matrices Q y R. R es triangular superior con rii ¡ 0 pero Q no es unitaria. Por lo tanto, no se trata de una factorizaci´on QR de A. Este fen´omeno es bastante frecuente con los algoritmos de Gram-Schmidt y lo analizaremos con m´as detalle en la u ´ltima secci´on de este Tema.
6.8.
Factorizaci´ on QR reducida y completa
Al igual que con los valores singulares existe una factorizaci´on reducida y una factorizaci´on completa de cualquier matriz A P Cmn (m ¥ n). La reducida ya hemos visto lo que es: una descomposici´on de A en dos factores Q P Cmn y R P Cnn , con Q una matriz cuyas columnas son ortonormales y R una matriz triangular superior con elementos diagonales n´ umeros reales positivos (si A es de rango completo) o cero (si A no es de rango completo). Para obtener una factorizaci´on QR completa de A se a˜ naden a Q m n columnas mm ˜ ortonormales para formar una matriz Q P C que sea unitaria. Y se a˜ naden a R m n filas de ceros para obtener la matriz ˜ R
R
0pmnqn
˜ R, ˜ pero esta factorizaci´on no tiene por qu´e ser u Claramente A Q ´nica, debido ˜ A cualquier a la arbitrariedad de la elecci´on de las u ´ltima m n columnas de Q. ˜yR ˜ con estas propiedades (i.e. Q ˜ unitaria y descomposici´on de A en dos factores Q ˜ P Cmn con rij 0 para i ¡ j y rjj ¥ 0) se le llama factorizaci´ R on QR completa de A. Los algoritmos de Gram-Schmidt proporcionan factorizaciones reducidas. En la pr´oxima lecci´on estudiaremos otro algoritmo que da una factorizaci´on completa. Ya hemos dicho c´omo pasar de una factorizaci´on reducida a una completa. El rec´ıproco ˜R ˜ es una factorizaci´on completa de A P Cmn y m ¥ n, entonces es obvio: si A Q ˜ sus u se obtiene una factorizaci´on reducida eliminando de Q ´ltimas m n columnas ˜ y de R sus u ´ltimas m n filas de ceros.
136
6.9.
Proyecciones y bases ortonormales
Un an´ alisis experimetal de la estabilidad de los algoritmos de Gram-Schmidt
El experimento que se expone a continuaci´on, hecho con MATLAB, tiene como finalidad explorar la diferencia en la estabilidad num´erica entre los algoritmos cl´asico y modificado de Gram-Schmidt. Primero constru´ımos una matriz A de tama˜ no 80 80 con vectores singulares aleatorios y valores singulares en progresi´on geom´etrica desde 21 hasta 280 de raz´on 21 : disp(’Escogidas U y V matrices ortogonales 80x80 aleatorias’) S=diag(2.^(-1:-1:-80)); A=U*S*V’; Y usamos los algoritmos cl´asico y modificado de Gram-Schmidt para calcular factorizaciones QR de A: [QC,RC]=clgs(A); [QM,RM]=mgs(A); Finalmente, dibujamos los elementos diagonales rjj de R producidos por ambos algoritmos con una escala logar´ıtmica: axis([1 80 10^(-25) 10^0]) semilogy(1:80, diag(RC),’*’) hold on semilogy(1:80,diag(RM),’o’); semilogy(1:80, 2.^(-1:-1:-80),’r.-’); plot(1:80,sqrt(eps).*ones(80),1:80,eps.*ones(80)) Esto, m´as algunas sentencias para presentar la gr´afica con el texto apropiado, proporciona la Figura 6.1, que interpretamos a continuaci´on. Por una parte rjj }vj }2 con vj pIm Qj Qj qaj , de modo que rjj nos da una idea del tama˜ no del vector proyecci´on de aj sobre a1 , . . . , aj 1 ¡K . Salta a la vista
6.9 Un an´alisis experimetal de la estabilidad de los algoritmos de Gram-Schmidt137
0
10
−5
10
eps1/2 −10
rjj
10
−15
10
eps
−20
10
2
−j
−25
10
10
20
30
40
j
50
60
70
80
Figura 6.1: Elementos diagonales de R en la descomposici´on QR de A calculados con MATLAB mediante los algoritmos cl´asico y modificado de G-S. El algoritmo cl´asico produce los n´ umeros representados por y el modificado, los representados por . en la figura que rjj es una funci´on decreciente de j y que este decrecimiento se hace, sobre todo al principio, a lo largo de la recta 2j . Esto es bastante razonable porque A U SV
T
80 ¸
si ui viT
21u1v1T
22 u2 v2T
T s80 u80 v80 ,
i 1
y la j-´esima columna de A es aj
21vj1u1
22 vj2 u2
280 vj,80 u80 .
Ahora bien, las filas de V son ortonormales as´ı que
}vj }2
80 ¸
|vji|2 1,
i 1
y, como V es aleatoria, podemos suponer que todos los |vji | son, aproximadamente, de la misma magnitud. Es decir,
|vji|
1 80
1{2
801{2 0,1.
138
Proyecciones y bases ortonormales
As´ı pues aj
801{2
21 u1
22 u2
280 u80 .
Todos los ui son vectores unitarios por lo que la componente de aj en la direcci´on de ui es el doble que su componente en la direcci´on de ui 1 . En particular, la componente m´as fuerte de aj es en la direcci´on de u1 . Pero al hacer la descomposici´on QR de A, q1 tiene la direcci´on de a1 , que es, aproximadamente, la de u1 . Es decir, q1 u1 y, por lo tanto, r11 21 801{2 . En el segundo paso, v2
a2 pq1a2qq1 a2 pu1 a2qu1 801{2
22 u2
280 u80
801{222u2,
as´ı que q2 u2 y r22 801{2 22 . Y as´ı sucesivamente. Esto justifica por qu´e los rjj decrecen aproximadamente pegados a la recta 2j . Lo segundo que llama la atenci´on es que este comportamiento de rjj de ir pe? gado a la recta 2j se detiene en un cierto punto. Aproximadamente en eps para el algotitmo cl´asico de Gram-Schmidt y en eps para el modificado. Esto es una consecuencia del redondeo. Con el algoritmo cl´asico los rjj nunca son menores de 108 mientras que con el modificado se reducen del orden de 8 d´ıgitos decimales m´as, hasta, aproximadamente, 1016 ; el orden del ´epsilon de la m´aquina, eps. Nada mejor es esperable. Algunos algoritmos son m´as sensibles que otros a los errores de redondeo. Es conocido que al algoritmo cl´asico le pueden afectar mucho. Por ello rara vez se usa.
6.10.
P´ erdida num´ erica de ortogonalidad
Hay otro tipo de inestabilidad que afecta a ambos algoritmos del que ya hemos visto un ejemplo m´as arriba: ambos pueden calcular matrices Q que est´an lejos de tener columnas ortonormales. Una medida de la ortogonalidad de las columnas de una matriz nos la da el siguiente procedimiento: Si las columnas de Q P Cmn son ortonormales entonces Q Q In . Si las columnas de Q son “casi” ortonormales entonces Q Q es “casi” In ; es decir, los elementos de la matriz Q Q In son casi cero. O equivalentemente }QQ In}2 0.
6.10 P´erdida num´erica de ortogonalidad
139
La p´erdida de ortogonalidad ocurre cuando A est´a pr´oxima a no tener rango completo; es decir, cuando σn 0. Y como en la mayor parte de los fen´omenos de inestabilidad, ´esta aparece incluso en matrices de dimensi´on peque˜ na. Vamos a presentar dos ejemplos, uno para hacer a mano y otro con MATLAB. Consideremos la siguiente matriz A
0,70000 0,70711 0,70001 0,70711
en un hipot´etico ordenador que redondea los resultados a 5 d´ıgitos de exactitud relativa. Los algoritmos cl´asico y modificado de Gram-Schmidt son id´enticos para matrices de dos columnas. Observemos, en primer lugar, que σ2 pAq 5,0252 106 es la distancia de A a las matrices de rango 1. Calculemos la descomposici´on QR que se obtendr´ıa en esta situaci´on (recordemosque lo hacemos redondeando a 5 d´ıgitos significativos). En el primer paso del algoritmo a1
0,70000 0,70001
ñ r11 }a1}2
ñ q1 a1{r11
a
0,72
0,700012
0,98996 ñ
0,70710 0,70711
En el segundo paso
0,70711 r12 q1 a2 0,70710 0,70711
0,70711
0,70711 v2 a2 pq a2 qq1 1
1,0000p045 . . .q
0,70710 1 0,70711 0,70711
0,00001 0,00000
Los errores de redondeo en este c´alculo de v2 son determinantes. De hecho q2 y la matriz Q obtenida es Q
1 0
0,70710 1 , 0,70711 0
que es una matriz lejos de ser ortogonal. Para serlo la primera columna deber´ıa ser 0 , muy lejos de q1 . De hecho 1
}QQ In}2 0,70710,
140
Proyecciones y bases ortonormales
un n´ umero muy grande. En un ordenador con 16 d´ıgitos de precisi´on como los habituales y usando MATLAB con esta misma matriz: >> A=[ 0.70000 0.70711; 0.70001 0.707111]; >> [Q,R]=mgs(A); >> norm(Q’*Q-eye(2)) produce ans= 2.3014e-11 un n´ umero peque˜ no pero muy grande comparado con el ´epsilon de la m´aquina, eps, que es del orden de 2.2204e-16. Si hacemos lo mismo con el algoritmo QR implementado por MATLAB y que estudiaremos en la pr´oxima lecci´on obtenemos >> [Q,R]=qr(A); >> norm(Q’*Q-eye(2)) ans = 2.3514e-16 que es del orden del ´epsilon de la m´aquina.
