SwifyPy Numerical Mathematics Library
Protocol-oriented numerical computing in Swift with BLAS-accelerated matrix operations
Overview
A numerical computing library in Swift implementing a full linear algebra stack from scratch — six specialized matrix types sharing a unified protocol, LU/QR/Hessenberg decompositions, eigenvalue solvers with Rayleigh quotient shifts and deflation, five numerical integration strategies, ODE solvers, and spectral clustering. Matrix multiplication delegates to Apple's Accelerate framework via a protocol-based BLAS dispatch that selects the correct routine at compile time with zero runtime branching.
Architecture
The core abstraction is MatrixProtocol with an associated MatrixScalar type that dispatches to BLAS. Six concrete matrix types compose by wrapping simpler types (SymTridiagonalMatrix wraps SymmetricBandMatrix wraps UpperBandMatrix), and LU decomposition naturally decomposes into the correct sub-types. The Transposable protocol is deliberately separate because the transpose of a lower-band matrix cannot be a lower-band matrix — a type-theoretic design choice.
Key Concepts
Protocol-Based BLAS Dispatch
The MatrixScalar protocol defines gemm, rotg, and rot as static requirements. Float and Double conform by calling cblas_sgemm/cblas_dgemm respectively. Matrix multiplication dispatches via Scalar.gemm(m1, m2) — the compiler selects the correct BLAS routine through protocol witness tables with no runtime branching.
Code Highlights
public protocol MatrixScalar: BinaryFloatingPoint {
static func gemm(_ m1: Matrix<Self>, _ m2: Matrix<Self>) -> Matrix<Self>
static func rotg(a: Self, b: Self) -> (Self, Self)
static func rot(_ v1: inout Vector<Self>, _ v2: inout Vector<Self>,
_ a: Self, _ b: Self)
}
extension Double: MatrixScalar {
public static func gemm(_ m1: Matrix<Self>, _ m2: Matrix<Self>)
-> Matrix<Self> {
dgemm(transA: false, transB: false, rowsA: m1.height,
colsB: m2.width, colsA: m1.width,
buffA: m1.buffer, buffB: m2.buffer)
}
}Highlights
- Protocol-oriented numerical computing: six matrix types sharing a MatrixProtocol interface with compile-time BLAS dispatch through protocol witness tables
- Sparsity-exploiting band matrix LU decomposition reducing O(n^3) to O(k*n^2), with tests verifying exact iteration counts
- Full eigenvalue solver: QR algorithm with Rayleigh quotient shifts, progressive deflation via LazyMatrixView, and eigenvector accumulation
- Self-referential architecture: Gauss-Legendre integration computes quadrature nodes as eigenvalues of the Jacobi tridiagonal matrix using the library's own solver
Related Projects
Music Classification with Deep Learning
IEEE-format research paper with novel noise-augmentation for robust music identification
Published IEEE-format research paper presenting a novel noise-augmentation technique for robust music identification — training for real-world degradation rather than clean-sample accuracy
Expanding-Shrinking Network Model
A novel random graph model combining growth and shrinkage dynamics to produce networks with real-world properties
Novel ES model: at each step, either merge two random nodes (shrink) or split a random node inheriting all edges (expand) — controlled by a single α parameter