Fix GPU Cholesky cache initialization for non-square matrices by ChrisRackauckas-Claude · Pull Request #857 · SciML/LinearSolve.jl

ChrisRackauckas-Claude · 2025-12-16T18:31:28Z

Summary

When using DefaultLinearSolver with non-square GPU matrices (e.g., for least squares problems in NonlinearSolve.jl), the init_cacheval function for CholeskyFactorization would fail with DimensionMismatch: matrix is not square because it tried to compute cholesky(A) on a non-square matrix.

Root Cause

The @generated init_cacheval for DefaultLinearSolver initializes caches for all possible algorithms upfront, including CholeskyFactorization. For non-square GPU matrices, this fails even though CholeskyFactorization would never actually be used (the defaultalg function correctly selects QRFactorization for non-square matrices).

Fix

The fix checks assumptions.issq before attempting Cholesky factorization in init_cacheval(::CholeskyFactorization, ::GPUArraysCore.AnyGPUArray, ...) and returns nothing for non-square matrices, allowing the cache initialization to succeed.

Changes

src/factorization.jl: Added square matrix check in init_cacheval for CholeskyFactorization with GPU arrays
test/gpu/cuda.jl: Added tests for non-square GPU matrices (overdetermined and underdetermined systems)

Test plan

LinearSolve package loads correctly
CI tests pass
New CUDA tests for non-square matrices pass

🤖 Generated with Claude Code

When using DefaultLinearSolver with non-square GPU matrices (e.g., for least squares problems), the init_cacheval function for CholeskyFactorization would fail because it tried to compute cholesky(A) on a non-square matrix. The fix checks assumptions.issq before attempting Cholesky factorization and returns nothing for non-square matrices, allowing the DefaultLinearSolver to properly use QRFactorization instead. Fixes SciML/NonlinearSolve.jl#746 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

Add check for assumptions.issq in the SparseArrays extension's NormalCholeskyFactorization init_cacheval to avoid DimensionMismatch errors when DefaultLinearSolver initializes caches for non-square GPU matrices. This is part of the fix for SciML/NonlinearSolve.jl#746 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

ChrisRackauckas-Claude · 2025-12-16T21:55:15Z

Added a second fix for the NormalCholeskyFactorization init_cacheval function in ext/LinearSolveSparseArraysExt.jl.

The previous commit fixed CholeskyFactorization in src/factorization.jl, but the GPU tests still failed because the NormalCholeskyFactorization dispatch in the SparseArrays extension was calling ArrayInterface.cholesky_instance(convert(AbstractMatrix, A)) directly on the matrix A without checking if it's square.

This commit adds the same assumptions.issq check to avoid the DimensionMismatch error during DefaultLinearSolver cache initialization for non-square GPU matrices.

🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

ChrisRackauckas-Claude · 2025-12-16T22:04:15Z

CI Status Update

The key downstream tests are passing:

✅ NonlinearSolve.jl/Core/lts
✅ SciMLSensitivity.jl/Core1/lts
✅ Spell Check with Typos

The failing tests appear to be pre-existing issues:

Format Check - The SciML format checker checks the entire repo, not just changed files. The changed files have been formatted with SciMLStyle.
Core alldeps (1.10) - Package resolution failure between ChainRulesCore (pinned to 1.25.0) and Zygote (requires >=1.25.1). This is a dependency issue, not related to this PR.

GPU tests run on Buildkite separately and will validate the fix for non-square GPU matrices.

The previous change incorrectly returned nothing for ALL non-square matrices (sparse and GPU), but NormalCholeskyFactorization actually works with non-square matrices because it factors A'*A which is square. The issue was only with GPU arrays where ArrayInterface.cholesky_instance fails for non-square. Sparse CPU arrays handle this correctly. Now the check only applies to GPUArraysCore.AnyGPUArray, allowing sparse CPU matrices to use the original cholesky_instance path. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

ChrisRackauckas-Claude · 2025-12-16T22:50:11Z

Fix for Non-Square Tests Failure

The previous change incorrectly returned nothing for ALL non-square matrices in NormalCholeskyFactorization init_cacheval. However, NormalCholeskyFactorization actually works with non-square matrices because it factors A' * A (which is always square).

