I have been trying to get my head around the famous(?) paper Quantum algorithm for linear systems of equations (Harrow, Hassidim & Lloyd, 2009) (more popularly known as the HHL09 algorithm paper) for some time, now.
On the very first page, they say:
We sketch here the basic idea of our algorithm and then discuss it in more detail in the next section. Given a Hermitian $N\times N$ matrix $A$, and a unit vector $\vec{b}$, suppose we would like to find $\vec{x}$ satisfying $A\vec{x} = \vec{b}$. (We discuss later questions of efficiency as well as how the assumptions we have made about $A$ and $\vec{b}$ can be relaxed.) First, the algorithm represents $\vec{b}$ as a quantum state $|b\rangle = \sum_{i=1}^{N}b_i|i\rangle$. Next, we use techniques of Hamiltonian simulation [3, 4] to apply $e^{iAt}$ to $|b_i\rangle$ for a superposition of different times $t$. This ability to exponentiate $A$ translates, via the well-known technique of phase-estimation [5–7] into the ability to decompose $|b\rangle$ in the eigenbasis of $A$ and to find the corresponding eigenvalues $\lambda_j$ Informally, the state of the system after this stage is close to $\sum_{j=1}^{j=N} \beta_j |u_j\rangle|\lambda_j\rangle$, where $u_j$ is the eigenvector basis of $A$ and $|b\rangle = \sum_{j=1}^{j=N} \beta_j|u_j\rangle$.
So far so good. As described in Nielsen & Chuang in the chapter "The quantum Fourier transform and its applications", the phase estimation algorithm is used to estimate $\varphi$ in $e^{i2\pi \varphi}$ which is the eigenvalue corresponding to an eigenvector $|u\rangle$ of the unitary operator $U$.
Here's the relevant portion from Nielsen & Chuang:
The phase estimation algorithm uses two registers. The first register contains $t$ qubits initially in the state $|0\rangle$. How we choose $t$ depends on two things: the number of digits of accuracy we wish to have in our estimate for $\varphi$, and with what probability we wish the phase estimation procedure to be successful. The dependence of $t$ on these quantities emerges naturally from the following analysis.
The second register begins in the state $|u\rangle$ and contains as many qubits as is necessary to store $|u\rangle$. Phase estimation is performed in two stages. First, we apply the circuit shown in Figure 5.2. The circuit begins by applying a Hadamard transform to the first register, followed by application of controlled - $U$ operations on the second register, with $U$ raised to successive powers of two. The final state of the first register is easily seen to be:
$$\frac{1}{2^{t/2}}\left(|0\rangle+\text{exp}(2\pi i 2^{t-1}\varphi)|1\rangle)(|0\rangle+\text{exp}(2\pi i 2^{t-2}\varphi)|1\rangle)...(|0\rangle+\text{exp}(2\pi i 2^{0}\varphi)|1\rangle\right)= \frac{1}{2^{t/2}}\sum_{k=0}^{2^{t}-1}\text{exp}(2\pi i \varphi k)|k\rangle$$
The second stage of phase estimation is to apply the inverse quantum Fourier transform on the first register. This is obtained by reversing the circuit for the quantum Fourier transform in the previous section (Exercise 5.5) and can be done in $\Theta (t^2)$ steps. The third and final stage of phase estimation is to read out the state of the first register by doing a measurement in the computational basis. We will show that this provides a pretty good estimate of $\varphi$. An overall schematic of the algorithm is shown in Figure 5.3.
To sharpen our intuition as to why phase estimation works, suppose $\varphi$ may be expressed exactly int bits, as $\varphi = 0.\varphi_1 ... \varphi_t$. Then the state (5.20) resulting from the first stage of phase estimation may be rewritten
$$\frac{1}{2^{t/2}}(|0\rangle + \exp(2\pi i 0.\varphi_t|1\rangle)(|0\rangle + \exp(2\pi i 0.\varphi_{t-1}\varphi_t|1\rangle)...(|0\rangle + \exp(2\pi i 0.\varphi_1...\varphi_t|1\rangle)$$
The second stage of phase estimation is to apply the inverse quantum Fourier transform. But comparing the previous equation with the product form for the Fourier transform, Equation (5.4), we see that the output state from the second stage is the product state $|\varphi_1 ...\varphi_t\rangle$. A measurement in the computational basis, therefore, gives us $\varphi$ exactly!
Summarizing, the phase estimation algorithm allows one to estimate the phase $\varphi$ of an eigenvalue of a unitary operator $U$, given the corresponding eigenvector $|u\rangle$. An essential feature at the heart of this procedure is the ability of the inverse Fourier transform to perform the transformation
$$\frac{1}{2^{t/2}}\sum_{j = 0}^{2^t-1}\exp(2\pi i \varphi j)|j\rangle |u\rangle \to |\tilde \varphi \rangle |u\rangle$$
Let's proceed from here. I found a nice circuit diagram for the HHL09 algorithm here[$\dagger$]:
Step 1 (Phase Estimation):
In the first step of the HHL09 algorithm the same concept (of the standard Quantum Phase Estimation algorithm as described in Nielsen and Chuang) is used. However, we must keep in mind that $A$ by itself isn't a unitary operator. However, if we assume that $A$ is Hermitian then the exponential $e^{iAt}$ is unitary (no worries, there's exists a workaround in case $A$ isn't Hermitian!).
Here, we can write $U=e^{iAt}$. There's another subtle point involved here. We do not know the eigenvectors $|u_j\rangle$ of $U$ beforehand (but we do know that for any unitary matrix of size $N\times N$ there exist $N$ orthonormal eigenvectors). Moreover, we need to remind ourselves that if the eigenvalues of $A$ are $\lambda_j$ then the eigenvalues of $e^{iAt}$ will be $e^{i \lambda_j t}$. If we compare this with the form of eigenvalues given in Nielsen and Chuang for $U$ i.e. if $e^{2\pi i \varphi} \equiv e^{ i \lambda_j t}$, we'd find $\varphi = \frac{\lambda_j t}{2\pi}$. In this case, we begin in the state $|b\rangle$ (which can be written as a superposition of the eigenvectors of $U$ i.e. $\sum_{j=1}^{j=N}\beta_j|u_j\rangle$) rather than any particular eigenvector $|u_j\rangle$ of $U$, as far as the second register of qubits is concerned. If we had begun in the state $|u\rangle \otimes (|0\rangle)^{\otimes t}$ we would have ended up with $|u\rangle \otimes |\tilde\varphi\rangle$ i.e. $|u_j\rangle \otimes |\tilde{\frac{\lambda_j t}{2\pi}}\rangle$ (considering that $\lambda_j$ is the eigenvalue associated with the eigenvector $|u_j\rangle$ of $A$). Now, instead if we begin in the superposition of eigenvectors $\sum_{j=1}^{j=N}\beta_j|u_j\rangle$ we should end up with $\sum_{j=1}^{j=N}\beta_j|u_j\rangle\otimes |\tilde{\frac{\lambda_j t}{2\pi}}\rangle$.
Question:
Part 1: In the HHL09 paper, they wrote about the state of the system after this Phase Estimation step is $\sum_{j=1}^{j=N}\beta_j|u_j\rangle\otimes |\tilde\lambda_j\rangle$. However, from what I wrote above it seems to me that the state of the system should rather be $\sum_{j=1}^{j=N}\beta_j|u_j\rangle\otimes |\tilde{\frac{\lambda_j t}{2\pi}}\rangle$.
What am I missing here? Where did the factor of $\frac{t}{2\pi}$ vanish in their algorithm?
Edit: Part 2 has been asked here to make the individual questions more focused.
I also have several confusions regarding Step 2 and Step 3 of the HHL09 algorithm too, but I decided to post them as separate question threads, as this one is becoming too long. I'll add the links to those question threads, on this post, once they are created.
[$\dagger$]: Homomorphic Encryption Experiments on IBM's Cloud Quantum Computing Platform Huang et al. (2016)
