HHL algorithm

I am going to study quantum algorithms in my own way using qiskit. Since this is a record of my personal study, I may have left out a lot of explanations.

I am following the qiskit website.

The reason why qiskit was categorized as a Rec (recommendation system) in my blog was all to understand HHL. Currently, I am interested in recommender systems and am developing them. Linear equations are used with high probability when solving mathematical models using computers, and it is important for recommendation systems to extract features with high user engagement from the User-Item matrix.

Therefore, I started reviewing quantum algorithms with the goal of using a quantum computer to solve simultaneous linear equations at high speed, and I have finally reached my goal.


  • The file in jupyter notebook format is here .

google colaboratory

  • If you want to run it on google colaboratory here

Author’s environment

ProductName: Mac OS X
ProductVersion: 10.14.6
BuildVersion: 18G103
Python -V
Python 3.8.5

Import the basic libraries and check their versions.

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

import qiskit
import json

import matplotlib.pyplot as plt
import numpy as np
import math

from qiskit import IBMQ, Aer, transpile, assemble
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister

from qiskit.visualization import plot_histogram

{'qiskit-terra': '0.17.4',
 'qiskit-aer': '0.8.2',
 'qiskit-ignis': '0.6.0',
 'qiskit-ibmq-provider': '0.13.1',
 'qiskit-aqua': '0.9.1',
 'qiskit': '0.26.2',
 'qiskit-nature': None,
 'qiskit-finance': None,
 'qiskit-optimization': None,
 'qiskit-machine-learning': None}

Conjugate gradient method

By way of review, we will review the classical algorithm, the conjugate gradient method. A linear equation with positive definite matrix, $A$, as coefficients, and

$$ A \boldsymbol{x}=\boldsymbol{b} $$

Using an iterative method, we can numerically solve for $x$. Since it is an iterative method, an error $(\epsilon)$ is required to determine the end of the calculation.

The matrix $A$,$x$,$b$ is as follows.

$$ A = \left(\begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1 n} \\ a_{21} & a_{22} & \cdots & a_{2 n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n 1} & a_{n 2} & \cdots & a_{n n} \end{array}\right),\quad x=\left(\begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{array}\right), \quad b=\left(\begin{array}{c} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \end{array}\right) $$

Written in matrix tabular form, it is as follows.

$$ \left(\begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1 n} \\ a_{21} & a_{22} & \cdots & a_{2 n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n 1} & a_{n 2} & \cdots & a_{n n} \end{array}\right)\left(\begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{array}\right)=\left(\begin{array}{c} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \end{array}\right) $$

Next, consider the function $f(x)$, defined as follows.

$$ f(\boldsymbol{x})=\frac{1}{2}(\boldsymbol{x}, A\boldsymbol{x})-(\boldsymbol{b}, \boldsymbol{x}) $$

$(-,-)$ is an operator to calculate the inner product of vectors.

$$ (\boldsymbol{x}, \boldsymbol{y})=\boldsymbol{x}^{T} \boldsymbol{y}=\sum_{i=1}^{n} \boldsymbol{x}{i} \boldsymbol{y}{i} $$

When displayed in components, it looks like the following.

$$ f(x)=\frac{1}{2} \sum_{i=1}^{n} \sum_{j=1}^{n} a_{i j} x_{i} x_{j}-\sum_{i=1}^{n} b_{i} x_{i} $$

Now, differentiating by $x_k$, we get

$$ \frac{\partial f(x)}{\partial x_{k}}=\frac{1}{2} \sum_{i=1}^{n} a_{i k} x_{i}+\frac{1}{2} \sum_{j=1}^{n} a_{k j} x_{j}-b_{k} $$

It becomes

Since $A$ is a Hermitian matrix.

$$ \frac{\partial f(x)}{\partial x_{i}}=\sum_{j=1}^{n} a_{i j} x_{j}-b_{i}=0 $$

This becomes

