

In [ ]:

    
from __future__ import print_function, division
import qutip as qt
import numpy as np
import scipy.linalg as sp
import math
import cmath
π = math.pi



In [ ]:

    
def preprocess(mat, vec):

    if mat.shape[0] != mat.shape[1] or mat.shape[0] != vec.shape[0]:
        raise Exception("Matrix A should be square and b should have a matching dimension!")
        
    if not mat.isherm:
        zero_block = np.zeros(mat.shape)
        mat = qt.Qobj(np.bmat([[zero_block, mat.full()] , 
                               [mat.dag().full(), zero_block]]))
        vec = qt.Qobj(np.hstack([b_init.full().flatten(), zero_block[0]]))
    
    ### Normalise b and remember normalisation factor
    if vec.norm() != 1:
        nf = vec.norm()
        vec = vec / nf
    else:
        nf = 1

    return mat, vec, nf



In [ ]:

    
def qft(N):
    mat = 1 / math.sqrt(N) * qt.Qobj([[cmath.exp(1j * 2 * π * i * j / N)
                                       for i in range(N)] for j in range(N)])
    return mat

def angle(k, t0, C):
    return math.asin(- C * t0 /(2 * π * k))    # Sign is for appropriate phases in the rotation, but should not affect the solution


def rot(k, t0, C):
    return qt.Qobj([[math.cos(angle(k, t0, C)), math.sin(angle(k, t0, C))],
                     [- math.sin(angle(k, t0, C)), math.cos(angle(k, t0, C))]])

# def inv(Quantobj):
#     mat = Quantobj.full()
#     matinv = np.linalg.inv(mat)
#     invQobj = qt.Qobj(matinv)
#     invQobj.dims = Quantobj.dims
#     return invQobj



In [ ]:

    
A_init = qt.Qobj([[0.2, 0.1],[0.1, 0.2]])
b_init = qt.Qobj([[0.2] , [0.7]])
prec = 1000    # Choose the dimension of the ancillary state

# Produce a Hermitian matrix A given the problem matrix A_int
# and a unit-length vector b given the problem vector b_int 
A, b, nf1 = preprocess(A_init, b_init)

eigenvals, eigenvecs = A.eigenstates()  # The eigendecomposition of A

# Condition number, of use for estimating constants
κ = eigenvals.max() / eigenvals.min()

# Additive error achieved in the estimation of |x>, of use for estimating constants
ϵ = 0.2

STEP 1: STATE PREPARATION OF b

[HHL: page 2, column 2, center]

Write b in the eigenbasis of A



In [ ]:

    
# Create the matrix that diagonalizes A
diagonalizer = qt.Qobj(np.array([eigenvecs[i].full().T.flatten()
                                 for i in range(len(eigenvals))]))
b = diagonalizer * b
A = diagonalizer.dag() * A * diagonalizer

STEP 2: QUANTUM SIMULATION AND QUANTUM PHASE ESTIMATION

[HHL: page 2, column 2, bottom]



In [ ]:

    
T = prec
t0 = κ / ϵ    # It should be O(κ/ϵ), whatever that means
ψ0 = qt.Qobj([[math.sqrt(2 / T) * math.sin(π * (τ + 0.5) / T)] 
                for τ in range(T)])

# Order is b, τ, and then ancilla
evo = qt.tensor(qt.identity(A.shape[0]), qt.ket2dm(qt.basis(T, 0))) 

for τ in range(1, T):
    evo += qt.tensor((1j * A * τ * t0 / T).expm(), qt.ket2dm(qt.basis(T, τ)))

ψev = evo * qt.tensor(b, ψ0)

ftrans = qt.tensor(qt.identity(b.shape[0]), qft(T))

ψfourier = ftrans * ψev    # This is Eq. 3 in HHL

STEP 2-1: Making true the assumption of perfect phase estimation



In [ ]:

    
# w = (ψfourier[:T] / b[0]).argmax()
# prj = qt.ket2dm(qt.basis(T, w))

# ψfourier = qt.tensor(qt.identity(b.shape[0]), prj) * ψfourier

STEP 3-1: Conditional Rotation of Ancilla



In [ ]:

    
total_state = qt.tensor(ψfourier, qt.basis(2, 0))   # Add ancilla for swapping

C = 1 / κ    # Constant, should be O(1/κ)

# Do conditional rotation only on τ and ancilla
rotation = qt.tensor(qt.ket2dm(qt.basis(T, 0)), qt.identity(2)) 

for τ in range(1, T):
    rotation += qt.tensor(qt.ket2dm(qt.basis(T, τ)), rot(τ, t0, C))
    
final_state = qt.tensor(qt.identity(b.shape[0]), rotation) * total_state

STEP 3-2: Post-selection on $\left|1\right\rangle$ on the ancilla register

We perform the post-selection by projecting onto the desired ancilla state and later tracing out the ancilla and the $\left|\psi_0\right\rangle$ registers.



In [ ]:

    
projector = qt.tensor(qt.identity(b.shape[0]), qt.identity(T), qt.ket2dm(qt.basis(2, 1)))
postsel = projector * final_state
prob1 = qt.expect(projector, final_state)

# Trace out ancilla and τ registers, leaving only the b register
finalstate = qt.ket2dm(postsel).ptrace([0]) / prob1
finalstate.eigenenergies()

$\left|finalstate\right\rangle$ is essentially a pure state (it should be if all the process was perfect), so now we just isolate that main part, that is our $\left|x\right\rangle$ (after a transformation if our original A was not Hermitian) in the basis that diagonalizes $A$.



In [ ]:

    
fsevls, fsevcs = finalstate.eigenstates()
x = math.sqrt(fsevls.max()) * fsevcs[fsevls.argmax()]
x

This state is supposed to be $\left|x\right\rangle=A^{-1}\left|b\right\rangle=\sum_j \beta_j\lambda_j^{-1}\left|u_j\right\rangle$, although I don't understand why it has complex numbers.

IMPORTANT NOTE: If A_init is not Hermitian, the output $\left|x\right\rangle$ has dimension double than the input b_init. It is to be interpreted as the representation of the vector (0, x) in the basis that diagonalizes $C = [[0,\,A], [A^{\dagger},\, 0]]$

Eq. (A26) in HHL is wrong, but Eq. (A27) is right. The global sign in Eq. (A29) is wrong, as well as the $\sqrt{2}$ (it should be in the denominator). The problem with the $\sqrt{2}$ carries until Eq. (A31)