Problem "Linear Equations"

Objective

Implement functions for solving linear equations, calculating matrix inverse and matrix determinant.

Introduction

To work with matrices you can use the matlib/matrix/matrix.cs[→] that implements some basic matrix operations. Alternatively, you can creaite your own matrix class the way you like it. At the bare minimum your class might look something like (read the manual about indexers in Csharp),

public class matrix{
	public readonly int size1,size2;
	private double[] data; // keep matrix elements in one big array
	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];      // getter: for example, aij=A[i,j];
		set => data[i+j*size1]=value;// setter: for example, A[i,j]=5;
		}
}

Tasks

1. (6 points) Solving linear equations using QR-decomposition by 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 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.

   1. Implement a function, say

      public static void QRGSdecomp(matrix A, matrix R)
      

      that performs in-place modified Gram-Schmidt orthogonalization of an n×m (where n≥m) matrix A: on exit matrix A is replaced with Q and the m×m matrix R holds the computed matrix R.

      Alternatively, you can make a class

      public class QRGS{
      	public matrix Q,R;
      	public QRGS(matrix A){/* run the Gram-Schmidt process */}
      	public vector solve(vector b){/* back-substitute Q^T*b */}
      	public matrix inverse(){/* return the inverse matrix (part B) */}
      	public double determinant(){/* return ΠiRii */}
      }
      

      Check that you function 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;

   2. Implement a function, say

      public static vector QRGSsolve(matrix Q, matrix R, vector b)
      

      or, if you like classes, implement the corresponding method

      public vector solve(vector b)
      

      that—given the matrices Q and R from QRGSdecomp—solves the equation QRx=b by applying Q^T to the vector b (saving the result in a new vector x) and then performing in-place back-substitution on x and returning it.

      Check that you function 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

   Implement a function, say

   public static matrix QRGSinverse(matrix Q, matrix R)
   

   or a method,

   public matrix inverse()
   

   that—given the matrices Q and R from QRGSdecomp—calculates the inverse of the matrix A into a new matrix 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 pyxplot/gnuplot and fit with N^3.