# Unveiling the Theoretical Pillars and Cutting-Edge Applications of Quantum Amplitude Estimation

*Editor: Wang, Xuan-Wei*
*Time: 2024.06.01*

## Abstract

In the last semester, we identified a critical issue in finance known as the "option pricing problem." We constructed a quantum circuit framework to solve the numerical integration problem.

However, when generating the computational option price, we need to verify and reduce the discrepancy between the computed option price and the actual one. Therefore, we may use Iterative Quantum Amplitude Estimation after completing the main quantum circuit.

In this semester's project, we elucidate the significance and primary challenges associated with the call option pricing problem. Specifically, we delve into the utilization of quantum amplitude estimation techniques to solve the convergence problem when addressing the PDE model, such as the Black-Scholes model.

## 1. Grover's Algorithm

Grover's algorithm is an important quantum algorithm for solving unstructured search problem. Comparing to classical algorithm, it can only use few quantum gate to search through an unsorted database in approximately $\sqrt{N}$ steps.

:::warning
**Search problem:**
*Input:* $f:\Sigma^n\rightarrow\Sigma,\ \text{where}\ \Sigma=\{0,1\}$
*Output:* a string $x\in\Sigma^n$ satisfying $f(x)=1$, or "no solution" if no such strings exist.
:::
※ $f:\{0,1,...,N-1\}\mapsto\{0,1\}$

Here, we have:
* $N=2^n$ : size of database
* $\mathcal{A}_0=\{x\in\Sigma^n:f(x)=0\}$: set of non-solutions, where $$|\mathcal{A}_0\rangle=\frac{1}{\sqrt{|\mathcal{A}_0}|}\sum_{x\in\mathcal{A}_0}|x\rangle$$

* $\mathcal{A}_1=\{x\in\Sigma^n:f(x)=1\}$: set of solutions, where $$|\mathcal{A}_1\rangle=\frac{1}{\sqrt{|\mathcal{A}_1}|}\sum_{x\in\mathcal{A}_1}|x\rangle$$

For deterministic algorithm, we can iterate through all the bits to find the value such that $f(x)=1$, which requires $O(N)$ queries.

Probabilistic algorithms offer minor improvements, but still require a number of quires linear in $N$.

Grover's algorithm is a quantum algorithm for search requiring $O(\sqrt{N})$ queries.

:::danger
***Grover's Algorithm***
1. *Initialize:* set $n$ qubits to the state $\mathcal{H}^{\otimes n}|0^n\rangle$
2. *Iterate:* apply the ==Grover's operator== $t$ times
3. *Measure:* a standard basis measurement yields a candidate solution.

**- - - - - - - -**

==Grover's Operator== $$G=\mathcal{H}^{\otimes n}\ Z_{OR}\ \mathcal{H}^{\otimes n}\ Z_f$$
:::
Here, for all $x\in\Sigma^n$, the operator $Z_{OR}$ and $Z_f$ are defined as: $$Z_{OR}\ |x\rangle=\begin{cases}
|x\rangle & \quad \text{if } x=0^n\\
-|x\rangle & \quad \text{if } x\neq0^n
\end{cases}$$ $$Z_f:\ |x\rangle\mapsto(-1)^{f(x)}|x\rangle$$

Next, we consider the action of the Grover's operator $G$ on the set of non-solutions and solutions.

We set the first initialized state as $$|u\rangle=\mathcal{H}^{\otimes n}|0^n\rangle=\frac{1}{\sqrt{N}}\sum_{x\in\Sigma^n}|x\rangle$$
Then this state is contained in the subspace spanned by $\mathcal{A}_0,\ \mathcal{A}_1$: $$|u\rangle=\sqrt{\frac{|\mathcal{A}_0|}{N}}|\mathcal{A}_0\rangle+\sqrt{\frac{|\mathcal{A_1}|}{N}}|\mathcal{A}_1\rangle$$

