From 0fbb786f1dd6adf5a2769e9c058556466961499f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jofre=20Vall=C3=A8s=20Muns?= <61060572+jofrevalles@users.noreply.github.com> Date: Fri, 8 Sep 2023 20:24:45 +0200 Subject: [PATCH] Implement generic `SU{N}` operator (#37) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Add Random in dependencies * Implement SU(N) operator * Refactor SU * Refactor SU{N} gate * Add tests for SU{N} gate * Upgrade to ComplexF64 the default precision type of rand(SU{N}, ...) * Apply @mofeing suggestions from code review Co-authored-by: Sergio Sánchez Ramírez <15837247+mofeing@users.noreply.github.com> * Fix code * Fix error * change data precision * Update src/Gate.jl Co-authored-by: Sergio Sánchez Ramírez <15837247+mofeing@users.noreply.github.com> * Fix naming * Update tests accordingly * Update docs * Apply @mofeing suggestions from code review Co-authored-by: Sergio Sánchez Ramírez <15837247+mofeing@users.noreply.github.com> * Fix code from last code review * Fix code from last code review * Fix tests * Fix tests * Fix tests --------- Co-authored-by: Sergio Sánchez Ramírez <15837247+mofeing@users.noreply.github.com> --- docs/src/api/gates.md | 9 +++++ src/Array.jl | 19 ++++++++-- src/Gate.jl | 23 ++++++++++- test/Array_test.jl | 6 +++ test/Gate_test.jl | 88 ++++++++++++++++++++++++++++++++++++++++++- 5 files changed, 140 insertions(+), 5 deletions(-) diff --git a/docs/src/api/gates.md b/docs/src/api/gates.md index 4b9dfb7..999dd25 100644 --- a/docs/src/api/gates.md +++ b/docs/src/api/gates.md @@ -101,3 +101,12 @@ Control ```@docs Swap ``` + +## Special Unitary gate + +!!! warn "Experimental interface" + This interface is experimental and may change in the future. + +```@docs +Quac.SU{N} +``` \ No newline at end of file diff --git a/src/Array.jl b/src/Array.jl index 25d6416..75a9d89 100644 --- a/src/Array.jl +++ b/src/Array.jl @@ -1,6 +1,6 @@ import Base: Matrix, Array import LinearAlgebra: Diagonal, diag, eigvals, eigvecs, eigen -using LinearAlgebra: Eigen, LinearAlgebra +using LinearAlgebra: Eigen, LinearAlgebra, qr # preferred representation function arraytype end @@ -13,7 +13,7 @@ for G in [I, Z, S, Sd, T, Td, Rz] @eval arraytype(::Type{<:Gate{$G}}) = Diagonal end -for G in [X, Y, H, Rx, Ry] +for G in [X, Y, H, Rx, Ry, SU] @eval arraytype(::Type{<:Gate{$G}}) = Matrix end @@ -22,7 +22,7 @@ end Matrix(x::Gate) = Matrix{ComplexF32}(x) Matrix(::Type{T}) where {T<:Gate} = Matrix{ComplexF32}(T) -for Op in [I, X, Y, Z, H, S, Sd, T, Td, Swap] +for Op in [I, X, Y, Z, H, S, Sd, T, Td, Swap, SU] @eval Matrix{T}(::G) where {T,G<:Gate{$Op}} = Matrix{T}(G) end @@ -119,6 +119,19 @@ function Matrix{T}(g::Gate{<:Control}) where {T} return M end +function Base.rand(::Type{SU{N}}, lanes::NTuple{M, Int}; eltype::Type = ComplexF64) where {N, M} + # keep unitary matrix Q from QR decomposition + q, _ = qr(rand(eltype, N, N)) + + SU{N}(lanes...; array = Matrix(q)) +end + +Base.rand(::Type{Gate{SU{N}}}, lanes::Integer...; kwargs...) where {N} = rand(SU{N}, lanes; kwargs...) + +function Matrix{T}(g::Gate{<:SU{N}}) where {T, N} + return g.array |> Matrix{T} +end + Array(x::Gate) = Array{ComplexF32}(x) Array(::Type{T}) where {T<:Gate} = Array{ComplexF32}(T) diff --git a/src/Gate.jl b/src/Gate.jl index d7c9f92..221151d 100644 --- a/src/Gate.jl +++ b/src/Gate.jl @@ -4,7 +4,7 @@ export X, Y, Z, H, S, Sd, T, Td export isparametric, parameters export Rx, Ry, Rz, U1, U2, U3, Hz export Rxx, Ryy, Rzz -export Control, Swap, FSim +export Control, Swap, FSim, SU export CX, CY, CZ, CRx, CRy, CRz export control, target, operator export Pauli, Phase @@ -206,6 +206,17 @@ abstract type Control{Op<:Operator} <: Operator{NamedTuple{(:target,),Tuple{Oper Base.length(::Type{Control{T}}) where {T<:Operator} = 1 + length(T) parameters(::Type{Control{Op}}) where {Op} = parameters(Op) +""" + SU(N)(lane_1, lane_2, ..., lane_log2(N)) + +The `SU{N}` multi-qubit general unitary gate that can be used to represent any unitary matrix that acts on +`log2(N)` qubits. A new random `SU{N}` can be created with `rand(SU{N}, lanes...)`, where `N` is the dimension + of the unitary matrix and `lanes` are the qubit lanes on which the gate acts. +""" +abstract type SU{N} <: Operator{NamedTuple{(:array,), Tuple{Matrix}}} end +Base.length(::Type{SU{N}}) where {N} = + ispow2(N) ? log2(N) |> Int : throw(DomainError(N, "N must be a power of 2")) + for Op in [:X, :Y, :Z, :Rx, :Ry, :Rz] @eval const $(Symbol("C" * String(Op))) = Control{$Op} end @@ -248,6 +259,15 @@ for Op in [:I, :X, :Y, :Z, :H, :S, :Sd, :T, :Td, :U2, :U3, :Rx, :Ry, :Rz, :Rxx, @eval $Op(lanes...; params...) = Gate{$Op}(lanes...; params...) end +function SU{N}(lanes...; array::Matrix) where {N} + ispow2(N) || throw(DomainError(N, "N must be a power of 2")) + 2^length(lanes) == N || throw(ArgumentError("SU{$N} requires $(log2(N) |> Int) lanes")) + size(array) == (N,N) || throw(ArgumentError("`array` must be a (N,N)-size matrix")) + isapprox(array * adjoint(array), Matrix{ComplexF64}(LinearAlgebra.I, N, N)) || throw(ArgumentError("`array` is not unitary")) + + Gate{SU{N}}(lanes...; array) +end + Base.sqrt(g::Gate{X}) = Rx(lanes(g)..., θ = π / 2) Base.sqrt(g::Gate{Y}) = Ry(lanes(g)..., θ = π / 2) Base.sqrt(g::Gate{Z}) = S(lanes(g)...) @@ -283,6 +303,7 @@ Base.getproperty(g::Gate{Op}, i::Symbol) where {Op} = i ∈ propertynames(g) ? p Base.adjoint(::Type{<:Gate{Op}}) where {Op} = Gate{adjoint(Op)} Base.adjoint(::Type{Gate{Op,N,P}}) where {Op,N,P} = Gate{adjoint(Op),N,P} Base.adjoint(g::Gate{Op}) where {Op} = Gate{Op'}(lanes(g)...; [key => -val for (key, val) in pairs(parameters(g))]...) +Base.adjoint(g::Gate{SU{N}}) where {N} = Gate{SU{N}}(lanes(g)...; array = adjoint(g.array)) # NOTE useful type piracy Base.rand(::Type{NamedTuple{N,T}}) where {N,T} = NamedTuple{N}(rand(type) for type in T.parameters) diff --git a/test/Array_test.jl b/test/Array_test.jl index 509d1be..6bfa835 100644 --- a/test/Array_test.jl +++ b/test/Array_test.jl @@ -25,6 +25,12 @@ @test size(Matrix(Gate{Op})) == (2^length(Op), 2^length(Op)) end + # Special case for SU{N} + for N in [2, 4, 8] + @test Matrix(rand(Gate{SU{N}}, 1:Int(log2(N))...)) isa Matrix{ComplexF32} + @test size(Matrix(rand(Gate{SU{N}}, 1:Int(log2(N))...))) == (N, N) + end + for Op in [ I, X, diff --git a/test/Gate_test.jl b/test/Gate_test.jl index 841ba7e..89fe7a1 100644 --- a/test/Gate_test.jl +++ b/test/Gate_test.jl @@ -1,5 +1,6 @@ @testset "Gate" begin using Quac: I + using LinearAlgebra: qr @testset "Base.length" begin for Op in [ @@ -28,6 +29,8 @@ Control{Swap}, Control{Control{Swap}}, Control{Control{Control{Swap}}}, + SU{2}, + SU{4}, ] @test length(Gate{Op}) === length(Op) end @@ -37,6 +40,15 @@ for Op in [I, X, Y, Z, H, S, Sd, T, Td, Rx, Ry, Rz, Rxx, Ryy, Rzz, U2, U3] @test Gate{Op}(range(1, length = length(Op))...; rand(parameters(Op))...) !== nothing end + + # Special case for SU{N} + for N in [2, 4, 8] + _lanes = range(1, length = log2(N) |> Int) + rand_matrix = rand(ComplexF32, N, N) + q, _ = qr(rand_matrix) + array = Matrix{ComplexF32}(q) + @test Gate{SU{N}}(_lanes...; array = array) !== nothing + end end @testset "Constructor aliases" begin @@ -74,6 +86,15 @@ P = parameters(Op) @test Op(1:N...; rand(P)...) isa Gate{Op,N,P} end + + # Special case for SU{N} + for N in [2, 4, 8] + _lanes = range(1, length = log2(N) |> Int) + rand_matrix = rand(ComplexF64, N, N) + q, _ = qr(rand_matrix) + array = Matrix{ComplexF64}(q) + @test SU{N}(_lanes...; array = array) isa Gate{SU{N},log2(N) |> Int, NamedTuple{(:array,),Tuple{Matrix}}} + end end @testset "lanes" begin @@ -109,6 +130,15 @@ ] @test lanes(Gate{Op}(1:length(Op)...; rand(parameters(Op))...)) === tuple(1:length(Op)...) end + + # Special case for SU{N} + for N in [2, 4, 8] + _lanes = 1:length(SU{N}) |> collect + rand_matrix = rand(ComplexF32, N, N) + q, _ = qr(rand_matrix) + array = Matrix{ComplexF32}(q) + @test lanes(Gate{SU{N}}(_lanes...; array = array)) === tuple(1:length(SU{N})...) + end end @testset "operator" begin @@ -145,6 +175,15 @@ # NOTE `operator(Gate{Op})` is implied @test operator(Gate{Op}(1:length(Op)...; rand(parameters(Op))...)) === Op end + + # Special case for SU{N} + for N in [2, 4, 8] + _lanes = 1:length(SU{N}) |> collect + rand_matrix = rand(ComplexF32, N, N) + q, _ = qr(rand_matrix) + array = Matrix{ComplexF32}(q) + @test operator(Gate{SU{N}}(_lanes...; array = array)) === SU{N} + end end @testset "parameters" begin @@ -183,6 +222,17 @@ params = rand(parameters(Op)) @test parameters(Gate{Op}(1:length(Op)...; params...)) === params end + + # Special case for SU{N} + for N in [2, 4, 8] + @test parameters(Gate{SU{N}}) === parameters(SU{N}) + + _lanes = 1:length(SU{N}) |> collect + rand_matrix = rand(ComplexF32, N, N) + q, _ = qr(rand_matrix) + array = Matrix{ComplexF32}(q) + @test parameters(Gate{SU{N}}(_lanes...; array = array)).array === array + end end @testset "adjoint" begin @@ -230,6 +280,17 @@ @test adjoint(Gate{Op}(1:length(Op)...; params...)) === Gate{adjoint(Op)}(1:length(Op)...; [key => -val for (key, val) in pairs(params)]...) end + + # Special case for SU{N} + for N in [2, 4, 8] + @test_throws MethodError adjoint(Gate{SU{N}}) + + _lanes = 1:length(SU{N}) |> collect + rand_matrix = rand(ComplexF64, N, N) + q, _ = qr(rand_matrix) + + @test adjoint(SU{N}(_lanes...; array = Matrix{ComplexF64}(q))).array == adjoint(Matrix{ComplexF64}(q)) + end end @testset "Base.rand" begin @@ -262,6 +323,9 @@ Control{Swap}, Control{Control{Swap}}, Control{Control{Control{Swap}}}, + SU{2}, + SU{4}, + SU{8}, ] N = length(Op) P = parameters(Op) @@ -299,6 +363,9 @@ Control{Swap}, Control{Control{Swap}}, Control{Control{Control{Swap}}}, + SU{2}, + SU{4}, + SU{8}, ] @test keys(parameters(Op)) ⊆ propertynames(Gate{Op}) @@ -314,7 +381,7 @@ @testset "targettype" begin using Quac: targettype - for Op in [I, X, Y, Z, H, S, Sd, T, Td, Rx, Ry, Rz, Rxx, Ryy, Rzz, U2, U3, Swap, FSim] + for Op in [I, X, Y, Z, H, S, Sd, T, Td, Rx, Ry, Rz, Rxx, Ryy, Rzz, U2, U3, Swap, FSim, SU{2}, SU{4}, SU{8}] @test targettype(Gate{Op}) === Op end @@ -342,6 +409,25 @@ end end + @testset "random unitary" begin + # test_throws on a non-unitary matrix + @test_throws ArgumentError SU{4}(1, 2; array = rand(ComplexF32, 4, 4)) + + # test_throws on a non-square matrix + @test_throws ArgumentError SU{4}(1, 2; array = rand(ComplexF32, 4, 2)) + + # test_throws on a matrix without size (N, N) + @test_throws ArgumentError SU{4}(1, 2; array = rand(ComplexF32, 2, 2)) + + # test_throws SU{N} with N not a power of 2 + @test_throws DomainError SU{3}(1, 2; array = rand(ComplexF32, 3, 3)) + + # test_throws when there are not log2(N) lanes + rand_matrix = rand(ComplexF32, 4, 4) + q, _ = qr(rand_matrix) + @test_throws ArgumentError SU{4}(1, 2, 3; array = Matrix{ComplexF32}(q)) + end + @testset "target" begin for Op in [I, X, Y, Z, H, S, Sd, T, Td, Rx, Ry, Rz, Rxx, Ryy, Rzz, U2, U3, Swap, FSim] @test_throws MethodError target(rand(Gate{Op}, 1:length(Op)...))