The issue was only with GPU arrays where ArrayInterface.cholesky_instance(A) fails for non-square A. Sparse CPU arrays handle cholesky_instance correctly for non-square matrices.

Fix: Changed the condition to only apply to GPUArraysCore.AnyGPUArray, allowing sparse CPU matrices to use the original cholesky_instance path.

# Before (incorrect - broke sparse CPU):
elseif !assumptions.issq
    nothing

# After (correct - only affects GPU):
elseif A isa LinearSolve.GPUArraysCore.AnyGPUArray && !assumptions.issq
    nothing

ChrisRackauckas-Claude · 2025-12-16T23:09:04Z

✅ All Core Tests Passing

After the fix to only apply the non-square check to GPU arrays (not sparse CPU arrays), all Core tests are now passing:

✅ 14/14 Tests (*, *, Core, *) - ALL PASSING
✅ NonlinearSolve.jl/Core/lts
✅ SciMLSensitivity.jl/Core1/lts

The only remaining failures are pre-existing infrastructure issues:

Format Check (repo-wide formatting)
test (Core, alldeps, 1.10) (ChainRulesCore/Zygote dependency resolution)

Summary of changes:

src/factorization.jl: Added assumptions.issq check for CholeskyFactorization init_cacheval with GPU arrays
ext/LinearSolveSparseArraysExt.jl: Added AnyGPUArray && !assumptions.issq check for NormalCholeskyFactorization init_cacheval (GPU arrays only, not sparse CPU)
test/gpu/cuda.jl: Added tests for non-square GPU matrices

The fix allows DefaultLinearSolver to properly initialize caches for non-square GPU matrices by returning nothing for Cholesky-based cachevals when the matrix is not square, allowing the solver to fall back to QRFactorization.

ChrisRackauckas · 2026-08-23T19:09:01Z

@AJ0070 look into rebasing this?

`solve(LinearProblem(A, b))` with a non-square GPU `A` threw `DimensionMismatch: matrix is not square`, in two different places. Cache initialization. `_init_default_cacheval` builds a cacheval for every algorithm slot before it knows which one it will use, and two of those slots call `cholesky` on `A` itself: - `CholeskyFactorization` on `AnyGPUArray`, which calls `cholesky(A; check = false)`. - `NormalCholeskyFactorization` on `AnyGPUArray`, which calls `ArrayInterface.cholesky_instance(A)`. That algorithm factorizes `A'A`, so a non-square `A` is a supported least-squares input, and the generic method already sidesteps the shape by instancing an empty `Symmetric`. Only the GPU method built the instance from `A`. The guard is scoped to GPU input: `cholesky_instance` handles a non-square sparse `A`, and `solve!` stores a real factorization into that slot afterwards. The wide solve. Past cache initialization, a wide `A` then failed inside `QRFactorization`: a GPU `qr` factors it, but solving with the result builds `UpperTriangular(R)` on an `R` that is not square. `MinNormQR` already exists for exactly this, used by the BandedMatrices and FastAlmostBandedMatrices extensions, so route wide GPU input through `Aᵀ` the same way. That also makes the answer the minimum-norm one, matching dense `\` on the CPU. Measured on a T4, against main: before tall 4x2 FAIL DimensionMismatch (4, 2) wide 2x4 FAIL DimensionMismatch (2, 4) after tall 4x2 OK relerr=2.47e-6 wide 2x4 OK relerr=5.62e-7 ‖x‖ 0.209165 vs 0.209165 square 3x3 OK relerr=6.91e-8 The wide test asserts `‖x‖` against LAPACK's, not just the residual: agreeing on the residual alone would not distinguish the minimum-norm solution from any other point on the solution manifold. Revives SciML#857, whose guard covered only the `CholeskyFactorization` slot. The `NormalCholeskyFactorization` slot moved into `sparsearrays.jl` in the eight months since, and fails first, so that fix alone no longer resolves the reported failure; verified on a T4 that it does not. Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>

* Solve a non-square GPU system instead of throwing `solve(LinearProblem(A, b))` with a non-square GPU `A` threw `DimensionMismatch: matrix is not square`, in two different places. Cache initialization. `_init_default_cacheval` builds a cacheval for every algorithm slot before it knows which one it will use, and two of those slots call `cholesky` on `A` itself: - `CholeskyFactorization` on `AnyGPUArray`, which calls `cholesky(A; check = false)`. - `NormalCholeskyFactorization` on `AnyGPUArray`, which calls `ArrayInterface.cholesky_instance(A)`. That algorithm factorizes `A'A`, so a non-square `A` is a supported least-squares input, and the generic method already sidesteps the shape by instancing an empty `Symmetric`. Only the GPU method built the instance from `A`. The guard is scoped to GPU input: `cholesky_instance` handles a non-square sparse `A`, and `solve!` stores a real factorization into that slot afterwards. The wide solve. Past cache initialization, a wide `A` then failed inside `QRFactorization`: a GPU `qr` factors it, but solving with the result builds `UpperTriangular(R)` on an `R` that is not square. `MinNormQR` already exists for exactly this, used by the BandedMatrices and FastAlmostBandedMatrices extensions, so route wide GPU input through `Aᵀ` the same way. That also makes the answer the minimum-norm one, matching dense `\` on the CPU. Measured on a T4, against main: before tall 4x2 FAIL DimensionMismatch (4, 2) wide 2x4 FAIL DimensionMismatch (2, 4) after tall 4x2 OK relerr=2.47e-6 wide 2x4 OK relerr=5.62e-7 ‖x‖ 0.209165 vs 0.209165 square 3x3 OK relerr=6.91e-8 The wide test asserts `‖x‖` against LAPACK's, not just the residual: agreeing on the residual alone would not distinguish the minimum-norm solution from any other point on the solution manifold. Revives #857, whose guard covered only the `CholeskyFactorization` slot. The `NormalCholeskyFactorization` slot moved into `sparsearrays.jl` in the eight months since, and fails first, so that fix alone no longer resolves the reported failure; verified on a T4 that it does not. Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com> * Instance an empty Cholesky for a non-square GPU A instead of returning nothing * Instance an empty Cholesky in the NormalCholesky GPU slot too * Compare the QR solve on the host, not against a device array --------- Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>