For the operator $Z_f$ : $$Z_f|\mathcal{A}_0\rangle=|\mathcal{A}_0\rangle$$$$Z_f|\mathcal{A}_1\rangle=-|\mathcal{A}_1\rangle$$
And we define $Z_{OR}\equiv 2|0^n\rangle\langle0^n|-\mathbb{1}$, then we have $$\mathcal{H}^{\otimes n}\ Z_{OR}\ \mathcal{H}^{\otimes n}=\mathcal{H}^{\otimes n}\ (2|0^n\rangle\langle0^n|-\mathbb{1})\ \mathcal{H}^{\otimes n}=2|u\rangle\langle u|-\mathbb{1}$$

Hence, the action of $G$ on $\text{span}\{|\mathcal{A}_0\rangle,|\mathcal{A}_0\rangle\}$ can be described by a $2\times2$ matrix:$$M=\begin{pmatrix}
\frac{|\mathcal{A}_0|-|\mathcal{A}_1|}{N} & -\frac{2\sqrt{|\mathcal{A}_0|\cdot|\mathcal{A}_1|}}{N} \\
-\frac{2\sqrt{|\mathcal{A}_0|\cdot|\mathcal{A}_1|}}{N} & \frac{|\mathcal{A}_0|-|\mathcal{A}_1|}{N} \\
\end{pmatrix}=\begin{pmatrix}
\sqrt{\frac{|\mathcal{A}_0|}{N}} & -\sqrt{\frac{|\mathcal{A}_1|}{N}} \\
-\sqrt{\frac{|\mathcal{A}_1|}{N}} & \sqrt{\frac{|\mathcal{A}_0|}{N}} \\
\end{pmatrix}^2$$
Note that we can transfer this matrix into a form of rotation matrix, we have $$M=\begin{pmatrix}
\cos(2\theta) & \sin(2\theta) \\
-\sin(2\theta) & \cos(2\theta) \\
\end{pmatrix},\ \text{where}\ \theta=\sin^{-1}(\sqrt{\frac{|\mathcal{A}_1|}{N}})$$

That is, for each time we apply $G$, the state is rotated by an angle $2\theta$. Therefore, after iterating $t$ times, we have $$G^t|u\rangle=\cos((2t+1)\theta)|\mathcal{A}_0\rangle+\sin((2t+1)\theta)|\mathcal{A}_1\rangle$$

Let's dive into the visualized interpretation.

First, $Z_f$ is a reflection about the line $L_1$ paralleling to $|\mathcal{A}_0\rangle$.
>

Second, $\mathcal{H}^{\otimes n}\ Z_{OR}\ \mathcal{H}^{\otimes n}$ is also a reflection about the line $L_2$ paralleling to $|u\rangle$.
>

Hence, we combines the two reflection together, and we have:
>

Finally, we need to decide the best choice of iteration time $t$ such that the output angle of the vector is close to $\frac{\pi}{2}$, which means the vector needs to be closer to our target state $|\mathcal{A}_1\rangle$.

Hence, measuring after $t$ iterations gives an outcome $x\in\mathcal{A}_1$ with probability $$\sin^2((2t+1)\theta)\approx\frac{\pi}{2}\Leftrightarrow t\approx\frac{\pi}{4\theta}-\frac{1}{2}\xrightarrow{\text{closest integer}}t=\lfloor\frac{\pi}{4\theta}\rfloor$$

Here are some important considerations:
1. $t$ must be an integer
2. $\theta$ depends on the number of solutions $s=|\mathcal{A}_1|$

Then, one can easily find out that unique search requires $O(\sqrt{N})$ queries, and that multple search requires $O(\sqrt{\frac{N}{s}})$ queries.

## 2. Quantum Amplitude Amplification and Estimation
*(Brassard et al. )*

For QAE problem, we aim to deal with multi-solutions problem; namely, we try to estimate $t=|\{x\in X|f(x)=1\}|$, which means the number of inputs on which $f:X\rightarrow \{0,1\}$ takes the value $1$.

Similarly, given a unitary transformation $\mathcal{A}$, let $|\Psi\rangle=\mathcal{A}|0\rangle$.

Write $|\Psi\rangle=|\Psi_1\rangle+|\Psi_0\rangle$ as the good and bad components of $|\Psi\rangle$. Then we aim to estimate $a=\langle\Psi_1||\Psi_1\rangle$.

