Applied Linear Algebra for Research

A skill for applying linear algebra to research computing, data analysis, and scientific modeling. Covers matrix decompositions, eigenvalue problems, least squares, dimensionality reduction, and practical implementation in NumPy/SciPy.

Essential Operations

Matrix Multiplication and Solving Systems

import numpy as np
from scipy import linalg


def solve_linear_system(A: np.ndarray, b: np.ndarray) -> dict:
    """
    Solve Ax = b and analyze the system.

    Args:
        A: Coefficient matrix (n x n)
        b: Right-hand side vector (n,)
    """
    n = A.shape[0]

    # Check condition number (sensitivity to perturbations)
    cond = np.linalg.cond(A)

    result = {
        "shape": A.shape,
        "rank": np.linalg.matrix_rank(A),
        "condition_number": cond,
        "well_conditioned": cond < 1e10,
    }

    if result["rank"] == n:
        x = np.linalg.solve(A, b)
        result["solution"] = x
        result["residual_norm"] = np.linalg.norm(A @ x - b)
    else:
        # Underdetermined or singular -- use least-squares
        x, residuals, rank, sv = np.linalg.lstsq(A, b, rcond=None)
        result["least_squares_solution"] = x
        result["note"] = "System is rank-deficient; least-squares solution returned"

    return result

Matrix Decompositions

LU Decomposition (Solving Multiple Systems)

def lu_factorization(A: np.ndarray) -> dict:
    """
    LU decomposition for efficiently solving Ax=b for multiple b.
    """
    P, L, U = linalg.lu(A)

    return {
        "P": P,  # Permutation matrix
        "L": L,  # Lower triangular
        "U": U,  # Upper triangular
        "usage": (
            "Once computed, solve for any new right-hand side b "
            "in O(n^2) instead of O(n^3). Use scipy.linalg.lu_solve()."
        )
    }

Singular Value Decomposition (SVD)

def svd_analysis(A: np.ndarray) -> dict:
    """
    SVD of matrix A = U S V^T and its applications.

    Args:
        A: Input matrix (m x n)
    """
    U, s, Vt = np.linalg.svd(A, full_matrices=False)

    return {
        "U_shape": U.shape,       # Left singular vectors (m x k)
        "singular_values": s,      # Sorted descending
        "Vt_shape": Vt.shape,     # Right singular vectors (k x n)
        "rank": np.sum(s > 1e-10),
        "condition_number": s[0] / s[-1] if s[-1] > 0 else float("inf"),
        "energy_ratio": np.cumsum(s ** 2) / np.sum(s ** 2),
        "applications": [
            "Low-rank approximation (truncated SVD)",
            "Principal Component Analysis (PCA)",
            "Pseudoinverse computation",
            "Latent Semantic Analysis (LSA) in text mining",
            "Image compression",
            "Noise reduction"
        ]
    }

Eigendecomposition

def eigen_analysis(A: np.ndarray) -> dict:
    """
    Eigenvalue decomposition of a square matrix.
    """
    eigenvalues, eigenvectors = np.linalg.eig(A)

    # Sort by magnitude
    idx = np.argsort(np.abs(eigenvalues))[::-1]

    return {
        "eigenvalues": eigenvalues[idx],
        "eigenvectors": eigenvectors[:, idx],
        "is_symmetric": np.allclose(A, A.T),
        "is_positive_definite": (
            np.all(np.real(eigenvalues) > 0)
            if np.allclose(A, A.T) else "N/A (not symmetric)"
        ),
        "spectral_radius": np.max(np.abs(eigenvalues)),
        "trace_check": (
            f"Sum of eigenvalues: {np.sum(eigenvalues):.4f}, "
            f"Trace of A: {np.trace(A):.4f}"
        )
    }

Research Applications

Principal Component Analysis

def pca_from_scratch(X: np.ndarray, n_components: int = 2) -> dict:
    """
    PCA using eigendecomposition of the covariance matrix.

    Args:
        X: Data matrix (n_samples x n_features), centered
        n_components: Number of principal components to retain
    """
    # Center the data
    X_centered = X - X.mean(axis=0)

    # Covariance matrix
    C = np.cov(X_centered, rowvar=False)

    # Eigendecomposition (symmetric matrix -> use eigh for stability)
    eigenvalues, eigenvectors = np.linalg.eigh(C)

    # Sort descending
    idx = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]

    # Select top components
    components = eigenvectors[:, :n_components]
    explained_variance = eigenvalues[:n_components]
    total_variance = eigenvalues.sum()

    # Project data
    X_projected = X_centered @ components

    return {
        "components": components,
        "explained_variance_ratio": explained_variance / total_variance,
        "cumulative_variance": np.cumsum(explained_variance) / total_variance,
        "projected_data": X_projected
    }

Least Squares Regression

def least_squares_fit(X: np.ndarray, y: np.ndarray) -> dict:
    """
    Solve the normal equations: beta = (X^T X)^{-1} X^T y
    """
    # Using the numerically stable QR decomposition
    Q, R = np.linalg.qr(X)
    beta = linalg.solve_triangular(R, Q.T @ y)

    y_hat = X @ beta
    residuals = y - y_hat

    return {
        "coefficients": beta,
        "r_squared": 1 - np.sum(residuals ** 2) / np.sum((y - y.mean()) ** 2),
        "residual_norm": np.linalg.norm(residuals),
        "method": "QR decomposition (more stable than normal equations)"
    }

Numerical Stability

Best Practices

1. Avoid explicitly computing matrix inverses:
   BAD:  x = np.linalg.inv(A) @ b
   GOOD: x = np.linalg.solve(A, b)

2. Use specialized routines for structured matrices:
   - Symmetric positive definite: Cholesky (linalg.cho_solve)
   - Sparse: scipy.sparse.linalg.spsolve
   - Banded: scipy.linalg.solve_banded

3. Check condition numbers before solving:
   - cond(A) > 10^10 suggests the solution may be unreliable
   - Consider regularization (Tikhonov/ridge) for ill-conditioned systems

4. Use appropriate precision:
   - float64 for most research computing
   - float32 for large-scale GPU computations (monitor for precision loss)

When working with very large matrices, leverage sparse matrix representations (scipy.sparse), iterative solvers (conjugate gradient, GMRES), and randomized algorithms (randomized SVD) to keep computation tractable.