Generalizing this, we get

$$ \nabla f(x)=\left(\begin{array}{c} \frac{\partial f}{\partial x_{1}} \\ \vdots \\ \frac{\partial f}{\partial x_{n}} \end{array}\right)=A\boldsymbol{x}-b = 0 $$

This shows that finding the minimum value of $x$ for the function $f(x)$ is the same as solving for $A\boldsymbol{x}-b = 0$$.


As mentioned above, the conjugate gradient method (CG method) boils down to minimizing the function $f(x)$.

The normal gradient method is less efficient because it searches in the vertical direction of the contour line.

$$ -\nabla f(x)=\left(\begin{array}{c} \frac{\partial f}{\partial x_{1}} \\ \vdots \\ \frac{\partial f}{\partial x_{n}} \end{array}\right)=-b + A\boldsymbol{x} $$

The $f(x)$ forms an ellipse, and the idea is that if we apply a linear transformation to the ellipse and convert the ellipse to a circle, we can find an efficient solution because the vertical direction of the contour line is the shortest distance to the minimum point.

To do so, starting from a certain $x^{(0)}$, find the $x$ that is the minimum value according to the following asymptotic formula.

$$ x^{(k+1)}=x^{(k)}+\alpha_{k} p^{(k)} $$

Here, $p^{(k)}$$ is the direction vector to search for the solution. Also, $r^{(k)}$$ is the residual vector with the solution.

$$ r^{(0)}=b-A x^{(0)}, \quad {p}^{(0)}={r}^{(0)} $$

Find the $\alpha_k$ that minimizes $f(x^{(k+1)})$ at the $k+1$th step.

$$ \begin{aligned} &f(x^{(k+1)}) \\ &=f\left(x^{(k)}+\alpha_{k} p^{(k)}\right) \\ &=\frac{1}{2}\alpha_{k}{ }^{2}\left(p^{(k)}, A p^{(k)}\right)-\alpha_{k}\left(p^{(k)}, b-A x^{(k)}\right)+f\left(x^{(k)}\right) \end{aligned} $$

Since this is a simple quadratic function, we can find $\alpha_k$.

$$ \alpha_{k}=\frac{\left(p^{(k)}, b-A x^{(k)}\right)}{\left(p^{(k)}, A p^{(k)}\right)}=\frac{\left(p^{(k)}, r^{(k)}\right)}{\left(p^{(k)}) )}, A p^{(k)}\right)} $$

Next, to update $r$, we need $r^{(k+1)}$. $$ \begin{aligned} &r^{(k+1)}=b-A x^{(k+1)} \\ &r^{(k)}=b-A x^{(k)} \\ &r^{(k+1)}-r^{(k)}=A x^{(k+1)}-A x^{(k)}=\alpha_{k} A p^{(k)} \end{aligned} $$


$$ r^{(k+1)}=r^{(k)}-\alpha_{k} A p^{(k)} $$

This becomes

If this $$r^{(k+1)}$$ is below a certain threshold, the search is terminated.

$$ \left|\boldsymbol{r}^{(k+1)}\right|<\varepsilon $$

If it is other than below the threshold, we then seek $p^{(k+1)}$ to update $p$. $p^{(k+1)}$ consists of a vector of residuals $r^{(k+1)}$ and a constant multiple of $p^{(k)}$.

$$ p^{(k+1)}=r^{(k+1)}+\beta_{k} p^{(k)} $$

The constant $\beta_k$ is taken so that $p^{(k)}$ and $p^{(k+1)}$$ are conjugate to $A$.

$$ \left(p^{(k+1)}, A p^{(k)}\right)=\left(r^{(k+1)}+\beta_{k} p^{(k)}, A p^{(k)}\right)=\left(r^{(k+1)}, A p^{(k)}\right)+\beta_{k}\left(p ^{(k)}, A p^{(k)}\right) $$