We take $X=\{0,1,\cdots, N-1\}$, if $N$ is a power of $2$, then we set $\mathcal{A}\equiv\mathcal{H}$; if not, we then set $\mathcal{A}\equiv F_M$, which represents the quantum Fourier transform that, for evry integer $M\geq 1$, is defined by $$F_M:|x\rangle\longmapsto\frac{1}{\sqrt{M}}\sum^{M-1}_{y=0}e^{2\pi ixy/M}|y\rangle\ \ \ (0\leq x<M)$$

Note that we have $a=t/N$. Thus an estimation for $a$ equals to an estimation for $t$. Here, we may set an operaor $Q$ for estimation of $a$ s.t. $$Q=-\mathcal{A}S_0\mathcal{A}^{-1}S_f$$, where
* $S_0(\phi)=\begin{cases}
e^{i\phi}|x\rangle & \quad \text{if } |x\rangle=|0\rangle\\
|x\rangle & \quad \text{if } |x\rangle\neq|0\rangle
\end{cases}$
* $S_f(\varphi):|x\rangle\longmapsto\begin{cases}
e^{i\varphi}|x\rangle & \quad \text{if } f(x)=1\\
|x\rangle & \quad \text{if } f(x)=0
\end{cases}$

Then the following algorithm computes an estimation for $a$, via an estimate for $\theta_a$:

:::danger
**Algorithm(Est_Amp($\mathcal{A},\chi, M$))**
1. Initialize two registers of appropriate sizes to the state $|0\rangle\mathcal{A}|0\rangle$
2. Apply $F_M$ to the first register
3. Apply $\Lambda_M(Q)$ where $Q=-\mathcal{A} S_0 \mathcal{A} ^{-1} S_f$
4. Apply $F^{-1}_M$ to the first register.
5. Measure the first register and denote the outcome $|y\rangle$
6. Output $\tilde{a}=\sin^2(\pi\frac{y}{M})$

**- - - - - - - -**
This algrithm can also be simplifly described as the unitary transformation $$((F_M^{-1}\otimes\mathbf{I})\Lambda_M(Q)(F_M\otimes\mathbf{I}))$$

:::

The following illustrates the quantum circuit:
>
We can see that there are $m$ ancilla qubits and $n+1$ state qubits.

However, although quantum Fourier transform(QFT) is believed to achieve exponential speedup -- such as in Shor's algorithm for factoring problem -- it will cost the computation complexity to do quantum Fourier transformation.

Hence, we may consider replace QFT by a set of Grover's iterations combined with a Maximum Likelihood Estimation(MLE), and is named as Maximum Likelihood Amplitude Estimation(MLAE).

## 3. Innovative QAE algorithm and experiments

### 1. QAE Simplified (QAES)

>Our algorithm for approximate counting mirrors a standard classical approach for the following problem: given a biased coin that is heads with probability $p$, estimate $p$. First, for $k = 0,1,2,...$ the coin is tossed $2k$ times, and stop at the $k$ where heads was observed at least once. This gives a rough guess for $p$, up to some multiplicative constant. Second, this rough guess is improved to the desired $1 +\epsilon$ approximation via more coin tosses.

Here uses Grover's algorithm to achieve quadratic speedup.

If $K$ out of $N$ items are marked, then define $θ := \arcsin (pK/N)$ to be the **Grover angle**. For any odd integer $r$, Grover’s algorithm lets us prepare a coin that lands heads with probability $p =\sin2(rθ)$, by making $O(r)$ queries to the oracle.

This, in turn, is done by finding some value of $r$ that distinguishes the two possibilities, by making $θ ≈ θ_{min}$ and $θ ≈ θ_{max}$ lead to two nearly-orthogonal quantum states that are easy to distinguish by a measurement.

>

