Generalize control gates to Qudits, introduce general spin site type

Author: mtfishman

ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl

@@ -5,7 +5,7 @@ using QuantumOperatorDefinitions:
   QuantumOperatorDefinitions, OpName, SiteType, StateName, has_fermion_string
 
 function QuantumOperatorDefinitions.SiteType(r::Index)
-  return SiteType(gettag(r, "sitetype", "Qudit"); length=Int(length(r)))
+  return SiteType(gettag(r, "sitetype", "Qudit"); dim=Int(length(r)))
 end
 
 function QuantumOperatorDefinitions.has_fermion_string(n::String, r::Index)

src/op.jl

@@ -1,16 +1,16 @@
+using LinearAlgebra: diag, qr
+
 struct OpName{Name,Params}
-  function OpName{Name,Params}(params::NamedTuple) where {Name,Params}
-    return new{Name,(; Params..., params...)}()
+  params::Params
+  function OpName{N}(params::NamedTuple) where {N}
+    return new{N,typeof(params)}(params)
   end
 end
 name(::OpName{Name}) where {Name} = Name
-params(::OpName{<:Any,Params}) where {Params} = Params
+params(n::OpName) = getfield(n, :params)
 
 Base.getproperty(n::OpName, name::Symbol) = getfield(params(n), name)
 
-OpName{Name,Params}(; kwargs...) where {Name,Params} = OpName{Name,Params}((; kwargs...))
-
-OpName{N}(params::NamedTuple) where {N} = OpName{N,params}()
 OpName{N}(; kwargs...) where {N} = OpName{N}((; kwargs...))
 
 # This compiles operator expressions, such as:
@@ -461,35 +461,84 @@ end
 alias(::OpName"√iSWAP") = √(OpName"iSWAP"())
 @op_alias "√iSwap" "√iSWAP"
 
-## TODO: Bring back these definitions.
-## function default_random_matrix(eltype::Type, s::Index...)
-##   n = prod(dim.(s))
-##   return randn(eltype, n, n)
-## end
-##
-## # Haar-random unitary
-## #
-## # Reference:
-## # Section 4.6
-## # http://math.mit.edu/~edelman/publications/random_matrix_theory.pdf
-## function op!(
-##   o::ITensor,
-##   ::OpName"RandomUnitary",
-##   ::SiteType"Generic",
-##   s1::Index,
-##   sn::Index...;
-##   eltype=ComplexF64,
-##   random_matrix=default_random_matrix(eltype, s1, sn...),
-## )
-##   s = (s1, sn...)
-##   Q, _ = NDTensors.qr_positive(random_matrix)
-##   t = itensor(Q, prime.(s)..., dag.(s)...)
-##   return settensor!(o, tensor(t))
-## end
-##
-## function op!(o::ITensor, ::OpName"randU", st::SiteType"Generic", s::Index...; kwargs...)
-##   return op!(o, OpName("RandomUnitary"), st, s...; kwargs...)
-## end
+# TODO: Generalize to more controlled operators, right now
+# there is only one controlled operator that is enabled when
+# all of the last values of each sites/Qudit is set. See:
+# https://docs.quantum.ibm.com/api/qiskit/qiskit.circuit.library.UCGate
+# https://arxiv.org/abs/2410.05122
+nsites(n::OpName"Controlled") = get(params(n), :ncontrol, 1) + nsites(n.arg)
+function (n::OpName"Controlled")(domain...)
+  # Number of target sites.
+  nt = nsites(n.arg)
+  # Number of control sites.
+  nc = get(params(n), :ncontrol, length(domain) - nt)
+  @assert length(domain) == nc + nt
+  d_control = prod(to_dim.(domain[1:nc]))
+  return cat(I(d_control), n.arg(domain[(nc + 1):end]...); dims=(1, 2))
+end
+@op_alias "CNOT" "Controlled" op = OpName"X"()
+@op_alias "CX" "Controlled" op = OpName"X"()
+@op_alias "CY" "Controlled" op = OpName"Y"()
+@op_alias "CZ" "Controlled" op = OpName"Z"()
+function alias(n::OpName"CPhase")
+  return controlled(OpName"Phase"(; params(n)...))
+end
+@op_alias "CPHASE" "CPhase"
+@op_alias "Cphase" "CPhase"
+function alias(n::OpName"CRx")
+  return controlled(OpName"Rx"(; params(n)...))
+end
+@op_alias "CRX" "CRx"
+function alias(::OpName"CRy")
+  return controlled(OpName"Ry"(; params(n)...))
+end
+@op_alias "CRY" "CRy"
+function alias(::OpName"CRz")
+  return controlled(OpName"Rz"(; params(n)...))
+end
+@op_alias "CRZ" "CRz"
+function alias(::OpName"CRn")
+  return controlled(; arg=OpName"Rn"(; params(n)...))
+end
+@op_alias "CRn̂" "CRn"
+
+@op_alias "CCNOT" "Controlled" ncontrol = 2 op = OpName"X"()
+@op_alias "Toffoli" "CCNOT"
+@op_alias "CCX" "CCNOT"
+@op_alias "TOFF" "CCNOT"
+
+@op_alias "CSWAP" "Controlled" ncontrol = 2 op = OpName"SWAP"()
+@op_alias "Fredkin" "CSWAP"
+@op_alias "CSwap" "CSWAP"
+@op_alias "CS" "CSWAP"
+
+@op_alias "CCCNOT" "Controlled" ncontrol = 3 op = OpName"X"()
+
+# Version of `sign` that returns one
+# if `x == 0`.
+function nonzero_sign(x)
+  iszero(x) && return one(x)
+  return sign(x)
+end
+
+function qr_positive(M::AbstractMatrix)
+  Q, R = qr(M)
+  Q′ = typeof(R)(Q)
+  signs = nonzero_sign.(diag(R))
+  Q′ = Q′ * Diagonal(signs)
+  R = Diagonal(conj.(signs)) * R
+  return Q′, R
+end
+
+# TODO: Store a random matrix or seed as a parameter
+# of the `OpName`?
+function (n::OpName"RandomUnitary")(domain...)
+  d = prod(to_dim.(domain))
+  elt = get(params(n), :eltype, Complex{Float64})
+  Q, _ = qr_positive(randn(elt, (d, d)))
+  return Q
+end
+@op_alias "randU" "RandomUnitary"
 
 # Expand the operator in a new basis.
 using LinearAlgebra: ⋅

