Homework "Linear Equations"

Implement functions to solve linear equations; calculate matrix inverse; and calculate matrix determinant using the (modified) Gram-Schmidt orthogonalization method.

The Gram-Schmidt orthogonalization process, even modified, is less stable and accurate than the Givens roation algorithm. On the other hand, the Gram-Schmidt process produces the j-th orthogonalized vector after the j-th iteration, while orthogonalization using Givens rotations produces all the vectors only at the end. This makes the Gram–Schmidt process applicable for iterative methods like the Arnoldi iteration. Also its ease of implementation makes it a useful exercise for the students.

• To work with matrices and vectors you may well use the matrix and vector classes from "../matlib/matrix". They implement some basic matrix/vector operations. Alternatively, you can creaite your own matrix/vector classes the way you like it. At the bare minimum your classes might look something like (read the manual about indexers in Csharp),

  public class vector{
     public double[] data;
     public int size => data.Length;
     public double this[int i]{   // indexer
        get => data[i];        // getter
        set => data[i]=value;  // setter
        }
     public vector(int n){        // constructor
        data=new double[n];
        }
  }
  	

  public class matrix{
  	public readonly int size1,size2;
  	private double[] data;  // to keep matrix elements
  	public matrix(int n,int m){      // constructor
  		size1=n; size2=m;
  		data = new double[size1*size2];
  		}
  	public double this[int i,int j]{     // indexer
  		get => data[i+j*size1];
  		set => data[i+j*size1]=value;
  		}
  }
  

• You can generate (pseudo-)random numbers, uniformly distributed from zero to one, like this,

  var rnd = new System.Random(1); /* or any other seed */
  double x = rnd.NextDouble();
  double y = rnd.NextDouble();
  double z = rnd.NextDouble();
  	

• QR-factorization of a matrix consists of two matrices, Q and R. There are several options for the interface to your routines:

  1. Return Q and R in a tuple (like scipy),

     public static class QRGS{
        public static (matrix,matrix) decomp(matrix A){
           matrix Q=A.copy();
           matrix R=new matrix(A.size2,A.size2);
           /* orthogonalize Q and fill-in R */
           return (Q,R);
           }
        public static vector solve(matrix Q, matrix R, vector b){ ... }
        public static double det(matrix R){ ... }
        public static matrix inverse(matrix Q,matrix R){ ... }
     }
     	

  2. Keep Q and R and an instance of the non-static class QRGS,

     public class QRGS{
        matrix Q,R;
        public QRGS(matrix A){ /* the above "decomp" is the constructor here */
           Q=A.copy();
           R=new matrix(A.size2,A.size2);
           /* orthogonalize Q and fill-in R */
           }
        public vector solve(vector b){ ... }
        public double det(){ ... }
        public matrix inverse(){ ... }
     }
     	

1. (6 points) Solving linear equations using QR-decomposition by modified Gram-Schmidt orthogonalization

   Implement a static (or, at your choice, non-static) class "QRGS" with functions "decomp", "solve", and "det" (as indicated above). In the non-static class "decomp" becomes a constructor and must be called "QRGS").

   The function "decomp" (or the constructor QRGS) should perform stabilized Gram-Schmidt orthogonalization of an n×m (where n≥m) matrix A and calculate R.

   The function/method "solve" should use the matrices Q and R from "decomp" and solve the equation QRx=b for the given right-hand-side "b".

   The function/method "det" should return the determinant of matrix R which is equal to the determinant of the original matrix. Determinant of a triangular matrix is given by the product of its diagonal elements.

   Check that your "decomp" works as intended:

   • generate a random tall (n>m) matrix A (of a modest size);
   • factorize it into QR;
   • check that R is upper triangular;
   • check that Q^TQ=1;
   • check that QR=A;

   Check that you "solve" works as intended:

   • generate a random square matrix A (of a modest size);
   • generate a random vector b (of the same size);
   • factorize A into QR;
   • solve QRx=b;
   • check that Ax=b;

2. (3 points) Matrix inverse by Gram-Schmidt QR factorization

   Add the function/method "inverse" to your "QRGS" class that, using the calculated Q and R, should calculate the inverse of the matrix A and returns it.

   Check that you function works as intended:

   • generate a random square matrix A (of a modest size);
   • factorize A into QR;
   • calculate the inverse B;
   • check that AB=I, where I is the identity matrix;

3. (1 point) Operations count for QR-decomposition

   Measure the time it takes to QR-factorize (with your implementation) a random NxN matrix as function of N. Convince yourself that it goes like O(N³): measure (using the POSIX time utility) the time it takes to QR-factorize an N×N matrix for several values of N, plot the time as function of N in gnuplot and fit with N³.

   You can build your timing data like this (you might want to read

   man bash | less +/"for name"
   man 1 seq
   man 1 time 

   ),

   out.times.data : main.exe
   	>$@
   	for N in $$(seq 100 20 200); do \
   		time --format "$$N %e" --output $@ --append \
   		mono $< -size:$$N 1>out 2>err ;\
   	done