6.3.2 Linear Algebra Computations

Overview

Linear algebra computations form the computational backbone of modern science, engineering, and machine learning. While theoretical linear algebra focuses on abstract vector spaces and structural properties, computational linear algebra emphasizes algorithmic efficiency, numerical stability, and practical implementation for finite-dimensional spaces.

This section covers the fundamental operations, core decomposition algorithms, complexity considerations, and modern computational frameworks used to solve linear systems, analyze matrices, and perform high-dimensional transformations.

Fundamental Operations

All advanced linear algebra computations reduce to a set of primitive operations. Understanding their mathematical definitions and computational costs is essential for algorithm design.

Vector Operations

Given vectors \(\mathbf{u}, \mathbf{v} \in \mathbb{R}^n\) and scalar \(\alpha \in \mathbb{R}\):

\mathbf{u} + \mathbf{v} = \begin{bmatrix} u_1 + v_1 \\ \vdots \\ u_n + v_n \end{bmatrix}, \quad \alpha\mathbf{u} = \begin{bmatrix} \alpha u_1 \\ \vdots \\ \alpha u_n \end{bmatrix}, \quad \mathbf{u} \cdot \mathbf{v} = \sum_{i=1}^{n} u_i v_i

The dot product has a geometric interpretation: \(\mathbf{u} \cdot \mathbf{v} = \|\mathbf{u}\| \|\mathbf{v}\| \cos\theta\), where \(\theta\) is the angle between vectors. Computationally, this requires \(O(n)\) operations.

Matrix Operations

For matrices \(A \in \mathbb{R}^{m \times n}\) and \(B \in \mathbb{R}^{n \times p}\):

(AB)_{ij} = \sum_{k=1}^{n} A_{ik} B_{kj}, \quad A^T \in \mathbb{R}^{n \times m} \text{ where } (A^T)_{ij} = A_{ji}

Matrix multiplication is associative but not commutative. The naive algorithm requires \(O(mnp)\) floating-point operations (FLOPs). For square \(n \times n\) matrices, this becomes \(O(n^3)\).

💡 Computational Note

Modern libraries (BLAS/LAPACK) optimize matrix multiplication using cache-blocking, vectorization, and multi-threading. The Strassen algorithm reduces complexity to \(O(n^{\log_2 7}) \approx O(n^{2.807})\), though constant factors limit practical use for small \(n\).

Core Algorithms & Decompositions

Direct computation of inverses is discouraged in numerical practice. Instead, matrix decompositions factorize \(A\) into structured components that simplify solving systems, computing eigenvalues, or reducing dimensionality.

Decomposition	Form	Complexity	Primary Use
LU	\(A = LU\)	\(O(n^3)\)	Solving linear systems \(Ax=b\)
QR	\(A = QR\)	\(O(n^3)\)	Least squares, eigenvalue algorithms
Cholesky	\(A = LL^T\)	\(O(n^3/3)\)	Positive definite systems
SVD	\(A = U\Sigma V^T\)	\(O(n^3)\)	Dimensionality reduction, pseudoinverse
Eig	\(A = V\Lambda V^{-1}\)	\(O(n^3)\)	Dynamical systems, PCA, stability

LU Decomposition & Gaussian Elimination

LU factorization expresses \(A\) as a product of a lower triangular matrix \(L\) and upper triangular matrix \(U\). It generalizes Gaussian elimination into a reusable factorization step.

\begin{bmatrix} 2 & 1 & 1 \\ 4 & -6 & 0 \\ -2 & 2 & 0 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & 0 & 1 \end{bmatrix} \begin{bmatrix} 2 & 1 & 1 \\ 0 & -8 & -2 \\ 0 & 0 & -2 \end{bmatrix}

Once factorized, solving \(Ax=b\) reduces to forward substitution for \(Ly=b\) and backward substitution for \(Ux=y\), both \(O(n^2)\). Partial pivoting \(PA=LU\) is standard for numerical stability.

Singular Value Decomposition (SVD)

SVD applies to any \(m \times n\) matrix and reveals its geometric structure. The singular values \(\sigma_1 \geq \sigma_2 \geq \dots \geq 0\) measure the matrix's action along orthogonal directions.

A = U \Sigma V^T, \quad U^TU = I_m, \; V^TV = I_n, \; \Sigma = \text{diag}(\sigma_1, \dots, \sigma_r)

Truncated SVD (keeping top \(k\) singular values) enables low-rank approximation with minimal error (Eckart-Young-Mirsky theorem). This is foundational for PCA, recommendation systems, and noise reduction.

Numerical Stability & Precision

Real-world computations use finite-precision arithmetic (IEEE 754). Rounding errors accumulate, and ill-conditioned problems amplify noise.

Condition Numbers

The condition number of a matrix \(A\) measures sensitivity to perturbations:

\kappa(A) = \|A\| \|A^{-1}\| = \frac{\sigma_{\max}}{\sigma_{\min}}

If \(\kappa(A) \approx 10^d\), expect loss of \(d\) decimal digits of accuracy. Matrices with \(\kappa(A) \gg 1\) are ill-conditioned and require regularization or higher precision.

Python / NumPy

import numpy as np
from numpy.linalg import cond, solve, lstsq

# Check conditioning before solving
A = np.array([[1, 1.0001], [1.0001, 1]])
b = np.array([2, 2.0001])

print(f"Condition number: {cond(A)}")
x = solve(A, b)
print(f"Solution: {x}")

# Least squares for over-determined systems
x_lstsq, residuals, rank, sv = lstsq(A, b, rcond=None)

Modern Computational Approaches

Contemporary linear algebra computation leverages hardware acceleration and algorithmic optimizations:

GPU Acceleration: CUDA/MAGMA libraries parallelize matrix kernels across thousands of cores, achieving 10-50x speedups for large \(n\).
Iterative Solvers: For sparse systems, conjugate gradient (CG) and GMRES avoid \(O(n^3)\) factorization, converging in \(O(kn)\) iterations with preconditioning.
Randomized Linear Algebra: Sketching and randomized SVD approximate large matrices in \(O(n^2 \log k)\) time, crucial for big data pipelines.
Auto-differentiation Integration: Frameworks like PyTorch and JAX embed linear algebra ops in computational graphs for gradient-based optimization.

Applications

Linear algebra computations power virtually every data-driven discipline:

Machine Learning: Neural network forward/backward passes, PCA, kernel methods, attention mechanisms
Computer Graphics: Homogeneous transformations, ray tracing, mesh processing
Physics & Engineering: FEM/FDM solvers, quantum state evolution, structural analysis
Signal Processing: Fourier transforms, filtering, compression (JPEG, MP3)
Finance: Portfolio optimization, risk modeling, option pricing PDEs