We use a rotation matrix U = (cs−sc) with c=cos(α), s=sin(α) and a real-valued angle α as an example. U has eigenvalues λ± = c±is = e±iα. Thus, ϕ=α/(2π). The corresponding eigenvectors are of the form u± = (1±i).
We use
6 t=
in the second register which allows us with probability 1−ϵ to get the correct phase up to t−⌈log(2+12ϵ)⌉ digits. Let us choose
1/4
epsilon <-## note the log in base-2
t-ceiling(log(2+1/(2*epsilon))/log(2))
digits <- digits
[1] 4
and therefore expect an error of less than
2^(-digits)
[1] 0.0625
We start with qubit 1 in state u+
S(1) * (H(1) * qstate(t+1, basis="")) x <-
and we define the gate corresponding to U
pi*3/7
alpha <- sin(alpha)
s <- cos(alpha)
c <-## note that R fills the matrix columns first
array(as.complex(c(c, -s, s, c)), dim=c(2,2))
M <- sqgate(bit=1, M=M, type=paste0("Uf")) Uf <-
Now we apply the Hadamard gate to qubits 2,,t+1
for(i in c(2:(t+1))) {
H(i) * x
x <- }
and the controlled Uf
for(i in c(2:(t+1))) {
cqgate(bits=c(i, 1),
x <-gate=sqgate(bit=1,
M=M, type=paste0("Uf", 2^(i-2)))) * x
M %*% M
M <-
}plot(x)
Next we apply the inverse Fourier transform
qft(x, inverse=TRUE, bits=c(2:(t+1)))
x <-plot(x)
x is now the state |˜φ⟩|u⟩. |˜φ⟩ is not necessarily a pure state. The next step is a projective measurement of |˜φ⟩
measure(x)
xtmp <- genStateNumber(which(xtmp$value==1)-1, t+1)
cbits <- sum(cbits[1:t]/2^(1:t))
phi <-
1:t] cbits[
[1] 0 0 1 1 1 0
phi
[1] 0.21875
Note that we can measure the complete state, because |u⟩ is not entangled to the rest. We find that usually
-alpha/(2*pi) phi
[1] 0.004464286
is indeed smaller than the maximal deviation 2−digits= 0.0625 we expect. The distribution of probabilities over the states in |˜φ⟩ is given as follows (factor 2 from dropping |u⟩)
plot(2*abs(x@coefs[seq(1,128,2)])^2, type="l",
ylab="p", xlab="state index")
The algorithm also works in case the specific eigenvector cannot be prepared. Starting with a random initial state |ψ⟩=∑ucu|u⟩, we may apply the very same algorithm and we will find the approximation to the phase φu with probability |cu|2(1−ϵ).
We prepare the second register in the state (11) = (1−i)u++(1+i)u−.
(H(1) * qstate(t+1, basis="")) x <-
This implies that we will find both φu with equal probability.
for(i in c(2:(t+1))) {
H(i) * x
x <-
} array(as.complex(c(c, -s, s, c)), dim=c(2,2))
M <-for(i in c(2:(t+1))) {
cqgate(bits=c(i, 1),
x <-gate=sqgate(bit=1,
M=M, type=paste0("Uf", 2^(i-2)))) * x
M %*% M
M <-
} qft(x, inverse=TRUE, bits=c(2:(t+1)))
x <-
function(x, t) {
measurephi <- measure(x)
xtmp <- genStateNumber(which(xtmp$value==1)-1, t+1)
cbits <- sum(cbits[1:t]/2^(1:t))
phi <-return(invisible(phi))
} measurephi(x, t=t)
phi <-2*pi*phi
[1] 1.374447
-c(+alpha, 2*pi-alpha)/2/pi phi
[1] 0.004464286 -0.566964286
We can draw the probability distribution again and observe the two peaks corresponding to the two eigenvalues
plot(abs(x@coefs)^2, type="l",
ylab="p", xlab="state index")
Let’s measure 1000 times, which is easily possible in our simulator
c()
phi <-for(i in c(1:N)) {
measurephi(x, t)
phi[i] <-
}hist(phi, breaks=2^t, xlim=c(0,1))
abline(v=c(alpha/2/pi, 1-alpha/2/pi), lwd=2, col="red")
The red vertical lines indicate the true values.