Linear Systems : Matrices, Vectors, eigen systems

Linear Systems : Matrices, Vectors, eigen systems#

In this module we will learn how to solve linear systems which are very common in engineering and applied physics.

Applications are numerous:

One particular common operation, the matrix multiplication, is still the subject of ongoing research:

Image Description
From: "https://www.quantamagazine.org/wp-content/uploads/2021/03/Matrix_multiplication_2880x1620_Lede.jpg"

Tips about matrix computing:

Furthermore, for eigen values, please check https://www.youtube.com/watch?v=PFDu9oVAE-g&t=0s

Eigen c++#

Eigen is a C++ template library for linear algebra: matrices, vectors, numerical solvers, and related algorithms.

Quick reference: https://eigen.tuxfamily.org/dox/group__QuickRefPage.html

Solving \(Ax = b\)#

This is one of the archetypical problems in numerical algebra. You migh try to implement it by hand, but it will probable be not efficient, stable, etc (Remember https://nhigham.com/2022/10/11/seven-sins-of-numerical-linear-algebra/). Better to use a numerical library to do the heavy task for you.

QR decomposition#

In this example we are going to use eigen to solve the system using HouseHolder QR decomposition as done in http://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html . Please create a file with the following code, compile it and run it. You are supposed to have already installed eigen (the computer room and the binder machine both have it).

# include <iostream>
# include <Eigen/Dense>

int main()
{
   Eigen::Matrix3d A;
   Eigen::Vector3d b;
   //std::cout.precision(16);
   //std::cout.setf(std::ios::scientific);
   
   A << 1,2,3,  4,5,6,  7,8,10;
   b << 3, 3, 4;
   std::cout << "Here is the matrix A:\n" << A << std::endl;
   std::cout << "Here is the vector b:\n" << b << std::endl;
   Eigen::Vector3d x = A.colPivHouseholderQr().solve(b);
   std::cout << "The solution is:\n" << x << std::endl;
   std::cout << (A*x - b).norm() << "\n";
   return 0;
}

If the library is installed on system paths, just compile it as

g++ -std=c++17 qr.cpp -o qr.x

or compile it optimized as

g++ -std=c++17 -O3 qr.cpp -o qr.x

and then run it

./qr.x

What do you get? is the optimized code faster than the original one? why? why not?

Installation and compilation tips#

Try the following:

  • Library installed at /usr/include/eigen3 or /usr/local/include/eigen3: Then you need to tell the compiler where to find the headers

    g++ -std=c++17 -O3 -I /usr/include/eigen3 qr.cpp -o qr.x

    or, if you installed it on a specific path

    g++ -std=c++17 -O3 -I /path/to/installdir/ qr.cpp -o qr.x

  • Install it : On systems like ubuntu, you could use

    sudo apt install libeigen3-dev

    In collab,

    apt install libeigen3-dev

    In the course provisioned binder, you do not need to install, it is already installed.

Exercise#

How to know that the solution is actually a solution of the original problem? design a criteria .

Exercise#

Now create a random matrix and vector, of arbitrary size N, and measure the time to compute the solution as a function of N. Use Eigen::MatrixXd A = Eigen::MatrixXd::Random(N, N); and std::chrono

LU decomposition#

THe LU decomposition os another tool we can use to solve the system (advantages? disadvantages?)

Now implement the following solution and compare with the previous one? The only change you need is

Eigen::MatrixXd x = A.fullPivLu().solve(b)

%%writefile lu.cpp
# YOUR CODE HERE
raise NotImplementedError()

Writing lu.cpp

Exercises#

Solving simple system#

Solve the system

../../_images/linear-example-01.png

Performance between QR and LU#

Compare the time between LU and QR decompositions on random matrices.

Solve simple system#

Simulating temperature#

  • Temperature discretized

  • System of equations

Laplace equation in 2D#

When solving the laplace equation \(\nabla V = 0\) to compute the electrostatic potential on a planar region, you can discretize the derivatives on a grid and then arrive to the following equation

(37)#\[\begin{equation} V_{i+1,j} + V_{i-1,j} + V_{i,j+1} + V_{i,j-1} - 4V_{i,j} = 0.0. \end{equation}\]

This can be written as a matrix problem. Solve the matrix problem for a square plate of lenght \(L\), with \(N\) points on each side. The boundary conditions, of Dirichlet type, are \(V(x, 0) = 5\sin(\pi x/L)\), \(V(x, L) = V(0, y) = V(L, y) = 0.0\).

Vandermonde determinant#

Compute the determinant of the [[https://en.wikipedia.org/wiki/Vandermonde_matrix][Vandermonde matrix]] of size \(N\times N\). Measure the time as a function of \(N\).

Condition number#

Compute the condition number for an arbitrary matrix \(A\). Apply it for the Vandermonde matrix. The condition number is defined as \(\kappa(A) = |A^{-1}| |A|\)

Rotation matrix#

Define a rotation matrix in 2d by an angle \(\theta\). Apply it to a given vector and check that it is actually rotated.

Coupled oscillators#

Imagine that you have two masses \(m_1, m_2\), coupled through springs in the form wall-spring-mass-spring-mass-spring-wall. All spring are linear with constant \(k\). Write the equations of motion, replace each solution with \(x_i(t) = A_i e^{i\omega t}\), and obtain a matrix representation to get the amplitudes. Compute the eigen values and eigen-vectors. Those are the [[https://en.wikipedia.org/wiki/Normal_mode][normal modes]] . Extend to n oscillators of the same mass.

Thick lens (Boas, 3.15.9)#

The next matrix is used when discussing a thick lens in air

(38)#\[\begin{equation} A = \begin{pmatrix} 1 & (n-1)/R_2\\ 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0\\ d/n & 1 \end{pmatrix} \begin{pmatrix} 1 & -(n-1)/R_1\\ 0 & 1 \end{pmatrix}, \end{equation}\]

where \(d\) is the thickness of the lens, \(n\) is the refraction index, and \(R_1\) and \(R_2\) are the curvature radius. Element \(A_{12}\) is equal to \(-1/f\), where \(f\) is the focal distance. Evaluate \(\det A\) and \(1/f\) as functions of \(n \in [1, 3]\).

Products production#

Teaching distribution#

Payments#

../../_images/problem-08.png

Eigen vectors and eigen values#

Eigen-values (\(\lambda\)) and eigen-vectors (\(\vec v\)) of a matrix \(A\) fulfill the following equation

(39)#\[\begin{equation} A\vec v = \lambda \vec v, \end{equation}\]

and are extensively use in physics, statistics, computational methods, etc (see https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors?useskin=vector). For example, they are use to compute (https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors?useskin=vector#Applications):

  • Eigen values of the Schroedinger equation

    Image Description
    From: "https://upload.wikimedia.org/wikipedia/commons/c/cf/HAtomOrbitals.png"

  • Diagonalization basis for matrices: Given a matrix \(A\), and the matrix \(P\) whose column vectors are the eigenvectors of \(A\), then one can apply a basis trnasformation, to move to the eigen-vector matrices, and obtain a diagonal matrix \(D\) (whose elements are the eigen values of \(A\)), as

    (40)#\[\begin{equation} D = P^{-1}AP \end{equation}\]
    Image Description
    From: "https://upload.wikimedia.org/wikipedia/commons/4/4e/Diagonalization_as_rotation.gif"

  • Normal mode of oscillatory systems

  • Principal component analysis in data science

https://towardsdatascience.com/wp-content/uploads/2022/02/1QinDfRawRskupf4mU5bYSA.png

  • Image analysis

    Image Description
    From: "https://upload.wikimedia.org/wikipedia/commons/6/67/Eigenfaces.png"

  • Moment of inertia

  • Stress tensor principal axis

Solving a general eigen value problem with eigen c++#

Following the documentation of the general eigen solver, https://eigen.tuxfamily.org/dox/classEigen_1_1EigenSolver.html, we can use the example to show how it works:

#include <iostream>
#include <string>
#include <Eigen/Dense>
#include <complex>

int main(int argc, char **argv) {

    // read the matrix ncols, to be NxN
    const int N = std::stoi(argv[1]);

    // create a random matrix
    // seed controlled by srand
    Eigen::MatrixXd A = Eigen::MatrixXd::Random(N,N);
    std::cout << "Here is a random matrix, A:" << std::endl 
            << A << std::endl << std::endl;

    // create an eigen solver with that matrix. 
    // Once created it solves everything inmediatly
    Eigen::EigenSolver<Eigen::MatrixXd> es(A);

    // print some results
    std::cout << "The eigenvalues of A are:" << std::endl 
            << es.eigenvalues() << std::endl;
    std::cout << "The matrix of eigenvectors, V, is:" << std::endl 
            << es.eigenvectors() << std::endl << std::endl;

    // Check the first eigen value: Av = lambda v
    std::complex<double> lambda = es.eigenvalues()[0];
    std::cout << "Consider the first eigenvalue, lambda = " << lambda << std::endl;
    Eigen::VectorXcd v = es.eigenvectors().col(0);
    std::cout << "If v is the corresponding eigenvector, then lambda * v = " << std::endl 
            << lambda * v << std::endl;
    std::cout << "... and A * v = " << std::endl 
            << A.cast<std::complex<double> >() * v << std::endl << std::endl;

    // diagnalization check
    Eigen::MatrixXcd D = es.eigenvalues().asDiagonal();
    Eigen::MatrixXcd V = es.eigenvectors();
    std::cout << "Finally, V * D * V^(-1) = " << std::endl 
            << V * D * V.inverse() << std::endl;
    
    return 0;
}

Please compile and run it. This is a very general example and cam be simplified for specific matrices (hermitian) to improve efficiency.

Matrix exponential from taylor expansion#

In scalar calculus, the exponential function \(e^x\) can be written as an infinite sum using its Taylor series (specifically, the Maclaurin series) centered at \(0\):

\[e^x = \sum_{k=0}^{\infty} \frac{x^k}{k!} = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \dots\]

We can extend this exact concept to linear algebra. If we have a square matrix \(A\) (of size \(n \times n\)), the matrix exponential, denoted as \(e^A\), is defined by substituting the matrix into the same Taylor expansion.

The Definition#

Since we are dealing with matrices, the scalar \(1\) becomes the identity matrix \(I\):

\[e^A = \sum_{k=0}^{\infty} \frac{A^k}{k!} = I + A + \frac{A^2}{2!} + \frac{A^3}{3!} + \dots\]

Where:

  • \(A^0 = I\) (Identity matrix)

  • \(A^2 = A \times A\)

  • \(k!\) is the factorial of \(k\) (\(k \times (k-1) \times \dots \times 1\))

Tip

Why is this important? Matrix exponentials are crucial for solving systems of linear differential equations. For example, the system \(\frac{d\mathbf{x}}{dt} = A\mathbf{x}\) has the elegant solution \(\mathbf{x}(t) = e^{At}\mathbf{x}(0)\).

How to Compute It (An Example)#

While the infinite sum is the definition, computing it directly by adding infinite matrices is usually impossible unless the matrix has special properties.

Scenario 1: Nilpotent Matrices (The Series Terminates)#

If a matrix becomes the zero matrix when raised to a certain power (called a nilpotent matrix), the Taylor series becomes finite and easy to calculate.

Let \(A = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}\). Let’s find \(A^2\): $\(A^2 = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} = \begin{pmatrix} 0 & 0 \\ 0 & 0 \end{pmatrix}\)$

Because \(A^2 = 0\), all higher powers (\(A^3, A^4, \dots\)) are also \(0\). The infinite Taylor series collapses into just the first two terms:

\[\begin{split}e^A = I + A = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} + \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix}\end{split}\]

Scenario 2: Diagonal Matrices#

If a matrix is diagonal, you can simply take the exponential of each diagonal element:

\[\begin{split}\text{If } D = \begin{pmatrix} a & 0 \\ 0 & b \end{pmatrix}, \text{ then } e^D = \begin{pmatrix} e^a & 0 \\ 0 & e^b \end{pmatrix}\end{split}\]

For more complex matrices, mathematicians usually diagonalize the matrix (\(A = PDP^{-1}\)) so that \(e^A = Pe^DP^{-1}\).

Key Properties to Keep in Mind#

  • Non-Commutativity: In scalar math, \(e^x e^y = e^{x+y}\). In matrix math, this is only true if the matrices commute (\(AB = BA\)). If they do not commute, \(e^A e^B \neq e^{A+B}\).

  • Inverse: The inverse of \(e^A\) is always \(e^{-A}\).

  • Determinant: \(\det(e^A) = e^{\text{trace}(A)}\).

Normal modes and matricws#

In physics and engineering, normal modes describe the natural, synchronized patterns of vibration in an interconnected system (like coupled pendulums, buildings during an earthquake, or atoms in a molecule).

Linear algebra—specifically matrices, eigenvalues, and eigenvectors—is the mathematical engine used to find these modes.

What is a Normal Mode?#

In a normal mode, all parts of the system vibrate at the exact same frequency and pass through their equilibrium positions at the exact same time. If you have a system with \(n\) degrees of freedom (e.g., \(n\) masses connected by springs), the system will have exactly \(n\) distinct normal modes. Any complex, chaotic-looking motion of the system is actually just a combination (superposition) of these fundamental normal modes.

How Matrices Represent the System#

Consider two identical masses, \(m\), connected to each other and to walls by identical springs with spring constant \(k\).

If we apply Newton’s second law (\(F = ma\)) to both masses, we get a system of coupled differential equations:

\[m \frac{d^2x_1}{dt^2} = -kx_1 - k(x_1 - x_2) = -2kx_1 + kx_2\]
\[m \frac{d^2x_2}{dt^2} = -kx_2 - k(x_2 - x_1) = kx_1 - 2kx_2\]

We can rewrite this entire system beautifully using a single matrix equation:

\[\begin{split}\begin{pmatrix} m & 0 \\ 0 & m \end{pmatrix} \begin{pmatrix} \ddot{x}_1 \\ \ddot{x}_2 \end{pmatrix} + \begin{pmatrix} 2k & -k \\ -k & 2k \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}\end{split}\]

More compactly, this is written as: $\(M\ddot{\mathbf{x}} + K\mathbf{x} = 0\)$

Where:

  • \(M\) is the Mass Matrix

  • \(K\) is the Stiffness Matrix (notice how the off-diagonal \(-k\) terms “couple” the motion of the two masses together).

Solving for Modes using Eigenvalues#

To find a normal mode, we look for a solution where all parts vibrate sinusoidally at the same frequency \(\omega\): $\(\mathbf{x}(t) = \mathbf{v} \cos(\omega t)\)$

If we take the second derivative of this guess (\(\ddot{\mathbf{x}} = -\omega^2 \mathbf{v} \cos(\omega t)\)) and plug it back into our matrix equation, the calculus vanishes, leaving us with a pure linear algebra problem:

\[K\mathbf{v} = \omega^2 M\mathbf{v}\]

If we assume \(M\) is the identity matrix (or multiply by \(M^{-1}\)), this becomes a classic Eigenvalue Problem:

\[A\mathbf{v} = \lambda\mathbf{v}\]

Where:

  • The matrix \(A = M^{-1}K\)

  • The eigenvalue \(\lambda = \omega^2\) gives the frequency of the vibration (\(\omega = \sqrt{\lambda}\)).

  • The eigenvector \(\mathbf{v}\) gives the shape of the vibration (the relative displacement of each mass).

The Physical Meaning of the Eigenvectors#

For the two-mass system described above, solving the matrix equation yields two eigenvalues and two corresponding eigenvectors:

Mode 1: Symmetric Motion (Low Frequency)#

  • Eigenvalue: \(\lambda_1 = \frac{k}{m} \implies \omega_1 = \sqrt{\frac{k}{m}}\)

  • Eigenvector: \(\mathbf{v}_1 = \begin{pmatrix} 1 \\ 1 \end{pmatrix}\)

  • Physical Meaning: The masses move in the exact same direction with the exact same amplitude. Because they move together, the middle spring never stretches or compresses, meaning the system vibrates at a lower frequency.

Mode 2: Antisymmetric Motion (High Frequency)#

  • Eigenvalue: \(\lambda_2 = \frac{3k}{m} \implies \omega_2 = \sqrt{\frac{3k}{m}}\)

  • Eigenvector: \(\mathbf{v}_2 = \begin{pmatrix} 1 \\ -1 \end{pmatrix}\)

  • Physical Meaning: The masses move in opposite directions (mirror images). This fiercely compresses and stretches the center spring, resulting in a much higher force and a higher frequency of vibration.

Connecting Back to the Matrix Exponential#

Remember how the matrix exponential \(e^{At}\) solves \(\frac{d\mathbf{x}}{dt} = A\mathbf{x}\)?

For second-order system equations like \(M\ddot{\mathbf{x}} + K\mathbf{x} = 0\), we can use Euler’s formula extensions with matrices (\(\cos(\sqrt{M^{-1}K}t)\)) or convert the system into a larger first-order system. Ultimately, the matrix exponential provides the analytical trajectory of the system over time, while the eigenvectors map out the geometric pathways (the normal modes) that the system naturally flows through.

Exercises for eigen systems#

Eigen values Hilbert matrix#

Compute the eigen values and the condition number for the https://en.wikipedia.org/wiki/Hilbert_matrix.

Power method and eigen values#

Apply the https://en.wikipedia.org/wiki/Power_iteration to compute the maximum eigen value of the Vandermonde matrix (or any other matrix).

Power of a matrix#

Given a matrix \(A\), one can compute easily \(A^k\) using diagonalization: \(A^k = (PDP^{-1})^k = PDP^{-1}PDP^{-1} ... PDP^{-1} = PD^{k}P^{-1}\), where \(D^k\) is just raising the eigen values to the given k power. Implement this in a function that receives an eigen matrix and a power k and returns the corresponding result.

Checking stability (Heath, scientific computing)#

Compute the eigen values of the matrix

(41)#\[\begin{equation} A = \begin{pmatrix} 9 & 4.5 & 3\\ -56 & -28 & -18\\ 60 & 30 & 19 \end{pmatrix} \end{equation}\]

Then do it again, but change the 19 to 18.95. What is the relative change in the eigen values magnitudes? Do the same changing it to 19.05. Can you draw any conclusion about this? can you devise a condition number to characterize this?

Markov chain (Heath, Scientific Computing)#

A Markov chain is a representation of a Markov process, which, in turns, is a random process with no memory: the current state depends only on the previous one and no more. The probability of a transition from state j to state i is given by \(a_{ij}\),where \(0 \le a_{ij} \le 1\) and \(\sum_{i=1}^n a_{ij} =1\). Let \(A\) denote the matrix of transition probabilities, and let \(x^{(k)}\) denote the probability that the system is in state \(i\) after transition \(k\). If the initial probability distribution vector is \(x^{(0)}\), then the probability distribution vector after \(k\) steps is given by

(42)#\[\begin{equation} x^{(k)} = A x^{(k-1)} = A^k x^{(0)}. \end{equation}\]

Consider the matrix

(43)#\[\begin{equation} A = \begin{pmatrix} 0.8 & 0.2 & 0.1\\ 0.1 & 0.7 & 0.3\\ 0.1 & 0.1 & 0.6\\ \end{pmatrix}, \end{equation}\]

and suposse the system is initially in state 1.

  • What is the probability distribution vector after three steps?

  • What is the long-term value of the probability distribution vector?

  • Does the long-term value of the probability distribution vector depend on the particular starting value ?

  • What is the value of \(\lim_{k \to \infty} A^k\) , and what is the rank of this matrix?

  • Explain your previous results in terms of the eigenvalues and eigenvectors of A

  • Must 1 always be an eigenvalue of the transition matrix of a Markov chain? Why?

  • A probability distribution vector x is said to be stationary if Ax = x. How can you determine such a stationary value x using the eigenvalues and eigenvectors of A?

  • How can you determine a stationary value x without knowledge of the eigenvalues and eigen- vectors of A?

Other eigen tools#