Variational Quantum Eigensolver for Ising Hamiltonians

Jupyter notebook here.

Variational Quantum Eigensolver (VQE) algorithm is a quantum realization of Ritz’s variational method.

A more lengthy explanation is given by Michał Stęchły. Here I will just work out how it operates.

Ising hamiltonians are the easiest to map to quantum circuits and learning how to encode them will help us with more complicated hamiltonians.

using LinearAlgebra: I, tr

# Computational basis
𝟎 = [1,0]
𝟏 = [0,1]

X = [0 1; 1 0]
Z = [1 0; 0 -1] # we measure in Z basis

H = [1 1; 1 -1]/√2
p = H * 𝟎 # = |+⟩
m = H * 𝟏 # = |-⟩

Rotate qubits to measure other Pauli matrices (than $\sigma_x$)

As an example, if we measure in the Z basis its not possible to distinguish between $|+\rangle$ and $|-\rangle$.

@info "⟨+|Z|+⟩ ≈ ⟨-|Z|-⟩", p'*Z*p ≈ m'* Z*m

┌ Info: ("⟨+|Z|+⟩ ≈ ⟨-|Z|-⟩", true)

Rotating the basis in a certain way allows us to unequivocally distinguish them.

Ry(θ) = [cos(θ/2) -sin(θ/2); sin(θ/2) cos(θ/2)]

@assert Ry(π/2)*p ≈ 𝟏
@assert Ry(π/2)*m ≈ 𝟎

p_rot = Ry(-π/2)*p
m_rot = Ry(-π/2)*m
@info "⟨+|Z|Ry(-π/2)|+⟩ ≈ ⟨-|Z|Ry(-π/2)|-⟩", p_rot'*Z*p_rot ≈ m_rot'* Z*m_rot
@info "Outcome of ⟨+|Z|Ry(-π/2)|+⟩ = ", p_rot'*Z*p_rot
@info "Outcome of ⟨-|Z|Ry(-π/2)|-⟩ = ", m_rot'* Z*m_rot

┌ Info: ("⟨+|Z|Ry(-π/2)|+⟩ ≈ ⟨-|Z|Ry(-π/2)|-⟩", false)
┌ Info: ("Outcome of ⟨+|Z|Ry(-π/2)|+⟩ = ", 0.9999999999999998)
┌ Info: ("Outcome of ⟨-|Z|Ry(-π/2)|-⟩ = ", -0.9999999999999998)

The solution is to apply

• $R_y(−\pi/2)$ if Hamiltonian has $\hat{X}$
• $R_x(\pi/2)$ if Hamiltonian has $\hat{Y}$
• $I$ if Hamiltonian has $\hat{Z}$

before measurement. Consequently, the circuit will be different for each Pauli matrix.

Divide the Hamiltonian into single Pauli terms

𝓗₁ = 2Z
𝓗₂ = X
𝓗₃ = I(2)
𝓗 = 𝓗₁ + 𝓗₂ + 𝓗₃
@assert 𝓗 ≈ [3 1; 1 -1]

Initial state

We use a parametrized initial wavefunction for reasons that will be explained later but you can already guess that we will do a variational search.

# Ansatz will be Ry(θ)*𝟎
ansatz(θ) = Ry(θ)*𝟎

A circuit for each term

function vqe(θ, verbose=false)

   verbose && @info "θ = $(θ[1])"
   # Circuit for 𝓗₁: ψ ---[ h₁ )===

   ψ = ansatz(θ[1])

   h₁ = 2ψ' * Z * ψ
   verbose && @info "  E₁ = $h₁"

   # Circuit for 𝓗₂: ψ ---[ Ry ]---[ h₂ )===

   ψ = Ry(-π/2)*ansatz(θ[1])

   h₂ = ψ' * Z * ψ
   verbose && @info "  E₂ = $h₂"

   # Circuit for 𝓗₃:

   h₃ = 1
   verbose && @info "  E₃ = $h₃"

   verbose && @info "  ⟨ψ($θ[1])|𝓗|ψ($θ[1])⟩ = $(h₁ + h₂ + h₃)"
   h₁ + h₂ + h₃
end

For two single parameter values we have

vqe([0],true)
vqe([π],true)

┌ Info: θ = 0
┌ Info:   E₁ = 2.0
┌ Info:   E₂ = 2.220446049250313e-16
┌ Info:   E₃ = 1
┌ Info:   ⟨ψ([0][1])|𝓗|ψ([0][1])⟩ = 3.0
┌ Info: θ = π
┌ Info:   E₁ = -2.0
┌ Info:   E₂ = -2.220446049250313e-16
┌ Info:   E₃ = 1
┌ Info:   ⟨ψ([π][1])|𝓗|ψ([π][1])⟩ = -1.0

using Optim

result = optimize(vqe, [0.0], BFGS())

 * Status: success

 * Candidate solution
    Final objective value:     -1.236068e+00

 * Found with
    Algorithm:     BFGS

 * Convergence measures
    |x - x'|               = 6.65e-06 ≰ 0.0e+00
    |x - x'|/|x'|          = 2.48e-06 ≰ 0.0e+00
    |f(x) - f(x')|         = 4.94e-11 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 4.00e-11 ≰ 0.0e+00
    |g(x)|                 = 1.37e-11 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    4
    f(x) calls:    10
    ∇f(x) calls:   10

θ_opt = Optim.minimizer(result)[1]
energy_opt = Optim.minimum(result)
@info "VQE: θ* = $θ_opt, E* = $energy_opt"

┌ Info: VQE: θ* = -2.6779450445907376, E* = -1.2360679774997898

More advanced hamiltonians

Fermionic operators can be transformed into Ising-type hamiltonians in several ways. The most common are the Jordan-Wigner and Bravyi-Kitaev transformations.

In this example we only had a single qubit but it can be extended to multiple qubits.