Quantum Fourier Transformation and Phase Estimation
Quantum Fourier Transformation
using Yao
# Control-R(k) gate in block-A
A(i::Int, j::Int, k::Int) = control([i, ], j=>shift(2π/(1<<k)))
# block-B
B(n::Int, i::Int) = chain(i==j ? kron(i=>H) : A(j, i, j-i+1) for j = i:n)
QFT(n::Int) = chain(n, B(n, i) for i = 1:n)
# define QFT and IQFT block.
num_bit = 5
qft = QFT(num_bit)
iqft = qft' # get the hermitian conjugate
Total: 5, DataType: Complex{Float64}
chain
├─ chain
│ └─ kron
│ └─ 5=>H gate
├─ chain
│ ├─ control(5)
│ │ └─ (4,)=>Phase Shift Gate:-1.5707963267948966
│ └─ kron
│ └─ 4=>H gate
├─ chain
│ ├─ control(5)
│ │ └─ (3,)=>Phase Shift Gate:-0.7853981633974483
│ ├─ control(4)
│ │ └─ (3,)=>Phase Shift Gate:-1.5707963267948966
│ └─ kron
│ └─ 3=>H gate
├─ chain
│ ├─ control(5)
│ │ └─ (2,)=>Phase Shift Gate:-0.39269908169872414
│ ├─ control(4)
│ │ └─ (2,)=>Phase Shift Gate:-0.7853981633974483
│ ├─ control(3)
│ │ └─ (2,)=>Phase Shift Gate:-1.5707963267948966
│ └─ kron
│ └─ 2=>H gate
└─ chain
├─ control(5)
│ └─ (1,)=>Phase Shift Gate:-0.19634954084936207
├─ control(4)
│ └─ (1,)=>Phase Shift Gate:-0.39269908169872414
├─ control(3)
│ └─ (1,)=>Phase Shift Gate:-0.7853981633974483
├─ control(2)
│ └─ (1,)=>Phase Shift Gate:-1.5707963267948966
└─ kron
└─ 1=>H gate
The basic building block - controled phase shift gate is defined as
In Yao, factory methods for blocks will be loaded lazily. For example, if you missed the total number of qubits of chain
, then it will return a function that requires an input of an integer. So the following two statements are equivalent
control([4, ], 1=>shift(-2π/(1<<4)))(5) == control(5, [4, ], 1=>shift(-2π/(1<<4)))
true
Both of then will return a ControlBlock
instance. If you missed the total number of qubits. It is OK. Just go on, it will be filled when its possible.
Once you have construct a block, you can inspect its matrix using mat
function. Let's construct the circuit in dashed box A, and see the matrix of $R_4$ gate
julia> a = A(4, 1, 4)(5)
Total: 5, DataType: Complex{Float64}
control(4)
└─ 1=>Phase Shift Gate:-0.39269908169872414
julia> mat(a.block)
2×2 Diagonal{Complex{Float64}}:
1.0+0.0im ⋅
⋅ 0.92388-0.382683im
Similarly, you can use put
and chain
to construct PutBlock
(basic placement of a single gate) and ChainBlock
(sequential application of MatrixBlock
s) instances. Yao.jl
view every component in a circuit as an AbstractBlock
, these blocks can be integrated to perform higher level functionality.
You can check the result using classical fft
# if you're using lastest julia, you need to add the fft package.
using FFTW: fft, ifft
using LinearAlgebra: I
using Test
@test chain(num_bit, qft, iqft) |> mat ≈ I
# define a register and get its vector representation
reg = rand_state(num_bit)
rv = reg |> statevec |> copy
# test fft
reg_qft = apply!(copy(reg) |>invorder!, qft)
kv = ifft(rv)*sqrt(length(rv))
@test reg_qft |> statevec ≈ kv
# test ifft
reg_iqft = apply!(copy(reg), iqft)
kv = fft(rv)/sqrt(length(rv))
@test reg_iqft |> statevec ≈ kv |> invorder
QFT and IQFT are different from FFT and IFFT in three ways,
- they are different by a factor of $\sqrt{2^n}$ with $n$ the number of qubits.
- the little end and big end will exchange after applying QFT or IQFT.
- due to the convention, QFT is more related to IFFT rather than FFT.
Phase Estimation
Since we have QFT and IQFT blocks we can then use them to realize phase estimation circuit, what we want to realize is the following circuit
In the following simulation, we use equivalent QFTBlock
in the Yao.Zoo
module rather than the above chain block, it is faster than the above construction because it hides all the simulation details (yes, we are cheating :D) and get the equivalent output.
using Yao
using Yao.Blocks
using Yao.Intrinsics
function phase_estimation(reg1::DefaultRegister, reg2::DefaultRegister, U::GeneralMatrixGate{N}, nshot::Int=1) where {N}
M = nqubits(reg1)
iqft = QFT(M) |> adjoint
HGates = rollrepeat(M, H)
control_circuit = chain(M+N)
for i = 1:M
push!(control_circuit, control(M+N, (i,), (M+1:M+N...,)=>U))
if i != M
U = matrixgate(mat(U) * mat(U))
end
end
# calculation
# step1 apply hadamard gates.
apply!(reg1, HGates)
# join two registers
reg = join(reg1, reg2)
# using iqft to read out the phase
apply!(reg, sequence(control_circuit, focus(1:M...), iqft))
# measure the register (on focused bits), if the phase can be exactly represented by M qubits, only a single shot is needed.
res = measure(reg; nshot=nshot)
# inverse the bits in result due to the exchange of big and little ends, so that we can get the correct phase.
breflect.(M, res)./(1<<M), reg
end
phase_estimation (generic function with 2 methods)
Here, reg1
($Q_{1-5}$) is used as the output space to store phase ϕ, and reg2
($Q_{6-8}$) is the input state which corresponds to an eigenvector of oracle matrix U
. The algorithm detials can be found here.
In this function, HGates
corresponds to circuit block in dashed box A
, control_circuit
corresponds to block in dashed box B
. matrixgate
is a factory function for GeneralMatrixGate
.
Here, the only difficult concept is focus
, focus
returns a FunctionBlock
, that will make focused bits the active bits. An operator sees only active bits, and operating active space is more efficient, most importantly, it becomes much easier to integrate blocks. However, it has the potential ability to change line orders, for safety consideration, you may also need safer Concentrator
.
r = rand_state(6)
apply!(r, focus(4,1,2)) # or equivalently using focus!(r, [4,1,2])
nactive(r)
3
Then we will have a check to above function
using LinearAlgebra: qr, Diagonal
rand_unitary(N::Int) = qr(randn(N, N)).Q
M = 5
N = 3
# prepair oracle matrix U
V = rand_unitary(1<<N)
phases = rand(1<<N)
ϕ = Int(0b11101)/(1<<M)
phases[3] = ϕ # set the phase of the 3rd eigenstate manually.
signs = exp.(2pi*im.*phases)
U = V*Diagonal(signs)*V' # notice U is unitary
# the state with phase ϕ
psi = U[:,3]
res, reg = phase_estimation(zero_state(M), register(psi), GeneralMatrixGate(U))
println("Phase is 2π * $(res[]), the exact value is 2π * $ϕ")
Phase is 2π * 0.5, the exact value is 2π * 0.90625