Universidad Politécnica de Madrid–Escuela Técnica Superior de Ingenieros Industriales Grado en Ingeniería en Tecnologías Industriales. Curso 2015-2016-3º Matemáticas de Especialidad–Ingeniería Eléctrica
Valores y vectores propios Valores singulares
José Luis de la Fuente O’Connor
[email protected] [email protected] Clase_valvec_2016.pdf
1/95
Índice
2/95
Cuestiones teóricas de valores y vectores propios
Localización de valores propios
Obtención numérica de valores y vectores propios
Método de Jacobi Método de la iteración de potencia Método de la iteración inversa Iteración mediante cociente de Rayleigh Deflación Iteración simultánea Iteración QR Subespacios de Krylov Comparación de los métodos
Cálculo de valores singulares
3/95
Los valores y vectores propios adquieren una relevancia destacada para analizar asuntos con oscilación y resonancia. Su conocimiento es básico en:
Sistemas eléctricos de corriente alterna.
Modos de vibración natural de estructuras.
Instrumentos musicales.
Mecánica cuántica.
Lasers.
Resonancia Magnética Nuclear (NMR).
:::
4/95
Su cálculo e interpretación es esencial para el análisis de sistemas de generación, transporte y demanda de energía eléctrica. Frequency range 1
Frequency range 2
electromagnetic
Phenomena
electromechanical
SVC
Frequency Interconnected System
G
1 MHz
1 kHz
1 Hz
1 µs
1 ms
1s
1 µs M Turbine Syngovernor chronous Torsion maschine AVR
On-load Current Breaker tap transformer Arc changer Saturation restriking Protective system
External system Line Cable HVDC FACTS
Converter VAR-Load VAR-Admittance Motor Load torque
100 µs
Travelling Transient wave phenomena, phenomena Switching overvoltages
10 ms
1 mHz
1 min
1s
Short-circuit Transient phenomena stability SSRphenomena
Harmonics, Transformer-saturation
1h
100 s Control phenomena with steam generators
Processtimes Calculation time steps
5/95 From the time domain subprogram of
NETOMAC®
9
NETOMAC eigenvalue calculation mode
System’s working point at t0, t>t0
Augmented state equations:
Linear model of a dynamic power system
x
Axx
Axz
x
bx
0
Azx
Azz
z
bz
Sparsitybased storage and matrix computation
u
State equation: A
x
3
(Dimentioned maximum 800 dynamic order)
Eigensystem solution
Network in RST Admittances by differential equations non-linearities
Single line network Complex admittances only fundamental frequency
Time range Instantaneous values ns...μs...ms...s
Time range Quasi steady-state values s...min
Electromagnetic and electromechanical phenomena, complete solution
Electromechanical phenomena, fundamental frequency
only Loadflow
Power method of implicit inverse iteration
QR transformation
Simulation Models for System Components, Machines, Controllers and Control Units
Loadflow Initial conditions
x
Possible ways of simulation
Eigenvalues and associated eigenvectors
Loadflow Operating point
Overview of modes
System components Linearization Coupling
Activities of different modes on one device
Activities of one mode on different devices
R
Loadflow for special requirements, e.g. homogeneous multiconductor system
Transfer function based dominant pole method
L
P
U1
U
U2
1
Frequency range all system-variables
Eigenvalue analysis
Small-signal characteristics Network, machines and control
Systemoscillation and -damping, Netreduction, Controller layout
Eigenvalue analysis & transfer function analysis
2
3
Legend: R: Right eigenvector (observability information) L: Left eigenvector (controllability information) P: Participation factors U: Transfer function residues U1: Mode activities in unit impulse response U2: Mode activities in unit step response
y(t): Unit impulse/step response
G(jw): Frequency response
4
SIGMA OMEGA (rad/sec)
1
SIGMA OMEGA (RAD/SEC) (HZ) -0,087 4,566
COMPLEX S-PLANE
ZETA (%)
17 out of 206 solved modes are selected all the selected modes are inside the displayed range
FREQ (Hz)
ELECTRIC MACHINES
14
1
-4,053
12,991
-29,8
2,060
NORMALIZED RIGHT EIGENVECTOR (MODE SHAPE) NORMALIZED RIGHT EIGENVECTOR (MODE SHAPE)
NORMALIZED RIGHT EIGENVECTOR (MODE SHAPE)
MODE DISTRIBUTION ON COMPLEX S-PLANE SELECTED MODES MACHINE SWING MODES
Hz
13
2
SYNCHRONOUS MACHINES 18 items COMPONENT NAME
ZETA
SIGMA OMEGA SIGMA OMEGA (RAD/SEC) (RAD/SEC) (HZ) (HZ) -0,087 4,566 -0,087 4,566
FACO (%) -1,9
ZETA ZETA
FREQUENCY RESPONSE OF TRANSFER FUNCTION
FACO FACO (%) (%) 0(s)=
-1,217
11,754
-10,3
1,071
-2,520
11,450
-21,5
1,022
12
3
4
-0,053
10,022
-0,5
1,595
-1,110
9,072
-11,2
1,571
COMPONENT NAME COMPONENT NAME
-13,6
1,494
-12,6
1,490
8
-1,992
9,233
-21,1
1,469
9
-1,285
9,388
-1,104
-0,900
9,359
9,146
-9,0
1,456
10
10
1 1 DACOLOUA 2X150 MVA DACOLOUA 2X150 MVA
2 2 MILLOA 2x2x 1818 MVA MILLOA MVA
-2,079
8,942
3
LOSQUI 2X 18+14+Adong
4
TEHUSEQUEL+JALA+BANELVOL
5
FALALFAL 2X95 MVA
6
TOLIZASAV, ZALSAV 10+3X32
7
PELRA 5X76 MVA
8
NASVENTA AT 18 KV
1 1,407
-1,207
8,700
-14,5
-0,550
8,440
9
NASVENTA AT 19.2 KV
10
1,399
5
-6,6
1,343
4
13
-1,137
7,062
-14,3
1,251
-0,503
7,033
-7,1
1,119
6,310
-6,1
1,004
5,718
-4,9
0,910
17
-0,007
4,566
-1,09
0,727
0,19 0,19
CAREN 2X58.8 MVA
0,45 0,45
NASVENTA 8 8 NASVENTA ATAT 1818 KVKV
16
1
0,43 0,43
NASVENTA 19.2 9 9 NASVENTA ATAT 19.2 KVKV
2
3
0,18 0,18
CAREN 2X58.8 MVA 1010CAREN 2X58.8 MVA 8
11
ENCHEPHEHU 2X283 MVA
9
ENCHEPHEHU 2X283 MVA 1111ENCHEPHEHU 2X283 MVA
0,06 0,06
BUNCOL 2X240 MVA 1212BUNCOL 2X240 MVA
0,05 0,05
13
CURAMACHI 2X53 MVA
CURAMACHI 2X53 MVA 1313CURAMACHI 2X53 MVA
0,04 0,04
14
COANTU 2X160 MVA
12
2 -0,383
-0,200
0,16 0,16
0,15 0,15
PELRA 5X76 MVA 7 7 PELRA 5X76 MVA 17
BUNCOL 2X240 MVA
0,5
3
15
0,14 0,14
5 5 FALALFAL 2X95 MVA FALALFAL 2X95 MVA
6 6 TOLIZASAV, ZALSAV 10+3X32 TOLIZASAV, ZALSAV 10+3X32
18
14
16
COANTU 2X160 MVA 1414COANTU 2X160 MVA
0,06 0,06
15
ROTOREL 4X105 MVA
ROTOREL 4X105 MVA 1515ROTOREL 4X105 MVA
0,06 0,06
16
NICOREA 5X21.5 MVA
NICOREA 5X21.5 MVA 1616NICOREA 5X21.5 MVA
0,22 0,22
1
MAGNITUDE SCALE : 6.31 (x18-4)
0
0
-6
-5
-4
-3
-2
-1
0
17
QUEPAN 500 MVA
QUEPAN 500 MVA 1717QUEPAN 500 MVA
18
TILLARNUCA 2X68 MVA
TILLARNUCA 2X68 MVA 1818TILLARNUCA 2X68 MVA
0,37 0,37
0,56 0,56
1 0 0
12
0,91 Hz
0,37 0,37
3 3 LOSQUI 2X2X 18+14+Adong LOSQUI 18+14+Adong
4 4 TEHUSEQUEL+JALA+BANELVOL TEHUSEQUEL+JALA+BANELVOL
8
6 -22,9
11
pusMVA
1,00 1,00
0,63 0,63
9
7
12
ROTOR SPEED
NASVENTA AT 18 KV
MECHANICAL TORQUE
NYQUIST DIAGRAMM
DACOLOUA 2X150 MVA
MILLOA 2x 18 MVA
0,63 0,63
1,5 6
7
NASVENTA AT 16 KV PV
V(s)
4
SYNCHRONOUS MACHINES SYNCHRONOUS MACHINES 1818 items items
1
2
11
5
Y(s)
-1,9 -1,9
ELECTRIC MACHINES ELECTRIC MACHINES
MODE OBSERVABILITY OF DEVIATION VARIABLES MOTOR SPEED PY ParLub the dynamic system and observe the mode different devices.
2 2
3
Mode distribution (Eigenvalues)
Phasor-Diagram
0,10,1
0,20,2
0,30,3
Bar diagram
0,40,4
0,50,5
0,60,6
0,70,7
0,80,8
0,90,9
1 1
Nyquist-Diagram
6/95
Caso histórico
Hundimiento del Puente Tacoma 1, Washington, EE.UU.
Hundimiento del Puente Tacoma 2, Washington, EE.UU.
Hoy
Aspectos teóricos y algunas propiedades de los valores y vectores propios
7/95
En general, los vectores propios de un operador matemático1 lineal T son los vectores no nulos que cuando son transformados por el operador dan lugar a un múltiplo escalar de sí mismos2: T .v/ D v. Ese escalar, , se denomina valor propio.
La formulación de su cálculo en el caso habitual de que ese operador lo caracterice una matriz, Dada una matriz A 2 Cnn, calcular un escalar y un x ¤ 0 tales que Ax D x: El escalar es un valor propio de A y x su correspondiente vector propio. 1
H
Transformación lineal u aplicación lineal de un espacio vectorial V en si mismo. Ej. la ecuación de Schrödinger E D E E. 2 No cambian su dirección.
Para que exista una solución distinta de la trivial x D 0 el valor propio deberá ser raíz del polinomio característico de grado n asociado a A, es decir,
8/95
Lo que es igual a n
g1n
1
det.A
I/ D 0:
C g 2 n
2
C . 1/ngn D 0:
Igual que cualquier matriz tiene asociado un polinomio característico, cualquier polinomio tiene asociado una matriz compañera. La matriz compañera de un polinomio mónico p.t / D c0 C c1t C C cn 1t n 1 C t n es 20 0 : : : 0 c 3 C .p/ D
610 4: ::
0 ::: 0 1 ::: 0 :: : : :: : : : 0 0 ::: 1
0
c1 c2 :: :
cn
1
7 5:
Los valores propios de esta matriz C .p/ son las raíces del polinomio p.t /.
El polinomio mínimo q.t / de la matriz A es el polinomio mónico único de grado mínimo tal que q.A/ D 0.
9/95
Los vectores propios de A pertenecen al subespacio nulo de Ax I, ker.Ax I/, y no están unívocamente determinados: Si v es un vector propio, ˛v también lo es.
Siempre existen n valores propios de A 2 Cnn, reales o complejos. No siempre existen n vectores propios.
La multiplicidad algebraica del valor propio de A es la multiplicidad de la raíz correspondiente del polinomio característico asociado a A.
La multiplicidad geométrica de es el número de vectores propios linealmente independientes que se corresponden con . Teorema La multiplicidad geométrica de un valor propio es menor o igual que su multiplicidad algebraica.
Por ejemplo, si A D I, D 1 es un valor propio con multiplicidad algebraica y geométrica n. El polinomio característico de A es p.z/ D .z 1/n y e i 2 Cn, i D 1; : : : ; n, sus vectores propios.
Si el valor propio tiene una multiplicidad geométrica menor que la algebraica, se dice defectuoso. Se dice que una matriz es defectuosa si tiene al menos un valor propio defectuoso. La matriz # " 210 A D 021 002
10/95
tiene un valor propio de multiplicidad algebraica 3 y multiplicidad geométrica 1. >> A=[2 1 0;0 2 1;0 0 2]; >> [V D]=eig(A) V = 1.0000 -1.0000 1.0000 0 0.0000 -0.0000 0 0 0.0000 D = 2 0 0 0 2 0 0 0 2
Si una matriz A 2 Cnn no es defectuosa dispone de n vectores propios linealmente independientes.
Si V 2 Cnn tiene como vectores columna los n vectores propios de A y D D diag.1; : : : ; n/, entonces Avi D i vi , i D 1; : : : ; n, es equivalente a AV D V D. Si A no es defectuosa, A D V DV
1
que es la descomposición en valores propios de A: o A diagonalizada por V . Definición Una matriz A se dice diagonalizable por semejanza si es semejante a una matriz diagonal.
Teorema Una matriz A 2 Cnn es diagonalizable si y sólo si tiene n vectores propios linealmente independientes.
11/95
12/95
Definición Dos matrices A; B 2 Cnn se dicen semejantes si existe una matriz inver-
tible P tal que A D P 1 BP.
Teorema Dos matrices semejantes tienen el mismo polinomio característico y, por consiguiente, los mismos valores propios.
Definición El espectro .A/ de A es el conjunto de sus valores propios: .A/ D f 2 C W det.A
I/ D 0g:
Definición El radio espectral, .A/, de la matriz A es el valor máximo de los módulos de sus valores propios: .A/ D max.i 2.A / ji j: Es el radio del menor círculo del plano complejo centrado en el origen que contiene todos sus valores propios.
1
1:1
1:3
1:7
2:5
4:1
7:3
13:7
26:5 13/95
Al aplicársele a cualquier vector la transformación que representa A, ese vector 4 7KHYHFWRUV [ [ [ DUHVKRZQLQ)LJ 7KHRWKHUYHFWRUVDUHJURZLQJ A ; : : : ; A 3 tiende a orientarse+RZHYHU hacia la OLQHVHJPHQWVDUHGUDZQVKRZLQJWKHGLUHFWLRQVRIWKRVH dirección del vector propio dominante de A. WRRORQJWRGLVSOD\
YHFWRUV ,QIDFW WKHGLUHFWLRQVRIWKHYHFWRUVDUHZKDWZHUHDOO\ZDQWWRVHH QRWWKHYHF 1;8 0;8 WRUVWKHPVHOYHV 7KHOLQHVVHHPWREHDSSURDFKLQJWKHOLQHUHSUHVHQWLQJWKHHLJHQVSDFH Ejemplo Estudiemos la matriz A D 0;2 1;2 cuyos valores propios son VSDQQHGE\ Y1 0RUHSUHFLVHO\ WKHDQJOHEHWZHHQWKHOLQH VXEVSDFH GHWHUPLQHGE\ 4 A1k [DDQGWKHOLQH HLJHQVSDFH GHWHUPLQHGE\ 2 y 2 D 1. El vector propio correspondiente a es Si se multiplica 1 k! 1 .1 Y1 JRHVWR]HURDV
A y sus potencias por el vector x D
0;5 1
el resultado es el de la figura.
x2
A4 x Ax x
A3x
A2 x 1
Espacio de
v1
v1 1
4
10
x1
'*(63& 'LUHFWLRQVGHWHUPLQHGE\ [ A[ A2 [; : : : ; A7 [ Si ese vector está en la dirección de alguno de los vectores propios de A, se expande o contrae por un factor que determina el correspondiente valor propio 3
!k
k
14/95
Otros ejemplos Ejemplo del apéndice "Definiciones, : : :
10 1 0 Si A D , 1 D 1; x 1 D ; 2 D 2 y x 2 D : 02 0 1
1;5 0;5 1 1 Si A D , 1 D 2; x 1 D ; 2 D 1 y x 2 D : 0;5 1;5 1 1 Si A D
01 1 , 1 D i; x 1 D ; 2 D 10 i
p i i y x2 D , donde i D 1. 1
15/95
Teorema La suma de los valores propios de A es igual a su traza: 1 C 2 C C n D
n X
ai i .
i D1
Teorema El producto de los valores propios de A es igual a su determinante.
Teorema Los valores propios de una matriz triangular son los coeficientes de su diagonal principal.
Teorema Una matriz es singular si y sólo si tiene un valor propio igual a 0.
Teorema Si los valores propios de una matriz A son i , 1 i n, los de A
son i
˛; 1 i n. Sus vectores propios son idénticos.
˛I
16/95
Teorema Los valores propios de las potencias de A son las potencias de los de A; los vectores propios son los mismos.
I Demostración. Si consideramos la definición, Ax D xI A 2 x D Ax D 2 xI A n x D n x S 1 AS D I S 1 AS S 1 AS D 2 ) S 1 A 2 S D 2 : Las potencias negativas también: Ax D xI A 1 Ax D A 1 x ) A 1 x D 1 x:
Corolario Los valores propios de A
1
son los recíprocos de los de A.
17/95
Los valores propios de multiplicidad m > 1 tienen un subespacio propio asociado de dimensión m. Ese subespacio es invariante4 por lo que:
4
Todos los vectores de un subespacio propio son vectores propios.
Los subespacios propios correspondientes a valores propios distintos sólo tienen en común el vector nulo.
Si para todo vector w 2 W se cumple que f .w/ 2 W .
El problema de calcular los valores y vectores propios es en general un problema no lineal en los valores propios y lineal en los vectores x.
El cálculo de los valores propios por las raíces del polinomio característico no es aconsejable por:
El trabajo de determinar los coeficientes y raíces del polinomio. La sensibilidad a errores de redondeo de los coeficientes del polinomio (recordemos los polinomios de Wilkinson).
Ejemplo Consideremos la matriz p que maq: K .
1 1 , donde es cualquier número menor
Los valores propios exactos de A son 1 C y 1 . El cálculo mediante el polinomio característico haría det.A
I/ D 2
2 C .1
2 / D 2
2 C 1I
las raíces serían (valores propios calculados) 1 y 1.
u
18/95
Valores propios de matrices destacadas
19/95
Ortogonales
Recordemos que estas matrices son aquellas que tienen como inversa su traspuesta: QQT D I.
Ejemplo
Todos los valores propios de una matriz ortogonal tienen módulo unidad.
p p 2=2 2=2 01 1 0 p p : ; ; 10 0 1 2=2 2=2
Hermíticas
Son aquellas matrices iguales a su traspuesta conjugada: A D AH : Son una generalización de las matrices simétricas al campo complejo.
Ejemplo La matriz
1 1
1Ci : i 2
20/95
Si A es hermítica, el producto x Ax es un número real. H
Los valores propios de una matriz hermítica, en consecuencia, son números reales. En efecto Ax D xI
x Ax D x ˜ H
H
R
R
2
2
˜
x D jjxjj2
En una matriz hermítica los vectores propios correspondientes a dos valores propios distintos son ortogonales entre sí. En efecto, ) H H H x 2 A x 1 D 1 x 2 x 1 Ax 1 D 1x 1 H . /x 1 2 2 x 1 D 0: H H H Ax 2 D 2x 2 x A x 1 D 2 x x 1 2
2
Como los valores propios 1 y 2 son distintos, xH 2 x 1 D 0: Si los vectores propios se normalizan, x H x D 1, la matriz de vectores propios se convierte en una matriz ortogonal.
21/95
Unitarias
Son matrices cuya inversa es su compleja conjugada: U H U D U U H D I:
Ejemplo La matriz
p p i p2=2 p2=2 : 2=2 i 2=2 Las matrices unitarias son una extensión de las matrices ortogonales al campo complejo. Todos los valores propios tienen módulo unidad. Una matriz unitaria no modifica ni los ángulos ni las normas:
.U x/H .U y/ D x H U H U y D x H y si y D x; jjU xjj2 D jjxjj2:
Normales 2Son 3las matrices que cumplen que AAH D AH A:
120 Ejemplo 40 1 25 ; 201
i 0 : 0 3 5i
22/95
Triangularización de Schur
Issai Schur, Alemania, 1875-1941.
Teorema Triangularización de Schur Para cualquier A 2 Cnn existe una matriz unitaria U y una triangular superior T tales que U H AU D T : Los valores propios de A son entonces los coeficientes de la diagonal principal de T .
Teorema Para cualquier matriz hermítica A 2 Cnn existe una unitaria U tal que donde D es una matriz diagonal.
U H AU D D,
1. Los valores propios de A son números reales. 2. Se pueden obtener vectores propios de A que sean ortonormales.
Corolario Si A 2 Rnn es simétrica, existe una matriz ortogonal Q y una diagonal D tales queQT AQ D D.
Teorema Los valores propios de una matriz hermítica definida positiva son todos positivos. Recíprocamente, si todos los valores propios de una matriz son positivos, debe ser definida positiva.
Teorema Forma canónica de Jordan Para una matriz A 2 Cnn existe una matriz T
regular tal que
T
1
2
6 AT D J D 6 4
J1
0
:::
3
0
Jn
7 7; 5
2 3 i 1 6 i 1 0 7 Ji D 6 7 4 5 0 1 i es una matriz de Jordan y los i son los valores propios de A. donde
Marie Ennemond Camille Jordan, Francia, 1838-1922.
23/95
24/95
Tabla que resume las posibles transformaciones por semejanza. A
T
Valores propios distintos Simétrica real Hermítica compleja Normal Real cualquiera Cualquiera Cualquiera
Regular Ortogonal Unitaria Unitaria Ortogonal Unitaria Regular
BDT
1
AT
Diagonal Diagonal real Diagonal real Diagonal Triangular en bloques (real) Triangular superior (Schur) Casi diagonal (Jordan)
25/95
Localización de valores propios
Si no se necesita calcular el valor numérico exacto de los valores propios, sino conocer grosso modo dónde se encuentran en el plano complejo, existen varias formas de hacerlo.
La más simple surge de tomar normas en la expresión Av D v: jjjjvjj D jjvjj D jjAvjj jjAjjjjvjj ) jj jjAjj; para cualquier norma matricial inducida por una norma vectorial. Por consiguiente:
Los valores propios de una matriz se localizan en el plano complejo, dentro del circulo centrado en el origen de radio jjAjj.
26/95
Semyon Aranovich Ger˘sgorin, Rusia, 1901-1933.
Teorema Ger˘sgorin Los valores propios de una matriz A 2 Cnn se encuentran en la
unión de los n discos de Gershgorin, cada uno de los cuales está centrado en akk , k D 1; : : : ; n, y tiene de radio n X jakj j rk D j D1
j ¤k
I Demostración. Sea un valor propio de A y x su vector propio asociado. De Ax D x y .I A/x D 0 se tiene que n X . akk /xk D akj xj ; k D 1; : : : ; n; j D1
j ¤k
donde xk es el componente k-ésimo del vector x. Si xi es el coeficiente de x más grande en valor absoluto, como jxj j=jxi j 1 para j ¤ i, se tiene que j
ai i j
n X j D1
j ¤i
Luego está contenido en el disco f W j
jxj j X jaij j: jxi j n
jaij j
ai i j ri g.
j D1
j ¤i
27/95
El siguiente programa de Matlab calcula los círculos o discos de Gershgorin y los dibuja. function C = Gershgorin(A) % % Se dibujan los círculos de Gershgorin de la matriz A. [m n] = size(A); d = diag(A); cx = real(d); cy = imag(d); B = A - diag(d); r = sum(abs(B’)); % Suma filas de A sin diagonal C = [cx cy r(:)]; t = 0:pi/100:2*pi; c = cos(t); s = sin(t); [v d] = eig(A); % eig calcula los valores propios de A d = diag(d); % En d valores propios u1 = real(d); v1 = imag(d); hold on, grid on, axis equal xlabel(’Re’), ylabel(’Im’) h1_line = plot(u1,v1,’or’); set(h1_line,’LineWidth’,1.5) for i=1:n % Se dibujan los círculos de Gershgorin x = zeros(1,length(t)); y = zeros(1,length(t)); x = cx(i)+r(i)*c; y = cy(i)+r(i)*s; h2_line = plot(x,y); set(h2_line,’LineWidth’,1.2) end hold off title(’Círculos de Gershgorin y valores propios de la matriz’) end
Ejemplo Calcular los discos de Gershgorin de la matriz 2 3 123 A D 43 4 95 : 111
28/95
Los radios de los discos son r1 D 5 r2 D 12 r3 D 2:
Los valores propios son: 7;3067; 0;6533 C 0;3473i y Con el programa anterior, se obtiene lo siguiente. Círculos de Gershgorin y valores propios de la matriz 10
5
Im
0
−5
−10 −10
−5
0
5 Re
10
15
0;6533
0;3473i .
29/95
Los gráficos de los discos de Gershgorin de otras matrices curiosas se pueden ver en esta otra gráfica. gersh(gallery(’lesp’,12))
gersh(gallery(’hanowa’,10))
20 5 10
0
0
−10 −5 −20 −40
−30
−20
−10
0
−5
gersh(gallery(’ipjfact’,6,1))
0
5
gersh(gallery(’smoke’,16,1))
0.5 2 1 0
0 −1 −2
−0.5 −0.2
0
0.2
0.4
0.6
0.8
−2
−1
0
1
2
30/95
Índice
Introducción teórica
Localización de valores propios
Obtención numérica de los valores y vectores propios
Método de Jacobi
Método de la iteración de potencia
Método de la iteración inversa
Iteración mediante cociente de Rayleigh
Deflación
Iteración simultánea
Iteración QR
Subespacios de Krylov
Cálculo de valores singulares
31/95
Obtención numérica Método de Jacobi
Carl Gustav Jacobi, Alemania (Prusia), 1804-1851.
Es un método formulado en 1846 para el cálculo de los valores y vectores propios de una matriz simétrica o compleja hermítica.
Utiliza transformaciones por semejanza basadas en rotaciones, idénticas a las de Givens, para hacer cero pares de elementos simétricamente dispuestos respecto a la diagonal principal.
32/95
Partiendo de A 0 D A, cada iteración conforma una transformación A kC1 D J Tk A k J k ; h i c s donde cada matriz J k D s c se calcula de tal manera que h ih i h i h c 2a 2csa C s2a i apq .c 2 s 2 / C cs.app aqq / pp pq qq c s app apq c s s c s c D a .c 2 s 2 / C cs.a apq aqq a / c 2 a C 2csa C s 2 a pq
pp
qq
qq
pq
pp
sea diagonal.
Para lograrlo, apq .c 2 D
s 2/ C cs.app
aqq app 2apq
aqq / ha de ser cero. Haciendo
y t D s=c;
tangente del ángulo de rotación
se obtiene la ecuación de segundo grado t 2 C 2 t
1 D 0:
De las dos posibles raíces, tD
˙
p
33/95
1 C 2;
se escoge p la más pequeña para que j j =4. Luego se obtienen c D 1= 1 C t 2 y s D c t . function J=Jac_Rot(A,b,d) % Cálculo de la rotación de Jacobi para anular un coef. de A de coordenadas (b,d) if A(b,d)~=0 tau=(A(d,d)-A(b,b))/2/A(b,d); if tau>=0 t=1/(tau+sqrt(1+tau^2)); else t=-1/(-tau+sqrt(1+tau^2)); end c=1/sqrt(1+t^2); s=c*t; else c=1; s=0; end J=[c s; -s c]; % Igual que Givens end
En Matlab:
Mediante unos “barridos” que apliquen sistemáticamente estas transformaciones a todos los coeficientes que no estén en la diagonal principal de la matriz que tengan un valor mayor que una tolerancia se conseguirá ir convirtiendo la matriz en una diagonal.
34/95
La convergencia del proceso es cuadrática.
v u n n uX X u El proceso termina cuando off .A/ D u aij2 t i D1 j D1 j ¤i
h
102 021 211
> t ol kAkF .
i . Apliquémosle el método de Jacobi.
Ejemplo Sea A 0 D
Anulemos para empezar los coeficientes .1; 3/ y .3; 1/. Con ese fin, definamos la rotación 2 3 y hagamos
0;707 0 J0 D 4 0 1 0;707 0 2
0;707 0 5 0;707
3 3 0;707 0 A 1 D J T0 A 0 J 0 D 40;707 2 0;7075 : 0 0;707 1
35/95
A continuación, anulemos .1; 2/ y .2; 1/ mediante la rotación 2
y hagamos
0;888 J 1 D 40;460 0
3 0;460 0 0;888 05 0 1
2 3 3;366 0 0;325 A 2 D J T1 A 1 J 1 D 4 0 1;634 0;6285 : 0;325 0;628 1
Luego anulamos el .3; 2/ y el .2; 3/ usando la rotación 2
y haciendo
1 0 J 2 D 40 0;975 0 0;226
3 0 0;2265 0;975
3 3;366 0;0735 0;317 0 5: A 3 D J T2 A 2 J 2 D 40;0735 1;780 0;317 0 1;145 2
36/95
Comenzando un nuevo “barrido”, hagamos cero los coeficientes .1; 3/ y .3; 1/. Usaremos la rotación 2 3 y luego
0;998 0 J3 D 4 0 1 0;070 0 2
0;070 0 5 0;998
3 3;388 0;0733 0 A 4 D J T3 A 3 J 3 D 40;0733 1;780 0;0051 5 : 0 0;0051 1;167
El proceso continuaría hasta llegar a conseguir una aproximación a los valores propios deseados.
u
37/95
El programa para poder calcular una rotación de Jacobi 2 2 y aplicarla luego a la matriz original completa, pre y post multiplicándola, es este. function J=jacrot(A,i,j) % Transf. de Jacobi del coeficiente (i,j) y (j,i) de A n=length(A); J1=Jac_Rot(A,i,j); J=eye(n); % Calcula qué rotación 2x2 elemental aplicar J([i j],[i j])=J1([1 2],[1 2]); % Se adapta a la propia A end function J=Jac_Rot(A,b,d) % Cálculo de rotación de Jacobi 2x2 en la matriz A. if A(b,d)~=0 tau=(A(d,d)-A(b,b))/2/A(b,d); if tau>=0 t=1/(tau+sqrt(1+tau^2)); else t=-1/(-tau+sqrt(1+tau^2)); end c=1/sqrt(1+t^2); s=c*t; else c=1; s=0; end J=[c s; -s c]; % Igual que Givens end
Ejemplo con Matlab
>> A=rand(4); >> A=A+A’ A = 1.7818 1.1086 1.3615 0.3352 1.1086 0.5150 1.0842 0.5054 1.3615 1.0842 1.8585 0.9660 0.3352 0.5054 0.9660 0.9466 >> % Anulemos el coeficiente (1,2) >> J=jacrot(A,1,2); A=J’*A*J A = 2.4252 0.0000 1.7218 0.5436 0.0000 -0.1284 0.2544 0.2688 1.7218 0.2544 1.8585 0.9660 0.5436 0.2688 0.9660 0.9466 >> % Ahora el (1,3) >> J=jacrot(A,1,3); A=J’*A*J A = 3.8868 0.1646 0.0000 1.0396 0.1646 -0.1284 0.1939 0.2688 0.0000 0.1939 0.3969 0.3847 1.0396 0.2688 0.3847 0.9466 >> % El (1,2) se ha hecho distinto de cero; el (1,4) >> J=jacrot(A,1,4); A=J’*A*J A = 4.2172 0.2383 0.1165 0 0.2383 -0.1284 0.1939 0.2063 0.1165 0.1939 0.3969 0.3666 -0.0000 0.2063 0.3666 0.6161
38/95
>> % Ahora el (2,3) >> J=jacrot(A,2,3); A=J’*A*J A = 4.2172 0.1899 0.1852 0 0.1899 -0.1922 -0.0000 0.0814 0.1852 -0.0000 0.4607 0.4127 -0.0000 0.0814 0.4127 0.6161 >> % El (2,4) >> J=jacrot(A,2,4); A=J’*A*J A = 4.2172 0.1890 0.1852 0.0188 0.1890 -0.2003 -0.0409 -0.0000 0.1852 -0.0409 0.4607 0.4107 0.0188 0.0000 0.4107 0.6243 >> % El (3,4) >> J=jacrot(A,3,4); A=J’*A*J A = 4.2172 0.1890 0.1312 0.1320 0.1890 -0.2003 -0.0316 -0.0260 0.1312 -0.0316 0.1237 0.0000 0.1320 -0.0260 -0.0000 0.9612 >> % Después de un barrido el valor de off(A) se ha reducido >> off(A) ans = 0.2845 >> % Desde >> of1 of1 = 3.1996 >> % Otro barrido >> J=jacrot(A,1,2); A=J’*A*J A = 4.3202 -0.0000 0.1580 0.0270 -0.0000 -1.3369 -0.0502 0.0230 0.1580 -0.0502 -0.5132 0 0.0270 0.0230 -0.0000 -0.7836
39/95
>> J=jacrot(A,1,3); A=J’*A*J A = 4.3254 -0.0016 -0.0000 -0.0016 -1.3369 -0.0501 -0.0000 -0.0501 -0.5183 0.0270 0.0230 -0.0009 >> J=jacrot(A,1,4); A=J’*A*J A = 4.3255 -0.0015 -0.0000 -0.0015 -1.3369 -0.0501 -0.0000 -0.0501 -0.5183 0.0000 0.0230 -0.0009 >> J=jacrot(A,2,3); A=J’*A*J A = 4.3255 -0.0015 0.0001 -0.0015 -1.3400 0.0000 0.0001 -0.0000 -0.5153 0.0000 0.0229 -0.0023 >> J=jacrot(A,2,4); A=J’*A*J A = 4.3255 -0.0015 0.0001 -0.0015 -1.3409 0.0001 0.0001 0.0001 -0.5153 -0.0001 -0.0000 -0.0023 >> J=jacrot(A,3,4); A=J’*A*J A = 4.3255 -0.0015 0.0001 -0.0015 -1.3409 0.0001 0.0001 0.0001 -0.5152 -0.0001 0.0000 -0.0000 >> eig(A) % Casi ha convergido ans = 4.3255 -1.3409 -0.5152 -0.7829
0.0270 0.0230 -0.0009 -0.7836
-0.0000 0.0230 -0.0009 -0.7838
-0.0000 0.0229 -0.0023 -0.7838
-0.0001 0 -0.0023 -0.7828
-0.0001 0.0000 0 -0.7829 en dos barridos
40/95
Programa completo de Matlab para realizar el proceso. function [V D it]=Jacobi_val_12_2(A) % Cálculo por Jacobi de valores y vectores propios de una MATRIZ SIMÉTRICA A tol=sqrt(eps)*norm(A,’fro’); D=A; n=length(A); V=eye(n); [m1 p]=max(triu(abs(D),1)); % En (p,q) elemento mayor valor no en diagonal [~, q]=max(m1); % Posición fila máximo valor no cero en L(A) p=p(q); it=0; % Posición columna máximo anterior while off(D)>tol % Procesos iterativo; necesita rutina off J=Jac_Rot(D,p,q); % Se hacen cero Dpq y Dqp (p debe ser < q) D([p q],:)=J’*D([p q],:); D(:,[p q])=D(:,[p q])*J; V(:,[p q])=V(:,[p q])*J; [m1 p]=max(triu(abs(D),1)); [~, q]=max(m1); p=p(q); it=it+1; end [D I]=sort(diag(D)); V=V(:,I); end function a=off(A) % Calcula off de la matriz cuadrada A: raiz cuadrada de la suma % de los cuadrados de los coeficientes de A no en la % diagonal principal; también sqrt(sum(sum(triu(a.^2,1)))). n=length(A); a=0; for k=1:n-1 a=a+sum(A(k,k+1:n).^2); end a=sqrt(a); end
Para el ejemplo hecho a mano, se obtiene este resultado. >> A=[1 0 2;0 2 1;2 1 1]; >> [v d it]=Jacobi_val_12_2(A) v = 0.6611 -0.4973 0.5618 0.2261 0.8460 0.4828 -0.7154 -0.1922 0.6718 d = -1.1642 1.7729 3.3914 it = 8
>> [V D]=eig(A) V = -0.6611 -0.4973 -0.2261 0.8460 0.7154 -0.1922 D = -1.1642 0 0 1.7729 0 0
-0.5618 -0.4828 -0.6718
41/95
La aproximación que hicimos iterando un par de “barridos” era bastante buena.
0 0 3.3914
Para la matriz de la sesión interactiva anterior. >> [v d it]=Jacobi_val_12_2(A) v = -0.3673 0.3703 -0.6272 0.9032 0.1289 -0.0831 -0.1590 -0.7022 0.2688 -0.1549 0.5944 0.7263 d = -0.2134 0.1238 0.9569 4.2346 it = 16 >> A*v-v*diag(d) ans = 1.0e-010 * -0.0096 0.3780 0.0061 -0.0034 0.2318 0.0042 0.0183 0.4117 0.0067 -0.0155 0.2006 0.0032
0.5785 0.4009 0.6399 0.3086
0.2292 0.0812 -0.4444 0.3862
>> [v d]=eig(A) v = 0.3673 -0.3703 -0.9032 -0.1289 0.1590 0.7022 0.1549 -0.5944 d = -0.2134 0 0 0.1238 0 0 0 0 >> A*v-v*d ans = 1.0e-014 * -0.0222 -0.0444 -0.0056 0.0097 -0.0479 0.0180 -0.0069 -0.0083
0.6272 0.0831 -0.2688 -0.7263
0.5785 0.4009 0.6399 0.3086
0 0 0.9569 0
0 0 0 4.2346
-0.0333 -0.0305 -0.0444 -0.0111
-0.0444 -0.0444 -0.1332 0
Método de la iteración de la potencia
42/95
Su objetivo es calcular el valor propio dominante y su vector propio asociado.
Se basa en que al aplicarle a cualquier vector la transformación que representa una matriz A, ese vector tiende a orientarse hacia la dirección del vector propio dominante de A. Recordad.
Partiendo de un x 0 de módulo unidad, opera mediante una iteración de punto fijo de fórmula de recurrencia yk 1 : y k 1 D Ax k 1; donde x k D jjy k 1jj1
La magnitud jjy k 1jj1 converge al valor propio dominante, 1, y el vector x i lo hace al vector propio dominante, v1.
Su convergencia está ligada a j2=1j : a menor valor mejor convergencia.
Ejemplo Partiendo de x T0 D Œ0; 1 calculemos el valor propio dominante de 1;5 0;5 : 0;5 1;5 Utilicemos una pequeña sesión de Matlab como esta
k
Los resultados son estos.
0 1 2 3 4 5 6 7 8 9 10
x Tk 0,0 0,333 0,600 0,778 0,882 0,939 0,969 0,984 0,992 0,996 0,998
1,0 1,0 1,0 1,0 1,0 1,0 1,0 1,0 1,0 1,0 1,0
jjx k jj1 1,5000 1,6667 1,8000 1,8889 1,9412 1,9697 1,9846 1,9922 1,9961 1,9981
>> x=[0;1]; >> A=[1.5 0.5;0.5 1.5]; >> for i=1:10 v=A*x m=max(abs(v)) x=v/m end
>> x=[0;1]; A=[1.5 0.5;0.5 1.5]; for i=1:10 v=A*x, m=max(abs(v)), x=v/m end v = 0.5000 1.5000 m = 1.5000 x = 0.3333 1.0000 v = 1.0000 1.6667 m = 1.6667 x = 0.6000 1.0000 . . . v = 1.9942 1.9981 m = 1.9981 x = 0.9980 1.0000
43/95
Geometric Interpretation
of power iteration El comportamientoBehavior de las iteraciones es el dedepicted la figurageometrically:
1.0
0.5
0.0
x..0 x.1 x..2 x...3...x 4... ....... ... ...... ...... ...... ....... ....... ....... ... . . ... ... ................... v 1 . v 2 ..... . .... . . . . .. ... ... ... .... .... .... .... ... ... ..... .................... ... . . . . .. ... .... ... ... ... ............. ... .... ..... ......................... ... ... .... .... ................... ... ... .... .... ...................... ... . . . ..... ... ... .... ............................ ... ............ ... ......... .... −1.0
Initial vector
El punto inicial
0.0
−0.5
"
#
"
#
0.5
"
0 1 −1 x0 = =1 +1 1 1 1
1.0
#
in two eigenvectors contains equalcomponents 0 1 1 (shownxby D dashed arrows) D1 C1 0 1
1
1
Repeated multiplication by A causes compo-
es una combinaciónnent lineal de eigenvector los dos vectores propios vto1 larger y v2 . in first (corresponding eigenvalue, 2) to dominate, and hence sequence of vectors converges to that eigenvector
La multiplicación sucesiva por A causa que el coeficiente 35 en el primer vector propio sea el que domine, por lo que la sucesión converge a ese vector propio.
44/95
45/95
El método puede fallar por diversas razones: 1. Porque haya más de un valor propio con el mismo módulo, en cuyo caso las iteraciones puede que converjan a una combinación lineal de los correspondientes vectores propios.
Este caso es bastante habitual pues esos dos valores propios pueden ser un par complejo conjugado.
2. El vector de partida puede que tenga coeficiente cero en el valor propio dominante.
El error de redondeo en la práctica puede que introduzca un pequeño valor, por lo que este peligro rara vez ocurre.
3. Para una matriz real y un punto de partida también real, puede que nunca se converja a un vector complejo.
46/95
Un programa de Matlab para el método de la potencia. function [lambda,V,cnt]= potencia(A,x,tol,max1) % Método de la potencia para la obtención de un valor % y un vector propios de A if nargin> x=[0;1]; A=[1.5 0.5;0.5 1.5]; for i=1:6 v=A\x, m=max(abs(v)), x=v/m end v = -0.2500 0.7500 m = 0.7500 x = -0.3333 1.0000 . . v = -0.8333 0.9444 m = 0.9444 x = -0.8824 1.0000 v = -0.9118 0.9706 m = 0.9706 x = -0.9394 1.0000 v = -0.9545 0.9848 m = 0.9848 x = -0.9692 1.0000
Converge al valor propio 1 y vector propio Π1; 1T .
49/95
50/95
En Matlab la iteración sería así. function [lambda,V,cnt]= ItInversa_2(A,x,tol,max1) % Método de la potencia inversa para la obtención % de un valor y vector propio de A if nargin> [lambda,V,cnt]=ItInversa(A,[0;1]) lambda = 1.0000 V = -0.7071 0.7071 cnt = 55
51/95
Mejora del método: desplazamiento
El desplazamiento es particularmente útil en el caso de la iteración inversa.
Si es el valor propio de A más próximo a , el valor propio de menor magnitud de A I será .
Con una elección apropiada de se puede calcular cualquier valor propio de A que esté próximo a esa . Si el desplazamiento está próximo a un valor propio, la convergencia es muy rápida.
Iteración mediante cociente de Rayleigh
John William Strutt, Lord Rayleigh, Reino Unido, 1842-1919.
T Definición El cociente de Rayleigh de un vector x 2 Rn es el escalar r.x/ D xx TAx x .
Dada una matriz A 2 Rnn y uno de sus vectores propios, x, la mejor estimación del correspondiente valor propio se puede considerar un problema de mínimos cuadrados n 1 como este: x Ax.
De sus ecuaciones normales, x Tx D x TAx, la solución que se obtiene es x TAx D T : x x
52/95
53/95
El cociente de Rayleigh puede utilizarse como desplazamiento para acelerar la convergencia de los métodos iterativos que hemos visto. El algoritmo de la iteración inversa quedaría así. Dado un punto de partida x 0 for k D 1, 2, : : :
x Tk 1 Ax k 1 k D T xk 1xk 1
Resolver .A k I/y k D x k yk xk D jjy k jj2 end
1
H El versión del cociente de Rayleigh para vectores complejos es xx HAx x .
Ejemplo (continuación)
El método de la iteración de potencia con cociente de Rayleigh aplicado paso a paso al problema que estudiamos: x Tk
k 0 1 2 3 4 5 6
54/95
0,000 0,333 0,600 0,778 0,882 0,939 0,969
1,0 1,0 1,0 1,0 1,0 1,0 1,0
jjy k jj1
x Tk Ax k =x Tk x k
0,150 1,667 1,800 1,889 1,941 1,970
1,500 1,800 1,941 1,985 1,996 1,999
k Partiendo de un punto aleatorio:
0 1 2
x Tk 0,807 0,924 1,000
k 0,397 1,000 1,000
1,896 1,998 2,000
En Matlab: function [lambda,V,cnt]= ItInvRayleigh_2(A,x,max1) % Método de la iteración inversa con cociente de Rayleigh if nargin1 x=x/norm(x); end tol=n*norm(A,1)*eps; sigma=x’*A*x; while err>tol y=(A-sigma*eye(n))\x; x=y/norm(y); sigma=x’*A*x; err=norm(A*x-sigma*x,1); cnt=cnt+1; end lambda=sigma; V=x; end
>> A=[1.5 0.5;0.5 1.5]; >> [lambda,V,cnt]=ItInvRayleigh_2(A,[0.5;0.55]) lambda = 2.0000 V = 0.7071 0.7071 cnt = 3
55/95
Deflación
Idea: Calculados un valor y un vector propios, obtener otros por deflación: Como cuando una vez conocida una de las raíces, x1, de un polinomio de grado n éste se divide por x x1 obteniéndose otro de grado n 1.
Si x es el vector propio asociado al valor propio dominante de A, 1, y H es una matriz de Householder, tal que H x D ˛e 1, la transformación de semejanza que define consigue transformar A así T b A 1 D H AH 1 D 1 ; 0 A2 donde A 2 es una matriz de orden n los restantes de A.
1 cuyos valores propios son 2; : : : ; n,
Luego, trabajando con A 2 y calculado 2, si y 2 es un vector propio asociado, el vector bT y 2 1 ˛ x2 D H ; ; donde ˛ D y2 2 1
56/95
es el vector propio correspondiente a 2 en la matriz original A, supuesto 1 ¤ 2 . Este método funciona bien para calcular algunos valores y vectores propios.
En Matlab:
function [l2 v2 B] = defl_1(A, v1) % B resulta de A por deflación de vec. propio v1. % Calcula el valor propio l2 y vec. propio v2. n = length(A); v1 = Housv(v1); C = Houspre(v1,A); B = zeros(n,n); for i=1:n B(:,i)=Housmvp(v1,C(i,:)); end + l1 = B(1,1); b = B(1,2:n); B = B(2:n,2:n); [l2 y] = ItInvRayleigh_2(B); if l1~=l2 a = b*y/(l2-l1); v2 = Housmvp(v1,[a;y]); else v2 = v1; end end
function u = Housv(x) % Transformación de Householder del vector x. m = max(abs(x)); u = x/m; if u(1) == 0, su = 1; else su = sign(u(1)); end u(1) = u(1)+su*norm(u); u = u/norm(u); u = u(:); end function P = Houspre(u, A) % Producto P = H*A, donde H is una transformación de Householder % definida por el vector u. n = length(A); v = u/norm(u); v = v(:); P = zeros(n,n); for j=1:n aj = A(:,j); P(:,j)=aj-2*(v’*aj)*v; end end function p = Housmvp(u, x) % Producto p = H*x, donde H es la transformación de Householder % definida por u. u = u(:); x = x(:); v = u/norm(u); p = x - 2*(v’*x)*v; end
57/95
Utilicémoslo: >> A=ones(5)+diag(ones(5,1)) A = 2 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 2 >> [lambda,V]= ItInvRayleigh_2(A) lambda = 6 V = 0.44721359549990 0.44721359550002 0.44721359549982 0.44721359550001 0.44721359550003 >> [l2,v2]=defl_1(A,V) l2 = 1 v2 = -0.89442719099992 0.22360679774998 0.22360679774998 0.22360679774998 0.22360679774998
% Generamos una matriz dada
% Cal. valor y vector propio dominantes
% Defl para obtener otro par val. A
>> [norm(A*V-lambda*V);norm(A*v2-l2*v2)] % Comprobar resultados ans = 1.0e-015 * 0.4441 0.6497 >> format short >> [l,v]=eig(A) % Comprobar con rutina Matlab eig l = 0.8333 -0.1667 -0.1667 0.2236 0.4472 -0.1667 0.8333 -0.1667 0.2236 0.4472 -0.1667 -0.1667 0.8333 0.2236 0.4472 -0.5000 -0.5000 -0.5000 0.2236 0.4472 0 0 0 -0.8944 0.4472 v = 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 6
58/95
Iteración simultánea o de subespacio
Para obtener varios valores/vectores propios, apliquemos la iteración de la potencia a los p vectores linealmente independientes que forman X 0; es decir: X 0 D matriz n p de rango p for k D 1, 2, : : : X k D AX k 1 end
La Im A x 1 A x 2 : : : A xp converge al subespacio invariante de los vectores propios de A correspondientes a sus p mayores valores propios: j1j > jp j > jpC1j jpC2j:
La velocidad de convergencia es O.pC1=p /.
k
k
k
59/95
Para optimizar prestaciones y conseguir una buena ortogonalidad se usa la factorización QR de X k , ortonormalizando así las columnas de X k . Se obtendría una Q como base de Im.X k /. X 0 D matriz n p de rango p for k D 1, 2, : : : Qk R k D X k 1 X k D AQk end
Las iteraciones convergen a una matriz triangular, si sólo hay valores propios reales, o a una triangular en bloques 22 en la diagonal principal cuando haya valores propios complejos conjugados.
60/95
Sesión en Matlab para matriz simétrica
>> A=randn(5); A=A*A’ A = 2.8486 -2.5688 0.1741 3.5081 -0.1959 -2.5688 2.8789 -0.3177 -4.0552 -2.2033 0.1741 -0.3177 2.3957 1.4227 -2.3486 3.5081 -4.0552 1.4227 6.2012 2.7463 -0.1959 -2.2033 -2.3486 2.7463 19.4886 >> A0=A; >> % Hagamos ahora 10 iteraciones de la iteración ortogonal >> for i=1:10, [Q R]=qr(A); A=R*Q; end >> A A = 20.8010 -0.0997 0.0000 0.0000 0.0000 -0.0997 10.6605 0.0000 0.0000 0.0000 0.0000 0.0000 2.2265 0.0000 0.0000 0.0000 -0.0000 0.0000 0.1085 0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0165
>> % Otras 5 iteraciones más >> for i=1:5, [Q R]=qr(A); A=R*Q; end >> A A = 20.8020 -0.0035 0.0000 0.0000 0.0000 -0.0035 10.6595 0.0000 0.0000 0.0000 0.0000 0.0000 2.2265 0.0000 0.0000 0.0000 -0.0000 0.0000 0.1085 0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0165 >> % Mejora; hagamos 10 más >> for i=1:10, [Q R]=qr(A); A=R*Q; end >> A A = 20.8020 -0.0000 0.0000 0.0000 0.0000 -0.0000 10.6595 -0.0000 0.0000 0.0000 0.0000 0.0000 2.2265 0.0000 0.0000 0.0000 -0.0000 0.0000 0.1085 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0165 >> % Comparemos todos con los de la matriz original >> eig(A0) ans = 0.0165 0.1085 2.2265 10.6595 20.8020
IterSimul obtiene los r valores propios máximos de A simétrica.
function [i r1]=IterSimul(A,r) % Máximos valores propios de A (simétrica); itera. simultánea % ortogonal n=length(A); if nargin> A=rand(20); A=A*A’; >> [i r]=IterSimul(A,6) Valores calculados con eig() 1.0e+002 * 1.006294230413485 0.052804894629075 0.039516583781399 0.032191801462925 0.031157323965477 0.026806664806072 0.022597606348368 0.019827755676060 0.013685589405813 0.012217263560720 0.009791335827920 0.005048187809920 0.004282145342196 0.003460064312945 0.002436587767534 0.001713272489506 0.001168137691697 0.000605614254568 0.000238729619383 0.000012007660353 i = 44 r = 1.0e+002 * 1.006294230413484 0.052804894629054 0.039516583695056 0.032053644988988 0.031291616159557 0.026806665421155
61/95
Iteración QR
62/95
G. GOLUB AND F. UHLIG
family until his death. At first, Kantorovich’s group developed analytic computational tools in Prorab for algebraic and trigonometric polynomials, for integer arithmetic and series, etc. When Vera joined, her task was to select and classify matrix operations that are useful in numerical linear algebra. Linear algebra subroutines were included in Prorab much later. This experience brought Vera close to numerical algebra and computation. In 1972 she obtained her secondary doctorate (Habilitation). More extended biographies of Vera’s life, achievements and long and productive career are available in Golub et al. (1990) and Konkova et al. (2003), for example.
Propuesta por John G.F. Francis y Vera Kublanovskaya. Downloaded from http://imajna.oxfordjournals.org at Università di Bologna - Sistema Bibliotecario d'Ateneo on March 30, 2010
474
John G.F. Francis, Reino Unido 1934 y Vera Kublanovskaya, Rusia 1920-2012. Vera Kublanovskaya in August 2008
In the Russian literature the QR algorithm was initially called the method of one-sided rotations. Vera Kublanovskaya started to develop her version of the QR algorithm in 1958 after reading Rutishauser’s LR algorithm paper (Rutishauser, 1958). She represented a matrix A = L · Q as the product of a lower triangular matrix L and an orthogonal matrix Q that she suggested to compute as the product of elementary rotations or reflections. Her eigenvalue algorithm factors A = A1 = L 1 · Q 1 , reverse order multiplies A2 = Q 1 · L 1 and then factors A2 again as A2 = L 2 · Q 2 , etc. In 1959 she performed numerical experiments with her LQ decomposition-based algorithm on an electromechanic ‘Mercedes’ calculator in low dimensions by hand since linear algebra computer codes had not been developed yet in Prorab. Vera’s hand computations indicated that her explicit factorization—reverse-order multiply algorithm— is convergent. However, D. K. Faddeev feared that the obtained results may only have been accidental. Vera’s QR-like algorithm was briefly mentioned in 1960 in the first edition of a monograph
Generalizando la iteración simultánea (ortogonal) se pueden calcular todos los valores propios de A, así como los correspondientes vectores propios. A 0 D AI U 0 D I for k D 1, 2, : : : Qk 1 R k 1 D A k 1 , obtener factorización QR A k D R k 1 Qk 1 ; U k D U k 1 Qk 1 end
Las A k convergerán a una matriz triangular, si todos los valores propios son reales, o con bloques 22 en la diagonal, si hay valores propios complejos. U k convergerá a los vectores propios.
63/95
Las matrices A k son ortogonalmente semejantes a las anteriores del proceso y a A: A 1 D R 0Q0 D QT0 A 0Q0 pues R 0 D QT0 A 0 A 2 D R 1Q1 D QT1 A 1Q1 D QT1 QT0 A 0Q0 Q1 D .Q0Q1/T A 0 .Q0Q1/
Es decir, todas las A k tendrán los valores propios de A D A 0.
Las columnas de Qk forman una base ortonormal del subespacio Im.A k /.
Si A es simétrica, las iteraciones conservan la simetría y las A k convergerán a una matriz triangular y simétrica: diagonal.
Trabajemos un poco “a mano” con el algoritmo.
64/95
D = diag([4 3 2 1]); rand(’seed’,0); format short e S=rand(4); S = (S - .5)*2; % Se perturba un poco A = S*D/S % A_0 = A = S*D*S^{-1} for i=1:6 [Q R] = qr(A); A = R*Q % Iteración QR end A =
A
A
A
A
A
4.1475e+00 -4.9214e-03 6.2706e-01 1.5691e-01 = 3.9304e+00 4.0097e-02 3.5144e-01 -2.7163e-02 = 3.9287e+00 -1.6308e-02 2.1041e-01 5.9415e-03 = 3.9514e+00 -5.7605e-02 1.1812e-01 -1.4093e-03 = 3.9696e+00 -7.1970e-02 6.3115e-02 3.4472e-04 = 3.9817e+00 -6.9838e-02 3.2744e-02 -8.5395e-05
-7.7246e-01 3.2731e+00 -1.9699e-01 -4.1310e-01
-7.3819e-01 -1.8990e-01 2.2190e+00 -6.3508e-01
-4.1172e+00 1.9440e+00 -6.0845e-01 3.6039e-01
-2.7682e-01 3.0595e+00 -3.8535e-02 8.4395e-02
1.5214e-01 -6.9040e-01 2.3156e+00 1.9457e-01
4.4847e+00 -2.1283e+00 -4.8346e-01 6.9447e-01
-1.6932e-01 3.0078e+00 -3.4259e-02 -2.3299e-02
4.1119e-01 -8.8862e-01 2.2044e+00 -8.1008e-02
-4.4113e+00 2.1571e+00 1.0227e+00 8.5905e-01
-1.3627e-01 2.9981e+00 -2.7149e-02 7.1443e-03
4.9740e-01 -9.6432e-01 2.1181e+00 3.8096e-02
4.3484e+00 -2.1072e+00 -1.3251e+00 9.3241e-01
-1.1610e-01 2.9977e+00 -1.8643e-02 -2.2854e-03
5.3907e-01 -9.8790e-01 2.0658e+00 -1.8738e-02
-4.3288e+00 2.0192e+00 1.4919e+00 9.6688e-01
-9.8200e-02 2.9984e+00 -1.2095e-02 7.4567e-04
5.6721e-01 -9.9032e-01 2.0362e+00 9.3684e-03
4.3356e+00 -1.9242e+00 -1.5816e+00 9.8360e-01
¿Alguna conclusión en términos de convergencia; trabajo?
65/95
Iteración QR con desplazamiento
La velocidad de convergencia se puede aumentar con desplazamientos: Qk 1 R k
1
D Ak
1
k 1I
A k D R k 1 Qk El desplazamiento k
1
1
C k 1I:
debe ser una buena aproximación de un valor propio.
Desplazamientos más usados:
Rayleigh. El valor de ann de A k 1: buena aproximación a priori de n. Wilkinson. Del bloque diagonal inferior 22 de A k 1 (si no es diagonal), el valor propio de esa submatriz que esté más cerca de ann.
66/95
Este script de Matlab es el método de la iteración QR con desplazamiento de Rayleigh para matrices con valores propios reales: simétricas, hermíticas, etc. function [lambda itt]= IteracionQR(A,tol) % Método de la iteración QR con desplazamiento % para matrices con valores propios reales if nargin tol sigma=A(k,k); [Q,R]=qr(A(1:k,1:k)-sigma*eye(k,k)); A(1:k,1:k)=R*Q+sigma*eye(k,k); it=it+1; end end lambda=sort(diag(A)); if nargout==2, itt=it; end end
67/95
Ejemplos
>> ab=ones(5)+eye(5) ab = 2 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 2 >> [lambda it]=IteracionQR(ab) lambda = 1.0000 1.0000 1.0000 1.0000 6.0000 it = 4
>> a=rand(20); >> [lambda it]=IteracionQR(a’*a) lambda = 0.0010 0.0175 0.1107 0.1496 0.2232 0.4124 0.7861 0.8132 0.9430 1.2765 1.4153 1.7930 1.9123 2.2906 2.6955 2.9732 3.7114 4.9037 5.3983 100.1106 it = 48
>> eig(a’*a) ans = 0.0010 0.0175 0.1107 0.1496 0.2232 0.4124 0.7861 0.8132 0.9430 1.2765 1.4153 1.7930 1.9123 2.2906 2.6955 2.9732 3.7114 4.9037 5.3983 100.1106
Sesión en Matlab para matriz sin estructura especial
>> A=randn(5) A = 0.1303 1.0169 -0.2750 0.7688 -0.3857 0.5583 0.2378 -0.2941 -0.5498 -0.1324 0.3125 0.8049 0.2869 -1.1517 -0.0024 -1.7196 -0.3503 0.4211 0.9847 -0.6632 >> A=hess(A) % VEREMOS por qué A = 0.1303 -0.1471 1.0342 -0.5765 0.8100 0.6165 -0.1303 0.5381 0 -1.0393 0.3560 0.3942 0 0 1.3485 -1.9025 0 0 0 0.4372 >> A0=A;
0.2815 1.2163 -0.2401 -0.4955 -1.3721
>> % Otras 10 iteraciones más >> for i=1:10, [Q R]=qr(A); A=R*Q; end >> A A = -1.6585 -0.7415 0.3436 -0.8538 0.1202 -2.0481 -0.4571 -0.6504 0 -0.0322 0.6956 -0.7333 0 0 1.0035 1.0146 0 0 0 0.0000
-0.5971 0.9983 0.6386 -0.0670 -1.2908
>> % Mejora; hagamos 50 más >> for i=1:50, [Q R]=qr(A); A=R*Q; end >> A A = -1.7423 -0.7746 -0.7995 -0.5685 0.0805 -1.9665 -0.6724 0.2789 0 -0.0000 1.0634 -0.8844 0 0 0.8436 0.6489 0 0 0 0.0000
>> % Hagamos ahora 10 iteraciones de la iteración QR >> for i=1:10, [Q R]=qr(A); A=R*Q; end >> A A = -1.9640 -1.2063 -0.1339 0.3196 -0.3104 0.1231 0.3343 -0.7129 1.0207 -0.4763 0 -1.5358 -0.9798 -0.1630 -0.8736 0 0 1.1031 0.6130 0.2523 0 0 0 0.0000 -0.0941
0.1063 -1.0568 -0.1472 -0.0113 -0.0941
>> % Intentemos unas 50 últimas >> for i=1:50, [Q R]=qr(A); A=R*Q; end >> A A = -1.8704 -0.7919 -0.1944 1.0214 0.0632 -1.8384 0.5215 0.3733 0 -0.0000 0.7034 -0.7224 0 0 1.0056 1.0089 0 0 0 0.0000
68/95
-0.2176 -1.0411 0.1361 0.0134 -0.0941
>> % Hay dos bloques 2x2 en la diagonal >> % El último parece claro que es -0.0941
-0.0297 -1.0632 -0.0561 0.1247 -0.0941
>> % Mejora, aunque no suficientemente; otras 50 >> for i=1:50, [Q R]=qr(A); A=R*Q; end >> A A = -1.8095 -0.7895 -0.9613 0.3250 -0.1298 0.0657 -1.8993 -0.1186 0.6697 -1.0556 0 -0.0000 0.7945 -1.0629 0.0612 0 0 0.6651 0.9179 0.1223 0 0 0 0.0000 -0.0941
>> % Calculemos los valores propios de esos bloques >> eig(A(1:2,1:2)) ans = -1.8544 + 0.2232i -1.8544 - 0.2232i >> eig(A(3:4,3:4)) ans = 0.8562 + 0.8385i 0.8562 - 0.8385i >> % Comparemos todos con los de la matriz original >> eig(A0) ans = 0.8562 + 0.8385i 0.8562 - 0.8385i -0.0941 -1.8544 + 0.2232i -1.8544 - 0.2232i
Iteración QR con doble desplazamiento
69/95
Si hay valores propios complejos, se incorpora el doble desplazamiento para tratar los bloques 22 irreducibles de la diagonal principal.
Se usan los dos valores propios, y N , de la submatriz 22 del bloque inferior, Qk 1 R k
1
D Ak
1
k 1 I
A k D R k 1 Qk
Qk R k D A k
1
C k 1 I
N k 1 I
A kC1 D R k Qk C N k 1 I
haciendo así un paso doble con desplazamiento en dos etapas.
Este doble paso se puede simplificar si M D .A k
1
k 1 I/ .A k
D .Qk 1 Qk / .R k R k 1 / D A 2k
1
sA k
1
siendo s la traza del bloque diagonal inferior 22 de A k
C tI; 1
1
N k 1 I/
y t su determinante.
Si se factorizase la matriz así, M D QR, mediante una sencilla manipulación algebraica resultaría que A kC1 D QT M Q, por lo que teóricamente, en lo que se denomina doble desplazamiento implícito, se evitaría un paso.
Calcular M requiere O.n3/ operaciones y su factorización ulterior QR es bastante inestable.
Con los años, los diferentes enfoques de esta variante han convergido en una estrategia de doble desplazamiento implícito tipo Francis, formulada en 1961 por John G.F. Francis.
70/95
John G.F. Francis, Reino Unido 1934
Dado el coste O.n3/, y su complejidad, hay que preparar el problema antes de usar este método.
Algoritmo QR: Transformaciones preliminares Definición Una matriz de Hessenberg es una matriz triangular excepto por una subdiagonal adyacente a la diagonal principal.
@ @ @
0
@ @ @ @
Cualquier matriz se puede reducir a la forma de Hessenberg mediante transformaciones de Householder o Givens.
Si la matriz original es simétrica, al reducirla a la forma de Hessenberg se obtendrá una tridiagonal.
71/95
72/95
Si se parte de una matriz con la forma de Hessenberg, el algoritmo QR más actual necesita un número de operaciones O.n2/, implementándose en dos fases: Matriz
1ª fase
Simétrica a tridiagonal General
2ª fase a diagonal
a Hessenberg a triangular
El trabajo necesario total de todo el proceso QR, preliminar más iterativo, es el que sigue. Matriz simétrica
Matriz general
4 3 n para valores propios sólo 10n3 para valores propios sólo 3 9n3 valores y vectores propios 25n3 valores y vectores propios
73/95
Un conjunto de rutinas sencillas para reducir cualquier matriz, mediante transformaciones de Householder, a la forma de Hessenberg seria este. function [A V] = Hessred(A) % Reducción de A a Hessenberg con Householder % En V vectores de transf. sucesivas de Householder [m,n] =size(A); if A == triu(A,-1), V = eye(m); return, end V = []; for k=1:m-2 x = A(k+1:m,k); v = Housv(x); A(k+1:m,k:m)=A(k+1:m,k:m)-2*v*(v’*A(k+1:m,k:m)); A(1:m,k+1:m)=A(1:m,k+1:m)-2*(A(1:m,k+1:m)*v)*v’; v = [zeros(k,1);v]; V = [V v]; end end function u = Housv(x) % Transf. Householder del x; vector en u. m = max(abs(x));u = x/m; if u(1) == 0, su = 1; else su = sign(u(1)); end u(1) = u(1)+su*norm(u); u = u/norm(u); u = u(:); end
74/95
Consideraciones finales en torno al algoritmo o método QR:
Para matrices muy grandes, el algoritmo es costosísimo. Si sólo se necesitan unos pocos valores y vectores propios, especialmente para n grandes, el método no saca partido de ello. Requiere mucha memoria de ordenador, aunque la matriz sea grande y dispersa. Las transformaciones por semejanza introducen muchos elementos no nulos por lo que se destruye una posible estructura de dispersidad.
75/95
Este es un programa del algoritmo QR para matrices simétricas con reducción inicial a forma Hessenberg y desplazamiento sencillo de Wilkinson. Se prueba con una matriz 20 20, simétrica, generada aleatoriamente. function [D itt]=It_QR_3(A) % Iteración QR con despla. de Wilkinson para una matriz simétrica n=length(A); tol=off(A)*1.e-8; k=n; D=zeros(n,1); it=0; A=Hessred(A); % Reducir a Hessenberg while k>1 while abs(A(k,k-1))>tol %Calcular desplazamiento Wilkinson s=eig2x2(A(k-1:k,k-1:k)); [i j]=min(abs(A(k,k)*[1 1]’-s)); % Mejor desp. Wilkinson [Q R]=qr(A(1:k,1:k)-s(j)*eye(k)); A(1:k,1:k)=R*Q+s(j)*eye(k); it=it+1; end k=k-1; end D=sort(diag(A),’ascend’); if nargout==2, itt=it; end end function [L]=eig2x2(a) tra=a(2,2)+a(1,1); sqtd=sqrt(tra*tra-4*(a(2,2)*a(1,1)-a(1,2)*a(2,1))); L(1)=(tra+sqtd)/2; L(2)=(tra-sqtd)/2; L=L(:); end
>> A=rand(20); >> A=(A+A’)/2; >> [D it]=It_QR_3(A) D = 10.5361 1.5883 1.5200 1.1285 1.0568 0.9498 0.6448 0.5156 0.2413 0.1139 -0.0201 -0.1461 -0.4263 -0.6220 -0.6692 -0.8182 -0.8895 -1.0250 -1.2618 -1.2714 it = 36
>> eig(A) ans = -1.2714 -1.2618 -1.0250 -0.8895 -0.8182 -0.6692 -0.6220 -0.4263 -0.1461 -0.0201 0.1139 0.2413 0.5156 0.6448 0.9498 1.0568 1.1285 1.5200 1.5883 10.5361
function [T it itd t] = qrstep_Francis_Bindel(H) % Computes eigenvalues of A using a implicit double-shift QR method % as Cornell’s Bindel notes (CS 620). % % Output: TT, vector containing eigenvalues; it, simple iterations; % itd, double iterations. global itr itc n = length(H); T=zeros(n,1); tol = norm(H,’fro’)*1e-8; H=hess(H); % Reduce A to Hessenberg form itr=0; itc=0; t=cputime; while n > 2 if abs(H(n,n-1)) < tol*(abs(H(n-1,n-1))+abs(H(n,n))) H(n,n-1) = 0; T(n)=H(n,n); % Deflate 1x1 block n = n-1; elseif abs(H(n-1,n-2)) < tol*(abs(H(n-2,n-2))+abs(H(n-1,n-1))) H(n-1,n-2) = 0; T(n-1:n)=eig2x2(H(n-1:n,n-1:n)); % Deflate 2x2 block n = n-2; else H(1:n,1:n) = qrstep(H(1:n,1:n)); % Main body end end T(1:2)=eig2x2(H(1:2,1:2)); T=sort(T,’descend’); t=cputime-t; it=itr; itd=itc; end function [NH] = qrstep(H) % Compute double-shift poly and initial column of H^2 + b*H + c*I [b c] = qrpoly(H); C1 = H(1:3,1:2)*H(1:2,1); C1(1:2) = C1(1:2) + b*H(1:2,1); C1(1) = C1(1) + c; % Apply a similarity associated with the first step of QR on C v = Housv(C1); H(1:3,:) = H(1:3,:)-2*v*(v’*H(1:3,:)); H(:,1:3) = H(:,1:3)-(H(:,1:3)*(2*v))*v’; % Do "bulge chasing" to return to Hessenberg form n = length(H); for j=1:n-2 k = min(j+3,n); v = Housv(H(j+1:k,j)); % Find W=I-2vv’ to put zeros below H(j+1,j), H=WHW’ H(j+1:k,:) = H(j+1:k,:)-2*v*(v’*H(j+1:k,:)); H(:,j+1:k) = H(:,j+1:k)-(H(:,j+1:k)*(2*v))*v’; H(k,j) = 0; end NH=triu(H,-1); end
76/95
function [b c] = qrpoly(H) % Compute b and c of z^2+b*z+c = (z-sigma)(z-conj(sigma)) % where sigma is the Wilkinson double shift: eigenvalue of % the 2x2 trailing submatrix closest to H(n,n) global itr itc HH = H(end-1:end,end-1:end); traHH = HH(1,1)+HH(2,2); detHH = HH(1,1)*HH(2,2)-HH(1,2)*HH(2,1); if traHH^2 > 4*detHH % Real eigenvalues % Use as double shift the eigenvalue closer to H(n,n) lambdaHH(1) = (traHH + sqrt(traHH^2-4*detHH))/2; lambdaHH(2) = (traHH - sqrt(traHH^2-4*detHH))/2; if abs(lambdaHH(1)-H(end,end))> [T it itd t]=qrstep_Francis_Bindel(A) T = -6.0770 5.2804 5.0190 3.2085 + 3.6195i 3.2085 - 3.6195i 0.0031 + 4.4372i 0.0031 - 4.4372i -4.2104 + 0.7456i -4.2104 - 0.7456i -1.3585 + 3.9497i -1.3585 - 3.9497i 1.2975 + 3.8800i 1.2975 - 3.8800i -2.4416 + 2.7407i -2.4416 - 2.7407i 2.8787 + 0.7120i 2.8787 - 0.7120i -1.0840 + 1.8066i -1.0840 - 1.8066i -1.2594 + 0.9430i -1.2594 - 0.9430i -0.3809 + 0.9062i -0.3809 - 0.9062i 0.6323 -0.1911 it = 7 itd = 27 t = 0.0624
>> sort(eig(A),’descend’) ans = -6.0770 5.2804 5.0190 3.2085 + 3.6195i 3.2085 - 3.6195i 0.0031 + 4.4372i 0.0031 - 4.4372i -4.2104 + 0.7456i -4.2104 - 0.7456i -1.3585 + 3.9497i -1.3585 - 3.9497i 1.2975 + 3.8800i 1.2975 - 3.8800i -2.4416 + 2.7407i -2.4416 - 2.7407i 2.8787 + 0.7120i 2.8787 - 0.7120i -1.0840 + 1.8066i -1.0840 - 1.8066i -1.2594 + 0.9430i -1.2594 - 0.9430i -0.3809 + 0.9062i -0.3809 - 0.9062i 0.6323 -0.1911
>> A=rand(1000); >> A = (A - .5)*2; >> [T it itd t]=qrstep_Francis_Bindel(A) T = 7.123846740573868 +17.026981501328613i 7.123846740573868 -17.026981501328613i 16.624570346272108 + 7.937272245254786i 16.624570346272108 - 7.937272245254786i -16.089705600685924 + 8.930438288025009i -16.089705600685924 - 8.930438288025009i -10.166595375856957 +15.258792634828080i -10.166595375856957 -15.258792634828080i 17.167325994176007 + 6.377944580510263i 17.167325994176007 - 6.377944580510263i 18.249349262178139 + 0.000000000000000i 7.805032835495162 +16.405785237414399i 7.805032835495162 -16.405785237414399i -8.822496260367785 +15.869321050483862i -8.822496260367785 -15.869321050483862i 16.001083815170652 + 8.555659101012008i 16.001083815170652 - 8.555659101012008i -3.803035356678313 +17.729394693350987i -3.803035356678313 -17.729394693350987i 4.837154691757405 +17.460853314474601i 4.837154691757405 -17.460853314474601i -17.976118771479754 + 2.032262837884847i -17.976118771479754 - 2.032262837884847i . . . -1.748103565292426 - 1.443294519531504i -2.229777171574543 + 0.000000000000000i 2.124127007088685 + 0.000000000000000i -2.000980926645378 + 0.000000000000000i -0.267964043130730 + 1.702030627215852i -0.267964043130730 - 1.702030627215852i 1.251743319582046 + 0.303853777326017i 1.251743319582046 - 0.303853777326017i -0.989726686530990 + 0.000000000000000i 0.222221915982704 + 0.687550759680182i 0.222221915982704 - 0.687550759680182i it = 93 itd = 1116 t = 41.563322418069106 >> tic, [V S U]=eig(A); toc Elapsed time is 1.323848 seconds.
78/95
Subespacios de Krylov
Alekxey Krylov, Rusia 1863-1945.
Recordamos que, para A 2 Cnn y 0 ¤ b 2 Cn1, al subespacio Kj D Genfb; Ab : : : A j 1bg se le denomina subespacio de Krylov.
Los métodos que trabajan en los sucesivos subespacios de Krylov que se crean sólo llevan a cabo multiplicaciones de matrices por vectores y van reduciendo la matriz a forma Hessenberg.
79/95
Si A tiene n valores propios distintos, 1; : : : ; n, con vectores propios asociados x 1; : : : ; x n ortonormales (base ortonormal de Rn), cualquier vector de Rn se puede escribir como combinación lineal de esos vectores propios; en particular, uno b D c1x 1 C c2x 2 C C cnx n. La matriz de Krylov, K j D b Ab : : : A j 1b , de dimensión n j , se puede escribir 2 3 j 1 1 1 61 6 2 K j D Œc1 x 1 c2 x 2 cn x n nn 6 :: :: 4: : 1 n
1 j 1 2 7 7 : : : : ::: 7 5 j 1 n nj
Cuando j D n, la matriz C n D K n AK n es tipo Hessenberg: 5
1
@ @
0
@ @ @
5
De hecho C n es la matrix compañera de A.
80/95
La matriz de Krylov K n suele estar mal condicionada. Para obtener una buena base de Im.K n/ se utiliza la factorización QnR n D K n de tal forma que H 1 C n D K n 1 AK n D R 1 QH n AQn R n ! Qn AQn D R n C n R n H (matriz de Hessenberg):
Si se igualan las columnas k de la expresión anterior escrita AQn D QnH se tiene que Aq k D h1k q 1 C C hkk q k C hkC1;k q kC1 expresión que relaciona el vector q kC1 con los anteriores q 1; : : : ; q k .
Si se premultiplica por qjH , teniendo en cuenta la ortonormalidad, hj k D qjH Aq k ;
j D 1; : : : ; k:
Estas expresiones dan lugar a la ya estudiada iteración de Arnoldi. Obtiene una matriz unitaria Qn y una Hessenberg H n, columna a columna, haciendo sólo productos de A por vectores y productos interiores de estos.
81/95
Ortogonalización de A por Iteración de Arnoldi
Dado x 0 cualquiera
q 1 D x 0 = kx 0 k2 for k D 1; 2; : : : uk D Aq k for j D 1 to k hj k D qjH uk uk D uk hj k qj
Según progresa el algoritmo, en la iteración k , si Qk D Œq 1 ; : : : ; q k , la matriz H k D QH k AQk es Hessenberg y sus valores propios, denominados valores de Ritz, por
Walther Ritz, Suiza, 1878-1909,
end
hkC1;k D kuk k2 if hkC1;k D 0 stop q kC1 D uk = hkC1;k
son muy buenas aproximaciones de los valores propios de A . Si y es un vector propio de H k , Qk y es un aproximación de un vector propio de A .
end
Si se requieren con precisión los valores propios de H k se pueden calcular por otro método –por ejemplo la iteración QR–, siendo una tarea menor si k n.
El coste en cálculos de la iteración de Arnoldi es elevado: cada nuevo q k se debe ortogonalizar respecto a todos los vectores columna de Qk , por lo que se reinicializa periódicamente desde un nuevo vector escogido adecuadamente.
Si A es simétrica o hermítica se utiliza la iteración de Lanczos que consigue una matriz tridiagonal.
Ortogonalización de A por Iteración de Lanczos
q 0 D 0; ˇ0 D 0 y x 0 (cualquiera) q 1 D x 0 = kx 0 k2 for k D 1; 2; : : : uk D Aq k ˛k D q H k uk uk D uk ˇk 1 q k 1 ˛k q k ˇk D kuk k2 if ˇk D 0 stop q kC1 D uk =ˇk
Cornelius Lanczos (Kornel Lowy), Hungría, 1893-1974.
end
Si ˇk D 0 los valores de Ritz son los valores propios exactos de A.
82/95
En Matlab:
83/95
function [U H flag] = arnoldi_W(A,ns,x) % Realiza ns iteraciones de Arnoldi con reortogonalización en A partiendo de x % % Si flag == 0, las columnas de U son ortonormales, H(ns+1,ns) matriz Hessenberg y % A*U(:,1:ns) = U(:,1:ns+1)*H. % % Si flag == j > 0, se ha obtenido subespacio invariante en j iter. y % A*U(:,1:j) = U(:,1:j)*H(1:j,1:j). % flag = 0; H = zeros(ns+1,ns); n = size(A,1); U = zeros(n,ns+1); U(:,1) = x/norm(x); delta = zeros(ns+1,1); for jj = 1:ns U(:,jj+1) = A*U(:,jj); % sólo se usa A aquí H(1:jj,jj) = U(:,1:jj)’*U(:,jj+1); U(:,jj+1) = U(:,jj+1) - U(:,1:jj)*H(1:jj,jj); % ortogonaliza delta(1:jj) = U(:,1:jj)’*U(:,jj+1); U(:,jj+1) = U(:,jj+1) - U(:,1:jj)*delta(1:jj); % reortogonaliza H(1:jj,jj) = H(1:jj,jj) + delta(1:jj); H(jj+1,jj) = norm(U(:,jj+1)); if H(jj+1,jj) == 0, flag = jj; return, end U(:,jj+1) = U(:,jj+1)/H(jj+1,jj); % arngo.m: script que demuestra las prestaciones de arnoldi_W end % end load west0479, A = west0479; % Matriz dispersa 479x479 n = size(A,1); ns = 30; % 30 iteraciones de Arnoldi randn(’state’,321) % mismo vector aleatorio cada vez x = randn(n,1); tic % vector de partida aleatorio [U,H,flag] = arnoldi_W(A,ns,x); if flag == 0 residual = norm(A*U(:,1:ns)-U*H) % comprobaciones precisión y ortogonalidad orthocheck = norm(eye(ns+1)-U’*U) ritz = eig(H(1:ns,1:ns),’nobalance’), toc, tic % valores de Ritz y tiempos lam = eig(full(A)); % cálculo de todos los val. propios de A toc, figure(1) plot(real(lam),imag(lam),’r+’,real(ritz),imag(ritz),’bo’) end
84/95
El resultado del ejemplo arngo.m con 30 iteraciones de Arnoldi es este. Las crucecitas rojas son los valores propios auténticos; los círculos azules los valores de Ritz que los aproximan. 2000 1500 1000 500 0 −500 −1000 −1500 −2000 −150
−100
−50
0
50
100
150
>> arngo residual = 8.432475148361903e-013 orthocheck = 1.251900907447703e-015 ritz = 1.0e+003 * 0.000009213557470 + 1.700662320900872i 0.000009213557470 - 1.700662320900872i -0.100885103219523 + 0.066606248354400i -0.100885103219523 - 0.066606248354400i -0.007240150911503 + 0.120672187071397i -0.007240150911503 - 0.120672187071397i 0.108125258880896 + 0.054065937658482i 0.108125258880896 - 0.054065937658482i -0.023312919813425 + 0.070655152714405i -0.023312919813425 - 0.070655152714405i 0.059637449793072 + 0.043628336812077i 0.059637449793072 - 0.043628336812077i 0.074697281694741 0.046089595120268 + 0.039266047684087i 0.046089595120268 - 0.039266047684087i 0.049326743009009 + 0.013259039094298i 0.049326743009009 - 0.013259039094298i 0.030463628327965 + 0.041942531078661i 0.030463628327965 - 0.041942531078661i 0.009119356708959 + 0.046886492625578i 0.009119356708959 - 0.046886492625578i -0.013773203484211 + 0.046711009454001i -0.013773203484211 - 0.046711009454001i -0.034755988534576 + 0.042746123558990i -0.034755988534576 - 0.042746123558990i -0.040459377125795 + 0.030311689517999i -0.040459377125795 - 0.030311689517999i -0.074653210517090 -0.047756691067739 + 0.011809670790806i -0.047756691067739 - 0.011809670790806i Elapsed time is 0.076879 seconds. Elapsed time is 0.226753 seconds.
85/95
Comparación de los métodos Cálculo de todos los valores propios y vectores propios Matrices generales reales o complejas. Reducción preliminar a forma de Hessenberg seguida de iteración QR. Matrices reales simétricas o complejas hermíticas. Reducción preliminar a matriz tridiagonal seguida de iteración QR. Cálculo de algunos valores y vectores propios Matrices simétricas de tamaño moderado. Reducción preliminar a matriz tridiagonal seguida de iteración inversa. Matrices de gran tamaño. Método con iteración de Arnoldi para matrices generales. Lanczos para matrices simétricas o complejas hermíticas.
86/95
Índice
Introducción teórica
Localización de valores propios
Obtención numérica
Cálculo de valores singulares
87/95
Cálculo de los valores singulares
El cálculo de los valores singulares de una matriz A, las raíces cuadradas no negativas de los valores propios de A T A se pueden obtener aplicando a A T A alguno de los métodos expuestos para obtener los valores propios de matrices simétricas.
Existen no obstante algoritmos más especializados basados en iteraciones QR de la matriz A. Veamos uno.
Algoritmo de Golub y Reinsch
88/95
Gene Howard Golub, EE.UU. 1932-2007.
Primera fase: Bidiagonalización
Consiste en transformar A 2 Cmn en triangular superior bidiagonal con transformaciones de Householder. Es decir, hacer donde
QB A˘ B D B D
B1 , 0
3 9 d1 f2 > > 7 > 6 d2 f3 0 = 7 > 6 6 ::: ::: 7 n 6 7 BD6 > 7 ::: f 7 > 6 0 > n7 > ; 6 4 5 dn m 0 2
y QB D Qn Q1 2 Rmm, ˘ B D ˘ 1 ˘ n
2
n
2 Rnn.
El esquema operativo que se sigue con una matriz A 64 es este: 0 0 0 0 0 0 0 0 0 0 Q 0 0 Q 0 0 Q 0 0 Q 0 0 0 1 0 ˘1 0 2 0 0 ˘2 0 0 3 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0000 0 0 0 0 0 0 0 0 0 0000 0 0 0 0 0 0 0 0 0
Esta rutina de Matlab consigue el resultado. function [Q1 A P1]= bidiag(A) % Se bidiagonaliza A por Householder. [m,n] = size(A); Q1=eye(m); P1=eye(n); for k = 1:min(m,n) % Introduce ceros debajo de la diagonal en la columna k. u = A(:,k); u(1:k-1) = 0; sigma = norm(u); if sigma ~= 0 if u(k) ~= 0, sigma = sign(u(k))*sigma; end u(k) = u(k) + sigma; rho = 1/(sigma’*u(k)); v = rho*(u’*A); q1=rho*(u’*Q1); A = A - u*v; Q1 = Q1 - u*q1; A(k+1:m,k) = 0; end % Introduce zeros a la derecha de la supradiagonal en la fila k. u = A(k,:); u(1:k) = 0; sigma = norm(u); if sigma ~= 0 if u(k+1) ~= 0, sigma = sign(u(k+1))*sigma; end u(k+1) = u(k+1) + sigma; rho = 1/(sigma’*u(k+1)); v = rho*(A*u’); p1=rho*(P1*u’); A = A - v*u; P1 = P1 - p1*u; A(k,k+2:n) = 0; end end end
89/95
Algoritmo de Golub y Reinsch. Segunda fase
90/95
Ya A bidiagonal, se diagonaliza mediante un algoritmo iterativo haciendo B kC1 D U Tk B k V k ;
k D 1; 2; : : : ;
con U k y V k ortogonales, y tal que lKımk!1 B k D ˙ D diag.1; : : : ; n/.
Cálculo práctico de la SVD (3)
La descomposición en valores singulares será U D QB diag.Ude y V D(3) ˘B V k . k ; I la m n/ culo práctico SVD
’
š kD:::;2;1
D de una matriz bidiagonal B, de tamaño n×n:
kD1;2:::
∗ 0 0 0 0
0 T=BTB.
En la figura se esquematiza cómo se hace. 6
∗ 0 0 0 ∗ ∗ ∗ 0 0 BG + 12 0 ∗ ∗ 0 → 0 0 0 ∗ ∗ 0 0 0 0 ∗ 0 ∗ 0 →60 0 0 BG 24
0 0 0 0 ∗ ∗ ∗ 0 0 G B 0 21 0 ∗ ∗ 0 → 0 0 0 ∗ ∗ 0 0 0 0 ∗ 0
∗ 0 0 0 ∗ ∗ ∗ 0 0 G B 0 43 0 ∗ ∗ 0 → 0 0 + ∗ ∗ 0 0 0 0 ∗ 0
∗ + 0 0 ∗ ∗ ∗ 0 0 BG 0 13 0 ∗ ∗ 0 → 0 0 0 ∗ ∗ 0 0 0 0 ∗ 0
∗ 0 0 0 ∗ ∗ ∗ 0 0 BG 0 35 0 ∗ ∗ + → 0 0 0 ∗ ∗ 0 0 0 0 ∗ 0
∗ 0 0 0 ∗ ∗ ∗ 0 0 G B 0 32 + ∗ ∗ 0 → 0 0 0 ∗ ∗ 0 0 0 0 ∗ 0
∗ 0 0 0 ∗ ∗ ∗ 0 0 G B 0 54 0 ∗ ∗ 0 → 0 0 0 ∗ ∗ 0 0 0 + ∗ 0
E
, donde
SVD de unaA matriz ˙ V TB, de tamaño n×n: D Ubidiagonal 0 ¾ Se puede aplicar el método QR con desplazamiento a la matriz T=BTB. ¾ También se pueden aplicar una serie (teóricamente infinita) de rotacion Givens a B para hacerla diagonal: ETSII-UPM
∗ ∗ 0 ∗ ∗ Se puede aplicar el método QR con desplazamiento a la matriz 0 0 ∗ de También se pueden aplicar una serie (teóricamente infinita) de rotaciones 0 0 0 Givens a B para hacerla diagonal: 0 0 0
0 0 ∗ 0 0 BG + 12 ∗ 0 → 0 ∗ ∗ 0 0 ∗ 0
∗ 0 0 0 ∗ ∗ ∗ + 0 BG 0 24 0 ∗ ∗ 0 → 0 0 0 ∗ ∗ 0 0 0 0 ∗ 0
0 0 0 0 ∗ ∗ ∗ 0 0 G B 0 21 0 ∗ ∗ 0 → 0 0 0 ∗ ∗ 0 0 0 0 ∗ 0
∗ 0 0 0 ∗ ∗ ∗ 0 0 G B 0 43 0 ∗ ∗ 0 → 0 0 + ∗ ∗ 0 0 0 0 ∗ 0
∗ + 0 0 ∗ ∗ ∗ 0 0 BG 0 13 0 ∗ ∗ 0 → 0 0 0 ∗ ∗ 0 0 0 0 ∗ 0
∗ 0 0 0 ∗ ∗ ∗ 0 0 BG 0 35 0 ∗ ∗ + → 0 0 0 ∗ ∗ 0 0 0 0 ∗ 0
∗ 0 0 0 ∗ ∗ ∗ 0 0 G B 0 32 + ∗ ∗ 0 → 0 0 0 ∗ ∗ 0 0 0 0 0 ∗
∗ 0 0 0 ∗ ∗ ∗ 0 0 G B 0 54 0 ∗ ∗ 0 → 0 0 0 ∗ ∗ 0 0 0 + ∗ 0
∗ 0 ∗ ∗
0 ∗
0 0 0 0
∗ 0 0 0 ∗ ∗ 0 0 0 ∗ ∗ 0 0 0 ∗ ∗ 0 0 0 ∗
∗ 0 0 0 ∗ ∗ 0 ¾0 Se comprueba que los elementos de la supra-diagonal tienden a cero 0 ∗ ∗ 0 n-1 hasta el primero) y los de la diagonal a los valores singulares. 0 0 ∗ ¾∗ La convergencia puede ser acelerada utilizando desplazamiento. 0 0 0 ∗
Recordemos las ideas el algoritmo de Jacobi para el cálculo de los valores propios de matrices simétricas, pues se hace algo parecido.
(d
91/95
Código de Matlab:
function [S U V it t t1 F]=svd_GolRein_1(A) % Se calcula al descomposición en valores singulares, [U S V]=U*S*V=A, % de una matriz cualquiera m x n. % Sigue casi exactamente el "paper" de Golub y Reinsch, "Singular Value % Decomposition and Least Squares Solutions", Numeri. Math., 14, 1970. c=1; s=1; % Comienzo iteración QR [m n]=size(A); tol=norm(A,’fro’)*sqrt(eps); itr=0; tic; for i=l+1:k if n>m, A=A’; [m n]=size(A); itr=1; end % Necesario m>=n g=e(i); y=q(i); h=s*g; g=c*g; % Bidiagonaliza A por Householder: U*B_n*V’ = A e(i-1)=sqrt(f*f+h*h); z=e(i-1); c=f/z; s=h/z; [U B_n V] = bidiag_Hansen(A); t1=toc; tic; f=x*c+g*s; g=-x*s+g*c; h=y*s; y=y*c; q=B_n(:,1); e=[0; B_n(1:n-1,2)]; U1=eye(m); V1=eye(n); V1(:,[i-1 i])=V1(:,[i-1 i])*[c -s;s c]; k=n; it=0; q(i-1)=sqrt(f*f+h*h); z=q(i-1); c=f/z; s=h/z; while k>0 % Bucle principal proceso itera. f=c*g+s*y; x=-s*g+c*y; caso=1; U1(:,[i-1 i])=U1(:,[i-1 i])*[c -s;s c]; for l=k:-1:1 end if abs(e(l)) U*S*V ans = 1.0000 2.0000 3.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000
>> A=[1 2 3;3 4 5;6 7 8]; >> [V S U]=svd(A) V = -0.2500 -0.8371 -0.4867 -0.4852 -0.3267 0.8111 -0.8379 0.4389 -0.3244 S = 14.5576 0 0 0 1.0372 0 0 0 0.0000 U = -0.4625 0.7870 0.4082 -0.5706 0.0882 -0.8165 -0.6786 -0.6106 0.4082 >> V*S*U’ ans = 1.0000 2.0000 3.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000
>> A=rand(1000); >> A = (A - .5)*2;
Tiempos y precisión
>> tic, [S U V]=svd_GolRein_1(A); toc Elapsed time is 13.971564 seconds. >> tic, [V S U]=svd(A); toc Elapsed time is 0.315793 seconds. >> norm(diag(S-S1)) ans = 1.1278e-08
92/95
93/95
1 0 0 0 2 Para una matriz complicada como M D 00 00 30 00 00 : 04000
>> M=[1 0 0 0 2;0 0 3 0 0;0 0 0 0 0;0 4 0 0 0] M = 1 0 0 0 2 0 0 3 0 0 0 0 0 0 0 0 4 0 0 0 >> [S U V]=svd_GolRein_1(M) S = 4.0000 0 0 0 0 0 3.0000 0 0 0 0 0 2.2361 0 0 0 0 0 0 0 U = 0 0 -1 0 0 -1 0 0 0 0 0 1 1 0 0 0 V = 0 1.0000 0 0 0 0 0 -1.0000 0 0 -0.4472 0 0 0 -0.8944 0 0 0 1.0000 0 0 0 0 0 1.0000
>> [U S V]=svd(M) U = 0 0 1 0 1 0 0 0 0 1 0 0 S = 4.0000 0 0 3.0000 0 0 0 0 V = 0 0 1.0000 0 0 1.0000 0 0 0 0
0 0 -1 0 0 0 2.2361 0
0 0 0 0
0 0 0 0
0.4472 0 0 0 0.8944
0 0 0 1.0000 0
-0.8944 0 0 0 0.4472
94/95
Algoritmo de Jacobi
Las raíces cuadradas no negativas de los valores propios de AA T son los valores singulares de A. Utilicemos Jacobi para calcular aquellos valores propios. >> A=[1 2 3;3 4 5;6 7 8]; >> [V D it]=Jacobi_val_12(A’*A) V = 0.7870 0.4082 -0.4625 0.0882 -0.8165 -0.5706 -0.6106 0.4082 -0.6786 D = 0.0000 1.0759 211.9241 it = 6 >> sqrt(D) ans = 0.0000 1.0372 14.5576
>> A=[1 2 3;3 4 5;6 7 8]; >> [V S U]=svd(A) V = -0.2500 -0.8371 -0.4867 -0.4852 -0.3267 0.8111 -0.8379 0.4389 -0.3244 S = 14.5576 0 0 0 1.0372 0 0 0 0.0000 U = -0.4625 0.7870 0.4082 -0.5706 0.0882 -0.8165 -0.6786 -0.6106 0.4082 >> V*S*U’ ans = 1.0000 2.0000 3.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000
Con la matriz M :
>> [V D it]=Jacobi_val_12(M’*M) V = 0.4472 0 0 0 -0.8944 0 1.0000 0 0 0 0 0 1.0000 0 0 0 0 0 1.0000 0 0.8944 0 0 0 0.4472 D = 0 0 5 9 16 it = 1 >> sqrt(D) ans = 0 0 2.2361 3.0000 4.0000
Algo más sofisticado, aunque en lo esencial igual: function [U S V] = svd_J(A,tol) % SVDJ Singular value decomposition using Jacobi algorithm. if nargin == 1, tol = eps; end % Tolerancia por defecto [M,N] = size(A); K = min(M,N); % K es el número de valores singulares On = 0; for c=A, On=On + sum(abs(c).^2); end On=On./N; % Suma coef.^2/N Previous_Off = Inf; V = eye(N); while true R = 0; % Contar las rotaciones for r = 1:N - 1 for c = r + 1:N % Calculate the three elements of the implicit matrix B that are % needed to calculate a Jacobi rotation. Since B is Hermitian, the % fourth element (b_cr) is not needed. b_rr = sum(abs(A(:,r)).^2); % Real value. b_cc = sum(abs(A(:,c)).^2); % Real value. b_rc = A(:,r)’ * A(:,c); % Same type as A. % Calculate a Jacobi rotation (four elements of G). The two values % that we calculate are a real value, C = cos(theta) and S, a value % of the same type as A, such that |S| = sin(theta). m = abs(b_rc); if m ~= 0 % If the off-diagonal element is zero, we don’t rotate. tau = (b_cc - b_rr)/(2*m); % tau is real and will be zero if % the two on-diagonal elements are % equal. In this case G will be an % identity matrix, and there is no % point in further calculating it. if tau ~= 0 R = R + 1; % Count the rotation we are about to perform. t = sign(tau)./(abs(tau) + sqrt(1+tau.^ 2)); C = 1./sqrt(1 + t.^ 2); S = (b_rc.* t.* C)./ m; % Initialize the rotation matrix, which is the same size as the % implicit matrix B. % We have to create an identity matrix here of the same type as A, % that is, quaternion if A is a quaternion, double if A is double. % To do this we use a function handle (q.v.) constructed from the % class type of A. This was done before the loop, since the type % of A is invariant. G = eye(N); G(r,r) = C; G(c,c) = C; G(r,c) = S; G(c,r) =-conj(S); A = A * G; V = V * G; end end end end
95/95
if R == 0, error(’No rotations performed during sweep.’), end % Calculate the sum of the off-diagonal elements of the matrix B. B = A’ * A; Off = sum(sum(abs(triu(B,1))))/(N.^2); % Normalise by the matrix size! if (Off/On) < tol, break; end % Off-diagonal sum is small enough to stop. if Previous_Off < Off warning(’QTFM:information’, ... ’Terminating sweeps: off diagonal sum increased on last sweep.’) break; end Previous_Off = Off; end % Extract and sort the singular values. The vector T may be longer than the % number of singular values (K) in cases where A is not square. [T,IX] = sort(sqrt(abs(diag(B))),’descend’); if nargout == 0 || nargout == 1 % .. only the singular values are needed. U = T(1:K); end if nargout == 3 % .. the singular vectors and singular values are needed. A = A(:, IX); % Map the columns of A and V into the same order as the V = V(:, IX); % singular values, using the sort indices in IX. % Construct the left singular vectors. These are in A but we need % to divide each column by the corresponding singular value. This % calculation is done by replicating T to make a matrix which can % then be divided into A element-by-element (vectorized division). U = A./ repmat(T’,M,1); S = diag(T); % Construct a diagonal matrix of singular values from % the vector T, because when there are three output % parameters, S is required to be a matrix. end end