Sparse Linear Algebra

AppleAccelerate wraps Apple's Sparse Solvers library for sparse matrix operations and direct solvers.

AASparseMatrix

A wrapper around Apple's SparseMatrix format, constructed from Julia's SparseMatrixCSC.

using AppleAccelerate, SparseArrays
import AppleAccelerate: AASparseMatrix, muladd!

A_jl = sprandn(100, 100, 0.05)
A = AASparseMatrix(A_jl)

x = randn(100)
y = A * x  # Sparse matrix-vector multiply

The constructor automatically detects symmetric and triangular structure and sets the appropriate Apple Accelerate attributes.

Matrix operations

Function	Description
`AASparseMatrix`	Construct from Julia sparse matrix
`A * x`	Sparse matrix-vector or matrix-matrix multiply
`alpha * A * x`	Scaled sparse multiply
`muladd!`	Multiply-add: `y += A * x` or `y += alpha * A * x`
`transpose(A)`	Transpose (sets flag, no copy)

Query functions

Function	Description
`size(A)`	Matrix dimensions
`eltype(A)`	Element type
`issymmetric(A)`	Check if symmetric
`istriu(A)`	Check if upper triangular
`istril(A)`	Check if lower triangular
`A[i, j]`	Element access

AAFactorization

Wraps Apple's SparseOpaqueFactorization. Lazy factorization wrapper: the factorization is computed on the first call to solve or by explicitly calling factor!.

import AppleAccelerate: AAFactorization, solve, solve!, factor!

A = sprandn(100, 100, 0.1) + 20I
f = AAFactorization(A)

# Factorization computed lazily on first solve
b = randn(100)
x = solve(f, b)

# Or explicitly
factor!(f)
x = solve(f, b)

Factorization types

Type	Use case
`SparseFactorizationQR`	Default for non-symmetric matrices
`SparseFactorizationCholesky`	Default for symmetric positive definite
`SparseFactorizationLDLT`	Symmetric indefinite (default LDLT)
`SparseFactorizationLDLTUnpivoted`	Symmetric indefinite, no pivoting
`SparseFactorizationLDLTSBK`	Symmetric indefinite, Bunch-Kaufman
`SparseFactorizationLDLTTPP`	Symmetric indefinite, threshold partial pivoting
`SparseFactorizationCholeskyAtA`	Cholesky of A'A (for least squares)

Solve functions

Function	Description
`AAFactorization`	Lazy factorization wrapper
`solve(f, b)`	Solve `Ax = b`, returns new vector/matrix
`solve!(f, xb)`	Solve in-place (`xb` is overwritten with solution)
`f \ b`	Equivalent to `solve(f, b)`
`ldiv!(f, xb)`	Equivalent to `solve!(f, xb)`
`ldiv!(x, f, b)`	Solve `Ax = b`, store result in `x`
`factor!(f)`	Explicitly compute the factorization
`factor!(f, type)`	Compute factorization with specific type
`factorize(A::AASparseMatrix)`	Create an `AAFactorization` from a sparse matrix

Factorization conveniences

For matrices wrapped as an AASparseMatrix, the usual LinearAlgebra factorization spellings build an AAFactorization of the requested kind:

import AppleAccelerate: AASparseMatrix

A = AASparseMatrix(sprandn(100, 100, 0.1) + 20I)

F = lu(A)        # AAFactorization (LU; requires macOS 15.5+)
G = qr(A)        # QR (also works for rectangular / least-squares)

cholesky(A) and ldlt(A) are available for symmetric/Hermitian matrices.

No type piracy

These methods dispatch on AppleAccelerate's own AASparseMatrix, not on SparseMatrixCSC, so they do not override Julia's built-in lu(::SparseMatrixCSC) / cholesky / qr / ldlt (UMFPACK/CHOLMOD/SPQR). You opt in to the Accelerate solvers by wrapping your matrix in AASparseMatrix first, e.g. lu(AASparseMatrix(A)).

Reusing a factorization: `refactor!`

When a matrix changes its values but not its sparsity pattern — the common case in Newton iterations, implicit time stepping, and parameter sweeps — refactor! recomputes the numeric factorization in place while reusing the existing symbolic factorization (the fill-reducing ordering and sparsity analysis). This is substantially cheaper than building a fresh AAFactorization.

import AppleAccelerate: AAFactorization, refactor!, solve

P = sprandn(100, 100, 0.05) + 20I
F = AAFactorization(P)
solve(F, randn(100))             # forces the first (full) factorization

# Same sparsity pattern, new values:
P2 = copy(P); P2.nzval .*= 1.5
refactor!(F, P2)                 # reuses the symbolic factorization
x = solve(F, randn(100))

F must already hold a completed factorization (call factor!/solve once first), and A must have the same number of stored nonzeros as the original matrix.

