For the VQE algorithm to work it is assumed that the Hamiltonian can be decomposed into a sum of Pauli constituents. Inspired by an implementation for 2 qubits on the quantum Slack channel I thought I'd try to implement a general decomposition algorithm in Julia.
From Wikipedia:
...the Pauli matrices are a set of three 2 × 2 complex matrices which are Hermitian and unitary.
...together with the identity matrix $I$ (sometimes considered as the zeroth Pauli matrix $σ_0$), the Pauli matrices form a basis for the real vector space of 2 × 2 Hermitian matrices. This means that any 2 × 2 Hermitian matrix can be written in a unique way as a linear combination of Pauli matrices, with all coefficients being real numbers.
First we need to define the Pauli matrices.
$$ \sigma_0 = I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}, \sigma_1 = \sigma_x = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \sigma_2 = \sigma_y = \begin{pmatrix} 0 &-i \\ i & 0 \end{pmatrix}, \sigma_3 = \sigma_z = \begin{pmatrix} 1 & 0 \\ 0 &-1 \end{pmatrix}. $$In Julia code this becomes:
σ₀=[1 0; 0 1]
σ₁=[0 1; 1 0]
σ₂=[0 -im; im 0]
σ₃=[1 0; 0 -1]
σ = [σ₀, σ₁, σ₂, σ₃];
Taking the wikipedia article at face value, for a 2x2 matrix $H$ the decomposition would be:
$$ H = \sum_{i=0,x,y,z} a_{i} \sigma_i, \quad a_{i,j} = \frac{1}{2} tr\left( \sigma_i H \right) $$Let's try it with a random matrix with real entries.
using LinearAlgebra # For the trace and kron operations
temp = rand(2,2)
H = 0.5(temp+temp')
a = [tr(σ[i]*H) for i in 1:4]/2
@assert H ≈ sum(a[i]*σ[i] for i in 1:4)
So it works!
What about a Hermitian matrix with complex entries?
temp = rand(ComplexF64, 2,2)
H = 0.5(temp+temp')
a = [tr(σ[i]*H) for i in 1:4]/2
@assert H ≈ sum(a[i]*σ[i] for i in 1:4)
That works too!
What about larger dimensions? From this blog post I found the formula:
$$ H = \sum_{i,j=0,x,y,z} a_{i,j} \left( \sigma_i \otimes \sigma_j \right), \quad a_{i,j} = \frac{1}{4} tr\left[\left( \sigma_i \otimes \sigma_j \right) H \right] $$This means that $a$ is a $4\times4$ matrix containing all the combinations of two Pauli matrices. This translates very easily into Julia.
⊗ = kron
temp = rand(ComplexF64, 4,4)
H = 0.5(temp+temp')
a = [tr((σ[i]⊗σ[j])*H) for i in 1:4, j in 1:4]/4
@assert H ≈ sum(a[i,j]*(σ[i]⊗σ[j]) for i in 1:4, j in 1:4)
So far so good. Now we want to generalise it to $n$ dimensions.
The pattern seems clear. To find the coefficients, $a_{1,2,3,...}$, calculate the tensor product of the $n$ combinations of Pauli matrices. Multiply by the matrix to be decomposed and then take the trace and divide by the normalisation factor, $2^n$.
The result will be an $n$ dimensional array with $4^n$ entries. Julia's CartesianIndices
type, gives us a way to iterate over all the dimensions in a linear way.
""" The tensor product of pauli matrices for the given index of arbitrary length [i,j,k,...]"""
pauli_product(idx) = foldr(⊗, σ[idx[j]] for j in 1:length(idx))
""" Returns multidimensional array of Pauli coefficients, depending on the size of H """
function pauli_decomposition(H)
@assert ishermitian(H)
n = Int(log(2, size(H, 1)))
a = Array{Float64}(undef, fill(4, n)...) # 4^n array
for idx in CartesianIndices(a)
a[idx] = real(tr(pauli_product(idx)*H)) # Ignore numerical error in imaginary part
end
a/2^n # Normalisatin factor
end
""" Sums the pauli terms with the given coefficients """
function pauli_sum(a)
sum(a[idx]*pauli_product(idx) for idx in CartesianIndices(a))
end;
Time to test it.
n=5
temp = rand(ComplexF64, 2^n,2^n)
H = 0.5(temp+temp')
@time a = pauli_decomposition(H)
@assert H ≈ pauli_sum(a)
0.422849 seconds (1.31 M allocations: 105.146 MiB, 4.67% gc time)
Works on my machine.. I've tested it up to n=7 and it works perfectly!
(After n=7 the matrices the calculation takes quite a while and consumes the memory on my laptop.)
We can also test it with a matrix with a known result: $$H = \frac{1}{2}(I\otimes I-\sigma_x\otimes \sigma_x-\sigma_y\otimes \sigma_y+\sigma_z\otimes \sigma_z)$$
H = [1 0 0 0;
0 0 -1 0;
0 -1 0 0;
0 0 0 1]
a = pauli_decomposition(H)
4×4 Array{Float64,2}: 0.5 0.0 0.0 0.0 0.0 -0.5 0.0 0.0 0.0 0.0 -0.5 0.0 0.0 0.0 0.0 0.5
Which matches perfectly what we expected :)