From this, we get $$ \beta_{k}=-\frac{\left(r^{(k+1)}, A p^{(k)}\right)}{\left(p^{(k)}, A p^{(k)}\right)} $$

which becomes

This is repeated until the residuals converge below a certain threshold.

Here, we have

$$ \left(p^{(i)}, A p^{(j)}\right)=0 \quad (i \neq j) $$

So, I guess the name of the gradient method is conjugate gradient method, because the vector with conjugate relation is the direction vector.

Also, for the residual vector

$$ \left(r^{(i)}, r^{(j)}\right)=0 \quad (i \neq j) $$

This is the orthogonal relationship.

Since the number of $r_k$ that are first-order independent is only $N$ at most, the calculation converges within $N$ times.

For more information on the preprocessed conjugate gradient method, Krylov subspaces, and convergence, please refer to the following.

Computational complexity

$s$ is the fraction of non-zero elements in matrix $A$, $\kappa$ is the ratio of the largest to the smallest eigenvalue of matrix $A$, $\displaystyle \left| \frac{\lambda_{max}}{\lambda_{min}}\right|$, and $\epsilon$ is the precision.

$$ O(N s \kappa \log (1 / \varepsilon)) $$

This is in the form proportional to $N$, so if $s\sim \log N$, it will be proportional to $N\log N$.

HHL Algorithm

Assumptions of HHL

HHL is an algorithm that is based on certain assumptions.

  • There are effective oracles that implement loading
  • Hamiltonian simulation and solution functions can be computed
  • $A$ is a Hermitian matrix
  • Can be computed in $\mathcal{O}\left(\log (N) s^{2} \kappa^{2} / \epsilon\right)$.
  • While the classical algorithm returns the complete solution, HHL only approximates the function that gives the solution vector.

Outline of HHL

Mapping to quantum circuits

To solve a simultaneous linear equation with a quantum algorithm, we need to map $Ax=b$ to a quantum circuit. One way to do this is to map the $i$th component of $b$ to the amplitude of the $i$th ground state of the quantum state $|b\rangle$. And, of course, this requires the normalization $\displaystyle \sum_i |b_i|^2 = 1$.

$$ Ax=b \rightarrow A|x\rangle = |b\rangle $$

This becomes.

Spectral decomposition

Since $A$ is a Hermitian matrix, it can be spectrally decomposed. We had implicitly assumed the Hermitianity of $A$. $$ A’=\left(\begin{array}{ll} 0 & A \\ A & 0 \end{array}\right) $$

Then, $A’$ is a Hermitian matrix, which is fine. Thus, $A$ can be expanded using the eigenvector $|u_i\rangle$ and its eigenvalue $\lambda_i$ as follows.

$$ A=\sum_{j=0}^{N-1}\lambda_{j}\left|u_{j}\right\rangle\left\langle u_{j}\right| $$

Therefore, the inverse matrix is as follows.

$$ A^{-1}=\sum_{j=0}^{N-1} \lambda_{j}^{-1}\left|u_{j}\right\rangle\left\langle u_{j}\right| $$

Since $u_i$ is an eigenvector of $A$, $|b\rangle$ can be represented by its superposition. Perhaps this is a strong motivation to use quantum computers.

$$ |b\rangle=\sum_{j=0}^{N-1} b_{j}\left|u_{j}\right\rangle $$

Normally, it is not possible to prepare $|b\rangle$ in this form without being able to calculate the eigenvectors of $A$, but a quantum computer can prepare this state automatically by reading $|b\rangle$.

In the end, we will use the quantum computer to obtain the following equation.

$$ |x\rangle=A^{-1}|b\rangle=\sum_{j=0}^{N-1}\lambda_{j}^{-1} b_{j}\left|u_{j}\right\rangle $$

1. Load the data

Load the amplitude of each data of the target $|b\rangle$ into the qubit $n_b$.

$$ |0\rangle_{n_{b}} \mapsto|b\rangle_{n_{b}} $$

