Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

stdlib: faster kronecker product between hermitian and symmetric matrices #53186

Merged
merged 10 commits into from
Apr 18, 2024
Merged

stdlib: faster kronecker product between hermitian and symmetric matrices #53186

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
c22b47a
faster kronecker product between hermitian matrices
araujoms Feb 5, 2024
e8cc43e
separate out kernel and add methods for triangular matrices
araujoms Feb 7, 2024
e758fa9
slightly faster triangular kernel
araujoms Feb 8, 2024
674fcf7
Merge branch 'master' into hermkron
araujoms Feb 9, 2024
6a32c6e
add tests for kron with partially initialized matrices
araujoms Feb 26, 2024
d3f5ec6
Merge branch 'master' into hermkron
araujoms Feb 26, 2024
8a60a4a
specialize methods to numeric matrices
araujoms Mar 3, 2024
f1326e3
Merge branch 'master' into hermkron
araujoms Mar 3, 2024
c356ede
Merge branch 'master' into hermkron
araujoms Mar 21, 2024
485eadf
Merge branch 'master' into hermkron
oscardssmith Apr 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
115 changes: 115 additions & 0 deletions stdlib/LinearAlgebra/src/symmetric.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -496,6 +496,121 @@ for (T, trans, real) in [(:Symmetric, :transpose, :identity), (:(Hermitian{<:Uni
end
end

function kron(A::Hermitian{T}, B::Hermitian{S}) where {T<:Union{Real,Complex},S<:Union{Real,Complex}}
resultuplo = A.uplo == 'U' || B.uplo == 'U' ? :U : :L
R = Hermitian(Matrix{promote_op(*, T, S)}(undef, _kronsize(A, B)), resultuplo)
return kron!(R, A, B)
end

function kron(A::Symmetric{T}, B::Symmetric{S}) where {T,S}
resultuplo = A.uplo == 'U' || B.uplo == 'U' ? :U : :L
R = Symmetric(Matrix{promote_op(*, T, S)}(undef, _kronsize(A, B)), resultuplo)
return kron!(R, A, B)
end

for (T, conj, real) in [(:Symmetric, :identity, :identity), (:(Hermitian{<:Union{Real,Complex}}), :conj, :real)]
Copy link
Contributor

@jishnub jishnub Feb 5, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The method should probably be restricted to Hermitian{<:Union{Real,Complex}, <StridedMatrix} and similarly for Symmetric, and a fallback implementation of the form kron!(C.data, A.data, B.data) may be added. This will ensure that sparse arrays (or other special array types) remain performant, if they provide a kron! implementation.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please forgive me my ignorance as I'm rather new to Julia. You are saying I should change the function declarations to

function kron(A::Hermitian{T, <:StridedMatrix}, B::Hermitian{S, <:StridedMatrix}) where {T<:Union{Real,Complex},S<:Union{Real,Complex}}
function kron(A::Symmetric{T, <:StridedMatrix}, B::Symmetric{S, <:StridedMatrix}) where {T,S}

in order to act only on dense matrices, not sparse ones. Fine, I get this. But then I should also give a fallback implementation, that only acts if the sparse type doesn't provide their own kron!? How? I thought the original function declaration was already doing this, as a more specialized function declaration for a specific type would take precedence.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm recommending

function kron!(C::Hermitian{TC, <:StridedMatrix}, A::Hermitian{T, <:StridedMatrix}, B::Hermitian{S, <:StridedMatrix}) where {T<:Union{Real,Complex},S<:Union{Real,Complex},TC<:Union{Real,Complex}}

as the signature of the method that is added here, and add another method

kron!(C::Hermitian{<:Union{Real,Complex}}, A::Hermitian{<:Union{Real,Complex}}, B::Hermitian{<:Union{Real,Complex}}) = Hermitian(kron!(C.data, A.data, B.data))

If the parent has an efficient kron! implementation, this would not loop over the entire array. This way, custom array types only need to specialize kron! for their own matrix types, and not for wrappers such as Hermitian{<:Union{Real,Complex}, MyArray}.

You're right that the original method already provides a more efficient implementation than the fully generic dense-array implementation that loops over the entire A and B. However, forwarding kron! to the parents allows us to access even more optimized methods that packages may provide. Ultimately, it's a tradeoff. You're right that a more specialized method kron!(::Hermitian{<:Union{Real,Complex}}, ::Hermitian{<:Union{Real,Complex}, MyArray}, ::Hermitian{<:Union{Real,Complex}, MyArray}) would take precedence, but it's better if packages don't need to define such methods for wrappers around their arrays, as that greatly increases the total number of methods that need to be defined. If we forward the operation to the parent, then the parent does not need to know about Hermitian at all.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see, the issue is that a generic Kronecker product over sparse arrays will still be vastly more efficient for Hermitian matrices than my dense version, so we don't want to capture that. Of course, such a sparse type can also implement its own specialization for Hermitian matrices, but that's a different issue.

The best I could come up with was

for T in [:Symmetric, :(Hermitian{<:Union{Real,Complex}})]
    @eval begin
        function kron!(C::$T, A::$T, B::$T)
            wrapper = Base.typename($T).wrapper
            return wrapper(kron!(C.data, A.data, B.data))
        end
    end
end

for the fallback and

for (T, conj, real) in [(:(Symmetric{<:Any,<:StridedMatrix}), :identity, :identity), (:(Hermitian{<:Union{Real,Complex}, <:StridedMatrix}), :conj, :real)]
    @eval begin
        function kron!(C::$T, A::$T, B::$T)

for the signature of the method itself. It's rather hideous. Is there a better solution?

Also, doesn't exactly the same issue applies to dot? Is there any package actually implementing faster versions of dot and kron for sparse arrays?

Copy link
Contributor

@jishnub jishnub Feb 6, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I take my comment back, as kron forwarded to the parent doesn't produce the same result:

julia> H1 = Hermitian(rand(2,2));

julia> H2 = Hermitian(rand(2,2));

julia> kron(H1, H2)
4×4 Matrix{Float64}:
 0.343295  0.436362  0.106336   0.135163
 0.436362  0.717831  0.135163   0.222349
 0.106336  0.135163  0.0616077  0.0783094
 0.135163  0.222349  0.0783094  0.128822

julia> Hermitian(kron(H1.data, H2.data))
4×4 Hermitian{Float64, Matrix{Float64}}:
 0.343295  0.436362  0.106336   0.135163
 0.436362  0.717831  0.219534   0.222349
 0.106336  0.219534  0.0616077  0.0783094
 0.135163  0.222349  0.0783094  0.128822

I think the current implementation is fine.

Copy link
Contributor

@jishnub jishnub Feb 16, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's just do the symmetric/hermitian for now

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤦 I wasted a lot of time doing the triangular case only because @dkarrasch asked. I have no use for it. Next time please only ask for features you actually want.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh. I'm in no way suggesting that we abandon the work if you've implemented it already. My impression was otherwise. Please continue with your work.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah ok. Well, it's done already, I've already pushed it.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well?

@eval begin
function kron!(C::$T, A::$T, B::$T)
size(C) == _kronsize(A, B) || throw(DimensionMismatch("kron!"))
if ((A.uplo == 'U' || B.uplo == 'U') && C.uplo != 'U') || ((A.uplo == 'L' && B.uplo == 'L') && C.uplo != 'L')
throw(ArgumentError("C.uplo must match A.uplo and B.uplo, got $(C.uplo) $(A.uplo) $(B.uplo)"))
end
n_A = size(A, 1)
n_B = size(B, 1)
@inbounds if A.uplo == 'U' && B.uplo == 'U'
for j = 1:n_A
jnB = (j - 1) * n_B
for i = 1:(j-1)
Aij = A.data[i, j]
inB = (i - 1) * n_B
for l = 1:n_B
for k = 1:(l-1)
C.data[inB+k, jnB+l] = Aij * B.data[k, l]
C.data[inB+l, jnB+k] = Aij * $conj(B.data[k, l])
end
C.data[inB+l, jnB+l] = Aij * $real(B[l, l])
end
end
Ajj = $real(A[j, j])
for l = 1:n_B
for k = 1:(l-1)
C.data[jnB+k, jnB+l] = Ajj * B.data[k, l]
end
C.data[jnB+l, jnB+l] = Ajj * $real(B[l, l])
end
end
elseif A.uplo == 'U' && B.uplo == 'L'
for j = 1:n_A
jnB = (j - 1) * n_B
for i = 1:(j-1)
Aij = A.data[i, j]
inB = (i - 1) * n_B
for l = 1:n_B
C.data[inB+l, jnB+l] = Aij * $real(B[l, l])
for k = (l+1):n_B
C.data[inB+l, jnB+k] = Aij * $conj(B.data[k, l])
C.data[inB+k, jnB+l] = Aij * B.data[k, l]
end
end
end
Ajj = $real(A[j, j])
for l = 1:n_B
C.data[jnB+l, jnB+l] = Ajj * $real(B[l, l])
for k = (l+1):n_B
C.data[jnB+l, jnB+k] = Ajj * $conj(B.data[k, l])
end
end
end
elseif A.uplo == 'L' && B.uplo == 'U'
for j = 1:n_A
jnB = (j - 1) * n_B
Ajj = $real(A[j, j])
for l = 1:n_B
for k = 1:(l-1)
C.data[jnB+k, jnB+l] = Ajj * B.data[k, l]
end
C.data[jnB+l, jnB+l] = Ajj * $real(B[l, l])
end
for i = (j+1):n_A
conjAij = $conj(A.data[i, j])
inB = (i - 1) * n_B
for l = 1:n_B
for k = 1:(l-1)
C.data[jnB+k, inB+l] = conjAij * B.data[k, l]
C.data[jnB+l, inB+k] = conjAij * $conj(B.data[k, l])
end
C.data[jnB+l, inB+l] = conjAij * $real(B[l, l])
end
end
end
else #if A.uplo == 'L' && B.uplo == 'L'
for j = 1:n_A
jnB = (j - 1) * n_B
Ajj = $real(A[j, j])
for l = 1:n_B
C.data[jnB+l, jnB+l] = Ajj * $real(B[l, l])
for k = (l+1):n_B
C.data[jnB+k, jnB+l] = Ajj * B.data[k, l]
end
end
for i = (j+1):n_A
Aij = A.data[i, j]
inB = (i - 1) * n_B
for l = 1:n_B
C.data[inB+l, jnB+l] = Aij * $real(B[l, l])
for k = (l+1):n_B
C.data[inB+k, jnB+l] = Aij * B.data[k, l]
C.data[inB+l, jnB+k] = Aij * $conj(B.data[k, l])
end
end
end
end
end
return C
end
end
end

(-)(A::Symmetric) = Symmetric(-A.data, sym_uplo(A.uplo))
(-)(A::Hermitian) = Hermitian(-A.data, sym_uplo(A.uplo))

Expand Down
32 changes: 32 additions & 0 deletions stdlib/LinearAlgebra/test/symmetric.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -467,6 +467,28 @@ end
@test dot(symblockml, symblockml) ≈ dot(msymblockml, msymblockml)
end
end

@testset "kronecker product of symmetric and Hermitian matrices" begin
for mtype in (Symmetric, Hermitian)
symau = mtype(a, :U)
symal = mtype(a, :L)
msymau = Matrix(symau)
msymal = Matrix(symal)
for eltyc in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
creal = randn(n, n)/2
cimag = randn(n, n)/2
c = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(creal, cimag) : creal)
symcu = mtype(c, :U)
symcl = mtype(c, :L)
msymcu = Matrix(symcu)
msymcl = Matrix(symcl)
@test kron(symau, symcu) ≈ kron(msymau, msymcu)
@test kron(symau, symcl) ≈ kron(msymau, msymcl)
@test kron(symal, symcu) ≈ kron(msymal, msymcu)
@test kron(symal, symcl) ≈ kron(msymal, msymcl)
end
end
end
end
end

Expand All @@ -481,6 +503,16 @@ end
@test dot(A, B) ≈ dot(Symmetric(A), Symmetric(B))
end

# let's make sure the analogous bug will not show up with kronecker products
@testset "kron Hermitian quaternion #52318" begin
A, B = [Quaternion.(randn(3,3), randn(3, 3), randn(3, 3), randn(3,3)) |> t -> t + t' for i in 1:2]
@test A == Hermitian(A) && B == Hermitian(B)
@test kron(A, B) ≈ kron(Hermitian(A), Hermitian(B))
A, B = [Quaternion.(randn(3,3), randn(3, 3), randn(3, 3), randn(3,3)) |> t -> t + transpose(t) for i in 1:2]
@test A == Symmetric(A) && B == Symmetric(B)
@test kron(A, B) ≈ kron(Symmetric(A), Symmetric(B))
end

#Issue #7647: test xsyevr, xheevr, xstevr drivers.
@testset "Eigenvalues in interval for $(typeof(Mi7647))" for Mi7647 in
(Symmetric(diagm(0 => 1.0:3.0)),
Expand Down