AMS 526 Final Project

1. AMS 526 Final Project By Bernard Moore and Xiaolong Huang

2. Problem (Part 1) • Write a C, C++, or Fortran program to generate a 200x200 strictly diagonally dominant symmetric matrix A using a random number generator. • Make a Cholesky decomposition of A. • Solve Ax=b, for a randomly generated vector b. • Find the inverse of A, via AA-1= I. • Results on the condition number of A.

3. Definition • Let A be a real n x nmatrix. A is said to be strictly (row) diagonally dominant provided

4. C Programming Option • We wrote a C program to generate a real symmetric strictly diagonally dominant matrix A by using the random number generator rand(). • We design the program so that A will have entries strictly between -1 and 1.

5. The Algorithm (Step 1) double R[N][N],b[N][1]; for(i=0;i<N;i++) for(j=0;j<N;j++) vvvR[i][j]=2*(double)rand()/(double)RAND_MAX-1 We begin by defining a matrix with random entries strictly between -1 and 1.

6. The Algorithm (Step 2) for(i=0;i<N;i++) S=0; for(j=0;j<N;j++) if(i!=j) S=S+fabs(R[i][j]); if(S>1) S=S/(N-1); for(k=0;k<N;k++) if(k!=i) R[i][k]=R[i][k]/(N-1) R[i][i]=(double)rand()*(1-S)/(double)RAND_MAX+S; Next we must guarantee that for each row of the matrix the row sum S (of the off-diagonal entries) is less than the diagonal entry in that row. We define the diagonal entries as random number between S and 1.

7. The Algorithm (Step 3) for(i=0;i<N;i++) for(j=i+1;j<N;j++) if(fabs(R[i][j])<fabs(R[j][i])) R[j][i]=R[i][j]; else R[i][j]=R[j][i]; This step guarantees that the matrix will be symmetric and strictly diagonally dominant by assingning the R[i][j]=min(R[i][j],R[j][i]), for all i and j.

8. Finding the Cholesky Decomposition double *cholesky(double *A, int n) { double *L = (double*)calloc(n * n, sizeof(double)); inti,j,k; for (i = 0; i < n; i++) for (j = 0; j < (i+1); j++) double s = 0; for (k = 0; k < j; k++) s += L[i * n + k] * L[j * n + k]; L[i * n + j] = (i == j) ? sqrt(A[i * n + i] - s) : (1.0 / L[j * n + j] * (A[i * n + j] - s)); return L; } • To find the Cholesky decomposition of A, we wrote the following function which returns the pointer *L, giving the entries of the lower triangular matrix L, where A=LLT.

9. Finding the Inverse of A (Step 1) double *syminv(double *R, int n) double *RT = (double*)calloc(n * n, sizeof(double)); inti,j,k; for(i=0;i<n;i++) for(j=0;j<n;j++) RT[j*n+i]=R[i*n+j]; Now we find the inverse of A by solving AA-1=I. To solve this equation we use the Cholesky decomposition of A, and solve LLTxi = eiand set LT xi = yi, then we must first solve the equation Lyi= eiand then solve LT xi = yi. We begin the algorithm by finding the transpose of the lower triangular matrix L from the Cholesky decomposition of A.

10. Finding the Inverse (Step 2) double *e = (double*)calloc(n*n, sizeof(double)); e[0]=1; for(i=0;i<n;i++) for(j=0;j<n;j++) if(i=j) e[n*i+j]=1; double *y = (double*)calloc(n*n, sizeof(double)); double s; for(k=0;k<n;k++) for(i=0;i<n;i++) s=0; for(j=0;j<i;j++) s=s+R[i*n+j]*y[j*n+k]; y[i*n+k]=(e[k*n+j]-s)/R[i+n*i]; Now we create the identity matrix (as a pointer), and solve the equation Lyi = eiby forward substitution.

