Numerical Linear Algebra

Linear algebra is a fundamental part of modern mathematics. It supports fields from calculus and statistics to geometry, control theory, and functional analysis. Most linear systems are well understood. Even nonlinear problems are often studied through linear approximations.

Numerical linear algebra extends these ideas to computation, enabling solutions of PDEs, optimization tasks, eigenvalue problems, and more. It addressing some of the hardest problems in physics and engineering.

Motivations from Physics¶

• Normal Modes: Vibrations near equilibrium reduce to generalized eigenvalue problems. Linear algebra therefore reveals resonance in materials, acoustics, and plasma waves.

• Quantum Mechanics: Described by the Schrödinger equation, quantum systems are inherently linear.

• Discretized PDEs: Discretizing PDEs yields large sparse linear systems. They can solved numerically by methods such as conjugate gradient.

• Nonlinear Problems: Nonlinear physics problems including turbulence are sometimes untrackable. Linearizing them with perturbation theory reduces them to sequences of linear systems.

Motivations from Computation¶

• Large-Scale Data: Modern sensors and simulations produce massive datasets. Matrix decompositions (e.g., SVD, PCA) provide compression, noise reduction, and feature extraction.

• Neural Networks: Core operations in training, i.e., backpropagation, is dominated by large matrix multiplications. Efficient linear algebra routines are therefore critical for scaling deep learning.

• Hardware Accelerators: GPUs and TPUs are optimized for matrix operations, making vectorized linear algebra essential for both neural networks and scientific computing.

Condition Number and Accuracy in Solving Linear Equations¶

When solving a linear system $A \mathbf{x} = \mathbf{b}$, accuracy depends not only on the algorithm but also on the conditioning of the matrix. The condition number describes how sensitive the solution is to small perturbations in the input.

Definition (2-norm)¶

• If $\kappa(A) \approx 1$, the system is well-conditioned.

• If $\kappa(A) \gg 1$, the system is ill-conditioned, and small changes in $A$ or $\mathbf{b}$ can produce large errors in $\mathbf{x}$.

Error Amplification¶

If $\mathbf{b}$ is perturbed by $\delta\mathbf{b}$, the solution changes by $\delta\mathbf{x}$. For small perturbations,

Thus the relative error in $\mathbf{x}$ can be amplified by a factor of up to $\kappa(A)$.

Interpretation¶

• $\kappa(A) \sim 10^0$--10²: problem is well-behaved; only a few digits may be lost.

• $\kappa(A) \sim 10^6$: up to 6 digits of accuracy may be lost (in double precision with $\sim 16$ digits).

• $\kappa(A) \sim 10^{16}$: essentially singular in double precision; solutions may be dominated by numerical noise.

In practice, $\log_{10}\kappa(A)$ estimates the number of accurate digits lost in the solution.

Direct Solvers¶

Direct methods are often the first approach taught for solving linear systems $A\mathbf{x} = \mathbf{b}$. They involve algebraic factorizations that can be computed in a fixed number of steps (roughly $\mathcal{O}(n^3)$) for an $n \times n$ matrix.

Gaussian Elimination¶

Gaussian Elimination transforms the system $A \mathbf{x} = \mathbf{b}$ into an equivalent upper-triangular form $U \mathbf{x} = \mathbf{c}$ by systematically applying row operations. Once in an upper-triangular form, one can perform back-substitution to solve for $\mathbf{x}$.

1. Row Operations

   • Subtract a multiple of one row from another to eliminate entries below the main diagonal.

   • Aim to create zeros in column $j$ below row $j$.

2. Partial Pivoting (optional)

   • When a pivot (diagonal) element is small (or zero), swap the current row with a row below that has a larger pivot element in the same column.

   • This step mitigates numerical instability by reducing the chance that small pivots lead to large rounding errors in subsequent operations.

3. Result

   • After eliminating all sub-diagonal entries, the matrix is in upper-triangular form $U$.

   • Solve $U\mathbf{x} = \mathbf{c}$ via back-substitution.

Here is an example:

The matrix is now in echelon form (also called triangular form):

