From dc76e149b585cd5c1d99f7577e209d1cadb983b3 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Tue, 8 Sep 2026 05:25:43 -0400 Subject: [PATCH] Add `isequal` for adjoint/transpose of sparse matrices #766 added stored-entry-only `isequal` for `SparseMatrixCSC` and for the `Adjoint`/`Transpose` wrappers of sparse vectors, but not for the wrappers of sparse matrices, which still fell back to the generic elementwise `AbstractArray` method. On two 10^5 x 10^5 matrices with one stored entry each, `isequal(A', B')` took ~14s; it now matches `A' == B'` at ~0.2ms. Parameterize `nzeq` on the elementwise predicate, mirror the `==` methods for `Adjoint`/`Transpose` matrices with `isequal`, and loosen the indexed argument of `nzeq` to `AbstractMatrix`. The latter also fixes a pre-existing `MethodError` in `A' == transpose(B)` for complex `A`, where the adjoint of a transpose does not collapse and yields a nested `Adjoint{<:Any,<:Transpose}`. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01LgBHUw9Hp7YW5ub29B4R5y --- src/sparsematrix.jl | 33 ++++++++++++++++++++++----------- test/sparsematrix_ops.jl | 35 +++++++++++++++++++++++++++++++++++ 2 files changed, 57 insertions(+), 11 deletions(-) diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index d4b19c64..51bdb5e8 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -2378,41 +2378,52 @@ Base.isequal(A1::AbstractSparseMatrixCSC, A2::AbstractSparseMatrixCSC) = _iseq(i ## Explicit efficient comparisons with transposed arrays # Check whether all nonzero elements of A are equal to the respective elements in B -function nzeq(A::AbstractSparseMatrixCSC, B::AbstractSparseMatrixCSCInclAdjointAndTranspose) +# under the elementwise predicate `eq` (`==` or `isequal`) +function nzeq(eq::F, A::AbstractSparseMatrixCSC, B::AbstractMatrix) where {F} @inbounds for j in axes(A,2) for k in nzrange(A, j) i = rowvals(A)[k] val = nonzeros(A)[k] - val ≠ B[i,j] && return false + eq(val, B[i,j]) || return false end end return true end # Peel off `Adjoint` and `Transpose` from first argument -nzeq(A::Adjoint{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}, - B::AbstractSparseMatrixCSCInclAdjointAndTranspose) = - nzeq(A', B') -nzeq(A::Transpose{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}, - B::AbstractSparseMatrixCSCInclAdjointAndTranspose) = - nzeq(transpose(A), transpose(B)) +# `B` may be a nested wrapper such as `Adjoint{<:Any,<:Transpose}` (from `A' == transpose(B)`), +# hence the loose `AbstractMatrix` bound: `B` is only ever indexed +nzeq(eq::F, A::Adjoint{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}, + B::AbstractMatrix) where {F} = + nzeq(eq, A', B') +nzeq(eq::F, A::Transpose{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}, + B::AbstractMatrix) where {F} = + nzeq(eq, transpose(A), transpose(B)) # Compare by walking both matrices # (We could further optimize the case `AbstractSparseMatrixCSC == # Adjoint(Transpose(AbstractSparseMatrixCSC))` more efficiently, i.e. # the case where the RHS is both adjoint and transposed, i.e. where it # is in CSC format again.) -function ==(A::AbstractSparseMatrixCSC, - B::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}) +function _iseq(eq::F, A::AbstractSparseMatrixCSC, + B::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}) where {F} # Different sizes are always different size(A) ≠ size(B) && return false # Compare nonzero elements - return nzeq(A, B) && nzeq(B, A) + return nzeq(eq, A, B) && nzeq(eq, B, A) end +==(A::AbstractSparseMatrixCSC, B::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}) = + _iseq(==, A, B) +Base.isequal(A::AbstractSparseMatrixCSC, B::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}) = + _iseq(isequal, A, B) # Peel off `Adjoint` and `Transpose` from first argument ==(A::Adjoint{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}, B::AbstractSparseMatrixCSCInclAdjointAndTranspose) = A' == B' ==(A::Transpose{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}, B::AbstractSparseMatrixCSCInclAdjointAndTranspose) = transpose(A) == transpose(B) +Base.isequal(A::Adjoint{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}, B::AbstractSparseMatrixCSCInclAdjointAndTranspose) = + isequal(A', B') +Base.isequal(A::Transpose{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}, B::AbstractSparseMatrixCSCInclAdjointAndTranspose) = + isequal(transpose(A), transpose(B)) ## Reductions diff --git a/test/sparsematrix_ops.jl b/test/sparsematrix_ops.jl index 46189e42..a0d941a3 100644 --- a/test/sparsematrix_ops.jl +++ b/test/sparsematrix_ops.jl @@ -74,6 +74,41 @@ end end end +@testset "isequal for adjoint/transpose of sparse matrices" begin + n = 10^5 + A = spzeros(n, n); A[1, 1] = 1 + B = copy(A) + for (L, R) in ((A', B'), (transpose(A), transpose(B)), (A, B'), (A', B), + (A, transpose(B)), (transpose(A), B), (A', transpose(B))) + @test isequal(L, R) + @test @elapsed(isequal(L, R)) < 0.1 + end + A[1, 2] = 1; B[2, 1] = 1 + @test isequal(A, B') && isequal(A', B) && !isequal(A', B') && !isequal(A, B) + @test !isequal(spzeros(2, 3)', spzeros(2, 3)) + # adjoint vs transpose of a complex matrix nests wrappers (`Adjoint{<:Any,<:Transpose}`) + C = sparse([1, 2], [2, 3], [1.0im, 2.0], 3, 3) + @test C' == transpose(conj(C)) && isequal(C', transpose(conj(C))) + @test C' != transpose(C) && !isequal(C', transpose(C)) + # isequal semantics for NaN, signed zeros and conjugation must match dense arrays + X = sparse([1, 2, 3, 1], [1, 1, 2, 3], [NaN, -0.0, 2.0, 1.0im], 3, 3) + for Y in (sparse([1, 2, 3, 1], [1, 1, 2, 3], [NaN, -0.0, 2.0, 1.0im], 3, 3), + sparse([1, 2, 3, 1], [1, 1, 2, 3], [NaN, -0.0, 2.0, -1.0im], 3, 3), + sparse([1, 1, 2, 3], [1, 2, 3, 1], [NaN, -0.0, 2.0, 1.0im], 3, 3), + sparse([1, 1, 2, 3], [1, 2, 3, 1], [NaN, -0.0, 2.0, -1.0im], 3, 3), + sparse([1, 3, 1], [1, 2, 3], [NaN, 2.0, 1.0im], 3, 3), + sparse([1, 2, 3, 1], [1, 1, 2, 3], [NaN, 0.0, 2.0, 1.0im], 3, 3), + sparse([1, 2, 3, 1, 3], [1, 1, 2, 3, 3], [NaN, -0.0, 2.0, 1.0im, 0.0], 3, 3), + sparse([1, 2, 3, 1], [1, 1, 2, 3], [1.0, -0.0, 2.0, 1.0im], 3, 3)) + for (L, R) in ((X', Y'), (transpose(X), transpose(Y)), (X, Y'), (X', Y), + (X, transpose(Y)), (transpose(X), Y), (X', transpose(Y))) + @test isequal(L, R) == isequal(Matrix(L), Matrix(R)) + @test isequal(R, L) == isequal(Matrix(R), Matrix(L)) + @test (L == R) == (Matrix(L) == Matrix(R)) + end + end +end + @testset "iszero specialization for SparseMatrixCSC" begin @test !iszero(sparse(I, 3, 3)) # test failure @test iszero(spzeros(3, 3)) # test success with no stored entries