2. Applying QPE

We use quantum phase estimation to estimate the phase of $|b\rangle$ of the unitary operator $U=e^{i A t}$.

The $|b\rangle$ is as follows, as described above.

$$ |b\rangle=\sum_{j=0}^{N-1} b_{j}\left|u_{j}\right\rangle $$

The $U$ can be expanded and organized as follows.

$$ \begin{aligned} &U=e^{i A t}=\sum_{k=0}^{\infty} \frac{\left(i A t\right)^{k}}{k !} \\ &=\sum_{k=0}^{\infty}\frac{\left(i A t\right)^{k}}{k !} \left(\sum_{j=0}^{N-1}\lambda_{j}|u_{j}\rangle\langle u_{j} |\right)^{k}\\ &=\sum_{k=0}^{\infty}\frac{(i t)^{k}}{k !} \sum_{j=0}^{N-1}\lambda_{j}^{k}\left|u_{j}\right\rangle \langle u_{j} |\\ &=\sum_{j=0}^{N-1}\left(\sum_{k=0}^{\infty}\frac{\left(i t\right)^{k}}{k_{i}}\lambda_{j}^{k}\right)|u_{j}\langle u_{j}| \\ &=\sum_{r=0}^{N-1} e^{i \lambda_{j}t}\left|u_{j}\right\rangle\left\langle u_{j}\right| \end{aligned} $$

If we let $U$ act on $|b\rangle$.

$$ \begin{aligned} U|b\rangle &=U\left(\sum_{j=0}^{N-1} b_{j}\left|u_{j}\right\rangle\right) \\ &=\sum_{j^{\prime}=0}^{N-1} e^{i \lambda j^{\prime} t}\left|h_{j^{\prime}}\right\rangle\left\langle h_{j^{\prime}}\right| \cdot\left(\ sum_{j=0}^{N-1} b_{j}\left|h_{j}\right\rangle\right) \\ &=\sum_{j=0}^{N-1} b_{j} e^{i \lambda_{j} t}\left|u_{j}\right\rangle \end{aligned} $$

Then, using quantum phase estimation, we can find the quantum state $|\tilde{\lambda_j}\rangle_{n_l}$ of $\lambda_j$.

The $\tilde{\lambda_{j}}$ is the $n_l$-bit binary approximation to $\displaystyle 2^{n_{l}}\frac{\lambda_{j} t}{2\pi}$.

Assuming $t=2\pi$ and $\lambda_l$ can be represented exactly in $n_l$ bits, quantum phase estimation can be expressed as follows.

$$ \operatorname{QPE}\left(e^{i A 2 \pi}, \sum_{j=0}^{N-1} b_{j}|0\rangle_{n_{l}}\left|u_{j}\right\rangle_{n_{b}}\right)=\sum_{j=0}^{N-1} b {j}\left|\lambda{j}\right\rangle_{n_{l}}\left|u_{j}\right\rangle_{n_{b}} $$

3. Use of auxiliary qubits

Although it is descent, we will consider how to use a control rotation gate to extract $|v_i\rangle$ as the amplitude of a qubit.

In this section, we refer to the following references.

We consider two qubits. The first is the control bit of the control rotation gate, and the second is the target bit of the control rotation gate.

Tangentially, input the state of $\left|b\left(\frac{1}{\pi}\cos ^{-1} v_{i}\right)\right\rangle$ to the control bit, and $|0\rangle$ to the target bit. and input $|0\rangle$ for the target bit. This target bit becomes the auxiliary qubit.

The $b(\cdots)$ indicates that it is in binary notation, and the $b_k(\cdots)$ represents the value of the $k$th bit in binary notation. To summarize, we have the following.

