Quantum Scientific Computation

Content

• Basic Quantum Algorithm

  • Phase Kickback

  • Grover’s Algorithm

  • Hadamard Test

  • Quantum Fourier Transform

  • Quantum Phase Estimation

  • Trotter Based Hamiltonian Simulation

• Matrix Computation over Quantum Computer

  • Block Encoding

  • Matrix Functions of Hermitian matrices

  • Quantum Singular Value Transformation

NOTE: More quantum algorithms can be found at Zoo of quantum algorithms.

Basic quantum algorithms

Some of them may be counterintuitive and sometimes it looks like we are trying to find applications of a particular quantum phenomena (formula).

phase kickback

Let $f: {0, 1}^n \rightarrow {0, 1}$ be a boolean function, which can be queried via the following unitary

i.e., the value of $f(x)$ is returned as a phase factor.

One of the quantum circuit to implement phase kickback is to use oracle as follows with an extra uncomputation ancilla qubit:

The concrete circuit depends on the function $f(x)$. Denote $\vert -\rangle = \frac{1}{\sqrt{2}}(\vert 0\rangle - \vert 1\rangle$, then we have

ignore the ancilla qubit $\vert -\rangle$ we have the phase kickback unitary $U_f$.

Remark. Combine phase kickback with Hadamard gates we could get the Deutsch’s algorithm designed for the judgement of balanced boolean functions.

Grover’s algorithm

Problem setting (Unstructured search): Given a boolean function $f:{0, 1}^n \rightarrow {0, 1}$, there exists a unique marked state $x_0$ such that $f(x_0)=1$ and we would like to find $x_0$.

Main idea: Use an oracle to implement a phase kick, which turn out to be a reflection for some particular states, then combine two reflection to form a rotation. With the rotation operator to move the initial state to desired $\vert x_0\rangle$.

The origin of quadratic speed up comes from the fact that classical probabilistic algorithms work with probability density, quantum algorithms work with wavefunction amplitudes, of which the square gives the probability density.

• Step 1: Oracle and phase kickback:

Given oracle

Let $\vert y\rangle=\vert -\rangle$, we have phase kickback with an ancilla qubit:

• Step 2: initial state preparation:

Our target is to rotate this state to close $\vert x_0\rangle$.

• Step 3: rotation operator:

Apply oracle $U_f$ to initial state we have

Thus $U_f$ are actually the reflection operator

In the same way one have reflection with $\psi_0$,

Thus we have Grover rotation operator $G=R_{\psi_0}R_{x_0}$ such that

where we define $\vert \psi_0\rangle = -\sin(\theta/2)\vert x_0\rangle + \cos(\theta/2)\vert \psi_0^\bot\rangle$.

Finally, for $\sin((2k+1)\theta/2) \approx 1$ we need $k\approx \frac{\pi}{2\theta} - \frac{1}{2} \approx \frac{\sqrt{N}\pi}{4}$.

Remark:

Grover’s algorithm is not restricted to the problem of unstructured search. One immediate application is amplitude amplification, i.e., for

we can amplify the amplitude of $\vert \psi_{\text{good}}\rangle$ to improve the success probability.

Hadamard Test

Hadamard test is a useful tool for computing the expectation value of a unitary operator $U$ (in general not Hermitian) with respect to a state $\psi$, i.e., $\langle \psi \vert U \vert \psi \rangle$.

The real part and imaginary part of the expectation value need to be measured separately.

From the Hadamard test for real part expectation value, we have

The probability of measuring the qubit 0 to be state $\vert 0\rangle$ is

Circuit model of Haramard test for $\text{Re} \langle \psi \vert U\vert \psi\rangle$:

Circuit model of Haramard test for $\text{Im} \langle \psi \vert U\vert \psi\rangle$:

Quantum Fourier transform

Definition:

For any $j$ in the computational basis, the (discrete) forward Fourier transform is defined as:

The (discrete) inverse Fourier transform is

The circuit implementation of QFT relies on the following formula:

where we use the $n$-bit fixed point representation, i.e., for $a = \frac{k}{2^n} \in [0, 1)$ with integer $k$,

Quantum phase estimation

Problem:

Let $U$ be a unitary, and $\vert \psi\rangle$ is an eigenvector, i.e.,

the goal is to find $\varphi$ up to certain precision.

Algorithm:

Suppose $\varphi = \frac{k}{2^d}$ for some $k\in [2^d]$, $d$ is the number of ancilla qubits, then we can use a pair of QFT and the power of $U$ to extract number $k$ by measuring ancilla qubits.

where $\mathcal{U}=\sum_{j\in[2^d]} \vert j\rangle \langle j\vert \otimes U^j $. $\vert \psi_0\rangle$ is the eigenvector corresponding to $\varphi$, which can be obtained in advanced or just used for immediate step during theoretic analysis.

Useful form:

Given a Hermitian matrix $A$ with eigenpair $A v = \lambda v$ satisfying $0 < \lambda < 1$ with an exact $d$-bit representation and $|v|=1$, then we have unitary operator $U_\text{QPE}= \mathrm{e}^{\mathrm{i}2\pi A}$ such that

where $\vert \lambda\rangle$ is the $d$-bit representation of $\lambda$.

Application: HHL algorithm for linear systems

Step 1: Convert $Ax=b$ into quantum linear system problem (QLSP)

We assume $b\in \mathbb{C}^N$ is a normalized vector and $\vert b\rangle=U_b \vert 0^n\rangle$.

Note that $ \vert 1\rangle\langle 0\vert =\begin{bmatrix} 0 & 0 \ 1 & 0 \end{bmatrix} $, if $A$ is not Hermitian, we can enlarge the system to get a Hermitian matrix (dilation method), i.e.,

and solve $\widetilde{A}\vert 0, x\rangle=\vert 1, b\rangle$.

The QLSP is to find a quantum state $\vert \widetilde{x}\rangle$ so that

Step 2: Solve QLSP with eigendecomposition of Hermitian matrix $A$

Suppose $A$ has the eigendecomposition $A\vert v_j\rangle=\lambda_i \vert v_j\rangle$. Then $\vert b\rangle = \sum_j \beta_j \vert v_j\rangle$, the unnormalized solution satisfies

Step 3: Use quantum circuit to express solution of QLSP

Assume $0 < \lambda_0 \leq \lambda_1 \leq \cdots \leq \lambda_{N_1} < 1 $ and all eigenvalues have an exact $d$-bit representation.

Consider unitary matrix $U =\mathrm{e}^{\mathrm{i}2\pi A}$ and quantum phase estimation,

With controlled rotation unitary

where each $\widetilde{\lambda}_j$ approximates $\lambda_j$, $C>0$ is a constant decided by user being a lower bound to $\lambda_0$ such that $0 < C/\lambda_j < 1$ for all j.

Then we have

Hence if the measurement of signal qubit is 1, we obtain

which is stored as a normalized vector

To retrieve the norm of the solution, we compute the success probability of measuring signal qubit

Step 0: complexity analysis

The success probability $p(1)$ is determined by

• the choice of the normalization constant $C$
• the norm of the true solution $|x|=|A^{-1}\vert b\rangle|$

For theoretical analysis, assume $C=\lambda_0$, if assume $\lambda_{N-1}=1$ further, one can connect with the condition number of matrix $A$, i.e.,

Futhermore, $|A^{-1}\vert b\rangle| \geq \frac{1}{|A|}|\vert b\rangle|=1$. Therefore $p(1)=\Omega(\kappa^{-2}).$

After connecting with the condition number, now we can perform the comparison with classical iterative linear system solvers and find applications that quantum computing can do better.

Note: implementation of CR gate

Given rotation angle in its $d$-bit fixed point representation, $\theta=.\theta_{d-1}\cdots\theta_{0}$,

then

However, to realize the desired controlled rotation for QLSP, one need to convert the coefficient into rotation angle, i.e.,

where $C > 0$ such that $0 < C/\lambda_j < 1$ is user-specific constant. The arcsin function can be implemented using classical arithmetics circuit.

Trotter based Hamiltonian simulation

Problem: given an initial state $\vert \psi_0\rangle$ and a time-independent Hamiltonian $H$, evaluate the quantum state at time $t$ via

Specifically, consider $H=H_1+H_2$, then one has Trotter splitting:

Let $\Delta t = t/n$, then we have first order Trotter method:

And further second order Trotter method

Matrix Computation over Quantum Computer

All the operations over quantum computer are unitary. Thus we have to encode the general matrix into unitary matrix and represent (approximate) their operations with unitary operations.

Block Encoding

Question:

How to read-in and read-out the information of a matrix $A\in \mathbb{C}^{N \times N}$, which is generally not unitary in a quantum computer?

One possible method is via $\mathrm{e}^{\mathrm{i}t A}$ (if $A$ is not Hermitian, we can use the dilation method).

Some more general methods will be discussed below.

Query model for matrix entries (read-out)

Assume $A$ is an $n$-qubit, square matrix and $\vert A \vert_{\max}:=\max_{ij} \vert A_{ij}\vert < 1$. The desired query oracle takes the following general form:

It is indeed a problem that how to implement this oracle. We will not get into the detail here.

Heuristic Illustration:

Find a $(n+1)$-qubit unitary matrix $U\in \mathbb{C}^{2N\times 2N}$ such that

In order to obtain $A\vert b\rangle$, consider

and

we only need to measure the qubit 0 and keep the state if it returns 0. The missing info of norm $|A\vert b\rangle|$ can be recovered via the success probability $p(0)$.

Definition (Block encoding)

Given an $n$-qubit matrix $A$, if we can find $\alpha, \epsilon\in\mathbb{R}_+$, and an $(m+n)$-qubit unitary matrix $U_A$ such that

then $U_A$ is called an $(\alpha, m, \epsilon)$-block-encoding of $A$, which consists the set denoted as $\text{BE}_{\alpha, m}(A, \epsilon)$.

Note: if

then $A=(\langle 0^m\vert \otimes I_n) U_A (\vert 0^m\rangle \otimes I_n)$ by some expansion of tensor products.

Example:

For any $n$-qubit matrix $A$ with $|A|_2 \leq 1$, the SVD of $A$ is denoted as $A=W\Sigma V^{\dagger}$, then we may construct

which is a $(1,1)$-block-encoding of $A$.

Note: Instead of first identifying $A$ and then finding its block encoding $U_A$, sometimes we reverse this though process: we first identifying unitary $U_A$ that is easy to implement on quantum computer, and then ask which matrix can be block encoded by $U_A$.

Matrix Functions of Hermitian matrices

Many tasks in scientific computation can be expressed in terms of matrix functions.

• Hamiltonian simulation: $f(A)=e^{\mathrm{i}At}$.
• Gibbs state preparation: $f(A)=e^{-\beta A}$.
• Solving linear systems of equation: $f(A)=A^{-1}$.
• Eigenstate filtering: $f(A)=\mathbf{1}_{(-\infty, 0)} (A-\mu I)$.

Matrix function of Hermitian matrices is defined by the function of its eigenvalues, given function $f:\mathbb{R}\rightarrow \mathbb{C}$ and matrix $A=V\Lambda V^{\dagger}$, we define

Target: Compute any matrix function over quantum computer

Linear combination of unitaries

Let $T=\sum_{i=0}^{K-1} \alpha_i U_i$ be the linear combination of $n$-by-$n$ unitary matrices $U_i$ (general matrix can be extended to a unitary matrices through block encoding). WLOG assume $K=2^a$, $\alpha_i > 0$.

Select oracle

Prepare oracle

LCU Lemma:

Define

Then for any $\vert \psi\rangle$,

where $\vert \widetilde{\bot}\rangle$ is an un-normalized state satisfying $(\vert 0^a\rangle \langle 0^a\vert \otimes I_n) \vert \widetilde{\bot}\rangle = 0.$

Example: Linear combination of two matrices

Qubitization of Hermitian matrices

Essential Fact: An arbitrarily unitary block encoding matrix $U_A$ can be partially block diagonalized into some subblocks of size $2\times 2$, with the transformed basis depending on the eigenvectors of $A$.

This fact gives a method to design the circuit for the Chebyshev polynomials of matrix $A$ since for rotation matrix $O$,

one has

where

are Chebyshev polynomials.

Similar, suppose $U_A$ is Hermitian, denote the eigen decomposition of $A$ as

Define

then basis $\mathcal{B}_i = {\vert 0^m\rangle\vert v_i\rangle, \vert \bot_i\rangle }$ forms a invariant subspace of $U_A$.

In order to turn $U_A$ as a rotation, introduce

hence,

Thus,

Note that $T_k(A)= \sum_i T_k(\lambda_i) \vert v_i\rangle \langle v_i\vert $, we then have

which gives the block encoding of the Chebyshev polynomial of matrix $A$.

If $U_A$ is not Hermitian, we have

Define

then we have

Let $\widetilde{O}=U_A^{\dagger} Z_{\Pi}U_A Z_{\Pi}$, then

Therefore,

Note that $U_A$ and $U_A Z_{\Pi} (U_A^{\dagger} Z_{\Pi}U_A Z_{\Pi})^k $ shares the same invariant subspaces, we know that

Quantum eigenvalue transformation

The construction of matrix function over quantum computer above is not practical due to the difficulty of oracle and complexity.

Below we generalize the block encoding for Chebyshev polynomial and find methods to make it more practical.

The block encoding of Chebyshev polynomial can be written as

Further, consider

For any $\widetilde{\Phi} := (\widetilde{\phi}_0, \cdots, \widetilde{\phi}_d ) \in \mathbb{R}^{d+1}$,

where $P\in \mathbb{C}[x]$ satisfying certain conditions.

Conversely, we have the conclusion known as quantum signal processing (QSP).

There exists a set of phase factors $\Phi:=(\phi_0, \cdots, \phi_d) \in \mathbb{R}^{d+1}$ such that

if and only if $P, Q\in\mathbb{C}[x]$ satisfy

• $\text{deg}(P) \leq d$, $\text{deg}(Q) \leq d - 1$,

• $P$ has parity $d \mod 2$ and $Q$ has parity $d-1 \mod 2$ and

• $\vert P(x)\vert ^2 + (1-x^2)\vert Q(x)\vert ^2 = 1$, $\forall x\in [-1, 1]$.

Note: if we are only interested in QSP for real polynomials $f(x)$, let $f(x)=\text{Re}(P(x))$, then the third condition above can be significantly relaxed.

Optimization method for finding phase factors

Given a target real polynomial $f(x)$ with degree $d$ and parity $d \mod 2$, $f(x)$ is entirely determined by its values on $\tilde{d}=\lceil\frac{d+1}{2} \rceil$ points $x_k$. The QSP problem can be equivalently solved via the following optimization problem:

where $g(x, \Phi):=\text{Re}(U(x, \Phi)).$

Quantum Singular Value Transformation

Generalized matrix function

Given SVD of matrix $A\in\mathbb{C}^{N\times N}$ as $A=W\Sigma V^{\dagger}$ or equivalently

Let $f:\mathbb{R}\rightarrow \mathbb{C}$ be a scalar function, the generalized matrix function is defined as

Besides,

Qubitization of general matrices

and $\widetilde{O}=U_A^{\dagger} Z_{\Pi}U_A Z_{\Pi}$, then

Therefore,

Remark: note that all singular values are non-negative, when performing the polynomial approximation, one can always extend the function $f$ to match the desired parity. For example,

QSVT with real polynomial

Under the same conditions as QET, we can find sequence of phase factors $\Phi\in\mathbb{R}^{d+1}$, so that for even $d$,

gives a $(1, m+1)$-block-encoding of $P^{\triangleright}(A)$ for some even polynomial $P\in\mathbb{R}[x]$; for add $d$,

gives a $(1, m+1)$-block-encoding of $P^{\diamond}(A)$ for some odd polynomial $P\in \mathbb{R}[x]$.

In general, the QSVT circuit takes the form

Remark:

This formulation gives the freedom to insert unitary matrices $P, Q$ and lead to the equivalent QSVT transformation

Assume $\widetilde{U}_A=QU_A P^{\dagger}$, then $U_A$ is the matrix representation of $\widetilde{U}_A$ with respect to the bases given by $P, Q$, respectively. This can be useful for certain application, e.g., fixed-point amplitude amplification.