Someone said that the best way of learning is to do; others said the best way to explain. So here I try to explain what I learned building VQE from scratch naively.

• If it's the first time you've heard of eigenvalue and Pauli matrices:

Variation algorithms are some hybrid cars, meaning part of the workload of the classic processor is replaced by quantum computers. The task the algorithm is supposed to fulfill is to calculate the minimum eigenvalue of an operator (sometimes this can have a physical meaning, like energy, for instance).

For starters, to determine the minim eigenvalue, we need the operator of which the minimum eigenvalue needs to be determined. For instance, consider the operator below:

H=-1.05 * np.kron(I,I) + 0.39 * np.kron(I,Z) - 0.39 * np.kron(Z,I) - 0.01 * np.kron(Z,Z) + 0.18 * np.kron(X,X)

When we composed the above operator, we composed it from a sum of Pauli matrices products.

Usually, we do not know the way our operator is decomposed, and as such we need to build a function to do it for us.

Think of how we would solve the problem if we were dealing with vectors. We would project the vectors on each axis. In this situation we have to do the exact same thing, except, we will use matrices. We need to start by building our matrix basis set, after which we need to project the matrix of our observable on the base elements, for which we will use Hilbert Schmidt product:

# Make Hilbert Schmidt product betwen matrices mat1 and mat2
def HS(mat1,mat2): 
    return(np.dot(mat1.conjugate().transpose(), mat2)).trace()

Of course, before this, we need to build the base elements. If, for example, our operator has a size of 2^3, each base element can be written as a tensor product of X, Z, Y matrices. The method below is meant to do exactly that: it builds the base and projects H on the basic elements.

# Decompose  observable in matrices create by kron. product of Pauli matrices (h-coeficient , h_label-kron product structure )
def decompose(O):
    size=len(O)
    nr_pauli=np.log2(size)
    norm_fact=1/(2**nr_pauli)

    elements=itertools.product(indice,repeat=int(nr_pauli))
    h_label=[]
    h=[]
    for i in elements:
        label=''
        matrice=pauli[i[0]]
        for  j in i :
            label=label+labels[indice[j]]
        for j in range(int(nr_pauli)-1):
            matrice=np.kron(matrice,pauli[i[j+1]])
        #print(matrice)
        h_label.append(label)
        h.append(norm_fact*HS(matrice,O))

    return h,h_label

Now we can write the expected value of H as a sum of the expected value of the projections

• $<\psi|H|\psi> = -1.05 *<\psi|I\otimes I|\psi>+0.39 *<\psi|I\otimes Z|\psi>-0.39 *<\psi|Z\otimes I|\psi>-0,01 *<\psi|Z\otimes Z|\psi>+0.18 *<\psi|X\otimes X|\psi>$
• $<\psi|H|\psi> = -1.05 *1+0.39 *<\psi|I\otimes Z|\psi>-0.39 *<\psi|Z\otimes I|\psi>-0,01 *<\psi|Z\otimes Z|\psi>+0.18 *<\psi|X\otimes X|\psi>$

The following critical element of VQE help us move through space; thus, we need to be extremely careful of choosing the ansatz because sometimes we do not wish to explore the entire space, so it would be better not to select ansatz points that are too hardware expensive. One ansatz example is RyRz ansatz.

# Ansatz RZRY -> these is just an example of ansatz 
# Look at ref [2]

def ryrz(H,parm,depth):
    
    nrq=int(np.log2(len(H)))
    ansatz=QuantumCircuit(nrq,nrq)
    it=iter(parm)
    
        
    for i in range(nrq):
        
        ansatz.ry(next(it),i)
        ansatz.rz(next(it),i)
            
            
    for g in range(depth):
        
        for l in range(nrq-1):
            ansatz.cx(l,l+1)
            
        for i in range(nrq):
            ansatz.ry(next(it),i)
            ansatz.rz(next(it),i)
          
        
            
    return ansatz

depth=2
pi=np.pi
parm=np.array([pi/2, pi/2, pi/2,pi/2,  pi/2, pi/2, pi/2,pi/2,pi/2, pi/2, pi/2,pi/2,pi/2, pi/2, pi/2,pi/2])
ansatz=ryrz(H,parm,depth=depth)
ansatz.draw('mpl')