NLA Engine

Numerical Linear Algebra engine with implementations in C++, Python, and R. Everything is built from scratch — no wrappers around LAPACK or scipy internals.

What's in here

Decompositions

• LU with partial pivoting (Doolittle variant)
• QR via Householder reflections
• Cholesky for symmetric positive definite matrices
• SVD via Jacobi eigendecomposition

Eigenvalue Solvers

• Power iteration (dominant eigenvalue)
• Inverse iteration (eigenvalue near a given shift)
• QR algorithm with Wilkinson shifts (all eigenvalues, symmetric case)
• Hessenberg reduction preprocessing

Iterative Solvers

• Conjugate Gradient (CG) for SPD systems
• GMRES with restarts for general systems
• BiCGSTAB for nonsymmetric systems
• Jacobi and Gauss-Seidel (classical methods)

Other

• Dense matrix class with column-major storage (BLAS/LAPACK convention)
• CSR sparse matrix with sparse matrix-vector product
• BLAS-like routines (axpy, gemv, gemm, trsv, trsm)
• Linear system solve, inverse, determinant, condition number
• Least squares via QR

Project structure

├── include/nla/       C++ header-only library
├── python/            Python implementation (numpy-backed, custom algorithms)
├── r/                 R package implementation
├── tests/             C++ tests
├── examples/          Usage examples
└── benchmarks/        Performance comparisons

Building the C++ library

mkdir build && cd build
cmake .. -DNLA_BUILD_TESTS=ON -DNLA_BUILD_EXAMPLES=ON
make -j$(nproc)
ctest  # run tests

The library is header-only so you can also just add include/ to your include path.

Python

cd python
pip install -e .
python tests/test_all.py

from nla_engine import DenseMatrix, lu, qr, solve, cg, eig_symmetric

A = DenseMatrix([[4, -1, 0], [-1, 4, -1], [0, -1, 3]])
b = DenseMatrix([[15], [10], [10]])
x = solve(A, b)

# eigenvalues of a symmetric matrix
vals, vecs = eig_symmetric(A)

# CG on a large SPD system
A_big = DenseMatrix.random_spd(500)
b_big = DenseMatrix.random(500, 1)
x, iters, residual, converged = cg(A_big, b_big)

R

source("r/R/decomposition.R")
source("r/R/solvers.R")
source("r/R/eigen.R")

A <- matrix(c(4, -1, 0, -1, 4, -1, 0, -1, 3), nrow=3, byrow=TRUE)
b <- c(15, 10, 10)
x <- nla_solve(A, b)

eig <- nla_eig_symmetric(A)
print(eig$values)

Algorithm references

Most implementations follow the standard algorithms from:

• Golub & Van Loan, Matrix Computations (4th ed)
• Trefethen & Bau, Numerical Linear Algebra
• Demmel, Applied Numerical Linear Algebra

See docs/algorithms.md for details on each algorithm.

License

MIT