Advanced Introduction to C++, Scientific Computing and Machine Learning




Claudius Gros, WS 2019/20

Institut für theoretische Physik
Goethe-University Frankfurt a.M.

Finding Eigenvalues

diagonalizing matrices

$$ A\ \to\ D \qquad\qquad D\, =\, \begin{pmatrix} d_1 & 0 & \cdots & 0 \\ 0 & d_2 & \ddots & 0 \\ \vdots &\ddots&\ddots& 0 \\ 0 & \cdots & 0 &d_n \end{pmatrix} $$

QR decomposition

reflection matrices



Householder reflection



Householder reflection for column vectors

QR decomposition with Householder reflections

$$ Q_1 A = R_1 $$
$$ Q_2 R_1 = Q_2\big( Q_1 A\big) = R_2 \quad\qquad R_2 = \begin{pmatrix} * & * & * & * \\ 0 & * & * & * \\ 0 & 0 & * & * \\ 0 & 0 & * & * \end{pmatrix} $$

eigenvalues & unitary transformations

eigenvalues of a (given) matrix A

eigenvalues by QR iteration


series of unitary transformations

long-term behavior of maps

standard QR factorization

standalone QR code

Copy Copy to clipboad
Downlaod Download
// QR ITERATION #############################################################
/*
QR iteration for finding eigenvalues
LANGUAGE:     C++0x
PRECISION:    double
AUTHOR:       Hendrik Wernecke <wernecke@itp.uni-frankfurt.de>
VERSION DATE: October 19, 2016
LICENCE:      free to share, distribute and modify
INCLUDES:     standard libraries (and possible bugs)
PRELIMINARIES
        works so far only with square matrices
        only valid for matrices with real eigenvalues
DESCRIPTION
 By means of Householder decomposition a m*m matrix A is decomposed
   A = Q*R
 where Q is a m*m orthogonal matrix and R is a m*m upper triangular matrix.
 Using the iteration
   A_n-1 = Q_n*R_n
 and
   A_n = R_n*Q_n
 with initially
  A0 = A
 one can -- under certain circumstances -- converged in the limit
   lim_{n\to\infty}A_n = A'
 towards an upper triangular matrix A' which has the same eigenvalues
 than the original matrix A.

*/
//###########################################################################

/* ### INCLUDES ### */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>

/* ### FUNCTIONS forward definitions ### */
// definitions below MAIN

double** newMatrix(int m, int n);
void matrixCopy(double** mat0, double** mat1, int m, int n);
void UnitMatrix(double** mat, int m, int n);
void printMatrix(double** mat, int m, int n);
void printArray(double* arr, int m);
double norm(double* arr, int m);
void getColumn(double** mat, int k, double* arr, int m, int n);
void multiplyMatrices(double** mat1, double** mat2, double** erg, int m, int n);
void updateAMatrix(double** Ah, double** Qk, int k, int m, int n);
void transposeMatrix(double** mat, double** tmat, int m, int n);
bool qrHouseholderDecomp(double** A0, double** Q, double** R, int m, int n, double tol);

/* ### MAIN ### */

int main(int argLength, char* argValues[])  
{
/* SETTINGS */
        int maxitnum = 1000;    // maximal number of iterations
        double tol = 1e-4;      // tolerance for QR decomposition
                                // and accuracy of eigenvalues

/*
// example 1) dimension 3*3
        const int m = 3, n = 3; // dimensions of matrix (n!=m not yet implemented)
        double A[m][n] = {      // matrix to be decomposed
                {12, -51,   4},
                { 6, 167, -68},
                {-4,  24, -41}
                };

*/
/*
// example 2) dimension 4*4
        const int m = 4, n = 4; // dimensions of matrix (n!=m not yet implemented)
        double A[m][n] = {      // matrix to be decomposed
                { 4,  1,  3, -2},
                { 1, -2,  4,  1},
                { 3,  4,  1,  2},
                {-2,  1,  2,  3}
                };
*/
// example 3) symmetric matrix 
int size = 8;
if (argLength==2)
  size = atof(argValues[1]);
const int m=size, n=size;
double A[m][n];
for (int i=0; i<m; i++) {
        for (int j=i; j<n; j++) {
                A[i][j] = rand()%100;
                A[j][i] = A[i][j];
        }
}

        /* COMPUTATIONS */      

        // map A to real 2D array
        double** A0 = newMatrix(m, n);
        for (int i=0; i<m; i++) {
                for (int j=0; j<n; j++) {
                        A0[i][j] = A[i][j];
                }
        }
        // allocate matrices for storage of results
        double** Aprev = newMatrix(m, n);
        matrixCopy(A0, Aprev, m, n);
        double** Anext = newMatrix(m, n);

        double** Q = newMatrix(m, n);
        double** R = newMatrix(m, n);

        // some auxiliar variables
        int itnum = 0;
        bool check = false;

        // get started with the iteration
        while ( itnum<maxitnum )
        {
//printf("A_%i = \n", itnum); printMatrix(Aprev, m, n);
                if (not qrHouseholderDecomp(Aprev, Q, R, m, n, tol)) {
                        printf("QR decomposition failed!\n");
                        return 1;
                }
                multiplyMatrices(R, Q, Anext, m, n);
                // check if converged (comparing diagonal entries)
                double diff = 0;
                for (int i=0; i<m && i<n; i++) {
                                diff += fabs(Aprev[i][i]-Anext[i][i]);
                }
//printf("diff = %f\n", diff);
                if ( diff<=tol ) { check = true; break; }
                matrixCopy(Anext, Aprev, m, n);
                itnum++;
        }
        
        /* RESULTS */
        // announce final result
        if (check) {
                printf("eigenvalues converged within %i steps of iteration\n", itnum);
                printf("original matrix A_0 = \n"); printMatrix(A0, m, n);
                printf("iterated matrix A_%i = \n", itnum); printMatrix(Anext, m, n);
                printf("eigenvalues are: ");
                for (int i=0; i<m; i++) { printf(" %f, ", Anext[i][i]); } printf("\n");
        } else {
                printf("iteration did not converge in %i steps\nresulting in\n", itnum);
                printf("A_%i\n", itnum);
                printMatrix(Anext, m, n);
        }

        /* CLEANING UP */
        delete [] A0;
        delete [] Aprev;
        delete [] Anext;
        delete [] Q;
        delete [] R;

        return 0;
}
// END of MAIN