$$ \begin{aligned} &b\left(\frac{1}{\pi} \cos ^{-1}\left(v_{i}\right)\right)=d_{0} d_{1} d_2 \cdot \cdots d_{m-1} \\ &b_{k}\left(\frac{1}{\pi} \cos ^{-1}\left(v_{i}\right)\right)=d_k \\ &\frac{2}{\pi} \cos ^{-1}\left(v_{i}\right)=\frac{d_{0}}{2}+\frac{d_{1}}{4}+\frac{d_{2}}{8} \cdots \frac{d_{m-1}}{2^{m}} \quad \left(0 \leqq \frac{2}{\pi} \cos ^{-1}\left(v_{i}\right) \leqq 1\right) \end{aligned} $$

From this, we get

$$ \begin{aligned} \frac{2}{\pi} \cos ^{-1}\left(v_{i}\right)&=\frac{1}{2} b_{0}\left(\frac{1}{\pi} \cos ^{-1}\left(v_{1}\right)\right)+\frac{1}{4} b_{1}\left(\frac{1}{\pi} \cos ^{-1}\left(v_{i}\right)\right)+ \cdots \\ &=\sum_{k=0}^{m-1} b_{k}\left(\frac{1}{\pi} \cos ^{-1}\left(v_{n}\right)\right) 2^{-k-1} \cdots (1) \end{aligned} $$

Apply a control rotation gate with $\displaystyle b_{k}\left(\frac{2}{\pi}\cos ^{-1}\left(v_{i}\right)\right)$ as the control bit, $R_{y}\left(2^{-k-1}\pi\right)$. Consider the following

$$ \begin{aligned} &\prod_{k=0}^{m-1} R_{y}\left(b_{k}\left(\frac{1}{\pi} \cos ^{-1}\left(v_{i}\right)\right) 2^{-k-1}\pi\right)|0\rangle \\ &=R_{y}\left(\sum_{k=0}^{m-1} b_{k}\left(\frac{1}{\pi} \cos ^{-1}\left(v_{i}\right)\right) 2^{-k-1}\pi\right)|0\rangle \\ \end{aligned} $$

The rotation gate is as follows.

$$ R_{y}(\theta)|0\rangle=\cos \frac{\theta}{2}|0\rangle+\sin \frac{\theta}{2}|1\rangle $$

Using this equation and (1), we get

$$ \begin{aligned} &\cos\frac{1}{2}\left(\sum_{k=0}^{m-1} b_{k}\left(\frac{1}{\pi} \cos ^{-1}\left(v_{i}\right)\right) 2^{-k-1}\pi\right) \\ &=\cos\left(\frac{1}{2}\times\frac{2}{\pi} \cos ^{-1}\left(v_{i}\right)\times \pi\right) = v_i \end{aligned} $$

and we can use this to get

$$ R_{y}\left(\sum_{k=0}^{m-1} b_{k}\left(\frac{1}{\pi} \cos ^{-1}\left(v_{i}\right)\right) 2^{-k-1}\pi\right)|0\rangle=v_{i}|0\rangle+\sqrt{1-v_{i}^{2}}|1\rangle $$

We can then obtain Now, by setting $\displaystyle v_i = \frac{1}{\lambda_i}$, we can set the auxiliary qubit to

$$ \frac{1}{\lambda_{j}}|0\rangle + \sqrt{1-\frac{1}{\lambda_{j}^{2}}}|1\rangle $$

This can be calculated as Now that we have $\frac{1}{\lambda_j}$ as an amplitude, we can use it to solve the linear equation.

4. Use inverse quantum phase estimation

The inverse of quantum phase estimation gives us the following.

$$ \sum_{j=0}^{N-1} b_{j}|0\rangle_{n_{l}}\left|u_{j}\right\rangle_{n_{b}}\left(\frac{1}{\lambda_{j}}|0\rangle+\sqrt{1-\frac{1}{\lambda_{j}^{2}}}|1\rangle\right) $$

5. Measure the auxiliary qubit

Measure the auxiliary qubit, and if $|0\rangle$ is measured, then