src/sitetype.jl

@@ -1,11 +1,12 @@
 struct SiteType{T,Params}
   params::Params
+  function SiteType{N}(params::NamedTuple) where {N}
+    return new{N,typeof(params)}(params)
+  end
 end
 params(t::SiteType) = getfield(t, :params)
-
 Base.getproperty(t::SiteType, name::Symbol) = getfield(params(t), name)
 
-SiteType{N}(params) where {N} = SiteType{N,typeof(params)}(params)
 SiteType{N}(; kwargs...) where {N} = SiteType{N}((; kwargs...))
 
 SiteType(s::AbstractString; kwargs...) = SiteType{Symbol(s)}(; kwargs...)

src/sitetypes/qubit.jl

@@ -86,62 +86,15 @@ or:
 alias(n::OpName"R") = OpName"Rn"(; θ=n.θ, ϕ=n.ϕ, λ=-n.ϕ)
 =#
 # https://docs.quantum.ibm.com/api/qiskit/qiskit.circuit.library.UGate
-function Base.AbstractArray(n::OpName"Rn", ::Tuple{SiteType"Qubit"})
+# TODO: Generalize to `"Qudit"`, see:
+# https://quantumcomputing.stackexchange.com/questions/16251/how-does-a-general-rotation-r-hatn-theta-related-to-u-3-gate
+# https://quantumcomputing.stackexchange.com/questions/4249/decomposition-of-an-arbitrary-1-qubit-gate-into-a-specific-gateset
+# https://en.wikipedia.org/wiki/List_of_quantum_logic_gates#Other_named_gates
+# https://en.wikipedia.org/wiki/Spin_(physics)#Higher_spins
+function (n::OpName"Rn")(::SiteType"Qubit")
   return [
     cos(n.θ / 2) -exp(im * n.λ)*sin(n.θ / 2)
     exp(im * n.ϕ)*sin(n.θ / 2) exp(im * (n.ϕ + n.λ))*cos(n.θ / 2)
   ]
 end
 @op_alias "Rn̂" "Rn"
