Add `cho_factor` and `cho_solve` to `dpnp.scipy.linalg` by antonwolfy · Pull Request #3116 · IntelPython/dpnp

antonwolfy · 2026-10-20T22:42:29Z

Add SciPy-compatible dpnp.scipy.linalg.cho_factor and dpnp.scipy.linalg.cho_solve, including batched input as in SciPy 1.18 (#3096).

Implementation

cho_factor uses the existing potrf and potrf_batch bindings. Only the referenced triangle is copied, in one pass, into a zeroed C-order buffer, which matches the SciPy 1.18 output. With overwrite_a=True, a 2-D, writable input of the result dtype that is C- or F-contiguous is factored in place. As in SciPy, a 0-D or 1-D input is viewed as a matrix.
cho_solve uses a new LAPACK potrs binding for a 2-D factor. For a batch of factors it calls a new BLAS trsm_batch binding twice, since oneMKL potrs_batch is much slower (MKLD-20341).
A single factor (2-D, or with size-1 batch dimensions) is not broadcast to a batched right-hand side: the batch is folded into the columns of one potrs call. For a real type on GPU with more than 16 folded columns and a batch of at most 32, the factor is shared by the batch in trsm_batch (stride 0) instead, because double-precision trsm on Xe GPUs slows down by up to 3.6x above 16 columns.
Batch dimensions broadcast as in lu_solve, and a broadcast right-hand side is copied without being materialized first.
lower can also be the batch-shaped array SciPy 1.18 returns, or a sequence, if all its values are equal.
With check_finite=True, both finiteness checks of cho_solve are submitted before the first host sync.

Differences from SciPy

Documented in the docstrings and covered by tests:

cho_factor returns lower as a Python bool; SciPy 1.18 returns an array tiled to the batch shape.
cho_solve supports a zero-size batch; SciPy raises.

As in the existing lu_factor and dpnp.linalg.cholesky, integer and boolean inputs are cast to the default floating type of the device; SciPy 1.18 uses float32 for boolean and up to 16-bit integer inputs.

Other changes

dpnp.linalg.cholesky: a See Also link to dpnp.scipy.linalg.cho_factor and a fixed example output for upper=True.
doc/reference/linalg.rst: removed the stale linalg.lu_factor and linalg.lu_solve entries.
Backported tests/cupyx_tests/linalg_tests/test_solve.py from CuPy, with TestChosolve adapted to dpnp.scipy.linalg.cho_solve and the invh tests skipped.

Testing

New tests in test_linalg.py, test_sycl_queue.py and test_usm_type.py, plus the CuPy tests.
On Panther Lake H (Level Zero GPU, OpenCL GPU and OpenCL CPU devices), the linalg, SYCL queue and USM type tests pass, and a comparison with SciPy 1.18.1 (6720 checks over dtypes including mixed ones, C and F layouts, lower, batching and broadcasting, overwrite_* and check_finite) finds no differences. Earlier revisions were also tested on PVC, Lunar Lake and Arc B580.
The build adds no compiler warnings.

Performance

cho_factor + cho_solve against dpnp.linalg.solve, median times, measured on an earlier revision of this branch:

Device	n=4096	batch (1000, 32, 32)	re-solve: `cho_solve` vs `lu_solve`
PVC Max 1100 GPU	3.1x faster (f64)	84x faster (f64)	1.0–2.0x faster
Panther Lake H Xe3 GPU	1.8–2.3x faster	10–21x faster	1.2–4.5x faster
Arc B580	1.7–2.2x faster	20–43x faster	1.07–3.2x faster

On client and Battlemage GPUs, the Cholesky path is slower than LU for small single matrices (n=256: 0.5–0.6x) and for batches of larger matrices such as (100, 128, 128) (0.24–0.5x), because of oneMKL issues in potrf and potrf_batch.

Have you provided a meaningful PR description?
Have you added a test, reproducer or referred to an issue with a reproducer?
Have you tested your changes locally for CPU and GPU devices?
Have you made sure that new changes do not introduce compiler warnings?
Have you checked performance impact of proposed changes?
Have you added documentation for your changes, if necessary?
Have you added your changes to the changelog?

github-actions · 2026-10-20T23:20:26Z

View rendered docs @ https://intelpython.github.io/dpnp/pull/3116/index.html

github-actions · 2026-10-20T23:30:06Z

Array API standard conformance tests for dpnp=0.22.0dev1=np2py314h8d9cdd5_21 ran successfully.
Passed: 1491
Failed: 0
Skipped: 6

coveralls · 2026-10-21T00:44:31Z

coverage: 84.21% (-0.2%) from 84.397% — cho-factor-cho-solve into master

Copilot

🟡 Changes recommended

Half-precision inputs can reach unsupported Cholesky dispatches, and overwriting a mixed-USM right-hand side can return a solution with the wrong USM type.

2 open issues

What changed in this PR

Adds SciPy-compatible Cholesky factorization and solving to dpnp.scipy.linalg, enabling positive-definite systems to use Cholesky rather than general LU solves.

Changes:

Adds cho_factor and cho_solve with batching, broadcasting, overwrite options, and finite-value checks.
Introduces oneMKL potrs and batched triangular-solve bindings.
Adds numerical, queue, and USM tests, plus API documentation and a changelog entry.

File	Description
`dpnp/tests/third_party/cupyx/linalg_tests/test_solve.py`	Adapts CuPy Cholesky-solve tests; skips unsupported `invh`.
`dpnp/tests/third_party/cupyx/linalg_tests/__init__.py`	Adds the test package marker.
`dpnp/tests/test_usm_type.py`	Checks result USM allocation types.
`dpnp/tests/test_sycl_queue.py`	Checks queue preservation and mixed-queue rejection.
`dpnp/tests/test_linalg.py`	Tests numerical results, batching, layouts, and validation.
`dpnp/scipy/linalg/_utils.py`	Implements factorization and optimized solve paths.
`dpnp/scipy/linalg/_decomp_cholesky.py`	Defines and documents the public Cholesky APIs.
`dpnp/scipy/linalg/__init__.py`	Exports both new functions.
`dpnp/linalg/dpnp_iface_linalg.py`	Adds a documentation cross-reference and corrects an example.
`dpnp/backend/extensions/lapack/types_matrix.hpp`	Defines supported `potrs` types.
`dpnp/backend/extensions/lapack/potrs.hpp`	Declares the `potrs` binding.
`dpnp/backend/extensions/lapack/potrs.cpp`	Implements validated oneMKL Cholesky solves.
`dpnp/backend/extensions/lapack/lapack_py.cpp`	Exposes `potrs` to Python.
`dpnp/backend/extensions/lapack/CMakeLists.txt`	Includes `potrs` in the build.
`dpnp/backend/extensions/blas/types_matrix.hpp`	Defines supported batched triangular-solve types.
`dpnp/backend/extensions/blas/trsm.hpp`	Declares the batched triangular-solve binding.
`dpnp/backend/extensions/blas/trsm_batch.cpp`	Implements batched solves with shared-factor support.
`dpnp/backend/extensions/blas/CMakeLists.txt`	Includes batched triangular solves in the build.
`dpnp/backend/extensions/blas/blas_py.cpp`	Exposes `trsm_batch` to Python.
`doc/reference/scipy_linalg.rst`	Lists the new Cholesky functions.
`doc/reference/linalg.rst`	Removes stale LU API entries.
`CHANGELOG.md`	Records the new public APIs.

_{Review effort: Balanced}

antonwolfy · 2026-10-30T18:34:34Z

`Ridge` benchmark with the new functions

To see what cho_factor / cho_solve give in a real workload, I benchmarked scikit-learn's Ridge.fit with dpnp inputs on a Panther Lake H Xe3 iGPU, using the scenario from scikit-learn issue scikit-learn/scikit-learn#35060 (array API support for solver="cholesky").

Setup:

scikit-learn 1.9.1, whose array API Ridge only supports solver="svd"; SciPy 1.18.1, NumPy 2.4.3, dpnp from this branch, oneAPI 2026.2.
Ridge(alpha=1.0, fit_intercept=True).fit with one target, standard normal data already on the device, fit time only.
As in the prototype from that issue, _solve_svd is swapped during Ridge(solver="svd").fit, so validation, centering and the intercept stay scikit-learn's own code:
- A: Gram (or kernel) matrix on the device, solved on the host with scipy.linalg.solve(assume_a="pos")
- B: dpnp.linalg.solve (LU) on the device
- C: dpnp.scipy.linalg.cho_factor + cho_solve on the device
A, B and C are medians of 6 interleaved rounds. NumPy (stock solver="cholesky" on the 16-core CPU) is the best of 3, and dpnp SVD a single run. The coefficients match NumPy to ~1e-6 (f32) and ~1e-14 (f64).

Fit time, seconds:

Shape	dtype	NumPy CPU	dpnp SVD (today)	A: host solve	B: device LU	C: device Cholesky	C vs B	C vs NumPy
500,000×100	f32	0.072	0.784	0.080	0.093	0.096	0.96x	0.75x
50,000×1000	f32	0.199	1.079	0.105	0.080	0.076	1.06x	2.6x
50,000×1000	f64	0.251	1.687	0.249	0.256	0.237	1.08x	1.06x
463,810×990	f32	1.194	-	0.771	0.758	0.765	0.99x	1.56x
2,500×20,000	f32	0.278	5.696	0.176	0.104	0.095	1.09x	2.9x
2,500×20,000	f64	0.723	6.595	0.498	0.464	0.444	1.05x	1.63x
5,000×40,000	f32	1.969	-	0.766	0.554	0.538	1.03x	3.7x
10,000×40,000	f32	7.210	-	2.386	1.834	1.473	1.25x	4.9x
100,000×5,000	f32	4.360	-	1.451	1.358	1.280	1.06x	3.4x

The last two shapes are not in that issue; they are added to show larger systems.

Solve step alone (the n×n system inside fit), ms:

n (shape)	dtype	B: LU	C: Cholesky	C vs B
100 (500,000×100)	f32	0.64	3.2	0.2x
1000 (50,000×1000)	f32	21.0	11.6	1.8x
1000 (50,000×1000)	f64	30.2	14.0	2.2x
2500 (2,500×20,000)	f32	41.0	19.6	2.1x
2500 (2,500×20,000)	f64	175.8	82.2	2.1x
5000 (5,000×40,000)	f32	165.2	95.6	1.7x
5000 (100,000×5,000)	f32	189.2	102.1	1.9x
10000 (10,000×40,000)	f32	588.4	282.4	2.1x

Takeaways:

Replacing SVD is the large gain (7–60x), and all three variants get it.
The Cholesky solve is about 2x faster than LU from n=1000 up. In fit the change is between −4% and +9% on the shapes from the issue, and +25% on 10,000×40,000, since the Gram matmul dominates; the gain grows with the size of the system.
Cholesky on the device is 1.4–1.9x faster than the host-solve variant A on the f32 dual shapes and on 50,000×1000 f32, and 1.05–1.1x in f64.
For n=100, Cholesky is slower than LU because of the small-n potrf latency on client GPUs (MKLD-20346).
For tall data, most of the fit time is scikit-learn's own validation and centering on dpnp, about 0.3 ms per MiB of X (517 ms of the 765 ms fit at 463,810×990). That is worth a separate look.
The remaining part of Add Cholesky-based linear solvers to dpnp.scipy.linalg: cho_factor, cho_solve, solve_triangular, solve #3096, solve(assume_a="pos"), is what a SciPy-style caller would use to reach this Cholesky path.

Copilot

🔵 Needs a closer look

The new asynchronous potrs/trsm_batch bindings need device-backed oneMath/CUDA/HIP validation beyond the supplied Intel results, which this source-only environment cannot provide.

No open issues 🎉

2 resolved since last review

_{Review effort: Balanced}

Implement SciPy-compatible `dpnp.scipy.linalg.cho_factor` and `dpnp.scipy.linalg.cho_solve` with batched input support (gh-3096). - `cho_factor` uses the existing `potrf` / `potrf_batch` bindings. Only the referenced triangle is copied into a zeroed C-order buffer, which matches the SciPy 1.18 output in a single pass. The input is factorized in place for `overwrite_a=True` when it is 2-D, writable, of the result dtype, and C- or F-contiguous. - `cho_solve` uses a new LAPACK `potrs` binding for 2-D input and a new BLAS `trsm_batch` binding (two triangular solves) for batched input, since oneMKL `potrs_batch` is much slower. Batch dimensions are broadcast as in `lu_solve`. - Add tests and API reference entries.

A C-order factor is passed to oneMKL transposed, so the solves are with conj(A), and `b` and the solution are conjugated to compensate. This is only needed for a complex factor. Decide it from the dtype of `c` rather than the result dtype, so a real factor with a complex right-hand side no longer conjugates `b` and the solution. Add a test for a real factor with complex right-hand sides.

Add `tests/cupyx_tests/linalg_tests/test_solve.py` from CuPy as `dpnp/tests/third_party/cupyx/linalg_tests/test_solve.py`, with `TestChosolve` adapted to `dpnp.scipy.linalg.cho_solve`: - pass `check_finite=False`, since dpnp checks the whole factor for infs and NaNs, as SciPy does, and the tests put NaNs in the unused triangle; - pass a batch of vectors with a trailing axis, since dpnp follows SciPy's `(..., n, k)` rule for a batched right-hand side. The `invh` tests are skipped, as `invh` is not supported.

- Add See Also link from dpnp.linalg.cholesky to dpnp.scipy.linalg.cho_factor - Fix the upper-triangular example output in dpnp.linalg.cholesky - Remove stale linalg.lu_factor and linalg.lu_solve entries from linalg.rst - Look up the order manager in _copy_rhs and _cho_solve_trsm instead of passing it

- Do not broadcast a 2-D factor `c` over a batched right-hand side: fold the batch into the RHS columns for one `potrs` call, or, for a real type on GPU with a small batch (batch * nrhs > 16, batch <= 32), where oneMKL potrs slows down above 16 columns, share the factor by all batch elements in `trsm_batch` - Let `_trsm_batch` take a 2-D triangular matrix shared by the batch (stride 0) - Submit both `check_finite` reductions before waiting for the first result - Accept the batch-shaped `lower` array returned by SciPy >= 1.18 in `cho_solve`, if all its values are equal

- `bool()` of a non-empty list or tuple is always True, so `cho_solve((c, [False]), b)` used the wrong triangle: convert a sequence like an array, as SciPy's `cho_solve` does, and do the same in `cho_factor` - Use finite values in `TestChoSolve.test_empty`: `check_finite` reads `c` or `b` when it is not empty, so uninitialized data could fail it - Add the changelog entry

…ciPy differences - `cho_factor` with `overwrite_a=True` and `cho_solve` with `overwrite_b=True` returned a `usm_ndarray` input as is: wrap it with `dpnp.get_result_array` - Reject a 0-d right-hand side in the `_potrs` binding instead of reading its shape - Document that `cho_factor` returns `lower` as a bool and that `cho_solve` supports a zero-size batch, unlike SciPy

- Move the dispatch rule for a 2-D factor with a batched right-hand side to `_fold_rhs_batch`, so tests can force each path on any device - Sample integers from [-2, 2) in `_make_hpd_np`: [-1, 1) truncated them to zeros, so integer test matrices were diagonal - Test a batched column-major factor, an empty `lower` array, the dispatch boundaries, a non-finite value in the triangle `cho_factor` does not reference, an in-place 1-D `b`, a 1-D `c` and mixed queues, and match the error messages - Drop the direct test of the `_potrs` binding

- List the `ValueError` cases of `cho_factor` next to `LinAlgError` - Add a Raises section to `cho_solve` - Mention the sync to read a dpnp array `lower` in `cho_solve`

- Catch `oneapi::mkl::exception` after the LAPACK one, so an error not from LAPACK (e.g. `unimplemented` or a failed allocation) also frees the scratchpad - Include the standard headers used and drop the unused `linalg_exceptions.hpp`

A `c` of shape (1, n, n) or (1, 1, n, n) with a batched right-hand side was broadcast to the batch and copied per batch element; solve with it as with a 2-D factor instead

- View a 0-D or 1-D `a` as a matrix in `cho_factor`, as SciPy does, and report its shape as `a1.shape` - Raise `ValueError` "The factored matrix c is not square." in `cho_solve` for `c` with fewer than two dimensions - Report the shapes of a single system for incompatible `c` and `b`, also for batched input

- Copy `b` into the solve buffer through a view in its batch shape, so a broadcast `b` is not first reshaped into a full copy - Restore `_align_lu_solve_broadcast` as on master: `cho_solve` checks the shapes before calling it - Include the standard headers used in `potrs.hpp`, `trsm.hpp` and `trsm_batch.cpp`

With `overwrite_b=True`, a right-hand side whose USM type differs from the coerced one of `c` and `b` was solved in place, so the solution kept the USM type of `b`. Copy it in that case, as without `overwrite_b`.

…patch - query the potrs scratchpad size inside the try and raise if its allocation fails - drop the unused unit_diag option of _trsm_batch - use the shared-factor trsm_batch path only on Intel GPUs - note the measured potrs_batch slowdown next to the internal ticket

ndgrigorian · 2026-11-05T15:50:04Z

+    sycl::event clean_up_event = exec_q.submit([&](sycl::handler &cgh) {
+        cgh.depends_on(potrs_event);
+        auto ctx = exec_q.get_context();
+        cgh.host_task([ctx, scratchpad]() {
+            dpnp::tensor::alloc_utils::sycl_free_noexcept(scratchpad, ctx);
+        });
+    });
+    host_task_events.push_back(clean_up_event);


With the allocation now skipped for a zero scratchpad size, this clean-up host task still runs in that case just to free a nullptr, so it's an extra host task per call for nothing. It could be submitted only when scratchpad != nullptr, with keep_args_alive in potrs() also depending on potrs_ev directly, as trsm_batch.cpp does:

host_task_events.push_back(potrs_ev); sycl::event args_ev = dpnp::utils::keep_args_alive( exec_q, {a_array, b_array}, host_task_events);

Done: the clean-up host task is now submitted only when a scratchpad was allocated, and keep_args_alive also depends on potrs_ev. With oneMKL 2026.2, potrs_scratchpad_size returns 0 for all four types (n up to 4096, nrhs up to 512) on CPU and GPU, over both OpenCL and Level Zero, so every call was paying for that extra host task.

oneMKL reports a zero potrs scratchpad, so the clean-up only freed a nullptr. keep_args_alive now depends on the potrs event directly.

antonwolfy added this to the 0.22.0 release milestone Oct 20, 2026

antonwolfy self-assigned this Oct 20, 2026

antonwolfy force-pushed the cho-factor-cho-solve branch from 410539a to 6b3035e Compare October 30, 2026 14:24

antonwolfy marked this pull request as ready for review October 30, 2026 16:31

antonwolfy requested review from ndgrigorian and vlad-perevezentsev as code owners October 30, 2026 16:31

antonwolfy requested a balanced review from Copilot October 30, 2026 17:17

Copilot started reviewing on behalf of antonwolfy October 30, 2026 17:18

View session

Copilot AI reviewed Oct 30, 2026

View reviewed changes

Comment thread dpnp/scipy/linalg/_utils.py Resolved

Comment thread dpnp/scipy/linalg/_utils.py Resolved

ndgrigorian reviewed Nov 3, 2026

View reviewed changes

Comment thread dpnp/backend/extensions/lapack/potrs.cpp Resolved

Comment thread dpnp/scipy/linalg/_utils.py Resolved

Comment thread dpnp/scipy/linalg/_utils.py Resolved

Comment thread dpnp/libtensor/extensions/blas/trsm_batch.cpp Resolved

antonwolfy requested a balanced review from Copilot November 3, 2026 14:02

Copilot started reviewing on behalf of antonwolfy November 3, 2026 14:03

View session

Copilot AI reviewed Nov 3, 2026

View reviewed changes

antonwolfy added 13 commits November 5, 2026 16:18

Skip dtype check in cho_solve CuPy test on devices without fp64

a8c86f6

Reference oneMKL issue MKLD-20341 for the batched cho_solve workaround

8711efe

Complete the Raises and Warnings sections of cho_factor and cho_solve

63d1263

- List the `ValueError` cases of `cho_factor` next to `LinAlgError` - Add a Raises section to `cho_solve` - Mention the sync to read a dpnp array `lower` in `cho_solve`

Say the potrf_batch issue behind a test ordering was seen under WSL only

e5409f8

antonwolfy added 5 commits November 5, 2026 16:18

Do not broadcast a factor with size-1 batch dimensions in cho_solve

1d1fe14

A `c` of shape (1, n, n) or (1, 1, n, n) with a batched right-hand side was broadcast to the batch and copied per batch element; solve with it as with a 2-D factor instead

Reuse b in cho_solve only if it has the USM type of the result

b254743

With `overwrite_b=True`, a right-hand side whose USM type differs from the coerced one of `c` and `b` was solved in place, so the solution kept the USM type of `b`. Copy it in that case, as without `overwrite_b`.

ndgrigorian reviewed Nov 5, 2026

View reviewed changes

antonwolfy force-pushed the cho-factor-cho-solve branch from 38ec63c to 20b0f45 Compare November 5, 2026 15:50

antonwolfy and others added 2 commits November 6, 2026 12:28

Skip the potrs clean-up host task when no scratchpad is allocated

e08e8e1

oneMKL reports a zero potrs scratchpad, so the clean-up only freed a nullptr. keep_args_alive now depends on the potrs event directly.

Merge branch 'master' into cho-factor-cho-solve

a17fa60

Navigation Menu

Repository navigation

Add `cho_factor` and `cho_solve` to `dpnp.scipy.linalg` - #3116

Add `cho_factor` and `cho_solve` to `dpnp.scipy.linalg`#3116
antonwolfy wants to merge 20 commits into
masterfrom
cho-factor-cho-solve

antonwolfy commented Oct 20, 2026 •

edited

Loading

Uh oh!

github-actions Bot commented Oct 20, 2026

Uh oh!

github-actions Bot commented Oct 20, 2026 •

edited

Loading

Uh oh!

coveralls commented Oct 21, 2026 •

edited

Loading

Uh oh!

Copilot AI left a comment

Uh oh!

Uh oh!

Uh oh!

antonwolfy commented Oct 30, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Copilot AI left a comment

Uh oh!

ndgrigorian Nov 5, 2026 •

edited

Loading

Uh oh!

antonwolfy Nov 6, 2026

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Conversation

antonwolfy commented Oct 20, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Implementation

Differences from SciPy

Other changes

Testing

Performance

Uh oh!

github-actions Bot commented Oct 20, 2026

Uh oh!

github-actions Bot commented Oct 20, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

coveralls commented Oct 21, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

🟡 Changes recommended

Uh oh!

Uh oh!

Uh oh!

antonwolfy commented Oct 30, 2026

Ridge benchmark with the new functions

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

🔵 Needs a closer look

Uh oh!

ndgrigorian Nov 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

antonwolfy Nov 6, 2026

Choose a reason for hiding this comment

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

antonwolfy commented Oct 20, 2026 •

edited

Loading

github-actions Bot commented Oct 20, 2026 •

edited

Loading

coveralls commented Oct 21, 2026 •

edited

Loading

`Ridge` benchmark with the new functions

ndgrigorian Nov 5, 2026 •

edited

Loading