Post on 06-Mar-2018
transcript
QR-Zerlegung mitHouseholder-Transformationen
Numerische Mathematik 1WS 2011/12
1/33
Orthogonales Eliminieren
Sei x 2 Rn ein Vektor x 6= 0.
Ziel: Ein orthogonales H 2 Rn;n bestimmen, sodass
Hx = �kxke1;
ein Vielfaches des ersten Einheitsvektors
e1 =
26664
10...0
37775
ist.
2/33
Householder-Transformation
Das kann man durch Drehungen (Givens-Rotationen) oderdurch Spiegelungen (Householder-Transformationen)erreichen.
x2
x1
x
3/33
Householder-TransformationSetzt man
~w := x + sgn(x1)kxke1; und w :=~w
k ~wk
dann ist der Orthogonalprojektor auf hwi := span(w) gegebendurch
P := wwT =
�~w
k ~wk
��~wT
k ~wk
�=
~w ~wT
~wT ~w:
x2 +kxke1
wPx
x1
hwi?x hwi
4/33
Householder-TransformationDamit ist der (duale) Orthogonalprojektor auf hwi? gegebendurch
Q := I � P = I � wwT ;
und entsprechend die Householder-Spiegelung (an hwi?)
H := I � 2P = I � 2wwT :
x2 +kxke1
wPx
x1
hwi?
Hx
Qx
x hwi
5/33
Alternative DarstellungDie Householder-Transformation H 2 Rn;n hat also dieDarstellung
H = I � 2wwT = I � 2wwT
wT w;
wobei w 2 Rn mit kwk = 1.
Ist v 2 Rn n f0g aber ein linear abhangiger Vektor v = � �w , mit� 2 R n f0g so gilt ebenfalls
I � 2vvT
vT v= I � 2
�2wwT
�2wT w
= I � 2wwT
wT w= H:
6/33
Normalisierte Darstellung
Ist x 6= 0 so ist der erste Eintrag von
~w := x + sgn(x1)kxke1
ebenfalls nicht 0, d.h. ~w1 6= 0.
Dann erfullt der Vektor
v :=1~w1
� ~w 2 Rn
die Bedingung v1 = 1.
7/33
Normalisierte Darstellung
Somit ist die Householder-Transformation gegeben durch
H = I � 2P = I � 2vvT
vT v= I �
2kvk2 vvT
= I � �vvT ;
wobei� :=
2kvk2 :
8/33
Berechnung von v und �
function [beta , v] = house(x)
n = length(x);
if( norm(x) < n*eps )
beta = 0;
v = eye(n,1);
else
% bereche w
w = ...
% berechne v
v = w/w(1);
% berechne beta
beta = ...
end
9/33
VorteilDie Householder-Transformation kann man also als
H = I � �vvT
schreiben, wobei � 2 R und v 2 Rn die Form
v =
26664
1v2...
vn
37775
hat.
H I � �vvT
Speicherbedarf: n2 n
10/33
H anwendenEs sei � 2 R, v 2 Rn und die Householder-Transformation
H = I � �vvT 2 Rn;n:
Dann ist fur X 2 Rn;` gerade
HX = (I � �vvT )X = X � �vvT X
= X � �v�
vT X�
| {z }=:X12R
1;` in O(n`)
= X � �vX1
= X � � (vX1)| {z }=:X22R
n;` in O(n`)
= X � �X2| {z }in O(n`)
=: Y :
H I � �vvT
Laufzeit, Matrix-Vektor-Multiplikation: O(n2) O(n)Laufzeit, Matrix-Matrix-Multiplikation: O(n3) O(n2)
11/33
H anwenden in Matlab
function Y = h_anwenden(v, beta , X)
Y = zeros(size(X));
% X1 in der ersten Zeile von Y speichern
% X2 in Y speichern
% Y ausrechnen
Y = ...
12/33
QR-Zerlegung (Formal)
Wir wollen eine QR-Zerlegung der Matrix
A =
26666664
a(1)11 a(1)
12 a(1)13
a(1)21 a(1)
22 a(1)23
a(1)31 a(1)
32 a(1)33
a(1)41 a(1)
42 a(1)43
a(1)51 a(1)
52 a(1)53
377777752 R5;3
mittels Householder-Transformationen berechnen.
13/33
QR-Zerlegung (Formal)
A =
26666664
a(1)11 a(1)
12 a(1)13
a(1)21 a(1)
22 a(1)23
a(1)31 a(1)
32 a(1)33
a(1)41 a(1)
42 a(1)43
a(1)51 a(1)
52 a(1)53
37777775
Schritt 1: Berechne v1 2 R5 und �1 sodass mit
H1 := I � �1v1vT1 2 R5;5
und
H1 := H1 gilt H1A =
26666664
a(2)11 a(2)
12 a(2)13
0 a(2)22 a(2)
230 a(2)
32 a(2)33
0 a(2)42 a(2)
430 a(2)
52 a(2)53
37777775:
14/33
QR-Zerlegung (Formal)
H1A =
26666664
a(2)11 a(2)
12 a(2)13
0 a(2)22 a(2)
230 a(2)
32 a(2)33
0 a(2)42 a(2)
430 a(2)
52 a(2)53
37777775
Schritt 2: Berechne v2 2 R4 und �2 sodass mit
H2 := I � �2v2vT2 2 R4;4
und
H2 :=
�1
H2
�gilt H2H1A =
26666664
a(3)11 a(3)
12 a(3)13
0 a(3)22 a(3)
230 0 a(3)
330 0 a(3)
430 0 a(3)
53
37777775:
15/33
QR-Zerlegung (Formal)
H2H1A =
26666664
a(3)11 a(3)
12 a(3)13
0 a(3)22 a(3)
230 0 a(3)
330 0 a(3)
430 0 a(3)
53
37777775
Schritt 3: Berechne v3 2 R3 und �3 sodass mit
H3 := I � �3v3vT3 2 R3;3
und
H3 :=
241
1H3
35 gilt H3H2H1A =
2666664
a(4)11 a(4)
12 a(4)13
0 a(4)22 a(4)
230 0 a(4)
330 0 00 0 0
3777775 :
16/33
QR-Zerlegung (Formal)
H3H2H1| {z }=:QT
A =
2666664
a(4)11 a(4)
12 a(4)13
0 a(4)22 a(4)
230 0 a(4)
330 0 00 0 0
3777775 =: R:
Dann ist Q orthogonal, R obere Dreiecksmatrix und somit
A = QReine QR-Zerlegung von A.
17/33
QR-Zerlegung (im Rechner)
Es werden die Hi , Hi , etc. nicht explizit gebildet.
Man rechnet nur sog. Elementarreflektoren v1, v2, ... aus undspeichert diese in den frei werdenden Null-Elementen (dieerste 1 muss nicht gespeichert werden).
Dabei braucht man ein extra Array fur die �1, �2, ...
18/33
QR-Zerlegung (im Rechner)Mit obigen Bezeichnungen durchlauft man die Abfolge
2666664
a(1)11 a(1)
12 a(1)13
a(1)21 a(1)
22 a(1)23
a(1)31 a(1)
32 a(1)33
a(1)41 a(1)
42 a(1)43
a(1)51 a(1)
52 a(1)53
3777775
� =
�0 0 0
�
2666664
a(2)11 a(2)
12 a(2)13
v1;2 a(2)22 a(2)
23
v1;3 a(2)32 a(2)
33
v1;4 a(2)42 a(2)
43
v1;5 a(2)52 a(2)
53
3777775
� =
��1 0 0
�
2666664
a(3)11 a(3)
12 a(3)13
v1;2 a(3)22 a(3)
23
v1;3 v2;2 a(3)33
v1;4 v2;3 a(3)43
v1;5 v2;4 a(3)53
3777775
� =
��1 �2 0
�
2666664
a(4)11 a(4)
12 a(4)13
v1;2 a(4)22 a(4)
23
v1;3 v2;2 a(4)33
v1;4 v2;3 v3;2v1;5 v2;4 v3;3
3777775
� =
��1 �2 �3
�
wobei vi;j := vi
���j
den j-ten Eintrag von vi bezeichnet.19/33
QR-Zerlegung (in LAPACK)SUBROUTINE DGEQRF( M, N, A, LDA , TAU , WORK , LWORK , INFO )
** .. Scalar Arguments ..
INTEGER INFO , LDA , LWORK , M, N* ..* .. Array Arguments ..
DOUBLE PRECISION A( LDA , * ), TAU( * ), WORK( * )* ..** Purpose* =======** DGEQRF computes a QR factorization of a real M-by-N matrix A:* A = Q * R.** Arguments* =========** M (input) INTEGER* The number of rows of the matrix A. M >= 0.** N (input) INTEGER* The number of columns of the matrix A. N >= 0.** A (input/output) DOUBLE PRECISION array , dimension (LDA ,N)* On entry , the M-by -N matrix A.* On exit , the elements on and above the diagonal of the array* contain the min(M,N)-by-N upper trapezoidal matrix R (R is* upper triangular if m >= n); the elements below the diagonal ,* with the array TAU , represent the orthogonal matrix Q as a* product of min(m,n) elementary reflectors (see Further* Details ).** LDA (input) INTEGER* The leading dimension of the array A. LDA >= max(1,M).** TAU (output) DOUBLE PRECISION array , dimension (min(M,N))* The scalar factors of the elementary reflectors (see Further* Details ).** WORK (workspace/output) DOUBLE PRECISION array , dimension (MAX(1,LWORK))* On exit , if INFO = 0, WORK (1) returns the optimal LWORK.** LWORK (input) INTEGER* The dimension of the array WORK. LWORK >= max(1,N).* For optimum performance LWORK >= N*NB, where NB is* the optimal blocksize.** If LWORK = -1, then a workspace query is assumed; the routine* only calculates the optimal size of the WORK array , returns* this value as the first entry of the WORK array , and no error* message related to LWORK is issued by XERBLA.** INFO (output) INTEGER* = 0: successful exit* < 0: if INFO = -i, the i-th argument had an illegal value** Further Details* ===============** The matrix Q is represented as a product of elementary reflectors** Q = H(1) H(2) . . . H(k), where k = min(m,n).** Each H(i) has the form** H(i) = I - tau * v * v**T** where tau is a real scalar , and v is a real vector with* v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),* and tau in TAU(i).** =====================================================================** .. Local Scalars ..
LOGICAL LQUERYINTEGER I, IB , IINFO , IWS , K, LDWORK , LWKOPT , NB,
$ NBMIN , NX* ..* .. External Subroutines ..
EXTERNAL DGEQR2 , DLARFB , DLARFT , XERBLA* ..* .. Intrinsic Functions ..
INTRINSIC MAX , MIN* ..* .. External Functions ..
INTEGER ILAENVEXTERNAL ILAENV
* ..* .. Executable Statements ..** Test the input arguments*
INFO = 0NB = ILAENV( 1, 'DGEQRF ', ' ', M, N, -1, -1 )LWKOPT = N*NBWORK( 1 ) = LWKOPTLQUERY = ( LWORK.EQ.-1 )IF( M.LT.0 ) THEN
INFO = -1ELSE IF( N.LT.0 ) THEN
INFO = -2ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
INFO = -7END IFIF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGEQRF ', -INFO )RETURN
ELSE IF( LQUERY ) THENRETURN
END IF** Quick return if possible*
K = MIN( M, N )IF( K.EQ.0 ) THEN
WORK( 1 ) = 1RETURN
END IF*
NBMIN = 2NX = 0IWS = NIF( NB.GT.1 .AND. NB.LT.K ) THEN
** Determine when to cross over from blocked to unblocked code.*
NX = MAX( 0, ILAENV( 3, 'DGEQRF ', ' ', M, N, -1, -1 ) )IF( NX.LT.K ) THEN
** Determine if workspace is large enough for blocked code.*
LDWORK = NIWS = LDWORK*NBIF( LWORK.LT.IWS ) THEN
** Not enough workspace to use optimal NB: reduce NB and* determine the minimum value of NB.*
NB = LWORK / LDWORKNBMIN = MAX( 2, ILAENV( 2, 'DGEQRF ', ' ', M, N, -1,
$ -1 ) )END IF
END IFEND IF
*IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
** Use blocked code initially*
DO 10 I = 1, K - NX, NBIB = MIN( K-I+1, NB )
** Compute the QR factorization of the current block* A(i:m,i:i+ib -1)*
CALL DGEQR2( M-I+1, IB , A( I, I ), LDA , TAU( I ), WORK ,$ IINFO )
IF( I+IB.LE.N ) THEN** Form the triangular factor of the block reflector* H = H(i) H(i+1) . . . H(i+ib -1)*
CALL DLARFT( 'Forward ', 'Columnwise ', M-I+1, IB ,$ A( I, I ), LDA , TAU( I ), WORK , LDWORK )
** Apply H**T to A(i:m,i+ib:n) from the left*
CALL DLARFB( 'Left', 'Transpose ', 'Forward ',$ 'Columnwise ', M-I+1, N-I-IB+1, IB,$ A( I, I ), LDA , WORK , LDWORK , A( I, I+IB ),$ LDA , WORK( IB+1 ), LDWORK )
END IF10 CONTINUE
ELSEI = 1
END IF** Use unblocked code to factor the last or only block.*
IF( I.LE.K )$ CALL DGEQR2( M-I+1, N-I+1, A( I, I ), LDA , TAU( I ), WORK ,$ IINFO )
*WORK( 1 ) = IWSRETURN
** End of DGEQRF*
END
20/33
QR-Zerlegung in Matlabfunction [QR , betas] = qr_zerlegung(A)
[n, m] = size(A);
min_nm = min([n-1, m]);
betas = zeros(1, min_nm );
for i=1: min_nm
[beta_i , v_i] = house (...);
% A^(i) --> A^(i+1)
... h_anwenden (...) ...
A(i+1:n,i) = v_i (2:end);
betas(i) = beta_i;
end
QR = A;
21/33
ProblemHat man die QR-Zerlegung in der Form2
666664
a(4)11 a(4)
12 a(4)13
v1;2 a(4)22 a(4)
23v1;3 v2;2 a(4)
33v1;4 v2;3 v3;2v1;5 v2;4 v3;3
3777775
� =��1 �2 �3
�gespeichert, so kann man
Q = H1H2H3
und ebenso (da Householder-Transformationen symmetrischsind)
QT = H3H2H1
nicht so einfach mit einem Vektor (oder mit einer Matrix)Multiplizieren.
22/33
Losungfunction Y = q_anwenden(QR , betas , X, transponiert)
if( transponiert )
% hier soll Q^T X ausgerechnet werden
for i=...
... h_anwenden (...) ...
end
else
% hier soll Q X ausgerechnet werden
for i=...
... h_anwenden (...) ...
end
end
23/33
Q-Anwenden (in LAPACK)SUBROUTINE DORMQR( SIDE , TRANS , M, N, K, A, LDA , TAU , C, LDC ,
$ WORK , LWORK , INFO )** Purpose* =======** DORMQR overwrites the general real M-by -N matrix C with** SIDE = 'L' SIDE = 'R'* TRANS = 'N': Q * C C * Q* TRANS = 'T': Q**T * C C * Q**T** where Q is a real orthogonal matrix defined as the product of k* elementary reflectors** Q = H(1) H(2) . . . H(k)** as returned by DGEQRF. Q is of order M if SIDE = 'L' and of order N* if SIDE = 'R'.** Arguments* =========** TRANS (input) CHARACTER *1* = 'N': No transpose , apply Q;* = 'T': Transpose , apply Q**T.** A (input) DOUBLE PRECISION array , dimension (LDA ,K)* The i-th column must contain the vector which defines the* elementary reflector H(i), for i = 1,2,...,k, as returned by* DGEQRF in the first k columns of its array argument A.* A is modified by the routine but restored on exit.** SIDE (input) CHARACTER *1* = 'L': apply Q or Q**T from the Left;* = 'R': apply Q or Q**T from the Right.** .. Scalar Arguments ..
CHARACTER SIDE , TRANSINTEGER INFO , K, LDA , LDC , LWORK , M, N
* ..* .. Array Arguments ..
DOUBLE PRECISION A( LDA , * ), C( LDC , * ), TAU( * ), WORK( * )* ..** M (input) INTEGER* The number of rows of the matrix C. M >= 0.** N (input) INTEGER* The number of columns of the matrix C. N >= 0.** K (input) INTEGER* The number of elementary reflectors whose product defines* the matrix Q.* If SIDE = 'L', M >= K >= 0;* if SIDE = 'R', N >= K >= 0.** LDA (input) INTEGER* The leading dimension of the array A.* If SIDE = 'L', LDA >= max(1,M);* if SIDE = 'R', LDA >= max(1,N).** TAU (input) DOUBLE PRECISION array , dimension (K)* TAU(i) must contain the scalar factor of the elementary* reflector H(i), as returned by DGEQRF.** C (input/output) DOUBLE PRECISION array , dimension (LDC ,N)* On entry , the M-by -N matrix C.* On exit , C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.** LDC (input) INTEGER* The leading dimension of the array C. LDC >= max(1,M).** WORK (workspace/output) DOUBLE PRECISION array , dimension (MAX(1,LWORK))* On exit , if INFO = 0, WORK (1) returns the optimal LWORK.** LWORK (input) INTEGER* The dimension of the array WORK.* If SIDE = 'L', LWORK >= max(1,N);* if SIDE = 'R', LWORK >= max(1,M).* For optimum performance LWORK >= N*NB if SIDE = 'L', and* LWORK >= M*NB if SIDE = 'R', where NB is the optimal* blocksize.** If LWORK = -1, then a workspace query is assumed; the routine* only calculates the optimal size of the WORK array , returns* this value as the first entry of the WORK array , and no error* message related to LWORK is issued by XERBLA.** INFO (output) INTEGER* = 0: successful exit* < 0: if INFO = -i, the i-th argument had an illegal value** =====================================================================** .. Parameters ..
INTEGER NBMAX , LDTPARAMETER ( NBMAX = 64, LDT = NBMAX +1 )
* ..* .. Local Scalars ..
LOGICAL LEFT , LQUERY , NOTRANINTEGER I, I1 , I2, I3, IB , IC, IINFO , IWS , JC , LDWORK ,
$ LWKOPT , MI , NB, NBMIN , NI , NQ , NW* ..* .. Local Arrays ..
DOUBLE PRECISION T( LDT , NBMAX )* ..* .. External Functions ..
LOGICAL LSAMEINTEGER ILAENVEXTERNAL LSAME , ILAENV
* ..* .. External Subroutines ..
EXTERNAL DLARFB , DLARFT , DORM2R , XERBLA* ..* .. Intrinsic Functions ..
INTRINSIC MAX , MIN* ..* .. Executable Statements ..** Test the input arguments*
INFO = 0LEFT = LSAME( SIDE , 'L' )NOTRAN = LSAME( TRANS , 'N' )LQUERY = ( LWORK.EQ.-1 )
** NQ is the order of Q and NW is the minimum dimension of WORK*
IF( LEFT ) THENNQ = MNW = N
ELSENQ = NNW = M
END IFIF( .NOT.LEFT .AND. .NOT.LSAME( SIDE , 'R' ) ) THEN
INFO = -1ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS , 'T' ) ) THEN
INFO = -2ELSE IF( M.LT.0 ) THEN
INFO = -3ELSE IF( N.LT.0 ) THEN
INFO = -4ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN
INFO = -5ELSE IF( LDA.LT.MAX( 1, NQ ) ) THEN
INFO = -7ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
INFO = -10ELSE IF( LWORK.LT.MAX( 1, NW ) .AND. .NOT.LQUERY ) THEN
INFO = -12END IF
*IF( INFO.EQ.0 ) THEN
** Determine the block size. NB may be at most NBMAX , where NBMAX* is used to define the local array T.*
NB = MIN( NBMAX , ILAENV( 1, 'DORMQR ', SIDE // TRANS , M, N, K,$ -1 ) )
LWKOPT = MAX( 1, NW )*NBWORK( 1 ) = LWKOPT
END IF*
IF( INFO.NE.0 ) THENCALL XERBLA( 'DORMQR ', -INFO )RETURN
ELSE IF( LQUERY ) THENRETURN
END IF** Quick return if possible*
IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) THENWORK( 1 ) = 1RETURN
END IF*
NBMIN = 2LDWORK = NWIF( NB.GT.1 .AND. NB.LT.K ) THEN
IWS = NW*NBIF( LWORK.LT.IWS ) THEN
NB = LWORK / LDWORKNBMIN = MAX( 2, ILAENV( 2, 'DORMQR ', SIDE // TRANS , M, N, K,
$ -1 ) )END IF
ELSEIWS = NW
END IF*
IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN** Use unblocked code*
CALL DORM2R( SIDE , TRANS , M, N, K, A, LDA , TAU , C, LDC , WORK ,$ IINFO )ELSE
** Use blocked code*
IF( ( LEFT .AND. .NOT.NOTRAN ) .OR.$ ( .NOT.LEFT .AND. NOTRAN ) ) THEN
I1 = 1I2 = KI3 = NB
ELSEI1 = ( ( K-1 ) / NB )*NB + 1I2 = 1I3 = -NB
END IF*
IF( LEFT ) THENNI = NJC = 1
ELSEMI = MIC = 1
END IF*
DO 10 I = I1, I2 , I3IB = MIN( NB, K-I+1 )
** Form the triangular factor of the block reflector* H = H(i) H(i+1) . . . H(i+ib -1)*
CALL DLARFT( 'Forward ', 'Columnwise ', NQ -I+1, IB, A( I, I ),$ LDA , TAU( I ), T, LDT )
IF( LEFT ) THEN** H or H**T is applied to C(i:m,1:n)*
MI = M - I + 1IC = I
ELSE** H or H**T is applied to C(1:m,i:n)*
NI = N - I + 1JC = I
END IF** Apply H or H**T*
CALL DLARFB( SIDE , TRANS , 'Forward ', 'Columnwise ', MI , NI ,$ IB , A( I, I ), LDA , T, LDT , C( IC, JC ), LDC ,$ WORK , LDWORK )
10 CONTINUEEND IFWORK( 1 ) = LWKOPTRETURN
** End of DORMQR*
END
24/33
Lineares Ausgleichsproblem
Sei A 2 Rn;m eine Matrix mit vollem Spaltenrang und b 2 Rn.Meistens ist auch n � m.
Ziel: Das lineare Ausgleichsproblem losen, d.h. ein x� 2 Rm
bestimmen mit
minx2Rm
kAx � bk = kAx� � bk :
25/33
Lineares AusgleichsproblemHat man eine QR-Zerlegung von A 2 Rn;m berechnet so ist
minx2Rm
kAx � bk = minx2Rm
kQRx � bk = minx2Rm
Q(Rx �QT b)
= minx2Rm
Rx �QT b ;
d.h. mit den Bezeichnungen R =:
�R10
�und QT b =:
�b1b2
�,
wobei R1 2 Rm;m, b1 2 R
m, b2 2 Rn�m, hat man
minx2Rm
kAx � bk2 = minx2Rm
�R10
�x �
�b1b2
� 2
= minx2Rm
kR1x � b1k2 + k0� b2k
2 = kb2k+ minx2Rm
kR1x � b1k2;
mit der Losungx� = R�1b1
und dem Residuum
kAx� � bk = kb2k: 26/33
Ausgleichsproblem losen (in Matlab)
function [x, res] = ausgleichsproblem(A, b)
[n, m] = size(A);
[QR, beta] = qr_zerlegung(A);
QTb = q_anwenden(QR , beta , b, true);
...
x = ...;
res = ...;
27/33
Asymptotische LaufzeitEs sei m 2 N fest, n 2 N mit n > m variabel, A 2 Rn;m undb 2 Rn.
Dann ist die asymptotische Laufzeit (fur n !1) zur Losungdes linearen Ausgleichsproblems
minx2Rm
kAx � bk ;
nach obiger MethodeO(n):
Wurde man erst orthogonales Q 2 Rn;n und eine obereDreiecksmatrix R 2 Rn;m berechnen mit QR = A, so wurdeschon die notwendige Matrix-Vektor-Multiplikation QT b
O(n2)
brauchen.28/33
Berechnung der SVD
Um die SVD einer Matrix der Form
A =
2666666664
� � � �� � � �� � � �� � � �� � � �� � � �� � � �
3777777775
zu berechnen, wendet man Householder-Transformationen(von links und rechts) an.
30/33
Berechnung der SVDMan durchlauft die Abfolge
A
2666666664
� � � �0 � � �0 � � �0 � � �0 � � �0 � � �0 � � �
3777777775
2666666664
� � 0 00 � � �0 � � �0 � � �0 � � �0 � � �0 � � �
3777777775
2666666664
� � 0 00 � � �0 0 � �0 0 � �0 0 � �0 0 � �0 0 � �
3777777775
2666666664
� � 0 00 � � 00 0 � �0 0 � �0 0 � �0 0 � �0 0 � �
3777777775
2666666664
� � 0 00 � � 00 0 � �0 0 0 �0 0 0 �0 0 0 �0 0 0 �
3777777775
2666666664
� � 0 00 � � 00 0 � �0 0 0 �0 0 0 00 0 0 00 0 0 0
3777777775=: B
31/33
Berechnung der SVDDann kann man die Eigenwerte von
BT B =
2664� 0 0 0 0 0 0� � 0 0 0 0 00 � � 0 0 0 00 0 � � 0 0 0
3775
2666666664
� � 0 00 � � 00 0 � �0 0 0 �0 0 0 00 0 0 00 0 0 0
3777777775
=
2664� � 0 0� � � 00 � � �0 0 � �
3775 ;
einer symmetrischen positiv semi-definiten Tridiagonalmatrix,berechnen.
32/33
Eigenwerte von BT B (in LAPACK)SUBROUTINE DSTEV( JOBZ , N, D, E, Z, LDZ , WORK , INFO )
** .. Scalar Arguments ..
CHARACTER JOBZINTEGER INFO , LDZ , N
* .. Array Arguments ..DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ , * )
** Purpose* =======** DSTEV computes all eigenvalues and , optionally , eigenvectors of a* real symmetric tridiagonal matrix A.** Arguments* =========** JOBZ (input) CHARACTER *1* = 'N': Compute eigenvalues only;* = 'V': Compute eigenvalues and eigenvectors.** N (input) INTEGER* The order of the matrix. N >= 0.** D (input/output) DOUBLE PRECISION array , dimension (N)* On entry , the n diagonal elements of the tridiagonal matrix* A.* On exit , if INFO = 0, the eigenvalues in ascending order.** E (input/output) DOUBLE PRECISION array , dimension (N-1)* On entry , the (n-1) subdiagonal elements of the tridiagonal* matrix A, stored in elements 1 to N-1 of E.* On exit , the contents of E are destroyed.** Z (output) DOUBLE PRECISION array , dimension (LDZ , N)* If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal* eigenvectors of the matrix A, with the i-th column of Z* holding the eigenvector associated with D(i).* If JOBZ = 'N', then Z is not referenced.** LDZ (input) INTEGER* The leading dimension of the array Z. LDZ >= 1, and if* JOBZ = 'V', LDZ >= max(1,N).** WORK (workspace) DOUBLE PRECISION array , dimension (max(1,2*N-2))* If JOBZ = 'N', WORK is not referenced.** INFO (output) INTEGER* = 0: successful exit* < 0: if INFO = -i, the i-th argument had an illegal value* > 0: if INFO = i, the algorithm failed to converge; i* off -diagonal elements of E did not converge to zero.** =====================================================================** .. Parameters ..
DOUBLE PRECISION ZERO , ONEPARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
* ..* .. Local Scalars ..
LOGICAL WANTZINTEGER IMAX , ISCALEDOUBLE PRECISION BIGNUM , EPS , RMAX , RMIN , SAFMIN , SIGMA , SMLNUM ,
$ TNRM* ..* .. External Functions ..
LOGICAL LSAMEDOUBLE PRECISION DLAMCH , DLANSTEXTERNAL LSAME , DLAMCH , DLANST
* ..* .. External Subroutines ..
EXTERNAL DSCAL , DSTEQR , DSTERF , XERBLA* ..* .. Intrinsic Functions ..
INTRINSIC SQRT* ..* .. Executable Statements ..** Test the input parameters.*
WANTZ = LSAME( JOBZ , 'V' )*
INFO = 0IF( .NOT.( WANTZ .OR. LSAME( JOBZ , 'N' ) ) ) THEN
INFO = -1ELSE IF( N.LT.0 ) THEN
INFO = -2ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
INFO = -6END IF
*IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DSTEV ', -INFO )RETURN
END IF** Quick return if possible*
IF( N.EQ.0 )$ RETURN
*IF( N.EQ.1 ) THEN
IF( WANTZ )$ Z( 1, 1 ) = ONE
RETURNEND IF
** Get machine constants.*
SAFMIN = DLAMCH( 'Safe minimum ' )EPS = DLAMCH( 'Precision ' )SMLNUM = SAFMIN / EPSBIGNUM = ONE / SMLNUMRMIN = SQRT( SMLNUM )RMAX = SQRT( BIGNUM )
** Scale matrix to allowable range , if necessary.*
ISCALE = 0TNRM = DLANST( 'M', N, D, E )IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
ISCALE = 1SIGMA = RMIN / TNRM
ELSE IF( TNRM.GT.RMAX ) THENISCALE = 1SIGMA = RMAX / TNRM
END IFIF( ISCALE.EQ.1 ) THEN
CALL DSCAL( N, SIGMA , D, 1 )CALL DSCAL( N-1, SIGMA , E( 1 ), 1 )
END IF** For eigenvalues only , call DSTERF. For eigenvalues and* eigenvectors , call DSTEQR.*
IF( .NOT.WANTZ ) THENCALL DSTERF( N, D, E, INFO )
ELSECALL DSTEQR( 'I', N, D, E, Z, LDZ , WORK , INFO )
END IF** If matrix was scaled , then rescale eigenvalues appropriately.*
IF( ISCALE.EQ.1 ) THENIF( INFO.EQ.0 ) THEN
IMAX = NELSE
IMAX = INFO - 1END IFCALL DSCAL( IMAX , ONE / SIGMA , D, 1 )
END IF*
RETURN** End of DSTEV*
END
33/33