QR-Zerlegung mit Householder-Transformationen · Householder-Transformation Das kann man durch...

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


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


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


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


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


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


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 = -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


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 = -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


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

QR-Zerlegung mit Householder-Transformationen · Householder-Transformation Das kann man durch...

Documents

Transcript of QR-Zerlegung mit Householder-Transformationen · Householder-Transformation Das kann man durch...