ChrisRackauckas · 2026-09-06T12:19:18Z

Closing as superseded: 679 commits behind master; the change is either on master already or no longer applies.

🤖 Posted by an AI agent — harness: Claude Code, model: claude-opus-5[1m]. Chris asked for these closures directly.

https://claude.ai/code/session_014FEzNTLFutCmTEAZ3zBg5R

ChrisRackauckas force-pushed the fix-gpu-nonsquare-cholesky-init branch from 91c9aec to 854d90f Compare December 16, 2025 21:17

Format LinearSolveSparseArraysExt.jl with SciMLStyle

9e1b11b

🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

AJ0070 mentioned this pull request Aug 24, 2026

Solve a non-square GPU system instead of throwing #1250

Merged

5 tasks

ChrisRackauckas closed this Sep 6, 2026

Navigation Menu

Uh oh!

Uh oh!

Fix GPU Cholesky cache initialization for non-square matrices - #857

Fix GPU Cholesky cache initialization for non-square matrices#857
ChrisRackauckas-Claude wants to merge 4 commits into
SciML:mainfrom
ChrisRackauckas-Claude:fix-gpu-nonsquare-cholesky-init

ChrisRackauckas-Claude commented Dec 16, 2025

Uh oh!

ChrisRackauckas-Claude commented Dec 16, 2025

Uh oh!

ChrisRackauckas-Claude commented Dec 16, 2025

Uh oh!

ChrisRackauckas-Claude commented Dec 16, 2025

Uh oh!

ChrisRackauckas-Claude commented Dec 16, 2025

Uh oh!

ChrisRackauckas commented Aug 23, 2026

Uh oh!

ChrisRackauckas commented Sep 6, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Uh oh!

Conversation

ChrisRackauckas-Claude commented Dec 16, 2025

Summary

Root Cause

Fix

Changes

Test plan

Uh oh!

ChrisRackauckas-Claude commented Dec 16, 2025

Uh oh!

ChrisRackauckas-Claude commented Dec 16, 2025

CI Status Update

Uh oh!

ChrisRackauckas-Claude commented Dec 16, 2025

Fix for Non-Square Tests Failure

Uh oh!

ChrisRackauckas-Claude commented Dec 16, 2025

✅ All Core Tests Passing

Uh oh!

ChrisRackauckas commented Aug 23, 2026

Uh oh!

ChrisRackauckas commented Sep 6, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants