Homework "Linear Equations"

Implement functions to i) solve linear equations; ii) calculate matrix inverse; and iii) calculate matrix determinant. The method to use is the modified Gram-Schmidt orthogonalization.

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 vectors only at the end of the process. This makes the Gram–Schmidt process applicable for iterative methods like the Arnoldir/Lanczos iteration. Also its ease of implementation makes it a useful exercise for the students.

• To work with matrices and vectors you may use the matrix and vector classes from "../matlib/matrix". They implement some basic matrix/vector operations. Alternatively, you can implement your own matrix/vector classes the way you like it (a good way to learn a language). 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(){ ... }
     } 

  3. You can add matrices Q and R and the methods decomp, solve, inverse, det directly to your matrix class (declare it partial then you can add things to it using a file in your homework directory).

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