Module 1b: Statistical Linear Algebra and Numerical Methods

Quick Navigation Guide¶

📚 Choose Your Learning Path¶

🎯 Navigation by Project Needs¶

Learning Objectives¶

By the end of this module, you will be able to:

• Determine whether a matrix is positive definite using multiple methods 🔴

• Construct and interpret covariance matrices from data 🔴

• Apply Cholesky decomposition for sampling and solving 🟡

• Debug numerical failures in statistical computations 🔴

• Connect linear algebra to MCMC, Gaussian Processes, and neural networks 🟡

• Implement robust algorithms that handle finite precision 🟢

Prerequisites from Module 0a¶

Module Overview¶

Building on Module 0a’s foundations, this module explores how linear algebra enables statistical computation and machine learning. We cover positive definite matrices that ensure physical validity (Project 4), covariance structures that encode correlations (Project 5), and advanced numerical methods that keep statistical computations stable despite finite precision (all projects). The module culminates by showing how these concepts unite classical physics with modern machine learning (Final Project).

Part 1: Positive Definite Matrices and Statistical Foundations¶

Priority: 🟡 Important - Essential for Projects 4-5

1.1 Quadratic Forms and Energy {#part-5-positive-definite}¶

A quadratic form is:

These appear as energy expressions throughout physics:

Physical Example: The kinetic energy of a rotating rigid body

where $\vec{\omega}$ is angular velocity and $I$ is the moment of inertia tensor. This quadratic form is always positive (energy can’t be negative), making $I$ positive definite. The eigenvectors of $I$ are the principal axes - spin around these and the object doesn’t wobble!

1.2 Positive Definiteness: Ensuring Physical Reality {#positive-definiteness}¶

A symmetric matrix is positive definite if its quadratic form is always positive. This property ensures physical validity—energies are positive, distances are non-negative, and covariance matrices represent valid uncertainties.

Four Equivalent Tests for Positive Definiteness:

1. ✅ All eigenvalues > 0

2. ✅ All leading principal minors > 0

3. ✅ Has Cholesky decomposition $A = LL^T$ (requires strict positive definiteness)

4. ✅ Can write as $A = B^T B$ for some invertible matrix $B$

Quick check: Positive definite matrices always have positive diagonal elements. Proof: For standard basis vector $\vec{e}_i$, we have $\vec{e}_i^T A \vec{e}_i = A_{ii} > 0$.

Important Distinction:

• Positive definite (all eigenvalues > 0): Standard Cholesky works

• Positive semi-definite (eigenvalues ≥ 0): May have rank deficiency (rank < dimension), requiring modified Cholesky or eigendecomposition

• Why it matters: Zero eigenvalues mean the matrix maps some directions to zero—information is lost

1.3 Covariance Matrices: The Bridge to Statistics {#covariance-matrices}¶

Priority: 🟡 Important - Foundation for Projects 4-5

Before diving into matrices, let’s understand what covariance actually tells us:

Intuitive Understanding:

• Positive covariance: Variables tend to increase together (mass and luminosity in stars)

• Negative covariance: One increases as the other decreases (stellar temperature and radius for giants)

• Zero covariance: No linear relationship (doesn’t mean independent!)

Mathematical Definition:

This measures the average product of deviations from the mean. When both variables are above (or below) their means together, the products are positive, giving positive covariance.

Relationship to Correlation:

Correlation is just normalized covariance! It ranges from -1 to +1, making it easier to interpret than raw covariance.

The Covariance Matrix:

For multiple variables $X_1, ..., X_n$, we organize all pairwise covariances into a matrix:

Structure and Properties:

• Diagonal elements: $\Sigma_{ii} = \text{Var}(X_i)$ (variances are “self-covariance”)

• Off-diagonal elements: $\Sigma_{ij}$ (covariances between different variables)

• Always symmetric: $\Sigma_{ij} = \Sigma_{ji}$ (order doesn’t matter)

• Always positive semi-definite: All eigenvalues ≥ 0 (why? see below)

Why Positive Semi-Definite?

For any linear combination of variables $Y = a^T X$:

Variance can’t be negative! This forces $\Sigma$ to be positive semi-definite. If it has negative eigenvalues, you could construct a linear combination with negative variance—physically impossible!

1.4 The Multivariate Gaussian Distribution {#multivariate-gaussian}¶

Priority: 🟡 Important - Core of Projects 4-5

The multivariate Gaussian extends the familiar bell curve to multiple dimensions. In 1D, you know the Gaussian as the bell-shaped $e^{-x^2}$ curve. In higher dimensions, it becomes an ellipsoidal cloud of probability.

The Formula Decoded:

Let’s understand each piece intuitively:

The Exponent: $-\frac{1}{2}(\vec{x}-\vec{\mu})^T\Sigma^{-1}(\vec{x}-\vec{\mu})$

• This is the Mahalanobis distance squared

• In 1D with variance $\sigma^2$: reduces to $-(x-\mu)^2/(2\sigma^2)$

• Measures “how many standard deviations away” but accounts for correlations

• If variables are correlated, being far in one direction might be “cheaper” than another

The Normalization: $(2\pi)^{n/2}|\Sigma|^{1/2}$

• Ensures total probability integrates to 1

• $|\Sigma|$ is the determinant—the “volume” of the uncertainty ellipsoid

• Larger determinant = more spread out = lower peak height (volume conserved)

Geometric Picture:

• Mean $\vec{\mu}$: Center of the probability cloud

• Covariance $\Sigma$: Shape and orientation of the ellipsoid

  • Eigenvalues: Lengths of ellipsoid axes

  • Eigenvectors: Directions of axes

• Contours of constant probability: Ellipsoids defined by $(\vec{x}-\vec{\mu})^T\Sigma^{-1}(\vec{x}-\vec{\mu}) = c$

1.5 Cholesky Decomposition: The Matrix Square Root¶

Priority: 🟡 Important - Essential for Project 5

Every positive definite matrix can be factored as:

where $L$ is lower triangular with positive diagonal entries.

Uniqueness Property: The Cholesky decomposition is unique — for any positive definite matrix $A$, there exists exactly one lower triangular matrix $L$ with positive diagonal entries such that $A = LL^T$.

This uniqueness is crucial for numerical methods: no matter which algorithm you use (whether the classical Cholesky-Banachiewicz or the Cholesky-Crout variant), you’ll always get the same $L$ matrix (up to numerical precision). This deterministic property makes Cholesky decomposition ideal for reproducible scientific computing.

Why uniqueness matters: When you use Cholesky in your MCMC proposals (Project 4) or GP implementations (Project 5), you’re guaranteed consistent results across different machines and libraries. This isn’t true for eigendecomposition, where eigenvectors can differ by sign or be arbitrarily chosen when eigenvalues are repeated.

Intuitive Understanding: Think of Cholesky as finding the “square root” of a matrix

• Just as $9 = 3 \times 3$, we have $A = L \times L^T$

• $L$ transforms uncorrelated unit variance variables into correlated variables with covariance $A$

Why Lower Triangular? The triangular structure means each variable depends only on previous ones:

• Variable 1: Independent

• Variable 2: Depends on variable 1

• Variable 3: Depends on variables 1 and 2

• And so on...

This sequential dependency structure makes computation efficient!

Two Key Applications:

1. Solving Linear Systems (faster than computing inverse):

# Problem: Solve Ax = b where A is positive definite

# Slow way: x = inv(A) @ b (NEVER DO THIS!)
# - Computing inverse: O(n³) operations
# - Numerically unstable
# - Loses sparsity structure

# Fast way using Cholesky:
L = np.linalg.cholesky(A)      # A = LL^T
# Now Ax = b becomes LL^T x = b

# Step 1: Solve Ly = b (forward substitution)
y = np.linalg.solve(L, b)

# Step 2: Solve L^T x = y (back substitution)
x = np.linalg.solve(L.T, y)

# Why faster? Triangular systems solve in O(n²) not O(n³)!

2. Generating Correlated Random Variables (Projects 4-5):

# Generate samples from N(mu, Sigma)

# The magic: If z ~ N(0, I), then Lz ~ N(0, LL^T) = N(0, Sigma)
L = np.linalg.cholesky(Sigma)
z = np.random.randn(n)          # Standard normal (uncorrelated)
x = mu + L @ z                  # Has mean mu, covariance Sigma

# Why does this work?
# - z has covariance I (independent unit variance)
# - Lz has covariance L·I·L^T = LL^T = Sigma
# - Adding mu shifts the mean without changing covariance

Progressive Problems: Positive Definiteness¶

Part 2: Advanced Topics for Your Projects¶

Priority: 🟢 Enrichment - Read as needed for specific projects

2.1 Singular Value Decomposition - The Swiss Army Knife {#svd-swiss-army}¶

Every matrix has a singular value decomposition:

where:

• $U$: Left singular vectors (orthonormal output directions)

• $\Sigma$: Diagonal matrix of singular values (stretching factors)

• $V^T$: Right singular vectors (orthonormal input directions)

SVD for Non-Square Matrices: For an $m \times n$ matrix $A$:

• $U$ is $m \times m$ (square, orthogonal)

• $\Sigma$ is $m \times n$ (rectangular!)

• $V^T$ is $n \times n$ (square, orthogonal)

Only $\min(m,n)$ singular values exist; the rest of $\Sigma$ is padded with zeros. This is crucial for understanding rank and dimensionality reduction.

Geometric Intuition: Any matrix transformation can be broken into three steps:

1. Rotate (by $V^T$): Align input to principal axes

2. Stretch (by $\Sigma$): Scale along each axis by σᵢ

3. Rotate (by $U$): Align to output space

This means ANY linear transformation—no matter how complex—is just rotate-stretch-rotate! This decomposition is unique (up to sign ambiguities) and always exists.

Understanding Rank Through SVD: The rank of a matrix equals the number of non-zero singular values. This tells you the true dimensionality:

• Full rank: All $σᵢ > 0$, no information lost

• Rank deficient: Some $σᵢ = 0$, transformation loses dimensions

• Numerical rank: Count $σᵢ >$ tolerance (e.g., 10⁻¹⁰) for finite precision

The Pseudoinverse and Least Squares:

The pseudoinverse $A^+ = V\Sigma^+ U^T$ where $\Sigma^+$ inverts non-zero singular values:

• For overdetermined systems (more equations than unknowns): Gives least-squares solution

• For underdetermined systems (more unknowns than equations): Gives minimum-norm solution

• Robust via truncation: Ignore singular values < threshold

# Computing pseudoinverse via SVD
U, s, Vt = np.linalg.svd(A, full_matrices=False)

# Threshold small singular values (regularization)
threshold = 1e-10 * s.max()
s_inv = np.where(s > threshold, 1/s, 0)

# Construct pseudoinverse
A_pseudo = Vt.T @ np.diag(s_inv) @ U.T

# Now x = A_pseudo @ b gives:
# - Least-squares solution if overdetermined
# - Minimum-norm solution if underdetermined

Why SVD Beats Eigendecomposition:

The Deep Connection: For any matrix $A$:

• $A^T A$ has eigenvalues $\lambda_i = \sigma_i^2$ and eigenvectors = columns of $V$

• $AA^T$ has eigenvalues $\lambda_i = \sigma_i^2$ and eigenvectors = columns of $U$

• This is why SVD always exists—symmetric matrices always have eigendecompositions!

2.2 Block Matrices and the Schur Complement¶

Large systems often have natural block structure:

Physical Motivation: Different parts of your system interact differently.

Concrete Example: Triple Star System Consider a close binary with a distant third star: $$\begin{pmatrix} F_{11} & F_{12} & F_{13} \ F_{21} & F_{22} & F_{23} \ F_{31} & F_{32} & F_{33} \end{pmatrix} =

The 2×2 upper-left block represents strong binary interactions, while off-diagonal blocks show weak coupling to the distant star.

The Schur Complement in Action:

Let’s solve $\begin{pmatrix} A & B \\ C & D \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} f \\ g \end{pmatrix}$

Instead of solving the full system, we can eliminate $y$:

1. From bottom equation: $y = D^{-1}(g - Cx)$

2. Substitute into top: $Ax + BD^{-1}(g - Cx) = f$

3. Rearrange: $(A - BD^{-1}C)x = f - BD^{-1}g$

The matrix $S = A - BD^{-1}C$ is the Schur complement—it’s the “effective $A$” after accounting for $y$’s influence through $D$.

Physical Intuition: Think of the Schur complement as “marginalizing out” variables. It’s like projecting 3D motion onto 2D while accounting for the 3rd dimension’s influence. When you eliminate variables, their effect doesn’t disappear—it gets encoded in the effective interactions between remaining variables.

Simple 2×2 Example:

Schur complement of $D$: $S = 3 - (1)(1/4)(2) = 3 - 0.5 = 2.5$

This 2.5 is the “effective resistance” in the first variable after the second adjusts optimally.

Why It Matters for Project 5 (GPs): When you add a new observation to a Gaussian Process:

• Old observations = block $A$

• Covariance with new point = blocks $B, C$

• New point variance = block $D$

• Prediction uses the Schur complement to efficiently update without re-inverting everything!

2.3 The Jacobian Matrix: Local Linear Approximation {#jacobian-matrix}¶

For vector function $\vec{f}: \mathbb{R}^n \to \mathbb{R}^m$:

Intuitive Understanding: Near any point, nonlinear functions look linear. The Jacobian is that linear approximation:

Concrete 2D Example: For the transformation $f(r, \theta) = (r\cos\theta, r\sin\theta)$ (polar to Cartesian):

The determinant $|J| = r$ tells you area scaling—this is why polar area elements are $r\,dr\,d\theta$!

In Project 2 (N-body Stability): Near an equilibrium configuration $\vec{s}_0$, perturbations $\delta\vec{s}$ evolve as:

The eigenvalues of $J$ determine the fate:

• All Re(λ) < 0: Stable (perturbations decay)

• Any Re(λ) > 0: Unstable (perturbations grow)

• Re(λ) = 0, Im(λ) ≠ 0: Neutral oscillations

Example: For a star orbiting at Lagrange point L4:

• Eigenvalues have small negative real parts → weakly stable

• Large imaginary parts → oscillates if perturbed

• This is why Trojan asteroids librate around L4/L5!

2.4 Matrix Exponentials: Solving Linear Evolution {#matrix-exponentials}¶

The matrix exponential solves any linear ODE system:

Three Ways to Understand $e^{At}$:

1. Series Definition (like scalar exponential):

2. Through Eigenvalues (if $A$ is diagonalizable): If $A = PDP^{-1}$ where $D$ is diagonal with eigenvalues:

3. Physical Interpretation: $e^{At}$ is the propagator—it evolves the system from time 0 to time $t$.

Practical Computation:

# Never use the series directly—it's inefficient and can be inaccurate
# Instead, use specialized algorithms:

import scipy.linalg as la

A = np.array([[0, 1], [-1, -0.1]])  # Damped oscillator
t = 5.0

# Method 1: Direct computation (best for dense matrices)
expAt = la.expm(A * t)

# Method 2: Via eigendecomposition (good for understanding)
eigvals, eigvecs = la.eig(A)
D = np.diag(np.exp(eigvals * t))
expAt_eig = eigvecs @ D @ la.inv(eigvecs)

# Method 3: For solving ODEs, use ODE solvers instead!
# They're more efficient than computing e^{At} explicitly

Example: Damped Harmonic Oscillator: For $\ddot{x} + \gamma\dot{x} + \omega_0^2 x = 0$, rewrite as first-order system:

The eigenvalues $\lambda = \frac{-\gamma \pm \sqrt{\gamma^2 - 4\omega_0^2}}{2}$ determine behavior:

• Underdamped ($\gamma < 2\omega_0$): Complex λ → oscillatory decay

• Critically damped ($\gamma = 2\omega_0$): Repeated real λ → fastest decay

• Overdamped ($\gamma > 2\omega_0$): Two real λ → slow decay

2.5 Matrix Norms: How Big is a Matrix?¶

Priority: 🟢 Enrichment - Crucial for understanding stability and convergence

Matrix norms measure “size” but different norms capture different aspects:

Frobenius Norm (total “energy”):

• Like treating matrix as a long vector

• All elements contribute equally

• Easy to compute but ignores structure

Spectral Norm (worst-case amplification):

• Maximum stretch factor for any unit vector

• Equals largest singular value

• Determines stability and convergence

Why Spectral Norm Matters:

1. Stability: System $x_{n+1} = Ax_n$ is stable iff $||A||_2 < 1$

2. Error amplification: Input error $\epsilon$ → output error ≤ $||A||_2 \cdot \epsilon$

3. Condition number: $\kappa = ||A||_2 \cdot ||A^{-1}||_2 = \sigma_{\max}/\sigma_{\min}$

4. Convergence rate: Iterations converge like $(||A||_2)^n$

Example: Why Deep Networks Are Hard to Train: Consider a 10-layer network where each layer multiplies by matrix $W$:

• If $||W||_2 = 1.1$: Gradients grow as $(1.1)^{10} \approx 2.6×$

• If $||W||_2 = 0.9$: Gradients shrink as $(0.9)^{10} \approx 0.35×$

• Need $||W||_2 \approx 1$ for stable training!

This is why techniques like batch normalization and careful initialization are crucial.

2.6: Numerical Implementation Examples¶

Part 3: Numerical Reality - When Mathematics Meets Silicon¶

Priority: 🔴 Essential - Critical for debugging all projects

7.1 Condition Numbers Revisited: Statistical Context {#condition-numbers}¶

Building on Module 0a’s introduction to condition numbers, let’s explore their specific implications for statistical computations:

Condition Numbers in Statistical Methods:

For covariance matrices, the condition number has special meaning:

This ratio tells you:

• Near 1: Variables have similar scales and low correlation

• Large: Either variables have vastly different scales OR near-perfect correlation

• Infinite: Perfect correlation (singular covariance)

MCMC Implications (Project 4):

# Poorly conditioned proposal covariance slows convergence
Sigma_bad = np.array([[1, 0.999], [0.999, 1]])  # κ ≈ 2000
# MCMC needs ~κ steps to explore the distribution

# Well-conditioned after reparameterization
L = np.linalg.cholesky(Sigma_bad)
# Sample in transformed space where κ = 1, then transform back

GP Implications (Project 5):

# Kernel matrices become ill-conditioned when:
# 1. Data points are very close (duplicate measurements)
# 2. Length scale is too small (overfitting)
# 3. Noise term is too small (numerical instability)

# Check your GP's condition number:
kappa = np.linalg.cond(K)
if kappa > 1e12:
    print("Warning: GP may give nonsense predictions!")
    # Solutions:
    # - Increase noise term σ²
    # - Remove duplicate points
    # - Adjust length scale

7.2 When Statistical Computations Fail¶

Understanding failure modes specific to statistical methods helps you recognize and fix problems:

Covariance Matrix Not Positive Semi-Definite:

# Computing sample covariance with insufficient data
n_samples = 3
n_features = 5  # More features than samples!

X = np.random.randn(n_samples, n_features)
Sigma = np.cov(X.T)  # Shape: (5, 5)

# This is rank-deficient: rank ≤ min(n_samples-1, n_features) = 2
print(f"Rank: {np.linalg.matrix_rank(Sigma)}")  # 2, not 5!

# Cholesky will fail:
try:
    L = np.linalg.cholesky(Sigma)
except np.linalg.LinAlgError:
    print("Covariance is singular!")

    # Fix: Regularization (add prior belief about variance)
    Sigma_reg = Sigma + 0.01 * np.eye(n_features)
    L = np.linalg.cholesky(Sigma_reg)  # Now works

Numerical Underflow in Likelihoods:

# Computing likelihood of many data points
log_likelihoods = []
for x in data:
    # BAD: Probability underflows to zero
    p = multivariate_normal.pdf(x, mu, Sigma)

    # GOOD: Work in log space
    log_p = multivariate_normal.logpdf(x, mu, Sigma)
    log_likelihoods.append(log_p)

# Sum log-likelihoods instead of multiplying probabilities
total_log_likelihood = np.sum(log_likelihoods)

Cholesky Decomposition Failures:

def safe_cholesky(A, max_tries=5):
    """
    Robust Cholesky with automatic regularization.
    """
    jitter = 1e-6

    for i in range(max_tries):
        try:
            L = np.linalg.cholesky(A + jitter * np.eye(len(A)))
            return L
        except np.linalg.LinAlgError:
            if i == max_tries - 1:
                # Last resort: eigenvalue thresholding
                eigvals, eigvecs = np.linalg.eigh(A)
                eigvals = np.maximum(eigvals, jitter)
                A_fixed = eigvecs @ np.diag(eigvals) @ eigvecs.T
                return np.linalg.cholesky(A_fixed)
            else:
                jitter *= 10  # Increase regularization

    return None

7.3 Statistical-Specific Numerical Guidelines¶

Part 4: Bridge to Machine Learning¶

Priority: 🟡 Important - Helps see the big picture

8.1 From Classical to Statistical to Learning¶

Linear algebra provides the mathematical continuity across your entire journey. Watch how the same mathematical objects evolve in meaning as you progress through different computational paradigms:

The Evolution of Mathematical Objects:

This isn’t coincidence—it’s the deep unity of mathematics. The same linear algebra that predicts planetary orbits also powers Google’s search algorithm and ChatGPT’s language understanding.

8.2 Matrix Calculus for Neural Networks {#matrix-calculus-nn}¶

Priority: 🟡 Important for Final Project

Your final project requires understanding how gradients flow through matrix operations. Here’s the essential matrix calculus:

Basic Matrix Derivatives:

For scalar function $f$ of matrix $W$:

Key Results You’ll Need:

1. Linear layer: $f(W) = W\vec{x}$

   $$\frac{\partial}{\partial W}(W\vec{x}) = \vec{x}^T$$

   (25)

2. Quadratic form: $f(W) = \vec{x}^T W \vec{x}$

   $$\frac{\partial}{\partial W}(\vec{x}^T W \vec{x}) = \vec{x}\vec{x}^T$$

   (26)

3. Frobenius norm: $f(W) = ||W||_F^2 = \text{trace}(W^T W)$

   $$\frac{\partial}{\partial W}||W||_F^2 = 2W$$

   (27)

Chain Rule for Matrices:

For composition $f(g(W))$:

This is what happens during backpropagation—gradients flow backward through the network via chain rule!

Example: Simple Two-Layer Network:

# Forward pass:
# Layer 1: z₁ = W₁x + b₁
# Activation: a₁ = relu(z₁)
# Layer 2: z₂ = W₂a₁ + b₂
# Loss: L = ||z₂ - y||²

# Backward pass (what JAX computes automatically):
# ∂L/∂z₂ = 2(z₂ - y)
# ∂L/∂W₂ = (∂L/∂z₂) @ a₁ᵀ
# ∂L/∂a₁ = W₂ᵀ @ (∂L/∂z₂)
# ∂L/∂z₁ = (∂L/∂a₁) ⊙ relu'(z₁)  # Element-wise product
# ∂L/∂W₁ = (∂L/∂z₁) @ xᵀ

# The pattern: gradients with respect to weights = (upstream gradient) @ (input)ᵀ

8.3 Connecting Everything: The Big Picture¶

The techniques in this module power every major simulation and discovery in modern astrophysics:

Your Learning Journey:

1. Classical (Projects 1-3): Deterministic systems, exact solutions

2. Statistical (Projects 4-5): Uncertainty quantification, probabilistic inference

3. Learning (Final Project): Pattern discovery, function approximation

The Mathematical Thread:

• Same linear algebra throughout

• Increasing complexity of applications

• Growing appreciation for numerical stability

• Deepening understanding of connections

Why This Matters: Every breakthrough in computational astrophysics builds on these foundations:

• Gaia: 1 billion stars → massive covariance matrices

• LIGO: Gravitational waves → matched filtering in high dimensions

• JWST: Spectroscopy → dimensionality reduction via SVD

• Vera Rubin: 20 TB/night → statistical methods at scale

Numerical Recipes Quick Reference¶

Glossary of Terms¶

All terms introduced in this module, arranged alphabetically:

• back substitution: Solving upper triangular systems from bottom to top, used in the second phase of solving systems via Cholesky decomposition

• block matrix: Matrix partitioned into submatrices, often reflecting natural system structure (e.g., strong/weak interactions in multi-body systems)

• Cholesky decomposition: Factorization of a positive definite matrix $A = LL^T$ where $L$ is lower triangular with positive diagonal entries; geometrically finds the “square root” of a positive definite matrix

• condition number: For statistical context: $\kappa(A) = \sigma_{\max}/\sigma_{\min}$, ratio of largest to smallest singular value; measures how much errors amplify in statistical computations

• covariance: Measure of how two variables change together; positive means they increase together, negative means one increases as the other decreases, zero means no linear relationship

• covariance matrix: Matrix $\Sigma$ containing all pairwise covariances between random variables; $\Sigma_{ij} = \text{Cov}(X_i, X_j)$; encodes all linear relationships in your data

• forward substitution: Solving lower triangular systems from top to bottom, used in the first phase of solving systems via Cholesky decomposition

• Frobenius norm: Matrix norm defined as $||A||_F = \sqrt{\sum_{i,j} |a_{ij}|^2}$; treats matrix like a vector and sums all squared elements

• Jacobian: Matrix of all first-order partial derivatives $J_{ij} = \partial f_i/\partial x_j$; represents the best linear approximation of a vector function at a point

• Mahalanobis distance: Scale-invariant distance that accounts for correlations: $d = \sqrt{(\vec{x}-\vec{\mu})^T\Sigma^{-1}(\vec{x}-\vec{\mu})}$; measures “how many standard deviations away” in correlated space

• matrix exponential: The operator $e^{At}$ that propagates linear systems forward in time; solution to $\dot{x} = Ax$ is $x(t) = e^{At}x(0)$

• multivariate Gaussian: Multi-dimensional generalization of the bell curve, defined by mean vector $\mu$ and covariance matrix $\Sigma$; probability density involves $\exp(-\frac{1}{2}(\vec{x}-\vec{\mu})^T\Sigma^{-1}(\vec{x}-\vec{\mu}))$

• positive definite: Property of a symmetric matrix where $\vec{x}^T A \vec{x} > 0$ for all $\vec{x} \neq \vec{0}$; ensures energies are positive and covariances are valid

• positive semi-definite: Property of a symmetric matrix where $\vec{x}^T A \vec{x} \geq 0$ for all $\vec{x}$; allows for zero eigenvalues (rank deficiency)

• precision matrix: Inverse of the covariance matrix: $\Sigma^{-1}$; encodes conditional independence structure (zero entries indicate conditional independence)

• pseudoinverse: Generalized inverse $A^+ = V\Sigma^+ U^T$ obtained via SVD; gives least-squares solution for overdetermined systems and minimum-norm solution for underdetermined systems

• quadratic form: Expression $Q(\vec{x}) = \vec{x}^T A \vec{x}$ where $A$ is symmetric; represents energies, distances, and other positive quantities in physics

• regularization: Adding small positive values (often to diagonal elements) to improve numerical stability or prevent singularities; balances stability with bias

• Schur complement: The “effective” matrix after eliminating some variables: $S = A - BD^{-1}C$; represents how eliminated variables influence remaining ones

• singular values: Non-negative values $\sigma_1 \geq \sigma_2 \geq ... \geq 0$ from SVD; measure importance/variance explained by each component

• spectral norm: Matrix norm defined as $||A||_2 = \sigma_{\max}$; equals maximum amplification factor for any unit vector

• SVD (Singular Value Decomposition): Universal matrix factorization $A = U\Sigma V^T$ that reveals the fundamental action of any matrix as rotate-stretch-rotate

Main Takeaways¶

Essential NumPy/SciPy Reference for Module 0b¶

Priority: 🔴 Essential - Keep this open while coding Projects 4-6

Quick Reference Card¶

With these foundations, you’re equipped to tackle the statistical and machine learning challenges ahead!