Problema de valores propios
Descrição do Produto
Cap´ıtulo 10
El problema de valores propios
10.1.
Introducci´ on
El objetivo de este cap´ıtulo es estudiar los algoritmos m´as importantes para calcular los valores propios de matrices generales; es decir, matrices que no son sim´etricas o herm´ıticas. Los algoritmos que veremos aqu´ı son aplicables tambi´en a matrices sim´etricas y herm´ıticas pero para ´estas hay, adem´as, otros algoritmos espec´ıficos. ´ En Algebra Lineal Num´erica se suele hablar de algoritmos directos y de algoritmos iterativos para distinguir entre los algoritmos que producen el resultado deseado tras un n´ umero finito de iteraciones (algoritmos directos como los de Gauss o QR) y aquellos que tienen como objetivo la construcci´on de una sucesi´on de objetos que converge al resultado deseado. Como veremos enseguida hay una raz´on muy poderosa para que los algoritmos para el c´alculo de valores y vectores propios sean del segundo tipo. Y siendo las cosas as´ı se habla de m´etodos directos e iterativos para el c´alculo de valores y vectores propios con un significado diferente a los anteriores. 215
216
El problema de valores propios
Los primeros buscan obtener todos los valores propios y (opcionalmente) los vectores propios; es decir, la forma de Schur de la matriz. Los m´etodos iterativos tratan de calcular un valor propio seleccionado y un vector propio, y de forma reiterada llegar a encontrar, si es posible, todos ellos. El algoritmo QR es el paradigma de los m´etodos directos y el M´etodo de las Potencias el de los m´etodos iterativos. En este Cap´ıtulo trataremos los m´etodos iterativos.
10.2.
El m´ etodo de las potencias
El objetivo del m´etodo de las potencias es calcular un vector propio de una matriz dada A P Fnn . Es extremadamente simple: Dado un vector no nulo x P Fn se construye la sucesi´on tAk xu, k 1, 2, . . .. Esta sucesi´on de vectores puede no converger a un vector, pero frecuentemente la sucesi´on de los subespacios generados por sus elementos,t Ak x ¡u, s´ı. Recordemos que para cada valor propio no hay un vector propio determinado un´ıvocamente, sino todo un subespacio de ellos (exceptuando el vector cero que no es un vector propio). Por lo tanto, lo que nos interesa es buscar un subespacio propio de A y despu´es escoger un vector propio particular; por ejemplo, un vector propio unitario. As´ı pues, debemos estudiar la convergencia de subespacios y no de vectores propiamente.
10.2.1.
Convergencia de subespacios
En primer lugar, al conjunto de subespacios de Fn de dimensi´on d se le llama Grassmanniana o Variedad de Grassman, y se representa mediante el s´ımbolo Grd pFn q: Grd pFn q tS Fn : S subespacio de dimensi´on du Para poder hablar de convergencia de subespacios, debemos dotar a Grd pFn q de una estructura topol´ogica. Hay dos formas cl´asicas de hacerlo que producen la misma topolog´ıa: definir una distancia entre subespacios, llamada distancia “gap”, o identificar el conjunto de subespacios de dimensi´on d con un conjunto cociente de matrices y definir en ´este la topolog´ıa cociente. Aunque la primera forma de proceder es m´as apropiada para el an´alisis num´erico porque proporciona cotas cuantitativas en el an´alisis del error, es tambi´en la m´as larga. Como, adem´as, no haremos un
´ 10.2. EL METODO DE LAS POTENCIAS
217
an´alisis del error muy riguroso, adoptaremos la segunda aproximaci´on, que es m´as r´apida. Dado un subespacio S de dimensi´on d en Fn y dada una base tx1 , . . . , xd u, llamaremos matriz base de S a la matriz X x1 xd . Las matrices base de los subespacios de dimensi´on d tienen la particularidad de ser matrices de tama˜ no n d y de rango d. Y toda matriz con estas propiedades genera un subespacio de dimensi´on d. Pero puede haber m´as de una matriz n d de rango d que genere el mismo subespacio: todas las matrices base deben estar relacionadas mediante una matriz de cambio de base. Es decir, X1 , X2 P Fnd generan el mismo subespacio de dimensi´on d si y s´olo si rang X1 rang X2 d y existe una matriz invertible P P Fdd tal que X1 X2 P . Ahora bien, la relaci´on X1 X2 si y s´olo si existe P P Fdd tal que X1 X2 P es una relaci´on de equivalencia en el conjunto Mn,d pFq de las matrices n d de rango d con elementos en F. Y cada subespacio de dimensi´on d puede identificarse con una clase de equivalencia. As´ı, identificaremos Grd pFn q M
pFq
Mn,d pFq Gld pFq
donde Gln,d representa el conjunto de clases de equivalencia por la relaci´on que d pFq acabamos de definir. En concecuencia, dada una matriz X P Mn,d pFq denotaremos por X ¡P Grd pFq el subespacio generado por las columnas de X. Ahora, Mn,d pFq es un conjunto abierto de Fnd que podemos considerar dotado de la topolog´ıa relativa. En este caso, U Mn,d pFq es abierto en Mn,d pFq si lo es en Fnd . Con esta topolog´ıa en Mn,d pFq podemos dotar a Grd pFn q de la topolog´ıa M pFq : Si cociente en Gln,d d pFq π : Mn,d pFq X
Ñ GrdpFnq Ñ X¡
es la proyecci´on can´onica, U Grd pFn q es abierto si y s´olo si π 1 pU q es abierto en Mn,d pFq. As´ı pues, la topolog´ıa de Grd pFn q es la topolog´ıa m´as fina que hace de π una aplicaci´on continua. Tenemos adem´as Lema 10.1 π es una aplicaci´on abierta. Demostraci´ on.- Hay que probar que si U Mn,d pFq es abierto entonces π pU q es abierto en Grd pFn q. Y, de acuerdo con la definici´on de la topolog´ıa definida en Grd pFn q, esto significa que π 1 pπ pU qq es abierto en Mn,d pFq. Ahora bien,
218
El problema de valores propios
Y P π 1 pπ pU qq si y s´olo si Y ¡P π pU q; y esto equivale a que alguna matriz X P U. Por lo tanto, π 1 pπ pU qq tXP : X
Y ¡ X ¡ para
P U y P P GldpFqu
Para demostrar que π 1 pπ pU qq es abierto tenemos que ver que si X P π 1 pπ pU qq existe un r ¡ 0 tal que si }X Y } r entonces Y P π 1 pπ pU qq. Tomamos, como norma, cualquier norma de matriz y sea X P π 1 pπ pU qq. Entonces XP P U para alguna matriz P P Gld pFq. Como U es abierto, existe s ¡ 0 tal que si }XP Z } s entonces Z P U. Sea s 0 r }P }
y supongamos }X Y } r. Como }XP Y P } ¤ }X Y }}P } resulta que }XP Y P } s; por lo que Y P P U. Es decir, Y P π 1 pπ pU qq tal y como se deseaba demostrar.
Hay algunas propiedades topol´ogicas importantes de Grd pFn q que utilizaremos pero que no vamos a demostrar. Son las siguientes:
n,dpFq es el subconjunto de Mn,dpFq de matrices cuyas columnas son (i) Si M ortonormales y Ud pFq representa el grupo matrices unitarias (ortogonales si n,dpFq{UdpFq son F R) de tama˜ no d d entonces Mn,d pFq{ Gld pFq y M n,dpFq est´a dotado de la topolog´ıa relativa homeomorfos. Se entiende que M n,dpFq{UdpFq es el conjunto cociente de como subconjunto de Mn,d pFq y M clases de matrices con columnas ortonormales que son matriz base del mismo n,dpFq est´an relacionadas si y s´olo si espacio vectorial. Es decir, Q1 , Q2 P M n,dpFq{UdpFq se existe U P Ud pFq tal que Q2 Q1 U . El espacio cociente M supone dotado con la correspondiente topolog´ıa cociente.
n,dpFq{UdpFq son homeomorfos se puede La prueba de que Mn,d pFq{ Gld pFq y M obtener como consecuencia de que los factores Q y R en la factorizaci´on QR de A (recordemos que por ser A de rango completo es u ´nica) son funciones continuas de A.
n,dpFq entonces M r1pπrpV qq tQZ : Q P V y Z P UdpFqu, π r es una aplicaci´on abierta. y se demuestra igual que en el Lema 10.1 que π
(ii) Si V
´ 10.2. EL METODO DE LAS POTENCIAS
219
(iii) Ya hemos dicho que a Grd pFn q se le puede dotar de una estructura de espacio m´etrico con la llamada m´etrica ”gap”. Somos ahora un poco m´as precisos en este punto: Dados dos subespacios S1 , S2 P Grd pFn q se define la distancia gap entre S1 y S2 como ΘpS1 , S2 q }P1 P2 }2 siendo Pi la proyecci´on ortogonal sobre Si . Esta definici´on es una generalizaci´on de un hecho bastante natural para subespacios de dimensi´on 1. Si dim S1 dim S2 1 entonces S1 x ¡ y S2 y ¡, digamos, con x 0, y 0. En este caso resulta (aunque no es trivial demostrarlo) que }P1 P2 }2 sen θ siendo θ el a´ngulo agudo que forman x, y. Es claro que en este caso ΘpS1 , S2 q sen θ define una m´etrica. En cualquier caso, Grd pFn q con la topolog´ıa derivada de la m´etrica gap y n,dpFq{UdpFq son homeomorfos. La los espacios cocientes Mn,d pFq{ Gld pFq y M demostraci´on de este resultado tampoco es trivial. (iv) Como todo espacio m´etrico es Hausdorff, Grd pFn q es Haussdorff.
n,dpFq es cerrado y acotado; por lo tanto com(v) Se demuestra f´acilmente que M pacto (por ser subespacio de Fnd y F R o C). Como la proyecci´on can´onica n,dpFq Ñ r:M π
pq
n p q es continua, Grd pF q es compacto.
Mn,d F Ud F
El siguiente resultado muestra que la convergencia de sucesiones de subespacios puede estudiarse a partir de la convergencia de matrices bases de los mismos. Teorema 10.2 Sea tSk uk0,1,2,... Grd pFn q una sucesi´on de subespacios y S P Grd pFq. Sea X P Mn,d pFq una matriz base de S y, para k 0, 1, 2, . . . Xk una matriz base de Sk . Entonces Sk Ñ S si y s´olo si para cada k 0, 1, 2, . . . existe una matriz invertible d d Pk tal que Xk Pk Ñ X. Demostraci´ on.- La demostraci´on es elemental. Supongamos Sk Ñ S y sea UX un entorno abierto de X. Como π es abierta VS π pUX q es un entorno abierto de S y como Sk Ñ S, existe un n´ umero natural N tal que si k ¥ N entonces Sk P VS . Dado que VS π pUX q, para cada k ¥ N existe una matriz Yk P UX tal que Sk Yk ¡. Es decir, para cada k ¥ N existe una matriz invertible Pk tal que Xk Pk P UX . Esto significa que Xk Pk Ñ X.
220
El problema de valores propios
El rec´ıproco se demuestra igual. Si VS es un entorno abierto de S entonces UX π 1 pVS q es un entorno abierto de X. Si existen matrices invertibles Pk tales que Xk Pk Ñ X entonces existe un entero positivo N ¡ 0 tal que si k ¥ N se tiene que Xk Pk P UX pero VS π pUX q. Por lo tanto, para k ¥ N se tiene que Sk π pXk Pk q P VS . Es decir, Sk Ñ S.
r es abierta, la demostraci´on del siguiente resultado es igual que la Usando que π del Teorema 10.2 Teorema 10.3 Sea tSk uk0,1,2,... Grd pFn q una sucesi´on de subespacios y S P n,dpFq una matriz base ortonormal de S y, para k 0, 1, 2, . . . Grd pFq. Sea Q P M Qk una matriz base ortonormal de Sk . Entonces Sk Ñ S si y s´olo si para cada k 0, 1, 2, . . . existe una matriz unitaria d d Zk tal que Qk Zk Ñ Q.
Con estos resultados sobre convergencia de subespacios podemos probar el siguiente resultado fundamental: Teorema 10.4 Sea X P Mn,d pFq y A P Fnn una matriz tal que rangpAk X q d para todo k 1, 2, . . .. Si la sucesi´on de subespacios t Ak X ¡u converge, lo hace a un subespacio S P Grd pFn q que es A-invariante. Demostraci´ on.- Supongamos que Ak X ¡Ñ S. Como rangpAk X q d para cada k ¥ 0 el subespacio Ak X ¡P Grd pFn q. Ahora bien, Grd pFn q es un espacio compacto y Hausdorff, por lo tanto cerrado. Si Ak X ¡ converge debe hacerlo a un subespacio de Grd pFn q. Veamos que, adem´as, ´este debe ser A-invariante. En efecto, como para k 0, 1, 2, . . . Ak X es una matriz base de Ak X ¡ aplicando el Teorema 10.2, para cada k 0, 1, 2, . . . existe una matriz invertible Pk tal que AK XPk Ñ Y , donde Y es una matriz base de S. Dado que A es lineal, y por lo tanto continua, ApAK XPk q Ñ AY . Aplicando de nuevo el Teorema 10.2 deducimos que Ak 1 X ¡Ñ AY ¡. Pero por hip´otesis Ak 1 X ¡Ñ Y ¡. Como Grd pFq es Haussdorf, el l´ımite es u ´nico. Por lo tanto, AY ¡ Y ¡ S. Ahora es f´acil deducir que AS S. En efecto, si x P S entonces x Y ξ y Ax AY ξ P AY ¡ S.
´ 10.2. EL METODO DE LAS POTENCIAS
10.2.2.
221
El algoritmo del m´ etodo de las potencias
Aplicamos el Teorema 10.4 al m´etodo de las potencias.
Corolario 10.5 Sea A P Fnn y x P Fn un vector no nulo. Si la sucesi´on de subespacios Ak x ¡ converge a un subespacio y ¡ no nulo entonces y es un vector propio de A.
Demostraci´ on.- En primer lugar, si Ak x ¡Ñ y ¡ con y 0 entonces A x 0 para todo k 1, 2, . . .. Por el Teorema 10.4 y ¡ es A-invariante, lo cual significa que Ay P y ¡. Es decir, Ay αy para alg´ un α P F. Por lo tanto, y es un vector propio de A asociado al valor propio α. k
N´otese que de acuerdo con la demostraci´on del corolario, si A y x son reales y
Ak x ¡ converge a un subespacio no nulo y ¡ entonces y es vector propio real asociado a un valor propio real.
De forma similar se demuestra
Corolario Sea A P Fnn y q P Fn un vector unitario. Si la sucesi´on de subes 10.6 Ak q pacios }Ak q}2 converge a un subespacio y ¡ no nulo entonces y es un vector propio unitario de A.
Este resultado es la base del m´etodo de las potencias cuya primera versi´on es la siguiente:
222
El problema de valores propios
M´ etodo de las Potencias (Primera Versi´ on) Dada A P Rnn Paso 1: El´ıjase x0 P Fn1 (F R o C) x0 Paso 2: q0 }x0}2 Paso 2: for j 0, 1, 2, 3, . . . hasta convergencia
Aqj xj 1 (vector propio unitario aproximado) 1 k xj 1 k ˜ j 1 q Aqj 1 (valor propio aproximado) λ j 1 xj qj
1
end for. Antes de analizar la convergencia de este algoritmo conviene hacer algunas observaciones: 1. En primer lugar, el algoritmo construye la sucesi´on de vectores unitarios: q0 , q1 , q2 , . . . donde Aqi1 qi }Aqi1}2 , i 1, 2, . . . Es muy f´acil ver que i
}AAiqq0} . 0 2 " k * A q0 Es decir, el algoritmo construye la sucesi´on }Ak q0}2 del Corolario 10.6. Por qi
lo tanto, si esta sucesi´on converge lo hace a un vector propio unitario de A. No obstante, esta sucesi´ ok n puede no converger y sin embargo s´ı hacerlo la sucesi´on A q0 de subespacios }Ak q0}2 . Esto es lo importante. En tal caso, de acuerdo con el Corolario 10.6, obtendremos un subespacio propio de A como l´ımite.
2. Se usa provisionalmente la expresi´on “hasta convergencia” para indicar que, en caso de que la sucesi´on de subespacios converja, la rutina for se realizar´ıa un n´ umero suficiente de veces. M´as adelante estudiaremos un criterio de terminaci´on de esta rutina.
´ 10.2. EL METODO DE LAS POTENCIAS
223
3. El algoritmo calcula vectores unitarios qj que cuando qj ¡Ñ q ¡ proporcionan un vector propio unitario, q, de A. Por eso llamamos vectores propios aproximados a los vectores qj . umero en F. Si qj fuera un 4. Con cada vector qj se calcula qj Aqj que es un n´ vector propio exacto entonces qj Aqj ser´ıa un valor propio exacto. En efecto: Aqj
λqj ñ qjAqj λqjqj λ
porque qj es unitario. Si qj es un vector pr´oximo a un vector propio entonces ˜ j q Aqj es, por continuidad, un n´ umero pr´oximo a un valor propio de A. λ j Por eso le llamamos ”valor propio aproximado”. Adem´as, si A es real y x0 es real, entonces todo se realiza en aritm´etica real. Es decir, los vectores y valores propios aproximados calculados son reales. 5. A´ un cuando la sucesi´on qj no converja, si qj ¡ converge lo hace a un subespacio propio y los vectores qj seguir´an siendo vectores propios aproximados. ˜ j converger´a a un valor proEn consecuencia, si qj ¡ converge, la sucesi´on λ pio de A. Es decir, si el algoritmo converge (en el sentido que qj ¡ converge) se obtiene un vector propio unitario y un valor propio tan aproximados como se quiera (dentro de la aritm´etica en punto flotante que se utilice) de A. Pero dado que los vectores propios pueden cambiar de direcci´on, los criterios para terminar el algoritmo deber´an tener en cuenta tanto a los vectores propios como a los valores propios.
10.2.3.
An´ alisis de la convergencia del m´ etodo de las potencias
La cuesti´on ahora es establecer condiciones suficientes lo m´as amplias posible sobre A y q para poder asegurar que qk converge. El siguiente teorema es el resultado fundamental: Teorema 10.7 Sea q0 P Fn un vector unitario y A P Fnn una matriz diagonalizable con λ1 , . . . , λn como valores propios. Sea tv1 , . . . , vn u una base de vectores propios unitarios de A (Avi λi vi y }vi }2 1, i 1, . . . , n) y supongamos que q0 α1 v1 αn vn . Supongamos tambi´en que
|λ1| ¡ |λ2| ¥ |λ3| ¥ ¥ |λn|
(10.1)
224
El problema de valores propios
Entonces qj ¡Ñ v1 ¡ si y s´olo si α1 0. M´as a´ un, en tal caso, para j 0, 1, 2, . . . existen n´ umeros complejos zj de m´odulo 1 tales que
j }zj qj v1}2 O λλ2 . 1
Recu´erdese que la expresi´on }zj qj v1 }2 O λλ
2 1
j
significa que para j sufi-
cientemente grande existe una constante positiva K tal que
j }zj qj v1}2 ¤ K λλ2 . 1
Por lo tanto, la segunda parte del Teorema nos dice que, cuando " *hay convergencia, ´esta se produce a la misma velocidad a la que la sucesi´on
converge a 0. λji vi y, en
0, como λ1 0 porque |λ1| ¡ |λ2| ¥ 0, podemos definir βj
} Aj q0 }2 .
Aj q0
Si α1
2 1
Demostraci´ on.- Es f´acil ver que si Avi consecuencia,
α1Aj v1
αn Aj vn
j λλ
λi vi entonces Aj vi
α1λj1v1
αn λjn vn .
j
0 y recordando que qj }AAj qq0} tenemos que 0 2 j j α2 λ2 αn λn βj qj v1 v2 vn . α1 λ1 α1 λ1 Como por hip´otesis λλ 1, tenemos que l´ım βj qj v1 . Por j Ñ8 concluimos que qj ¡Ñ v1 ¡.
α1 λj1
As´ı βj
2 1
(10.2) el Teorema 10.2
Rec´ıprocamente, supongamos qj ¡Ñ v1 ¡. Por el Teorema 10.2 existen n´ umeros βj tales que βj qj Ñ v1 . Sea P la proyecci´on sobre v1 ¡ a lo largo de v2, . . . , vn ¡. Entonces P pβj qj q Ñ P v1 v1. Pero P pβj qj q
βj βj j j j P p A q q 0 }Aj q0}2 }Aj q0}2 P pα1λ1v1 α2λ2v2
αn λjn vn q
βj j }Aj q0}2 α1λ1v1.
´ 10.2. EL METODO DE LAS POTENCIAS
225
βj α1 λj1 v1 Ñ v1 . Como v1 0 debe ser α1 0. A la proyecci´on sobre v1 ¡ j }A q0}2 a lo largo de v2 , . . . , vn ¡ se le llama proyecci´on espectral sobre el subespacio propio v1 ¡. As´ı
Para demostrar la segunda parte del teorema partimos de la expresi´on (10.2) y ponemos j j j α2 λ2 α3 λ3 α n λn xj v2 v3 vn . α1 λ1 α1 λ1 α 1 λ1 Entonces, teniendo en cuenta que }vi }2
1, j j j α λ α λ 2 2 3 3 }xj }2 ¤ α λ }v2}2 α λ }v3}2 ααn λλn }vn}2 1 j 1 1 j 1 j 1 1 λ α λ α λ2 α2 α3 λ3 ααn λλn 1 1 1 2 1 2
La sucesi´on
α2 α1
j α3 λ3 α1 λ2
j αn λn α1 λ2
es mon´otona decreciente (o mejor, no creciente) y acotada superiormente. Sea K una cota superior (por ejemplo, K puede tomarse el elemento de la sucesi´on correspondiente a j 0). As´ı j }xj }2 ¤ K λλ2 . 1 Por otra parte, como βj qj v1 xj y l´ım xj 0 tenemos que l´ım βj qj v1 . j
Ñ8
Adem´as, la norma es una funci´on continua, de modo que l´ım |βj | j
Ñ8
j
Ñ8
1 dado que
}qj }2 }v1}2 1. M´as a´un, j } A q 0 }2 ||βj | 1| j 1 1 j }Aj q0}2 |α1λj1| |α1λ1| |α1λ1| 1 j }A q0}2 }α1λj1v1}2 ¤ |α1λj1| ¤ 1 j }Aj q0 α1λj1v1}2 }xj }2. |α1λ1| As´ı, para j 0, 1, 2, . . . 1 }xj }2 ¤ |βj | ¤ 1 }xj }2 .
226
El problema de valores propios
y como xj Ñ 0, para j suficientemente grande 1 }xj }2 ¡ 0. Supondremos de ahora en adelante que j es suficientemente grande como para que se cumpla esta condici´on. β Ahora si zj |βjj | entonces |zj | 1 y zj qj
1 |βj | pv1
xj q ¤
1 pv1 1 } xj } 2
xj q
1
}xj }2 v 1 1 }xj }2
1 xj . 1 }xj }2
Por lo tanto, teniendo en cuenta que }v1 }2
1, }zj qj v1}2 ¤ 1 2}x}jx}2} j 2 Para j suficientemente grande 1 }xj }2 ¥ 12 , por lo que }zj qj v1}2 ¤ 4}xj }2. j λ2 Finalmente, como }xj }2 ¤ K , conclu´ımos que λ1 j }zj qj v1}2 ¤ 4K λλ2 1 para j suficientemente grande, tal y como se deseaba demostrar. La condici´on suficiente de convergencia del Teorema 10.7 requiere que la matriz A tenga un valor propio dominante en m´odulo y que para el vector inicial x0 la componente de ´este en la direcci´on del subespacio propio asociado a dicho valor propio dominante sea no nula. Ambas condiciones son gen´ericas (es decir, casi todas las matrices y vectores las cumplen). Ya sabemos que el conjunto de las matrices que tienen valores propios distintos es abierto y denso, de modo que la condici´on de tener un valor propio dominante en m´odulo es gen´erica salvo que la matriz sea real y tenga valores propios complejos. Esta situaci´on la analizamos enseguida. En el supuesto de que la matriz dada tenga un valor propio dominante en m´odulo, la elecci´on de un vector x0 que no tenga componente nula en la direcci´on del correspondiente subespacio propio no presenta problemas. Pr´acticamente cualquier vector cumple esta condici´on. Pensemos, por ejemplo, en la siguiente situaci´on: sabemos que alguien ha dibujado dos vectores linealmente independientes en R2 con origen en p0, 0q pero desconocemos qu´e vectores son. ¿Qu´e posibilidades tenemos de dar un vector que sea colineal con uno de ellos? Convendremos inmediatamente que es muy improbable
´ 10.2. EL METODO DE LAS POTENCIAS
227
que escogido un vector al azar, ´este sea colineal con uno de los vectores dados. De la misma forma, si escogemos x0 al azar es casi seguro que su componente en la direcci´on del subespacio propio asociado al valor propio dominante sea no nula. Una buena pr´actica consiste en NO tomar x0 un vector con estructura (por ejemplo, x0 p1, 0, . . . , 0q) sino al azar. Si la matriz tiene valores propios complejos a´ un cuando sus valores propios sean distintos pueden tener igual m´odulo. Esto es claro si la matriz es real con valores propios complejos, pero puede suceder, incluso, con matrices complejas. Por ejemplo, la matriz A = 1.0000 + 1.0000i 0 0 - 1.0000i
0 + 2.0000i 1.0000 - 1.0000i -3.0000
0 0 1.0000
tiene como valores propios: >> eig(A) ans = 1.0000 1.0000 + 1.0000i 1.0000 - 1.0000i Aplicando el m´etodo de las potencias con vector >> x0=rand(3,1) x0 = 0.8147 0.9058 0.1270 produce 4 subsucesiones de “valores propios” que parecen estabilizarse en los valores 1.0000 - 0.6119i, 1.5931 - 0.2490i, 1.0000 - 0.1563i, 0.4069 - 0.2490i
228
El problema de valores propios
Veamos que ´este no es un hecho casual. Supongamos que A es diagonalizable y que |λ1 | |λ2 | ¡ |λ3 | ¥ . Supongamos tambi´en que λ1 λ2 y que las componentes, α1 y α2 , de x en las direcciones de los vectores propios v1 y v2 no son cero. Supongamos finalmente que λ1 eiθ λ2 , 0 θ 2π. Entonces Ak x λk2
k
eikθ α1 v1
α2 v 2
λ3 λ2
α3
v3
k
αn
λn λ2
vn
Si existe un n´ umero racional p t{s tal que θ 2π entonces hay t subsucesiones de p ( k A x ¡ que convergen a los subespacios generados por los vectores e En el ejemplo de arriba 1
2ksπi t
α1 v1
ie
2πi 4
α2 v 2 ,
k
1, . . . , t.
p1 iq, de modo que θ 2π4 y t 4.
Hay otros casos en los que no cumpli´endose la condici´on (10.1) todav´ıa se puede asegurar la convergencia. Recordemos que A es diagonalizable y supongamos que λ1 es un valor propio semisimple pero λ1 λ2 λr ¡ |λr 1 | ¥ con r n o λ1 λn con |λn | ¡ 0 (el caso A diagonalizable y λi 0 para i 1, . . . , n implica que A 0). En estas condiciones
r ¸ k
Ak x λ1
i 1
αi vi
n ¸
i r 1
k
αi
λi λ1
vi
siendo x1 α1 v1 αr vr la proyecci´on de x sobre el subespacio v1 , . . . .vr ¡ a lo largo de vr 1 , . . . , vn ¡. Como el valor propio λ1 es semisimple de multiplicidad r, resulta que v1 , . . . .vr ¡ Kerpλ1 In Aq. Por lo tanto, si la proyecci´on de x sobre Kerpλ1 In Aq a lo largo de vr 1 , . . . , vn ¡ no es cero, para k 0, 1, 2, . . . 1 existen escalares βk k (n´otese que λ1 0 porque o bien |λ1 | ¡ ¡ |λn | ¥ 0 λ1
r λn con |λn| ¡ 0) tales que βk Ak x Ñ ° αivi. Es decir, Ak x ¡Ñ i1 °r °r α v ¡. Pero α v P Kerpλ I Aq, de modo que ´este es un vector propio de
o λ1
i 1
i i
i i
1 n
i 1
A asociado al valor propio λ1 . Un u ´ltimo aspecto a considerar es el de la velocidad de convergencia cuando ´esta se da. Tal y como comentamos antes de la demostraci´on del Teorema 10.7, la segunda parte de ´este nos indica que cuando se dan las condiciones de convergencia,
´ 10.2. EL METODO DE LAS POTENCIAS
229
la velocidad a la que los subespacios Aj x0 ¡ convergen al subespacio propio correspondiente al valor propio de mayor m´odulo es aproximadamente la misma λ2 j que la velocidad de convergencia de la sucesi´on . Si el cociente |λ2 |{|λ1 | es λ1 pr´oximo a 1 debemos esperar una convergencia lenta y cuanto m´as peque˜ no sea esta raz´on, la convergencia ser´a m´as r´apida. En cualquier caso se trata de una raz´on de convergencia lineal porque depende del cociente |λ2 |{|λ1 | y no de su cuadrado, o su cubo, etc. Cuando se pretende calcular un u ´nico valor y vector propio, el m´etodo de las potencia puede ser un primer algoritmo para conseguir en unas pocas iteraciones un vector y un valor propio aproximados que puedan servir como valores iniciales para otros m´etodos que convergen m´as r´apidamente.
10.2.4.
Criterios de terminaci´ on del algoritmo
Recordemos que al escribir el algoritmo fuimos muy imprecisos sobre cuando deber´ıa terminar. Dec´ıamos que se ejecutar´ıa un bucle para calcular los sucesivos vectores unitarios propios aproximados y los correspondientes valores propios aproximados “hasta la convergencia”(cuando la haya). Ahora vamos a ser m´as precisos. Un posible criterio de finalizaci´on del algoritmo podr´ıa ser cuando la distancia entre dos subespacios pr´oximos propios aproximados sea menor que una cantidad previamente especificada. Para ello bastar´ıa almacenar los vectores propios aproximados calculados en dos pasos sucesivos, digamos qj y qj 1 , y calcular la distancia entre ellos utilizando la m´etrica “gap”. Esto es sencillo porque sabemos que si Xj qj ¡ y Xj 1 qj 1 ¡ entonces ΘpXj , Xj
1
q }PX PX }2 j
j 1
siendo PX la proyecci´on ortogonal sobre X . Ahora bien como qj y qj PXj
qj qj
y PX j
1
qj
1
son unitarios
1 qj 1 .
As´ı pues, s´olo habr´ıa que modificar el algoritmo para que se ejecutara el c´alculo de qj hasta que }qj qj qj 1 qj 1 }2 fuera menor que una cantidad previamente especificada. Nosotros, sin embargo, utilizaremos otro criterio para terminar el algoritmo de las potencias y todos los que se derivan de ´el. Este criterio tiene su fundamento en el siguiente resultado:
230
El problema de valores propios
Proposici´ on 10.8 Sea A P Fnn , x P Cn un vector unitario y µ P C. Sea r Ax µx. Existe una matriz E P Cnn tal que k E kF }r}2 y pA E qx µx.
Podemos interpretar este resultado como una especie de “estabilidad hacia atr´as”. En efecto, esta proposici´on dice que los valores-vectores propios aproximados de A son valores-vectores propios exactos de matrices pr´oximas a A. La demostraci´on se deja como ejercicio.
Debe observarse que el error relativo de tomar A
E por A es
}E }F }r}2 }A}F }A}F El c´alculo del residuo r es muy f´acil en el proceso que define el algoritmo de las }F potencias y consideraremos obtenida la convergencia cuando el error relativo }}E A}F sea menor que una cantidad predeterminada. Es decir, cuando hayamos calculado r tan pr´oxima a A como necesitemos. un valor-vector propio exactos de una matriz A Dado que puede no haber convergencia, habr´a que tomar las precauciones necesarias por si el error relativo nunca alcanza a ser menor que la cantidad prefijada. En conclusi´on, fijaremos un ¡ 0 y un n´ umero m´aximo de iteraciones y el algoritmo terminar´a cuando el error relativo sea menor que o cuando se sobrepase el n´ umero de iteraciones m´aximo. Con todo esto en mente, la versi´on final del algoritmo de las potencias es la siguiente:
´ ´ INVERSA 10.3. EL METODO DE ITERACION
231
M´ etodo de las Potencias (Versi´ on Final) Datos: A P Fnn , x0 P Fn , ¡ 0 e itermax Objetivo: calcular un par pλ, xq valor-vector propio aproximado de A. El algoritmo devolver´a, adem´as, el n´ umero de iteraciones iter realizado as´ı como el residuo en la u ´ltima iteraci´on. x
x0 }x0}2 res=1 iter=0 normA=}A}F while res ¡ normA e iter itermax y Ax (siguiente potencia) λ x y (valor propio aproximado) res=}y λx}2 (residuo) y (normalizaci´on) x }y}2 iter=iter+1 end while
10.3.
El m´ etodo de iteraci´ on inversa
El objetivo del m´etodo de iteraci´on inversa es mejorar la velocidad de convergencia del m´etodo de las potencias y posibilitar la obtenci´on de un valor-vector propio que no se tenga que corresponder necesariamente con el valor propio dominante. Se parte de la siguiente observaci´on: si λ0 es un valor propio de A P Fnn y ppλq pd λd pd1 λd1 p1 λ p0 es un polinomio cualquiera entonces ppλ0 q es un valor propio de la matriz ppAq pd Ad
pd1 Ad1
p1 A
p0 In .
En efecto, si v es un vector propio de A para λ0 entonces Av Por lo tanto ppAqv
pd A d v pdλd0 v ppλ0qv.
pd1 Ad1 v pd1 λd01 v
p1 Av p1 λ0 v
λ0v y Aj v λj0v.
p0 v p0 v
232
El problema de valores propios
Es decir, v es tambi´en vector propio de ppAq para ppλ0 q. Se prueba de forma similar que si ppAq es invertible entonces ppλ0 q1 es valor propio de ppAq1 y v es vector propio de esta matriz para dicho valor propio:
ppλ0qv ô v ppλ0qppAq1v ô ppλ0q1v ppAq1v. En particular si σ R ΛpAq entonces σIn A es invertible y λ0 P ΛpAq si y s´olo si 1 P ΛppσIn Aq1q. Por lo tanto, si λ1,. . . ,λn son los valores propios de A σ λ0 entonces los valores propios de pσIn Aq1 son ppAqv
µ1
σ 1 λ
1
, . . . µn
σ 1 λ
n
El m´etodo de iteraci´on inversa es el m´etodo de las potencias aplicado a la matriz
pσIn Aq1 para alg´un n´umero σ. A este n´umero se le suele llamar desplazamiento
o precondicionamiento y se suele usar para aumentar la velocidad de convergencia del m´etodo de las potencias y para obtener un valor-propio particular de la matriz A cuando se tiene un cierto conocimiento de ellos. Veamos como se hace. Supongamos r0 λi de A (por que por el medio que sea conocemos un valor propio aproximado λ ejemplo con una o dos cifras decimales de precisi´on) pero λi no tiene por qu´e ser el r0 de modo que σ λi y σ es muy valor propio dominante de A. Tomamos σ λ diferente de los dem´as valores propios de A. Esto significa que
1 1 |µ1| σ λ ¡¡ µj σ λ i j
(¡¡ significa ”mucho mayor”) para todo valor propio λj de A distinto de λi . Si aplicamos el m´etodo de las potencias a pσIn Aq1 tenemos que si |µ2 | ¥ |µ3 | ¥ k ¥ |µn| entonces la velocidad de convergencia ser´ıa del orden de µµ2 , tanto 1 m´as grande cuanto m´as pr´oximo est´e σ a λi . En estas condiciones el m´etodo de las potencias aplicado directamente a pσIn Aq1 proporcionar´ıa, supuesto que converja, un valor-vector propio de esta matriz; digamos, µ1 y v1 . Pero, tal y como hemos visto m´as arriba v1 ser´ıa tambi´en un vector propio de A para λi y podemos usar este vector propio para calcular directamente el valor propio de A: λi v1 Av1 . En otras palabras, modificaremos el m´etodo de las potencias para, haciendo uso de los vectores propios unitarios aproximados qj , obtener valores propios aproximados de A en la forma qj Aqj .
´ ´ INVERSA 10.3. EL METODO DE ITERACION
233
Hay otra modificaci´on en la implementaci´on del m´etodo de iteraci´on inversa respecto del de las potencias. El primer paso en el algoritmo del m´etodo de las potencias consiste en calcular qj 1 Aqj . Si lo aplicamos directamente a la matriz pσIn Aq1 habr´ıa que hacer qj 1 pσIn Aq1qj . En vez de calcular la inversa, lo que se hace es resolver el sistema
pσIn Aqqj 1 qj . Ambas operaciones son te´oricamente equivalentes pero la segunda requiere menos operaciones. De hecho, gen´ericamente σIn A admite una descomposici´on LU que conviene calcular de una vez por todas y usar s´olo sustitusci´on hacia delante y hacia atr´as para resolver el sistema pσIn Aqqj 1 qj . En nuestra implementaci´on del algoritmo, y puesto que practicaremos con MATLAB, dejaremos que sea ´este quien se encargue de resolver este sistema. Para ello usaremos la orden qj 1 pσIn Aqzqj . Se podr´ıa argumentar que en cualquier caso la mariz σIn A est´a tanto m´as cerca de ser singular cuanto m´as cerca σ de un valor propio exacto de A, y que esto hace que los errores de redondeo se magnifiquen tanto si se calcula la inversa como si se resuelve el sistema. En efecto es as´ı y a pesar de todo, el algoritmo converge. ¿Por qu´e? La respuesta a esta cuesti´on est´a en el an´alisis del error en el algoritmo que se emplee para resolver los sistemas linales. Si este algoritmo es estable hacia atr´as (como por ejemplo lo es para casi todas las matrices el de la eliminaci´on gaussiana), la soluci´on calculada para el sistema (pσIn Aqy x es la soluci´on exacta de alg´ un sistema pσIn A E qy x conde el error relativo }E }F {}A}F es del orden del epsilon de la m´aquina. Si el valor-vector propio (“exacto”) que estamos buscando est´a bien condicionado, ser´a muy parecido al de pσIn A E q y la sucesi´on de valores-vectors propios calculados seguir´a convergiendo hacia una valor-vector propio de σIn A. Dicho de otra manera m´as intuitiva, la mayor parte de las veces la propagaci´on del error juega a nuestro favor porque cuando σ est´a muy cerca de un valor propio de A, el vector y se calcula, en efecto, de manera inexacta pero, tal y como vimos en la demostraci´on del Teorema 10.7,
pσIn Aqk q0 pσ α1λ qk v1 pσ α2λ qk v2 pσ αnλ qk vn. 1 2 n Esto significa que cuanto m´as pr´oximo est´e σ a λ1 mayor ser´a la componente de pσIn Aqk q0 en la direcci´on de v1, de manera que los errores en el redondeo redundan
234
El problema de valores propios
en una potenciaci´on de la direcci´on en la que queremos que converjan los subespacios propios aproximados. De hecho, si σ y λ1 son casi iguales, una sola iteraci´on nos proporciona un vector que es ya casi un vector propio porque la proyecci´on espectral de ´este sobre el subespacio v2 , . . . , vn ¡ es insignificante en comparaci´on con la proyecci´on espectral sobre v1 ¡. En resumen, el algoritmo para el m´etodo de iteraci´on simult´anea es el siguiente:
M´ etodo de Iteraci´ on Inversa (Primera Versi´ on) Datos: A P Fnn , x0 P Fn , σ, ¡ 0 e itermax Objetivo: calcular un par pλ, xq valor-vector propio aproximado de A. x
x0 }x0}2 res=1 iter=0 normA=}A}F while res ¡ normA e iter itermax y pσIn Aqzx (ô pA σIn qy x) y x }y }2 (vector propio unitario aproximado) λ x Ax (valor propio aproximado) res=}λx Ax}2 (nuevo residuo) iter=iter+1 end while
Hay algunas mejoras que se pueden introducir en el algoritmo a fin optimizar el n´ umero de operaciones evitando sucesivas multiplicaciones de A por x, que es la parte m´as costosa. Concretamente, en la notaci´on del algoritmo representamos con el mismo s´ımbolo x dos n´ umeros diferentes en las dos primeras l´ıneas del bucle while. Vamos a mantener el s´ımbolo x para el primero; i.e. y
p su actualizaci´on: y con x
pσIn Aq1x p x
y }y}2 .
´ ´ INVERSA 10.3. EL METODO DE ITERACION
De esta forma
pσIn Aqxp pσIn}y} Aqy }yx} 2
Pongamos w
235
}yx}
.
2
. Entonces
2
pAxp σ xppσIn Aqxp σ ρ λx con
ppσIn Aqxp xpw. ρx
Finalmente el residuo es r
λxp Axp pσIn Aqxp ρxp w ρxp.
Con todo esto en mente el algoritmo de iteraci´on inversa queda como sigue:
M´ etodo de Iteraci´ on Inversa (Versi´ on Final) nn n Datos: A P F , x0 P F , σ, ¡ 0 e itermax Objetivo: calcular un par pλ, xq valor-vector propio aproximado de A. x
x0 }x0}2 res=1 iter=0 normA=}A}F while res ¡ normA e iter y pσIn Aqzx x w }y}2 y x }y}2 ρ x w λσρ res=}w ρx}2 iter=iter+1 end while
itermax
236
10.4.
El problema de valores propios
El m´ etodo del cociente de Rayleigh
El m´etodo del cociente de Rayleigh es exactamente igual que el de iteraci´on inversa pero con un desplazamiento que var´ıa en cada iteraci´on. La cuesti´on es ¿c´omo elegir el desplazamiento para hacer que el residuo sea lo m´as peque˜ no posible en cada iteraci´on? La respuesta es sencilla porque se trata de calcular la soluci´on de un problema de m´ınimos cuadrados lineal. Recordemos que el residuo es r
}λx Ax}
donde λ y x son el valor y vector unitario propios aproximados que se calculan en cada iteraci´on. Puesto que λ es un valor propio aproximado (que est´a variando en cada iteraci´on) bien podemos hacer uso de ´el en cada iteraci´on. La pregunta es ¿cu´al es el valor de λ que hace m´ınimo el residuo? Nos preguntamos por el problema de hallar m´ın }αx Ax}2
P
α F
Debe observarse que en este problema la matriz de coeficientes es x (n 1) y el vector de t´erminos independientes es Ax (tambi´en n 1). Recordemos que una forma resolver el problema de m´ınimos cuadrados m´ınyPFn }By c}2 con B P Fmn y c P Fm1 es resolviendo el sistema de ecuaciones normales: B By0 B c. En nuestro caso el sistema de ecuaciones normales es x xα0 Por lo tanto α0
xAx
xxAx x
es el n´ umero que minimiza el residuo.
Definici´ on 10.9 Para A n´ umero
P Fnn se llama cociente de Rayleigh de A en x P Fn al Rpxq
x Ax x x
As´ı pues, lo que hemos demostrado m´as arriba es
´ 10.4. EL METODO DEL COCIENTE DE RAYLEIGH
237
Teorema 10.10 Si A P Fnn y x P Fn entonces el cociente de Rayleigh de A en x es el n´ umero donde se alcanza el m´ınimo del residuo }αx Ax}2 . Es decir m´ın }αx Ax}2
P
α F
}Rpxqx Ax}2.
En la pr´actica, el algoritmo del cociente de Rayleigh consiste en sustituir σ por λ en el algoritmo de iteraci´on inversa iniciando el primer desplazamiento con el cociente de Rayleigh de la matriz dada en el vector inicial elegido. N´otese que puesto que los vectores propios aproximados son unitarios, el denominador en el cociente de Rayleigh es siempre 1.
M´ etodo del cociente de Rayleigh Datos: A P Fnn , x0 P Fn , ¡ 0 e itermax Objetivo: calcular un par pλ, xq valor-vector propio aproximado de A. x
x0 }x0}2 λ x Ax res=1 iter=0 normA=}A}F while res ¡ normA e iter y pλIn Aqzx x w }y}2 y x }y}2 ρ x w λλρ res=}w ρx}2 iter=iter+1 end while
itermax
Aunque en esta forma de implementar el algoritmo no lo parezca, resulta que λ Rpxq en cada iteraci´on. Para verlo basta deshacer los cambios que hicimos en el algoritmo de iteraci´on inversa. En efecto, al igual que en aquel caso denotamos con
238
El problema de valores propios
x p y λp los valores actualizados y con x y λ los que entran en la iteraci´on. Debemos p Rpxpq xpAxp. Para ello, en la iteraci´on λp λ ρ. Entonces ver que λ
p λ
y x p w λ }y} }y} λρλx 2 2 λy y y Ay λ xp Axp. y y yy
λ
y pλIn Aqy yy
En los ejercicios veremos que, en general, la convergencia del m´etodo del cociente de Rayleigh sigue siendo lineal, pero si la matriz es sim´etrica es, al menos, cuadr´atica. Ello es debido al siguiente resultado, que se conoce con el nombre de propiedad estacionaria de los vectores propios de las matrices sim´etricas. Teorema 10.11 Sea A P Rnn una matriz sim´etrica. Entonces pλ0 , x0 q es un valorvector propio de A si y s´olo si x0 es un punto estacionario o cr´ıtico de Rpxq y Rpx0 q λ0 . Esto es, ∇Rpx0 q 0 y Rpx0 q λ0 donde ∇Rpxq es el gradiente de Rpxq.
Demostraci´ on: Debemos demostrar que pλ0 , x0 q es un valor-vector propio de xT Ax B R p x0 q 0, i 1, . . . , n y Rpx0 q λ0 . En el caso real Rpxq T . A si y s´olo si B xi x x Bp xT Axq T Calculamos primero Bxi y luego haremos lo mismo con x x. En primer lugar x Ax T
n ¸
ajk xj xk ,
j,k 1
as´ı que
n n ¸ BpxT Axq ¸ aik xk aji xj . B xi j 1 k1 Como A es sim´etrica, aji aij y entonces n BpxT Axq 2 ¸ aij xj 2pAxqi B xi j 1
´ 10.4. EL METODO DEL COCIENTE DE RAYLEIGH
239
siendo pAxqi la i-´esima componente del vector Ax. De la misma forma xT x
°n
x2j y
j 1
BpxT xq 2x . i B xi Por lo tanto
BR 2pxT xqpAxqi 2pxT Axqxi 2 pAxq xT Ax x i i B xi pxT xq2 xT x xT x xT2x ppAxqi Rpxqxiq xT2x pAx Rpxqxqi
donde pAx Rpxqxqi q es la i-´esima componente del vector Ax Rpxqx. Ahora es claro que ∇Rpx0 q 0 y Rpx0 q λ0 si y s´olo si Ax0 λ0 x0 0; i.e., pλ0 , x0 q es un valor-vector propio de A. Tal y como hab´ıamos anunciado una consecuencia interesante de este Teorema es
Corolario 10.12 Si para una matriz sim´etrica A y un vector q0 el algoritmo del cociente de Rayleigh converge, lo hace, a partir de un cierto momento, cuadr´aticamente.
Demostraci´ on: Ya sabemos que si el algoritmo del cociente de Rayleigh converge, lo hace a un valor-vector unitario propios, pλ0 , x0 q, de A. Es decir, el algoritmo produce una secuencia pλk , αk qk q de valores-vectores propios aproximados de A que converge a pλ0 , x0 q. Pero λk Rpαk qk q y desarrolando Rpxq en serie de Taylor alrededor de x0 tenemos que Rpxq Rpx0 q ∇Rpx0 qpx x0 q
Op}x x0 }2 q Op}x x0 }2 q
por ser x0 un vector propio de A. Esto implica que para k suficientemente grande Rpαk qk q Rpx0 q Op}qk x0 }2 q. Como qk ¡ converge a x0 ¡ linealmente (al menos), el residuo converge a cero, a partir de ese k suficientemente grande, cuadr´aticamente (al menos).
240
El problema de valores propios
10.5.
El algoritmo QR con desplazamiento expl´ıcito
Todos los algoritmos tienen como objetivo u ´ltimo reducir una matriz dada A P
F , a una forma de Schur (real o compleja) para extraer de ah´ı la informaci´on n n
correspondiente a los valores y vectores propios. La finalidad de los que hemos visto hasta ahora era obtener un vector propio unitario (siempre que se den las condiciones apropiadas) y a partir de ah´ı, usando el cociente de Rayleigh, el correspondiente valor propio. Este es el primer paso para obtener la forma de Schur compleja porque, tal y como mostramos en la demostraci´on del Teorema de Schur, una vez obtenidos, r , para digamos pλ1 , q1 q, se puede construir una matriz unitaria, digamos Q q Q la que λ 1 b Q AQ . 0 B A continuaci´on se aplica el algoritmo correspondiente a B y as´ı sucesivamente. Este proceso de ir aplicando el algoritmo a matrices cada vez m´as peque˜ nas se llama deflaci´ on. Si se dan las condiciones apropiadas (valores propios con m´odulos diferentes y proyecci´on espectral no nula sobre los correspondientes subespacios propios) el resultado de este proceso de deflaci´on es una matriz triangular superior unitariamente semejante a A; i.e., una forma de Schur de A. El algoritmo QR comparte el mismo objetivo: reducir A a una forma de Schur. La forma en que este proceso se lleva a cabo, sin embargo, parece no tener nada que ver con los algoritmos que hemos visto hasta ahora. Esta es la primera versi´on del algoritmo QR con desplazamiento expl´ıcito:
Dada A P Fnn A0 A for k 0, 1, 2, . . . hasta convergencia El´ıjase un desplazamiento σk Calc´ ulese la descomposici´on QR de Ak σk In : Ak σk In Revi´ertanse los factortes: Ak 1 Rk Qk σk In end for
Qk Rk
En resumen, el algoritmo construye una sucesi´on de matrices tAk u a partir de
10.5. EL ALGORITMO QR CON DESPLAZAMIENTO EXPL´ICITO
241
una dada A A0 calculando la factorizaci´on QR de Ak σk In (para un σk que se elige expl´ıtamente en cada paso) y definiendo la siguiente matriz de la sucesi´on, Ak 1 , multiplicando los factores Q y R al rev´es y a˜ nadi´endo el desplazamiento. ¿Por qu´e hacer tal cosa produce algo significativo en relaci´on con el c´alculo de los valores y vectores propios de A? La historia empieza con John Francis y Eva Kublanovskaya que lo desarrollaron independientemente hacia 1959 en Inglaterra y la extinta Uni´on Sovi´etica. Se tard´o muchos a˜ nos, sin embargo, en comprenderlo bien. La raz´on de por qu´e este algoritmo calcula valores y vectores propios de A es, esencialmente, porque es una forma elegante de implementar el algoritmo del cociente de Rayleigh y, cuando se toma desplazamiento σk 0 (algoritmo QR a secas, sin desplazamiento), entonces calcula las mismas matrices que el algoritmo de iteraci´on simult´anea (u ortogonal) que es una generalizaci´on del algoritmo de las potencias. En este curso s´olo veremos la relaci´on entre el algoritmo QR con desplazamiento expl´ıcito y el del cociente de Rayleigh. Pero antes de nada es importante observar que todas las matrices en la sucesi´on tAk u que construye el algoritmo son unitariamente semejantes. Lema 10.13 Las matrices Ak que construye el algoritmo QR con desplazamiento expl´ıcito son todas unitariamente semejantes. Demostraci´ on: Para un k 0, 1, 2, . . . arbitario, supongamos dada Ak y veamos que Ak 1 es unitariamente semejante a Ak . Esto demuestra que todas las matrices de la sucesi´on son unitariamente semejantes. Sea Ak σk In Qk Rk la factorizaci´on QR de Ak σk In . Ak
1
por lo que Ak
Rk Qk σk In Qk pAk σk In qQk σk In Qk Ak Qk σk Qk Qk σk In Qk Ak Qk , 1
pRk Qk pAk σk Inq
y Ak son unitariamente semejantes.
1, 2, . . . Ak Qk1 Ak1 Qk1 Qk1 Qk2 Ak2 Qk2 Qk1 Qk1 Q0 A0 Q0 Qk1 Este Lema tiene una consecuencia importante: Para cada k
Es decir, el algoritmo QR no s´olo proporciona una sucesi´on de matrices unitariamente semejantes a A sino que nos da tambi´en la matriz unitaria que conecta A con cada una de las matrices de la sucesi´on: es el producto de los factores Q en la factorizaci´on QR de Ak σk In . Si para alg´ un valor de k la matriz Ak es (casi) una forma de Schur de A tendremos adem´as la matriz unitaria que convierte A en dicha (casi) forma de Schur.
242
10.5.1.
El problema de valores propios
El algoritmo QR y el cociente de Rayleigh
Para ver que el algoritmo QR converge en las mismas condiciones que el algoritmo del cociente de Rayleigh vamos a ver que ambos son esencialmente lo mismo. En realidad el algoritmo QR aplica el m´etodo del cociente de Rayleigh para calcular vectores propios por la izquierda en vez de hacerlo para calcular vectores propios por la derecha. ¿C´omo ser´ıa el m´etodo del cociente de Rayleigh en este caso? Igual que el ya visto pero operando con vectores filas: Dada A, un vector fila inicial x0 el algoritmo del cociente de Rayleigh consta de 4 pasos para j 1, 2, . . .: 1. Elegir un desplazamiento σj , que es el cociente de Rayleigh de la iteraci´on anterior, 2. Resolver el sistema yj pA σj In q xj 1 , 3. Normalizar xj
}yyj}
(vector propio aproximado por la izquierda), y
j 2
4. Calcular el nuevo cociente de Rayleigh λj ximado)
xj Axj (nuevo valor propio apro-
Cuando el algoritmo converge lo hace a un valor-vector propio unitario por la iz- p q quierda pλ, q q de A. Con este q podemos formar una matriz unitaria Q Q de forma que (v´ease la Secci´on ?? del Cap´ıtulo ?? sobre vectores propios por la izquierda) B b Q AQ 0 λ Ahora proceder´ıamos por deflaci´on aplicando el algoritmo por filas a B y as´ı sucesivamente. Supongamos ahora que hacemos una sola iteraci´on del algoritmo del cociente de Rayleigh por filas con vector inicial x0 en . El desplazamiento inicial ser´ıa en Aen
σ.
Sea q el vector unitario que obtenemos al aplicar esta iteraci´on (quitamos los sub´ındices porque s´olo hablamos de una iteraci´on y por comodidad): q
σIn q1 }eenppAAσI q1} n
n
(10.3) 2
10.5. EL ALGORITMO QR CON DESPLAZAMIENTO EXPL´ICITO
243
Analicemos una iteraci´on del algoritmo QR con el mismo desplazamiento σ en Aen . Calculamos la descomposici´on QR de A σIn : A σIn Entonces y
Q
QR.
RpA σInq1
en Q
enRpA σInq1 rnnenpA σInq1. ´ltima fila de Q . Denot´emosla por qr . Este vector es unitario Ahora bien, en Q es la u y acabamos de ver que
qr
rnnenpA σInq1.
De aqu´ı deducimos (recordamos que los elementos diagonales de R son n´ umeros reales positivos) 1 rnn }enpA σInq1}2 , y en consecuencia, comparando con (10.3) q˜ q. En palabras: la u ´ltima fila de la matriz Q en la factorizaci´on QR de A σIn con σ en Aen coincide con el vector propio por la izquierda unitario aproximado, q, de A que calcula el m´etodo del cociente de Rayleigh cuando se emplea en como vector inicial en la primera iteraci´on . Es decir,
r q Q Q
(10.4)
r El valor propio aproximado que nos proporciona el cociente de Rayleigh ser´a λ q Aq y el residuo (recordemos que estamos calculando vectores propios por la izquierrq. Estos c´alculos concluyen la iteraci´on del m´etodo del cociente da) r q A λ de Rayleigh. Respecto al algoritmo QR con desplazamiento, la segunda y u ´ltima etapa consiste en calcular la siguiente matriz en la sucesi´on revirtiendo los factores Q y R de la factorizaci´on QR de A y sumando el desplazamiento:
p RQ A
σIn
QpA σInqQ
σIn
QAQ.
244
El problema de valores propios
Escribiendo Q como en (10.4)
p A
QrAQr QrAq r q r Q q AQ q Aq
r Q q
A
r y qA r λrq. Por ello, teniendo en cuenta que q es Ahora bien, q Aq λ r ortogonal a todas las columnas de Q, r rQr q AQ
rqQr rQr λ
y
(10.5) }qAQr}22 }rQr}22 rQrQrr. r}2 }r}2. Para ello debemos recordar Nuestro objetivo es demostrar que }q AQ r es la soluci´on del problema de el Teorema 10.10. De acuerdo con este Teorema, λ m´ınimos cuadrados consistente en calcular el n´ umero complejo, α, que hace m´ınima la distancia αq q A. Es decir,
}r}2 }λrq qA}2 m´ ın }αq q A}2 . αP C
Recordemos tambi´en (Teorema ?? del Cap´ıtulo 7) que para el problema de m´ınimos cuadrados m´ınn }By c}2 , y0 es una soluci´on si y s´olo si By c P pIm B qK . En nuestro
P
y C
rq qA P pIm qqK o equivalentemente qr 0. Es caso, debe suceder que r λ r Con esto en mente, podemos demostrar que }qAQr}2 }r}2. En decir, r P Im Q. efecto, por (10.5) }qAQr}22 rQrQrr, rQr es la proyecci´on ortogonal sobre Im Qr y r P Im Q; r por lo tanto QrQrr pero Q r}22 rr }r}22 tal y como dese´abamos demostrar. En definitiva y }q AQ
r QrAq p QrAQ A r λr q AQ
r
r}2 }r}2 siendo r qA λrq es el residuo del algoritmo del cociente y }q AQ de Rayleigh aplicado por filas a la matriz A con vector inicial en . En conclusi´on Teorema 10.14 Para una matriz A P Fnn el algoritmo QR con desplazamiento expl´ıto σk en Ak en produce sucesi´on de matrices tAk uk¥0 tales que si pa pkuna p kq q p kq p kq ra k 0, 1, 2, . . ., an an1 ann1 ann es la u ´ltima fila de Ak , Ak
10.5. EL ALGORITMO QR CON DESPLAZAMIENTO EXPL´ICITO
pk1q
245
pkq pk1q
Qk1 Ak1 Qk1 y qn es la u ´ltima fila de Qk1 entonces pann , qn q es el valorvector propio unitario por la izquierda aproximados que produce el algoritmo del cociente Rayleigh aplicado a Ak1 con vector inicial en . Adem´as la norma del de kq vector apn1 apnnkq1 es igual a la norma del residuo producido por el algoritmo del cociente de Rayleigh.
Este resultado nos indica que si el algoritmo del cociente de Rayleigh converge para la matriz A con vector inicial en entonces usando como desplazamiento expl´ıcito el elemento en la posici´on pn, nq de la matriz obtenida en cada iteraci´on, el algoritmo QR con estos desplazamientos produce una sucesi´on de matrices que converge a una matriz con la forma
B c 0 λ
siendo λ un valor propio de A. Adem´as , el vector formado por los elementos diagonales de la u ´ltima fila de la matriz producida en cada iteraci´on tiene la misma norma que el residuo producido por el algoritmo del cociente de Rayleigh. Podemos usar como criterio para terminar el proceso iterativo que la norma de este vector sea menor que una cantidad previamente especificada. Una vez obtenido un vector suficientemente peque˜ no sustituir´ıamos los elementos no diagonales de la u ´ltima fila por cero y proceder´ıamos a deflactar la matriz hasta obtener, si hay convergencia, una matriz en forma de Schur. El proceso, adem´as, produce las matrices unitarias de paso. El algoritmo en cada etapa ser´ıa el siguiente:
246
El problema de valores propios
Una etapa del algoritmo QR Datos: A P Fnn , ¡ 0 e itermax Objetivo: Obtener una matriz T , semejante a A, cuyos elementos en la u ´ltima fila sean todos cero excepto el de la diagonal, y la matriz unitaria de semejanza Q. normA=}A}F T A, Q In d T pn, 1 : n 1q, iter=0 while }d}2 ¡ normA e iter itermax s T pn, nq (desplazamiento) T sIn Q1 R1 (factorizaci´on QR de T sIn ) T R1 Q1 sIn (revertimos los factores) d T pn, 1 : n 1q (nuevos elementos no diagonales) Q QQ1 (actualizamos la matriz unitaria de paso) iter=iter+1 end while if iter itermax T pn, 1 : n 1q 0 end if
10.5.2.
(para que quede bonito)
An´ alisis de la velocidad de convergencia del algoritmo QR
Para analizar la velocidad de convergencia del algoritmo QR vamos a determinar pk 1q pk q una cota de }an }2 en t´erminos de }an }2 . Para no arrastrar los sub´ındices todo p Ak 1. La notaci´on ser´a similar al referirnos el rato, vamos a poner A Ak y A a los diversos elementos y matrices relacionados con Ak y Ak 1 . As´ı, pondremos p σIn RQ. A σIn QR y A Consideremos A, Q y R particionadas como sigue: A σIn Entonces
B σIn1 h g µσ
P f S r e π 0 ρ QR.
p ph S r P f B σIn1 p A σIn p σ 0 ρ e π . gp µ
(10.6)
(10.7)
10.5. EL ALGORITMO QR CON DESPLAZAMIENTO EXPL´ICITO
247
El objetivo es acotar gp en funci´on de g. En primer lugar, como Q es unitaria, la norma de cada una de sus filas y columnas es 1. En particular, }e}22 |π |2 }f }22 |π |2 , de donde concluimos que
}e}2 }f }2. Por otra parte, en (10.6)
g
eS. Suponiendo que S es no singular y δ }S 1 }2 , tenemos que }e}2 ¤ δ}g}2. Ahora como Q es unitaria, despejando R en (10.6) obtenemos
S r 0 ρ
n1 Pf π¯e B gσI
h
µσ
.
De aqu´ı sacamos que ρ f h π ¯ pµ σ q. Como |π | ¤ 1, }e}2 y, por la desigualdad de Cauchy-Schwarz, |f h| ¤ }f }2 }h}2 :
}f }2, }e}2 ¤ δ}g}2
|ρ| ¤ δ}g}2}h}2 |µ σ|. Finalmente de (10.7) obtenemos gp ρe . Poniendo todo junto y teniendo en cuenta, otra vez, que }e}2 ¤ δ }g }2 conclu´ımos }gp}2 ¤ δ2}h}2}g}22 δ|µ σ|}g}2, o, con los sub´ındices restaurados,
}gk 1}2 ¤ δk2}hk }2}gk }22
δk |µk σk |}gk }2 .
Notemos finalmente que para todo k 1, 2, . . . }hk }2 ¤ }A}2 . En efecto, todas las matrices Ak que se van construyendo mediante el algoritmo QR son unitariamente semejantes (Ak 1 Qk Ak Qk ), de modo que }Ak }2 }A}2 . Pero
}hk }22 ¤ }hk }22 |µ|2 }Ak en}22 ¤ }Ak }22}en}22 }Ak }22. En conclusi´on, existe un n´ umero real η ¡ 0 tal que }hk }2 ¤ η para todo k 1, 2, . . ..
Entonces
}gk 1}2 ¤ ηδk2}gk }22
δk |µk σk |}gk }2 .
(10.8)
248
El problema de valores propios
Lo primero que sugiere esta expresi´on es, tal y como ya hemos deducido al comparar los algoritmos QR y cociente de Rayleigh, que el desplazamiento apropiado en cada iteraci´on es σk µk porque esta elecci´on proporciona una cota de gk 1 en t´erminos del cuadrado de }gk }2 . En tal caso, la cota queda
}gk 1}2 ¤ ηδk2}gk }22 Vamos a suponer que existe un n´ umero real δ ¡ 0 tal que para todo k 1, 2, . . . }Sk1}2 ¤ δ. Se puede demostrar que, de hecho, esta propiedad es verdadera para valores de gk suficientemente peque˜ nos. En estas condiciones
}gk 1}2 ¤ δ2η}gk }22.
As´ı pues, cuando el algoritmo QR converge a un valor propio simple y }Sk1 }2 ¤ δ para todo k, la norma de gk 1 est´a acotada por una cantidad que es proporcional al cuadrado de la norma de gk . Se dice que la sucesi´on t}gk }2 u converge cuadr´aticamente, al menos, a cero. Tal y como ya hemos dicho la condici´on }Sk1 }2 ¤ δ se cumple para valores de gk suficientemente peque˜ nos, as´ı que no podemos decir que la convergencia es siempre cuadr´atica. Al principio puede no serlo, pero llegar´a un momento en el que lo ser´a con seguridad. Se dice que la convergencia del algoritmo QR es localmente cuadr´ atica. La convergencia cuadr´atica es muy conveniente porque, una vez que empieza, es muy r´apida. Para hacernos una idea, si δ 2 η 1 y }g0 }2 101 , entonces
}g1}2 ¤ 102 }g2}2 ¤ 104 }g3}2 ¤ 108 }g4}2 ¤ 1016 Cuatro iteraciones bastan para reducir el error por debajo de la unidad de redondeo para la doble precisi´on. Si A fuera herm´ıtica, entonces tambi´en lo ser´ıan cada una de las Ak (Recordemos que Ak Qk1 Ak1 Qk1 ). En particular tendr´ıamos que hk gk , de modo que
10.6. CONSIDERACIONES FINALES
249
podemos reemplazar η por }gk }2 y la acotaci´on ser´ıa
}gk 1}2 ¤ δ2}gk }32. Este tipo de convergencia se llama c´ ubica. Observaciones 10.15 El an´alisis de la convergencia que hemos hecho es local; esto es, se supone que gk es suficientemente peque˜ no (o, equivalentemente, µk suficientemente pr´oximo a un valor propio). No hay teoremas para la convergencia global del algoritmo QR aplicables a cualquier matriz. El comportamiento t´ıpico del algoritmo parece ser el siguiente: se consumen unas cuantas iteraciones en las que la convergencia puede ser muy lenta o parecer, incluso, que no hay convergencia y en un momento dado emp`ıeza la convergencia de manera muy r´apida. Para el c´alculo de los siguientes valores propios se necesitan cada vez menos iteraciones. La raz´on de esta u ´ltima caracter´ıstica del comportamiento del algoritmo QR se puede explicar, al menos parcialmente, por la relaci´on entre este algoritmo y el del m´etodo de iteraci´on ortogonal o simult´anea que no abordaremos en este curso.
10.6.
Consideraciones finales
Hay muchas m´as cosas que decir sobre el algoritmo QR. En particular, su relaci´on con el algoritmo de iteraci´on simult´anea, que es una generalizaci´on del m´etodo de las potencias, es muy importante porque sirve para establecer condiciones suficientes de convergencia. Otro aspecto importante es la relaci´on entre el algoritmo QR y la reducci´on a forma Hessenberg y, en general, la implementaci´on pr´actica del algoritmo QR (algoritmo QR con desplazamiento impl´ıcito). Tambi´en es importante la elecci´on de los desplazamientos y la reducci´on a forma de Schur real cuando la matriz es real. Gen´ericamente las matrices complejas tienen valores propios distintos en m´odulo, pero esto no es verdad para las matrices reales en las que cuando aparecen valores propios complejos ´estos vienen conjugados a pares. Es por eso que la reducci´on a forma de Schur real es tan importante para estas matrices. De vez en cuando surgen tipos de matrices para las que el algoritmo QR con los desplazamientos conocidos no converge. Lo habitual en tales casos es a˜ nadir nuevos desplazamientos a la lista de los ya conocidos. Sin embargo, por ahora no hay un algoritmo universal eficiente para todas las matrices.
250
El problema de valores propios
Finalmente, los algoritmos para el c´alculo de los valores propios no terminan con el algoritmo QR aunque ´este sea el algoritmo de prop´osito general que m´as se utiliza. Para matrices con estructura existen algoritmos especiales. Por ejemplo, para matrices sim´etricas o herm´ıticas que son unitariamente diagonalizables adem´as del algoritmo QR existen otros que aprovechan la estructura de estas matrices. Otro caso es el de las matrices con muchos elementos iguales a cero (sparse matrices en ingl´es). Los algoritmos para el c´alculo de valores propios han sido y sigue siendo un campo de mucha actividad investigadora.
Lihat lebih banyak...
Comentários
