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

Claudius Gros, WS 2021/22

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

# linear equations

$$\left.\begin{array}{cccccc} 2x&-3 y& -z & + 2w& =&4\\ 4x&-4 y& -z & + 4w& =&4\\ 2x&-5 y& -3z& + 3w& =&9\\ \end{array}\right\} \qquad\quad A\,\vec x=\vec b \qquad\quad \vec x = \left(\begin{array}{c} x \\ y \\ z \\ w \end{array}\right)$$
• $3\times 4$ matrix $\ \ A$
• augmented matrix $\ \ (A|\vec b)$
• $$\left( \begin{array}{rrrr|r} 2&-3&-1& 2&4\\ 4&-4&-1& 4&4\\ 2&-5&-3& 3&9\\ \end{array} \right)$$
• solution via iterative elimination

# Gauss elimination

$$\left( \begin{array}{rrrr|r} 2&-3&-1& 2&4\\ 4&-4&-1& 4&4\\ 2&-5&-3& 3&9\\ \end{array} \right)$$
• solve for $x$ in first equation and substitute into the other equations
$2x=4+3y+z-2w$
• substract (multiples of) first row
$$\left( \begin{array}{rrrr|r} 2&-3&-1& 2&4\\ 0& 2& 1& 0&-4\\ 0&-2&-2& 1&5\\ \end{array} \right)$$
• solve for $y$ in second equation and substitute into third and forth equation
$2y=-4-z$
• substract (multiples of) second row
$$\left( \begin{array}{rrrr|r} 2&-3&-1& 2&4\\ 0& 2& 1& 0&-4\\ 0& 0&-1& 1&1\\ \end{array} \right)$$
• in matrix form
$$\tilde{A}\,\vec x = \tilde{b} \qquad\quad \left( \begin{array}{rrrr} 2&-3&-1& 2\\ 0& 2& 1& 0\\ 0& 0&-1& 1\\ \end{array} \right) \left(\begin{array}{c} x \\ y \\ z \\ w \end{array}\right) = \left(\begin{array}{c} 4 \\ -4 \\ 1 \end{array}\right)$$

# back substitution

$$\left( \begin{array}{rrrr} 2&-3&-1& 2\\ 0& 2& 1& 0\\ 0& 0&-1& 1\\ \end{array} \right) \left(\begin{array}{c} x \\ y \\ z \\ w \end{array}\right) = \left(\begin{array}{c} 4 \\ -4 \\ 1 \end{array}\right)$$
• last equation
$$-z+w=1$$
• take $w$ as free parameter and substitute back
$$\begin{array}{rcl} z&=&w-1\\ 2y+z=-4\qquad\quad y &=&\frac{-3-w}{2} \\ 2x-3y-z+2w=4\qquad\quad x &=&\frac{-3-5w}{4} \end{array}$$

with

$$2x = 4+3y+z-2w= 4-\frac{9}{2}-\frac{3}{2}w+w-1-2w =-\frac{3}{2} -\frac{5}{2}w$$
• final (one-dimensional) solution
$$\vec x= \left(\begin{array}{c} x \\ y \\ z \\ w \end{array}\right) = \left(\begin{array}{c} -(3+5w)/4 \\ -(3-w)/2 \\ w-1 \\ w \end{array}\right) \qquad\qquad \mbox{for any}\ \ w$$

# pivoting

• Gauss elimination
$$\left( \begin{array}{cccc} *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ 0&*&*&*\\ 0&*&*&*\\ 0&*&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ 0&*&*&*\\ 0&0&*&*\\ 0&0&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ 0&*&*&*\\ 0&0&*&*\\ 0&0&0&*\\ \end{array} \right)$$
• pivoting
select largest element $[*]$ in row
$$\left( \begin{array}{cccc} *&[*]&*&*\\ *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ *&0&*&[*]\\ *&0&*&*\\ *&0&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ *&0&*&*\\ *&0&[*]&0\\ *&0&*&0\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ *&0&*&*\\ *&0&*&0\\ *&0&0&0\\ \end{array} \right)$$
• one needs to divide by the pivot selected,
pivot need to be non-vanishing
• numerical stable if pivot is large
• partial pivoting
only largest element in given row/column
• complete pivoting
largest element in matrix