Below is a simple code for naive (no pivoting) Gaussian Elimination in python. Although normally we want to avoid for loops in python for performance, let’s stick with for loop this time so we can directly implement the algorithm we just described.

Let’s also compare with numpy’s solver.

Let’s set the (0,0) element to a small value.

Let’s now improve the above naive (no pivoting) Gaussian elimination by adding pivoting.

$LU$ Decomposition¶

$LU$ Decomposition is a systematic way to express $A$ as $A = P L U$, where:

• $P$ is a permutation matrix

• $L$ is lower triangular (with 1s on the diagonal following the standard convention).

• $U$ is upper triangular.

Once $A$ is factored as $P L U$, solving $A\mathbf{x} = \mathbf{b}$ becomes:

1. $L \mathbf{y} = P^t \mathbf{b}$ (forward substitution)

2. $U \mathbf{x} = \mathbf{y}$ (back-substitution)

Gaussian elimination essentially constructs the $L$ and $U$ matrices behind the scenes:

• The multipliers used in the row operations become the entries of $L$.

• The final upper-triangular form is $U$.

Matrix Inverse¶

Although the matrix inverse $A^{-1}$ is a central theoretical concept, explicitly forming $A^{-1}$ just to solve $A \mathbf{x} = \mathbf{b}$ is almost always unnecessary and can degrade numerical stability. Instead, direct approaches (like Gaussian Elimination or LU factorization) find $\mathbf{x}$ with fewer operations and less error accumulation. Both methods cost about $\mathcal{O}(n^3)$, but computing the entire inverse introduces extra steps and can magnify floating-point errors.

In rare cases, you might actually need $A^{-1}$.

For example, when:

• Multiple Right-Hand Sides: If you must solve $A \mathbf{x} = \mathbf{b}_i$ for many different $\mathbf{b}_i$, you might form or approximate $A^{-1}$ for convenience.

• Inverse-Related Operations: Some advanced algorithms (e.g., computing discrete Green’s functions or certain control theory design methods) explicitly require elements of the inverse.

In the demonstration below, we show how to compute the inverse using np.linalg.inv(). However, again, it is generally safer and more efficient to solve individual systems via a factorization method in practice.

Iterative Solvers for Large/Sparse Systems¶

Direct methods are powerful but quickly become impractical when dealing with very large or sparse systems. In such cases, iterative solvers provide a scalable alternative.

• Large-Scale Problems: PDE discretizations in 2D/3D can easily lead to systems with millions of unknowns, making $\mathcal{O}(n^3)$ direct approaches infeasible.

• Sparse Matrices: Many physical systems produce matrices with a significant number of zero entries. Iterative methods may be implemented to access only nonzero elements each iteration, saving time and memory.

• Scalability: On parallel architectures (clusters, GPUs, etc), iterative methods often scale better than direct factorizations.

Jacobi Iteration¶

1. Idea: Rewrite $A\mathbf{x} = \mathbf{b}$ as $\mathbf{x} = D^{-1}(\mathbf{b} - R\mathbf{x})$, where $D$ is the diagonal of $A$ and $R$ is the remainder.

2. Update Rule:

   $$\begin{align} x_i^{(k+1)} = \frac{1}{a_{ii}} \Big(b_i - \sum_{j \neq i} a_{ij} x_j^{(k)}\Big). \end{align}$$

   (20)

3. Pros/Cons:

   • Easy to implement: each iteration updates $\mathbf{x}$ using only values from the previous iteration.

   • Slow convergence unless $A$ is well-conditioned (e.g., diagonally dominant).

Gauss-Seidel Iteration¶

Gauss-Seidel iteration is closely related to Jacobi but usually converges faster, particularly if the matrix is strictly or diagonally dominant.

1. Idea:

   • Instead of using only the old iteration values (from step $k$), Gauss-Seidel uses the most recent updates within the same iteration.

   • This often converges faster because you incorporate newly computed values immediately rather than waiting for the next iteration.

2. Update Rule:

   $$\begin{align} x_i^{(k+1)} = \frac{1}{a_{ii}} \Big( b_i - \sum_{j < i} a_{ij} x_j^{(k+1)} - \sum_{j > i} a_{ij} x_j^{(k)} \Big). \end{align}$$

   (21)

   Note that $\mathbf{x}^{(k+1)}$ is already partially updated (for $j < i$).

