Variational Quantum Linear Solver
====================================

This notebook is an implementation of the pennylane notebook into myqlm notebook developed  by Andrea Mari (Xanadu) who implemented "the Variational quantum linear solver" (VQLS) algorithm  originaly introduced by Bravo-Prieto et al. (2019) <https://arxiv.org/abs/1909.05820>.
In the following we have reproduced Andrea's explanations from her original notebook.
We thanks Jean-christophe Jaskula (ATOS) for his help.
In this notebook, unstead of using a gradient descent based method for optimizing the objective function as in the orginal implementation, we use COBYLA.



    

Introduction
------------------

We first define the problem and the general structure of a VQLS.
As a second step, we consider a particular case and we solve it explicitly with PennyLane.

## The problem


We are given a $2^n \times 2^n$ matrix $A$ which can be expressed as a linear
combination of $L$ unitary matrices $A_0, A_1, \dots A_{L-1}$, i.e.,

\begin{align}A = \sum_{l=0}^{L-1} c_l A_l,\end{align}

where $c_l$ are arbitrary complex numbers. Importantly, we assume that each of the
unitary components $A_l$ can be efficiently implemented with a quantum circuit
acting on $n$ qubits.

We are also given a normalized complex vector in the physical form of a quantum
state $|b\rangle$, which can be generated by a unitary operation $U$
applied to the ground state of $n$ qubits. , i.e.,

\begin{align}|b\rangle = U |0\rangle,\end{align}

where again we assume that $U$ can be efficiently implemented with a quantum circuit.

The problem that we aim to solve is that of preparing a quantum state $|x\rangle$, such that
$A |x\rangle$ is proportional to $|b\rangle$ or, equivalently, such that

\begin{align}|\Psi\rangle :=  \frac{A |x\rangle}{\sqrt{\langle x |A^\dagger A |x\rangle}} \approx |b\rangle.\end{align}


## Variational quantum linear solver


The approach used in a VQLS is to approximate the solution $|x\rangle$ with a variational quantum
circuit, i.e., a unitary circuit $V$ depending on a finite number of classical real parameters
$w = (w_0, w_1, \dots)$:

\begin{align}|x \rangle = V(w) |0\rangle.\end{align}

The parameters should be optimized in order to maximize the overlap between the quantum states
$|\Psi\rangle$ and $|b\rangle$. This suggests to define the following cost function:

\begin{align}C_G = 1- |\langle b | \Psi \rangle|^2,\end{align}

such that its minimization with respect to the variational parameters should lead towards the problem solution.

In her original implementation Andrea Mari  propose to replace $\color{blue}{|0\rangle\langle 0|}$ in the original objective fonction $C_G$:

\begin{align}C_G = 1- \frac{ \sum_{l, l'}  c_l c_{l'}^* \langle 0|  V^\dagger A_{l'}^\dagger U \color{blue}{|0\rangle \langle 0|} U^\dagger A_l  V |0\rangle}
    {\sum_{l,l'} c_l c_{l'}^* \langle 0| V^\dagger A_{l'}^\dagger A_l V |0\rangle} .\end{align}
    
with the following positive operator:

\begin{align}\color{blue}{P} =  \frac{1}{2} + \frac{1}{2n}\sum_{j=0}^{n-1} Z_j,\end{align}

where $Z_j$ is the Pauli $Z$ operator locally applied to the $j\rm{th}$ qubit. This gives a new cost function:

\begin{align}C_L = 1- \frac{ \sum_{l, l'}  c_l c_{l'}^* \langle 0|  V^\dagger A_{l'}^\dagger U \color{blue}{P} U^\dagger A_l  V |0\rangle}
    {\sum_{l,l'} c_l c_{l'}^* \langle 0| V^\dagger A_{l'}^\dagger A_l V |0\rangle},\end{align}

which satisfies

\begin{align}C_G \rightarrow 0   \Leftrightarrow C_L \rightarrow 0,\end{align}

and so we can solve our problem by minimizing $C_L$ instead of $C_G$.

Substituting the definition of $P$ into the expression for $C_L$ we get:

\begin{align}C_L
    &=& \frac{1}{2} - \frac{1}{2n} \frac{ \sum_{j=0}^{n-1} \sum_{l, l'}  c_l c_{l'}^* \langle 0|  V^\dagger A_{l'}^\dagger U Z_j U^\dagger A_l  V |0\rangle}
    {\sum_{l,l'} c_l c_{l'}^* \langle 0| V^\dagger A_{l'}^\dagger A_l V |0\rangle} \\
    &&\\
    &=& \frac{1}{2} - \frac{1}{2n} \frac{ \sum_{j=0}^{n-1} \sum_{l, l'}  c_l c_{l'}^* \mu_{l,l',j}}
    {\sum_{l,l'} c_l c_{l'}^* \mu_{l,l',-1}},\end{align}

which can be computed whenever we are able to measure the following coefficients

\begin{align}\mu_{l, l', j} = \langle 0|  V^\dagger A_{l'}^\dagger U Z_j U^\dagger A_l  V |0\rangle,\end{align}

where we used the convention that if $j=-1$,  $Z_{-1}$ is replaced with the identity.

In our implementation we don't use the hadamard test  as in the original notebook but take advantage of myQLM to provide the expectaion value of the Observable we have defined:  $ A_{l'}^\dagger U Z_j U^\dagger A_l$



### A simple example


In this notebook we reproduce the  simple example based on a system of 3 qubits  proposed in Andrea Mari's implementation.

\begin{align}
    A  &=  c_0 A_0 + c_1 A_1 + c_2 A_2 = \mathbb{I} + 0.2 X_0 Z_1 + 0.2 X_0, \\
    |b\rangle &= U |0 \rangle = H_0  H_1  H_2 |0\rangle,
\end{align}


where $Z_j, X_j, H_j$ represent the Pauli $Z$, Pauli $X$ and Hadamard gates applied to the qubit with index $j$.

This problem is computationally quite easy since a single layer of local rotations is enough to generate the
solution state, i.e., we can use the following simple ansatz:

\begin{align}
|x\rangle = V(w) |0\rangle = \Big [  R_y(w_0) \otimes  R_y(w_1) \otimes  R_y(w_2) \Big ]  H_0  H_1  H_2 |0\rangle.
\end{align}


In the code presented below we solve this particular problem by minimizing the local cost function $C_L$.
Eventually we will compare the quantum solution with the classical one.


```python
import time
import numpy as np
import importlib.util
from scipy.optimize import minimize

from qat.lang.AQASM import Program,QRoutine,RY,RZ,H,CNOT,Z,PH
from qat.core import Observable, Term
from qat.plugins import ObservableSplitter
from qat.qpus import get_default_qpu

class Timer():
    text: str = "Elapsed time: {:0.4f} seconds"
    _start_time: float = 0

    def start(self) -> None:
        """Start a new timer"""        
        self._start_time = time.perf_counter()

    def stop(self) -> float:
        """Stop the timer, and report the elapsed time"""
        elapsed_time = time.perf_counter() - self._start_time
        self._start_time = None

        return elapsed_time
    
    def __enter__(self):
        """Start a new timer as a context manager"""
        self.start()
        return self

    def __exit__(self, *exc_info):
        """Stop the context manager timer"""
        print(self.text.format(self.stop()))


qpu = ObservableSplitter() | get_default_qpu()

def U_b(qprog,qubits,n_qubits):
    """Unitary matrix rotating the ground state to the problem vector |b> = U_b |0>."""
    for idx in range(n_qubits):
        qprog.apply(H,qubits[idx])

H_mat = np.array([[1,1],[1,-1]])/np.sqrt(2)
X_m = np.array([[0,1],[1,0]])
Z_m = np.array([[1,0],[0,-1]])

def apply_n(f, m, n):
    r = m
    for i in range(n-1):
        r = f(r,m)
    return r

def variational_block(qprog,qubits,n_qubits,weights):
    """Variational circuit mapping the ground state |0> to the ansatz state |x>."""
    # We first prepare an equal superposition of all the states of the computational basis.
    for idx in range(n_qubits):
        qprog.apply(H,qubits[idx])

    # A very minimal variational circuit.
    for idx, element in enumerate(weights):
        var = qprog.new_var(float, '\\theta'+str(idx))
        qprog.apply(RY(var), qubits[idx])

def psi_norm(c):
    """Returns the normalization constant <psi|psi>, where |psi> = A |x>."""

    return abs(qpu.submit(c.to_job(job_type='OBS', observable=A*A)).value)

def bind(c, w):
    #print({'\\theta'+str(idx):e for idx, e in enumerate(w)})
    return c.bind_variables({'\\theta'+str(idx):e for idx, e in enumerate(w)})

# Expectation value of the Observable  and corresponding cost function
def cost_loc(weights: np.array,
             O:       Observable
            )->float:
    c = bind(parametrized_circuit, weights)
    num = qpu.submit(c.to_job(job_type='OBS', observable=O, nbshots=0)).value
    res = 0.5 - 0.5 * abs(num) / (n_qubits * psi_norm(c))
    cost_history.append(res)
    return res

n_qubits = 3       # Number of system qubits.
n_shots = 10 ** 2  # Number of quantum measurements.
steps = 30         # Number of optimization steps
eta = 0.8          # Learning rate
q_delta = 0.001    # Initial spread of random quantum weights
rng_seed = 0       # Seed for random number generator

np.random.seed(rng_seed)
w0 = np.random.randn(n_qubits) * q_delta 

# define the observables 
A = Observable(n_qubits, pauli_terms=[Term(0.2, "XZ", [0,1]),Term(0.2, "X", [0])], constant_coeff=1)
U = Observable(n_qubits, matrix=apply_n(np.kron, H_mat, 3)/(2**(3-1)))
Z = Observable(n_qubits, pauli_terms=[Term(1, "Z", [i]) for i in range(n_qubits)])
# note: by defining O will give a ~4x speedup in the algotihm
O=A*U*Z*U*A


# Main program
qprog=Program()
qubits=qprog.qalloc(n_qubits)
variational_block(qprog,qubits,n_qubits,w0)
parametrized_circuit = qprog.to_circ()
cost_history = []
with Timer():
    w = minimize(cost_loc,w0,O, method='COBYLA',
                           constraints=(), tol=0.0001, callback=None, 
                           options={'rhobeg': 1.0, 'maxiter': 50, 'disp': True, 'catol': 0.001})

# Preparation of the quantum solution
# ------------------------------------
# Given the variational weights ``w`` that we have previously optimized,
# we can generate the quantum state $|x\rangle$. By measuring $|x\rangle$
# in the computational basis we can estimate the probability of each basis state.
# For this task, we initialize a new PennyLane device and define the associated *qnode* circuit.
# To estimate the probability distribution over the basis states we first take ``n_shots``
# samples and then compute the relative frequency of each outcome.
def prob_dist(r):
    n = 2**len(r.raw_data[0].state.bitstring)
    p = np.zeros(n)
    for j, e in enumerate(r.raw_data):
        p[sum([int(s)*2**int(i) for i,s in enumerate(e.state.bitstring)])] = e.probability
        
    return p

q_probs = prob_dist(qpu.submit(bind(parametrized_circuit,w.x).to_job()))
q_probs_init = prob_dist(qpu.submit(bind(parametrized_circuit,w0).to_job()))

# Classical algorithm
# ----------------------
# To solve the problem in a classical way, we use the explicit matrix representation in
# terms of numerical NumPy arrays.c = np.array([1.0, 0.2, 0.2])
c = np.array([1.0, 0.2, 0.2])
Id = np.identity(2)
Z_mat = np.array([[1, 0], [0, -1]])
X_mat = np.array([[0, 1], [1, 0]])

A_0 = np.identity(8)
A_1 = np.kron(np.kron(X_mat, Z_mat), Id)
A_2 = np.kron(np.kron(X_mat, Id), Id)

A_num = c[0] * A_0 + c[1] * A_1 + c[2] * A_2
b = np.ones(8) / np.sqrt(8)

# We can print the explicit values of  𝐴  and  𝑏 :
print("A = \n", A_num)
print("b = \n", b)

# The solution can be computed via a matrix inversion:
A_inv = np.linalg.inv(A_num)
x = np.dot(A_inv, b)

# Finally, in order to compare x with the quantum state $|x\rangle$, we normalize
# and square its elements.
c_probs = (x / np.linalg.norm(x)) ** 2

# Comparison
# Let us print the classical result.
print("classical result")
print("x_n^2 =\n", c_probs)

# The previous probabilities should match the following quantum state probabilities.
print("VQLS result")
print("|<x|n>|^2=\n", q_probs)
```

    Elapsed time: 5.9850 seconds
    A = 
     [[1.  0.  0.  0.  0.4 0.  0.  0. ]
     [0.  1.  0.  0.  0.  0.4 0.  0. ]
     [0.  0.  1.  0.  0.  0.  0.  0. ]
     [0.  0.  0.  1.  0.  0.  0.  0. ]
     [0.4 0.  0.  0.  1.  0.  0.  0. ]
     [0.  0.4 0.  0.  0.  1.  0.  0. ]
     [0.  0.  0.  0.  0.  0.  1.  0. ]
     [0.  0.  0.  0.  0.  0.  0.  1. ]]
    b = 
     [0.35355339 0.35355339 0.35355339 0.35355339 0.35355339 0.35355339
     0.35355339 0.35355339]
    classical result
    x_n^2 =
     [0.08445946 0.08445946 0.16554054 0.16554054 0.08445946 0.08445946
     0.16554054 0.16554054]
    VQLS result
    |<x|n>|^2=
     [0.08447514 0.08445027 0.16557213 0.16552339 0.08446807 0.0844432
     0.16555827 0.16550953]



```python

```

c80a548b98dc556839212879eccb0ca8d55fec2c
ipynb
Jupyter Notebook
notebooks/vqls/VQLS_notebook.ipynb
myQLM/myqlm-contrib
6e00c349b857a68b9462e4d8f18de7b391b3910f
Apache-2.0
2020-09-09T09:44:33.000Z
2021-09-16T07:10:53.000Z
notebooks/vqls/VQLS_notebook.ipynb
myQLM/myqlm-contrib
6e00c349b857a68b9462e4d8f18de7b391b3910f
Apache-2.0
notebooks/vqls/VQLS_notebook.ipynb
myQLM/myqlm-contrib
6e00c349b857a68b9462e4d8f18de7b391b3910f
Apache-2.0
2020-07-06T10:01:09.000Z
2021-09-17T04:57:22.000Z
Qwen/Qwen-72B
1. YES
2. YES