# Quantum algorithm for linear systems of equations (HHL09): Step 2 - What is $|\Psi_0\rangle$?

In the paper: Quantum algorithm for linear systems of equations (Harrow, Hassidim & Lloyd, 2009), what's written up to the portion

The next step is to decompose $|b\rangle$ in the eigenvector basis, using phase estimation [5–7]. Denote by $|u_j\rangle$ the eigenvectors of $A$ (or equivalently, of $e^{iAt}$), and by $\lambda_j$ the corresponding eigenvalues.

on page $2$ makes some sense to me (the confusions up till there have been addressed in the previous posts linked above). However, the next portion i.e. the $R(\lambda^{-1})$ rotation seems a bit cryptic.

Let $$|\Psi_0\rangle := \sqrt{\frac{2}{T}}\sum_{\tau =0}^{T-1} \sin \frac{\pi(\tau+\frac{1}{2})}{T}|\tau\rangle$$

for some large $T$. The coefficients of $|\Psi_0\rangle$ are chosen (following [5-7]) to minimize a certain quadratic loss function which appears in our error analysis (see  for details).

Next, we apply the conditional Hamiltonian evolution $\sum_{\tau = 0}^{T-1}|\tau\rangle \langle \tau|^{C}\otimes e^{iA\tau t_0/T}$ on $|\Psi_0\rangle^{C}\otimes |b\rangle$, where $t_0 = \mathcal{O}(\kappa/\epsilon)$.

## Questions:

1. What exactly is $|\Psi_0\rangle$? What do $T$ and $\tau$ stand for? I've no idea from where this gigantic expression $$\sqrt{\frac{2}{T}}\sum_{\tau =0}^{T-1} \sin \frac{\pi(\tau+\frac{1}{2})}{T}|\tau\rangle$$ suddenly comes from and what its use is.

2. After the phase estimation step, the state of our system is apparently:

$$\left(\sum_{j=1}^{j=N}\beta_j|u_j\rangle\otimes |\tilde\lambda_j\rangle\right)\otimes |0\rangle_{\text{ancilla}}$$

This surely cannot be written as $$\left(\sum_{j=1}^{j=N}\beta_j|u_j\rangle\right)\otimes \left(\sum_{j=1}^{j=N}|\tilde\lambda_j\rangle\right)\otimes |0\rangle_{\text{ancilla}}$$ i.e.

$$|b\rangle\otimes \left(\sum_{j=1}^{j=N}|\tilde\lambda_j\rangle\right)\otimes |0\rangle_{\text{ancilla}}$$

So, it is clear that $|b\rangle$ is not available separately in the second register. So I've no idea how they're preparing a state like $|\Psi_0\rangle^{C}\otimes |b\rangle$ in the first place! Also, what does that $C$ in the superscript of $|\Psi_0\rangle^{C}$ denote?

3. Where does this expression $\sum_{\tau = 0}^{T-1}|\tau\rangle \langle \tau|^{C}\otimes e^{iA\tau t_0/T}$ suddenly appear from ? What's the use of simulating it? And what is $\kappa$ in $\mathcal{O}(\kappa/\epsilon)$ ?

# 1. Definitions

Names and symbols used in this answer follow the ones defined in Quantum linear systems algorithms: a primer (Dervovic, Herbster, Mountney, Severini, Usher & Wossnig, 2018). A recall is done below.

## 1.1 Register names

Register names are defined in Figure 5. of Quantum linear systems algorithms: a primer (Dervovic, Herbster, Mountney, Severini, Usher & Wossnig, 2018) (reproduced below):

• $S$ (1 qubit) is the ancilla register used to check if the output is valid or not.
• $C$ ($n$ qubits) is the clock register, i.e. the register used to estimate the eigenvalues of the hamiltonian with quantum phase estimation (QPE).
• $I$ ($m$ qubits) is the register storing the right-hand side of the equation $Ax = b$. It stores $x$, the result of the equation, when $S$ is measured to be $\left|1\right>$ at the end of the algorithm. # 2. About $\left|\Psi_0\right>$:

1. What exactly is $\left|\Psi_0\right>$?

$\left|\Psi_0\right>$ is one possible initial state of the clock register $C$.

2. What do $T$ and $\tau$ stand for?

$T$ stands for a big positive integer. This $T$ should be as large as possible because the expression of $\left|\Psi_0\right>$ asymptotically minimise a given error for $T$ growing to infinity. In the expression of $\left|\Psi_0\right>$, $T$ will be $2^n$, the number of possible states for the quantum clock $C$.

$\tau$ is just the summation index

3. Why such a gigantic expression for $\left|\Psi_0\right>$?

See DaftWullie's post for a detailed explanation.

Following the citations in Quantum algorithm for linear systems of equations (Harrow, Hassidim & Lloyd, 2009 v3) we end up with:

1. The previous version of the same paper Quantum algorithm for linear systems of equations (Harrow, Hassidim & Lloyd, 2009 v2). The authors revised the paper 2 times (there are 3 versions of the original HHL paper) and version n°3 does not include all the informations provided in the previous versions. In the V2 (section A.3. starting at page 17), the authors provide a detailed analysis of the error with this special initial state.
2. Optimal Quantum Clocks (Buzek, Derka, Massar, 1998) where the expression of $\left|\Psi_0\right>$ is given as $\left|\Psi_{opt}\right>$ in Equation 10. I don't have the knowledge to understand fully this part, but it seems like this expression is "optimal" in some sense.

# 3. Preparation of $\left|\Psi_0\right>$:

As said in the previous part, $\left|\Psi_0\right>$ is an initial state. They do not prepare $\left|\Psi_0\right>$ after the phase estimation procedure. The sentence ordering is not really optimal in the paper. The phase estimation procedure they use in the paper is a little bit different from the "classic" phase estimation algorithm represented in the quantum circuit linked in part 1, and that is why they explain it in details.

Their phase estimation algorithm is:

1. Prepare the $\left|\Psi_0\right>$ state in the register $C$.
2. Apply the conditional Hamiltonian evolution to the registers $C$ and $I$ (which are in the state $\left|\Psi_0\right>\otimes \left|b\right>$).
3. Apply the quantum Fourier transform to the resulting state.

Finally, the $C$ in $\left| \Psi_0 \right>^C$ means that the state $\left| \Psi_0 \right>$ is stored in the register $C$. This is a short and convenient notation to keep track of the registers used.

# 4. Hamiltonian simulation:

First of all, $\kappa$ is the condition number (Wikipedia page on "condition number") of the matrix $A$.

$\sum_{\tau = 0}^{T-1}|\tau\rangle \langle \tau|^{C}\otimes e^{iA\tau t_0/T}$ is the mathematical representation of a quantum gate.

The first part in the sum $|\tau\rangle \langle \tau|^{C}$ is a control part. It means that the operation will be controlled by the state of the first quantum register (the register $C$ as the exponent tells us).

The second part is the "Hamiltonian simulation" gate, i.e. a quantum gate that will apply the unitary matrix given by $e^{iA\tau t_0/T}$ to the second register (the register $I$ that is in the initial state $\left|b\right>$).

The whole sum is the mathematical representation of the controlled-U operation in the quantum circuit of "1. Definitions", with $U = e^{iA\tau t_0/T}$.

$\newcommand{\bra}{\left\langle#1\right|}\newcommand{\ket}{\left|#1\right\rangle}\newcommand{\proj}{|#1\rangle\langle#1|}\newcommand{\half}{\frac12}$In answer to your first question, I wrote myself some notes some time ago about my understanding of how it worked. The notation is probably a bit different (I've tried to bring it more into line, but it's easy to miss bits), but attempts to explain that choice of the state $|\Psi_0\rangle$. There also seem to be some factors of $\frac12$ floating around in places.

When we first study phase estimation, we're usually thinking about it in respect to use in some particular algorithm, such as Shor's algorithm. This has a specific goal: getting the best $t$-bit approximation to the eigenvalue. You either do, or you don't, and the description of phase estimation is specifically tuned to give as high a success probability as possible.

In HHL, we are trying to produce some state $$\ket{\phi}=\sum_j\frac{\beta_j}{\lambda_j}\ket{\lambda_j},$$ where $\ket{b}=\sum_j\beta_j\ket{\lambda_j}$, making use of phase estimation. The accuracy of the approximation of this will depend far more critically on an accurate estimation of the eigenvalues that are close to 0 rather than those that are far from 0. An obvious step therefore, is to attempt to modify the phase estimation protocol so that rather than using `bins' of fixed width $2\pi/T$ for approximating the phases of $e^{-iAt}$ ($T=2^t$ and $t$ is number of qubits in phase estimation register), we might rather specify a set of $\phi_y$ for $y\in\{0,1\}^t$ to act as the centre of each bin so that we can have vastly increased accuracy close to 0 phase. More generally, you might specify a trade-off function for how tolerant you might be of errors as a function of the phase $\phi$. The precise nature of this function can then be tuned to a given application, and the particular figure of merit which you will use to determine success. In the case of Shor's algorithm, our figure of merit was simply this binning protocol -- we were successful if the answer was in the correct bin, and unsuccessful outside it. This is not going to be the case in HHL, whose success is more reasonably captured by a continuous measure such as the fidelity. So, for the general case, we shall designate a cost function $C(\phi,\phi')$ which specifies a penalty for answers $\phi'$ if the true phase is $\phi$.

Recall that the standard phase estimation protocol worked by producing an input state that was the uniform superposition of all basis states $\ket{x}$ for $x\in\{0,1\}^t$. This state was used to control the sequential application of multiple controlled-$U$ gates, which are followed up by an inverse Fourier transform. Imagine we could replace the input state with some other state $$\ket{\Psi_0}=\sum_{x\in\{0,1\}^t}\alpha_x\ket{x},$$ and then the rest of the protocol could work as before. For now, we will ignore the question of how hard it is to produce the new state $\ket{\Psi_0}$, as we are just trying to convey the basic concept. Starting from this state, the use of the controlled-$U$ gates (targeting an eigenvector of $U$ of eigenvalue $\phi$), produces the state $$\sum_{x\in\{0,1\}^t}\alpha_xe^{i\phi x}\ket{x}.$$ Applying the inverse Fourier transform yields $$\frac{1}{\sqrt{T}}\sum_{x,y\in\{0,1\}^t}e^{ix\left(\phi-\frac{2\pi y}{M}\right)}\alpha_x\ket{y}.$$ The probability of getting an answer $y$ (i.e. $\phi'=2\pi y/T$) is $$\frac{1}{T}\left|\sum_{x\in\{0,1\}^t}e^{ix\left(\phi-\frac{2\pi y}{T}\right)}\alpha_x\right|^2$$ so the expected value of the cost function, assuming a random distribution of the $\phi$, is $$\bar C=\frac{1}{2\pi T}\int_0^{2\pi}d\phi\sum_{y\in\{0,1\}^t}\left|\sum_{x\in\{0,1\}^t}e^{ix\left(\phi-\frac{2\pi y}{T}\right)}\alpha_x\right|^2C(\phi,2\pi y/T),$$ and our task is to select the amplitudes $\alpha_x$ that minimise this for any specific realisation of $C(\phi,\phi')$. If we make the simplifying assumption that $C(\phi,\phi')$ is only a function of $\phi-\phi'$, then we can make a change of variable in the integration to give $$\bar C=\frac{1}{2\pi}\int_0^{2\pi}d\phi\left|\sum_{x\in\{0,1\}^t}e^{ix\phi}\alpha_x\right|^2C(\phi),$$ As we noted, the most useful measure is likely to be a fidelity measure. Consider we have a state $\ket{+}$ and we wish to implement the unitary $U_\phi=\proj{0}+e^{i\phi}\proj{1}$, but instead we implement $U_{\phi'}=\proj{0}+e^{i\phi'}\proj{1}$. The fidelity measures how well this achieves the desired task, $$F=\left|\bra{+}U_{\phi'}^\dagger U\ket{+}\right|^2=\cos^2\left(\frac{\phi-\phi'}{2}\right),$$ so we take $$C(\phi-\phi')=\sin^2\left(\frac{\phi-\phi'}{2}\right),$$ since in the ideal case $F=1$, so the error, which is what we want to minimise, can be taken as $1-F$. This will certainly be the correct function for evaluating any $U^t$, but for the more general task of modifying the amplitudes, not just the phases, the effects of inaccuracies propagate through the protocol in a less trivial manner, so it is difficult to prove optimality, although the function $C(\phi-\phi')$ will already provide some improvement over the uniform superposition of states. Proceeding with this form, we have $$\bar C=\frac{1}{2\pi}\int_0^{2\pi}d\phi\left|\sum_{x\in\{0,1\}^t}e^{ix\phi}\alpha_x\right|^2\sin^2\left(\half\phi\right),$$ The integral over $\phi$ can now be performed, so we want to minimise the function $$\half\sum_{x,y=0}^{T-1}\alpha_x\alpha_y^\star(\delta_{x,y}-\half\delta_{x,y-1}-\half\delta_{x,y+1}).$$ This can be succinctly expressed as $$\min\bra{\Psi_0}H\ket{\Psi_0}$$ where $$H=\half\sum_{x,y=0}^{T-1}(\delta_{x,y}-\half\delta_{x,y-1}-\half\delta_{x,y+1})\ket{x}\bra{y}.$$ The optimal choice of $\ket{\Psi_0}$ is the minimum eigenvector of the matrix $H$, $$\alpha_x=\sqrt{\frac{2}{T+1}}\sin\left(\frac{(x+1)\pi}{T+1}\right),$$ and $\bar C$ is the minimum eigenvalue $$\bar C=\half-\half\cos\left(\frac{\pi}{T+1}\right).$$ Crucially, for large $T$, $\bar C$ scales as $1/T^2$ rather than the $1/T$ that we would have got from the uniform coupling choice $\alpha_x=1/\sqrt{T}$. This yields a significant benefit for the error analysis.

If you want to get the same $|\Psi_0\rangle$ as reported in the HHL paper, I believe you have to add the terms $-\frac14\left(\ket{0}\bra{T-1}+\ket{T-1}\bra{0}\right)$ to the Hamiltonian. I have no justification for doing so, however, but this is probably my failing.