upgrade (#8)

* upgrade

* new QSVD

* polish QSVD

Author: GiggleLiu

examples/QAOA.jl

@@ -4,7 +4,7 @@ using NLopt
 include("maxcut_gw.jl")
 
 HB(nbit::Int) = sum([put(nbit, i=>X) for i=1:nbit])
-tb = TimeEvolution(HB(3), 0.1)
+tb = TimeEvolution(HB(3) |> cache, 0.1; check_hermicity=false)
 function HC(W::AbstractMatrix)
     nbit = size(W, 1)
     ab = Any[]
@@ -22,11 +22,11 @@ function qaoa_circuit(W::AbstractMatrix, depth::Int; use_cache::Bool=false)
     hc = HC(W)
 	use_cache && (hb = hb |> cache; hc = hc |> cache)
     c = chain(nbit, [repeat(nbit, H, 1:nbit)])
-    append!(c, [chain(nbit, [time_evolve(hc, 0.0, tol=1e-5), time_evolve(hb, 0.0, tol=1e-5)]) for i=1:depth])
+    append!(c, [chain(nbit, [time_evolve(hc, 0.0, tol=1e-5, check_hermicity=false), time_evolve(hb, 0.0, tol=1e-5, check_hermicity=false)]) for i=1:depth])
 end
 
 
-function cobyla_optimize(circuit::MatrixBlock{N}, hc::MatrixBlock; niter::Int) where N
+function cobyla_optimize(circuit::AbstractBlock{N}, hc::AbstractBlock; niter::Int) where N
     function f(params, grad)
         reg = zero_state(N) |> dispatch!(circuit, params)
         loss = expect(hc, reg) |> real
@@ -60,7 +60,7 @@ using Random, Test
 	@test -expect(hc, product_state(5, bmask(sets[1]...)))+sum(W)/4 ≈ exact_cost
 
 	# build a QAOA circuit and start training.
-	qc = qaoa_circuit(W, 20; use_cache=false)
+	qc = qaoa_circuit(W, 20; use_cache=true)
 	opt_pl = nothing
 	opt_cost = Inf
 	for i = 1:10

examples/VQE.jl

@@ -2,7 +2,7 @@ using Yao
 using QuAlgorithmZoo
 using KrylovKit
 
-function ed_groundstate(h::MatrixBlock)
+function ed_groundstate(h::AbstractBlock)
     E, V = eigsolve(h |> mat, 1, :SR, ishermitian=true)
     println("Ground State Energy is $(E[1])")
     ArrayReg(V[1])

src/CircuitBuild.jl

@@ -115,7 +115,7 @@ function rand_single_gate()
 end
 
 """
-    rand_gate(nbit::Int, mbit::Int, [ngate::Int]) -> MatrixBlock
+    rand_gate(nbit::Int, mbit::Int, [ngate::Int]) -> AbstractBlock
 
 random nbit gate.
 """

src/Diff.jl

@@ -3,39 +3,40 @@ import Yao: expect, content, chcontent, mat, apply!
 using StatsBase
 
 ############# General Rotor ############
-const Rotor{N, T} = Union{RotationGate{N, T}, PutBlock{N, <:Any, <:RotationGate, <:Complex{T}}}
+const Rotor{N, T} = Union{RotationGate{N, T}, PutBlock{N, <:Any, <:RotationGate{<:Any, T}}}
 const CphaseGate{N, T} = ControlBlock{N,<:ShiftGate{T},<:Any}
 const DiffBlock{N, T} = Union{Rotor{N, T}, CphaseGate{N, T}}
 """
-    generator(rot::Rotor) -> MatrixBlock
+    generator(rot::Rotor) -> AbstractBlock
 
 Return the generator of rotation block.
 """
 generator(rot::RotationGate) = rot.block
 generator(rot::PutBlock{N, C, GT}) where {N, C, GT<:RotationGate} = PutBlock{N}(generator(rot|>content), rot |> occupied_locs)
 generator(c::CphaseGate{N}) where N = ControlBlock(N, c.ctrol_locs, ctrl_config, control(2,1,2=>Z), c.locs)
 
-abstract type AbstractDiff{GT, N, T} <: TagBlock{GT, N, T} end
+abstract type AbstractDiff{GT, N, T} <: TagBlock{GT, N} end
 Base.adjoint(df::AbstractDiff) = Daggered(df)
 
 istraitkeeper(::AbstractDiff) = Val(true)
 
 #################### The Basic Diff #################
 """
-    QDiff{GT, N, T} <: AbstractDiff{GT, N, Complex{T}}
+    QDiff{GT, N} <: AbstractDiff{GT, N, T}
     QDiff(block) -> QDiff
 
 Mark a block as quantum differentiable.
 """
-mutable struct QDiff{GT, N, T} <: AbstractDiff{GT, N, Complex{T}}
+mutable struct QDiff{GT, N, T} <: AbstractDiff{GT, N, T}
     block::GT
     grad::T
     QDiff(block::DiffBlock{N, T}) where {N, T} = new{typeof(block), N, T}(block, T(0))
 end
 content(cb::QDiff) = cb.block
 chcontent(cb::QDiff, blk::DiffBlock) = QDiff(blk)
 
-@forward QDiff.block mat, apply!
+@forward QDiff.block apply!
+mat(::Type{T}, df::QDiff) where T = mat(T, df.block)
 Base.adjoint(df::QDiff) = QDiff(content(df)')
 
 function YaoBlocks.print_annotation(io::IO, df::QDiff)
@@ -52,19 +53,19 @@ Mark a block as differentiable, here `GT`, `PT` is gate type, parameter type.
 Warning:
     please don't use the `adjoint` after `BPDiff`! `adjoint` is reserved for special purpose! (back propagation)
 """
-mutable struct BPDiff{GT, N, T, PT} <: AbstractDiff{GT, N, T}
+mutable struct BPDiff{GT, N, T} <: AbstractDiff{GT, N, T}
     block::GT
-    grad::PT
+    grad::T
     input::AbstractRegister
-    BPDiff(block::MatrixBlock{N, T}, grad::PT) where {N, T, PT} = new{typeof(block), N, T, typeof(grad)}(block, grad)
+    BPDiff(block::AbstractBlock{N}, grad::T) where {N, T} = new{typeof(block), N, T}(block, grad)
 end
-BPDiff(block::MatrixBlock) = BPDiff(block, zeros(parameters_eltype(block), nparameters(block)))
+BPDiff(block::AbstractBlock) = BPDiff(block, zeros(parameters_eltype(block), nparameters(block)))
 BPDiff(block::DiffBlock{N, T}) where {N, T} = BPDiff(block, T(0))
 
 content(cb::BPDiff) = cb.block
-chcontent(cb::BPDiff, blk::MatrixBlock) = BPDiff(blk)
+chcontent(cb::BPDiff, blk::AbstractBlock) = BPDiff(blk)
 
-@forward BPDiff.block mat
+mat(::Type{T}, df::BPDiff) where T = mat(T, df.block)
 function apply!(reg::AbstractRegister, df::BPDiff)
     if isdefined(df, :input)
         copyto!(df.input, reg)
@@ -105,7 +106,7 @@ autodiff(mode::Symbol, block::AbstractBlock) = autodiff(Val(mode), block)
 autodiff(::Val{:BP}, block::DiffBlock) = BPDiff(block)
 autodiff(::Val{:BP}, block::AbstractBlock) = block
 # Sequential, Roller and ChainBlock can propagate.
-function autodiff(mode::Val{:BP}, blk::Union{ChainBlock, Roller, Sequential})
+function autodiff(mode::Val{:BP}, blk::Union{ChainBlock, Sequential})
     chsubblocks(blk, autodiff.(mode, subblocks(blk)))
 end
 
@@ -148,11 +149,11 @@ Numeric differentiation.
 end
 
 """
-    opdiff(psifunc, diffblock::AbstractDiff, op::MatrixBlock)
+    opdiff(psifunc, diffblock::AbstractDiff, op::AbstractBlock)
 
 Operator differentiation.
 """
-@inline function opdiff(psifunc, diffblock::AbstractDiff, op::MatrixBlock)
+@inline function opdiff(psifunc, diffblock::AbstractDiff, op::AbstractBlock)
     r1, r2 = _perturb(()->expect(op, psifunc()) |> real, diffblock, π/2)
     diffblock.grad = (r2 - r1)/2
 end
@@ -215,14 +216,14 @@ end
 end
 
 """
-    backward!(δ::AbstractRegister, circuit::MatrixBlock) -> AbstractRegister
+    backward!(δ::AbstractRegister, circuit::AbstractBlock) -> AbstractRegister
 
 back propagate and calculate the gradient ∂f/∂θ = 2*Re(∂f/∂ψ*⋅∂ψ*/∂θ), given ∂f/∂ψ*.
 
 Note:
 Here, the input circuit should be a matrix block, otherwise the back propagate may not apply (like Measure operations).
 """
-backward!(δ::AbstractRegister, circuit::MatrixBlock) = apply!(δ, circuit')
+backward!(δ::AbstractRegister, circuit::AbstractBlock) = apply!(δ, circuit')
 
 """
     gradient(circuit::AbstractBlock, mode::Symbol=:ANY) -> Vector

src/Grover.jl

@@ -38,21 +38,21 @@ num_grover_step(psi::ArrayReg, oracle) = _num_grover_step(prob_match_oracle(psi,
 _num_grover_step(prob::Real) = Int(round(pi/4/sqrt(prob)))-1
 
 """
-    GroverIter{N, T}
+    GroverIter{N}
 
-    GroverIter(oracle, ref::ReflectBlock{N, T}, psi::ArrayReg, niter::Int)
+    GroverIter(oracle, ref::ReflectBlock{N}, psi::ArrayReg, niter::Int)
 
 an iterator that perform Grover operations step by step.
 An Grover operation consists of applying oracle and Reflection.
 """
-struct GroverIter{N, T}
+struct GroverIter{N}
     psi::ArrayReg
     oracle
-    ref::ReflectBlock{N, T}
+    ref::ReflectBlock{N}
     niter::Int
 end
 
-groveriter(psi::ArrayReg, oracle, ref::ReflectBlock{N, T}, niter::Int) where {N, T} = GroverIter{N, T}(psi, oracle, ref, niter)
+groveriter(psi::ArrayReg, oracle, ref::ReflectBlock{N}, niter::Int) where {N} = GroverIter{N}(psi, oracle, ref, niter)
 groveriter(psi::ArrayReg, oracle, niter::Int) = groveriter(psi, oracle, ReflectBlock(psi |> copy), niter)
 groveriter(psi::ArrayReg, oracle) = groveriter(psi, oracle, ReflectBlock(psi |> copy), num_grover_step(psi, oracle))
 
@@ -68,17 +68,17 @@ end
 Base.length(it::GroverIter) = it.niter
 
 """
-    groverblock(oracle, ref::ReflectBlock{N, T}, niter::Int=-1)
+    groverblock(oracle, ref::ReflectBlock{N}, niter::Int=-1)
     groverblock(oracle, psi::ArrayReg, niter::Int=-1)
 
 Return a ChainBlock/Sequential as Grover Iteration, the default `niter` will stop at the first optimal step.
 """
-function groverblock(oracle::MatrixBlock{N, T}, ref::ReflectBlock{N, T}, niter::Int=-1) where {N, T}
+function groverblock(oracle::AbstractBlock{N}, ref::ReflectBlock{N}, niter::Int=-1) where {N}
     if niter == -1 niter = num_grover_step(ref.psi, oracle) end
     chain(N, chain(oracle, ref) for i = 1:niter)
 end
 
-function groverblock(oracle, ref::ReflectBlock{N, T}, niter::Int=-1) where {N, T}
+function groverblock(oracle, ref::ReflectBlock{N}, niter::Int=-1) where {N}
     if niter == -1 niter = num_grover_step(ref.psi, oracle) end
     sequence(sequence(oracle, ref) for i = 1:niter)
 end

src/HHL.jl

@@ -1,15 +1,15 @@
 export hhlcircuit, hhlproject!, hhlsolve, HHLCRot
 
 """
-    HHLCRot{N, NC, T} <: PrimitiveBlock{N, Complex{T}}
+    HHLCRot{N, NC} <: PrimitiveBlock{N}
 
 Controlled rotation gate used in HHL algorithm, applied on N qubits.
 
     * cbits: control bits.
     * ibit:: the ancilla bit.
     * C_value:: the value of constant "C", should be smaller than the spectrum "gap".
 """
-struct HHLCRot{N, NC, T} <: PrimitiveBlock{N, Complex{T}}
+struct HHLCRot{N, NC, T} <: PrimitiveBlock{N}
     cbits::Vector{Int}
     ibit::Int
     C_value::T

src/HadamardTest.jl

@@ -3,15 +3,15 @@ export hadamard_test, hadamard_test_circuit, swap_test_circuit
 """
 see WiKi.
 """
-function hadamard_test_circuit(U::MatrixBlock{N}, ϕ::Real) where N
+function hadamard_test_circuit(U::AbstractBlock{N}, ϕ::Real) where N
     chain(N+1, put(N+1, 1=>H),
         put(N+1, 1=>Rz(ϕ)),
         control(N+1, 1, 2:N+1=>U),  # get matrix first, very inefficient
         put(N+1, 1=>H)
         )
 end
 
-function hadamard_test(U::MatrixBlock{N}, reg::AbstractRegister, ϕ::Real) where N
+function hadamard_test(U::AbstractBlock{N}, reg::AbstractRegister, ϕ::Real) where N
     c = hadamard_test_circuit(U, ϕ::Real)
     reg = join(reg, zero_state(1))
     expect(put(N+1, 1=>Z), reg |> c)

src/Miscellaneous.jl

@@ -16,13 +16,12 @@ function inverselines(nbit::Int; n_reg::Int=nbit)
     c
 end
 
-function singlet_block(::Type{T}, nbit::Int, i::Int, j::Int) where T
+function singlet_block(nbit::Int, i::Int, j::Int)
     unit = chain(nbit)
-    push!(unit, put(nbit, i=>chain(XGate{T}(), HGate{T}())))
-    push!(unit, control(nbit, -i, j=>XGate{T}()))
+    push!(unit, put(nbit, i=>chain(X, H)))
+    push!(unit, control(nbit, -i, j=>X))
 end
 
-singlet_block(nbit::Int, i::Int, j::Int) = singlet_block(ComplexF64, nbit, i, j)
 singlet_block() = singlet_block(2,1,2)
 
 Yao.mat(ρ::DensityMatrix{1}) = dropdims(state(ρ), dims=3)

src/QFT.jl

@@ -8,11 +8,11 @@ CRk(i::Int, j::Int, k::Int) = control([i, ], j=>shift(2π/(1<<k)))
 CRot(n::Int, i::Int) = chain(n, i==j ? kron(i=>H) : CRk(j, i, j-i+1) for j = i:n)
 QFTCircuit(n::Int) = chain(n, CRot(n, i) for i = 1:n)
 
-struct QFTBlock{N} <: PrimitiveBlock{N,ComplexF64} end
-mat(q::QFTBlock{N}) where N = applymatrix(q)
+struct QFTBlock{N} <: PrimitiveBlock{N} end
+mat(::Type{T}, q::QFTBlock{N}) where {T, N} = T.(applymatrix(q))
 
 apply!(reg::DefaultRegister{B}, ::QFTBlock) where B = (reg.state = ifft!(invorder_firstdim(reg |> state), 1)*sqrt(1<<nactive(reg)); reg)
-apply!(reg::DefaultRegister{B}, ::Daggered{N, T, <:QFTBlock}) where {B,N,T} = (reg.state = invorder_firstdim(fft!(reg|>state, 1)/sqrt(1<<nactive(reg))); reg)
+apply!(reg::DefaultRegister{B}, ::Daggered{N, <:QFTBlock}) where {B,N} = (reg.state = invorder_firstdim(fft!(reg|>state, 1)/sqrt(1<<nactive(reg))); reg)
 
 # traits
 ishermitian(q::QFTBlock{N}) where N = N==1
@@ -26,7 +26,7 @@ function print_block(io::IO, pb::QFTBlock{N}) where N
     printstyled(io, "QFT(1-$N)"; bold=true, color=:blue)
 end
 
-function print_block(io::IO, pb::Daggered{N, T,<:QFTBlock}) where {N, T}
+function print_block(io::IO, pb::Daggered{N,<:QFTBlock}) where {N, T}
     printstyled(io, "IQFT(1-$N)"; bold=true, color=:blue)
 end
 

src/QSVD.jl

@@ -0,0 +1,66 @@
+# Quantum SVD
+# Reference: https://arxiv.org/abs/1905.01353
+export QuantumSVD, circuit_qsvd, train_qsvd!, readout_qsvd
+
+"""
+Quantum singular value decomposition algorithm.
+
+    * `reg`, input register (A, B) as the target matrix to decompose,
+    * `circuit_a`, U matrix applied on register A,
+    * `circuit_b`, V matrix applied on register B,
+    * `optimizer`, the optimizer, normally we use `Adam(lr=0.1)`,
+    * `Nc`, log2 number of singular values kept,
+    * `maxiter`, the maximum number of iterations.
+"""
+function train_qsvd!(reg, circuit_a::AbstractBlock{Na}, circuit_b::AbstractBlock{Nb}, optimizer; Nc::Int=min(Na, Nb), maxiter::Int=100) where {Na, Nb}
+    nbit = Na+Nb
+    c = circuit_qsvd(circuit_a, circuit_b, Nc) |> autodiff(:QC)   # construct a differentiable circuit for training
+
+    obs = -mapreduce(i->put(nbit, i=>Z), +, (1:Na..., Na+Nc+1:Na+Nb...))
+    params = parameters(c)
+    for i = 1:maxiter
+        grad = opdiff.(() -> copy(reg) |> c, collect_blocks(AbstractDiff, c), Ref(obs))
+        QuAlgorithmZoo.update!(params, grad, optimizer)
+        println("Iter $i, Loss = $(Na+expect(obs, copy(reg) |> c))")
+        dispatch!(c, params)
+    end
+end
+
+"""build the circuit for quantum SVD training."""
+function circuit_qsvd(circuit_a::AbstractBlock{Na}, circuit_b::AbstractBlock{Nb}, Nc::Int) where {Na, Nb}
+    nbit = Na+Nb
+    cnots = chain(control(nbit, i+Na, i=>X) for i=1:Nc)
+    c = chain(concentrate(nbit, circuit_a, 1:Na), concentrate(nbit, circuit_b, Na+1:nbit), cnots)
+end
+
+"""read QSVD results"""
+function readout_qsvd(reg::AbstractRegister, circuit_a::AbstractBlock{Na}, circuit_b::AbstractBlock{Nb}, Nc::Int) where {Na, Nb}
+    reg = copy(reg) |> concentrate(Na+Nb, circuit_a, 1:Na) |> concentrate(Na+Nb, circuit_b, Na+1:Na+Nb)
+    _S = [select(reg, b|b<<Na).state[] for b in basis(Nc)]
+    S = abs.(_S)
+    order = sortperm(S, rev=true)
+    S, _S = S[order], _S[order]
+    mat(circuit_a)[order,:]'.*transpose(_S./S), S, transpose(mat(circuit_b)[order,:])
+end
+
+"""
+    QuantumSVD(M; kwargs...)
+Quantum SVD.
+    * `M`, the matrix to decompose, size should be (2^Na × 2^Nb), the sum of squared spectrum must be 1.
+kwargs includes
+    * `Nc`, log2 number of singular values kept,
+    * `circuit_a` and `circuit_b`, the circuit ansatz for `U` and `V` matrices,
+    * `maxiter`, maximum number of iterations,
+    * `optimizer`, default is `Adam(lr=0.1)`.
+"""
+function QuantumSVD(M::AbstractMatrix; Nc::Int=log2i(min(size(M)...)),
+        circuit_a=random_diff_circuit(log2i(size(M, 1)), 5, pair_ring(log2i(size(M, 1)))),
+        circuit_b=random_diff_circuit(log2i(size(M, 2)), 5, pair_ring(log2i(size(M, 2)))),
+        maxiter=200, optimizer=Adam(lr=0.1))
+
+    dispatch!(circuit_a, :random)
+    dispatch!(circuit_b, :random)
+    reg = ArrayReg(vec(M))
+    train_qsvd!(reg, circuit_a, circuit_b, optimizer, Nc=Nc, maxiter=maxiter)
+    readout_qsvd(reg, circuit_a, circuit_b, Nc)
+end