> Theorem 1.
> Let $S\subset [N]$ be a nonempty set of marked item, and let $K=|S|$. Given access to a membership oracle to $S$ and $\epsilon,\ \delta>0$, there exists a quantum algorithm that outputs an estimate $\hat{K}$ satisfies $$K(1-\epsilon)<\hat{K}<K(1+\epsilon)$$ with prob. at least $1-\delta$. There exists a function $Q(N,K,\epsilon,\delta)\in O(\sqrt{N/K}\epsilon^{-1}\log(\delta^{-1}))$ such that the algorithm makes fewer than **$Q(N,K,\epsilon,\delta)$** queries whenever the estimate is accurate. The algorithm needs $O(\log N)$ qubits of space.

>

Lemma 2. provide the possibility and feasibility to find $r$.

> Theorem 3.
> Given $\epsilon,\ \delta>0$ as well as access to an $(n+1)$-qubits unitary $U$ satisfying $$U|0^n\rangle|0\rangle=a|\phi\rangle|0\rangle+\sqrt{1-a^2}|\tilde{\phi}\rangle|1\rangle$$, where $|\phi\rangle$ and $|\tilde{\phi}\rangle$ are arbitary $n$-qubit states and $0<a<1$, there exists an algorithm that outputs an estimate $\hat{a}$ that satisfies $$a(1-\epsilon)<\hat{a}<a(1+\epsilon)$$ and uses $O(a^{-1}\epsilon^{-1}\log(\delta^{-1}))$ applications of $U$ and $U^{\dagger}$ with prob. at least $1-\delta$.

Here, Aaronson et al.(2021) derived a way to complete amplitude estimation. They try to generalize the algorithm for approximate counting to amplified estimation. Here is the algorithm:

>

### 2. Maximum Likelihood Amplified Estimation (MLQAE)
(Suzuki et al.(2020))(on IBM 5.1)

Suzuki et al. provided a way to do amplitude estimation via **Maximum Likelihood Estimation** instead of QFT.

>The key idea of MLAE is post-processing the measurements of the quantum circuits, $Q^k\mathcal{A}$, with different values of $k$.

>

Note that
* $N_k$ : the number of measurements(shots)
* $h_k$ : the number of measuring good states for the state $Q^{m_k}|\Psi\rangle$
* $Q$ : the amplitude amplification operator
* $m_k$ : the number of application of $Q$ to qubits

> 

Fig 1. above shows the process for the whole algorithm:
* Left: Preparing states $Q^{m_k}|\Psi\rangle$, we obtain the number of $h_k$, $\forall k=1,...,M,\ M\in\mathbb{N}$
* Center: Construct likelihood functions $L_k(h_k;\theta_a)$
* Right: Combining likelihood functions to get $L(h;\theta_a)$, we then output $$\hat{\theta_a}=\arg\max_{\theta_a} L(h;\theta_a)=\arg\max_{\theta_a} \ln L(h;\theta_a)$$

We need to designing the sequences $\{m_k,N_k\}$ so that the resulting ML estimate $\hat{\theta_a}$ outperforms the classical limit $O(N_q^{-1/2})$ and hopefully achieves the Heisenberg limit $O(N_q^{-1})$, i.e., the quadratic speedup.

Then, $\forall k=1,...,M$we consider three ways to determine our $\{m_k\}$:
* $m_k=0$: this means we conduct classical random sampling (numerical simulation)
* $m_k=k$ and $N_k=N_{shots}$: Linearly incremental sequence (**LIS**)
* $m_k=2^k$ and $N_k=N_{shots}$: Exponentially incremental sequence (**EIS**)

The following picture represents the experiment results by Suzuki et al.:

>

The figure shows the relationship between number of queries and errors for the target prob. $a=\sin^2(\theta_a)=\frac{2}{3},\ \frac{1}{3},\ \frac{1}{6},\ \frac{1}{12},\ \frac{1}{24},\ \frac{1}{48}$ with $N_{shots}=100$.

From the figure, we can see that the errors by taking **EIS**(black) is less than by taking **LIS**(red) and **numerical estimation**(blue) under the same queries. Moreover, with the number of queries increasing, the errors are monotonically decreasing.

Hence, Suzuki et al. compared the query complexity and computational complexity for the three different choosen of $\{m_k\}$ as following table have shown:
> 

Moreover, we can compare the complexity of the quantum circuit for conventional QAE and MLQAE by Suzuki et al., in the case of $n=2$ with single $Q$ operation:
> 

