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

Claudius Gros, SS 2023

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

# 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}$$
• finding the eigenvealues $\ d_n\$ numerically

### QR decomposition

• series of successive transformations
here:: $\ QR\$ iteration
$R\$ (upper right) triangular matrix
$Q\$ orthogonal matrix $$\quad Q Q^T\,=\,I\qquad\qquad Q^{-1}\,=\,Q^T$$
• determine the transformation matrix $\ Q$
here:: using Householder transformations

# reflection matrices

• plane defined by $\quad \mathbf{x}\cdot\mathbf{v}=0$
• reflection matrix
$$Q \equiv I-2\mathbf{v}\mathbf{v}^T \qquad\quad Q_{ij} = \delta_{ij} -2v_iv_j$$
with
$$\mathbf{v}\cdot\mathbf{v}=1 \quad\quad \mathbf{x}=\mathbf{x}_0 + (\mathbf{x}\cdot\mathbf{v}) \mathbf{v} \quad\quad \mathbf{x}_0 \cdot\mathbf{v}=0$$
we have $$\begin{array}{rcl} Q_{ij}(x_0)_j &=& \big(\delta_{ij} -2v_iv_j\big)(x_0)_j \ =\ (x_0)_i\\[0.5ex] Q_{ij}v_j &=& \big(\delta_{ij} -2v_iv_j\big)v_j \ =\ v_i-2v_i \ =\ -v_i \end{array}$$
$$Q\,\mathbf{x} = \mathbf{x}_0-(\mathbf{x}_0\cdot\mathbf{v})\mathbf{v}$$
• orthogonality for $\ Q^T = Q\$ (symmetric) $$\begin{array}{rcl} \big(QQ^T\big)_{ik} &=& \sum_j\big(\delta_{ij}-2v_iv_j\big) \big(\delta_{jk}-2v_jv_k\big) \\ &=& \delta_{ik}-2v_iv_k-2v_iv_k +4 v_i(\mathbf{v}\cdot\mathbf{v})v_k \\[0.5ex] &=& \delta_{ik} \ \equiv\ (I)_{ik} \end{array}$$

# Householder reflection

• given a direction $$\mathbf{a}\,=\,|\mathbf{a}|\mathbf{e} \qquad\quad |\mathbf{e}|=1 \qquad\quad \qquad\quad \qquad\quad$$ then $$\mathbf{u} = \mathbf{x}-|\mathbf{x}|\mathbf{e} \qquad\quad \mathbf{v} = \mathbf{u}/|\mathbf{u}| \qquad\quad \qquad\quad \qquad\quad$$ defines a rotation matrix $$H_{\mathbf{e}}\,=\,Q_{\mathbf{v}}\,=\, I-2\mathbf{v}\mathbf{v}^T \qquad\quad \qquad\quad \qquad\quad$$ that rotates $\ \mathbf{x}\$ to $\ \mathbf{e}$
• Householder matrix $\ H_{\mathbf{e}}$

# Householder reflection for column vectors

• column vectors $$A = \begin{pmatrix} a_{11} & a_{12} & \cdots \\ a_{21} & a_{22} & \cdots \\ \vdots & \vdots & \ddots \\ \end{pmatrix} \equiv \big(\mathbf{a}_1,\mathbf{a}_2,\dots\big) \quad\qquad \mathbf{a}_n = \begin{pmatrix} a_{1n} \\ a_{2n} \\ \vdots \end{pmatrix}$$
• first step: zero first column $$Q_1 A = R_1 \quad\qquad R_1 = \begin{pmatrix} * & * & * & * \\ 0 & * & * & * \\ 0 & * & * & * \\ 0 & * & * & * \end{pmatrix} \equiv \big(Q_1\mathbf{a}_1,\, Q_1\mathbf{a}_2,\, Q_1\mathbf{a}_3,\, Q_1\mathbf{a}_4\big)$$ by rotating $\ \mathbf{a}_1\$ to $\ \mathbf{e}_1\$ (direction) $$\fbox{\displaystyle\phantom{\large|} \mathbf{a}_1 \ \to\ \mathbf{e}_1 \phantom{\large|}}$$
• many transformations possible
here Hauseholder rotation $\quad \mathbf{a}_1\,\hat{=}\,\mathbf{x} \quad \mathbf{e}_1\,\hat{=}\,\mathbf{e}$

$$Q_1 = I-2\mathbf{v}_1\mathbf{v}_1^T \qquad\quad \mathbf{v}_1 = \mathbf{u}_1/|\mathbf{u}_1| \qquad\quad \mathbf{u}_1 = \mathbf{a}_1-|\mathbf{a}_1|\mathbf{e}_1$$

# QR decomposition with Householder reflections

$$Q_1 A = R_1$$
• iteration
$$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}$$
• repeat first step for $\ R_1'$, with $$R_1 = \begin{pmatrix} * & * & * & * \\ 0 & & & \\ \vdots & & R_1' & \\ 0 & & & \end{pmatrix} \qquad\quad R_1' = \begin{pmatrix} * & * & * \\ * & * & * \\ * & * & * \end{pmatrix}$$ and $$Q_2 = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & & & \\ \vdots & & Q_2' & \\ 0 & & & \end{pmatrix} \qquad\quad Q_2' = \begin{pmatrix} * & * & * \\ * & * & * \\ * & * & * \end{pmatrix}$$ obtaining
$$R = Q_n\cdots Q_2Q_1A \quad\qquad Q R =A \quad\qquad Q = Q_1^T\cdots Q_n^T$$
• numerical effort:: $\ O(n^3)$