$$ \left(\sqrt{\frac{1}{\sum_{j=0}^{N-1}\left|b_{j}\right|^{2} /\left|\lambda_{j}\right|^{2}}}\right) \sum_{j=0}^{N-1} \frac{b_{j}}{\lambda_{j}}|0\rangle_{n_{l}}\left|u_{j}\right\rangle_{n_{b}} $$

and the solution is in the form.

Computational complexity comparison

Quantum algorithm

$$ O(s \kappa \operatorname{poly} \log (s \kappa / \varepsilon))) $$

For a matrix $A$, if we can assume sparsity $(s \sim O(\operatorname{poly} \log N)))$, then

$$ O(s \kappa \operatorname{poly} \log (s \kappa / \varepsilon))) \sim O(s \kappa \operatorname{poly} \log N \operatorname{poly} \log (s \kappa / \varepsilon)) varepsilon)) $$

It becomes.

Conjugate gradient method

$$ O(N s \kappa \log (1 / \varepsilon)) $$

From this, we can expect an exponential speedup for the quantum algorithm.

Implement with qiskit

Let’s try to implement it according to the qiskit website.

As it turns out, my environment doesn’t give the results as per the site: I had a DepricateWarning, and I thought that might be it, so I tried various things, but the results didn’t match, so I’ll try to find the detailed cause later.

from qiskit import Aer
from qiskit.circuit.library import QFT
from qiskit.aqua import QuantumInstance, aqua_globals
from qiskit.quantum_info import state_fidelity
from qiskit.aqua.algorithms import HHL, NumPyLSsolver
from qiskit.aqua.components.eigs import EigsQPE
from qiskit.aqua.components.reciprocals import LookupRotation
from qiskit.aqua.operators import MatrixOperator
from qiskit.aqua.components.initial_states import Custom
import numpy as np

from qiskit import aqua

{'qiskit-terra': '0.17.4',
 'qiskit-aer': '0.8.2',
 'qiskit-ignis': '0.6.0',
 'qiskit-ibmq-provider': '0.13.1',
 'qiskit-aqua': '0.9.1',
 'qiskit': '0.26.2',
 'qiskit-nature': None,
 'qiskit-finance': None,
 'qiskit-optimization': None,
 'qiskit-machine-learning': None}
def create_eigs(matrix, num_ancillae, num_time_slices, negative_evals):
    ne_qfts = [None, None].
    if negative_evals:
        num_ancillae += 1
        ne_qfts = [QFT(num_ancillae - 1), QFT(num_ancillae - 1).inverse()].

    return EigsQPE(MatrixOperator(matrix=matrix),
                   evo_time=None, # np.pi*3/4, #None, # This is t, can set to: np.pi*3/4
def fidelity(hhl, ref):
    solution_hhl_normed = hhl / np.linalg.norm(hhl)
    solution_ref_normed = ref / np.linalg.norm(ref)
    fidelity = state_fidelity(solution_hhl_normed, solution_ref_normed)
    print("Fidelity:\t\t %f" % fidelity)
matrix = [[1, -1/3], [-1/3, 1]]
vector = [1, 0]]
orig_size = len(vector)
matrix, vector, truncate_powerdim, truncate_hermitian = HHL.matrix_resize(matrix, vector)

# Initialize eigenvalue finding module
eigs = create_eigs(matrix, 3, 100, False)
num_q, num_a = eigs.get_register_sizes()

# Initialize initial state module
init_state = Custom(num_q, state_vector=vector)

# Initialize reciprocal rotation module
reci = LookupRotation(negative_evals=eigs._negative_evals, evo_time=eigs._evo_time)

algo = HHL(matrix, vector, truncate_powerdim, truncate_hermitian, eigs,
           init_state, reci, num_q, num_a, orig_size)