>

By compare only the numbers of CNOT gates and numbers of qubits, we can easily observe how simple the MLQAE circuit are comparing to conventional quantum algorithm.

Finally, we denote that MLQAE needs $O(\frac{1}{\epsilon}(\sqrt{\frac{K}{N}}))$ oracle queries to estimate $K$ good states in $N$ samples.

### 3. Iterative QAE
(Grinko et al.)

The iteravive QAE conducted by Grinko et al. aims to find a confidence interval $[\theta_l,\theta_u]\subset[0,\pi/2]$, and then choose $\theta_a=\frac{\theta_l+\theta_l}{2}$.

However, they try to estimate $\cos((4k+2)\theta_a)$ instead of estimate $\sin^2((2k+1)\theta_a)$, since cosine function is monotonic in either $[0,\pi]$ and $[\pi,2\pi]$. Thus, they aim to find the largest $k$, namely, the numbers of application of amplitude amplification operator $Q$ to quantum states, such that the scaled interval $[(4k+2)\theta_l,(4k+2)\theta_u]_{\mod2\pi}$ is fully contained in $[0,\pi]$ or $[\pi,2\pi]$.

Finally, once they find out $\theta_l,\ \theta_u$ such that $\theta_u-\theta_l<2\epsilon$ for some $\epsilon>0$, they generate an interval $[a_l,a_u]=[\sin^2(\theta_l),\sin^2(\theta_u)]$. We ultimately satisfy the following issue:

>Let $\mathcal{A}$ be $(n+1)$-qubit quantum circuit satisfying $$\mathcal{A}|0\rangle_{n}|0\rangle=\sqrt{1-a}|\psi_0\rangle_n|0\rangle+\sqrt{a}|\psi_1\rangle_n|1\rangle$$where $a\in[0,1]$ is unknown and $|\psi_0\rangle,\ |\psi_1\rangle$ are two normalized states, and $\mathcal{A}$ can be constructed with $N_{\mathcal{A}}\in\mathbb{N}$ number of elementary gates.
>Then, the following holds:
>The IQAE algorithm outputs a confidence interval $[a_l,a_u]$ that satisfies $$a\notin[a_l,a_u]$$ with probability at most $\alpha$, where $0\leq a_u-a_l\leq 2\epsilon$.
>In particular, the estimator $\hat{a}:=\frac{(a_u-a_l)}{2}$ of a satisfies $$|a-\hat{a}|<\epsilon$$ where $a\in[a_l,a_u]$ with probability $(1-\alpha)$.

For the correctness of IQAE, we have following theorem 1:

> **Theorem 1.(Correctness of IQAE)**
> Suppose a confidence level $1-a\in(0,1)$, a target accuracy $\epsilon>0$, and a number of shots $N_{shots}\in\{1,...,N_{max}(\epsilon,a)\}$, where $$N_{max}(\epsilon,a)=\frac{32}{(1-2\sin(\pi/14))^2}\log(\frac{2}{a}\log_2(\frac{\pi}{4\epsilon}))$$
> In this case, IQAE terminates after a miximum number of $\lceil\log_2(\pi/8\epsilon)\rceil$ rounds, where we define one round as a set of iterations with the same $k_i$, and each round consists of at most $N_{max}(\epsilon,a)/N_{shots}$ iterations. IQAE returns $[a_l,a_u]$ with $a_u-a_l<2\epsilon$ and $$P(a\notin[a_l,a_u])\leq a$$
> Thus, $\tilde{a}=(a_l+a_u)/2$ leads to an estimate for $a$ with $|a-\tilde{a}|\leq \epsilon$ with a confidence of $1-a$.
> Furthermore, for the total number of $Q$-applications, $N_{oracle}$, it holds that $$N_{oracle}<\frac{50}{\epsilon}\log(\frac{2}{a}\log_2(\frac{\pi}{4\epsilon}))$$.

The check work for theorem 1. and algorithm details have been done by Grinko et al.(2021) in their paper. Hence, we aim to demonstrate their comparison experiment results here.