# eigenvalues & unitary transformations

• starting with $$A = QR \qquad\quad R = Q^{-1}A = Q^TA$$ we find
$$B \ \equiv\ RQ \ =\ Q^TAQ$$
• eigenvalues are invariant under unitary transformations $$\begin{array}{rcl} \lambda\mathbf{x} &=& A\mathbf{x} \ =\ QQ^T AQQ^T\mathbf{x}\\[0.5ex] \lambda\big(Q^T\mathbf{x}\big) & =& Q^TAQ\big(Q^T\mathbf{x}\big) \ =\ B\big(Q^T\mathbf{x}\big) \end{array}$$ $B\$ has hence an eigenvalue $\ \lambda\$ if and only if $\ A\$ has.

### eigenvalues of a (given) matrix A

• the eigenvalues $\ \lambda_n\$ of an R (upper triangular) matrix are on the diagonal $$R = \begin{pmatrix} \lambda_1 & * & * \\ 0 & \lambda_2 & * \\ 0 & 0 & \lambda_3 \end{pmatrix} \qquad\quad R-\lambda I = \begin{pmatrix} \lambda_1-\lambda & * & * \\ 0 & \lambda_2-\lambda & * \\ 0 & 0 & \lambda_3-\lambda \end{pmatrix}$$ having the characteristic polynomial $\ det(R-\lambda I)=\prod_n (\lambda_n-\lambda)\$
• $B\$ is however -not- diagonal!
• task: find a orthogonal matrix $\ Q^*$
such that $\ B^*\to R^*\$ becomes a upper right matrix $\ R^*$ $$B^*\quad \to\quad \fbox{\displaystyle\phantom{\large|} R^* \ =\ (Q^*)^TAQ^* \phantom{\large|}}$$
• iterative approximation

# eigenvalues by QR iteration

### series of unitary transformations

• select algorithm for a QR decomposition for the mapping $$A\quad \to\quad B\,=\,Q^TAQ \qquad\qquad (A\,=\,QR)$$ here:: Householder
task:: $B$ upper right matrix
• now consider a map in the space of orthogonal matrices defined by $$Q_k\ \to\ Q_{k+1} \qquad\quad \fbox{\displaystyle\phantom{\large|} Q_{k+1}R_{k+1}\ \equiv\ AQ_{k} \phantom{\large|}}$$
• starting $\ Q_0=I$
$$Q_{1}R_{1}\,=\,A, \qquad\quad Q_{2}R_{2}\,=\, AQ_{1}, \qquad\quad \ldots$$
• question: behavior of map for $\ k\to\infty$

# long-term behavior of maps

• dynamical system theory:
maps can converge (here:: in the phase space of othogonal matrices) to
• fixpoints
• limit cycles
• chaotic atractors
see complex-systems lecture course
• convergence (to a fixpoint) $$\lim_{k\to\infty} Q_k\ =\ Q^* \qquad\quad \lim_{k\to\infty} R_k\ =\ R^* \qquad\quad \fbox{\displaystyle\phantom{\large|} Q^*R^*\ \equiv\ AQ^* \phantom{\large|}}$$ leads to
$$R^* \ =\ (Q^*)^TAQ^*$$
• $Q^*\$ orthononal
:: eigenvalues of $\ A\$ and $\ R^*\$ identical
• $R^*\$ upper right matrix
:: eigenvalues of $\ R^*\$ on diagonal
$\color{red}\Longrightarrow\quad$ eigenvalue problem solved

# standard QR factorization

• avoiding confusion: $\ W_k$ for WR (QR) decompostion
• starting $\ A_0=A,\qquad A_0=W_0R_0,\quad A_1=R_0W_0$
• map (unitary transformation) $$A_k\quad \to\quad A_{k+1}\ \equiv\ R_kW_k \ =\ W_k^TA_k W_k \qquad\qquad A_k \ =\ W_kR_k$$ viz $$A_{k+1} \ =\ W_k^T\cdots W_0^T A_0 W_0\cdots W_k$$
• reformulation $$\tilde Q_k \equiv W_0W_1\cdots W_k \qquad\qquad \tilde Q_k A_{k+1} \ =\ A \tilde Q_k$$
• equivalence of fixpoints $$\fbox{\displaystyle\phantom{\large|} Q_{k+1}R_{k+1}\ \equiv\ AQ_{k} \phantom{\large|}} \qquad\qquad \fbox{\displaystyle\phantom{\large|} \tilde Q_k A_{k+1} \ =\ A \tilde Q_k \phantom{\large|}}$$ when $$\lim_{k\to\infty} A_{k} \ \to\ \tilde R^* \qquad\quad \lim_{k\to\infty} \tilde Q_{k} \ \to\ \tilde Q^* \qquad\quad \tilde R^* \ =\ (\tilde Q^*)^TA\tilde Q^*$$ converging to a upper triangular matrix
• care when
• close-by eigenvalues
• care when complex eigenvalues
:: all matrices are real
• the eigenvectors $$\begin{array}{rcl} \lambda\mathbf{x}_\lambda^* & =& \tilde R^*\mathbf{x}_\lambda^* \ = \ (\tilde Q^*)^TA\tilde Q^* \mathbf{x}_\lambda^* \\[0.5ex] \lambda\tilde Q^*\mathbf{x}_\lambda^* & =& A\tilde Q^* \mathbf{x}_\lambda^* \end{array}$$ of $\ A$ are $\quad \mathbf{x}_\lambda = \tilde Q^* \mathbf{x}_\lambda^*$

# standalone QR code

// 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 ###