# swapping rows

### order of equations irrelevant

• swapping rows:: pivot of first column to first row $$\left( \begin{array}{rrr|r} 2&-3&-1& 4\\ \underline{4}&-4&-1& 4\\ 2&-5&-3& 9\\ \end{array} \right) \qquad \equiv \qquad \left( \begin{array}{rrr|r} \underline{4}&-4&-1& 4\\ 2&-3&-1& 4\\ 2&-5&-3& 9\\ \end{array} \right)$$
• Gauss + pivot of second column to second row $$\left( \begin{array}{rrr|r} \underline{4}&-4&-1& 4\\ 2&-3&-1& 4\\ 2&-5&-3& 9\\ \end{array} \right) \quad \Rightarrow\quad \left( \begin{array}{rrr|r} 4&-4&-1& 4\\ 0&-1&-1/2& 2\\ 0&\underline{-3}&-5/2& 7\\ \end{array} \right) \quad\equiv\quad \left( \begin{array}{rrr|r} 4&-4&-1& 4\\ 0&\underline{-3}&-5/2& 7\\ 0&-1&-1/2& 2\\ \end{array} \right)$$
• Gauss $$\left( \begin{array}{rrr|r} 4&-4&-1& 4\\ 0&\underline{-3}&-5/2& 7\\ 0&-1&-1/2& 2\\ \end{array} \right) \quad \Rightarrow\quad \left( \begin{array}{rrr|r} 4&-4&-1& 4\\ 0&-3&-5/2& 7\\ 0& 0&1/3& -1/3\\ \end{array} \right)$$
• back substitution $$z=-1 \quad \qquad -3y=7+5z/2=9/2$$ $$y = -3/2 \qquad\quad 4x=4+4y+z=-3 \qquad\quad x =-3/4$$
• final (compare previous result with $w=0$) $$\vec x= \left(\begin{array}{c} x \\ y \\ z \end{array}\right) = \left(\begin{array}{c} -3/4 \\ -3/2 \\ -1 \end{array}\right)$$

# Gauss elimination with pivoting

• with swaping of pivot rows
• testing of result
• input from file
• note (&A) in function argument:
the reference of an array contains information about the length, a pointer not
int  intNumber    = 10;
int  intArray[]   = { 1, 2, 3, 4, 5 };

int* intPointer_1 = &intNumber;      // ok
int* intPointer_2 = &intArray;       // not allowed : cannot convert ‘int (*)[5]’ to ‘int*’
int* intPointer_3 = &(intArray[0]);  // ok

• note that calling of the array-length templated function gaussElimination()
with dynamic arrays would not be possible
#include <iostream>
#include <stdio.h>     // for printf
#include <string>      // c_str()
#include <cmath>       // for abs
#include <fstream>     // file streams
using namespace std;

// --- Gauss elimination of Ax=b with pivoting.
// --- On input quadratic matrix A and vector b.
// --- On ouput solution x.

template<int N>
void gaussElimination(double (&A)[N][N], double (&b)[N], double (&x)[N])
{
double temp;                       // temp variable
const double smallNumber = 1e-10;  // a small number
//
//--- loop over all columns A[][col];  A[][col]
for (int col=0; col<N; col++)
{
//--- find pivot row; swap pivot with actualy row
int pivotRow = col;
for (int row=col+1;row<N;row++)
if ( abs(A[row][col]) > abs(A[pivotRow][col]) )
pivotRow = row;

//--- swap row A[col][] with A[pivotRow][]
for (int i=0;i<N;i++)
{
temp = A[col][i];
A[col][i] = A[pivotRow][i];
A[pivotRow][i] = temp;
}

//--- swap elements b[col] with b[pivotRow]
temp = b[col];
b[col] = b[pivotRow];
b[pivotRow] = temp;

//--- throw a warning if matrix is singular or nearly singular
if (abs(A[col][col]) <= smallNumber)
cerr << "Matrix is singular or nearly singular" << endl;

//--- Gauss with pivot within A and b
for (int row=col+1;row<N;row++)
{
double alpha = A[row][col] / A[col][col];  // divide by pivot
b[row] -= alpha*b[col];
for (int j=col;j<N;j++)
A[row][j] -= alpha * A[col][j];
}
}  // end of loop over all columns

//--- back substitution
for (int i=N-1;i>=0;i--)
{
double sum = 0.0;
for (int j=i+1;j<N;j++)
sum += A[i][j] * x[j];
x[i] = (b[i] - sum) / A[i][i];   // pivots on diagonal after swapping
}
}  // end of gaussElimination()

