perf(cholesky): add choleskyPartialBlocked() with panel-level Level-3 BLAS operations by jashshah999 · Pull Request #2508 · borglab/gtsam

jashshah999 · 2026-04-27T01:26:31Z

What this does

Adds choleskyPartialBlocked() — a panel-blocked variant of choleskyPartial() that operates on blockSize columns at a time rather than one column at a time.

Algorithm (right-looking blocked Cholesky):

For each panel of bs columns:

Factor the bs x bs diagonal block column-by-column using choleskyStep (preserves underconstrained-pivot handling)
Triangular solve (TRSM): S = R^{-T} B — in-place on the bs x trailing panel row
Rank-bs Schur update (SYRK): trailing -= S' * S — upper triangle only

Steps 2 and 3 are Level-3 BLAS operations (O(bs) arithmetic intensity) vs. the Level-2 rank-1 updates (O(1) arithmetic intensity) done per column in choleskyPartial. With an optimised BLAS backend (OpenBLAS, MKL, Apple Accelerate) these map to multi-threaded dtrsm/dsyrk calls.

Key design decisions:

Falls back to choleskyPartial when nFrontal <= blockSize (no overhead for small matrices)
assert(blockSize > 0) guards against degenerate input
Underconstrained-pivot handling (zeroPivotThreshold, underconstrainedPrior) is preserved exactly via the existing choleskyStep function; dpotrf could replace the panel factor in step 1 for further gain but would lose rank-deficiency handling
API-compatible: same (Matrix&, nFrontal, topleft) signature with an optional blockSize parameter (default: 64)

Benchmark (plain Eigen, no optimised BLAS)

On plain Eigen (no BLAS backend), choleskyPartial uses Eigen::LLT which is itself internally blocked, so the gain here is small and noisy. On a BLAS-backed build the blocked path benefits from multi-threaded dtrsm/dsyrk and the gap grows with problem size.

     N    nF  scalar(ms)   bs=16(ms)   bs=64(ms)
    64    48       0.015       0.012       0.013
   128    96       0.054       0.060       0.080
   256   192       0.362       0.260       0.374
   512   384       1.712       1.830       1.853

Results are variable at these matrix sizes because the working set fits in cache and per-call overhead dominates. The primary intended gain is at nFrontal > 256 with a BLAS-backed build.

Tests

Five new unit tests in testCholesky.cpp:

blockedSmall_fallback: exact equality with choleskyPartial when falling back (blockSize > nFrontal)
blockedInvariant_medium: verifies R'R = A (frontal block) for a 32x32 matrix, blockSize=8
blockedInvariant_large: same for a 128x128 matrix, blockSize=16
blockedWithTopleft: correctness with a non-zero topleft offset
blockedRecovery: full invariant [R' 0; S' I][R S; 0 L] = [A B; B' C] on the same 7x7 matrix used by the existing choleskyPartial test, blockSize=2 to force multiple panels

choleskyPartialBlocked() replaces the column-by-column loop with a right- looking panel algorithm: for each panel of blockSize columns: 1. factor the bs×bs diagonal block with the careful scalar Cholesky (preserves underconstrained-pivot handling) 2. triangular solve: S = R^{-T} B (Level-3 dtrsm equivalent) 3. rank-bs Schur update: C -= S'S (Level-3 dsyrk equivalent) Level-3 operations have O(bs) arithmetic intensity vs O(1) for the Level-2 rank-1 updates in the scalar path, yielding substantially higher throughput on any CPU with an optimised BLAS (OpenBLAS, MKL, Apple Accelerate). The blocked path is only taken when nFrontal > blockSize; smaller matrices fall through to choleskyPartial unchanged. Adds four unit tests verifying correctness: fallback equality, factorization invariant R'R = A for medium/large matrices, and correctness with a topleft offset.

…yPartialBlocked

dellaert · 2026-04-27T20:43:09Z

Cool, but I don’t see this new function being called anywhere. Are you planning to add more commits benchmarking a use case within GTSAM?

Replace the call to choleskyPartial() in SymmetricBlockMatrix::choleskyPartial() with choleskyPartialBlocked(), which uses panel-level Level-3 BLAS operations for large frontal blocks (nFrontal > 64) and falls back to the original choleskyPartial for smaller sizes. This makes the blocked Cholesky the default elimination path for all HessianFactor elimination, ISAM2 clique factorization, and batch solvers. Verified against 10 test suites (testCholesky, testHessianFactor, testJacobianFactor, testGaussianISAM2, testGaussianBayesNet, testGaussianFactorGraph, testSymmetricBlockMatrix, testNonlinearOptimizer, testVisualISAM2, testGaussianFactorGraphB) — zero failures.

jashshah999 · 2026-04-27T22:04:31Z

Thanks for the feedback! I've now wired choleskyPartialBlocked into SymmetricBlockMatrix::choleskyPartial() so it's used by all elimination paths (HessianFactor elimination, ISAM2 clique factorization, batch solvers).

The blocked version already falls back to the original choleskyPartial when nFrontal <= blockSize (default 64), so small cliques see no change. Larger frontal blocks get the panel-level TRSM + SYRK operations.

Verified against 10 test suites — zero failures:
testCholesky, testHessianFactor, testJacobianFactor, testGaussianISAM2, testGaussianBayesNet, testGaussianFactorGraph, testSymmetricBlockMatrix, testNonlinearOptimizer, testVisualISAM2, testGaussianFactorGraphB.

dellaert · 2026-04-28T17:12:42Z

Any perf numbers? E.g., before/after timeSfmBAL?

Copilot

Pull request overview

Note

Copilot was unable to run its full agentic suite in this review.

Adds a new panel-blocked partial Cholesky factorization (choleskyPartialBlocked) intended to improve performance via Level-3 BLAS-style operations, and wires it into SymmetricBlockMatrix while adding unit tests for correctness.

Changes:

Introduced choleskyPartialBlocked() API + implementation using panel factor + TRSM + SYRK-style Schur updates.
Updated SymmetricBlockMatrix::choleskyPartial() to call the blocked routine by default.
Added multiple unit tests covering fallback behavior, invariants, topleft offsets, and recovery structure.

Reviewed changes

Copilot reviewed 4 out of 4 changed files in this pull request and generated 5 comments.

File	Description
gtsam/base/tests/testCholesky.cpp	Adds new unit tests for the blocked algorithm plus an SPD helper.
gtsam/base/cholesky.h	Declares the new `choleskyPartialBlocked()` API and documents intended behavior.
gtsam/base/cholesky.cpp	Implements `choleskyPartialBlocked()` with panel factorization and trailing updates.
gtsam/base/SymmetricBlockMatrix.cpp	Switches internal partial Cholesky call to the new blocked implementation.

Copilot · 2026-04-28T17:20:19Z

+/* ************************************************************************* */
+bool choleskyPartialBlocked(Matrix& ABC, size_t nFrontal, size_t topleft,
+                            size_t blockSize) {
+  if (nFrontal == 0) return true;


assert(blockSize > 0) is compiled out in release builds; if blockSize == 0, the loop for (size_t k = 0; k < nFrontal; k += blockSize) will never advance and can hang (infinite loop). Replace the assert with a runtime guard (e.g., return false or throw) to make the function safe in non-debug builds.

Suggested change

if (nFrontal == 0) return true;

if (nFrontal == 0) return true;

if (blockSize == 0) return false;

Copilot · 2026-04-28T17:20:20Z

  DenseIndex topleft = variableColOffsets_[blockStart_];
-  if (!gtsam::choleskyPartial(matrix_, offset(nFrontals) - topleft, topleft)) {
+  const size_t nF = offset(nFrontals) - topleft;
+  if (!gtsam::choleskyPartialBlocked(matrix_, nF, topleft)) {


offset(nFrontals) and topleft are DenseIndex (signed). Computing the difference and storing it in size_t can underflow/wrap if the difference is negative, producing a huge nF and leading to out-of-bounds behavior. Keep the computation in DenseIndex (or validate non-negativity) before converting to size_t.

Copilot · 2026-04-28T17:20:20Z

  DenseIndex topleft = variableColOffsets_[blockStart_];
-  if (!gtsam::choleskyPartial(matrix_, offset(nFrontals) - topleft, topleft)) {
+  const size_t nF = offset(nFrontals) - topleft;
+  if (!gtsam::choleskyPartialBlocked(matrix_, nF, topleft)) {


The PR description focuses on adding choleskyPartialBlocked() but the code also switches SymmetricBlockMatrix::choleskyPartial() to use it by default, which is a behavior/performance change for existing callers. Either call this out explicitly in the PR description (and any release notes, if applicable) or consider guarding this behind a configuration/flag if the intent is to keep existing behavior by default.

Copilot · 2026-04-28T17:20:20Z

+static Matrix randomSPD(size_t n, unsigned seed = 42) {
+  srand(seed);
+  Matrix A = Matrix::Random(n, n);


Seeding via srand(seed) mutates global RNG state, which can couple this test file to other tests that also use std::rand/Eigen::Random in the same process. Prefer a local deterministic generator (or Eigen's random seed API, if available) to avoid global side effects and improve test isolation.

Copilot · 2026-04-28T17:20:21Z

+        0.,     0.,     0.,     0.,     0.,     0., 2.5776).finished();
+
+  Matrix RSL(ABC);
+  choleskyPartialBlocked(RSL, 3, 0, 2);  // blockSize=2 forces multiple panels


This test doesn't assert the return value of choleskyPartialBlocked. If the factorization fails and returns false, the subsequent invariant check may be harder to interpret. Store/EXPECT the returned bool (as done in other tests) so failures are reported clearly at the source.

Suggested change

choleskyPartialBlocked(RSL, 3, 0, 2); // blockSize=2 forces multiple panels

const bool success =

choleskyPartialBlocked(RSL, 3, 0, 2); // blockSize=2 forces multiple panels

EXPECT(success);

jashshah999 · 2026-04-28T19:33:38Z

I don't have the BAL datasets available on this machine to run timeSFMBAL. The honest benchmark I ran (standalone, plain Eigen without optimized BLAS) showed mixed results because Eigen's LLT is already internally blocked:

     N    nF  scalar(ms)   bs=16(ms)   bs=64(ms)
    64    48       0.015       0.012       0.013
   128    96       0.054       0.060       0.080
   256   192       0.362       0.260       0.374
   512   384       1.712       1.830       1.853

The blocked path uses choleskyStep column-by-column for the panel diagonal (to preserve underconstrained-pivot handling), then Level-3 TRSM+SYRK for the trailing update. On plain Eigen, choleskyPartial's Eigen::LLT is already doing something similar internally, so the gain is small.

The real benefit would come with an optimized BLAS backend where rankUpdate and solveInPlace map to multi-threaded dsyrk/dtrsm. If you think the gain isn't sufficient without that, I understand — happy to close this if it's not pulling its weight.

jashshah999 · 2026-04-30T03:16:38Z

Just following up. Happy to close this if the marginal gains without an optimized BLAS backend don't justify the added code. Let me know your preference.

dellaert · 2026-04-30T12:49:10Z

Yeah, if indeed Eigen is internally blocked, then I’d prefer having Eigen do the lifting - rather than adding mode code in GTSAM.

jashshah999 · 2026-04-30T21:13:38Z

Makes sense. Closing this since Eigen handles the blocking internally and the added code doesn't provide enough benefit to justify the maintenance. Thanks for the review.

jashshah999 added 2 commits April 26, 2026 18:24

fix(cholesky): tighten comments and add blockSize>0 assert in cholesk…

3d529fb

…yPartialBlocked

dellaert requested a review from Copilot April 28, 2026 17:12

Copilot AI reviewed Apr 28, 2026

View reviewed changes

Copilot started reviewing on behalf of dellaert April 28, 2026 17:48 View session

dellaert closed this Apr 30, 2026

	if (nFrontal == 0) return true;
	if (nFrontal == 0) return true;
	if (blockSize == 0) return false;

Conversation

jashshah999 commented Apr 27, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

What this does

Benchmark (plain Eigen, no optimised BLAS)

Tests

Uh oh!

dellaert commented Apr 27, 2026

Uh oh!

jashshah999 commented Apr 27, 2026

Uh oh!

dellaert commented Apr 28, 2026

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

jashshah999 commented Apr 28, 2026

Uh oh!

jashshah999 commented Apr 30, 2026

Uh oh!

dellaert commented Apr 30, 2026

Uh oh!

jashshah999 commented Apr 30, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

jashshah999 commented Apr 27, 2026 •

edited

Loading