macldlt

Python wrapper for Apple Accelerate's sparse LDL^T factorization.

macldlt exposes the symmetric indefinite sparse solver from Apple's Accelerate framework to Python via pybind11. It accepts SciPy sparse matrices and NumPy arrays directly, with no manual conversion needed.

macOS only — requires macOS 13.0+.

Published wheels target macosx_13_0_arm64, which means macOS 13.0 is the minimum supported version; they are intended to run on macOS 13 and newer on Apple Silicon.

Prebuilt wheels are published for Apple Silicon (macOS arm64) only.

Installation

pip install macldlt

Or install from source:

git clone https://github.com/bodono/macldlt.git
cd macldlt
pip install -e ".[test]"

Building from source requires:

• macOS 13.0+
• A C++17 compiler (Xcode command-line tools)
• Python >= 3.10
• pybind11 >= 2.12

Quick start

import numpy as np
import scipy.sparse as sp
from macldlt import LDLTSolver

# Build a symmetric positive-definite matrix
A = sp.csc_matrix(np.array([
    [ 4.0, 1.0],
    [ 1.0, 3.0],
]))

solver = LDLTSolver(A)
b = np.array([1.0, 2.0])
x = solver.solve(b)
print(x)  # [0.09090909 0.63636364]

API reference

LDLTSolver(A, triangle="upper", ordering="amd", factorization="ldlt")

Perform symbolic analysis and numeric factorization of A. The solver is immediately ready to call solve().

Parameters:

Notes:

• Symmetry is assumed but not checked. Only the specified triangle is read; the other triangle is ignored. If you pass a full symmetric matrix, set triangle to whichever triangle contains the data you want used.
• The solver is not thread-safe. Do not call methods concurrently on the same instance from multiple threads.
• For a new sparsity pattern, create a new solver.

Example:

# Upper triangle of a 3x3 symmetric matrix
A_upper = sp.csc_matrix(np.array([
    [4.0, 1.0, 0.0],
    [0.0, 3.0, 2.0],
    [0.0, 0.0, 5.0],
]))
solver = LDLTSolver(A_upper, triangle="upper")

solver.refactor(values)

Reuse the existing symbolic analysis and numeric factorization workspace for new values with the same sparsity pattern. This calls Accelerate's SparseRefactor, which can be faster than a full factor().

Parameters:

# In a tight loop, just pass new values
for new_values in value_generator:
    solver.refactor(new_values)
    x = solver.solve(b)

solver.solve(rhs, inplace=False)

Solve Ax = rhs.

By default, allocates and returns a new NumPy array. With inplace=True, overwrites rhs in place and returns it (avoids allocation).

Parameters:

Returns: numpy.ndarray — the solution, same shape as rhs.

# Allocating solve
x = solver.solve(b)

# In-place solve (no allocation)
solver.solve(b, inplace=True)
# b now contains the solution

solver.inertia()

Return the inertia of the factored matrix as a tuple (num_negative, num_zero, num_positive), where:

• num_negative — number of negative pivots
• num_zero — number of zero pivots
• num_positive — number of positive pivots

The sum num_negative + num_zero + num_positive equals n.

neg, zero, pos = solver.inertia()
if zero > 0:
    print("Matrix is singular")
if neg == 0 and zero == 0:
    print("Matrix is positive definite")

solver.info()

Return a dictionary with solver state and workspace information.

Returns: dict with keys:

Properties

Typical workflow

solver = LDLTSolver(A) ──► solve()              # one-shot usage
              │
              └── refactor(values) ──► solve()   # same pattern, new values

solver = LDLTSolver(A_new) ──► solve()           # new sparsity pattern

1. One-shot solve: Pass A to the constructor, then call solve().
2. Repeated solves, same pattern: Call refactor(new_vals) then solve(). The symbolic analysis from the constructor is reused.
3. New sparsity pattern: Create a new LDLTSolver with the new matrix.

Triangle conventions

Accelerate's symmetric solver reads only one triangle of the matrix. You specify which triangle is stored via the triangle parameter:

import numpy as np
import scipy.sparse as sp

A_full = np.array([
    [4.0, 1.0],
    [1.0, 3.0],
])

# If you store the upper triangle:
A_upper = sp.csc_matrix(np.triu(A_full))
solver = LDLTSolver(A_upper, triangle="upper")

# If you store the lower triangle:
A_lower = sp.csc_matrix(np.tril(A_full))
solver = LDLTSolver(A_lower, triangle="lower")

If you pass a full symmetric matrix with triangle="upper", only the upper triangle entries are used — the lower triangle is ignored (and vice versa).

License

MIT