3. Convergence Properties:

   • Gauss-Seidel is shown to converge for strictly diagonally dominant matrices, and often outperforms Jacobi in practice.

   • Still relatively slow for large, ill-conditioned systems.

Eigenvalue Problems¶

Eigenvalue problems show up in many physics applications, including normal mode analysis, quantum mechanics, and stability analyses.

1. Normal Modes

   • Vibrational analyses of structures or molecules reduce to $K \mathbf{x} = \omega^2 M \mathbf{x}$. The solutions are eigenpairs $(\omega^2, \mathbf{x})$.

   • Each eigenvector $\mathbf{x}$ describes a mode shape; the eigenvalue $\omega^2$ is the squared frequency.

2. Quantum Mechanics

   • The Schrodinger equation in matrix form $\hat{H}|E\rangle = E|E\rangle$ yields eigenvalues $E$ (energy levels) and eigenstates $|E\rangle$.

   • Collecting all eigenvectors (as columns) in a matrix diagonalizes $\hat{H}$ if it is Hermitian.

3. Stability Analysis

   • A linearized system near equilibrium produces a Jacobian $J$. Eigenvalues of $J$ show growth/decay rates of perturbations.

   • If real parts of eigenvalues are negative, the system is stable; if positive, it’s unstable.

Given an $n \times n$ matrix $A$, the eigenvalue problem seeks scalars $\lambda$ (eigenvalues) and nonzero vectors $\mathbf{v}$ (eigenvectors) such that:

• Eigenvalues can be real or complex.

• Eigenvectors identify the directions in which $A$ acts as a simple scale transformation ($\lambda$).

Power Method for the Dominant Eigenpair¶

The power method is a straightforward iterative technique for finding the eigenpair associated with the largest-magnitude eigenvalue:

1. Algorithm:

   • Begin with an initial vector $\mathbf{x}^{(0)}$.

   • Apply $A$ repeatedly and normalize:

     $$\begin{align} \mathbf{x}^{(k+1)} = \frac{A \mathbf{x}^{(k)}}{\|A \mathbf{x}^{(k)}\|}. \end{align}$$

     (23)

   • This converges to the eigenvector of the eigenvalue with the largest magnitude, assuming it’s distinct.

2. Limitations

   • Only computes one eigenvalue/eigenvector.

   • Convergence can be slow if other eigenvalues have magnitude close to the dominant one.

QR Algorithm for Dense Matrices¶

For dense matrices of moderate size, the $QR$ algorithm (or related methods) is the workhorse for computing all eigenvalues and optionally eigenvectors.

1. Idea:

   • Repeatedly factor $A$ as $QR$, where $Q$ is orthonormal and $R$ is upper triangular.

   • Set $A \leftarrow RQ$.

   • After enough iterations, $A$ becomes near-upper-triangular, revealing the eigenvalues on its diagonal.

2. Practical Approach

   • In Python, use numpy.linalg.eig or numpy.linalg.eigh (for Hermitian matrices).

   • These functions wrap optimized LAPACK routines implementing $QR$ or related transformations, e.g., divide-and-conquer, multiple relatively robust representations (MRRR).

Application: Coupled Harmonic Oscillators¶

Harmonic oscillators are a classic problem in physics, in this hands-on, we will:

1. Derive or reference the analytical solution for two coupled oscillators.

2. Numerically solve the same system (using an eigenvalue approach).

3. Generalize to $n$ (and even $n \times n$) coupled oscillators, visualizing the mode shapes.

Two Coupled Oscillators--Analytical Solution¶

Consider two masses $m$ connected by three springs (constant $k$), arranged in a line and connected to two walls:

|--k--[m]--k--[m]--k--|

If each mass can move only horizontally, the equations of motion form a $2 \times 2$ eigenvalue problem. Let:

• $x_1(t)$ be the horizontal displacement of Mass 1 from its equilibrium position.