/Users/hiroshi.wayama/anaconda3/lib/python3.8/site-packages/qiskit/aqua/components/initial_states/custom.py:79: DeprecationWarning: The Custom class is deprecated as of Aqua 0.9 and will be removed no earlier than 3 months after the release date. Instead, all algorithms and circuits accept a plain QuantumCircuit. Custom(state_vector=vector) is the same as a circuit where the ``initialize(vector/np.linalg.norm(vector))`` method has been called. method has been called.
  super(). __init__()
result = algo.run(QuantumInstance(Aer.get_backend('statevector_simulator')))
print("Solution:\t\t", np.round(result['solution'], 5))

result_ref = NumPyLSsolver(matrix, vector).run()
print("Classical Solution:\t", np.round(result_ref['solution'], 5))

print("Probability:\t\t %f" % result['probability_result'])
fidelity(result['solution'], result_ref['solution'])
/Users/hiroshi.wayama/anaconda3/lib/python3.8/site-packages/qiskit/aqua/components/initial_states/custom.py:151: DeprecationWarning The StateVectorCircuit class is deprecated as of Qiskit Aqua 0.9.0 and will be removed no earlier than 3 months after the release. Initialize a circuit, use the QuantumCircuit.initialize or QuantumCircuit.isometry methods. For a parameterized initialization, try the qiskit.ml RawFeatureVector class.
  svc = StateVectorCircuit(self._state_vector)

Solution: [ 0.66576-0.j -0.38561+0.j].
Classical Solution: [1.125 0.375].
Probability: 0.211527
Fidelity: 0.438807
print("circuit_width:\t", result['circuit_info']['width'])
print("circuit_depth:\t", result['circuit_info']['depth'])
print("CNOT gates:\t", result['circuit_info']['operations']['cx'])
circuit_width: 7
circuit_depth: 101
CNOT gates: 54

On the qiskit site, it looks like

Solution: [1.13586-0.j 0.40896+0.j].
Classical Solution: [1.125 0.375].

This time, the result is

Solution: [ 0.66576-0.j -0.38561+0.j].
Classical Solution: [1.125 0.375].

So, it seems that we are not getting good results.

As a note, I’ll leave you with the Custom and QuantumCircuit.initialize docs from the Depricate log.


init signature:
    num_qubits: int,
    state: str = 'zero',
    state_vector: Union[numpy.ndarray, qiskit.aqua.operators.state_fns.state_fn.StateFn, NoneType] = None,
    Circuit: Union[qiskit.circuit.quantumcircuit.QuantumCircuit, NoneType] = None,
) -> None

    Args: params (str or list or int):
        Args: params (str or list or int): * str: labels of basis states of the Pauli eigenstates Z, X, Y. See
            Args: params (str or list or int): * str: labels of basis states of the Pauli eigenstates Z, X, Y. See :meth:`~qiskit.quantum_info.states.statevector.Statevector.from_label`.
            Notice the order of the labels is reversed with respect to the qubit index to be applied to.
            Example label '01' initializes the qubit zero to `|1>` and the
            qubit one to `|0>`.
        * list: vector of complex amplitudes to initialize to.
        * int: an integer that is used as a bitmap indicating which qubits to initialize to `|1>`.
           * list: vector of complex amplitudes to initialize to. Example: setting params to 5 would initialize qubit 0 and qubit 2
           qubit 0 and qubit 2 would initialize qubit 0 to `|1>` and qubit 1 to `|0>`.
    qubits (QuantumRegister or int):
        * QuantumRegister: A list of qubits to be initialized [Default: None].
        * int: Index of qubits to be initialized [Default: None].

    qiskit.circuit.Instruction: a handle to the instruction that was just initialized.


Based on the HHL algorithm, a quantum recommendation algorithm and a classical algorithm inspired by it have been presented.

In the inspired classical algorithm, the User-Item matrix used in the recommendation system often has a low-rank approximation available (the huge number of user categories is much smaller than the number of users), which can be used to sample the solution state quickly. That’s the gist of it.

I will summarize the contents of this article in the near future.