11. Finding the Inverse (Step 3) double *x = (double*)calloc(n*n, sizeof(double)); double t; for(k=0;k<n;k++) for(i=n-1;i>=0;i--) t=0; for(j=n-1;j>i;j--) t=t+RT[i*n+j]*x[j*n+k]; x[i*n+k]=(y[n*i+k]-t)/RT[i+n*i]; return x; We now generate the solution to Ax=b, by solving the equation LTxi = yi by row-oriented backward substitution.

12. Calculating the Condition Number of A (Part 1) • First we calculate the norms of A. • As it happens, for a symmetric matrix A we have, |A|1 = |A|∞. • We begin by finding the ∞-norm of A, by first defining the maximum of two numbers, and running a for loop to find the maximum row sum. #define max(a,b) ((a) > (b) ? (a) : (b)) double r, maxr; maxr=0; for(k=0;k<N;k++) { r=0; for(j=0;j<N;j++) { r=r+fabs(R[k][j]); } maxr=max(r,maxr); }

13. Calculating the Condition Number of A (Part 2) • We now use the same algorithm to compute the norms of A-1. • Recall that the condition number with respect to the p-norm is given by kp (A) = |A|p|A-1|p • Symmetric diagonally dominant n x n matrices are well conditioned especially for large values of n. • In particular, for entries |aij|< 1, we have observed that for large values of n, • k1(A) = k∞(A) ≈ 4, • and • k2(A) ≈ 2.

14. Calculating the Condition Number of A (Part 3)

15. Results on the Condition Number of AChatterjee and Diaconis [1]

16. Results on the Condition Number of A Varah [2] Let Varah showed that if A is strictly diagonally dominant then the following inequality holds We have observed that as n approaches ∞, the value of L(A), which is bounded above by 1, tends to approach 0 making the upper bound on |A-1|∞ very large.

17. Results on the Condition Number of A • From our observations this inequality, like the one given in the Watkins [3] seems to only be useful for small values of n. • In our case where we have the entries|aij|<1, we encounter a more convergent nature of the condition number. We expect there is a bit nicer bound for the norms of A and A-1 for our specialized matrices, which remains to be seen, and is an interesting topic open for further investigation.

18. Problem (Part 2) For the 200x200 strictly diagonally dominant symmetric matrix A as in Part 1 • Perform similarity transformation to A and reduce it to a tridiagonal matrix. • Find the eigenvalues of A. • Analyze the efficiency of the algorithm used for the previous similarity transformation. • Prove that such a matrix is positive definite.

19. Householder Reduction • For each step k, use the column vector A(k+1:n, k) to generate the reflector. • After left and right multiplication Q’*A*Q, A(k+2:n, k) is eliminated • In the end, the lower half matrix below the first sub-diagonal is zero, which means the matrix is in Hessenberg form. • Total flops: (10/3) n3

20. Symmetric Householder Reduction • Simple version of algorithm • Only use half of the matrix entries • The matrix in the end becomes tri-diagonal • Total flops: (4/3) n3 • In programming, use one-dimensional array to record only the lower half of the matrix, and then perform Householder Reduction. • (A(i,j)=L(i*(i+1)/2+j))

21. QR algorithm (For tri-diagonal matrix) • Start with k=n • In each step • Use a 2x2 rotator to create the bulge (with proper shift) • Apply the rotator from left and right to chase the bulge • Iteration stops when |A(k-1,k)|<10^(-16) • Record A(k,k) as one eigenvalue • Move to A(1:k-1,1:k-1), perform the same algorithm • Stop when k=2. Then solve the remaining two eigenvalues from A(1:2,1:2)

22. Shift in QR algorithm • Wilkinson Shift: - one of the two eigenvalues of A(k-1:k,k-1:k) whichever is closer to A(k,k) - cubic convergence • Rayleigh Quotient Shift: A(k,k) - cubic convergence, with exceptions

23. Proposition

24. Proof of Proposition (Part 1)

25. Proof of Proposition (Part 2)

26. Proof of Proposition (Part 3)

27. Proof of Proposition (Part 4)

28. Proof of Proposition (Part 5)

29. Proof of Proposition (Part 6)

30. References [1] [2]