• $x_2(t)$ be the horizontal displacement of Mass 2. We assume small oscillations, so Hooke’s law applies linearly.

• Mass 1 experiences:

  • A restoring force $-k \,x_1$ from the left wall spring.

  • A coupling force from the middle spring: if $x_2 > x_1$, that spring pulls Mass 1 to the right; if $x_1 > x_2$, it pulls Mass 1 to the left. The net contribution is $-k (x_1 - x_2)$. Summing forces (Newton’s second law) gives:

  $$\begin{align} m \ddot{x}_1 = -k x_1 - k (x_1 - x_2). \end{align}$$

  (24)

• Mass 2 experiences:

  • A restoring force $-k\,x_2$ from the right wall spring.

  • The coupling force from the middle spring: $-k(x_2 - x_1)$. Hence,

  $$\begin{align} m \ddot{x}_2 = -k x_2 - k (x_2 - x_1). \end{align}$$

  (25)

Rewrite each equation:

We can write $\mathbf{x} = \begin{pmatrix}x_1 \\ x_2\end{pmatrix}$ and express the system as:

where

Equivalently,

This is a second-order linear system describing small oscillations.

We look for solutions of the form

where $\mathbf{x}(0)$ is the initial condition and $\omega$ is the (angular) oscillation frequency.

Plugging into $m \ddot{\mathbf{x}} + K \mathbf{x} = 0$ gives:

Nontrivial solutions exist only if

which is the eigenvalue problem for $\omega^2$.

Explicitly, $K - m \omega^2 I$ is:

The determinant is $(2k - m \omega^2)^2 - k^2$. Setting this to zero results

Hence, we get two solutions for $\omega^2$:

1. $\omega_+^2$: taking the $+$ sign:

   $$\begin{align} 2k - m \omega_+^2 = k \quad \Longrightarrow \quad \omega_+^2 = \frac{k}{m} \quad \Longrightarrow \quad \omega_1 = \sqrt{\frac{k}{m}}. \end{align}$$

   (35)

2. $\omega_-^2$: taking the $-$ sign:

   $$\begin{align} 2k - m \omega_-^2 =-k \quad \Longrightarrow \quad \omega_-^2 = \frac{3k}{m} \quad \Longrightarrow \quad \omega_2 = \sqrt{\frac{3k}{m}}. \end{align}$$

   (36)

For each of the normal modes:

• Lower Frequency $\omega_+ = \sqrt{k/m}$: Plug $\omega_+^2 = k/m$ back into $(K - m\,\omega_+^2 I)\mathbf{X} = 0$. For instance,

  $$\begin{align} \begin{pmatrix} 2k - k & -k \\ -k & 2k - k \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} k & -k \\ -k & k \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \mathbf{0}. \end{align}$$

  (37)

  This implies $x_1 = x_2$. Physically, the in-phase mode has both masses moving together.

• Higher Frequency $\omega_- = \sqrt{3k/m}$:

  $$\begin{align} \begin{pmatrix} 2k - 3k & -k \\ -k & 2k - 3k \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} -k & -k \\ -k & -k \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \mathbf{0}. \end{align}$$

  (38)

  This yields $x_1 = -x_2$. Physically, the out-of-phase mode has the two masses moving in opposite directions.

We can compute the position of these coupled oscillators according to the analytical solutions.

Plot multiple frames:

Two Coupled Oscillators--Semi-Analytical/Numerical Solution¶

Instead of solving the coupled oscillator problem analytically, we can at least solve the eigenvalue part of the problem numerically.

Future Topics¶

Singular Value Decomposition (SVD)¶

SVD factorizes $A = U \Sigma V^T$. It reveals how a matrix stretches vectors in orthogonal directions. It works for non-square matrices, provides optimal low-rank approximations, is numerically stable, and is widely used for pseudoinverses, compression, and noise reduction.

Principal Component Analysis (PCA)¶

PCA applies SVD to centered data to find directions of maximum variance. It reduces dimensionality, highlights dominant modes, and simplifies noisy datasets, making it essential in physics, data modeling, and machine learning.

SVD and PCA are cores of dimensionality reduction. We may revisit them when we connect linear algebra to machine learning.