-
-# TODO: Generalize to `"Qudit"` and other SiteTypes.
-# https://docs.quantum.ibm.com/api/qiskit/qiskit.circuit.library.UCGate
-nsites(n::OpName"Controlled") = get(params(n), :ncontrol, 1) + nsites(n.arg)
-function Base.AbstractArray(n::OpName"Controlled", ts::Tuple{Vararg{SiteType"Qubit"}})
-  # Number of target qubits.
-  nt = nsites(n.arg)
-  # Number of control qubits.
-  nc = get(params(n), :ncontrol, length(ts) - nt)
-  @assert length(ts) == nc + nt
-  return [
-    I(2^nc) falses(2^nc, 2^nt)
-    falses(2^nt, 2^nc) AbstractArray(n.arg, ts[(nc + 1):end])
-  ]
-end
-@op_alias "CNOT" "Controlled" op = OpName"X"()
-@op_alias "CX" "Controlled" op = OpName"X"()
-@op_alias "CY" "Controlled" op = OpName"Y"()
-@op_alias "CZ" "Controlled" op = OpName"Z"()
-function alias(n::OpName"CPhase")
-  return controlled(OpName"Phase"(; params(n)...))
-end
-@op_alias "CPHASE" "CPhase"
-@op_alias "Cphase" "CPhase"
-function alias(n::OpName"CRx")
-  return controlled(OpName"Rx"(; params(n)...))
-end
-@op_alias "CRX" "CRx"
-function alias(::OpName"CRy")
-  return controlled(OpName"Ry"(; params(n)...))
-end
-@op_alias "CRY" "CRy"
-function alias(::OpName"CRz")
-  return controlled(OpName"Rz"(; params(n)...))
-end
-@op_alias "CRZ" "CRz"
-function alias(::OpName"CRn")
-  return controlled(; arg=OpName"Rn"(; params(n)...))
-end
-@op_alias "CRn̂" "CRn"
-
-@op_alias "CCNOT" "Controlled" ncontrol = 2 op = OpName"X"()
-@op_alias "Toffoli" "CCNOT"
-@op_alias "CCX" "CCNOT"
-@op_alias "TOFF" "CCNOT"
-
-@op_alias "CSWAP" "Controlled" ncontrol = 2 op = OpName"SWAP"()
-@op_alias "Fredkin" "CSWAP"
-@op_alias "CSwap" "CSWAP"
-@op_alias "CS" "CSWAP"
-
-@op_alias "CCCNOT" "Controlled" ncontrol = 3 op = OpName"X"()

src/sitetypes/qudit.jl

@@ -1,2 +1,3 @@
-Base.length(t::SiteType"Qudit") = t.length
-alias(::SiteType"Boson") = SiteType"Qudit"()
+Base.length(t::SiteType"Qudit") = t.dim
+alias(t::SiteType"Boson") = SiteType"Qudit"(; dim=t.dim)
+alias(t::SiteType"S") = SiteType"Qudit"(; dim=Int(2t.spin + 1))

src/state.jl

@@ -2,12 +2,13 @@ using Random: randstring
 
 struct StateName{Name,Params}
   params::Params
+  function StateName{N}(params::NamedTuple) where {N}
+    return new{N,typeof(params)}(params)
+  end
 end
 params(n::StateName) = getfield(n, :params)
-
 Base.getproperty(n::StateName, name::Symbol) = getfield(params(n), name)
 
-StateName{N}(params) where {N} = StateName{N,typeof(params)}(params)
 StateName{N}(; kwargs...) where {N} = StateName{N}((; kwargs...))
 
 StateName(s::AbstractString; kwargs...) = StateName{Symbol(s)}(; kwargs...)

test/test_basics.jl

@@ -1,5 +1,5 @@
 using QuantumOperatorDefinitions: OpName, SiteType, ⊗, expand, op, opexpr, state, nsites
-using LinearAlgebra: Diagonal
+using LinearAlgebra: Diagonal, I
 using Test: @test, @testset
 
 const real_elts = (Float32, Float64)
@@ -54,7 +54,7 @@ const elts = (real_elts..., complex_elts...)
         [0 1; 0 0],
       ),
       (
-        SiteType("Qudit"; length=3),
+        SiteType("Qudit"; dim=3),
         3,
         √2 * [0 1 0; 1 0 1; 0 1 0],
         √2 * [0 -im 0; im 0 -im; 0 im 0],
@@ -161,6 +161,19 @@ const elts = (real_elts..., complex_elts...)
       @test op(name, SiteType("Electron")) == Matrix(OpName(name), SiteType("Electron"))
     end
   end
+  @testset "Random unitary" begin
+    U = op("RandomUnitary")
+    @test eltype(U) === Complex{Float64}
+    @test U * U' ≈ I(2)
+
+    U = op("RandomUnitary"; eltype=Float32)
+    @test eltype(U) === Float32
+    @test U * U' ≈ I(2)
+
+    U = op("RandomUnitary", 3)
+    @test eltype(U) === Complex{Float64}
+    @test U * U' ≈ I(3)
+  end
   @testset "state" begin
     @test state(1) == [1, 0]
     @test state("0") == [1, 0]
@@ -176,4 +189,10 @@ const elts = (real_elts..., complex_elts...)
     @test state("|0⟩ + 2|+⟩") == state("0") + 2 * state("+")
     @test state("|0⟩ ⊗ |+⟩") == kron(state("0"), state("+"))
   end
+  @testset "Boson and spin" begin
+    for t in (SiteType("Qudit"; dim=3), SiteType("Boson"; dim=3), SiteType("S"; spin=1))
+      @test length(t) == 3
+      @test op("X", t) == op("X", 3)
+    end
+  end
 end