diff --git a/docs/linear_solvers/hhl_tutorial.ipynb b/docs/linear_solvers/hhl_tutorial.ipynb
new file mode 100644
index 00000000..c4d87e1c
--- /dev/null
+++ b/docs/linear_solvers/hhl_tutorial.ipynb
@@ -0,0 +1,637 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "tags": [
+ "remove_cell"
+ ]
+ },
+ "source": [
+ "# Solving linear systems of equations using HHL and its Qiskit implementation\n",
+ "\n",
+ "The implementation is based on [https://dl.acm.org/doi/10.1145/3490631](https://dl.acm.org/doi/10.1145/3490631)[1](#enhancing)."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 1. Introduction \n",
+ "\n",
+ "Systems of linear equations arise naturally in many real-life applications in a wide range of areas, such as in the solution of Partial Differential Equations, the calibration of financial models, fluid simulation or numerical field calculation. The problem can be defined as, given a matrix $A\\in\\mathbb{C}^{N\\times N}$ and a vector $\\vec{b}\\in\\mathbb{C}^{N}$, find $\\vec{x}\\in\\mathbb{C}^{N}$ satisfying $A\\vec{x}=\\vec{b}$\n",
+ "\n",
+ "For example, take $N=2$, \n",
+ "\n",
+ "$$A = \\begin{pmatrix}1 & -1/3\\\\-1/3 & 1 \\end{pmatrix},\\quad \\vec{x}=\\begin{pmatrix} x_{1}\\\\ x_{2}\\end{pmatrix}\\quad \\text{and} \\quad \\vec{b}=\\begin{pmatrix}1 \\\\ 0\\end{pmatrix}$$\n",
+ "\n",
+ "Then the problem can also be written as find $x_{1}, x_{2}\\in\\mathbb{C}$ such that\n",
+ "$$\\begin{cases}x_{1} - \\frac{x_{2}}{3} = 1 \\\\ -\\frac{x_{1}}{3} + x_{2} = 0\\end{cases} $$\n",
+ "\n",
+ "A system of linear equations is called $s$-sparse if $A$ has at most $s$ non-zero entries per row or column. Solving an $s$-sparse system of size $N$ with a classical computer requires $\\mathcal{ O }(Ns\\kappa\\log(1/\\epsilon))$ running time using the conjugate gradient method [2](#conjgrad). Here $\\kappa$ denotes the condition number of the system and $\\epsilon$ the accuracy of the approximation.\n",
+ "\n",
+ "The HHL is a quantum algorithm to estimate a function of the solution with running time complexity of $\\mathcal{ O }(\\log(N)s^{2}\\kappa^{2}/\\epsilon)$[3](#hhl) when $A$ is a Hermitian matrix under the assumptions of efficient oracles for loading the data, Hamiltonian simulation and computing a function of the solution. This is an exponential speed up in the size of the system, however one crucial remark to keep in mind is that the classical algorithm returns the full solution, while the HHL can only approximate functions of the solution vector."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 2. The HHL algorithm\n",
+ "\n",
+ "### A. Some mathematical background\n",
+ "The first step towards solving a system of linear equations with a quantum computer is to encode the problem in the quantum language. By rescaling the system, we can assume $\\vec{b}$ and $\\vec{x}$ to be normalised and map them to the respective quantum states $|b\\rangle$ and $|x\\rangle$. Usually the mapping used is such that $i^{th}$ component of $\\vec{b}$ (resp. $\\vec{x}$) corresponds to the amplitude of the $i^{th}$ basis state of the quantum state $|b\\rangle$ (resp. $|x\\rangle$). From now on, we will focus on the rescaled problem\n",
+ "\n",
+ "$$ A|x\\rangle=|b\\rangle$$\n",
+ "\n",
+ "Since $A$ is Hermitian, it has a spectral decomposition\n",
+ "$$\n",
+ "A=\\sum_{j=0}^{N-1}\\lambda_{j}|u_{j}\\rangle\\langle u_{j}|,\\quad \\lambda_{j}\\in\\mathbb{ R }\n",
+ "$$\n",
+ "where $|u_{j}\\rangle$ is the $j^{th}$ eigenvector of $A$ with respective eigenvalue $\\lambda_{j}$. Then,\n",
+ "$$\n",
+ "A^{-1}=\\sum_{j=0}^{N-1}\\lambda_{j}^{-1}|u_{j}\\rangle\\langle u_{j}|\n",
+ "$$\n",
+ "and the right hand side of the system can be written in the eigenbasis of $A$ as\n",
+ "$$\n",
+ "|b\\rangle=\\sum_{j=0}^{N-1}b_{j}|u_{j}\\rangle,\\quad b_{j}\\in\\mathbb{ C }\n",
+ "$$\n",
+ "It is useful to keep in mind that the goal of the HHL is to exit the algorithm with the readout register in the state\n",
+ "$$\n",
+ "|x\\rangle=A^{-1}|b\\rangle=\\sum_{j=0}^{N-1}\\lambda_{j}^{-1}b_{j}|u_{j}\\rangle\n",
+ "$$\n",
+ "Note that here we already have an implicit normalisation constant since we are talking about a quantum state."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### B. Description of the HHL algorithm \n",
+ "\n",
+ "The algorithm uses three quantum registers, all of them set to $|0\\rangle $ at the beginning of the algorithm. One register, which we will denote with the subindex $n_{l}$, is used to store a binary representation of the eigenvalues of $A$. A second register, denoted by $n_{b}$, contains the vector solution, and from now on $N=2^{n_{b}}$. There is an extra register, for the auxiliary qubits. These are qubits used as intermediate steps in the individual computations but will be ignored in the following description since they are set to $|0\\rangle $ at the beginning of each computation and restored back to the $|0\\rangle $ state at the end of the individual operation.\n",
+ "\n",
+ "The following is an outline of the HHL algorithm with a high-level drawing of the corresponding circuit. For simplicity all computations are assumed to be exact in the ensuing description, and a more detailed explanation of the non-exact case is given in the [Qiskit Textbook](https://qiskit.org/textbook/ch-applications/hhl_tutorial.html).\n",
+ "\n",
+ "![](\"hhlcircuit.png\")
\n",
+ "\n",
+ "1. Load the data $|b\\rangle\\in\\mathbb{ C }^{N}$. That is, perform the transformation\n",
+ " $$ |0\\rangle _{n_{b}} \\mapsto |b\\rangle _{n_{b}} $$\n",
+ "2. Apply Quantum Phase Estimation (QPE) with\n",
+ "\t$$\n",
+ "\tU = e ^ { i A t } := \\sum _{j=0}^{N-1}e ^ { i \\lambda _ { j } t } |u_{j}\\rangle\\langle u_{j}|\n",
+ "\t$$\n",
+ "\tThe quantum state of the register expressed in the eigenbasis of $A$ is now\n",
+ "\t$$\n",
+ "\t\\sum_{j=0}^{N-1} b _ { j } |\\lambda _ {j }\\rangle_{n_{l}} |u_{j}\\rangle_{n_{b}}\n",
+ "\t$$\n",
+ " where $|\\lambda _ {j }\\rangle_{n_{l}}$ is the $n_{l}$-bit binary representation of $\\lambda _ {j }$.\n",
+ " \n",
+ "3. Add an auxiliary qubit and apply a rotation conditioned on $|\\lambda_{ j }\\rangle$,\n",
+ "\t$$\n",
+ "\t\\sum_{j=0}^{N-1} b _ { j } |\\lambda _ { j }\\rangle_{n_{l}}|u_{j}\\rangle_{n_{b}} \\left( \\sqrt { 1 - \\frac { C^{2} } { \\lambda _ { j } ^ { 2 } } } |0\\rangle + \\frac { C } { \\lambda _ { j } } |1\\rangle \\right)\n",
+ "\t$$\n",
+ "\twhere $C$ is a normalisation constant, and, as expressed in the current form above, should be less than the smallest eigenvalue $\\lambda_{min}$ in magnitude, i.e., $|C| < \\lambda_{min}$.\n",
+ " \n",
+ "4. Apply QPE$^{\\dagger}$. Ignoring possible errors from QPE, this results in\n",
+ "\t$$\n",
+ "\t\\sum_{j=0}^{N-1} b _ { j } |0\\rangle_{n_{l}}|u_{j}\\rangle_{n_{b}} \\left( \\sqrt { 1 - \\frac {C^{2} } { \\lambda _ { j } ^ { 2 } } } |0\\rangle + \\frac { C } { \\lambda _ { j } } |1\\rangle \\right)\n",
+ "\t$$\n",
+ " \n",
+ "5. Measure the auxiliary qubit in the computational basis. If the outcome is $1$, the register is in the post-measurement state\n",
+ "\t$$\n",
+ "\t\\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}}|u_{j}\\rangle_{n_{b}}\n",
+ "\t$$\n",
+ "\twhich up to a normalisation factor corresponds to the solution.\n",
+ "\n",
+ "6. Apply an observable $M$ to calculate $F(x):=\\langle x|M|x\\rangle$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 4. Qiskit Implementation"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now that we have analytically solved the problem from the example we are going to use it to illustrate how to run the HHL on a quantum simulator and on the real hardware. For the quantum simulator, Qiskit already provides an implementation of the HHL algorithm requiring only the matrix $A$ and $|b\\rangle$ as inputs in the simplest example. Although we can give the algorithm a general Hermitian matrix and an arbitrary initial state as NumPy arrays, in these cases the quantum algorithm will not achieve an exponential speedup. This is because the default implementation is exact and therefore exponential in the number of qubits (there is no algorithm that can prepare exactly an arbitrary quantum state using polynomial resources in the number of qubits or that can perform exactly the operation $e^{iAt}$ for some general Hermitian matrix $A$ using polynomial resources in the number of qubits). If we know an efficient implementation for a particular problem, the matrix and/or the vector can be given as `QuantumCircuit` objects. Alternatively, there's already an efficient implementation for tridiagonal Toeplitz matrices and in the future there might be more."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The interface for all algorithms to solve the linear system problem is `LinearSolver`. The problem to be solved is only specified when the `solve()` method is called:\n",
+ "```python\n",
+ "LinearSolver(...).solve(matrix, vector)\n",
+ "```\n",
+ "\n",
+ "The simplest implementation takes the matrix and the vector as NumPy arrays. Below we also create a `NumPyLinearSolver` (the classical algorithm) to validate our solutions."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "from qiskit_research.linear_solvers.numpy_linear_solver import NumPyLinearSolver\n",
+ "from qiskit_research.linear_solvers.hhl import HHL\n",
+ "\n",
+ "matrix = np.array([[1, -1 / 3], [-1 / 3, 1]])\n",
+ "vector = np.array([1, 0])\n",
+ "naive_hhl_solution = HHL().solve(matrix, vector)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "For the classical solver we need to rescale the right hand side (i.e. `vector / np.linalg.norm(vector)`) to take into account the renormalisation that occurs once `vector` is encoded in a quantum state within HHL."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "classical_solution = NumPyLinearSolver().solve(matrix, vector / np.linalg.norm(vector))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The `linear_solvers` package contains a folder called `matrices` intended to be a placeholder for efficient implementations of particular types of matrices. At the time of writing the only truly efficient implementation it contains (i.e. complexity scaling polynomially in the number of qubits) is the `TridiagonalToeplitz` class. Tridiagonal Toeplitz symmetric real matrices are of the following form \n",
+ "$$A = \\begin{pmatrix}a & b & 0 & 0\\\\b & a & b & 0 \\\\ 0 & b & a & b \\\\ 0 & 0 & b & a \\end{pmatrix}, a,b\\in\\mathbb{R}$$\n",
+ "(note that in this setting we do not consider non symmetric matrices since the HHL algorithm assumes that the input matrix is Hermitian).\n",
+ "\n",
+ "Since the matrix $A$ from our example is of this form we can create an instance of `TridiagonalToeplitz(num_qubits, a, b)` and compare the results to solving the system with an array as input."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from qiskit_research.linear_solvers.matrices.tridiagonal_toeplitz import (\n",
+ " TridiagonalToeplitz,\n",
+ ")\n",
+ "\n",
+ "tridi_matrix = TridiagonalToeplitz(1, 1, -1 / 3)\n",
+ "\n",
+ "tridi_solution = HHL().solve(tridi_matrix, vector)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Recall that the HHL algorithm can find a solution exponentially faster in the size of the system than their classical counterparts (i.e. logarithmic complexity instead of polynomial). However the cost for this exponential speedup is that we do not obtain the full solution vector.\n",
+ "Instead, we obtain a quantum state representing the vector $x$ and learning all the components of this vector would take a linear time in its dimension, diminishing any speedup obtained by the quantum algorithm.\n",
+ "\n",
+ "Therefore, we can only compute functions from $x$ (the so called observables) to learn information about the solution.\n",
+ "This is reflected in the `LinearSolverResult` object returned by `solve()`, which contains the following properties\n",
+ "- `state`: either the circuit that prepares the solution or the solution as a vector\n",
+ "- `euclidean_norm`: the euclidean norm if the algorithm knows how to calculate it \n",
+ "- `observable`: the (list of) calculated observable(s)\n",
+ "- `circuit_results`: the observable results from the (list of) circuit(s)\n",
+ "\n",
+ "Let's ignore `observable` and `circuit_results` for the time being and check the solutions we obtained before.\n",
+ "\n",
+ "First, `classical_solution` was the result from a classical algorithm, so if we call `.state` it will return an array:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "classical state: [1.125 0.375]\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(\"classical state:\", classical_solution.state)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Our other two examples were quantum algorithms, hence we can only access to the quantum state. This is achieved by returning the quantum circuit that prepares the solution state:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "naive state:\n",
+ " ┌────────────┐┌──────┐ ┌─────────┐\n",
+ " q4: ┤ circuit-85 ├┤3 ├────────┤3 ├\n",
+ " └────────────┘│ │┌──────┐│ │\n",
+ "q5_0: ──────────────┤0 ├┤2 ├┤0 ├\n",
+ " │ QPE ││ ││ QPE_dg │\n",
+ "q5_1: ──────────────┤1 ├┤1 ├┤1 ├\n",
+ " │ ││ 1/x ││ │\n",
+ "q5_2: ──────────────┤2 ├┤0 ├┤2 ├\n",
+ " └──────┘│ │└─────────┘\n",
+ " q6: ──────────────────────┤3 ├───────────\n",
+ " └──────┘ \n",
+ "tridiagonal state:\n",
+ " ┌─────────────┐┌──────┐ ┌─────────┐\n",
+ " q26: ┤ circuit-298 ├┤3 ├────────┤3 ├\n",
+ " └─────────────┘│ │┌──────┐│ │\n",
+ "q27_0: ───────────────┤0 ├┤2 ├┤0 ├\n",
+ " │ QPE ││ ││ QPE_dg │\n",
+ "q27_1: ───────────────┤1 ├┤1 ├┤1 ├\n",
+ " │ ││ 1/x ││ │\n",
+ "q27_2: ───────────────┤2 ├┤0 ├┤2 ├\n",
+ " └──────┘│ │└─────────┘\n",
+ " q28: ───────────────────────┤3 ├───────────\n",
+ " └──────┘ \n"
+ ]
+ }
+ ],
+ "source": [
+ "print(\"naive state:\")\n",
+ "print(naive_hhl_solution.state)\n",
+ "print(\"tridiagonal state:\")\n",
+ "print(tridi_solution.state)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Recall that the Euclidean norm for a vector $\\mathbf{x}=(x_1,\\dots,x_N)$ is defined as $||\\mathbf{x}||=\\sqrt{\\sum_{i=1}^N x_i^2}$. Therefore, the probability of measuring $1$ in the auxiliary qubit from Step $5$ in Section B is the squared norm of $\\mathbf{x}$. This means that the HHL algorithm can always calculate the euclidean norm of the solution and we can compare the accuracy of the results:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "classical Euclidean norm: 1.1858541225631423\n",
+ "naive Euclidean norm: 1.1858541225631407\n",
+ "tridiagonal Euclidean norm: 1.185854122563139\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(\"classical Euclidean norm:\", classical_solution.euclidean_norm)\n",
+ "print(\"naive Euclidean norm:\", naive_hhl_solution.euclidean_norm)\n",
+ "print(\"tridiagonal Euclidean norm:\", tridi_solution.euclidean_norm)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Comparing the solution vectors componentwise is more tricky, reflecting again the idea that we cannot obtain the full solution vector from the quantum algorithm. However, for educational purposes we can check that indeed the different solution vectors obtained are a good approximation at the vector component level as well. \n",
+ "\n",
+ "To do so first we need to use `Statevector` from the `quantum_info` package and extract the right vector components, i.e. those corresponding to the ancillary qubit (bottom in the circuits) being $1$ and the work qubits (the two middle in the circuits) being $0$. Thus, we are interested in the states `1000` and `1001`, corresponding to the first and second components of the solution vector respectively."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "naive raw solution vector: [ 5.74337504e-17+8.93013426e-18j -1.12425450e-17+5.23599292e-16j]\n",
+ "tridi raw solution vector: [-1.92535646e-16+9.81307787e-17j -2.24504108e-17+6.80168960e-16j]\n"
+ ]
+ }
+ ],
+ "source": [
+ "from qiskit.quantum_info import Statevector\n",
+ "\n",
+ "naive_sv = Statevector(naive_hhl_solution.state).data\n",
+ "tridi_sv = Statevector(tridi_solution.state).data\n",
+ "\n",
+ "# Extract the right vector components. 1000 corresponds to the index 8 and 1001 corresponds to the index 9\n",
+ "naive_full_vector = np.array([naive_sv[8], naive_sv[9]])\n",
+ "tridi_full_vector = np.array([tridi_sv[8], tridi_sv[9]])\n",
+ "\n",
+ "print(\"naive raw solution vector:\", naive_full_vector)\n",
+ "print(\"tridi raw solution vector:\", tridi_full_vector)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "At a first glance it might seem that this is wrong because the components are complex numbers instead of reals. However note that the imaginary part is very small, most likely due to computer accuracy, and can be disregarded in this case."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "naive_full_vector = np.real(naive_full_vector)\n",
+ "tridi_full_vector = np.real(tridi_full_vector)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Next, we will divide the vectors by their respective norms to suppress any constants coming from the different parts of the circuits. The full solution vector can then be recovered by multiplying these normalised vectors by the respective Euclidean norms calculated above:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "full naive solution vector: [ 1.16376749 -0.22780523]\n",
+ "full tridi solution vector: [-1.17787369 -0.13734469]\n",
+ "classical state: [1.125 0.375]\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(\n",
+ " \"full naive solution vector:\",\n",
+ " naive_hhl_solution.euclidean_norm\n",
+ " * naive_full_vector\n",
+ " / np.linalg.norm(naive_full_vector),\n",
+ ")\n",
+ "print(\n",
+ " \"full tridi solution vector:\",\n",
+ " tridi_solution.euclidean_norm\n",
+ " * tridi_full_vector\n",
+ " / np.linalg.norm(tridi_full_vector),\n",
+ ")\n",
+ "print(\"classical state:\", classical_solution.state)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "It should not come as a surprise that `naive_hhl_solution` is exact because all the default methods used are exact. However, `tridi_solution` is exact only in the $2\\times 2$ system size case. For larger matrices it will be an approximation, as shown in the slightly larger example below."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "classical euclidean norm: 1.237833351044751\n",
+ "naive euclidean norm: 1.2099806231119126\n",
+ "tridiagonal euclidean norm: 1.209457721870593\n"
+ ]
+ }
+ ],
+ "source": [
+ "from scipy.sparse import diags\n",
+ "\n",
+ "num_qubits = 2\n",
+ "matrix_size = 2**num_qubits\n",
+ "# entries of the tridiagonal Toeplitz symmetric matrix\n",
+ "a = 1\n",
+ "b = -1 / 3\n",
+ "\n",
+ "matrix = diags([b, a, b], [-1, 0, 1], shape=(matrix_size, matrix_size)).toarray()\n",
+ "vector = np.array([1] + [0] * (matrix_size - 1))\n",
+ "\n",
+ "# run the algorithms\n",
+ "classical_solution = NumPyLinearSolver().solve(matrix, vector / np.linalg.norm(vector))\n",
+ "naive_hhl_solution = HHL().solve(matrix, vector)\n",
+ "tridi_matrix = TridiagonalToeplitz(num_qubits, a, b)\n",
+ "tridi_solution = HHL().solve(tridi_matrix, vector)\n",
+ "\n",
+ "print(\"classical euclidean norm:\", classical_solution.euclidean_norm)\n",
+ "print(\"naive euclidean norm:\", naive_hhl_solution.euclidean_norm)\n",
+ "print(\"tridiagonal euclidean norm:\", tridi_solution.euclidean_norm)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now we can return to the topic of observables and find out what the `observable` and `circuit_results` properties contain.\n",
+ "\n",
+ "The way to compute functions of the solution vector $\\mathbf{x}$ is through giving the `.solve()` method a `LinearSystemObservable` as input. There are are two types of available `LinearSystemObservable` which can be given as input:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from qiskit_research.linear_solvers.observables import AbsoluteAverage, MatrixFunctional"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "For a vector $\\mathbf{x}=(x_1,...,x_N)$, the `AbsoluteAverage` observable computes $|\\frac{1}{N}\\sum_{i=1}^{N}x_i|$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "quantum average: 0.7499999999999976\n",
+ "classical average: 0.75\n",
+ "quantum circuit results: (0.49999999999999706+0j)\n"
+ ]
+ }
+ ],
+ "source": [
+ "num_qubits = 1\n",
+ "matrix_size = 2**num_qubits\n",
+ "# entries of the tridiagonal Toeplitz symmetric matrix\n",
+ "a = 1\n",
+ "b = -1 / 3\n",
+ "\n",
+ "matrix = diags([b, a, b], [-1, 0, 1], shape=(matrix_size, matrix_size)).toarray()\n",
+ "vector = np.array([1] + [0] * (matrix_size - 1))\n",
+ "tridi_matrix = TridiagonalToeplitz(1, a, b)\n",
+ "\n",
+ "average_solution = HHL().solve(tridi_matrix, vector, AbsoluteAverage())\n",
+ "classical_average = NumPyLinearSolver().solve(\n",
+ " matrix, vector / np.linalg.norm(vector), AbsoluteAverage()\n",
+ ")\n",
+ "\n",
+ "print(\"quantum average:\", average_solution.observable)\n",
+ "print(\"classical average:\", classical_average.observable)\n",
+ "print(\"quantum circuit results:\", average_solution.circuit_results)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The `MatrixFunctional` observable computes $\\mathbf{x}^T B \\mathbf{x}$ for a vector $\\mathbf{x}$ and a tridiagonal symmetric Toeplitz matrix $B$. The class takes the main and off diagonal values of the matrix for its constuctor method."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "quantum functional: 1.8281249999999882\n",
+ "classical functional: 1.828125\n",
+ "quantum circuit results: [(0.6249999999999963+0j), (0.49999999999999706+0j), (0.12499999999999926+0j)]\n"
+ ]
+ }
+ ],
+ "source": [
+ "observable = MatrixFunctional(1, 1 / 2)\n",
+ "\n",
+ "functional_solution = HHL().solve(tridi_matrix, vector, observable)\n",
+ "classical_functional = NumPyLinearSolver().solve(\n",
+ " matrix, vector / np.linalg.norm(vector), observable\n",
+ ")\n",
+ "\n",
+ "print(\"quantum functional:\", functional_solution.observable)\n",
+ "print(\"classical functional:\", classical_functional.observable)\n",
+ "print(\"quantum circuit results:\", functional_solution.circuit_results)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Therefore, `observable` contains the final value of the function on $\\mathbf{x}$, while `circuit_results` contains the raw values obtained from the circuit and used to process the result of `observable`.\n",
+ "\n",
+ "This 'how to process the result' is better explained by looking at what arguments `.solve()` takes. The `solve()` method accepts up to five arguments: \n",
+ "```python\n",
+ "def solve(self, matrix: Union[np.ndarray, QuantumCircuit],\n",
+ " vector: Union[np.ndarray, QuantumCircuit],\n",
+ " observable: Optional[Union[LinearSystemObservable, BaseOperator,\n",
+ " List[BaseOperator]]] = None,\n",
+ " post_rotation: Optional[Union[QuantumCircuit, List[QuantumCircuit]]] = None,\n",
+ " post_processing: Optional[Callable[[Union[float, List[float]]],\n",
+ " Union[float, List[float]]]] = None) \\\n",
+ " -> LinearSolverResult:\n",
+ "```\n",
+ "The first two are the matrix defining the linear system and the vector right hand side of the equation, which we have already covered. The remaining parameters concern the (list of) observable(s) to be computed out of the solution vector $x$, and can be specified in two different ways. One option is to give as the third and last parameter a (list of) `LinearSystemObservable`(s). Alternatively, we can give our own implementations of the `observable`, `post_rotation` and `post_processing`, where\n",
+ "- `observable` is the operator to compute the expected value of the observable and can be e.g. a `PauliSumOp`\n",
+ "- `post_rotation` is the circuit to be applied to the solution to extract information if additional gates are needed.\n",
+ "- `post_processing` is the function to compute the value of the observable from the calculated probabilities.\n",
+ "\n",
+ "In other words, there will be as many `circuit_results` as `post_rotation` circuits, and `post_processing` is telling the algorithm how to use the values we see when we print `circuit_results` to obtain the value we see when we print `observable`.\n",
+ "\n",
+ "Finally, the `HHL` class accepts the following parameters in its constructor method:\n",
+ "- error tolerance : the accuracy of the approximation of the solution, the default is `1e-2`\n",
+ "- expectation : how the expectation values are evaluated, the default is `PauliExpectation`\n",
+ "- quantum instance: the `QuantumInstance` or backend, the default is a `Statevector` simulation"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "1.1858541225631407\n",
+ "1.1858541225631423\n"
+ ]
+ }
+ ],
+ "source": [
+ "from qiskit import BasicAer\n",
+ "\n",
+ "backend = BasicAer.get_backend(\"qasm_simulator\")\n",
+ "hhl = HHL(1e-3, quantum_instance=backend)\n",
+ "\n",
+ "accurate_solution = hhl.solve(matrix, vector)\n",
+ "classical_solution = NumPyLinearSolver().solve(matrix, vector / np.linalg.norm(vector))\n",
+ "\n",
+ "print(accurate_solution.euclidean_norm)\n",
+ "print(classical_solution.euclidean_norm)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 9. References"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "1. A. Carrera Vazquez, A. Frisch, D. Steenken, H. S. Barowski, R. Hiptmair, and S. Woerner, “Enhancing Quantum Linear System Algorithm by Richardson Extrapolation,” ACM Trans. Quantum Comput. 3 (2022).\n",
+ "2. J. R. Shewchuk. An Introduction to the Conjugate Gradient Method Without the Agonizing Pain. Technical Report CMU-CS-94-125, School of Computer Science, Carnegie Mellon University, Pittsburgh, Pennsylvania, March 1994. \n",
+ "3. A. W. Harrow, A. Hassidim, and S. Lloyd, “Quantum algorithm for linear systems of equations,” Phys. Rev. Lett. 103.15 (2009), p. 150502."
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "contribute",
+ "language": "python",
+ "name": "contribute"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.8.12"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/linear_solvers/hhlcircuit.png b/docs/linear_solvers/hhlcircuit.png
new file mode 100644
index 00000000..fc17727a
Binary files /dev/null and b/docs/linear_solvers/hhlcircuit.png differ
diff --git a/qiskit_research/linear_solvers/__init__.py b/qiskit_research/linear_solvers/__init__.py
new file mode 100644
index 00000000..97cc43a5
--- /dev/null
+++ b/qiskit_research/linear_solvers/__init__.py
@@ -0,0 +1,92 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""
+Linear solvers (:mod:`qiskit_research.linear_solvers`)
+=========================================================
+It contains classical and quantum algorithms to solve systems of linear equations such as
+:class:`~qiskit_research.linear_solvers.HHL`.
+Although the quantum algorithm accepts a general Hermitian matrix as input, Qiskit's default
+Hamiltonian evolution is exponential in such cases and therefore the quantum linear solver will
+not achieve an exponential speedup.
+Furthermore, the quantum algorithm can find a solution exponentially faster in the size of the
+system than their classical counterparts (i.e. logarithmic complexity instead of polynomial),
+meaning that reading the full solution vector would kill such speedup (since this would take
+linear time in the size of the system).
+Therefore, to achieve an exponential speedup we can only compute functions from the solution
+vector (the so called observables) to learn information about the solution.
+Known efficient implementations of Hamiltonian evolutions or observables are contained in the
+following subfolders:
+
+`Matrices`_
+ A placeholder for efficient implementations of the Hamiltonian evolution of particular types of
+ matrices.
+
+`Observables`_
+ A placeholder for efficient implementations of functions that can be computed from the solution
+ vector to a system of linear equations.
+
+.. currentmodule:: qiskit_research.linear_solvers
+
+Linear Solvers
+==============
+
+.. autosummary::
+ :toctree: ../stubs/
+ :nosignatures:
+
+ LinearSolver
+ LinearSolverResult
+ HHL
+ NumPyLinearSolver
+
+Matrices
+========
+
+.. autosummary::
+ :toctree: ../stubs/
+ :nosignatures:
+
+ LinearSystemMatrix
+ NumPyMatrix
+ TridiagonalToeplitz
+
+Observables
+===========
+
+.. autosummary::
+ :toctree: ../stubs/
+
+ LinearSystemObservable
+ AbsoluteAverage
+ MatrixFunctional
+
+"""
+
+from .hhl import HHL
+from .numpy_linear_solver import NumPyLinearSolver
+from .linear_solver import LinearSolver, LinearSolverResult
+from .matrices import LinearSystemMatrix, NumPyMatrix, TridiagonalToeplitz
+from .observables import LinearSystemObservable, AbsoluteAverage, MatrixFunctional
+
+__all__ = [
+ "HHL",
+ "NumPyLinearSolver",
+ "LinearSolver",
+ "LinearSolverResult",
+ "LinearSystemMatrix",
+ "NumPyMatrix",
+ "TridiagonalToeplitz",
+ "LinearSystemObservable",
+ "AbsoluteAverage",
+ "MatrixFunctional",
+]
diff --git a/qiskit_research/linear_solvers/hhl.py b/qiskit_research/linear_solvers/hhl.py
new file mode 100644
index 00000000..4ffe9ebd
--- /dev/null
+++ b/qiskit_research/linear_solvers/hhl.py
@@ -0,0 +1,559 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""The HHL algorithm."""
+
+from typing import Optional, Union, List, Callable, Tuple
+import numpy as np
+
+from qiskit.circuit import QuantumCircuit, QuantumRegister, AncillaRegister
+from qiskit.circuit.library import PhaseEstimation
+from qiskit.circuit.library.arithmetic.piecewise_chebyshev import PiecewiseChebyshev
+from qiskit.circuit.library.arithmetic.exact_reciprocal import ExactReciprocal
+from qiskit.opflow import (
+ Z,
+ I,
+ StateFn,
+ TensoredOp,
+ ExpectationBase,
+ CircuitSampler,
+ ListOp,
+ ExpectationFactory,
+ ComposedOp,
+)
+from qiskit.providers import Backend
+from qiskit.utils import QuantumInstance
+
+from .linear_solver import LinearSolver, LinearSolverResult
+from .matrices.numpy_matrix import NumPyMatrix
+from .observables.linear_system_observable import LinearSystemObservable
+
+
+class HHL(LinearSolver):
+ r"""Systems of linear equations arise naturally in many real-life applications in a wide range
+ of areas, such as in the solution of Partial Differential Equations, the calibration of
+ financial models, fluid simulation or numerical field calculation. The problem can be defined
+ as, given a matrix :math:`A\in\mathbb{C}^{N\times N}` and a vector
+ :math:`\vec{b}\in\mathbb{C}^{N}`, find :math:`\vec{x}\in\mathbb{C}^{N}` satisfying
+ :math:`A\vec{x}=\vec{b}`.
+
+ A system of linear equations is called :math:`s`-sparse if :math:`A` has at most :math:`s`
+ non-zero entries per row or column. Solving an :math:`s`-sparse system of size :math:`N` with
+ a classical computer requires :math:`\mathcal{ O }(Ns\kappa\log(1/\epsilon))` running time
+ using the conjugate gradient method. Here :math:`\kappa` denotes the condition number of the
+ system and :math:`\epsilon` the accuracy of the approximation.
+
+ The HHL is a quantum algorithm to estimate a function of the solution with running time
+ complexity of :math:`\mathcal{ O }(\log(N)s^{2}\kappa^{2}/\epsilon)` when
+ :math:`A` is a Hermitian matrix under the assumptions of efficient oracles for loading the
+ data, Hamiltonian simulation and computing a function of the solution. This is an exponential
+ speed up in the size of the system, however one crucial remark to keep in mind is that the
+ classical algorithm returns the full solution, while the HHL can only approximate functions of
+ the solution vector.
+
+ Examples:
+
+ .. jupyter-execute::
+
+ import numpy as np
+ from qiskit import QuantumCircuit
+ from qiskit_research.linear_solvers.hhl import HHL
+ from qiskit_research.linear_solvers.matrices import TridiagonalToeplitz
+ from qiskit_research.linear_solvers.observables import MatrixFunctional
+
+ matrix = TridiagonalToeplitz(2, 1, 1 / 3, trotter_steps=2)
+ right_hand_side = [1.0, -2.1, 3.2, -4.3]
+ observable = MatrixFunctional(1, 1 / 2)
+ rhs = right_hand_side / np.linalg.norm(right_hand_side)
+
+ # Initial state circuit
+ num_qubits = matrix.num_state_qubits
+ qc = QuantumCircuit(num_qubits)
+ qc.isometry(rhs, list(range(num_qubits)), None)
+
+ hhl = HHL()
+ solution = hhl.solve(matrix, qc, observable)
+ approx_result = solution.observable
+
+ References:
+
+ [1]: Harrow, A. W., Hassidim, A., Lloyd, S. (2009).
+ Quantum algorithm for linear systems of equations.
+ `Phys. Rev. Lett. 103, 15 (2009), 1–15. `_
+
+ [2]: Carrera Vazquez, A., Hiptmair, R., & Woerner, S. (2020).
+ Enhancing the Quantum Linear Systems Algorithm using Richardson Extrapolation.
+ `arXiv:2009.04484 `_
+
+ """
+
+ def __init__(
+ self,
+ epsilon: float = 1e-2,
+ expectation: Optional[ExpectationBase] = None,
+ quantum_instance: Optional[Union[Backend, QuantumInstance]] = None,
+ ) -> None:
+ r"""
+ Args:
+ epsilon: Error tolerance of the approximation to the solution, i.e. if :math:`x` is the
+ exact solution and :math:`\tilde{x}` the one calculated by the algorithm, then
+ :math:`||x - \tilde{x}|| \le epsilon`.
+ expectation: The expectation converter applied to the expectation values before
+ evaluation. If None then PauliExpectation is used.
+ quantum_instance: Quantum Instance or Backend. If None, a Statevector calculation is
+ done.
+ """
+ super().__init__()
+
+ self._epsilon = epsilon
+ # Tolerance for the different parts of the algorithm as per [1]
+ self._epsilon_r = epsilon / 3 # conditioned rotation
+ self._epsilon_s = epsilon / 3 # state preparation
+ self._epsilon_a = epsilon / 6 # hamiltonian simulation
+
+ self._scaling = None # scaling of the solution
+
+ self._sampler = None
+ self.quantum_instance = quantum_instance
+
+ self._expectation = expectation
+
+ # For now the default reciprocal implementation is exact
+ self._exact_reciprocal = True
+ # Set the default scaling to 1
+ self.scaling = 1
+
+ @property
+ def quantum_instance(self) -> Optional[QuantumInstance]:
+ """Get the quantum instance.
+
+ Returns:
+ The quantum instance used to run this algorithm.
+ """
+ return None if self._sampler is None else self._sampler.quantum_instance
+
+ @quantum_instance.setter
+ def quantum_instance(
+ self, quantum_instance: Optional[Union[QuantumInstance, Backend]]
+ ) -> None:
+ """Set quantum instance.
+
+ Args:
+ quantum_instance: The quantum instance used to run this algorithm.
+ If None, a Statevector calculation is done.
+ """
+ if quantum_instance is not None:
+ self._sampler = CircuitSampler(quantum_instance)
+ else:
+ self._sampler = None
+
+ @property
+ def scaling(self) -> float:
+ """The scaling of the solution vector."""
+ return self._scaling
+
+ @scaling.setter
+ def scaling(self, scaling: float) -> None:
+ """Set the new scaling of the solution vector."""
+ self._scaling = scaling
+
+ @property
+ def expectation(self) -> ExpectationBase:
+ """The expectation value algorithm used to construct the expectation measurement from
+ the observable."""
+ return self._expectation
+
+ @expectation.setter
+ def expectation(self, expectation: ExpectationBase) -> None:
+ """Set the expectation value algorithm."""
+ self._expectation = expectation
+
+ def _get_delta(self, n_l: int, lambda_min: float, lambda_max: float) -> float:
+ """Calculates the scaling factor to represent exactly lambda_min on nl binary digits.
+
+ Args:
+ n_l: The number of qubits to represent the eigenvalues.
+ lambda_min: the smallest eigenvalue.
+ lambda_max: the largest eigenvalue.
+
+ Returns:
+ The value of the scaling factor.
+ """
+ formatstr = "#0" + str(n_l + 2) + "b"
+ lambda_min_tilde = np.abs(lambda_min * (2**n_l - 1) / lambda_max)
+ # floating point precision can cause problems
+ if np.abs(lambda_min_tilde - 1) < 1e-7:
+ lambda_min_tilde = 1
+ binstr = format(int(lambda_min_tilde), formatstr)[2::]
+ lamb_min_rep = 0
+ for i, char in enumerate(binstr):
+ lamb_min_rep += int(char) / (2 ** (i + 1))
+ return lamb_min_rep
+
+ def _calculate_norm(self, qc: QuantumCircuit) -> float:
+ """Calculates the value of the euclidean norm of the solution.
+
+ Args:
+ qc: The quantum circuit preparing the solution x to the system.
+
+ Returns:
+ The value of the euclidean norm of the solution.
+ """
+ # Calculate the number of qubits
+ nb = qc.qregs[0].size
+ nl = qc.qregs[1].size
+ na = qc.num_ancillas
+
+ # Create the Operators Zero and One
+ zero_op = (I + Z) / 2
+ one_op = (I - Z) / 2
+
+ # Norm observable
+ observable = one_op ^ TensoredOp((nl + na) * [zero_op]) ^ (I ^ nb)
+ norm_2 = (~StateFn(observable) @ StateFn(qc)).eval()
+
+ return np.real(np.sqrt(norm_2) / self.scaling)
+
+ def _calculate_observable(
+ self,
+ solution: QuantumCircuit,
+ ls_observable: Optional[LinearSystemObservable] = None,
+ observable_circuit: Optional[QuantumCircuit] = None,
+ post_processing: Optional[
+ Callable[[Union[float, List[float]], int, float], float]
+ ] = None,
+ ) -> Tuple[float, Union[complex, List[complex]]]:
+ """Calculates the value of the observable(s) given.
+
+ Args:
+ solution: The quantum circuit preparing the solution x to the system.
+ ls_observable: Information to be extracted from the solution.
+ observable_circuit: Circuit to be applied to the solution to extract information.
+ post_processing: Function to compute the value of the observable.
+
+ Returns:
+ The value of the observable(s) and the circuit results before post-processing as a
+ tuple.
+ """
+ # Get the number of qubits
+ nb = solution.qregs[0].size
+ nl = solution.qregs[1].size
+ na = solution.num_ancillas
+
+ # if the observable is given construct post_processing and observable_circuit
+ if ls_observable is not None:
+ observable_circuit = ls_observable.observable_circuit(nb)
+ post_processing = ls_observable.post_processing
+
+ if isinstance(ls_observable, LinearSystemObservable):
+ observable = ls_observable.observable(nb)
+
+ # in the other case use the identity as observable
+ else:
+ observable = I ^ nb
+
+ # Create the Operators Zero and One
+ zero_op = (I + Z) / 2
+ one_op = (I - Z) / 2
+
+ is_list = True
+ if not isinstance(observable_circuit, list):
+ is_list = False
+ observable_circuit = [observable_circuit]
+ observable = [observable]
+
+ expectations: Union[ListOp, ComposedOp] = []
+ for circ, obs in zip(observable_circuit, observable):
+ circuit = QuantumCircuit(solution.num_qubits)
+ circuit.append(solution, circuit.qubits)
+ circuit.append(circ, range(nb))
+
+ ob = one_op ^ TensoredOp((nl + na) * [zero_op]) ^ obs
+ expectations.append(~StateFn(ob) @ StateFn(circuit))
+
+ if is_list:
+ # execute all in a list op to send circuits in batches
+ expectations = ListOp(expectations)
+ else:
+ expectations = expectations[0]
+
+ # check if an expectation converter is given
+ if self._expectation is not None:
+ expectations = self._expectation.convert(expectations)
+ # if otherwise a backend was specified, try to set the best expectation value
+ elif self._sampler is not None:
+ if is_list:
+ op = expectations.oplist[0]
+ else:
+ op = expectations
+ self._expectation = ExpectationFactory.build(
+ op, self._sampler.quantum_instance
+ )
+
+ if self._sampler is not None:
+ expectations = self._sampler.convert(expectations)
+
+ # evaluate
+ expectation_results = expectations.eval()
+
+ # apply post_processing
+ result = post_processing(expectation_results, nb, self.scaling)
+
+ return result, expectation_results
+
+ def construct_circuit(
+ self,
+ matrix: Union[List, np.ndarray, QuantumCircuit],
+ vector: Union[List, np.ndarray, QuantumCircuit],
+ neg_vals: Optional[bool] = True,
+ ) -> QuantumCircuit:
+ """Construct the HHL circuit.
+
+ Args:
+ matrix: The matrix specifying the system, i.e. A in Ax=b.
+ vector: The vector specifying the right hand side of the equation in Ax=b.
+ neg_vals: States whether the matrix has negative eigenvalues. If False the
+ computation becomes cheaper.
+
+ Returns:
+ The HHL circuit.
+
+ Raises:
+ ValueError: If the input is not in the correct format.
+ ValueError: If the type of the input matrix is not supported.
+ """
+ # State preparation circuit - default is qiskit
+ if isinstance(vector, QuantumCircuit):
+ nb = vector.num_qubits
+ vector_circuit = vector
+ elif isinstance(vector, (list, np.ndarray)):
+ if isinstance(vector, list):
+ vector = np.array(vector)
+ nb = int(np.log2(len(vector)))
+ vector_circuit = QuantumCircuit(nb)
+ # pylint: disable=no-member
+ vector_circuit.isometry(
+ vector / np.linalg.norm(vector), list(range(nb)), None
+ )
+
+ # If state preparation is probabilistic the number of qubit flags should increase
+ nf = 1
+
+ # Hamiltonian simulation circuit - default is Trotterization
+ if isinstance(matrix, QuantumCircuit):
+ matrix_circuit = matrix
+ elif isinstance(matrix, (list, np.ndarray)):
+ if isinstance(matrix, list):
+ matrix = np.array(matrix)
+
+ if matrix.shape[0] != matrix.shape[1]:
+ raise ValueError("Input matrix must be square!")
+ if np.log2(matrix.shape[0]) % 1 != 0:
+ raise ValueError("Input matrix dimension must be 2^n!")
+ if not np.allclose(matrix, matrix.conj().T):
+ raise ValueError("Input matrix must be hermitian!")
+ if matrix.shape[0] != 2**vector_circuit.num_qubits:
+ raise ValueError(
+ "Input vector dimension does not match input "
+ "matrix dimension! Vector dimension: "
+ + str(vector_circuit.num_qubits)
+ + ". Matrix dimension: "
+ + str(matrix.shape[0])
+ )
+ matrix_circuit = NumPyMatrix(matrix, evolution_time=2 * np.pi)
+ else:
+ raise ValueError(f"Invalid type for matrix: {type(matrix)}.")
+
+ # Set the tolerance for the matrix approximation
+ if hasattr(matrix_circuit, "tolerance"):
+ matrix_circuit.tolerance = self._epsilon_a
+
+ # check if the matrix can calculate the condition number and store the upper bound
+ if (
+ hasattr(matrix_circuit, "condition_bounds")
+ and matrix_circuit.condition_bounds() is not None
+ ):
+ kappa = matrix_circuit.condition_bounds()[1]
+ else:
+ kappa = 1
+ # Update the number of qubits required to represent the eigenvalues
+ # The +neg_vals is to register negative eigenvalues because
+ # e^{-2 \pi i \lambda} = e^{2 \pi i (1 - \lambda)}
+ nl = max(nb + 1, int(np.ceil(np.log2(kappa + 1)))) + neg_vals
+
+ # check if the matrix can calculate bounds for the eigenvalues
+ if (
+ hasattr(matrix_circuit, "eigs_bounds")
+ and matrix_circuit.eigs_bounds() is not None
+ ):
+ lambda_min, lambda_max = matrix_circuit.eigs_bounds()
+ # Constant so that the minimum eigenvalue is represented exactly, since it contributes
+ # the most to the solution of the system. -1 to take into account the sign qubit
+ delta = self._get_delta(nl - neg_vals, lambda_min, lambda_max)
+ # Update evolution time
+ matrix_circuit.evolution_time = (
+ 2 * np.pi * delta / lambda_min / (2**neg_vals)
+ )
+ # Update the scaling of the solution
+ self.scaling = lambda_min
+ else:
+ delta = 1 / (2**nl)
+ print("The solution will be calculated up to a scaling factor.")
+
+ if self._exact_reciprocal:
+ reciprocal_circuit = ExactReciprocal(nl, delta, neg_vals=neg_vals)
+ # Update number of ancilla qubits
+ na = matrix_circuit.num_ancillas
+ else:
+ # Calculate breakpoints for the reciprocal approximation
+ num_values = 2**nl
+ constant = delta
+ a = int(round(num_values ** (2 / 3)))
+
+ # Calculate the degree of the polynomial and the number of intervals
+ r = 2 * constant / a + np.sqrt(np.abs(1 - (2 * constant / a) ** 2))
+ degree = min(
+ nb,
+ int(
+ np.log(
+ 1
+ + (
+ 16.23
+ * np.sqrt(np.log(r) ** 2 + (np.pi / 2) ** 2)
+ * kappa
+ * (2 * kappa - self._epsilon_r)
+ )
+ / self._epsilon_r
+ )
+ ),
+ )
+ num_intervals = int(np.ceil(np.log((num_values - 1) / a) / np.log(5)))
+
+ # Calculate breakpoints and polynomials
+ breakpoints = []
+ for i in range(0, num_intervals):
+ # Add the breakpoint to the list
+ breakpoints.append(a * (5**i))
+
+ # Define the right breakpoint of the interval
+ if i == num_intervals - 1:
+ breakpoints.append(num_values - 1)
+
+ reciprocal_circuit = PiecewiseChebyshev(
+ lambda x: np.arcsin(constant / x), degree, breakpoints, nl
+ )
+ na = max(matrix_circuit.num_ancillas, reciprocal_circuit.num_ancillas)
+
+ # Initialise the quantum registers
+ qb = QuantumRegister(nb) # right hand side and solution
+ ql = QuantumRegister(nl) # eigenvalue evaluation qubits
+ if na > 0:
+ qa = AncillaRegister(na) # ancilla qubits
+ qf = QuantumRegister(nf) # flag qubits
+
+ if na > 0:
+ qc = QuantumCircuit(qb, ql, qa, qf)
+ else:
+ qc = QuantumCircuit(qb, ql, qf)
+
+ # State preparation
+ qc.append(vector_circuit, qb[:])
+ # QPE
+ phase_estimation = PhaseEstimation(nl, matrix_circuit)
+ if na > 0:
+ qc.append(
+ phase_estimation, ql[:] + qb[:] + qa[: matrix_circuit.num_ancillas]
+ )
+ else:
+ qc.append(phase_estimation, ql[:] + qb[:])
+ # Conditioned rotation
+ if self._exact_reciprocal:
+ qc.append(reciprocal_circuit, ql[::-1] + [qf[0]])
+ else:
+ qc.append(
+ reciprocal_circuit.to_instruction(),
+ ql[:] + [qf[0]] + qa[: reciprocal_circuit.num_ancillas],
+ )
+ # QPE inverse
+ if na > 0:
+ qc.append(
+ phase_estimation.inverse(),
+ ql[:] + qb[:] + qa[: matrix_circuit.num_ancillas],
+ )
+ else:
+ qc.append(phase_estimation.inverse(), ql[:] + qb[:])
+ return qc
+
+ def solve(
+ self,
+ matrix: Union[List, np.ndarray, QuantumCircuit],
+ vector: Union[List, np.ndarray, QuantumCircuit],
+ observable: Optional[
+ Union[
+ LinearSystemObservable,
+ List[LinearSystemObservable],
+ ]
+ ] = None,
+ observable_circuit: Optional[
+ Union[QuantumCircuit, List[QuantumCircuit]]
+ ] = None,
+ post_processing: Optional[
+ Callable[[Union[float, List[float]], int, float], float]
+ ] = None,
+ ) -> LinearSolverResult:
+ """Tries to solve the given linear system of equations.
+
+ Args:
+ matrix: The matrix specifying the system, i.e. A in Ax=b.
+ vector: The vector specifying the right hand side of the equation in Ax=b.
+ observable: Optional information to be extracted from the solution.
+ Default is the probability of success of the algorithm.
+ observable_circuit: Optional circuit to be applied to the solution to extract
+ information. Default is `None`.
+ post_processing: Optional function to compute the value of the observable.
+ Default is the raw value of measuring the observable.
+
+ Raises:
+ ValueError: If an invalid combination of observable, observable_circuit and
+ post_processing is passed.
+
+ Returns:
+ The result object containing information about the solution vector of the linear
+ system.
+ """
+ # verify input
+ if observable is not None:
+ if observable_circuit is not None or post_processing is not None:
+ raise ValueError(
+ "If observable is passed, observable_circuit and post_processing cannot be set."
+ )
+
+ solution = LinearSolverResult()
+ solution.state = self.construct_circuit(matrix, vector)
+ solution.euclidean_norm = self._calculate_norm(solution.state)
+
+ if isinstance(observable, List):
+ observable_all, circuit_results_all = [], []
+ for obs in observable:
+ obs_i, circ_results_i = self._calculate_observable(
+ solution.state, obs, observable_circuit, post_processing
+ )
+ observable_all.append(obs_i)
+ circuit_results_all.append(circ_results_i)
+ solution.observable = observable_all
+ solution.circuit_results = circuit_results_all
+ elif observable is not None or observable_circuit is not None:
+ solution.observable, solution.circuit_results = self._calculate_observable(
+ solution.state, observable, observable_circuit, post_processing
+ )
+
+ return solution
diff --git a/qiskit_research/linear_solvers/linear_solver.py b/qiskit_research/linear_solvers/linear_solver.py
new file mode 100644
index 00000000..3bcec7b7
--- /dev/null
+++ b/qiskit_research/linear_solvers/linear_solver.py
@@ -0,0 +1,138 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""An abstract class for linear systems solvers."""
+
+from abc import ABC, abstractmethod
+from typing import Union, Optional, List, Callable
+import numpy as np
+
+from qiskit import QuantumCircuit
+from qiskit.algorithms.algorithm_result import AlgorithmResult
+
+from .observables.linear_system_observable import LinearSystemObservable
+
+# pylint: disable=too-few-public-methods
+
+
+class LinearSolverResult(AlgorithmResult):
+ """A base class for linear systems results.
+
+ The linear systems algorithms return an object of the type ``LinearSystemsResult``
+ with the information about the solution obtained.
+ """
+
+ def __init__(self) -> None:
+ super().__init__()
+
+ # Set the default to None, if the algorithm knows how to calculate it can override it.
+ self._state: Optional[Union[QuantumCircuit, np.ndarray]] = None
+ self._observable: Optional[Union[float, List[float]]] = None
+ self._euclidean_norm: Optional[float] = None
+ self._circuit_results: Optional[
+ Union[complex, List[complex], List[Union[complex, List[complex]]]]
+ ] = None
+
+ @property
+ def observable(self) -> Union[float, List[float]]:
+ """return the (list of) calculated observable(s)"""
+ return self._observable
+
+ @observable.setter
+ def observable(self, observable: Union[float, List[float]]) -> None:
+ """Set the value(s) of the observable(s).
+
+ Args:
+ observable: The new value(s) of the observable(s).
+ """
+ self._observable = observable
+
+ @property
+ def state(self) -> Union[QuantumCircuit, np.ndarray]:
+ """return either the circuit that prepares the solution or the solution as a vector"""
+ return self._state
+
+ @state.setter
+ def state(self, state: Union[QuantumCircuit, np.ndarray]) -> None:
+ """Set the solution state as either the circuit that prepares it or as a vector.
+
+ Args:
+ state: The new solution state.
+ """
+ self._state = state
+
+ @property
+ def euclidean_norm(self) -> float:
+ """return the euclidean norm if the algorithm knows how to calculate it"""
+ return self._euclidean_norm
+
+ @euclidean_norm.setter
+ def euclidean_norm(self, norm: float) -> None:
+ """Set the euclidean norm of the solution.
+
+ Args:
+ norm: The new euclidean norm of the solution.
+ """
+ self._euclidean_norm = norm
+
+ @property
+ def circuit_results(
+ self,
+ ) -> Union[complex, List[complex], List[Union[complex, List[complex]]]]:
+ """return the results from the circuits"""
+ return self._circuit_results
+
+ @circuit_results.setter
+ def circuit_results(
+ self,
+ results: Union[complex, List[complex], List[Union[complex, List[complex]]]],
+ ):
+ self._circuit_results = results
+
+
+class LinearSolver(ABC):
+ """An abstract class for linear system solvers in Qiskit."""
+
+ @abstractmethod
+ def solve(
+ self,
+ matrix: Union[np.ndarray, QuantumCircuit],
+ vector: Union[np.ndarray, QuantumCircuit],
+ observable: Optional[
+ Union[
+ LinearSystemObservable,
+ List[LinearSystemObservable],
+ ]
+ ] = None,
+ observable_circuit: Optional[
+ Union[QuantumCircuit, List[QuantumCircuit]]
+ ] = None,
+ post_processing: Optional[
+ Callable[[Union[float, List[float]], int, float], float]
+ ] = None,
+ ) -> LinearSolverResult:
+ """Solve the system and compute the observable(s)
+
+ Args:
+ matrix: The matrix specifying the system, i.e. A in Ax=b.
+ vector: The vector specifying the right hand side of the equation in Ax=b.
+ observable: Optional information to be extracted from the solution.
+ Default is the probability of success of the algorithm.
+ observable_circuit: Optional circuit to be applied to the solution to extract
+ information. Default is ``None``.
+ post_processing: Optional function to compute the value of the observable.
+ Default is the raw value of measuring the observable.
+
+ Returns:
+ The result of the linear system.
+ """
+ raise NotImplementedError
diff --git a/qiskit_research/linear_solvers/matrices/__init__.py b/qiskit_research/linear_solvers/matrices/__init__.py
new file mode 100644
index 00000000..a62dfdef
--- /dev/null
+++ b/qiskit_research/linear_solvers/matrices/__init__.py
@@ -0,0 +1,19 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""System matrices for Qiskit's linear solvers."""
+
+from .linear_system_matrix import LinearSystemMatrix
+from .numpy_matrix import NumPyMatrix
+from .tridiagonal_toeplitz import TridiagonalToeplitz
+
+__all__ = ["LinearSystemMatrix", "NumPyMatrix", "TridiagonalToeplitz"]
diff --git a/qiskit_research/linear_solvers/matrices/linear_system_matrix.py b/qiskit_research/linear_solvers/matrices/linear_system_matrix.py
new file mode 100644
index 00000000..ee0a4500
--- /dev/null
+++ b/qiskit_research/linear_solvers/matrices/linear_system_matrix.py
@@ -0,0 +1,130 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""An abstract class for matrices input to the linear systems solvers in Qiskit."""
+
+from abc import ABC, abstractmethod
+from typing import Tuple
+
+from qiskit import QuantumCircuit
+from qiskit.circuit.library import BlueprintCircuit
+
+
+class LinearSystemMatrix(BlueprintCircuit, ABC):
+ """Base class for linear system matrices."""
+
+ def __init__(
+ self,
+ num_state_qubits: int,
+ tolerance: float,
+ evolution_time: float,
+ name: str = "ls_matrix",
+ ) -> None:
+ """
+ Args:
+ num_state_qubits: the number of qubits where the unitary acts.
+ tolerance: the accuracy desired for the approximation
+ evolution_time: the time of the Hamiltonian simulation
+ name: The name of the object.
+ """
+ super().__init__(name=name)
+
+ # store parameters
+ self._num_state_qubits = num_state_qubits
+ self._reset_registers(num_state_qubits)
+ self.tolerance = tolerance
+ self.evolution_time = evolution_time
+
+ @property
+ def num_state_qubits(self) -> int:
+ r"""The number of state qubits representing the state :math:`|x\rangle`.
+
+ Returns:
+ The number of state qubits.
+ """
+ return self._num_state_qubits
+
+ @num_state_qubits.setter
+ def num_state_qubits(self, num_state_qubits: int) -> None:
+ """Set the number of state qubits.
+
+ Note that this may change the underlying quantum register, if the number of state qubits
+ changes.
+
+ Args:
+ num_state_qubits: The new number of qubits.
+ """
+ if num_state_qubits != self._num_state_qubits:
+ self._invalidate()
+ self._num_state_qubits = num_state_qubits
+ self._reset_registers(num_state_qubits)
+
+ @property
+ def tolerance(self) -> float:
+ """Return the error tolerance"""
+ return self._tolerance
+
+ @tolerance.setter
+ def tolerance(self, tolerance: float) -> None:
+ """Set the error tolerance
+ Args:
+ tolerance: The new error tolerance.
+ """
+ self._tolerance = tolerance
+
+ @property
+ def evolution_time(self) -> float:
+ """Return the time of the evolution."""
+ return self._evolution_time
+
+ @evolution_time.setter
+ def evolution_time(self, evolution_time: float) -> None:
+ """Set the time of the evolution.
+
+ Args:
+ evolution_time: The new time of the evolution.
+ """
+ self._evolution_time = evolution_time
+
+ @abstractmethod
+ def eigs_bounds(self) -> Tuple[float, float]:
+ """Return lower and upper bounds on the eigenvalues of the matrix."""
+ raise NotImplementedError
+
+ @abstractmethod
+ def condition_bounds(self) -> Tuple[float, float]:
+ """Return lower and upper bounds on the condition number of the matrix."""
+ raise NotImplementedError
+
+ @abstractmethod
+ def _reset_registers(self, num_state_qubits: int) -> None:
+ """Reset the registers according to the new number of state qubits.
+
+ Args:
+ num_state_qubits: The new number of qubits.
+ """
+ raise NotImplementedError
+
+ @abstractmethod
+ def power(self, power: int, matrix_power: bool = False) -> QuantumCircuit:
+ """Build powers of the circuit.
+
+ Args:
+ power: The power to raise this circuit to.
+ matrix_power: If True, the circuit is converted to a matrix and then the
+ matrix power is computed. If False, and ``power`` is a positive integer,
+ the implementation defaults to ``repeat``.
+
+ Returns:
+ The quantum circuit implementing powers of the unitary.
+ """
+ raise NotImplementedError
diff --git a/qiskit_research/linear_solvers/matrices/numpy_matrix.py b/qiskit_research/linear_solvers/matrices/numpy_matrix.py
new file mode 100644
index 00000000..0d1d9c63
--- /dev/null
+++ b/qiskit_research/linear_solvers/matrices/numpy_matrix.py
@@ -0,0 +1,215 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Hamiltonian simulation of matrices given as numpy arrays."""
+
+from typing import Tuple
+import numpy as np
+import scipy as sp
+
+from qiskit import QuantumCircuit, QuantumRegister
+
+from .linear_system_matrix import LinearSystemMatrix
+
+
+class NumPyMatrix(LinearSystemMatrix):
+ """Class of matrices given as a numpy array.
+
+ Examples:
+
+ .. jupyter-execute::
+
+ import numpy as np
+ from qiskit import QuantumCircuit
+ from qiskit_research.linear_solvers.matrices.numpy_matrix import NumPyMatrix
+
+ matrix = NumPyMatrix(np.array([[1 / 2, 1 / 6, 0, 0], [1 / 6, 1 / 2, 1 / 6, 0],
+ [0, 1 / 6, 1 / 2, 1 / 6], [0, 0, 1 / 6, 1 / 2]]))
+ power = 2
+
+ num_qubits = matrix.num_state_qubits
+ # Controlled power (as used within QPE)
+ pow_circ = matrix.power(power).control()
+ circ_qubits = pow_circ.num_qubits
+ qc = QuantumCircuit(circ_qubits)
+ qc.append(matrix.power(power).control(), list(range(circ_qubits)))
+ """
+
+ def __init__(
+ self,
+ matrix: np.ndarray,
+ tolerance: float = 1e-2,
+ evolution_time: float = 1.0,
+ name: str = "np_matrix",
+ ) -> None:
+ """
+ Args:
+ matrix: The matrix defining the linear system problem.
+ tolerance: The accuracy desired for the approximation.
+ evolution_time: The time of the Hamiltonian simulation.
+ name: The name of the object.
+ """
+
+ # define internal parameters
+ self._num_state_qubits = None
+ self._tolerance = None
+ self._evolution_time = None # makes sure the eigenvalues are contained in [0,1)
+ self._matrix = None
+
+ super().__init__(
+ num_state_qubits=int(np.log2(matrix.shape[0])),
+ tolerance=tolerance,
+ evolution_time=evolution_time,
+ name=name,
+ )
+
+ # store parameters
+ self.num_state_qubits = int(np.log2(matrix.shape[0]))
+ self.tolerance = tolerance
+ self.evolution_time = evolution_time
+ self.matrix = matrix
+
+ @property
+ def num_state_qubits(self) -> int:
+ r"""The number of state qubits representing the state :math:`|x\rangle`.
+
+ Returns:
+ The number of state qubits.
+ """
+ return self._num_state_qubits
+
+ @num_state_qubits.setter
+ def num_state_qubits(self, num_state_qubits: int) -> None:
+ """Set the number of state qubits.
+
+ Note that this may change the underlying quantum register, if the number of state qubits
+ changes.
+
+ Args:
+ num_state_qubits: The new number of qubits.
+ """
+ if num_state_qubits != self._num_state_qubits:
+ self._invalidate()
+ self._num_state_qubits = num_state_qubits
+ self._reset_registers(num_state_qubits)
+
+ @property
+ def tolerance(self) -> float:
+ """Return the error tolerance"""
+ return self._tolerance
+
+ @tolerance.setter
+ def tolerance(self, tolerance: float) -> None:
+ """Set the error tolerance
+ Args:
+ tolerance: The new error tolerance.
+ """
+ self._tolerance = tolerance
+
+ @property
+ def evolution_time(self) -> float:
+ """Return the time of the evolution."""
+ return self._evolution_time
+
+ @evolution_time.setter
+ def evolution_time(self, evolution_time: float) -> None:
+ """Set the time of the evolution.
+
+ Args:
+ evolution_time: The new time of the evolution.
+ """
+ self._evolution_time = evolution_time
+
+ @property
+ def matrix(self) -> np.ndarray:
+ """Return the matrix."""
+ return self._matrix
+
+ @matrix.setter
+ def matrix(self, matrix: np.ndarray) -> None:
+ """Set the matrix.
+
+ Args:
+ matrix: The new matrix.
+ """
+ self._matrix = matrix
+
+ def eigs_bounds(self) -> Tuple[float, float]:
+ """Return lower and upper bounds on the eigenvalues of the matrix."""
+ matrix_array = self.matrix
+ lambda_max = max(np.abs(np.linalg.eigvals(matrix_array)))
+ lambda_min = min(np.abs(np.linalg.eigvals(matrix_array)))
+ return lambda_min, lambda_max
+
+ def condition_bounds(self) -> Tuple[float, float]:
+ """Return lower and upper bounds on the condition number of the matrix."""
+ matrix_array = self.matrix
+ kappa = np.linalg.cond(matrix_array)
+ return kappa, kappa
+
+ def _check_configuration(self, raise_on_failure: bool = True) -> bool:
+ """Check if the current configuration is valid."""
+ valid = True
+
+ if self.matrix.shape[0] != self.matrix.shape[1]:
+ if raise_on_failure:
+ raise AttributeError("Input matrix must be square!")
+ return False
+ if np.log2(self.matrix.shape[0]) % 1 != 0:
+ if raise_on_failure:
+ raise AttributeError("Input matrix dimension must be 2^n!")
+ return False
+ if not np.allclose(self.matrix, self.matrix.conj().T):
+ if raise_on_failure:
+ raise AttributeError("Input matrix must be hermitian!")
+ return False
+
+ return valid
+
+ def _reset_registers(self, num_state_qubits: int) -> None:
+ """Reset the quantum registers.
+
+ Args:
+ num_state_qubits: The number of qubits to represent the matrix.
+ """
+ qr_state = QuantumRegister(num_state_qubits, "state")
+ self.qregs = [qr_state]
+
+ def _build(self) -> None:
+ """If not already built, build the circuit."""
+ if self._is_built:
+ return
+
+ super()._build()
+
+ self.compose(self.power(1), inplace=True)
+
+ def inverse(self):
+ return NumPyMatrix(self.matrix, evolution_time=-1 * self.evolution_time)
+
+ def power(self, power: int, matrix_power: bool = False) -> QuantumCircuit:
+ """Build powers of the circuit.
+
+ Args:
+ power: The power to raise this circuit to.
+ matrix_power: If True, the circuit is converted to a matrix and then the
+ matrix power is computed. If False, and ``power`` is a positive integer,
+ the implementation defaults to ``repeat``.
+
+ Returns:
+ The quantum circuit implementing powers of the unitary.
+ """
+ qc = QuantumCircuit(self.num_state_qubits)
+ evolved = sp.linalg.expm(1j * self.matrix * self.evolution_time)
+ # pylint: disable=no-member
+ qc.unitary(evolved, qc.qubits)
+ return qc.power(power)
diff --git a/qiskit_research/linear_solvers/matrices/tridiagonal_toeplitz.py b/qiskit_research/linear_solvers/matrices/tridiagonal_toeplitz.py
new file mode 100644
index 00000000..81c08b7f
--- /dev/null
+++ b/qiskit_research/linear_solvers/matrices/tridiagonal_toeplitz.py
@@ -0,0 +1,492 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Hamiltonian simulation of tridiagonal Toeplitz symmetric matrices."""
+
+from typing import Tuple, List
+import numpy as np
+from scipy.sparse import diags
+
+from qiskit.circuit import QuantumCircuit, QuantumRegister, AncillaRegister
+from qiskit.circuit.library import UGate, MCMTVChain
+
+from .linear_system_matrix import LinearSystemMatrix
+
+
+class TridiagonalToeplitz(LinearSystemMatrix):
+ r"""Class of tridiagonal Toeplitz symmetric matrices.
+
+ Given the main entry, :math:`a`, and the off diagonal entry, :math:`b`, the :math:`4\times 4`
+ dimensional tridiagonal Toeplitz symmetric matrix is
+
+ .. math::
+
+ \begin{pmatrix}
+ a & b & 0 & 0 \\
+ b & a & b & 0 \\
+ 0 & b & a & b \\
+ 0 & 0 & b & a
+ \end{pmatrix}.
+
+ Examples:
+
+ .. jupyter-execute::
+
+ import numpy as np
+ from qiskit import QuantumCircuit
+ from qiskit_research.linear_solvers.matrices import TridiagonalToeplitz
+
+ matrix = TridiagonalToeplitz(2, 1, -1 / 3)
+ power = 3
+
+ # Controlled power (as within QPE)
+ num_qubits = matrix.num_state_qubits
+ pow_circ = matrix.power(power).control()
+ circ_qubits = pow_circ.num_qubits
+ qc = QuantumCircuit(circ_qubits)
+ qc.append(matrix.power(power).control(), list(range(circ_qubits)))
+ """
+
+ def __init__(
+ self,
+ num_state_qubits: int,
+ main_diag: float,
+ off_diag: float,
+ tolerance: float = 1e-2,
+ evolution_time: float = 1.0,
+ trotter_steps: int = 1,
+ name: str = "tridi",
+ ) -> None:
+ """
+ Args:
+ num_state_qubits: the number of qubits where the unitary acts.
+ main_diag: the main diagonal entry
+ off_diag: the off diagonal entry
+ tolerance: the accuracy desired for the approximation
+ evolution_time: the time of the Hamiltonian simulation
+ trotter_steps: the number of Trotter steps
+ name: The name of the object.
+ """
+ # define internal parameters
+ self._main_diag = None
+ self._off_diag = None
+ self._tolerance = None
+ self._evolution_time = None # makes sure the eigenvalues are contained in [0,1)
+ self._trotter_steps = None
+
+ # store parameters
+ self.main_diag = main_diag
+ self.off_diag = off_diag
+ super().__init__(
+ num_state_qubits=num_state_qubits,
+ tolerance=tolerance,
+ evolution_time=evolution_time,
+ name=name,
+ )
+ self.trotter_steps = trotter_steps
+
+ @property
+ def num_state_qubits(self) -> int:
+ r"""The number of state qubits representing the state :math:`|x\rangle`.
+
+ Returns:
+ The number of state qubits.
+ """
+ return self._num_state_qubits
+
+ @num_state_qubits.setter
+ def num_state_qubits(self, num_state_qubits: int) -> None:
+ """Set the number of state qubits.
+
+ Note that this may change the underlying quantum register, if the number of state qubits
+ changes.
+
+ Args:
+ num_state_qubits: The new number of qubits.
+ """
+ if num_state_qubits != self._num_state_qubits:
+ self._invalidate()
+ self._num_state_qubits = num_state_qubits
+ self._reset_registers(num_state_qubits)
+
+ @property
+ def main_diag(self) -> float:
+ """Return the entry in the main diagonal."""
+ return self._main_diag
+
+ @main_diag.setter
+ def main_diag(self, main_diag: float) -> None:
+ """Set the entry in the main diagonal.
+ Args:
+ main_diag: The new entry in the main diagonal.
+ """
+ self._main_diag = main_diag
+
+ @property
+ def off_diag(self) -> float:
+ """Return the entry in the off diagonals."""
+ return self._off_diag
+
+ @off_diag.setter
+ def off_diag(self, off_diag: float) -> None:
+ """Set the entry in the off diagonals.
+ Args:
+ off_diag: The new entry in the main diagonal.
+ """
+ self._off_diag = off_diag
+
+ @property
+ def tolerance(self) -> float:
+ """Return the error tolerance"""
+ return self._tolerance
+
+ @tolerance.setter
+ def tolerance(self, tolerance: float) -> None:
+ """Set the error tolerance.
+ Args:
+ tolerance: The new error tolerance.
+ """
+ self._tolerance = tolerance
+
+ @property
+ def evolution_time(self) -> float:
+ """Return the time of the evolution."""
+ return self._evolution_time
+
+ @evolution_time.setter
+ def evolution_time(self, evolution_time: float) -> None:
+ """Set the time of the evolution and update the number of Trotter steps because the error
+ tolerance is a function of the evolution time and the number of trotter steps.
+
+ Args:
+ evolution_time: The new time of the evolution.
+ """
+ self._evolution_time = evolution_time
+ # Update the number of trotter steps. Max 7 for now, upper bounds too loose.
+ self.trotter_steps = int(
+ np.ceil(
+ np.sqrt(
+ ((evolution_time * np.abs(self.off_diag)) ** 3) / 2 / self.tolerance
+ )
+ )
+ )
+
+ @property
+ def trotter_steps(self) -> int:
+ """Return the number of trotter steps."""
+ return self._trotter_steps
+
+ @trotter_steps.setter
+ def trotter_steps(self, trotter_steps: int) -> None:
+ """Set the number of trotter steps.
+ Args:
+ trotter_steps: The new number of trotter steps.
+ """
+ self._trotter_steps = trotter_steps
+
+ @property
+ def matrix(self) -> np.ndarray:
+ """Returns the tridiagonal Toeplitz matrix built according to the main and off diagonal
+ entries."""
+ matrix = diags(
+ [self.off_diag, self.main_diag, self.off_diag],
+ [-1, 0, 1],
+ shape=(2**self.num_state_qubits, 2**self.num_state_qubits),
+ ).toarray()
+ return matrix
+
+ def eigs_bounds(self) -> Tuple[float, float]:
+ """Return lower and upper bounds on the absolute eigenvalues of the matrix."""
+ n_b = 2**self.num_state_qubits
+
+ # Calculate minimum and maximum of absolute value of eigenvalues
+ # according to the formula for Toeplitz 3-diagonal matrices
+
+ # For maximum it's enough to check border points of segment [1, n_b]
+ candidate_eig_ids = [1, n_b]
+
+ # Trying to add candidates near the minimum value of absolute eigenvalues
+ # function abs(main_diag - 2 * off_diag * cos(i * pi / (nb + 1))
+ if abs(self.main_diag) < 2 * abs(self.off_diag):
+ optimal_index = int(
+ np.arccos(self.main_diag / 2 / self.off_diag) / np.pi * (n_b + 1)
+ )
+
+ def add_candidate_index_if_valid(index_to_add: int) -> None:
+ if 1 <= index_to_add <= n_b:
+ candidate_eig_ids.append(index_to_add)
+
+ add_candidate_index_if_valid(optimal_index - 1)
+ add_candidate_index_if_valid(optimal_index)
+ add_candidate_index_if_valid(optimal_index + 1)
+
+ candidate_abs_eigs = np.abs(
+ [
+ self.main_diag - 2 * self.off_diag * np.cos(eig_id * np.pi / (n_b + 1))
+ for eig_id in candidate_eig_ids
+ ]
+ )
+
+ lambda_min = np.min(candidate_abs_eigs)
+ lambda_max = np.max(candidate_abs_eigs)
+ return lambda_min, lambda_max
+
+ def condition_bounds(self) -> Tuple[float, float]:
+ """Return lower and upper bounds on the condition number of the matrix."""
+ matrix_array = self.matrix
+ kappa = np.linalg.cond(matrix_array)
+ return kappa, kappa
+
+ def _check_configuration(self, raise_on_failure: bool = True) -> bool:
+ """Check if the current configuration is valid."""
+ valid = True
+
+ if self.trotter_steps < 1:
+ valid = False
+ if raise_on_failure:
+ raise AttributeError(
+ "The number of trotter steps should be a positive integer."
+ )
+ return False
+
+ return valid
+
+ def _reset_registers(self, num_state_qubits: int) -> None:
+ """Reset the quantum registers.
+
+ Args:
+ num_state_qubits: The number of qubits to represent the matrix.
+ """
+ qr_state = QuantumRegister(num_state_qubits, "state")
+ self.qregs = [qr_state]
+ self._ancillas: List[AncillaRegister] = []
+ self._qubits = qr_state[:]
+
+ if num_state_qubits > 1:
+ qr_ancilla = AncillaRegister(max(1, num_state_qubits - 1))
+ self.add_register(qr_ancilla)
+
+ def _build(self) -> None:
+ """If not already built, build the circuit."""
+ if self._is_built:
+ return
+
+ super()._build()
+
+ self.compose(self.power(1), inplace=True)
+
+ def _main_diag_circ(self, theta: float = 1) -> QuantumCircuit:
+ """Circuit implementing the matrix consisting of entries in the main diagonal.
+
+ Args:
+ theta: Scale factor for the main diagonal entries (e.g. evolution_time/trotter_steps).
+
+ Returns:
+ The quantum circuit implementing the matrix consisting of entries in the main diagonal.
+ """
+ theta *= self.main_diag
+ qc = QuantumCircuit(self.num_state_qubits, name="main_diag")
+ qc.x(0)
+ qc.p(theta, 0)
+ qc.x(0)
+ qc.p(theta, 0)
+
+ # pylint: disable=unused-argument
+ def control(num_ctrl_qubits=1, label=None, ctrl_state=None):
+ qc_control = QuantumCircuit(self.num_state_qubits + 1, name="main_diag")
+ qc_control.p(theta, 0)
+ return qc_control
+
+ qc.control = control
+ return qc
+
+ def _off_diag_circ(self, theta: float = 1) -> QuantumCircuit:
+ """Circuit implementing the matrix consisting of entries in the off diagonals.
+
+ Args:
+ theta: Scale factor for the off diagonal entries (e.g. evolution_time/trotter_steps).
+
+ Returns:
+ The quantum circuit implementing the matrix consisting of entries in the off diagonals.
+ """
+ theta *= self.off_diag
+
+ qr = QuantumRegister(self.num_state_qubits)
+ if self.num_state_qubits > 1:
+ qr_ancilla = AncillaRegister(max(1, self.num_state_qubits - 2))
+ qc = QuantumCircuit(qr, qr_ancilla, name="off_diags")
+ else:
+ qc = QuantumCircuit(qr, name="off_diags")
+ qr_ancilla = None
+
+ qc.u(-2 * theta, 3 * np.pi / 2, np.pi / 2, qr[0])
+
+ for i in range(0, self.num_state_qubits - 1):
+ q_controls = []
+ qc.cx(qr[i], qr[i + 1])
+ q_controls.append(qr[i + 1])
+
+ # Now we want controlled by 0
+ qc.x(qr[i])
+ for j in range(i, 0, -1):
+ qc.cx(qr[i], qr[j - 1])
+ q_controls.append(qr[j - 1])
+ qc.x(qr[i])
+
+ # Multicontrolled rotation
+ if len(q_controls) > 1:
+ ugate = UGate(-2 * theta, 3 * np.pi / 2, np.pi / 2)
+ qc.append(
+ MCMTVChain(ugate, len(q_controls), 1),
+ q_controls[:] + [qr[i]] + qr_ancilla[: len(q_controls) - 1],
+ )
+ else:
+ qc.cu(-2 * theta, 3 * np.pi / 2, np.pi / 2, 0, q_controls[0], qr[i])
+
+ # Uncompute
+ qc.x(qr[i])
+ for j in range(0, i):
+ qc.cx(qr[i], qr[j])
+ qc.x(qr[i])
+ qc.cx(qr[i], qr[i + 1])
+
+ # pylint: disable=unused-argument
+ def control(num_ctrl_qubits=1, label=None, ctrl_state=None):
+ qr_state = QuantumRegister(self.num_state_qubits + 1)
+ if self.num_state_qubits > 1:
+ qr_ancilla = AncillaRegister(max(1, self.num_state_qubits - 1))
+ qc_control = QuantumCircuit(qr_state, qr_ancilla, name="off_diags")
+ else:
+ qc_control = QuantumCircuit(qr_state, name="off_diags")
+ qr_ancilla = None
+ # Control will be qr[0]
+ q_control = qr_state[0]
+ qr = qr_state[1:]
+ qc_control.cu(-2 * theta, 3 * np.pi / 2, np.pi / 2, 0, q_control, qr[0])
+
+ for i in range(0, self.num_state_qubits - 1):
+ q_controls = []
+ q_controls.append(q_control)
+ qc_control.cx(qr[i], qr[i + 1])
+ q_controls.append(qr[i + 1])
+
+ # Now we want controlled by 0
+ qc_control.x(qr[i])
+ for j in range(i, 0, -1):
+ qc_control.cx(qr[i], qr[j - 1])
+ q_controls.append(qr[j - 1])
+ qc_control.x(qr[i])
+
+ # Multicontrolled x rotation
+ if len(q_controls) > 1:
+ ugate = UGate(-2 * theta, 3 * np.pi / 2, np.pi / 2)
+ qc_control.append(
+ MCMTVChain(ugate, len(q_controls), 1).to_gate(),
+ q_controls[:] + [qr[i]] + qr_ancilla[: len(q_controls) - 1],
+ )
+ else:
+ qc_control.cu(
+ -2 * theta, 3 * np.pi / 2, np.pi / 2, 0, q_controls[0], qr[i]
+ )
+
+ # Uncompute
+ qc_control.x(qr[i])
+ for j in range(0, i):
+ qc_control.cx(qr[i], qr[j])
+ qc_control.x(qr[i])
+ qc_control.cx(qr[i], qr[i + 1])
+ return qc_control
+
+ qc.control = control
+ return qc
+
+ def inverse(self):
+ return TridiagonalToeplitz(
+ self.num_state_qubits,
+ self.main_diag,
+ self.off_diag,
+ evolution_time=-1 * self.evolution_time,
+ )
+
+ def power(self, power: int, matrix_power: bool = False) -> QuantumCircuit:
+ """Build powers of the circuit.
+
+ Args:
+ power: The power to raise this circuit to.
+ matrix_power: If True, the circuit is converted to a matrix and then the
+ matrix power is computed. If False, and ``power`` is a positive integer,
+ the implementation defaults to ``repeat``.
+
+ Returns:
+ The quantum circuit implementing powers of the unitary.
+ """
+ qc_raw = QuantumCircuit(self.num_state_qubits)
+
+ # pylint: disable=unused-argument
+ def control(num_ctrl_qubits=1, label=None, ctrl_state=None):
+ qr_state = QuantumRegister(self.num_state_qubits + 1, "state")
+ if self.num_state_qubits > 1:
+ qr_ancilla = AncillaRegister(max(1, self.num_state_qubits - 1))
+ qc = QuantumCircuit(qr_state, qr_ancilla, name="exp(iHk)")
+ else:
+ qc = QuantumCircuit(qr_state, name="exp(iHk)")
+ qr_ancilla = None
+ # Control will be qr[0]
+ q_control = qr_state[0]
+ qr = qr_state[1:]
+ # A1 commutes, so one application with evolution_time*2^{j} to the last qubit is enough
+ qc.append(
+ self._main_diag_circ(self.evolution_time * power).control().to_gate(),
+ [q_control] + qr[:],
+ )
+
+ # Update trotter steps to compensate the error
+ trotter_steps_new = int(np.ceil(np.sqrt(power) * self.trotter_steps))
+
+ # exp(iA2t/2m)
+ qc.u(
+ self.off_diag * self.evolution_time * power / trotter_steps_new,
+ 3 * np.pi / 2,
+ np.pi / 2,
+ qr[0],
+ )
+ # for _ in range(power):
+ for _ in range(0, trotter_steps_new):
+ if qr_ancilla:
+ qc.append(
+ self._off_diag_circ(
+ self.evolution_time * power / trotter_steps_new
+ )
+ .control()
+ .to_gate(),
+ [q_control] + qr[:] + qr_ancilla[:],
+ )
+ else:
+ qc.append(
+ self._off_diag_circ(
+ self.evolution_time * power / trotter_steps_new
+ )
+ .control()
+ .to_gate(),
+ [q_control] + qr[:],
+ )
+ # exp(-iA2t/2m)
+ qc.u(
+ -self.off_diag * self.evolution_time * power / trotter_steps_new,
+ 3 * np.pi / 2,
+ np.pi / 2,
+ qr[0],
+ )
+ return qc
+
+ qc_raw.control = control
+ return qc_raw
diff --git a/qiskit_research/linear_solvers/numpy_linear_solver.py b/qiskit_research/linear_solvers/numpy_linear_solver.py
new file mode 100644
index 00000000..4e7b2d77
--- /dev/null
+++ b/qiskit_research/linear_solvers/numpy_linear_solver.py
@@ -0,0 +1,105 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""The Numpy LinearSolver algorithm (classical)."""
+
+from typing import List, Union, Optional, Callable
+import numpy as np
+
+from qiskit import QuantumCircuit
+from qiskit.quantum_info import Operator, Statevector
+
+from .linear_solver import LinearSolverResult, LinearSolver
+from .observables.linear_system_observable import LinearSystemObservable
+
+# pylint: disable=too-few-public-methods
+
+
+class NumPyLinearSolver(LinearSolver):
+ """The Numpy Linear Solver algorithm (classical).
+
+ This linear system solver computes the exact value of the given observable(s) or the full
+ solution vector if no observable is specified.
+
+ Examples:
+
+ .. jupyter-execute::
+
+ import numpy as np
+ from qiskit_research.linear_solvers import NumPyLinearSolver
+ from qiskit_research.linear_solvers.matrices import TridiagonalToeplitz
+ from qiskit_research.linear_solvers.observables import MatrixFunctional
+
+ matrix = TridiagonalToeplitz(2, 1, 1 / 3, trotter_steps=2)
+ right_hand_side = [1.0, -2.1, 3.2, -4.3]
+ observable = MatrixFunctional(1, 1 / 2)
+ rhs = right_hand_side / np.linalg.norm(right_hand_side)
+
+ np_solver = NumPyLinearSolver()
+ solution = np_solver.solve(matrix, rhs, observable)
+ result = solution.observable
+ """
+
+ def solve(
+ self,
+ matrix: Union[np.ndarray, QuantumCircuit],
+ vector: Union[np.ndarray, QuantumCircuit],
+ observable: Optional[
+ Union[
+ LinearSystemObservable,
+ List[LinearSystemObservable],
+ ]
+ ] = None,
+ observable_circuit: Optional[
+ Union[QuantumCircuit, List[QuantumCircuit]]
+ ] = None,
+ post_processing: Optional[
+ Callable[[Union[float, List[float]], int, float], float]
+ ] = None,
+ ) -> LinearSolverResult:
+ """Solve classically the linear system and compute the observable(s)
+
+ Args:
+ matrix: The matrix specifying the system, i.e. A in Ax=b.
+ vector: The vector specifying the right hand side of the equation in Ax=b.
+ observable: Optional information to be extracted from the solution.
+ Default is the probability of success of the algorithm.
+ observable_circuit: Optional circuit to be applied to the solution to extract
+ information. Default is ``None``.
+ post_processing: Optional function to compute the value of the observable.
+ Default is the raw value of measuring the observable.
+
+ Returns:
+ The result of the linear system.
+ """
+ # Check if either matrix or vector are QuantumCircuits and get the array from them
+ if isinstance(vector, QuantumCircuit):
+ vector = Statevector(vector).data
+ if isinstance(matrix, QuantumCircuit):
+ if hasattr(matrix, "matrix"):
+ matrix = matrix.matrix
+ else:
+ matrix = Operator(matrix).data
+
+ solution_vector = np.linalg.solve(matrix, vector)
+ solution = LinearSolverResult()
+ solution.state = solution_vector
+ if observable is not None:
+ if isinstance(observable, list):
+ solution.observable = []
+ for obs in observable:
+ solution.observable.append(
+ obs.evaluate_classically(solution_vector)
+ )
+ else:
+ solution.observable = observable.evaluate_classically(solution_vector)
+ solution.euclidean_norm = float(np.linalg.norm(solution_vector))
+ return solution
diff --git a/qiskit_research/linear_solvers/observables/__init__.py b/qiskit_research/linear_solvers/observables/__init__.py
new file mode 100644
index 00000000..1dba66cf
--- /dev/null
+++ b/qiskit_research/linear_solvers/observables/__init__.py
@@ -0,0 +1,19 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Observables for Qiskit's linear solvers."""
+
+from .linear_system_observable import LinearSystemObservable
+from .absolute_average import AbsoluteAverage
+from .matrix_functional import MatrixFunctional
+
+__all__ = ["LinearSystemObservable", "AbsoluteAverage", "MatrixFunctional"]
diff --git a/qiskit_research/linear_solvers/observables/absolute_average.py b/qiskit_research/linear_solvers/observables/absolute_average.py
new file mode 100644
index 00000000..281d50fe
--- /dev/null
+++ b/qiskit_research/linear_solvers/observables/absolute_average.py
@@ -0,0 +1,129 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""The absolute value of the average of a linear system of equations solution."""
+
+from typing import Union, List
+import numpy as np
+
+from qiskit import QuantumCircuit
+from qiskit.opflow import I, Z, TensoredOp
+from qiskit.quantum_info import Statevector
+
+from .linear_system_observable import LinearSystemObservable
+
+
+class AbsoluteAverage(LinearSystemObservable):
+ r"""An observable for the absolute average of a linear system of equations solution.
+
+ For a vector :math:`x=(x_1,...,x_N)`, the absolute average is defined as
+ :math:`\abs{\frac{1}{N}\sum_{i=1}^{N}x_i}`.
+
+ Examples:
+
+ .. jupyter-execute::
+
+ import numpy as np
+ from qiskit import QuantumCircuit
+ from qiskit_research.linear_solvers.observables.absolute_average import \
+ AbsoluteAverage
+ from qiskit.opflow import StateFn
+
+ observable = AbsoluteAverage()
+ vector = [1.0, -2.1, 3.2, -4.3]
+
+ init_state = vector / np.linalg.norm(vector)
+ num_qubits = int(np.log2(len(vector)))
+
+ qc = QuantumCircuit(num_qubits)
+ qc.isometry(init_state, list(range(num_qubits)), None)
+ qc.append(observable.observable_circuit(num_qubits), list(range(num_qubits)))
+
+ # Observable operator
+ observable_op = observable.observable(num_qubits)
+ state_vec = (~StateFn(observable_op) @ StateFn(qc)).eval()
+
+ # Obtain result
+ result = observable.post_processing(state_vec, num_qubits)
+
+ # Obtain analytical evaluation
+ exact = observable.evaluate_classically(init_state)
+ """
+
+ def observable(self, num_qubits: int) -> Union[TensoredOp, List[TensoredOp]]:
+ """The observable operator.
+
+ Args:
+ num_qubits: The number of qubits on which the observable will be applied.
+
+ Returns:
+ The observable as a sum of Pauli strings.
+ """
+ zero_op = (I + Z) / 2
+ return TensoredOp(num_qubits * [zero_op])
+
+ def observable_circuit(
+ self, num_qubits: int
+ ) -> Union[QuantumCircuit, List[QuantumCircuit]]:
+ """The circuit implementing the absolute average observable.
+
+ Args:
+ num_qubits: The number of qubits on which the observable will be applied.
+
+ Returns:
+ The observable as a QuantumCircuit.
+ """
+ qc = QuantumCircuit(num_qubits)
+ qc.h(qc.qubits)
+ return qc
+
+ def post_processing(
+ self, solution: Union[float, List[float]], num_qubits: int, scaling: float = 1
+ ) -> float:
+ """Evaluates the absolute average on the solution to the linear system.
+
+ Args:
+ solution: The probability calculated from the circuit and the observable.
+ num_qubits: The number of qubits where the observable was applied.
+ scaling: Scaling of the solution.
+
+ Returns:
+ The value of the absolute average.
+
+ Raises:
+ ValueError: If the input is not in the correct format.
+ """
+ if isinstance(solution, list):
+ if len(solution) == 1:
+ solution = solution[0]
+ else:
+ raise ValueError(
+ "Solution probability must be given as a single value."
+ )
+
+ return np.real(np.sqrt(solution / (2**num_qubits)) / scaling)
+
+ def evaluate_classically(
+ self, solution: Union[np.ndarray, QuantumCircuit]
+ ) -> float:
+ """Evaluates the given observable on the solution to the linear system.
+
+ Args:
+ solution: The solution to the system as a numpy array or the circuit that prepares it.
+
+ Returns:
+ The value of the observable.
+ """
+ # Check if it is QuantumCircuits and get the array from them
+ if isinstance(solution, QuantumCircuit):
+ solution = Statevector(solution).data
+ return np.abs(np.mean(solution))
diff --git a/qiskit_research/linear_solvers/observables/linear_system_observable.py b/qiskit_research/linear_solvers/observables/linear_system_observable.py
new file mode 100644
index 00000000..123d3fcf
--- /dev/null
+++ b/qiskit_research/linear_solvers/observables/linear_system_observable.py
@@ -0,0 +1,81 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""An abstract class for linear systems solvers in Qiskit's aqua module."""
+
+from abc import ABC, abstractmethod
+from typing import Union, List
+import numpy as np
+
+from qiskit import QuantumCircuit
+from qiskit.opflow import TensoredOp
+
+
+class LinearSystemObservable(ABC):
+ """An abstract class for linear system observables in Qiskit."""
+
+ @abstractmethod
+ def observable(self, num_qubits: int) -> Union[TensoredOp, List[TensoredOp]]:
+ """The observable operator.
+
+ Args:
+ num_qubits: The number of qubits on which the observable will be applied.
+
+ Returns:
+ The observable as a sum of Pauli strings.
+ """
+ raise NotImplementedError
+
+ @abstractmethod
+ def observable_circuit(
+ self, num_qubits: int
+ ) -> Union[QuantumCircuit, List[QuantumCircuit]]:
+ """The circuit implementing the observable.
+
+ Args:
+ num_qubits: The number of qubits on which the observable will be applied.
+
+ Returns:
+ The observable as a QuantumCircuit.
+ """
+ raise NotImplementedError
+
+ @abstractmethod
+ def post_processing(
+ self, solution: Union[float, List[float]], num_qubits: int, scaling: float = 1
+ ) -> float:
+ """Evaluates the given observable on the solution to the linear system.
+
+ Args:
+ solution: The probability calculated from the circuit and the observable.
+ num_qubits: The number of qubits where the observable was applied.
+ scaling: Scaling of the solution.
+
+ Returns:
+ The value of the observable.
+ """
+ raise NotImplementedError
+
+ @abstractmethod
+ def evaluate_classically(
+ self, solution: Union[np.ndarray, QuantumCircuit]
+ ) -> float:
+ """Calculates the analytical value of the given observable from the solution vector to the
+ linear system.
+
+ Args:
+ solution: The solution to the system as a numpy array or the circuit that prepares it.
+
+ Returns:
+ The value of the observable.
+ """
+ raise NotImplementedError
diff --git a/qiskit_research/linear_solvers/observables/matrix_functional.py b/qiskit_research/linear_solvers/observables/matrix_functional.py
new file mode 100644
index 00000000..cd7c37e1
--- /dev/null
+++ b/qiskit_research/linear_solvers/observables/matrix_functional.py
@@ -0,0 +1,183 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""The matrix functional of the vector solution to the linear systems."""
+
+from typing import Union, List
+import numpy as np
+from scipy.sparse import diags
+
+from qiskit import QuantumCircuit
+from qiskit.quantum_info import Statevector
+from qiskit.opflow import I, Z, TensoredOp
+
+from .linear_system_observable import LinearSystemObservable
+
+
+class MatrixFunctional(LinearSystemObservable):
+ """A class for the matrix functional of the vector solution to the linear systems.
+
+ Examples:
+
+ .. jupyter-execute::
+
+ import numpy as np
+ from qiskit import QuantumCircuit
+ from qiskit_research.linear_solvers.observables.matrix_functional import \
+ MatrixFunctional
+ from qiskit.transpiler.passes import RemoveResetInZeroState
+ from qiskit.opflow import StateFn
+
+ tpass = RemoveResetInZeroState()
+
+ vector = [1.0, -2.1, 3.2, -4.3]
+ observable = MatrixFunctional(1, -1 / 3)
+
+ init_state = vector / np.linalg.norm(vector)
+ num_qubits = int(np.log2(len(vector)))
+
+ # Get observable circuits
+ obs_circuits = observable.observable_circuit(num_qubits)
+ qcs = []
+ for obs_circ in obs_circuits:
+ qc = QuantumCircuit(num_qubits)
+ qc.isometry(init_state, list(range(num_qubits)), None)
+ qc.append(obs_circ, list(range(num_qubits)))
+ qcs.append(tpass(qc.decompose()))
+
+ # Get observables
+ observable_ops = observable.observable(num_qubits)
+ state_vecs = []
+ # First is the norm
+ state_vecs.append((~StateFn(observable_ops[0]) @ StateFn(qcs[0])).eval())
+ for i in range(1, len(observable_ops), 2):
+ state_vecs += [(~StateFn(observable_ops[i]) @ StateFn(qcs[i])).eval(),
+ (~StateFn(observable_ops[i + 1]) @ StateFn(qcs[i + 1])).eval()]
+
+ # Obtain result
+ result = observable.post_processing(state_vecs, num_qubits)
+
+ # Obtain analytical evaluation
+ exact = observable.evaluate_classically(init_state)
+ """
+
+ def __init__(self, main_diag: float, off_diag: float) -> None:
+ """
+ Args:
+ main_diag: The main diagonal of the tridiagonal Toeplitz symmetric matrix to compute
+ the functional.
+ off_diag: The off diagonal of the tridiagonal Toeplitz symmetric matrix to compute
+ the functional.
+ """
+ self._main_diag = main_diag
+ self._off_diag = off_diag
+
+ def observable(self, num_qubits: int) -> Union[TensoredOp, List[TensoredOp]]:
+ """The observable operators.
+
+ Args:
+ num_qubits: The number of qubits on which the observable will be applied.
+
+ Returns:
+ The observable as a list of sums of Pauli strings.
+ """
+ zero_op = (I + Z) / 2
+ one_op = (I - Z) / 2
+ observables = []
+ # First we measure the norm of x
+ observables.append(I ^ num_qubits)
+ for i in range(num_qubits):
+ j = num_qubits - i - 1
+
+ # TODO this if can be removed once the bug in Opflow is fixed where
+ # TensoredOp([X, TensoredOp([])]).eval() ends up in infinite recursion
+ if i > 0:
+ observables += [
+ (I ^ j) ^ zero_op ^ TensoredOp(i * [one_op]),
+ (I ^ j) ^ one_op ^ TensoredOp(i * [one_op]),
+ ]
+ else:
+ observables += [(I ^ j) ^ zero_op, (I ^ j) ^ one_op]
+
+ return observables
+
+ def observable_circuit(
+ self, num_qubits: int
+ ) -> Union[QuantumCircuit, List[QuantumCircuit]]:
+ """The circuits to implement the matrix functional observable.
+
+ Args:
+ num_qubits: The number of qubits on which the observable will be applied.
+
+ Returns:
+ The observable as a list of QuantumCircuits.
+ """
+ qcs = []
+ # Again, the first value in the list will correspond to the norm of x
+ qcs.append(QuantumCircuit(num_qubits))
+ for i in range(0, num_qubits):
+ qc = QuantumCircuit(num_qubits)
+ for j in range(0, i):
+ qc.cx(i, j)
+ qc.h(i)
+ qcs += [qc, qc]
+
+ return qcs
+
+ def post_processing(
+ self, solution: Union[float, List[float]], num_qubits: int, scaling: float = 1
+ ) -> float:
+ """Evaluates the matrix functional on the solution to the linear system.
+
+ Args:
+ solution: The list of probabilities calculated from the circuit and the observable.
+ num_qubits: The number of qubits where the observable was applied.
+ scaling: Scaling of the solution.
+
+ Returns:
+ The value of the absolute average.
+
+ Raises:
+ ValueError: If the input is not in the correct format.
+ """
+ if not isinstance(solution, list):
+ raise ValueError("Solution probabilities must be given in list form.")
+
+ # Calculate the value from the off-diagonal elements
+ off_val = 0.0
+ for i in range(1, len(solution), 2):
+ off_val += (solution[i] - solution[i + 1]) / (scaling**2)
+ main_val = solution[0] / (scaling**2)
+ return np.real(self._main_diag * main_val + self._off_diag * off_val)
+
+ def evaluate_classically(
+ self, solution: Union[np.ndarray, QuantumCircuit]
+ ) -> float:
+ """Evaluates the given observable on the solution to the linear system.
+
+ Args:
+ solution: The solution to the system as a numpy array or the circuit that prepares it.
+
+ Returns:
+ The value of the observable.
+ """
+ # Check if it is QuantumCircuits and get the array from them
+ if isinstance(solution, QuantumCircuit):
+ solution = Statevector(solution).data
+
+ matrix = diags(
+ [self._off_diag, self._main_diag, self._off_diag],
+ [-1, 0, 1],
+ shape=(len(solution), len(solution)),
+ ).toarray()
+
+ return np.dot(solution.transpose(), np.dot(matrix, solution))
diff --git a/test/linear_solvers/test_linear_solvers.py b/test/linear_solvers/test_linear_solvers.py
new file mode 100644
index 00000000..5f62db43
--- /dev/null
+++ b/test/linear_solvers/test_linear_solvers.py
@@ -0,0 +1,333 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Test the quantum linear system solver algorithm."""
+
+import unittest
+from scipy.linalg import expm
+import numpy as np
+from ddt import ddt, idata, unpack
+from qiskit import BasicAer, QuantumCircuit
+from qiskit_research.linear_solvers.hhl import HHL
+from qiskit_research.linear_solvers.matrices.tridiagonal_toeplitz import (
+ TridiagonalToeplitz,
+)
+from qiskit_research.linear_solvers.matrices.numpy_matrix import NumPyMatrix
+from qiskit_research.linear_solvers.observables.absolute_average import AbsoluteAverage
+from qiskit_research.linear_solvers.observables.matrix_functional import (
+ MatrixFunctional,
+)
+from qiskit.circuit.library.arithmetic.exact_reciprocal import ExactReciprocal
+from qiskit.quantum_info import Operator, partial_trace
+from qiskit.opflow import I, Z, StateFn
+from qiskit.utils import QuantumInstance
+from qiskit import quantum_info
+
+
+@ddt
+class TestMatrices(unittest.TestCase):
+ """Tests based on the matrices classes.
+
+ This class tests
+ * the constructed circuits
+ """
+
+ @idata(
+ [
+ [TridiagonalToeplitz(2, 1, -1 / 3)],
+ [TridiagonalToeplitz(3, 2, 1), 1.1, 3],
+ [
+ NumPyMatrix(
+ np.array(
+ [
+ [1 / 2, 1 / 6, 0, 0],
+ [1 / 6, 1 / 2, 1 / 6, 0],
+ [0, 1 / 6, 1 / 2, 1 / 6],
+ [0, 0, 1 / 6, 1 / 2],
+ ]
+ )
+ )
+ ],
+ ]
+ )
+ @unpack
+ def test_matrices(self, matrix, time=1.0, power=1):
+ """Test the different matrix classes."""
+ matrix.evolution_time = time
+
+ num_qubits = matrix.num_state_qubits
+ pow_circ = matrix.power(power).control()
+ circ_qubits = pow_circ.num_qubits
+ qc = QuantumCircuit(circ_qubits)
+ qc.append(matrix.power(power).control(), list(range(circ_qubits)))
+ # extract the parts of the circuit matrix corresponding to TridiagonalToeplitz
+ zero_op = (I + Z) / 2
+ one_op = (I - Z) / 2
+ proj = Operator(
+ (zero_op ^ pow_circ.num_ancillas) ^ (I ^ num_qubits) ^ one_op
+ ).data
+ circ_matrix = Operator(qc).data
+ approx_exp = partial_trace(
+ np.dot(proj, circ_matrix), [0] + list(range(num_qubits + 1, circ_qubits))
+ ).data
+
+ exact_exp = expm(1j * matrix.evolution_time * power * matrix.matrix)
+ np.testing.assert_array_almost_equal(approx_exp, exact_exp, decimal=2)
+
+ @idata(
+ [
+ [TridiagonalToeplitz(2, 1.5, 2.5)],
+ [TridiagonalToeplitz(4, -1, 1.6)],
+ ]
+ )
+ @unpack
+ def test_eigs_bounds(self, matrix):
+ """Test the capability of TridiagonalToeplitz matrix class
+ to find accurate absolute eigenvalues bounds."""
+
+ matrix_lambda_min, matrix_lambda_max = matrix.eigs_bounds()
+
+ numpy_matrix = matrix.matrix
+ eigenvalues, _ = np.linalg.eig(numpy_matrix)
+ abs_eigenvalues = np.abs(eigenvalues)
+ exact_lambda_min = np.min(abs_eigenvalues)
+ exact_lambda_max = np.max(abs_eigenvalues)
+
+ np.testing.assert_almost_equal(matrix_lambda_min, exact_lambda_min, decimal=6)
+ np.testing.assert_almost_equal(matrix_lambda_max, exact_lambda_max, decimal=6)
+
+
+@ddt
+class TestObservables(unittest.TestCase):
+ """Tests based on the observables classes.
+
+ This class tests
+ * the constructed circuits
+ """
+
+ @idata(
+ [
+ [AbsoluteAverage(), [1.0, -2.1, 3.2, -4.3]],
+ [AbsoluteAverage(), [-9 / 4, -0.3, 8 / 7, 10, -5, 11.1, 13 / 11, -27 / 12]],
+ ]
+ )
+ @unpack
+ def test_absolute_average(self, observable, vector):
+ """Test the absolute average observable."""
+ init_state = vector / np.linalg.norm(vector)
+ num_qubits = int(np.log2(len(vector)))
+
+ qc = QuantumCircuit(num_qubits)
+ qc.isometry(init_state, list(range(num_qubits)), None)
+ qc.append(observable.observable_circuit(num_qubits), list(range(num_qubits)))
+
+ # Observable operator
+ observable_op = observable.observable(num_qubits)
+ state_vec = (~StateFn(observable_op) @ StateFn(qc)).eval()
+
+ # Obtain result
+ result = observable.post_processing(state_vec, num_qubits)
+
+ # Obtain analytical evaluation
+ exact = observable.evaluate_classically(init_state)
+
+ np.testing.assert_almost_equal(result, exact, decimal=2)
+
+ @idata(
+ [
+ [MatrixFunctional(1, -1 / 3), [1.0, -2.1, 3.2, -4.3]],
+ [
+ MatrixFunctional(2 / 3, 11 / 7),
+ [-9 / 4, -0.3, 8 / 7, 10, -5, 11.1, 13 / 11, -27 / 12],
+ ],
+ ]
+ )
+ @unpack
+ def test_matrix_functional(self, observable, vector):
+ """Test the matrix functional class."""
+ from qiskit.transpiler.passes import RemoveResetInZeroState
+
+ tpass = RemoveResetInZeroState()
+
+ init_state = vector / np.linalg.norm(vector)
+ num_qubits = int(np.log2(len(vector)))
+
+ # Get observable circuits
+ obs_circuits = observable.observable_circuit(num_qubits)
+ qcs = []
+ for obs_circ in obs_circuits:
+ qc = QuantumCircuit(num_qubits)
+ qc.isometry(init_state, list(range(num_qubits)), None)
+ qc.append(obs_circ, list(range(num_qubits)))
+ qcs.append(tpass(qc.decompose()))
+
+ # Get observables
+ observable_ops = observable.observable(num_qubits)
+ state_vecs = []
+ # First is the norm
+ state_vecs.append((~StateFn(observable_ops[0]) @ StateFn(qcs[0])).eval())
+ for i in range(1, len(observable_ops), 2):
+ state_vecs += [
+ (~StateFn(observable_ops[i]) @ StateFn(qcs[i])).eval(),
+ (~StateFn(observable_ops[i + 1]) @ StateFn(qcs[i + 1])).eval(),
+ ]
+
+ # Obtain result
+ result = observable.post_processing(state_vecs, num_qubits)
+
+ # Obtain analytical evaluation
+ exact = observable.evaluate_classically(init_state)
+
+ np.testing.assert_almost_equal(result, exact, decimal=2)
+
+
+@ddt
+class TestReciprocal(unittest.TestCase):
+ """Tests based on the reciprocal classes.
+
+ This class tests
+ * the constructed circuits
+ """
+
+ @idata([[2, 0.1, False], [3, 1 / 9, True]])
+ @unpack
+ def test_exact_reciprocal(self, num_qubits, scaling, neg_vals):
+ """Test the ExactReciprocal class."""
+ reciprocal = ExactReciprocal(num_qubits + neg_vals, scaling, neg_vals)
+
+ qc = QuantumCircuit(num_qubits + 1 + neg_vals)
+ qc.h(list(range(num_qubits)))
+ # If negative eigenvalues, set the sign qubit to 1
+ if neg_vals:
+ qc.x(num_qubits)
+ qc.append(reciprocal, list(range(num_qubits + 1 + neg_vals)))
+
+ # Create the operator 0
+ state_vec = quantum_info.Statevector.from_instruction(qc).data[
+ -(2**num_qubits) :
+ ]
+
+ # Remove the factor from the hadamards
+ state_vec *= np.sqrt(2) ** num_qubits
+
+ # Analytic value
+ exact = []
+ for i in range(0, 2**num_qubits):
+ if i == 0:
+ exact.append(0)
+ else:
+ if neg_vals:
+ exact.append(-scaling / (1 - i / (2**num_qubits)))
+ else:
+ exact.append(scaling * (2**num_qubits) / i)
+
+ np.testing.assert_array_almost_equal(state_vec, exact, decimal=2)
+
+
+@ddt
+class TestLinearSolver(unittest.TestCase):
+ """Tests based on the linear solvers classes.
+
+ This class tests
+ * the constructed circuits
+ """
+
+ @idata(
+ [
+ [
+ TridiagonalToeplitz(2, 1, 1 / 3, trotter_steps=2),
+ [1.0, -2.1, 3.2, -4.3],
+ MatrixFunctional(1, 1 / 2),
+ ],
+ [
+ np.array(
+ [
+ [0, 0, 1.585, 0],
+ [0, 0, -0.585, 1],
+ [1.585, -0.585, 0, 0],
+ [0, 1, 0, 0],
+ ]
+ ),
+ [1.0, 0, 0, 0],
+ MatrixFunctional(1, 1 / 2),
+ ],
+ [
+ [
+ [1 / 2, 1 / 6, 0, 0],
+ [1 / 6, 1 / 2, 1 / 6, 0],
+ [0, 1 / 6, 1 / 2, 1 / 6],
+ [0, 0, 1 / 6, 1 / 2],
+ ],
+ [1.0, -2.1, 3.2, -4.3],
+ MatrixFunctional(1, 1 / 2),
+ ],
+ [
+ np.array([[82, 34], [34, 58]]),
+ np.array([[1], [0]]),
+ AbsoluteAverage(),
+ 3,
+ ],
+ [
+ TridiagonalToeplitz(3, 1, -1 / 2, trotter_steps=2),
+ [-9 / 4, -0.3, 8 / 7, 10, -5, 11.1, 13 / 11, -27 / 12],
+ AbsoluteAverage(),
+ ],
+ ]
+ )
+ @unpack
+ def test_hhl(self, matrix, right_hand_side, observable, decimal=1):
+ """Test the HHL class."""
+ if isinstance(matrix, QuantumCircuit):
+ num_qubits = matrix.num_state_qubits
+ elif isinstance(matrix, (np.ndarray)):
+ num_qubits = int(np.log2(matrix.shape[0]))
+ elif isinstance(matrix, list):
+ num_qubits = int(np.log2(len(matrix)))
+
+ rhs = right_hand_side / np.linalg.norm(right_hand_side)
+
+ # Initial state circuit
+ qc = QuantumCircuit(num_qubits)
+ qc.isometry(rhs, list(range(num_qubits)), None)
+
+ hhl = HHL()
+ solution = hhl.solve(matrix, qc, observable)
+ approx_result = solution.observable
+
+ # Calculate analytical value
+ if isinstance(matrix, QuantumCircuit):
+ exact_x = np.dot(np.linalg.inv(matrix.matrix), rhs)
+ elif isinstance(matrix, (list, np.ndarray)):
+ if isinstance(matrix, list):
+ matrix = np.array(matrix)
+ exact_x = np.dot(np.linalg.inv(matrix), rhs)
+ exact_result = observable.evaluate_classically(exact_x)
+
+ np.testing.assert_almost_equal(approx_result, exact_result, decimal=decimal)
+
+ def test_hhl_qi(self):
+ """Test the HHL quantum instance getter and setter."""
+ hhl = HHL()
+ self.assertIsNone(hhl.quantum_instance) # Defaults to None
+
+ # First set a valid quantum instance and check via getter
+ qinst = QuantumInstance(backend=BasicAer.get_backend("qasm_simulator"))
+ hhl.quantum_instance = qinst
+ self.assertEqual(hhl.quantum_instance, qinst)
+
+ # Now set quantum instance back to None and check via getter
+ hhl.quantum_instance = None
+ self.assertIsNone(hhl.quantum_instance)
+
+
+if __name__ == "__main__":
+ unittest.main()