From 6cef73ce6766cfa1a1f57b66c3ca19f875781748 Mon Sep 17 00:00:00 2001 From: GiggleLiu Date: Sat, 7 Oct 2023 10:19:44 +0800 Subject: [PATCH 1/4] fast overlap and fast add --- src/auto_diff.jl | 20 ++++++++++++++++++-- src/utils.jl | 20 +++++++++++++++++--- test/runtests.jl | 18 ++++++++++++++++++ 3 files changed, 53 insertions(+), 5 deletions(-) diff --git a/src/auto_diff.jl b/src/auto_diff.jl index eaca00c..b40863b 100644 --- a/src/auto_diff.jl +++ b/src/auto_diff.jl @@ -38,11 +38,27 @@ function Yao.AD.backward_params!(st::Tuple{<:MajoranaReg,<:MajoranaReg}, out, outδ = st ham = Yao.AD.generator(block) majoranaham = yaoham2majoranasquares(ham) - g = outδ.state ⋅ (majoranaham * out.state) / 4 - pushfirst!(collector, g) + # multiply majoranaham with out.state and inner product with outδ.state + # i.e. g = outδ.state ⋅ (majoranaham * out.state) / 4 + pushfirst!(collector, fast_overlap(outδ.state, majoranaham, out.state)/4) return nothing end +# compute tr(y' * A * x) +function fast_overlap(y::AbstractVecOrMat{T}, A::SparseMatrixCSC{T}, x::AbstractVecOrMat{T}) where T + @assert size(x, 1) == size(A, 2) && size(y, 1) == size(A, 1) && size(x, 2) == size(y, 2) "Dimension mismatch" + g = zero(T) + @inbounds for j = 1:size(A, 2) + for k in SparseArrays.nzrange(A, j) + i = A.rowval[k] + for m = 1:size(y, 2) + g += conj(y[i, m]) * A.nzval[k] * x[j, m] + end + end + end + return g +end + function Yao.AD.backward_params!(st::Tuple{<:MajoranaReg,<:MajoranaReg}, block::TimeEvolution, collector) out, outδ = st diff --git a/src/utils.jl b/src/utils.jl index fa27e81..003adc6 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -252,15 +252,16 @@ Convert a hamiltonian written as a YaoBlock into the corresponding """ function yaoham2majoranasquares(::Type{T}, yaoham::Add{2}) where {T<:Real} ham = zeros(T, 2nqubits(yaoham), 2nqubits(yaoham)) - for k in yaoham + @inbounds for k in yaoham if k isa Scale - ham += k.alpha * yaoham2majoranasquares(T, k.content) + fast_add!(ham, rmul!(yaoham2majoranasquares(T, k.content), k.alpha)) + #ham += k.alpha * yaoham2majoranasquares(T, k.content) elseif k isa KronBlock i1, i2 = kron2majoranaindices(k) ham[i1,i2] += 2 ham[i2,i1] -= 2 elseif k isa Add - ham += yaoham2majoranasquares(T, k) + fast_add!(ham, yaoham2majoranasquares(T, k)) elseif k isa PutBlock{2,1,ZGate} i1, i2 = kron2majoranaindices(k) ham[i1,i2] += 2 @@ -277,6 +278,19 @@ function yaoham2majoranasquares(::Type{T}, yaoham::Add{2}) where {T<:Real} return ham end +function fast_add!(A::AbstractMatrix{T}, B::SparseMatrixCSC{T}) where T + @assert size(A, 1) == size(B, 1) && size(A, 2) == size(B, 2) "Dimension mismatch" + for j = 1:size(B, 2) + for k in SparseArrays.nzrange(B, j) + @inbounds A[B.rowval[k],j] += B.nzval[k] + end + end + return A +end +function fast_add!(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T + return A .+= B +end + function yaoham2majoranasquares(::Type{T}, yaoham::KronBlock{2}) where {T<:Real} ham = spzeros(T, 2nqubits(yaoham), 2nqubits(yaoham)) i1, i2 = kron2majoranaindices(yaoham) diff --git a/test/runtests.jl b/test/runtests.jl index 8d90dca..5744ca1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -42,6 +42,24 @@ function ising_hamiltonian(nq, J, h) return hamiltonian end +@testset "fast_overlap" begin + nq = 4 + x = randn(ComplexF64, 10, 10) + y = randn(ComplexF64, 10, 10) + A = FLOYao.sprand(ComplexF64, 10, 10, 0.3) + r = FLOYao.fast_overlap(y, A, x) + @test isapprox(r, tr(y' * A * x), atol=1e-7) + @test isapprox(r, y ⋅ (A * x), atol=1e-7) +end + +@testset "fast_add!" begin + A = randn(10, 10) + B = FLOYao.sprand(10, 10, 0.2) + C = copy(A) + @test FLOYao.fast_add!(C, B) ≈ A + B + @test FLOYao.fast_add!(C, Matrix(B)) ≈ A + B +end + @testset "MajoranaRegister" begin nq = 2 mreg = FLOYao.zero_state(nq) From 450a3f0226a676a72a2b015c02b57b71447a12c0 Mon Sep 17 00:00:00 2001 From: GiggleLiu Date: Sat, 7 Oct 2023 10:31:31 +0800 Subject: [PATCH 2/4] fix tests --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index 5744ca1..4eb3281 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -57,6 +57,7 @@ end B = FLOYao.sprand(10, 10, 0.2) C = copy(A) @test FLOYao.fast_add!(C, B) ≈ A + B + C = copy(A) @test FLOYao.fast_add!(C, Matrix(B)) ≈ A + B end From f54a022504832ccc60291c37563bd5289c37e32b Mon Sep 17 00:00:00 2001 From: GiggleLiu Date: Thu, 12 Oct 2023 20:42:50 +0800 Subject: [PATCH 3/4] reply comments --- src/auto_diff.jl | 17 +---------------- src/utils.jl | 47 ++++++++++++++++++++++++++++++++++------------- 2 files changed, 35 insertions(+), 29 deletions(-) diff --git a/src/auto_diff.jl b/src/auto_diff.jl index b40863b..7055f21 100644 --- a/src/auto_diff.jl +++ b/src/auto_diff.jl @@ -44,27 +44,12 @@ function Yao.AD.backward_params!(st::Tuple{<:MajoranaReg,<:MajoranaReg}, return nothing end -# compute tr(y' * A * x) -function fast_overlap(y::AbstractVecOrMat{T}, A::SparseMatrixCSC{T}, x::AbstractVecOrMat{T}) where T - @assert size(x, 1) == size(A, 2) && size(y, 1) == size(A, 1) && size(x, 2) == size(y, 2) "Dimension mismatch" - g = zero(T) - @inbounds for j = 1:size(A, 2) - for k in SparseArrays.nzrange(A, j) - i = A.rowval[k] - for m = 1:size(y, 2) - g += conj(y[i, m]) * A.nzval[k] * x[j, m] - end - end - end - return g -end - function Yao.AD.backward_params!(st::Tuple{<:MajoranaReg,<:MajoranaReg}, block::TimeEvolution, collector) out, outδ = st ham = block.H majoranaham = yaoham2majoranasquares(ham) - g = outδ.state ⋅ (majoranaham * out.state) / 2 + g = fast_overlap(outδ.state, majoranaham, out.state) / 2 pushfirst!(collector, g) return nothing end diff --git a/src/utils.jl b/src/utils.jl index 003adc6..be3fbce 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -278,19 +278,6 @@ function yaoham2majoranasquares(::Type{T}, yaoham::Add{2}) where {T<:Real} return ham end -function fast_add!(A::AbstractMatrix{T}, B::SparseMatrixCSC{T}) where T - @assert size(A, 1) == size(B, 1) && size(A, 2) == size(B, 2) "Dimension mismatch" - for j = 1:size(B, 2) - for k in SparseArrays.nzrange(B, j) - @inbounds A[B.rowval[k],j] += B.nzval[k] - end - end - return A -end -function fast_add!(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T - return A .+= B -end - function yaoham2majoranasquares(::Type{T}, yaoham::KronBlock{2}) where {T<:Real} ham = spzeros(T, 2nqubits(yaoham), 2nqubits(yaoham)) i1, i2 = kron2majoranaindices(yaoham) @@ -362,3 +349,37 @@ function random_orthogonal_matrix(::Type{T}, n) where {T} end random_orthogonal_matrix(n) = random_orthogonal_matrix(Float64, n) + +# Fast sparse matrix operations +# # compute A .+= B +function fast_add!(A::AbstractMatrix, B::SparseMatrixCSC) + @assert size(A, 1) == size(B, 1) && size(A, 2) == size(B, 2) "Dimension mismatch" + for j = 1:size(B, 2) + for k in SparseArrays.nzrange(B, j) + @inbounds A[B.rowval[k],j] += B.nzval[k] + end + end + return A +end +function fast_add!(A::AbstractMatrix, B::AbstractMatrix) + return A .+= B +end + +# compute tr(y' * A * x) +function fast_overlap(y::AbstractVecOrMat{T1}, A::SparseMatrixCSC{T2}, x::AbstractVecOrMat{T3}) where {T1,T2,T3} + @assert size(x, 1) == size(A, 2) && size(y, 1) == size(A, 1) && size(x, 2) == size(y, 2) "Dimension mismatch" + g = zero(promote_type(T1, T2, T3)) + @inbounds for j = 1:size(A, 2) + for k in SparseArrays.nzrange(A, j) + i = A.rowval[k] + for m = 1:size(y, 2) + g += conj(y[i, m]) * A.nzval[k] * x[j, m] + end + end + end + return g +end + +function fast_overlap(y::AbstractVecOrMat{T1}, A::AbstractMatrix{T2}, x::AbstractVecOrMat{T3}) where {T1,T2,T3} + return y ⋅ (A * x) +end From 72c94a4f4bb46d78ccc666cb4a096294cceab277 Mon Sep 17 00:00:00 2001 From: janlukas Date: Thu, 12 Oct 2023 17:11:21 +0100 Subject: [PATCH 4/4] Bumped version number and added docstrings to fast_add and fast_overlap --- Project.toml | 2 +- src/auto_diff.jl | 5 ++--- src/utils.jl | 18 +++++++++++++++--- 3 files changed, 18 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index 7a1775c..22ba715 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "FLOYao" uuid = "6d9310a3-f1d0-41b7-8edb-11c1cf57cd2d" authors = ["janlukas.bosse and contributors"] -version = "1.4.0" +version = "1.4.1" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/auto_diff.jl b/src/auto_diff.jl index 7055f21..d734d5e 100644 --- a/src/auto_diff.jl +++ b/src/auto_diff.jl @@ -38,9 +38,8 @@ function Yao.AD.backward_params!(st::Tuple{<:MajoranaReg,<:MajoranaReg}, out, outδ = st ham = Yao.AD.generator(block) majoranaham = yaoham2majoranasquares(ham) - # multiply majoranaham with out.state and inner product with outδ.state - # i.e. g = outδ.state ⋅ (majoranaham * out.state) / 4 - pushfirst!(collector, fast_overlap(outδ.state, majoranaham, out.state)/4) + g = fast_overlap(outδ.state, majoranaham, out.state) / 4 + pushfirst!(collector, g) return nothing end diff --git a/src/utils.jl b/src/utils.jl index be3fbce..a769daa 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -350,8 +350,15 @@ end random_orthogonal_matrix(n) = random_orthogonal_matrix(Float64, n) -# Fast sparse matrix operations -# # compute A .+= B + +# ------------------------------------------- +# Utilities for fast sparse matrix operations +# ------------------------------------------- +""" + fast_add!(A::AbstractMatrix, B::SparseMatrixCSC) + +Fast implementation of `A .+= B` for sparse `B`. +""" function fast_add!(A::AbstractMatrix, B::SparseMatrixCSC) @assert size(A, 1) == size(B, 1) && size(A, 2) == size(B, 2) "Dimension mismatch" for j = 1:size(B, 2) @@ -361,11 +368,16 @@ function fast_add!(A::AbstractMatrix, B::SparseMatrixCSC) end return A end + function fast_add!(A::AbstractMatrix, B::AbstractMatrix) return A .+= B end -# compute tr(y' * A * x) +""" + fast_overlap(y::AbstractVecOrMat, A::SparseMatrixCSC, x::AbstractVecOrMat) + +Fast implementation of `tr(y' * A * x)` for sparse `A`. +""" function fast_overlap(y::AbstractVecOrMat{T1}, A::SparseMatrixCSC{T2}, x::AbstractVecOrMat{T3}) where {T1,T2,T3} @assert size(x, 1) == size(A, 2) && size(y, 1) == size(A, 1) && size(x, 2) == size(y, 2) "Dimension mismatch" g = zero(promote_type(T1, T2, T3))