Specialize Diagonal products and dot for adjoint/transpose of sparse matrices - #770
Open
ViralBShah wants to merge 1 commit into
Open
Specialize Diagonal products and dot for adjoint/transpose of sparse matrices#770ViralBShah wants to merge 1 commit into
Diagonal products and dot for adjoint/transpose of sparse matrices#770ViralBShah wants to merge 1 commit into
Conversation
…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
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #770 +/- ##
==========================================
+ Coverage 84.19% 84.23% +0.04%
==========================================
Files 13 13
Lines 9403 9428 +25
==========================================
+ Hits 7917 7942 +25
Misses 1486 1486 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
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:
A' * DandD * A'withD::DiagonalMaterialize the adjoint (O(nnz)) and reuse the existing CSC-times-
Diagonalkernels, mirroring howA' * Bis already handled for sparseB. The result stays aSparseMatrixCSC.On released Julia 1.11 the generic fallback is ~300x slower than
A * D(N=1000, density 0.1):A * DA' * DD * A'On current nightly the generic path already materializes the adjoint, so the gain there is only ~10%; the explicit methods make the fast path independent of LinearAlgebra internals and are worth backporting.
#627:
dot(A', B)for sparseA,BThe existing wrapper method walks the stored entries of
Band does a binary search intoA'for each one. Replace it with a merge: walk the columns ofBin order and keep one cursor per column ofparent(A). Since the row index being sought in each column of the parent is nondecreasing, the cursors only move forward, giving 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 existingconj(dot(A', B)).N=1000, density 0.1)dot(copy(A'), A)dot(A', A)Tests cover both wrappers, real and complex eltypes, mixed eltypes, stored zeros, empty columns, non-square shapes, dimension errors, and timing guards. Full test suite passes on nightly and
Test.detect_ambiguities(SparseArrays)remains empty.🤖 Generated with Claude Code
https://claude.ai/code/session_01LgBHUw9Hp7YW5ub29B4R5y