// ---
// --- print system
// ---
template<int N>
void printSystem(double A[N][N], double b[N], double x[], string message)
{
printf("\n---------------- %s\n",message.c_str());
for (int col=0;col<N;col++)
printf("%7d",col);
printf(" |%7s   :: %7s\n","b","x");
//
for (int row=0;row<N;row++)
{
for (int col=0;col<N;col++)
printf("%7.3f",A[row][col]);
printf(" |%7.3f   :: %7.3f\n",b[row],x[row]);
}
} // end of printSystem()

// ---
// --- main
// ---
int main()
{
/*
const int N = 2;
double A0[N][N], A[N][N], b0[N], b[N], x[N];
A0[0][0] = 1.0;  A0[0][1] = 3.0;   b0[0] = 1;
A0[1][0] = 2.0;  A0[1][1] = 4.0;   b0[1] = 1;
*/
const int N = 3;
double A0[N][N], A[N][N], b0[N], b[N], x[N];
//
A0[0][0] = 2.0;  A0[0][1] =-3.0; A0[0][2] = -1.0;   b0[0] =  4.0;
A0[1][0] = 4.0;  A0[1][1] =-4.0; A0[1][2] = -1.0;   b0[1] =  4.0;
A0[2][0] = 2.0;  A0[2][1] =-5.0; A0[2][2] = -3.0;   b0[2] =  9.0;
//
// an example for finding an EES (evolutionay stable strategy; game theory)
//
A0[0][0] = 6.0;  A0[0][1] = 5.0; A0[0][2] =  5.0;   b0[0] = 28.0;
A0[1][0] = 7.0;  A0[1][1] = 4.0; A0[1][2] =  3.0;   b0[1] = 28.0;
A0[2][0] = 7.0;  A0[2][1] = 5.0; A0[2][2] =  2.0;   b0[2] = 28.0;
//
for (int row=0;row<N;row++)
{
b[row] = b0[row];
for (int col=0;col<N;col++)
A[row][col] = A0[row][col];
}
//
printSystem(A,b,x,"at the start");
gaussElimination(A,b,x);
printSystem(A,b,x,"after elimination (A,b)");
printSystem(A0,b0,x,"after elimination (A0,B0)");
//
printf("\n");
for (int row=0;row<N;row++)  // testing the result
{
double sum = 0.0;
for (int col=0;col<N;col++)
sum += A0[row][col]*x[col];
printf("# row %d %8.3f | %8.3f\n",row, sum, b0[row]);
}
//
return 1;
}  // end of main()


# matrices