For the correctness of theorem 1, we use **Chernoff-Hoeffding(CH) bound** to generate the application times of operator $Q$, which could be a loose upper bound for the counting time. If we need more presice result, we can use **Clopprt-Pearson's confident inteveral for Bernoulli distribution(CP)** instead, which is said to lower the vonstant overhead in $N_{oracle}$ by a factor of $3$. But it is more complex to analyze analytically.

Here, we put some analysis results from Grinko et al. to demonstrate the different between CH and CP bounding methods:
>
>
>

## Conclusion

Grinko et al. have analyzed the performance for several types of quantum amplitude amplification and estimation algorithm:
>
>

They have shown that MLAE(i.e. MLQAE) and IQAE with both CH bounding method and CP bounding method performed better then conventional QAE, QAES, and classical Monte Carlo simulation.

Note that we need to clarify the following things about the experiment above:
1. how many qubits did they use
2. how were the performance / efficiency of the simulation hardware and source code?

Moreover, Rao et al.(2020) analyzed the performance between MLQAE and IQAE on IBM Q device, Vigo, and on IBM Q simulator respectively. And they performed with $2,\ 3$ and $10$ qubits for each experiment as following:
>
>
>

We can observe some key points from the above results:
* For each experiments, they conducted the number of shots from $2^4$ to $2^{10}$, and each run corredponding to a shot number is repeated $30$ times. This procedure leads to a statistically converged value for estimated $a$, the amplitude of the "good" state.
* For fig. 6 and 7, we can observe the gap for relative error for $a$ between IBM Vigo and IBM Q Simulator.
* For fig. 8, since the number of qubits is $10$, they can only use IBM Q Simulator to conduct the experiment.

Note:
1. The most important parameter for MLQAE is $m$, which means the number of application of $Q$. When $m$ changes from $3$ to $4$, the system will generate noise by the reason of decoherence. Hence, Rao et al. have set $m=3$ while performing the experiments.
2. Again, we consider the most vital parameter fot IQAE is $\epsilon$, which refers to the error tolerance and determines $k$, the power of $Q$. They have taken $\epsilon=0.01$.

>**Conclusion**
>* IQAE seems to perform better as shown in Figs. 6, 7, and 9. The main reason for the inaccuracy of MLQAE on the quantum device is that the fourth circuit, $(Q^4\mathcal{A}|0\rangle_{n+1})$, has
long circuit depth, which is proportional to the number of $Q$, while the IQAE runs have at most $Q^2$ even though they have more iterations than four.
>* When we tested IQAE with error tolerance $0.005$ on devices, the runs readily needed $Q^4$ and even $Q^7$. Thus, the results was so noisy that they did not have meaning. However, in the absence of quantum noise, MLQAE with $m=3$ and IQAE with $0.01$ error tolerance show similar performance as shown in Fig. 10.
> 
> 

# Future work and open problems
1. Discuss the first step of QFT applied in QAE algorithm
2. Deduce the performance of IQAE algorithm for both CP and CH bounding method
3. Reproduce the comparison results of Grinko et al. (2019) by using Qiskit
4. Connect the result to the option pricing problem, or to solve a relatively smaller problem

## References
1. Brassard, G., Hoyer, P., Mosca, M., and Tapp, A., *"Quantum amplitude amplication and estimation,"* Contemporary Mathematics 305, 53-74 (2002).
2. Suzuki, Y., Uno, S., Raymond, R., Tanaka, T., Onodera, T., and Yamamoto, N., *"Amplitude estimation without phase estimation,"* Quantum Information Processing 19(2), 75 (2020).
3. Grinko, D., Gacon, J., Zoufal, C., and Woerner, S., *"Iterative quantum amplitude estimation,"* arXiv preprint arXiv:1912.05559 (2019).
4. Pooja, R., Kwangmin, Y., Hyunkyung, L., Dasol, J., and Deokkyu, C., "Quantum amplitude estimation algorithms on IBM quantum devices", Aug(2020)
5. Aaronson, S. & Rall, P. *"Quantum Approximate Counting, Simplified."*, (2019)
6. IBM Quantum Learning, *"Grover's Algorithm"*