Specialize Diagonal products and dot for adjoint/transpose of sparse matrices - #770
Merged
Merged
Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #770 +/- ##
==========================================
+ Coverage 92.59% 92.61% +0.02%
==========================================
Files 12 12
Lines 8780 8848 +68
==========================================
+ Hits 8130 8195 +65
- Misses 650 653 +3 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Contributor
There was a problem hiding this comment.
🟢 Approval recommended
The implementations are consistent with existing sparse kernels and comprehensively tested.
Pull request overview
Adds sparse specializations for lazy adjoint/transpose operations, resolving performance regressions in diagonal multiplication and Frobenius dot products.
Changes:
- Materializes sparse wrappers before diagonal scaling.
- Implements a linear-time cursor-based sparse
dotkernel. - Adds correctness, edge-case, and performance tests.
File summaries
| File | Description |
|---|---|
src/linalg.jl |
Adds specialized multiplication and dot methods. |
test/linalg.jl |
Tests diagonal products and performance. |
test/linalg_products.jl |
Tests sparse wrapper dot products and edge cases. |
Review details
- Files reviewed: 3/3 changed files
- Comments generated: 0
- Review effort level: Balanced
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
ViralBShah
marked this pull request as draft
September 9, 2026 13:09
ViralBShah
marked this pull request as ready for review
September 10, 2026 09:25
ViralBShah
added this pull request to stack #815
September 11, 2026 10:39
…rse matrices Fixes #619: `A' * D` and `D * A'` for a sparse `A` and `Diagonal` `D` fell through to the generic `AbstractMatrix` product, ~300x slower than `A * D` on Julia 1.11. Materialize the adjoint (O(nnz)) and reuse the existing CSC-times-Diagonal kernels, mirroring how `A' * B` is handled for sparse `B`. Fixes #627: `dot(A', B)` for sparse `A`, `B` walked the stored entries of `B` and did a binary search into `A'` for each, ~50x slower than `dot(copy(A'), B)` on Julia 1.11. Add a merge that walks the columns of `B` in order while keeping one cursor per column of `parent(A)`, so it runs in O(nnz(A) + nnz(B) + n) time with O(n) extra memory and no O(nnz) temporary. `dot(B, A')` reaches the same kernel through the existing `conj(dot(A', B))`. Tests cover both wrappers, real and complex eltypes, mixed eltypes, stored zeros, empty columns, non-square shapes, dimension errors, and timing guards. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LgBHUw9Hp7YW5ub29B4R5y - Add `mul!` kernels for adjoint/transpose of a sparse matrix with a `Diagonal` - Address review: search when the other `dot` operand is far denser, alias check by storage
…the `dot` walk The adjoint kernels formed the transpose in the destination eltype before scaling, so a `Float32` or `Int` destination lost precision or threw where dense and the CSC kernels compute the product first. The direct path now requires the destination eltype to be the product's, and `alpha == 0` ignores `A` in all four `Diagonal` kernels, as for dense. A fixed destination whose pattern contains the product's is filled through the merge branch instead of failing in `copyinds!`, and one lacking an entry throws an `ArgumentError` with the destination untouched. The `dot` walk now picks the operand with fewer stored entries plus columns, since the walk visits every column. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…e numbers from comments Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…ypes of `C` and `A` differ Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
ViralBShah
force-pushed
the
vs/adjoint-diagonal-dot
branch
from
September 19, 2026 13:28
a7d3611 to
da8de50
Compare
ViralBShah
added a commit
that referenced
this pull request
Sep 19, 2026
Stacked on #770. Every sparse `dot` seeded its accumulator with `zero` of the entry types, either directly or through `dot(zero(T1), zero(T2))`, so any sparse operand with matrix-valued entries failed with `no method matching zero(::Type{Matrix{Float64}})`, while the dense `dot` in LinearAlgebra handles such entries. Seed with the zero of the *result* type instead, `zero(promote_op(dot, Ts...))`, which is a scalar even when the entries themselves have no `zero`. For number entries this is the same value as before. Covered: the Frobenius `dot` of two CSC matrices, of a dense or wrapped matrix with a CSC matrix, and of a lazy adjoint/transpose with a CSC matrix; sparse-vector `dot` with dense and sparse vectors; and the three-argument forms with a CSC matrix or a `Diagonal`. The sparse-vector/CSC/sparse-vector form also computed `dot(x, A) * y` per column, which is only `dot(x, A, y)` for scalar entries, so it now applies `dot(x, a, y)` per matching entry; on `ComplexF64` at `n=5000`, density 0.01 and 0.3, that is 2.25 ms against 2.83 ms before. The two `Diagonal` forms no longer call `first` on the operands, so they also work for empty vectors. The multiplication-count assertions from #770 lose the one multiplication the old seed performed. Symmetric/Hermitian three-argument `dot` with matrix-valued entries is unchanged: its kernels sum matrix products before the final `dot`, which needs a different seed. 🤖 Generated with [Claude Code](https://claude.com/claude-code) https://claude.ai/code/session_015VF52nADauBDqAQaUHjoNV Co-authored-by: Viral B. Shah <ViralBShah@users.noreply.github.com>
This was referenced Sep 19, 2026
ViralBShah
added a commit
that referenced
this pull request
Sep 19, 2026
Follow-up to #770 for #627. `dot` of an `Adjoint` with a `Transpose` of a sparse matrix (either order) still reaches LinearAlgebra's generic O(mn) `getindex` loop, because LinearAlgebra only strips matching wrappers. ```julia using SparseArrays, LinearAlgebra A = sprand(ComplexF64, 10^4, 10^4, 1e-3); C = sprand(ComplexF64, 10^4, 10^4, 1e-3) dot(A', transpose(C)) ``` Both wrappers transpose positions, so the stored entries of the parents line up column by column and only the elementwise operation differs. The column-merge loop of `dot(::AbstractSparseMatrixCSC, ::AbstractSparseMatrixCSC)` moves into `_dot_walk(f, A, B)`, and the two new methods call it on the parents with `dot(adjoint(a), transpose(b))` and its mirror. This touches dispatch: two new `dot` methods, no new ambiguities with LinearAlgebra. Minimum of 10 runs, nightly 1.14.0-DEV.3300, 0 allocations before and after: | N, density | `A' ⋅ transpose(C)` main | PR | `A ⋅ C` | |---|---|---|---| | 1000, 0.1 | 36.6 ms | 0.66 ms | 0.66 ms | | 10⁴, 1e-3 | 1469 ms | 0.65 ms | 0.65 ms | Tested against dense with complex entries, stored zeros and a non-square shape, plus a multiplication count proving the kernel is reached; `test/linalg_products.jl` passes on nightly. Not covered: `Symmetric`/view operands against a lazy adjoint, which remain on the fallback. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Fable 5.1 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Fixes #619 and #627. Both are cases where a lazy
Adjoint/Transposeof aSparseMatrixCSCfell through to a generic method that ignores sparsity.#619:
mul!of an adjoint/transpose of a sparse matrix with aDiagonalLinearAlgebra's
Diagonalkernel visits every element of the destination, somul!(C, S', D),mul!(C, D, S')and their 5-argument forms were O(m·n). New 5-argumentmul!methods form the adjoint directly inCwith onehalfperm!and scale it in place, falling back to a materialized copy whenbetais nonzero or whenCshares storage with the parent, is fixed, or has another index type.*reaches these throughmatop_dest, which now hands the adjoint of a sparse matrix an empty writable destination (with a size hint) and a fixed sparse matrix a fixed destination with its own structure. As a resultD * FandF' * Dfor a fixedFno longer fail with "Can't change ReadOnly";F * DandD * F'already worked.On nightly (
N=1000, density 0.1, min of 7):mul!(C, A', D)mul!(C, D, A')mul!(C, A', D, 2, 3)A' * D,D * A'*is unchanged on nightly, where the generic path already materializes the adjoint; the explicit methods make it independent of LinearAlgebra internals.#627:
dot(A', B)for sparseA,BThe existing wrapper method walks the stored entries of
Band does a binary search intoA'for each one. The new method walks the sparser operand and keeps one cursor per column of the other, so it runs in O(nnz(A) + nnz(B) + n) time with O(n) extra memory and no O(nnz) temporary. Since the cursors sweep every stored entry of the other operand, it switches to a binary search per entry once the other operand holds over 32x more entries (or columns) than the walked one; the measured crossover is a ratio of 20-50. It returns early when either operand has no stored entries, so two nearly empty operands with 10^7 rows cost nothing.dot(B, A')reaches the same kernel through the existingconj(dot(A', B)).N=1000, min of 7)dot(A', B), both density 0.1dot(copy(A'), B)dot(P', Bd), nnz(P)=2e3, nnz(Bd)=5e5dot(Bd', P)Semantics shared with dense and with the existing CSC kernels
Review of the first version found the direct adjoint path diverging from
mul!on dense arrays and from the plain CSC kernels, so the fourDiagonalkernels now agree with both:Float32orIntdestination goes through the materialized copy and the CSC kernel, somul!(spzeros(Float32,1,1), sparse([1e40;;])', Diagonal([1e-40]))is1f0and anIntdestination takes0.5 * 2.0without anInexactError, as onmain.alpha == 0ignoresAandbeta == 0ignoresC, so nonfinite entries in either do not leak into the result. This was already wrong in the plain CSC kernels onmain, which returnedNaNformul!(C, sparse([Inf;;]), D, 0, 0); they are fixed here too so that the adjoint and plain paths agree.mainonly handled the adjoint case, through the generic kernel, and threw fromcopyinds!for a plain operand. A fixed destination lacking an entry throws anArgumentErrornaming the entry, with the destination untouched, wheremainthrew aMethodErrorfrominsert!.dotwalk picks the operand with fewer stored entries plus columns, because the walk visits every column of the operand it walks. A 1×10⁶ parent with one stored entry against a 10⁶×1 operand with two went from 0.2 µs onmainto 627 µs with the first version, and is back to 0.2 µs.Dispatch is touched: three
matop_destmethods and themul!anddotmethods above are new. Not a backport candidate, since it adds methods.Tests
Wall-clock guards are replaced by a multiplication-counting eltype (
test/util/mulcount.jl): the kernels perform exactlynnzmultiplications, anddotonly multiplies where both operands store an entry, from either side of the walk. Also covered: fixed operands, the 3- and 5-argumentmul!forms, a destination that aliases or shares storage with the parent, another index type, bothdotbranches, stored zeros, empty columns, non-square shapes and dimension errors.Test.detect_ambiguities(SparseArrays; recursive=true)is empty. The dense-consistency cases above are tested for the plain, adjoint and transposed operand in both orders;test/linalg.jl,test/linalg_products.jlandtest/fixed.jlpass on nightly.🤖 Generated with Claude Code
https://claude.ai/code/session_015VF52nADauBDqAQaUHjoNV