/* ### FUNCTIONS DEFINITIONS ### */

bool qrHouseholderDecomp(double** A0, double** Q, double** R, const int m, const int n, double tol)
// performs a QR decomposition base on Householder reflection
// see <https://en.wikipedia.org/wiki/QR_decomposition#Using_Householder_reflections>
// takes the matrix <A0> to be decomposed and the matrices <Q>, <R> (resulting composites)
// all matrices are <m>*<n> matrices
// for selfconsistent check A0=Q*R the tolerance <tol> is encountered
{
        // auxiliar matrix H
        double** H = newMatrix(m, n);
        UnitMatrix(H, m, n);

        // create local copy
        double** Ak = newMatrix(m, n);
        matrixCopy(A0, Ak, m, n);

        for (int k = 0; k<m-1 && k<n ; k++) {
//printf("before k = %i, Ah = \n", k); printMatrix(Ak, m, n);

                // get first column
                double x[m];
                getColumn(Ak, k, x, m, n);
//printf("x = \n"); printArray(x, m);

                // compute alpha as norm of first column
                double alpha = norm(x, m);
//printf("alpha = %f\n", alpha);

                double e[m];
                for (int i=0; i<m; i++) { e[i] = 0; } e[k] = 1;      // k-th unit vector
                double u[m];
                for (int i=0; i<m; i++) { u[i] = x[i] - alpha * e[i]; }
                double unorm = norm(u, m);
//printf("u = \n"); printArray(u, m);
                double v[m];
                for (int i=0; i<m; i++) { v[i] = u[i]/unorm; }
//printf("v = \n"); printArray(v, m);

                // assemble l-th Householder matrix Qk = 1-2 v*v^T
                double** Qk = newMatrix(m, n);
                for (int i=0; i<m; i++) {
                        for (int j=0; j<n; j++) {
                                if (i==j) {     // diagonal
                                        Qk[i][j] = 1 - 2*v[i]*v[j];
                                } else {        // non-diagonal
                                        Qk[i][j] = - 2*v[i]*v[j];
                                }
                        }
                }

//printf("Qk = \n"); printMatrix(Qk, m, n);

                // updating Ah matrix
                updateAMatrix(Ak, Qk, k, m, n);
//printf("A_k = \n"); printMatrix(Ak, m, n);

                // store product of Q_k (H=Q_n*Q_n-1*...*Q_1)
                double** Hh = newMatrix(m, n);
                matrixCopy(H, Hh, m, n);
                multiplyMatrices(Qk, H, Hh, m, n);
                H = newMatrix(m, n);
                matrixCopy(Hh, H, m, n);
//printf("H = \n"); printMatrix(H, m, n);
                delete [] Hh;
                delete [] Qk;
        }
        // Q = (Hn*...*H1)^T = H^T
        UnitMatrix(Q, m, n);
        transposeMatrix(H, Q, m, n);
//printf("Q = \n"); printMatrix(Q, m, n);

        // R = Q^T*A = H*A
        UnitMatrix(R, m, n);
        multiplyMatrices(H, A0, R, m, n);
//printf("R = \n"); printMatrix(R, m, n);

        // check that A = Q*R
        double** Ah = newMatrix(m, n);
        multiplyMatrices(Q, R, Ah, m, n);
//printf("A = \n"); printMatrix(Ah, m, n);

        // check deviation
        double dev = 0;
        for (int i=0; i<m; i++) {
                for (int j=0; j<n; j++) {
                        dev += fabs(Ah[i][j]-A0[i][j]);
                }
        }
        if (dev<=tol) {return true;}
        else {return false;}
}


double** newMatrix(int m, int n)
// returns pointer to new <m>*<n> matrix
{
        // allocate new matrix
        double** mat = (double**) malloc(m*sizeof(double));
        for (int i=0; i<m; i++) { mat[i] = (double*) malloc(n*sizeof(double)); }
        return mat;
}


void matrixCopy(double** mat0, double** mat1, int m, int n)
// copy content of matrix <mat0> to matrix <mat1>, both <m>*<n>
{
        // copy entries
        for (int i=0; i<m; i++) {
                for (int j=0; j<n; j++) {
                        mat1[i][j] = mat0[i][j];
                }
        }
        return;
}


void UnitMatrix(double** mat, int m, int n)
// fills the <m>*<n> matrix <mat> with a unity matrix
{
        // write entries
        for (int i=0; i<m; i++) {
                for (int j=0; j<n; j++) {
                        if (i==j) { mat[i][j] = 1; }
                        else { mat[i][j] = 0; }
                }
        }
        return;
}


void printMatrix(double** mat, int m, int n)
// prints the matrix <mat> with <m> rows and <n> columns
{
        printf("\n");
        for (int i=0; i<m; i++) {
                for (int j=0; j<n; j++) {
                        printf("  %8.3f", mat[i][j]);
                }
                printf("\n");
        }
        printf("\n");
        return;
}


void printArray(double* arr, int m)
// prints the array <arr> with length <m>
{
        printf("\n");
        for (int i=0; i<m; i++) {
                printf("  %8.3f", arr[i]);
        }
        printf("\n");
        return;
}


double norm(double* arr, int m)
// computes the norm of an array <arr> with length <m>
{
        double h = 0;
        for (int i=0; i<m; i++) { h += arr[i]*arr[i]; }
        return sqrt(h);
}


void getColumn(double** mat, int k, double* arr, int m, int n)
// extracts the <k>-th column of a <m>*<n> matrix <mat> and saves it to array <arr>
{
        if ( k<0 or k>=n ) { arr = NULL; return; }
        for (int i=0; i<m; i++) { arr[i] = mat[i][k]; }
        return;
}


void multiplyMatrices(double** mat1, double** mat2, double** erg, int m, int n)
// multiplies the matrices <mat1> and <mat2>
{
        for (int i=0; i<m; i++) { for (int j=0; j<n; j++) { erg[i][j] = 0; } }
        // compute product
        for (int i=0; i<m; i++) {
                for (int j=0; j<n; j++) {
                        for (int k=0; k<n; k++) {
                                erg[i][j] += mat1[i][k]*mat2[k][j];
                        }
                }
        }
        return;
}


void updateAMatrix(double** Ah, double** Qk, int k, int m, int n)
// updateds the new matrix <Ah> for next Householder step
// takes the Householder matrix <Qk>, both <m>*<n> matrices
// and the Householder iteration step number <k>
{
        // copy Qk*A
        double** QA = newMatrix(m, n);
        multiplyMatrices(Qk, Ah, QA, m, n);
        // create dummy
        double** H = newMatrix(m, n);

        /* assemble new matrix */

        // k+1 zero columns
        for (int i=0; i<k+1; i++) { for (int j=0; j<n; j++) { H[i][j] = 0; } }
        // k+1 zero rows
        for (int i=0; i<m; i++) { for (int j=0; j<k+1; j++) { H[i][j] = 0; } }
        // k+1 diagonal elements set to one
        for (int i=0; i<k+1; i++) { H[i][i] = 1; }
        // rest of the matrix filled with values from QA
        for (int i=k+1; i<m; i++) { for (int j=k+1; j<n; j++) { H[i][j] = QA[i][j]; } }

        // copy result to Ah
        for (int i=0; i<m; i++) { for (int j=0; j<n; j++) { Ah[i][j] = H[i][j]; } }

        // clean up
        delete [] QA;
        delete [] H;

        return;
}


void transposeMatrix(double** mat, double** tmat, int m, int n)
// transposes the <m>*<n> matrix <mat>
// and stores the transpose in the <n>*<m> matrix  <tmat>
{
        // transpose    
        for (int j=0; j<n; j++) {
                for (int i=0; i<m; i++) {
                        tmat[j][i] = mat[i][j];
                }
        }
        return;
}

// ### END OF FILE ###