Loading...

1. A matrix Q ∈ Rn×n is called an orthogonal matrix if Qt Q = I. An orthogonal matrix Q and an upper triangular matrix R is called a QR factorization of A if A = QR. The Frobenius norm of a matrix A, denoted by ∥A∥F , is defined as ∥A∥F =

n n ∑ (∑

|Ai,j |2

)1/2 .

i=1 j=1

The n × n Hilbert H matrix is defined as the matrix with entries Hij = 1/(i + j − 1)

for

i, j = 1, 2, . . . , n.

(i) Suppose A is a non-singular n × n matrix with columns denoted by ui . Apply Gram–Schmidt orthogonalization to the ui ’s to obtain an orthonormal basis vi . Let Q be the matrix with columns consisting of the vectors vi and let R = Qt A. Prove the matrix Q is orthogonal and that R is upper triangular. Clearly Q ∈ Rn×n . Since the vi ’s form an orthonormal basis then { 1 for i = j (vi , vj ) = 0 for i ̸= j. It follows that

v1t

. [ t Q Q = .. v1

(v1 , vn ) .. = I. .

··· .. . (vn , v1 ) · · ·

(v1 , v1 ) ] .. vn = .

···

vnt

(vn , vn )

Therefore Q is an orthogonal matrix. To see that R is upper triangular recall the Gram–Schmidt algorithm ensures that span{ ui : i = 1, . . . , k } = span{ vi : i = 1, . . . , k }

for

k = 1, . . . , n.

Therefore (vi , uj ) = 0 for i > j. The fact that R is upper triangular now follows from

v1t

. [ t . R=QA= . u1

···

(v1 , u1 ) ] .. un = . (vn , u1 )

vnt 1

··· .. . ···

(v1 , un ) .. . . (vn , un )

(∗)

Math/CS 466/666 Programming Project 2 (ii) Suppose Q is an orthogonal matrix. Show that QQt = I. Since Q ∈ Rn×n then Qt Q = I implies the columns of Q are orthonormal vectors. Since there are n columns, the columns of Q must form an orthonormal basis of Rn . Suppose y ∈ Rn . As the columns of Q form an orthonormal basis, there is x ∈ Rn such that y = Qx. It follows that Qt y = Qt Qx = Ix = x. Thus, Qt y = 0 implies x = 0 and consequently that y = Qx = 0. Therefore Qt y = 0 if any only if y = 0. Now, since I = Qt Q then Qt I = Qt = IQt = (Qt Q)Qt = Qt (QQt ). Consequently Qt (I − QQt )x = 0

x ∈ Rn .

for every

As y = 0 is the only value of y such that Qt y = 0, it follows that (I − QQt )x = 0 for every x ∈ Rn . This is identical as saying QQt = I. (iii) [Extra Credit and Math/CS 666] Show that H is non-singular for any n. Define

∫

1

(f, g) =

f (x)g(x)dx

and

∥f ∥ =

√

(f, f ).

0

Note for any two polynomials p and q that ∥p − q∥ = 0 if and only if p = q. Therefore, the above definition yields a norm and inner product on the space of polynomials. Let a ∈ Rn and consider the polynomial p(x) = a1 + a2 x + · · · + an xn−1 =

n ∑

ai xi−1 .

i=1

Now

∫ ∥p∥ = 2

∫

1

0

=

1

p(x)p(x)dx = 0

n ∑ n ∫ 1 ∑ i=1 j=1

n (∑

ai xk−1 )

i=1

ai aj xi+j−2 dx =

0

n (∑

) aj xj−1 dx

j=1 n ∑ n ∑ i=1 j=1

ai aj dx = a · Ha. i+j−1

Therefore, if Ha = 0 then ∥p∥ = 0. Since this is a norm, that would imply p = 0, which in turn implies a = 0. It follows that H is injective and consequently invertible.

2

Math/CS 466/666 Programming Project 2 2. Let A be a matrix with columns given by ui for i = 1, . . . , n. Consider the following modification of the Gram–Schmidt orthogonalization algorithm: t1,1 t1,2 t2,2 t1,3 t2,3 t3,3 t1,4 t2,4 t3,4 t4,4

v1 = t1,1 /∥t1,1 ∥

= u1 = u2 = t1,2 − (v1 , t1,2 )v1 = u3 = t1,3 − (v1 , t1,3 )v1 = t2,3 − (v2 , t2,3 )v2 = u4 = t1,4 − (v1 , t1,4 )v1 = t2,4 − (v2 , t2,4 )v2 = t3,4 − (v3 , t3,4 )v3 .. .

v2 = t2,2 /∥t2,2 ∥

v3 = t3,3 /∥t3,3 ∥

v4 = t4,4 /∥t4,4 ∥

t1,n = un t2,n = t1,n − (v1 , t1,n )v1 t3,n = t2,n − (v2 , t2,n )v2 ... tn,n = tn−1,n − (vn−1 , tn−1,n )vn−1 Let Q be the matrix with columns vi and let (vi , ti,j ) Ri,j = ∥ti,i ∥ 0

vn = tn,n /∥tn,n ∥.

R be the matrix with entries for i < j for i = j for i > j.

(i) Show that ∥ti,i ∥ = (vi , ti,i ) and (vi , ti,j ) = (vi , uj ) for i ≤ j. Since vi = ti,i /∥ti,i ∥ by definition it follows that (vi , ti,i ) =

(ti,i , ti,i ) ∥ti,i ∥2 = = ∥ti,i ∥. ∥ti,i ∥ ∥ti,i ∥

Suppose there exists i < j such that vi · vj ̸= 0. Let j0 to be the smallest value for which such an i exists and further choose i0 to be the smallest value of i such that vi · vj0 ̸= 0. Write j0 = i0 + k. Since ( ) (vi , ti+1,i+1 ) = vi , ti,i+1 − (vi , ti,i+1 )vi = (vi , ti,i+1 ) − (vi , ti,i+1 )(vi , vi ) = 0, then the fact that vi+1 is proportional to ti+1,i+1 implies (vi , vi+1 ) = 0. We proceed using induction on k to show that (vi , vi+k ) = 0 for i = 1, . . . , n − k when k = 1, . . . , n − 1. Suppose (vi , vi+j ) = 0 for all i = 1, . . . , n − j where j = 1, . . . , k and k < n − 1. By hypothesis ( ) (vi , ti+j+1,i+k+1 ) = vi+j , ti+j,i+k+1 − (vi+j , ti+j,i+k+1 )vi+j = (vi , ti+j,i+k+1 ) − (vi , ti+j,i+k+1 )(vi , vi+j ) = (vi , ti+j,i+k+1 ) 3

Math/CS 466/666 Programming Project 2 holds for j = 1, . . . , k where i = 1, . . . , n − k − 1. Therefore (vi , ti+k+1,i+k+1 ) = (vi , ti+1,i+k+1 ). Since

( ) (vi , ti+1,i+k+1 ) = vi , ti,i+k+1 − (vi , ti,i+k+1 )vi = (vi , ti,i+k+1 ) − (vi , ti,i+k+1 )(vi , vi ) = 0

the fact that ti+k+1,i+k+1 is proportional to vi+k+1 implies (vi+k+1 , vi ) = 0, which completes the induction. In particular, the vi ’s are orthonormal. Since t1,j = uj , then (v1 , t1,j ) = (v1 , uj ) for j = 1, . . . , n. Now, if 1 < k ≤ i ≤ j then ( ) (vi , tk,j ) = vi , tk−1,j − (vk−1 , tk−1,j )vk−1 = (vi , tk−1,j ) − (vk−1 , tk−1,j )(vi , vk−1 ) = (vi , tk−1,j ). Proceeding inductively obtains that (vi , ti,j ) = (vi , t1,j ) = (vi , uj ) for every i ≤ j.

4

Math/CS 466/666 Programming Project 2 (ii) Prove the matrix Q defined above is orthogonal and that A = QR. We have already shown that that vi is an orthonormal basis of Rn and that Q is an orthogonal matrix. In fact, since (vi , ti,j ) = (vi , uj ), then the modified Gram–Schmidt algorithm is mathematically identical to the original Gram–Schmidt algorithm. Therefore the Q obtained by the modified algorithm is the the same as for the original algorithm. We also have ∥ti,i ∥ = (vi , ti,i ) = (vi , uj ). Thus, { Ri,j =

(vi , uj ) for i ≤ j 0 for i > j.

From (∗) it follows that R = Qt A. Finally, the fact that QQt = I implies QR = Q(Qt A) = (QQt )A = IA = A. (iii) Let H be the Hilbert matrix and write a program that computes Q and R such that H = QR using the modified Gram–Schmidt algorithm described above. We begin with a simple matrix library similar to the one created in the computing lab. Note that some functions have been renamed and other have been added. The library header file is 1 2

/* p2lib.h--Simple Matrix Libraries for Programming Project 2 Written November 15, 2016 by Eric Olson for Math/CS 466/666 */

3 4 5

#ifndef P2LIB_H #define P2LIB_H

6 7 8 9 10 11 12 13 14 15 16 17

extern extern extern extern extern extern extern extern extern extern extern

void matprintA(int n,double A[n][n]); void matprintAt(int n,double A[n][n]); double vecXtY(int n,double X[n],double Y[n]); void matCisABt(int n,double C[n][n],double A[n][n],double B[n][n]); void matCisAtB(int n,double C[n][n],double A[n][n],double B[n][n]); void matCisAB(int n,double C[n][n],double A[n][n],double B[n][n]); double frobI(int n,double A[n][n]); double frobD(int n,double A[n][n],double B[n][n]); void matset(int n,double A[n][n],double B[n][n]); void mattrans(int n,double A[n][n]); void matRQisA(int n,double R[n][n],double Q[n][n],double A[n][n]);

18 19

#endif

while the code for the library is 1 2

/* p2lib.c--Simple Matrix Libraries for Programming Project 2 Written November 15, 2016 by Eric Olson for Math/CS 466/666 */

3

5

Math/CS 466/666 Programming Project 2 4 5 6 7

#include #include #include #include

8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

void matprintA(int n,double A[n][n]) { for(int i=0;i

6

Math/CS 466/666 Programming Project 2 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75

double r=0; for(int i=0;i

76 77

/*

Since the C Programming Language stores matrices in row-major order, we switch the order of the indices and work with tranposed matrices in the code. In this way the column operations of the Gram–Schmidt algorithm become row operations for the transposed matrices. Thus, the code below actually performs the factorization Qt Rt = H t or equivalently RQ = H. */

83 84 85 86 87 88 89 90 91 92

for(int i=0;i

7

Math/CS 466/666 Programming Project 2 } r=sqrt(vecXtY(n,Q[i],Q[i])); R[i][i]=r; for(int j=0;j

93 94 95 96

}

97 98

}

The main program now consists of a routine to generate the Hilbert matrix and perform the QR decomposition for the desired values of n. 1 2

/* p2iii.c--Compute QR decomposition of the Hilbert Matrix Written November 15, 2016 by Eric Olson for Math/CS 466/666 */

3 4 5 6 7

#include #include #include #include

8 9 10 11 12 13

void mathilb(int n,double H[n][n]) { for(int i=0;i

14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

#define N 4 double A[N][N],Q[N][N],R[N][N]; int main(){ printf("#p2iii\n"); mathilb(N,A); mattrans(N,A); matRQisA(N,R,Q,A); printf("Q=[\n"); matprintAt(N,Q); printf("]\nR=[\n"); matprintAt(N,R); printf("]\n"); return 0; }

8

Math/CS 466/666 Programming Project 2 (iv) Use the above program to find an orthogonal matrix Q and an upper triangular matrix R such that H = QR when n = 4. The output of the program is #p2iii Q=[ 8.3811635492e-01, -5.2264837396e-01, 1.5397276152e-01, -2.6306682088e-02; 4.1905817746e-01, 4.4171332392e-01, -7.2775380737e-01, 3.1568018506e-01; 2.7937211831e-01, 5.2882138625e-01, 1.3950552218e-01, -7.8920046265e-01; 2.0952908873e-01, 5.0207166632e-01, 6.5360920576e-01, 5.2613364176e-01; ] R=[ 1.1931517553e+00, 6.7049308394e-01, 4.7493260112e-01, 3.6983547090e-01; 0.0000000000e+00, 1.1853326749e-01, 1.2565509463e-01, 1.1754199276e-01; 0.0000000000e+00, 0.0000000000e+00, 6.2217740601e-03, 9.5660929494e-03; 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 1.8790487206e-04; ]

9

Math/CS 466/666 Programming Project 2 (v) For the matrices Q and R found using the above program compute ∥H − QR∥F ,

∥Qt Q − I∥F

and

∥QQt − I∥F

when n = 4, 6, 8, 10, 12 and 20. Slight changes to the code from the previous section resulted in the program 1 2

/* p2iv.c--Compute QR decomposition of the Hilbert Matrix Written November 15, 2016 by Eric Olson for Math/CS 466/666 */

3 4 5 6 7

#include #include #include #include

8 9 10 11 12 13

void mathilb(int n,double H[n][n]) { for(int i=0;i

14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

int main(){ printf("#p2iv\n"); int nvals[]={4, 6, 8, 10, 12, 20}; printf("#%3s %22s %22s %22s\n","n","|QtQ-I|_F","|QQt-I|_F","|H-QR|_F"); for(int i=0;i

10

Math/CS 466/666 Programming Project 2 The resulting output was #p2iv # n 4 6 8 10 12 20

|QtQ-I|_F 4.05023225521046e-13 2.72353370005502e-10 7.33787259380259e-07 2.55016714252341e-04 4.38888436635874e-01 2.37604463513680e+00

|QQt-I|_F 4.05051714551953e-13 2.72353367319302e-10 7.33787259425052e-07 2.55016714252345e-04 4.38888436635874e-01 2.37604463513680e+00

11

|H-QR|_F 2.77555756156289e-17 6.35960131078450e-17 1.11022302462516e-16 1.16936318589005e-16 1.08388987728284e-16 1.69079929654581e-16

Math/CS 466/666 Programming Project 2 3. The LAPACK routine DGEQRFP computes the QR factorization of a matrix. After factorization, the upper triangular part of the matrix A is overwritten with R. The matrix Q may be obtained from DORGQR routine. (i) Use these subroutines or equivalent ones to find an orthogonal matrix Q and an upper triangular matrix R such that H = QR when n = 4. The program 1 2

/* p3i.c--Compute QR decomposition of the Hilbert Matrix using LAPACK Written November 15, 2016 by Eric Olson for Math/CS 466/666 */

3 4 5 6 7 8

#include #include #include #include #include

9 10 11 12 13 14

void mathilb(int n,double H[n][n]) { for(int i=0;i

15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

#define N 4 double Q[N][N],R[N][N],TAU[N]; int main(){ printf("#p3i\n"); mathilb(N,Q); mattrans(N,Q); LAPACKE_dgeqrfp(LAPACK_COL_MAJOR,N,N,Q[0],N,TAU); for(int i=0;i

12

Math/CS 466/666 Programming Project 2 produces the output #p3i Q=[ 8.3811635492e-01, -5.2264837396e-01, 1.5397276152e-01, -2.6306682088e-02; 4.1905817746e-01, 4.4171332392e-01, -7.2775380737e-01, 3.1568018506e-01; 2.7937211831e-01, 5.2882138625e-01, 1.3950552218e-01, -7.8920046265e-01; 2.0952908873e-01, 5.0207166632e-01, 6.5360920576e-01, 5.2613364176e-01; ] R=[ 1.1931517553e+00, 6.7049308394e-01, 4.7493260112e-01, 3.6983547090e-01; 0.0000000000e+00, 1.1853326749e-01, 1.2565509463e-01, 1.1754199276e-01; 0.0000000000e+00, 0.0000000000e+00, 6.2217740601e-03, 9.5660929494e-03; 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 1.8790487206e-04; ]

(ii) Is the QR factorization found using LAPACK the same as the factorization found using the modified Gram–Schmidt algorithm in the previous question? The output appears to be the same.

13

Math/CS 466/666 Programming Project 2 (iii) For the matrices Q and R found using LAPACK compute ∥H − QR∥F ,

∥Qt Q − I∥F

and

∥QQt − I∥F

when n = 4, 6, 8, 10, 12 and 20. The program 1 2

/* p3iii.c--Compute QR decomposition of the Hilbert Matrix using LAPACK Written November 15, 2016 by Eric Olson for Math/CS 466/666 */

3 4 5 6 7 8 9

#include #include #include #include #include #include

10 11 12 13 14 15

void mathilb(int n,double H[n][n]) { for(int i=0;i

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

int main(){ printf("#p3iii\n"); int nvals[]={4, 6, 8, 10, 12, 20}; printf("#%3s %22s %22s %22s\n","n","|QtQ-I|_F","|QQt-I|_F","|H-QR|_F"); for(int i=0;i

14

Math/CS 466/666 Programming Project 2 matCisAB(N,T,Q,R); double d3=frobD(N,T,A); printf("%4d %22.14e %22.14e %22.14e\n",N,d1,d2,d3);

39 40 41

} return 0;

42 43 44

}

produces the output #p3iii # n 4 6 8 10 12 20

|QtQ-I|_F 1.37677167724569e-15 1.24564258906621e-15 9.30497652132227e-16 1.00435845246017e-15 2.65012986360284e-15 2.69041862515974e-15

|QQt-I|_F 1.29470816722684e-15 1.16912542563685e-15 8.82796670053343e-16 1.06611484391197e-15 2.55560000106286e-15 2.61831256177488e-15

|H-QR|_F 3.10316769155909e-16 1.86190061493545e-16 2.45915026490370e-16 4.53742119551398e-16 7.21311288194086e-16 4.10085107424621e-16

(iv) Compare the accuracy of the results obtained using LAPACK with the accuracy of the results obtained in using the modified Gram–Schmidt algorithm. The accuracy of the LAPACK routines are significantly better than the modified Gram– Schmidt algorithm.

15

Math/CS 466/666 Programming Project 2 (v) [Extra Credit and Math/CS 666] Look up the the algorithm used by LAPACK and explain how it works. What is special about the Hilbert matrix? Compare the speed and accuracy of the two methods when finding the QR factorization for random n × n matrices where n = 10, 50, 100, 500, 1000 and 2000. The LAPACK routine is based on the method of Householder reflections. A description of this algorithm can be found in [1] Wikipedia, The Free Encyclopedia, “QR Decomposition,” accessed on November 19, 2016, http://en.wikipedia.org/wiki/QR decomposition. [2] First Steps in Numerical Analysis, R.J. Hosking, Hodder Arnold, 1978. The basic idea is to construct a sequence of n Householder matrices of the form Hk = I − 2wk wkt such that

∥w∥ = 1

where

α 0 H1 A = ...

∗

··· A1

∗

0 where α = −signum(A1,1 )∥u1 ∥ and u1 is the first column of A. Note that w1 ∈ Rn is a unit vector of length n with n degrees of freedom which needs to be chosen to create n − 1 zeros upon multiplication. The sign of α is chosen by multiplying w1 by ±1 as needed. Note that the choice of sign for α increases the stability of the algorithm by maximizing the entry in the upper-left corner of the matrix H1 . This procedure is then repeated inductively to find a vector wk+1 ∈ {0}k × Rn−k such that the the last n − k components of wk+1 denoted by w ˜k+1 ∈ Rn−k satisfy ∥w ˜k+1 ∥ = 1 t and such that (I − 2w ˜k+1 w ˜k+1 )Ak is an (n − k) × (n − k) matrix with zeros in all but the top corner of the first column. The QR decomposition may then be obtained by taking Q = Hn Hn−1 · · · H2 H1 . Some special things about the Hilbert matrix were found in [3] Wikipedia, The Free Encyclopedia, “Hilbert matrix,” accessed on November 19, 2016, http://en.wikipedia.org/wiki/Hilbert matrix. In particular, the Hilbert matrix is known to be non-singular, symmetric, positive definite and to have an inverse matrix consisting only of integer entries. Moreover, the condition number of the n × n Hilbert matrix is such that √ (1 + 2)4n √ cond(H) ∼ K . n Thus, the Hilbert matrix is ill-conditioned for relatively small values of n. This is why the Gram–Schmidt process failed to find a useful QR decomposition for n = 20. The code 16

Math/CS 466/666 Programming Project 2 1 2

/* p4.c--Compute QR decomposition of a Random Matrix in two ways Written November 15, 2016 by Eric Olson for Math/CS 466/666 */

3 4 5 6 7 8 9 10 11

#include #include #include #include #include #include #include #include

12 13 14 15 16 17 18 19 20 21 22 23

static double tic_time; void tic() { struct timeval tp; gettimeofday(&tp,0); tic_time=tp.tv_sec+tp.tv_usec*1.e-6; } double toc() { struct timeval tp; gettimeofday(&tp,0); return tp.tv_sec+tp.tv_usec*1.e-6-tic_time; }

24 25 26 27 28 29 30

void matrand(int n,double A[n][n]) { srand(1234); for(int i=0;i

31 32 33 34 35 36 37 38 39 40 41 42 43 44

int main(){ printf("#p4\n"); setrlimit(RLIMIT_STACK,&(const struct rlimit){ RLIM_INFINITY,RLIM_INFINITY}); int nvals[]={10,50,100,500,1000,2000}; printf("#%3s %18s %18s %18s %18s\n", "n","GS sec","GS |QtQ-I|_F", "LAPACK sec","LAPACK |QtQ-I|_F"); for(int i=0;i

17

Math/CS 466/666 Programming Project 2 tic(); matRQisA(N,R,Q,A); double t1=toc(); matCisABt(N,T,Q,Q); double d1=frobI(N,T);

45 46 47 48 49

// QtQ since column major

50

tic(); memcpy(Q,A,sizeof(double)*N*N); LAPACKE_dgeqrfp(LAPACK_COL_MAJOR,N,N,Q[0],N,TAU); for(int i=0;i

51 52 53 54 55 56 57 58 59 60 61 62

} return 0;

63 64 65

}

produces the output #p4 # n 10 50 100 500 1000 2000

GS sec 5.9604644775e-06 2.6392936707e-04 1.7039775848e-03 2.9938697815e-01 2.2616238594e+00 1.8388961077e+01

GS |QtQ-I|_F 1.7127857513e-14 1.8314443720e-14 5.5393586389e-13 8.5169821271e-13 3.6279737321e-12 9.3367459484e-12

LAPACK sec 2.6202201843e-04 4.3702125549e-04 1.2831687927e-03 6.5408945084e-02 4.3951296806e-01 3.2365288734e+00

LAPACK |QtQ-I|_F 1.2667751215e-15 4.3156599481e-15 8.0321821104e-15 2.9375320935e-14 4.9525951276e-14 8.5636611542e-14

For reference, the above timings were performed using an AMD Athlon 64 X2 Dual Core 5400+ Processor with DDR2 800 RAM. While the Modified Gram–Schmidt Method is faster for n = 10 and 50, for larger values of n the LAPACK routines are faster. In particular, LAPACK is 5.7 times faster when finding the QR decomposition for n = 2000. Although the rounding error is not so severe when working with random matrices and both methods yielded usable QR decompositions for all values of n, LAPACK achieved greater accuracy in all cases.

18

Math/CS 466/666 Programming Project 2 QR Factorization Your work should be presented in the form of a typed report using clear and properly punctuated...