For LU/Cholesky/LDLᵀ factorizations you can equivalently use the LinearAlgebra-style spellings that mirror how SparseArrays exposes symbolic reuse — lu!(F, A), cholesky!(F, A), ldlt!(F, A). Each requires F to already hold a factorization of the matching kind and delegates to refactor!. QR reuse has no stdlib spelling, so use refactor! directly for it.

Round-tripping to `SparseMatrixCSC`

An AASparseMatrix can be materialized back to a standard Julia SparseMatrixCSC. The transpose/adjoint/symmetric/Hermitian/triangular attribute bits are honored, so the result equals the logical matrix the wrapper represents:

A = AASparseMatrix(sprandn(50, 50, 0.1))
B = SparseMatrixCSC(A)
@assert B ≈ SparseMatrixCSC(A)   # round-trips the logical matrix

AppleAccelerate.AASparseMatrix — Type

Matrix wrapper, containing the Apple sparse matrix struct and the pointed-to data. Construct from a SparseMatrixCSC.

Multiplication (*) and multiply-add (muladd!) with both Vector and Matrix objects are working. transpose creates a new matrix structure with the opposite transpose flag, that references the same CSC data.

source

AppleAccelerate.AAFactorization — Type

Factorization object.

Create via f = AAFactorization(A::SparseMatrixCSC{T, Int64}). Calls to solve, ldiv, and their in-place versions require explicitly passing in the factorization object as the first argument. On construction, the struct stores a placeholder yet-to-be-factored object: the factorization is computed upon the first call to solve, or by explicitly calling factor!. If the matrix is symmetric, it defaults to a Cholesky factorization; otherwise, it defaults to QR.

source

AppleAccelerate.muladd! — Function

Computes y += A*x in place. Note that this modifies its LAST argument.

source

Computes y += alphaAx in place. Note that this modifies its LAST argument.

source

AppleAccelerate.factor! — Function

factor!(f::AAFactorization, [type::SparseFactorization_t])

Explicitly compute the factorization. If type is not specified, the default is chosen from the matrix's kind attribute:

Hermitian (real symmetric or complex Hermitian) → Cholesky
square, non-Hermitian → LU (macOS 15.5+; falls back to QR on older systems)
rectangular → QR

Called automatically by solve if the factorization has not yet been computed.

Complex symmetric (not Hermitian)

When the user leaves the factorization type unspecified and passes a complex symmetric (but not Hermitian) matrix, we default to LU factorization.

On older versions of MacOS, Apple's SparseFactorizationLDLT errors on complex symmetric (non-Hermitian) matrices, reporting "Cannot perform Hermitian matrix factorization of non-Hermitian matrix." On newer versions, it accepts the call and performs a true LDLᵀ (not LDLᴴ). This behavior isn't documented by Apple; the exact version threshold is unknown.

source

AppleAccelerate.solve — Function

solve(f::AAFactorization, b::StridedVecOrMat)

Solve the linear system Ax = b using Apple's Sparse Solvers, returning the solution x. The factorization is computed lazily on the first call if not already factored. Equivalent to f \ b.

source

AppleAccelerate.solve! — Function

solve!(f::AAFactorization, xb::StridedVecOrMat)

Solve the linear system Ax = b in-place, overwriting xb with the solution. On input xb contains the right-hand side b; on output it contains the solution x. Equivalent to ldiv!(f, xb).

source

AppleAccelerate.refactor! — Function

refactor!(f::AAFactorization, A::AASparseMatrix)
refactor!(f::AAFactorization, A::SparseMatrixCSC)

Recompute the numeric factorization stored in f using the values of A, reusing the existing symbolic factorization (the fill-reducing ordering and sparsity analysis). A must have the same sparsity pattern as the matrix f was originally factored from; only its numeric values may differ.

This is substantially cheaper than building a fresh AAFactorization whenever a matrix changes values but not structure — the common case in Newton iterations, implicit time stepping, and parameter sweeps.

For LU/Cholesky/LDLᵀ factorizations you can equivalently use the LinearAlgebra-style spellings lu!(f, A), cholesky!(f, A), and ldlt!(f, A) (matching SparseArrays' symbolic-reuse API); they require f to already hold a factorization of the matching kind and delegate here. QR reuse has no stdlib spelling, so use refactor! directly for it.

f must already hold a completed factorization (call factor! or solve at least once first); otherwise an ArgumentError is thrown. The factorization object f is mutated in place and returned. After refactor!, f references A's data, so keep A alive as long as f is used.

Warning

Apple's libSparse does not validate that the new pattern matches the old one. Passing a matrix with a different sparsity pattern is undefined behavior. The dimensions and stored nonzero count are checked here as a cheap guard, but an identical count with a different pattern will not be caught.

source