• diagonal
• $$\left( \begin{array}{ccccc} *&0&0&0&0\\ 0&*&0&0&0\\ 0&0&*&0&0\\ 0&0&0&*&0\\ 0&0&0&0&*\\ \end{array} \right)$$
• tri-diagonal
• $$\left( \begin{array}{ccccc} *&*&0&0&0\\ *&*&*&0&0\\ 0&*&*&*&0\\ 0&0&*&*&*\\ 0&0&0&*&*\\ \end{array} \right)$$
• block-diagonal (blocks are independent
• $$\left( \begin{array}{ccccc} *&*&0&0&0\\ *&*&0&0&0\\ 0&0&*&*&0\\ 0&0&*&*&0\\ 0&0&0&0&*\\ \end{array} \right)$$
• L (lower left) and R (upper right) triangular matrices
• $$L = \left( \begin{array}{ccccc} *&0&0&0&0\\ *&*&0&0&0\\ *&*&*&0&0\\ *&*&*&*&0\\ *&*&*&*&*\\ \end{array} \right) \qquad\quad R = \left( \begin{array}{ccccc} *&*&*&*&*\\ 0&*&*&*&*\\ 0&0&*&*&*\\ 0&0&0&*&*\\ 0&0&0&0&*\\ \end{array} \right)$$
• sparse (few non-zero entries, there are special algorithms)
• $$\left( \begin{array}{ccccc} 0&*&0&0&0\\ *&0&0&*&0\\ 0&*&0&0&0\\ *&0&0&0&*\\ 0&0&*&0&0\\ \end{array} \right)$$

# LR decomposition

• decompose matrix into product of $L$ and $R$ matrices
$$A = L\,R \qquad\qquad A\,\vec x =\vec b$$
• always possible (Gauss, Hausholder, ...)
number of variables: $\quad 2 N(N+1)/2 = N^2+N$
$N$ free variables
$$L\,R\,\vec x =\vec b \qquad\qquad L\,\vec y =\vec b \qquad\qquad \vec y = R\,\vec x$$
• forward substitution: solve for $\ \vec y$
$$L\,\vec y = \vec b \qquad\qquad L = \left( \begin{array}{cccc} *&0&0&0\\ *&*&0&0\\ *&*&*&0\\ *&*&*&*\\ \end{array} \right)$$
• backward substitution: solve for $\ \vec x$
$$R\,\vec x = \vec y \qquad\qquad R = \left( \begin{array}{cccc} *&*&*&*\\ 0&*&*&*\\ 0&0&*&*\\ 0&0&0&*\\ \end{array} \right)$$

# LR decomposition with Gauss elimination

• elimination is a linear process
:: equivalent to a matrix multiplication
here: no pivoting $$A = \left( \begin{array}{ccc} 1&1&1 \\ 2&3&5 \\ 4&6&8 \\ \end{array} \right)\qquad\quad L_1A = \left( \begin{array}{ccc} 1&0&0 \\ -2&1&0 \\ -4&0&1 \\ \end{array} \right) \left( \begin{array}{ccc} \underline{1}&1&1 \\ 2&3&5 \\ 4&6&8 \\ \end{array} \right) = \left( \begin{array}{ccc} 1&1&1 \\ 0&1&3 \\ 0&2&4 \\ \end{array} \right)$$
• Gauss for second column $$L_2\big(L_1A\big) = \left( \begin{array}{ccc} 1&0&0 \\ 0&1&0 \\ 0&-2&1 \\ \end{array} \right) \left( \begin{array}{ccc} 1&1&1 \\ 0&\underline{1}&3 \\ 0&2&4 \\ \end{array} \right) = \left( \begin{array}{ccc} 1&1&1 \\ 0&1&3 \\ 0&0&-2 \\ \end{array} \right) \equiv R$$
• inversion $$L_2L_1A=R\qquad\quad A=(L_1)^{-1}(L_2)^{-1}R\equiv LR \qquad\quad L=(L_1)^{-1}(L_2)^{-1}$$ with $$(L_1)^{-1}= \left( \begin{array}{ccc} 1&0&0 \\ 2&1&0 \\ 4&0&1 \\ \end{array} \right) \qquad\quad (L_1)^{-1}L_1= \left( \begin{array}{ccc} 1&0&0 \\ 2&1&0 \\ 4&0&1 \\ \end{array} \right) \left( \begin{array}{ccc} 1&0&0 \\ -2&1&0 \\ -4&0&1 \\ \end{array} \right) = \left( \begin{array}{ccc} 1&0&0 \\ 0&1&0 \\ 0&0&1 \\ \end{array} \right)$$
• the left-matrix $\ L\$ has therefore one's on the diagonals
(the $N$ degrees of freedom): $$L=(L_1)^{-1}(L_2)^{-1}= \left( \begin{array}{ccc} 1&0&0 \\ 2&1&0 \\ 4&0&1 \\ \end{array} \right) \left( \begin{array}{ccc} 1&0&0 \\ 0&1&0 \\ 0&2&1 \\ \end{array} \right) = \left( \begin{array}{ccc} 1&0&0 \\ 2&1&0 \\ 4&2&1 \\ \end{array} \right)$$
• check: $$LR= \left( \begin{array}{ccc} 1&0&0 \\ 2&1&0 \\ 4&2&1 \\ \end{array} \right) \left( \begin{array}{ccc} 1&1&1 \\ 0&1&3 \\ 0&0&-2 \\ \end{array} \right) = \left( \begin{array}{ccc} 1&1&1 \\ 2&3&5 \\ 4&6&8 \\ \end{array} \right) =A$$