Cap´ıtulo 7
Reflexiones de Householder y el Algoritmo QR
7.1.
Introducci´ on
Hay otros dos m´etodos importantes para calcular la descomposici´on QR de una matriz: las reflexiones de Householder y las rotaciones de Givens. Ambos comparten la propiedad de ser procesos que conducen a una “triangularizaci´on de la matriz dada mediante transformaciones ortogonales” (caso real) “o unitarias” (caso complejo). En este curso estudiaremos el primero de los dos m´etodos (el de las reflexiones de Householder). Una buena referencia para el m´etodo de las rotaciones de Givens es el libro de Golub y van Loan [9]. En realidad, estos dos m´etodos, son el resultado de aplicar un teorema general de a´lgebra que establece que toda transformaci´on ortogonal (o unitaria) es un producto de simetr´ıas (reflexiones o giros). Los m´etodos de Householder y Givens difieren de los de Gram-Schmidt en un hecho b´asico fundamental. Para ver esta diferencia analicemos, de nuevo, el proecedi141
142
Reflexiones de Householder y el Algoritmo QR
miento de ortonormalizaci´on de Gram-Schmidt: dada A a1 a2 1 a1 con r11 }a1 }2 . As´ı, si (F=R o C) definimos q1 r11 R1
Diag
tenemos que
AR1 A continuaci´on, ponemos r22 }a2 }2 y A2
q2
Si definimos
1 0 R2 0 . .. 0
tenemos que A3
an
P Fm n
1 , 1, . . . , 1 r11
1 1 p a2 pq1 a2 qq1 q r22 r22
q 1 a2
a2
r12 a1 r11
0 0 1 .. . . . . 0
rr
12 22
1 r22
0 .. .
0
A2R2 A1R1R2
an .
r1 a2 r r12r 22
a1
11 22
0 0 0 .. . 1
q 1 q 2 a3
an .
Procediendo de esta forma sucesivamente podemos construir, para la matriz dada A, matrices R1 , R2 ,. . . , Rn de modo que AR1 . . . Rn
Q,
ˆ siendo Q una matriz cuyas columnas son ortonormales. Se pone, entonces, R ˆ la descomposici´on QR, reducida, de A. Rn1 . . . R11 y se obtiene A QR,
El m´etodo de Householder; y tambi´en el de Givens, cambia completamente el punto de vista: se trata de realizar sobre la matriz A una serie de transformaciones unitarias (u ortogonales, en el caso real) Q1 , Q2 . . . , Qn tales que Qn . . . Q1 A R sea una matriz triangular superior. En lo sucesivo hablaremos de transformaciones y matrices unitarias pero se debe entender que si A es real las transformaciones y matrices son ortogonales.
7.1 Introducci´on
143
Lo primero que hay que observar es que mientras los m´etodos de ortonormalizaci´on de Gram-Schmidt produce una factorizaci´on QR reducida, el de Householder produce una factorizaci´on completa. En aquel caso Q es de tama˜ no m n y R es n n. En el m´etodo de Householder cada Qi es unitaria (o, insistimos por u ´ltima vez, ortogonal en el caso real);por lo tanto, de tama˜ no m m. Y R ser´a de tama˜ no m n. Sucintamente, el m´etodo de Householder consiste en construir matrices unitarias elementales (reflexiones) que operan sucesivamente sobre las columnas de A para reducir esta matriz a una forma triangular superior. Esquem´aticamente el proceso ser´ıa el siguiente:
x x x
x x x x x x x x x
x x x
Q Ñ 1
x 0 0 0 0
x x x x x x x x x x Q1 A
Q Ñ 2
x x x x x 0 x 0 x 0 x Q2 Q1 A
Q Ñ 3
x x x x x x 0 0 Q3 Q2 Q1 A
Q1 opera sobre las filas 1 a 5 para hacer ceros en todas las posiciones de la primera columna excepto en la posici´on de la diagonal. A continuaci´on se construye una matriz Q2 , unitaria u ortogonal en el caso real, que opera sobre las filas 2 a 5 de la matriz obtenida en el paso anterior, para hacer ceros en todas las posiciones de la segunda columna por debajo de la diagonal principal. Y as´ı sucesivmente. Es decir, para k 1, 2, . . . , n, Qk opera sobre Ak p1 : k, 1 : nq haciendo ceros en las posiciones pk 1, kq, pk 2, kq,. . . , pm, kq, donde A1 A y para k ¡ 1, Ak Qk1Ak1. Se debe observar que Qk no act´ ua sobre las primeras k 1 filas de Ak1 ; es decir, las deja intactas, por lo que tendr´a la siguiente forma: Qk
Ik1 0 ˜k 0 Q
˜ k una matriz unitaria cuyo objetivo es reducir el vector formado por los siendo Q elementos en las posiciones pk, k q, pk 1, k q,. . . , pm, k q de Ak1 a un vector cuya u ´nica componente posiblemente no nula sea la primera:
x
0
x x x 0 ˜k Q .. .. . . .
144
Reflexiones de Householder y el Algoritmo QR
La primera cuesti´on surge de forma natural: ¿existe siempre, para cualquier vector, una matriz unitaria que lo reduce a otro cuya u ´nica componente no nula sea, posiblemente, la primera?. Como veremos enseguida la respuesta es afirmativa y es la base del m´etodo de Householder.
7.2.
Las reflexiones de Householder
El objetivo es resolver el siguiente problema: dado un vector x P Cn1 construir una matriz unitaria Q P Cnn tal que Qx αe1
siendo e1 el vector p1, 0, . . . , 0q.
Antes de nada debemos observar que, como Q es unitaria,
|α|2 }Qx}22 xQQx xx }x}22, de donde deducimos que }x}2 |α|.
La idea de Householder fu´e construir tal matriz como la matriz de una reflexi´on: una simetr´ıa respecto de un hiperplano. El caso particular de R2 nos servir´a como gu´ıa de lo que queremos conseguir: dado un vector x P R2 debemos encontrar una hiperplano (en este caso, una recta que pase por el origen) tal que el sim´etrico de x respecto de este hiperplano sea αe1 (Figura 7.1).
H
v1 x v
||x||2 e1
Suponiendo que α ¡ 0, dicha recta (H en la figura) es la bisectriz del ´angulo que forma los vector x y e1 . Sea v αe1 x }x}2e1 x (este vector es el de la Figura 7.1 pero trasladado al origen). Como el tri´angulo formado por los vectores x, }x}2e1 y v (trasladado desde el origen al extremo de x) es is´osceles, resulta que H es perpendicular a v. Es decir, H
Figura 7.1: Reflexi´on sobre el eje de abscisas
v ¡K
¿Cu´al es la transformaci´on que lleva x sobre }x}2 e1 ?
7.2 Las reflexiones de Householder
145
Enseguida notamos que v1 es la proyecci´on ortogonal de x sobre v a lo largo de v H. Por lo tanto, si q en la direcci´on de v, tenemos que }v}2 es el vector unitario la proyecci´on ortogonal sobre v es Pq qq ; y
Pq x qqx. Como, finalmente, v 2v1 y x v }x}2 e1 conclu´ımos que pI 2qqqx }x}2e1. As´ı pues, la matriz Q que representa la reflexi´on respecto de v ¡K y que hace que Qx }x}2 e1 es Q I 2qq v donde q }v} y v }x}2e1 x. v1
2
Hay otra forma de reflejar x sobre el eje de abscisas: hacerlo sobre la direcci´on negativa de ´este. Es decir, la reflexi´on ser´ıa Q1 x }x}2 e1 .
H−
Razonando como antes, pero respecto de la Figura 7.2, ser´ıa Q1 I 2q1 q1 con q1 }uu}2 y u }x}2 e1 x. Todo lo anterior nos permite asegurar que en el caso R2 , para un vector x dado, siempre podemos encontrar una transformaci´on ortogonal Q de forma que Qx }x}2e1. Adem´as, nos dice que la matriz Q es una reflexi´on. Es decir, Q I 2q1 q1 con Q1 }uu}2 y u }x}2 e1 x. esta es nuestra gu´ıa para demostrar que lo mismo es posible en Cn .
H+
u1 u −||x||2 e1
x
v +||x||2 e1
Figura 7.2: Las dos reflexi´on sobre el eje de abscisas
x1 .. Lema 7.1 Sea x . P Cn y escribamos x1 |x1 |eiθ . Existen vectores unitarios xn n u1 , u2 P C tales que si Qui pIn 2ui ui q, i 1, 2, entonces Qu1 }x}2 eiθ e1 y Qu2 }x}2 eiθ e1 .
146
Reflexiones de Householder y el Algoritmo QR
Observaci´ on: En el caso complejo no se puede asegurar, en general, que la primera componente de Qui x vaya a ser real. Demostraci´ on.- Elegimos los vectores u1 y u2 siguiendo la orientaci´on marcada por el caso R2 . Definimos v1
}x}2eiθ e1 x,
y u1
v2
}vv1} ,
}x}2eiθ e1 x,
u2
1
}vv2} . 2
Veremos que, en efecto, Qu1 x }x}2 eiθ e1 . La demostraci´on de que Qu2 se hace igual. Qu1 x pIn 2u1 u1 qx x 2pu1 xqu1
}x}2eiθ e1
x 2 pvvv1 xq1{2 pvvv1q1{2 x 2 vv1vx v1. 1 1
1 1
1 1
Ahora bien
p}x}2eiθ eT1 xqx }x}2eiθ x1 xx }x}2eiθ eiθ |x1| }x}22 }x}2p|x1| }x}2q.
v1 x Y
v1 v1
p}x}2eiθ eT1 xqp}x}2eiθ e1 xq }x}22eT1 e1 }x}2eiθ eT1 x }x}2eiθ xe1 xx }x}22 }x}2eiθ x1 }x}2eiθ x¯1 }x}22 2}x}22 }x}2 eiθ |x1 |eiθ }x}2 eiθ |x1 |eiθ 2}x}22 2|x1 |}x}2 2}x}2 p}x}2 |x1 |q.
Por lo tanto Qu1 x x 2
}x}2p|x1| }x}2q v x 1 2}x}2 p}x}2 |x1 |q
v1
x }x}2eiθ e1 x }x}2eiθ e1,
tal y como se quer´ıa mostrar. Antes de continuar veamos un ejemplo en R4 . Ejemplo 7.2 Sea x p3, 1, 5, 1q. Encontrar una reflexi´on Qu tal que Qu x pα, 0, 0, 0q.
7.2 Las reflexiones de Householder
147
Las operaciones las haremos con ayuda de MATLAB. Seg´ un el Lema 7.1 debemos poner v Ahora
}x}2eiθ e1 x
}x}2
?
9
1
y u
v }v}2 .
1 6 y x1
25
3.
Por lo tanto "
v
5, 1q v1 p6, 0, 0, 0q p3, 1, 5, 1q pp3,9,1, 1, 5, 1q v2
Entonces u1
16 p3, 1, 5, 1q
As´ı
y v2
? 1 p9, 1, 5, 1q. 116
Qu1
1{2 1{6 1 { 6 17 pI4 2u1u1 q 5{6 5{{18 18 1{6 1{18
Entonces
5 {6 1{6 5{18 1{18 . 7{18 5{18 5{18 17{18
6
Qu1 x
0 . 0
0 De la misma forma
Qu2
1{2 1{6 5{6 1{6 5{54 1{54 . pI4 2u2u2 q 51{{66 535{{54 54 29{54 5{54 1{6 1{54 5{54 53{54
Y Qu2 x 6e1 . Ahora podemos describir el proceso de triangularizaci´on de Householder. La inducci´on del siguiente teorema lo pone claramente de manifiesto.
148
Reflexiones de Householder y el Algoritmo QR
Teorema 7.3 Sea A P Cmn y m ¥ n. Existen vectores unitarios u1 , . . . , un tales que para i 1, . . . , n, ui P Cmi 1 y si Qui Imi 1 2ui ui y Qi entonces
Ii1 0 0 Qui
P Cmm
Qn . . . Q1 A DR
siendo D P C una matriz diagonal unitaria y R P Cmn una matriz triangular m m
superior cuyos elementos diagonales son n´ umeros reales.
Demostraci´ on.- Procederemos por inducci´on sobre n. Si n 1 entonces A P Cm1 . Si A 0 escogemos cualquier vector unitario u P Cm1 y ponemos D Im y R 0 P Cm1 . Claramente, Qu A DR. Si A 0, por el Lema 7.1 existe u1 P Cm1 , unitario, tal que Qu1 A }A}2 eiθ e1 (con se quiere decir que se puede elegir u1 para que Qu1 A }A}2 eiθ e1 y se puede elegir u1 para que Qu1 A }A}2 eiθ e1 . Arrastraremos esta notaci´on), siendo θ el argumento de la primera componente de A. Entonces Qu1 A DR
con D
Diagpeiθ , Im1q
}A}2
0 y R .. . 0
.
Supongamos el teorema verdadero para matrices hasta con n 1 columnas y sea m n A P C . Escribamos A a1 A1 . Por el primer paso de la inducci´on, existe u1 P Cm1 , unitario, tal que Qu1 a1 }a1 }2 eiθ e1 . As´ı
}a1}2eiθ Qu A 0 1
xT B
con x P Cpn1q1 y B P Cpm1qpn1q . Aplicamos la hip´otesis de inducci´on a B. Existen vectores unitarios u2 , . . . , un tales que para i 2, . . . , n, ui P Cpmi 1q1 y si Qui Imi 1 2ui ui y I 0 i 2 ˜i Q 0 Qui
7.2 Las reflexiones de Householder
˜2 (Q
Qu ) entonces
149
˜n . . . Q ˜ 2C Q
2
˜ D˜ R,
˜ P Cpm1qpm1q una matriz diagonal unitaria y R ˜ P Cpm1qpn1q una matriz siendo D triangular superior cuyos elementos diagonales son n´ umeros reales. Sea Q1 Qu1 y Qi Entonces
Poniendo
1 0 ˜i 0 Q
P Cmm,
i 2, . . . , n.
}a1}2eiθ Qn . . . Q2 Q1 A 0 ˜ q, D Diagpe , D iθ
y
xT ˜R ˜ . D
}a1}2 R 0
xT ˜ R
obtenemos el resultado deseado. Observaciones 7.4 1. Las matrices de la forma Qu In 2uu , con u P Cn1 unitario, se llaman matrices de Householder o transformaciones de Householder o reflexiones de Householder, por ser Householder quien puso de manifiesto su inter´es en el c´alculo num´erico. Estas matrices tienen algunas propiedades importantes: (i) Son matrices herm´ıticas (sim´etricas en el caso real). En efecto, Qu
pIn 2uuq In 2uu In 2uu Qu
(ii) Son unitarias: Qu Qu
(iii)
pIn 2uuqpIn 2uuq In 4uu 4uuuu In donde hemos utilizado que u u 1 porque u es unitario. 1 De las dos propiedades anteriores se deduce que Q u Qu Qu . Es decir, la inversa de una matriz de Householder es ella misma.
2. Para que una matriz D P Cmm sea diagonal y unitaria debe ser de la forma D Diagpeiθ1 , . . . , eiθm q. Por lo tanto, si A P Rmn es real en el Teorema 7.3, debemos concluir que D Im y que las matrices Qui y R son reales. Es decir, en el caso real las matrices de Householder son ortogonales.
150
Reflexiones de Householder y el Algoritmo QR
3. En el Teorema 7.3 no se ha fijado el signo de los elementos diagonales de la matriz R. Ello es debido a que dos elecciones de vectores ui es posible. Es decir, es posible elegir tanto vi }ai }2 eiθ e1 ai como vi }ai }2 eiθ e1 ai . Para ambas elecciones, definiendo ui vi {}vi }2 , las matrices Qui I 2ui ui permiten llevar la matriz A a forma triangular. El signo de los elementos diagonales depende de la elecci´on realizada. Ahora bien, esta elecci´on no es inocua desde el punto de vista del c´alculo num´erico. El siguiente ejemplo, hecho con MATLAB, lo pone de manifiesto: >> A=[2 3; 10^(-7) 4] A = 2.000000000000000 0.000000100000000
3.000000000000000 4.000000000000000
>> v1=norm(A(:,1))*[1;0]-A(:,1), v2=-norm(A(:,1))*[1;0]-A(:,1) v1 = 1.0e-07 * 0.000000022204460 -1.000000000000000 v2 = -4.000000000000002 -0.000000100000000 >> u1=v1/norm(v1), u2=v2/norm(v2) u1 = 0.000000022204460 -1.000000000000000 u2 = -1.000000000000000 -0.000000025000000 >> Q1=eye(2)-2*u1*u1’, Q2=eye(2)-2*u2*u2’ Q1 = 0.999999999999999 0.000000044408921 0.000000044408921 -0.999999999999999 Q2 = -0.999999999999999 -0.000000050000000 -0.000000050000000 0.999999999999999
7.2 Las reflexiones de Householder
>> Q1*Q1’, Q2*Q2’ ans = 1.000000000000000 -0.000000000000000 ans = 1.000000000000001 0.000000000000000 >> R1=Q1*A, R2=Q2*A R1 = 2.000000000000002 -0.000000011182158 R2 = -2.000000000000003 -0.000000000000000
151
-0.000000000000000 1.000000000000000 0.000000000000000 1.000000000000000
3.000000177635681 -3.999999866773233
-3.000000199999997 3.999999849999995
La primera columna, a1 , de A est´a “casi” en la direcci´on de e1 p1, 0q. Si elegimos v1 }a1 }2 e1 a1 la matriz Q1 es ortogonal pero el elemento de la posici´on p2, 1q de R1 Q1 A es del orden de 108 , un orden m´as que el elemento p2, 1q de A. Sin embargo, si elegimos v2 }a1}2e1 a1 entonces Q2 es un “poco menos” ortogonal que Q1 , pero R2 Q2 A es perfectamente triangular (hasta la precisi´on con la que trabajamos 1016 ). Este ejemplo no es especialmente rebuscado, sucede muy a menudo. Si los vector x P Rm1 y }x}2 e1 son “casi” iguales (i.e., x1 }x}2 , x2 , x3 ,. . . , xm son “casi” cero), al hacer la diferencia
v
}x}2 x1 x 2
}x}e1 x
.. . xm
tenemos que hacer la resta }x}2 x1 de dos cantidades casi iguales. Este problema, tal y como se ve en los problemas del Cap´ıtulo 4, est´a mal condicionado y los errores de redondeo pueden dar lugar a c´alculos inexactos. Esto no sucede
152
Reflexiones de Householder y el Algoritmo QR
si elegimos
v
}x}2 x1 x2
}x}e1 x
.. . xm
porque entonces lo que tenemos que hacer es sumar |x}2 problema bien condicionado De la misma forma si x y }x}2e1 x.
x1 , y la suma es un
}x}2e1 son “casi” iguales y debemos elegir v
En conclusi´on, la elecci´on apropiada para v es
signpx1q}x}2eiθ e1 x donde signpx1 q es el signo de x1 ; es decir, signpx1 q 1 si x1 ¡ 0 y signpx1 q 1 si x1 0. Y para ser completamente exahustivos en este punto debemos considerar el caso en que x1 0. En este caso se define signpx1 q 1. Se v
debe observar que MATLAB no considera el caso en que x1 pueda ser cero y hay que definirlo expl´ıcitamente. Debe observarse tambi´en que puesto que Qu I 2uu se puede, igualmente, elegir v
signpx1q}x}2eiθ e1
x,
que es la elecci´on habitual. 1 4. Teniendo en cuenta que Q un el Teorema 7.3 para cada u Qu y que, seg´ mn APC existen matrices de Householder Qui , i 1, . . . , n tales que
Qun . . . Qu1 A DR, tenemos que
A Qu1 . . . Qun DR.
Si ponemos Q Qu1 . . . Qun D, resulta que A QR
es una factorizaci´on QR de A. Esta es una factorizaci´on QR completa; para obtener una factorizaci´on reducida basta suprimir de Q las u ´ltimas m n columnas y de R sus u ´ltima m n filas (que son filas cero).
7.2 Las reflexiones de Householder
153
5. Vamos a hacer expl´ıcito el proceso constructivo de factorizaci´on QR que est´a detr´as de la inducci´on en el Teorema 7.3. El objetivo es formular un algoritmo que produzca una factorizaci´on QR de una matriz arbitraria A P Rmn . Por ello, nos centraremos en el caso real. Supongamos A P Rmn dada. El Teorema 7.3 nos dice que existen vectores unitarios u1 , u2 ,. . . , un , ui P Rmi 1 tales que si Qui Imi 1 2ui ui y Qi DiagpIi1 , Qui q entonces
r11
Qi1 . . . Q1 A 0 0
Como Qi
... 0 0
ri1 i1
Ai
DiagpIi1Qu q i
r11
Qi Qi1 . . . Q1 A 0 0
..
0 0
.
ri1i1
Qui Ai
Es decir la acci´on de Qi sobre la matriz obtenida hasta el paso i 1 s´olo afecta a la submatriz Ai . Y Qui Ai
pImi 1 2uiui qAi Ai 2uiui Ai.
Esta es la operaci´on que tenemos que realizar en cada paso de el proceso iterativo sobre A. Algo similar sucede con la matriz Q. Recordemos que Q r i Q1 Q2 qi , i 1, . . . n tenemos que ponemos Q r Q r i1 Qi
Q1Q2 Qn. Si
Ii1 0 . 0 Qui
Por lo tanto, en el paso i s´olo quedan afectadas las m i u ´ltimas columnas de r As´ı pues Qi. r i p:, i : mqQu Q i
Qri1p:, i : mq 2Qri1p:, i : mquiui .
154
Reflexiones de Householder y el Algoritmo QR
Algoritmo QR de Householder (Caso real) Dada A P Rmn R A, Q Im for k 1 to n x Rrpk, . . . , mq|pk qs u signpx1 q}x}2 e1 x u u }u}2 Rrpk, . . . , mq|pk, . . . , nqs
Rrpk, . . . , mq|pk, . . . , nqs 2uu Rrpk, . . . , mq|pk, . . . , nqs Qp:, k : mq Qp:, k : mq 2Qp:, k : mquu end for
La salida de este algoritmo debe ser: una matriz ortogonal Q y una matriz triangular superior R tales que A QR. Se trata, como ya se ha dicho, de una factorizaci´on QR completa de A. Ya se ha indicado m´as arriba como conseguir una factorizaci´on reducida a partir de ´esta. Una observaci´on final: El algoritmo QR de Householder, en el caso real, devuelve una matriz triangular superior R cuyos elementos diagonales son n´ umeros reales pero que pueden ser positivos o negativos. Los algoritmos de Gram–Schmidt devuelven matrices R con n´ umeros reales positivos en la diagonal. Se puede conseguir que el algoritmo de Householder tambi´en proporcione matrices R con n´ umeros reales positivos (o negativos) en la diagonal. Para ello, si el elemento en la posici´on pk, k q es, digamos, negativo, multiplicamos la k-´esima fila de R por 1 y tambi´en la k´esima columna de Q por 1. Es decir, si Ek Diagp1, . . . , 1, . . . , 1q (1 en la k-´esima posici´on) entonces A QEk Ek R. Como Ek es unitaria, tambi´en lo es QEk . Y tambi´en Ek R es triangular superior. En el algoritmo bastar´ıa a˜ nadir las siguientes l´ıneas de c´odigo justo antes de finalizar el bucle for:
if Rpk, k q 0 Rpk, k : nq Rpk, k : nq Qp:, k q Qp:, k q end if
7.3 N´ umero de operaciones
7.3.
155
N´ umero de operaciones
Estudiamos a continuaci´on el coste del algoritmo QR por transformaciones de Householder. Antes de nada debemos hacer una observaci´on que es pertinente tambi´en para la siguiente secci´on donde estudiaremos la estabilidad de este algoritmo: tal y como se ha expuesto aqu´ı el algoritmo (recuadro con el Algoritmo de Householder (Caso real)) la matriz Q se forma expl´ıcitamente. Sin embargo, esto no es estrictamente necesario ya que Q es producto de reflexiones de Householder y ´estas quedan determinadas completamente por los sucesivos vectores unitarios u. Podr´ıamos pensar que lo que devuelve el algoritmo no son las matrices R y Q sino la matriz R y los n vectores unitarios u que permiten (si se desea) construir la matriz Q. Teniendo en cuenta esto, las operaciones del algoritmo se encuentran en el bucle (1) (2) (3) (4)
x Rrpk, . . . , mq|pk qs u signpx1 q}x}2 e1 x u u }u}2 Rpk : m, k : nq Rpk : m, k : nq 2uu Rpk : m, k : nq
que se realiza para k 1 : n. El n´ umero m k 1 aparece muy frecuentemente porque es el tama˜ no del vector x. Pongamos ` m k 1 para abreviar la notaci´on. La sentencia (1) es una asignaci´on y por lo tanto en ella no se realizan flops. En la sentencia (2) se hacen: ` multiplicaciones: x21 , x22 ,. . . , x2` . ` 1 sumas: x21 1 ra´ız cuadrada:
a
x21
x2` .
x2`
1 multiplicaci´o en }x}2 e1 . 1 suma para la primera componente de x. En total 2`
2
2pm k
2q flops. Como k
1 : n, el n´umero total de flops
156
Reflexiones de Householder y el Algoritmo QR
necesarios para implementar la sentencia (2) es n °
k 1
2pm k
2q
pm n 1 kq 2npm n 1q 2 npn2 1q k 1 2mn n2 potencias menores en m y n. 2
n °
Para implementar la sentencia (3) se necesita un n´ umero similar de flops. Es en la sentencia (4) donde se acumula la mayor parte de las flops como veremos ahora. Analicemos cu´antas son necesarias para cada j k : n: Rpk : m, j q Rpk : m, j q 2upu Rpk : m, j qq Los par´entesis en upu Rpk : m, j qq indican la forma que asociamos para realizar los productos. N´otese que u Rpk : m, j q es un producto escalar y, por tanto, un n´ umero cuya formaci´on requiere ` productos y ` 1 sumas. En total, 2` 1 flops. Ahora, si ponemos x u Rpk : m, j q tenemos que multiplicar el n´ umero x por 2u; es decir, 2x requiere un producto y 2xu requiere ` productos para un total de ` 1 flops. Finalmente, si v 2xu tenemos que hacer la resta Rpk : m, j q v; lo cual requiere ` flops. As´ı pues, el n´ umero total de flops para cada j k : n es p2` 1q p` 1q ` 4`. Y el n´umero total de flops involucrados en la sentencia (4) para cada k 1 : n son 4`pn k 1q 4pm k 1qpn k 1q. . As´ı pues, el n´ umero total de flops necesarios para la sentencia (4) es n °
k 1
4pm k
1qpn k
1q
pm n kqk 4pm nq ° k 4 ° k2 k 1 k 1 k1 npn 1qp2n 1q npn 1q 4pm nq 2 4 6 2 3 2 2mn 3 n potencias menores en m y n.
4
n °
n
Teorema 7.5 El algoritmo QR de Householder (caso real) requiere flops para calcular la factorizaci´on QR de una matriz m n.
n
2mn2 32 n3
2 N´otese que este algoritmo requiere del orden de n3 flops menos que los algoritmos 3 de Gram-Schmidt.
7.4 Estabilidad del algoritmo de Householder
7.4.
157
Estabilidad del algoritmo de Householder
Estrictamente hablando el algoritmo de Householder para la factorizaci´on QR no puede ser estable hacia atr´as si se entiende que el objetivo del mismo es calcular para cada matriz A matrices Q, unitaria (ortogonal), y R, triangular superior con n´ umeros reales en la diagonal. La raz´on es la misma que la que impide que ning´ un algoritmo para el c´alculo de la descomposici´on en valores singulares de una matriz sea estable hacia atr´as (Problema 11 de la Lecci´on 5). Dicho lo anterior, y recordando (ver Secci´on 7.3) que la formaci´on de la matriz Q en el algoritmo de Householder se hace a partir de los vectores unitarios que se construyen en los pasos u signpx1 q}x}2 e1 u u }u}2
x
del algoritmo, podemos interpretar que el resultado que se busca no es la matriz Q expl´ıcitamente, sino cada uno de los vectores unitarios v1 , v2 ,. . . , vn que, una vez normalizados, se usan para construir las reflexiones cuyo producto permite obtener la matriz unitaria (ortogonal) Q: vi Q Q1 Q2 Qn , Qi Im 2ui ui , ui }vi} .
Si se entiende que el problema de la factorizaci´on QR consiste en, para cada A P Fmn (n ¤ m), obtener vectores v1 , v2 ,. . . , vn en Fm1 tales que si Qi Im 2ui ui , vi ui }vi} , y Q Q1Q2 Qn entonces A QR con R triangular superior con n´umeros reales en la diagonal, entonces el algoritmo de Householder para la factorizaci´on QR es estable hacia atr´as para cualquier matriz A. En efecto, tal algoritmo aplicado r con n´ a A P Fmn producir´ıa una matriz triangular superior R, umeros reales en la vri r ri diagonal, y vectores vri , i 1, . . . , tales que definiendo u }vri} , Qi Im 2uriuri , rQ r1 Q r n entonces Q rR r es una factorizaci´ i 1, . . . , n, y Q on QR exacta de una r matriz A pr´oxima a A. El siguiente teorema, que damos sin demostraci´on, establece la estabilidad hacia atr´as del algoritmo de Householder rigurosamente: Teorema 7.6 Supongamos que el algoritmo QR por reflexiones de Householder se aplica en un ordenador que cumple los dos axiomas de la aritm´etica de punto flotante
158
Reflexiones de Householder y el Algoritmo QR
ryR r como (5.2) y (5.4) de la Lecci´on 5. Sea A P Fmn una matriz arbitraria. Sean Q r A δA tal que m´as arriba. Entonces existe una matriz A
}δA} Op q M }A} rR rA yQ
δA.
