Cholesky Factorization

Zero-allocation Cholesky decomposition and related linear algebra operations for symmetric positive definite (SPD) matrices, implemented in Julia.

Built for MATH 473 (Numerical Linear Algebra). The full writeup is in Writeup.pdf.

What's Implemented

• Cholesky factorization (cholesky) -- in-place, zero-allocation decomposition using the sub-matrix formulation with BLAS syr! for rank-1 updates
• Cholesky-Crout (cholesky_crout) -- vectorized reference implementation (not my own algorithm)
• Trefethen-Bau (cholesky_naive) -- direct implementation of Algorithm 23.1
• Determinant (det) -- via Cholesky factor; faster than stdlib for SPD matrices
• Solve (solve) -- forward/backward substitution on the Cholesky factor
• Inverse (inv) -- via triangular inverse of the Cholesky factor
• Least squares (solve_ls) -- solve via normal equations for non-SPD systems

Performance

The final algorithm achieves zero heap allocations by operating in-place using @views and BLAS routines. At 256x256 it runs within 2x of the LAPACK-backed stdlib, and at 16x16 within 1.4x.

Method	n=256	n=1024	n=4096
Stdlib (LAPACK)	155 us	4.7 ms	243 ms
Mine (vectorized)	316 us	15.2 ms	4.04 s
Crout	781 us	40.6 ms	2.53 s
Trefethen-Bau	9.53 ms	1.24 s	124.8 s

Usage

# In the Julia REPL, press ] to enter Pkg mode
pkg> activate .
pkg> resolve

# Run the test suite (1566 tests)
pkg> test

# Press backspace to return to Julia REPL
julia> include("src/ops.jl")

julia> using LinearAlgebra, BenchmarkTools

# Factorize a 256x256 SPD matrix
julia> A = Ops.rand_spd(256)
julia> L = Ops.cholesky(A)

# Verify zero allocations
julia> @benchmark Ops.cholesky(A) setup=(A=Ops.rand_spd(256)) evals=1

# Solve Ax = b
julia> b = rand(256)
julia> x = Ops.solve(A, b)

# Determinant
julia> Ops.det(A)

Benchmarks

julia> include("bench/benchmark.jl")
julia> tune!(suite)
julia> results = run(suite, verbose=true, seconds=5)

Project Structure

src/ops.jl          -- All algorithms (Cholesky, solve, det, inv)
src/assignment1.jl  -- Module wrapper
test/runtests.jl    -- Test suite
bench/benchmark.jl  -- BenchmarkTools suite
Writeup.pdf         -- Formal writeup with analysis
figures/            -- Performance plots

Testing

The test suite validates all operations against Julia's stdlib (LinearAlgebra.cholesky). Tests run across matrix sizes 3 to 1024, including both well-conditioned and ill-conditioned random SPD matrices. All results must satisfy ||x - x_hat||_2 < 10^-8.