Using the HHL algorithm to compute $A |b \rangle$ instead of $A^{-1} |b \rangle$

In the paper "Compiling basic linear algebra subroutines for quantum computers" here, they discuss (page 2 bottom right) using the HHL algorithm to multiply a vector by a matrix. However, after having read the HHL09 paper, what is being estimated is the state $$|x \rangle = A^{-1} |b \rangle$$, where $$A |x \rangle = |b \rangle$$ for some Hermitian matrix $$A$$ of dimension $$n$$, and vectors $$b, x$$ of size $$n$$, where $$b$$ and $$A$$ are known.

Is it possible to use the HHL algorithm to compute $$A |b \rangle$$ knowing $$A$$ and $$b$$. For instance, running the algorithm in reverse seems not possible (?), because not all steps in the algorithm are unitary.

Note: the graphics have been generated with the LaTeX code available here. Credits to @Niel de Beaudrap.

Yes it is possible!

The HHL algorithm can be schematically depicted as

Let's split down the parts:

1. The first part aims at computing an approximation of the eigenvalues of $$H$$, $$H = A$$ if $$A$$ is hermitian, else $$H = \begin{pmatrix} 0 & A \\ A^\dagger & 0 \end{pmatrix}$$.
2. The second part invert these eigenvalues and adds a phase proportional to these eigenvalues in front of your state (see the original HHL09 paper for the exact expression).
3. The third part uncompute the first part, leaving you with a state like $$\sum_{j=1}^{N}\beta_j \vert u_j \rangle \left( \sqrt{1 - \frac{C^2}{\lambda_j^2}} \vert 0 \rangle + \frac{C}{\lambda_j}\vert 1 \rangle \right).$$
4. The measure project the resulting state. If the measure succeed (measured value is $$1$$), then this means that your state is now $$\sum_{j=1}^{N}\beta_j \vert u_j \rangle \frac{C}{\lambda_j}\vert 1 \rangle = \left(\sum_{j=1}^{N}C\frac{\beta_j}{\lambda_j} \vert u_j \rangle \right) \vert 1 \rangle$$ which is the solution of the linear system considered.

Everything apart being equal, if you do not invert the eigenvalues in step 2, you will end up with a state like $$\sum_{j=1}^{N}\beta_j \vert u_j \rangle \left( \sqrt{1 - C^2\lambda_j^2} \vert 0 \rangle + C\lambda_j\vert 1 \rangle \right).$$ for step 3 and a result like $$\left(C\sum_{j=1}^{N}\beta_j\lambda_j \vert u_j \rangle \right) \vert 1 \rangle$$ if your measurement succeeded.

This is the state you are looking for: the result of the operation $$A\times b$$, scaled by $$C$$ in order to be of unit-norm.

• Isn't this rather circular since you're assuming you have $U$? Is it not the case that $U|b\rangle=A|b\rangle$? So it's a basic assumption that you can already compute that product. – DaftWullie Jun 18 at 8:08
• Or is $U=e^{iHt}$? You could still extract the action fairly easily with measurement on an ancilla, up to some higher order terms. – DaftWullie Jun 18 at 8:11
• $U = e^{iHt}$ for different values of $t$ (as I said, the drawing is schematic, in the real algorithm you apply controlled $U$, $U^2$, ...) so you are right, $U\vert b \rangle = e^{iHt} \vert b \rangle = C\sum_{j=1}^{N}\beta_je^{i\lambda_j t} \vert u_j \rangle$. We don't have the result, but something close. The phase estimation step is here to change this $e^{i\lambda_jt}$ in $\lambda_j$. – Nelimee Jun 18 at 8:32
• But why not simply apply $e^{iX\otimes H\delta t}$ on $|0\rangle|b\rangle$? If $\delta t$ is small, and you find the first qubit to be in the $|1\rangle$ state, you've got $H|b\rangle$ with high probability. – DaftWullie Jun 18 at 10:40
• The standard Pauli X matrix – DaftWullie Jun 18 at 18:19

Let's assume that you have a Hermitian matrix $$H=\left(\begin{array}{cc} 0 & A^\dagger \\ A & 0 \end{array}\right).$$ Let $$|b\rangle$$ be the state that we want to apply $$A$$ to, extended to work on the space that $$H$$ acts on. So, our aim is to implement $$H|b\rangle$$.

Let $$X$$ be the standard Pauli $$X$$ matrix. If we implement a unitary evolution $$U=e^{-i X\otimes H\delta t},$$ then the action of this on an input state $$|0\rangle|b\rangle$$ can be expressed through a Taylor expansion, assuming $$\delta t\|H\|\ll 1$$: $$U|0\rangle|b\rangle=|0\rangle|b\rangle-i\delta t |1\rangle(H|b\rangle)-\frac{\delta t^2}{2}|0\rangle(H^2|b\rangle)+\ldots.$$ So, if we measure the ancilla qubit in the standard basis, then with probability $$\delta t^2\|H|b\rangle\|^2$$, we find the ancilla to be in state $$|1\rangle$$, and the other register is in the desired state $$H|b\rangle$$, up to an accuracy $$O(\delta t^2)$$ (coming from the $$O(\delta t^4)$$ term in the Taylor expansion). If you don't get the $$|1\rangle$$ answer to your measurement, then to accuracy $$O(\delta t^2)$$, you've stil got the original state $$|0\rangle|b\rangle$$, so you can just repeat until success.

By that simple method, the overall error is very bad - you have to repeat $$O(1/\delta t^2)$$ times so that the overall error becomes $$O(1)$$. I'm sure this could be improved by replacing the $$X$$ with a cycle operation on a larger dimensional Hilbert space. However, note that this method is very similar to part of the HHL algorithm (where they do the controlled-rotation from the phase estimation register onto an ancilla) and, at least there, it is claimed that it can be improved using amplitude amplification. So I assume (without having looked at the details) that something similar is possible here.

• Ok, did not think about the Taylor expansion of the exponential, you are right it is simpler to implement than HHL. – Nelimee Jun 20 at 7:02