We introduce the four main postulates of quantum mechanics related to this course. For more details, we refer readers to [1, Section 2.2]. All postulates concern closed quantum systems (i.e., systems isolated from environments) only.

2.1.1 State space postulate

2.1.2 Quantum operator postulate

2.1.3 Quantum measurement postulate

2.1.4 Tensor product postulate

So far all quantum states encountered can be described by a single \( \ket{\psi}\in\mathcal{H}\) , called the pure state. More generally, if a quantum system is in one of a number of states \( \ket{\psi_i}\) with respective probabilities \( p_i\) , then \( \{p_i, \ket{\psi_i}\}\) is an ensemble of pure states. The density operator of the quantum system is

For a pure state \( \ket{\psi}\) , we have

is a rank-\( 1\) matrix.

Consider a quantum observable in (19) associated with the projectors \( \{P_m\}\) . For a pure state, it can be verified that the probability result of returning \( \lambda_m\) , and the expectation value of the measurement are respectively,

The expression (39) also holds for general density operators \( \rho\) .

An operator \( \rho\) is the density operator associated to some ensemble \( \{p_i, \ket{\psi_i}\}\) if and only if (1) \( \operatorname{Tr} \rho=1\) (2) \( \rho \succeq 0\) , i.e., \( \rho\) is a positive semidefinite matrix (also called a positive operator). All postulates in (2.1) can be stated in terms of density operators (see [1, Section 2.4.2]). Note that a pure state satisfies \( \rho^2=\rho\) . In general we have \( \rho^2\preceq \rho\) . If \( \rho^2\prec \rho\) , then \( \rho\) is called (the density operator of) a mixed state. Furthermore, an ensemble of admissible density operators is also a density operator.

A quantum operator \( U\) that transforms \( \ket{\psi}\) to \( U\ket{\psi}\) also transforms the density operator according to

However, not all quantum operations on density operators need to be unitary! See [1, Section 8.2] for more general discussions on quantum operations.

Most of the discussions in this course will be restricted to pure states, and unitary quantum operations. Even in this restricted setting, the density operator formalism can still be convenient, particularly for describing a subsystem of a composite quantum system. Consider a quantum system of \( (n+m)\) -qubits, partitioned into a subsystem \( A\) with \( n\) qubits (the state space is \( \mathcal{H}_A=\mathbb{C}^{2^n}\) ) and a subsystem \( B\) with \( m\) qubits (the state space is \( \mathcal{H}_B=\mathbb{C}^{2^m}\) ) respectively. The quantum state is a pure state \( \ket{\psi}\in \mathbb{C}^{2^{n+m}}\) with density operator \( \rho_{AB}\) . Let \( \ket{a_1},\ket{a_2}\) be two state vectors in \( \mathcal{H}_A\) , and \( \ket{b_1},\ket{b_2}\) be two state vectors in \( \mathcal{H}_B\) . Then the partial trace over system \( B\) is defined as

Since we can expand the density operator \( \rho_{AB}\) in terms of the basis of \( \mathcal{H}_A,\mathcal{H}_B\) , the definition of (41) can be extended to define the reduced density operator for the subsystem \( A\)

The reduced density operator for the subsystem \( B\) can be similarly defined. The reduced density operators \( \rho_A,\rho_B\) are generally mixed states.

If a quantum observable is defined only on the subsystem \( A\) , i.e., \( M=M_A\otimes I\) where \( M_A\) has the decomposition (19), then the success probability of returning \( \lambda_m\) , and the expectation value are respectively

Nearly all quantum algorithms operate on multi-qubit quantum systems. When quantum operators operate on two or more qubits, writing down quantum states in terms of its components as in (29) quickly becomes cumbersome. The quantum circuit language offers a graphical and compact manner for writing down the procedure of applying a sequence of quantum operators to a quantum state. For more details see [1, Section 4.2, 4.3].

In the quantum circuit language, time flows from the left to right, i.e., the input quantum state appears on the left, and the quantum operator appears on the right, and each “wire” represents a qubit i.e.,

Here are a few examples:

which is a graphical way of writing

The relation between these states can be expressed in terms of the following diagram

Also verify that

which is a graphical way of writing

Note that the input state can be general, and in particular does not need to be a product state. For example, if the input is a Bell state (32), we just apply the quantum operator to \( \ket{00}\) and \( \ket{11}\) , respectively and multiply the results by \( 1/\sqrt{2}\) and add together. To distinguish with other symbols, these single qubit gates may be either written as \( X,Y,Z,H\) or (using the roman font) \( \mathrm{X,Y,Z,H}\) .

The quantum circuit for the CNOT gate is

Here the “dot” means that the quantum gate connected to the dot only becomes active if the state of the qubit \( 0\) (called the control qubit) is \( a=1\) . This justifies the name of the CNOT gate (controlled NOT).

Similarly,

is the controlled \( U\) gate for some unitary \( U\) . Here \( U^a=I\) if \( a=0\) . The CNOT gate can be obtained by setting \( U=X\) .

Another commonly used two-qubit gate is the SWAP gate, which swaps the state in the \( 0\) -th and the \( 1\) -st qubits.

Quantum operators applied to multiple qubits can be written in a similar manner:

For a multi-qubit quantum circuit, unless stated otherwise, the first qubit will be referred to as the qubit 0, and the second qubit as the qubit 1, etc.

When the context is clear, we may also use a more compact notation for the multi-qubit quantum operators:

\[ \Leftrightarrow \]

\[ \Leftrightarrow \]

One useful multiple qubit gate is the Toffoli gate (or controlled-controlled-NOT, CCNOT gate).

We may also want to apply a \( n\) -qubit unitary \( U\) only when certain conditions are met

where the empty circle means that the gate being controlled only becomes active when the value of the control qubit is \( 0\) . This can be used to write down the quantum “if” statements, i.e., when the qubits \( 0,1\) are at the \( \ket{1}\) state and the qubit \( 2\) is at the \( \ket{0}\) state, then apply \( U\) to \( \ket{x}\) .

A set of qubits is often called a register (or quantum variable). For example, in the picture above, the main quantum state of interest (an \( n\) qubit quantum state \( \ket{x}\) ) is called the system register. The first \( 3\) qubits can be called the control register.

One of the most striking early results of quantum computation is the no-cloning theorem (by Wootters and Zurek, as well as Dieks in 1982), which forbids generic quantum copy operations (see also [1, Section 12.1]). The no-deleting theorem is a consequence of linearity of quantum mechanics.

Assume there is a unitary operator \( U\) that acts as the copy operations, i.e.,

for any black-box state \( x\) , and a chosen target state \( \ket{s}\) (e.g. \( \ket{0^n}\) ). Then take two states \( \ket{x_1},\ket{x_2}\) , we have

Taking the inner product of the two equations, we have

which implies \( \braket{x_1|x_2}=0\) or \( 1\) . When \( \braket{x_1|x_2}=1\) , \( \ket{x_1},\ket{x_2}\) refer to the same physical state. Therefore a cloning operator \( U\) can at most copy states which are orthogonal to each other, and a general quantum copy operation is impossible.

Given the ubiquity of the copy operation in scientific computing like \( y=x\) , the no-cloning theorem has profound implications. For instance, all classical iterative algorithms for solving linear systems require storing some intermediate variables. This operation is generally not possible, or at least cannot be efficiently performed.

There are two notable exceptions to the range of applications of the no-cloning theorem. The first is that we know how a quantum state is prepared, i.e., \( \ket{x}=U_x\ket{s}\) for a known unitary \( U_x\) and some \( \ket{s}\) . Then we can of course copy this specific vector \( \ket{x}\) via

The second is the copying of classical information. This is an application of the CNOT gate.

i.e.,

The same principle applies to copying classical information from multiple qubits. (1) gives an example of copying the classical information stored in 3 bits.

In general, a multi-qubit CNOT operation can be used to perform the classical copying operation in the computational basis. Note that in the circuit model, this can be implemented with a depth 1 circuit, since they all act on different qubits.

The quantum no-cloning theorem implies that there does not exist a unitary \( U\) that performs the deleting operation, which resets a black-box state \( \ket{x}\) to \( \ket{0^n}\) . This is because such a deleting unitary can be viewed as a copying operation

Then take \( \ket{x_1},\ket{x_2}\) that are orthogonal to each other, apply the deleting gate, and compute the inner products, we obtain

which is a contradiction.

A more common way to express the no-deleting theorem is in terms of the time reversed dual of the no-cloning theorem: in general, given two copies of some arbitrary quantum state, it is impossible to delete one of the copies. More specifically, there is no unitary \( U\) performing the following operation using known states \( \ket{s},\ket{s'}\) ,

for an arbitrary unknown state \( \ket{x}\) ((7)).

The quantum measurement applied to any qubit, by default, measures the outcome in the computational basis. For example,

outputs \( 0\) or \( 1\) each w.p. \( 1/2\) . We may also measure some of the qubits in a multi-qubit system.

\[ \equiv \]

There are two important principles related to quantum measurements: the principle of deferred measurement, and the principle of implicit measurement. At a first glance, both principles may seem to be counterintuitive.

The principle of deferred measurement states that measurement operations can always be moved from an intermediate stage of a quantum circuit to the end of the circuit. This is because even if a measurement is performed as an intermediate step in a quantum circuit, and the result of the measurement is used to conditionally control subsequent quantum gates, such classical controls can always be replaced by quantum controls, and the result of the quantum measurement is postponed to later.

Example 9 (Deferring quantum measurements)

Consider the circuit

Here the double line denotes the classical control operation. The outcome is that qubit \( 0\) has probability \( 1/2\) of outputting \( 0\) , and the qubit \( 1\) is at state \( \ket{0}\) . Qubit \( 0\) also has probability \( 1/2\) of outputting \( 1\) , and the qubit \( 1\) is at state \( \ket{1}\) .

However, we may replace the classical control operation after the measurement by a quantum controlled \( X\) (which is CNOT), and measure the qubit \( 0\) afterwards:

It can be verified that the result is the same. Note that CNOT acts as the classical copying operation. So qubit \( 1\) really stores the classical information (i.e., in the computational basis) of qubit \( 0\) .

Example 10 (Deferred measurement requires extra qubits)

The procedure of deferring quantum measurements using CNOTs is general, and important. Consider the following circuit:

The probability of obtaining \( 0,1\) is \( 1/2\) , respectively. However, if we simply “defer” the measurement to the end by removing the intermediate measurement, we obtain

The result of the measurement is deterministically \( 0\) ! The correct way of deferring the intermediate quantum measurement is to introduce another qubit

Measuring the qubit \( 0\) , we obtain \( 0\) or \( 1\) w.p. \( 1/2\) , respectively. Hence when deferring quantum measurements, it is necessary to store the intermediate information in extra (ancilla) qubits, even if such information is not used afterwards.

The principle of implicit measurements states that at the end of a quantum circuit, any unmeasured qubit may be assumed to be measured. More specifically, assume the quantum system consists of two subsystems \( A\) and \( B\) . If qubits \( A\) are to be measured at the end of the circuits, the results of the measurements does not depend on whether the qubits \( B\) are measured or not. Recall from (44) that a measurement on the subsystem \( A\) only depends on the reduced density matrix \( \rho_A\) . So we only need to show that \( \rho_A\) does not depend on the measurement in \( B\) . To see why this is the case, let \( \{P_i\}\) be the projectors onto the computational basis of \( B\) . Before the measurement, the density operator is \( \rho\) . If we measure the subsystem \( B\) , the resulting density operator is transformed into

Then it can be verified that

This proves the principle of implicit measurements.

By definition, the output of all quantum algorithms must be obtained through measurements, and hence the measurement outcome is probabilistic in general. If the goal is to compute the expectation value of a quantum observable \( M_A\) acting on a subsystem \( A\) , then its variance is

The number of samples \( \mathcal{N}\) needed to estimate \( \operatorname{Tr}\) to additive precision \( \epsilon\) satisfies

which only depends on \( \rho_A\) .

If a quantum algorithm denoted by a unitary \( U\) can be decomposed into a series of simpler unitaries as \( U=U_{K}… U_1\) , and if we can implement each \( U_i\) to precision \( \epsilon\) , then what is the global error? We now introduce a simple technique connecting the local error with the global error. In the context of quantum computation, this is often referred to as the “hybrid argument”.

In other words, if we can implement each local unitary to precision \( \epsilon\) , the global error grows at most linearly with respect to the number of gates and is bounded by \( K\epsilon\) . The telescoping series (65), as well as the hybrid argument can also be seen as a discrete analogue of the variation of constants method (also called Duhamel’s principle).

As a special case, consider \( B(t)=E(t)\widetilde{U}(t)\) , then (68) becomes

and

This is a direct analogue of the hybrid argument in the continuous setting.

In classical computation, there are many universal gate sets, in the sense that any classical gate can be represented as a combination of gates from the set. For example, the NAND gate (“Not AND”) alone forms a universal gate set [1, Section 3.1.2]. The NOR gate (“Not OR”) is also a universal gate set.

In the quantum setting, any unitary operator on \( n\) qubits can be implemented using \( 1\) - and \( 2\) -qubit gates[1, Section 4.5]. It is desirable to come up with a set of discrete universal gates, but this means that we need to give up the notion that the unitary \( U\) can be exactly represented. Instead, a set of quantum gates \( \mathcal{S}\) is universal if given any unitary operator \( U\) and desired precision \( \epsilon\) , we can find \( U_1,…,U_m\in \mathcal{S}\) such that

Here \( \left\lVert A\right\rVert =\sup_{\braket{\psi|\psi}=1} \left\lVert A\ket{\psi}\right\rVert \) is the operator norm (also called the spectral norm) of \( A\) , and \( \left\lVert \ket{\psi}\right\rVert =\sqrt{\braket{\psi|\psi}}\) is the vector \( 2\) -norm). There are many possible choices of universal gate sets, e.g. \( \{H,T, \operatorname{CNOT} \}\) . Another universal gate set is \( \{H, \operatorname{Toffoli} \}\) , which only involves real numbers.

Are some universal gate sets better than others? The Solovay-Kitaev theorem states that all choices of universal gate sets are asymptotically equivalent (see e.g.  [8, Chapter 2]):

Another natural question is about the computational power of quantum computers. Perhaps surprisingly, it is very difficult to prove that quantum computer is more powerful than classical computer. But is quantum computer at least as powerful as classical computers? The answer is yes! More specifically, any classical circuit can also be asymptotically efficiently implemented using a quantum circuit.

The proof rests on that the classical universal gate can be efficiently simulated using quantum circuits. Note that this is not a straightforward process: NAND, and other classical gates (such as AND, OR etc.) are not reversible gates! Hence the first step is to perform classical computation with reversible gates. More specifically, any irreversible classical gate \( x\mapsto f(x)\) can be made into a reversible classical gate

In particular, we have \( (x,0)\mapsto (x,f(x))\) computed in a reversible way. The key idea is to store all intermediate steps of the computation (see [1, Section 3.2.5] for more details).

On the quantum computer, storing all intermediate computational steps indefinitely creates two problems: (1) tremendous waste of quantum resources (2) the intermediate results stored in some extra qubits are still entangled to the quantum state of interest. So if the environments interfere with intermediate results, the quantum state of interest is also affected.

Fortunately, both problems can be solved by a step called “uncomputation”. In order to implement a Boolean function \(f :\{0,1\}^{n} \rightarrow\{0,1\}^{m}\), we assume there is an oracle

where \(\ket{0^m}\) comes from a \(m\)-qubit output register. The oracle is often further implemented with the help of a working register (a.k.a. “garbage” register) such that

From the no-deleting theorem, there is no generic unitary operator that can set a black-box state to \(\ket{0^w}\). In order to set the working register back to \(\ket{0^w}\) while keeping the input and output state, we introduce yet another \(m\)-qubit ancilla register initialized at \(\ket{0^m}\). Then we can use an \( n\) -qubit CNOT controlled on the output register and obtain

It is important to remember that in the operation above, the multi-qubit CNOT gate only performs the classical copying operation in the computational basis, and does not violate the no-cloning theorem.

Recall that \( U_f^{-1}=U_f^{†}\) , we have

Finally we apply an \( n\) -qubit SWAP operator on the ancilla and output registers to obtain

After this procedure, both the ancilla and the working register are set to the initial state. They are no longer entangled to the input or output register, and can be reused for other purposes. This procedure is called uncomputation. The circuit is shown in (2).

Using the technique of uncomputation, if the map \( x\mapsto f(x)\) can be efficiently implemented on a classical computer, then we can implement this map efficiently on a quantum computer as well. To do this, we first turn it into a reversible map (73). All reversible single-bit and two-bit classical gates can be implemented using single-qubit and two-qubit quantum gates. So the reversible map can be made into a unitary operator

on a quantum computer. This proves that a quantum computer is at least as powerful as classical computers.

The unitary transformation \( U_f\) in (81) can be applied to any superposition of states in the computational basis, e.g.

This does not necessarily mean that we can efficient implement the map \( \ket{x}\mapsto \ket{f(x)}\) . However, if \( f\) is a bijection, and we have access to the inverse of the reversible circuit for computing \( f^{-1}\) , then we may use the technique of uncomputation to implement such a map ((9)).

Let \([N]=\set{0,1,…,N-1}\). Any integer \(k\in[N]\) where \(N=2^n\) can be expressed as an \(n\)-bit string as \(k=k_{n-1}… k_0\) with \(k_i\in\{0,1\}\). This is called the binary representation of the integer \( k\) . It should be interpreted as

The number \( k\) divided by \( 2^m\) (\( 0\le m\le n\) ) can be written as (note that the decimal is shifted to be after \( k_m\) ):

The most common case is \( m=n\) , where

Sometimes we may also write \( a=0.k_{1}… k_{n}\) , so that \( k_i\) is the \( i\) -th decimal of \( a\) in the binary representation. For a given floating number \( 0\le a<1\) written as

the number \( (0.k_1… k_n)\) is called the \( n\) -bit fixed point representation of \( a\) . Therefore to represent \( a\) to additive precision \( \epsilon\) , we will need \( n=\lceil\log_2 \epsilon\rceil\) qubits. If the sign of \( a\) is also important, we may reserve one extra bit to indicate its sign.

Together with the reversible computational model, we can perform classical arithmetic operations, such as \( (x,y)\mapsto x+y\) , \( (x,y)\mapsto xy\) , \( x\mapsto x^{\alpha}\) , \( x\mapsto \cos(x)\) etc. using reversible quantum circuits. The number of ancilla qubits, and the number of elementary gates needed for implementing such quantum circuits is \( \mathcal{O}(\operatorname{poly}(n))\) (see [4, Chapter 6] for more details).

It is worth commenting that while quantum computer is theoretically as powerful as classical computers, there is a very significant overhead in implementing reversible classical circuits on quantum devices, both in terms of the number of ancilla qubits and the circuit depth.

All previous discussions assume that quantum operations can be perfectly performed. Due to the immense technical difficulty for realizing quantum computers, both quantum gates and quantum measurements may involve (significant) errors, particularly on near-term quantum devices. However, the threshold theorem states that if the noise in individual quantum gates is below a certain constant threshold (around \( 10^{-4}\) or above), it is possible to efficiently perform an arbitrarily large quantum computation with any desired precision (see [1, Section 10.6]). This procedure requires quantum error correction protocols.

This course will not discuss any details on quantum error corrections. We always assume fault-tolerant protocols have been implemented, and all errors come from either approximation errors at the mathematical level, or Monte Carlo errors in the readout process due to the probabilistic nature of the measurement process.

Let \( n\) be the number of qubits needed to represent the input. A quantum algorithm is efficient if the number of gates in the quantum circuit is \( \mathcal{O}(\operatorname{poly}(n))\) . Due to the probabilistic nature of the measurement outcome, we are typically satisfied if a quantum algorithm can produce the correct answer with sufficiently high probability \( p\) . For a decision problem that asks for a binary answer \( 0\) or \( 1\) , we require \( p>2/3\) (or at least \( p > 1/2+1/\operatorname{poly}(n)\) ). For other problems that we have an efficient procedure to check the correctness of the answer, we require \( p=\Omega(1)\) . Repeating this process many times and apply the Chernoff bound [1, Box 3.4], we can make the probability of outputting an incorrect answer vanishingly small.

In quantum algorithms, the computational cost is often measured in terms of the query complexity. Assume that we have access to black-box unitary operator \( U_f\) (e.g. the one used in the reversible computation), which is often called a quantum oracle. Our goal is to perform a given task using as few queries as possible to \( U_f\) .

The concept of query complexity hides the implementation details of \( U_f\) , and in some cases we can prove lower bounds on the number of queries to solve a certain problem, e.g. in the case of Grover’s search (proving a lower bound of the number of gates among all quantum algorithms can be much harder). Furthermore, once we have access to the number of elementary gates needed to implement \( U_f\) , we obtain immediately the gate complexity of the total algorithm. However, some queries can be (provably) difficult to implement, and then there can be a large gap between the query complexity and gate complexity. In order to obtain a meaningful query complexity analysis, one should also make sure that other components of the quantum algorithm will not end up dominating the total gate complexity, when all factors are taken into account.

Another important measure of the complexity is the circuit depth, i.e., the maximum number of gates along any path from an input to an output. Since quantum gates can be performed in parallel, the circuit depth is approximately equivalent to the concept of “wall-clock time” in classical computation, i.e., the real time needed for a quantum computer to carry out a certain task. Since quantum states can only be preserved for a short period of time (called the coherence time), the circuit depth also provides an approximate measure of whether the quantum algorithm exceeds the coherence limit of a given quantum computer. In many scenarios, the maximum coherence time is the most severe limiting factor. When possible, it is often desirable to reduce the circuit depth, even if it means that the quantum circuit needs to be carried out many more times.

Let us summarize the basic components of a typical quantum algorithm: the set of qubits can be separated into system registers (storing quantum states of interest) and ancilla registers (auxiliary registers needed to implement the unitary operation acting on system registers). Starting from an initial state, apply a series of one-/two-qubit gates, and perform measurements. Uncomputation should be performed whenever possible. Within the ancilla registers, if a register can be “freed” after the uncomputation, it is called a working register. Since working registers can be reused for other purposes, the cost of working registers is often not (explicitly) factored into the asymptotic cost analysis in the literature.

We use \( \|\cdot\|\) to denote vector or matrix 2-norm: when \( v\) is a vector we denote by \( \|v\|\) its 2-norm, and when \( A\) is matrix we denote by \( \|A\|\) its operator norm. Other matrix and vector norms will be introduced when needed. Unless otherwise specified, a vector \( v\in\mathbb{C}^N\) is an unnormalized vector, and a normalized vector (stored as a quantum state) is denoted by \( \ket{v}=v/\left\lVert v\right\rVert \) . A vector \( v\) can be expressed in terms of \( j\) -th component as \( v=(v_j)\) or \( (v)_j=v_j\) . We use a \( 0\) -based indexing, i.e., \( j=0,…,N-1\) or \( j\in [N]\) . When \( 1\) -based indexing is used, we will explicitly write \( j=1,…,N\) . We use the following asymptotic notations besides the usual \( \mathcal{O}\) (or “big-O”) notation: we write \( f=\Omega(g)\) if \( g=\mathcal{O}(f)\) ; \( f=\Theta(g)\) if \( f=\mathcal{O}(g)\) and \( g=\mathcal{O}(f)\) ; \( f=\widetilde{\mathcal{O}}(g)\) if \( f=\mathcal{O}(g\operatorname{polylog}(g))\) .

Assume we have two boxes, each of them may contain either an apple or an orange. We would like to answer: whether the two boxes contain the same type of fruit (but do not need to answer whether it is apple or orange).

This seems to be a weird question. If the content of a box can only be checked by opening it, then to answer the question we would need to open two boxes and check what is inside. It is impossible to answer whether the fruit types are the same without knowing the types! Nevertheless, this is precisely the question to be addressed by Deutsch’s algorithm.

Mathematically, consider a boolean function \( f:\{0,1\}\to\{0,1\}\) . The question is whether \( f(0)=f(1)\) or \( f(0)\ne f(1)\) ? A quantum fruit-checker assumes the access to \( f\) via the following quantum oracle:

A classical fruit-checker can only query \( U_f\) in the computational basis as

After these two queries, we can measure qubit \( 1\) with a deterministic outcome, and answer whether \( f(0)=f(1)\) . However, a quantum fruit-checker can apply \( U_f\) to a linear combination of states in the computational basis.

Let us first check that \( U_f\) is unitary:

which gives \( U_f^{†}U_f=I\) .

The idea behind Deutsch’s algorithm is to convert the oracle (90) into a phase kickback in (87). Take \( \ket{y}=\ket{-}=\frac{1}{\sqrt{2}}(\ket{0}-\ket{1})\) . Then

Note that \( \ket{y}=H X\ket{0}\) , (93) can also be interpreted as

The application of \( XH\) can be viewed as the step of uncomputation. Neglecting the qubit 1 which does not change before and after the application, we can focus on the first qubit only, which effectively defines a unitary

Hence the information of \( f(x)\) is stored as a phase factor (\( 0\) or \( \pi\) ). Recall that using a Hadamard gate \( H\ket{0}=\ket{+}, H\ket{1}=\ket{-}\) , the quantum circuit of Deutsch’s algorithm is

or in the commonly seen form in (3).

The answer is embedded in the measurement outcome of qubit \( 0\) . To verify this:

So if \( f(0)=f(1)\) , the final state is \( \pm\ket{0,-}\) . Measuring qubit \( 0\) returns \( 0\) deterministically (the globally phase factor is irrelevant). Similarly if \( f(0)\ne f(1)\) , the final state is \( \pm\ket{1,-}\) . Measuring qubit \( 0\) returns \( 1\) deterministically. In summary, only one query to \( U_f\) is sufficient to answer whether the two boxes contain the same type of fruit. The procedure is equally counterintuitive. Note that a classical fruit checker as implemented in (91) naturally receives the information by measuring qubit 1. On the other hand, Deutsch’s algorithm only uses qubit 1 as a signal qubit, and all the information is retrieved by measuring qubit 0, which, at least from the classical perspective of (91), seems to contain no information at all!

Now we have seen that in terms of the query complexity, a quantum fruit-checker is clearly more efficient. However, it is a fair question how to implement the oracle \( U_f\) , especially in a way that somehow does not already reveal the values of \( f(0),f(1)\) . We give the implementation of some cases of \( U_f\) in (14). In general, proving the query complexity alone may not be convincing enough that quantum computers are better than classical computers, and the gate complexity matters. This will not be the last time we hide away such “implementation details” of quantum oracles.

Example 14 (Qiskit example of Deutsch’s algorithm)

For \( f(0)=f(1)=1\) (constant case), we can use \( U_f=I\otimes X\) . For \( f(0)=0,f(1)=1\) (balanced case), we can use \( U_f= \operatorname{CNOT} \) .

Assume we have \( N=2^{n}\) boxes, and we are given the promise that only one of the boxes contains an orange, and each of the remaining boxes contains an apple. The goal is to find the box that contains the orange.

Mathematically, given a boolean function \( f:\{0,1\}^n\to\{0,1\}\) and the promise that there exists a unique marked state \( x_0\) that \( f(x_0)=1\) , we would like to find \( x_0\) . This is called an unstructured search problem. Classically, there is no simpler methods than opening \( (N-1)\) boxes in the worst case to determine \( x_0\) .

The quantum algorithm below, called Grover’s algorithm, relies on access to an oracle

and can find \( x_0\) using \( \mathcal{O}(\sqrt{N})\) queries. A classical computer again can only query \( U_f\) in the computational basis, and Grover’s algorithm achieves a quadratic speedup in terms of the query complexity.

The origin of the quadratic speedup can be summarized as follows: while classical probabilistic algorithms work with probability densities, quantum algorithms work with wavefunction amplitudes, of which the square gives the probability densities. More specifically, we start from a uniform superposition of all states as the initial state

This state can be prepared using Hadamard gates as

We would like to amplify the desired amplitude corresponding to \( \ket{x_0}\) from \( 1/\sqrt{N}\) to \( \sqrt{p}=\Omega(1)\) . We demonstrate below that this requires \( \mathcal{O}(\sqrt{N})\) queries to \( U_f\) . After this procedure, by measuring the final state in the computational basis, we obtain some output state \( \ket{x}\) . We can check whether \( x=x_0\) by applying another query of \( U_f\) according to \( U_f \ket{x,0}=\ket{x,f(x)}\) . The probability of obtaining \( f(x)=1\) is \( p\) . If \( f(x)\ne 1\) , we repeat the process. Then after \( \mathcal{O}(1/p)\) times of repetition, we can obtain \( x_0\) with high probability.

The first step of Grover’s algorithm is to turn the oracle (97) into a phase kickback. For this we take \( \ket{y}=\ket{-}\) , and for any \( x\in\{0,1\}^n\) ,

Any quantum state \( \ket{\psi}\) can be decomposed as

where \( \ket{\psi^{\perp}}\) is the component of \( \ket{\psi}\) orthogonal to \( \ket{x_0}\) , i.e., \( \braket{\psi^{\perp}|x_0}=0\) . We have

Here the minus sign is gained through the phase kickback. Discarding the \( \ket{-}\) which is unchanged by applying \( U_f\) , we obtain an \( n\) -qubit unitary

Therefore \( R_{x_0}\) is a reflection operator across the hyperplane orthogonal to \( \ket{x_0}\) , i.e., the Householder reflector

Let us write

with \( \theta=2\arcsin \frac{1}{\sqrt{N}}\approx \frac{2}{\sqrt{N}}\) , and \( \ket{\psi_0^{\perp}}=\frac{1}{\sqrt{N-1}}\sum_{x\ne x_0} \ket{x}\) is a normalized state orthogonal to \( \ket{x_0}\) . Then

So \( \operatorname{span} \{\ket{x_0},\ket{\psi_0^{\perp}}\}\) is an invariant subspace of \( R_{x_0}\) .

The next key step is to consider another Householder reflector with respect to \( \ket{\psi_0}\) . For later convenience we add a global phase factor \( -1\) (which is irrelevant to the physical outcome):

Direct computation shows

So define \( G=R_{\psi_0}R_{x_0}\) as the product of the two reflection operators (called the Grover operator), then it amplifies the angle from \( \theta/2\) to \( 3\theta/2\) . The geometric picture is in fact even clearer in (4) and the conclusion can be observed without explicit computation.

Applying the Grover operator \( k\) times, we obtain

So 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}\) . This proves that Grover’s algorithm can solve the unstructured search problem with \( \mathcal{O}(\sqrt{N})\) queries to \( U_f\) .

Another derivation of the Grover method is to focus on the operator, instead of the initial vector \( \ket{\psi_0}\) at each step of the calculation. In the orthonormal basis \( \mathcal{B}=\{\ket{x_0},\ket{\psi_0^{\perp}}\}\) , the matrix representation of the reflector \( R_{x_0}=I-2\ket{x_0}\bra{x_0}\) is

The matrix representation for the Grover diffusion operator \( R_{\psi_0}=2\ket{\psi_0}\bra{\psi_0}-I\) is

Here \( \sin(\theta/2)=a=1/\sqrt{N}\) . Therefore for the matrix representation of the Grover iterate \( G=R_{\psi_0}R_{x_0}\) is

i.e., \( G\) is a rotation matrix restricted to the two-dimensional space \( \mathcal{H}= \operatorname{span} \mathcal{B}\) .

The initial vector satisfies

so Grover’s search can be applied as before, via \( G^k\) for \( k\approx \frac{\pi\sqrt{N}}{4}\) times.

To draw the quantum circuit of Grover’s algorithm, we need an implementation of \( R_{\psi_0}\) . Note that

This can be implemented via the following circuit using one ancilla qubit: }

Here the controlled-NOT gate is an \( n\) -qubit controlled-\( X\) gate, and is only active if the system qubits are in the \( 0^n\) state. Discarding the signal qubit, we obtain and implementation of \( R_{\psi_0}\) . Since the signal qubit \( \ket{-}\) only changes up to a sign, it can be reused for both \( R_{\psi_0}\) and \( R_{x_0}\) .

The reflector \( R_{\psi_0}\) can also be implemented without using the ancilla qubit as (use a 3-qubit system as an example)

Grover’s algorithm is not restricted to the problem of unstructured search. One immediate application is called amplitude amplification (AA), which is used ubiquitously as a subroutine to achieve quadratic speedups.

Let \(\ket{\psi_0}\) be prepared by an oracle \( U_{\psi_0}\) , i.e., \( U_{\psi_0}\ket{0^n}=\ket{\psi_0}\) . We have the knowledge that

and \( p_0\ll 1\) . Here \( \ket{\psi_{\mathrm{bad}}}\) is an orthogonal state to the desired state \( \ket{\psi_{\mathrm{good}}}\) . We cannot get access to \( \ket{\psi_{\mathrm{good}}}\) directly, but would like to obtain a state that has a large overlap with \( \ket{\psi_{\mathrm{good}}}\) , i.e., amplify the amplitude of \( \ket{\psi_{\mathrm{good}}}\) .

In the problem of unstructured search, \( \ket{\psi_{\mathrm{good}}}=\ket{x_0}\) , and \( p_0=1/N\) . Although we do not have access to the answer \( \ket{x_0}\) , we assume access to the reflection oracle \( R_{x_0}\) . Here we also assume access to the reflection oracle

From \( U_{\psi_0}\) , we can construct the reflection with respect to the initial state

via the \( n\) -qubit controlled-\( X\) gate. So following exactly the same procedure as the unstructured search problem, we can construct the Grover iterate

Applying \( G^k\) to \( \ket{\psi_0}\) for some \( k=\mathcal{O}(1/\sqrt{p_0})\) , we obtain a state that has \( \Omega(1)\) overlap with \( \ket{\psi_{\mathrm{good}}}\) .

Recall that the unstructured search problem tries to find a marked state \( x_0\in[N]\) , using a reflection oracle

Grover’s algorithm can find \( x_0\) with constant probability (e.g. at least \( 1/2\) ) by making \( \mathcal{O}(\sqrt{N})\) times querying \( R_{x_0}\) . It turns out that this is asymptotically optimal, i.e., no quantum algorithm can perform this task using fewer than \( \Omega(\sqrt{N})\) access to \( R_{x_0}\) .

Any quantum search algorithm that starts from a universal initial state \( \ket{\psi_0}\) and queries \( R_{x_0}\) for \( k\) steps can be written in the following form:

for some unitaries \( U_1,…,U_k\) . For simplicity we assume no ancilla qubits are used, and the proof can be generalized to the case in the presence of ancilla qubits. The superscript \( x_0\) indicates that the state depends on the marked state \( x_0\) . Specifically, by “solving” the search problem, it means that there exists a \( k\) so that for each marked state \( x_0\) ,

In other words, measuring \( \ket{\psi^{x_0}_k}\) in the computational basis, the probability of obtaining \( \ket{x_0}\) is at least \( 1/2\) .

To prove the lower bound, we compare the action of \( \mathcal{U}^{x_0}_k\) with another “fake algorithm” \( \mathcal{U}_k\) , defined as

In particular, \( \ket{\psi_k}\) does not involve any information of the solution \( x_0\) and therefore cannot possibly solve the search problem.

For a set of vectors \( \{f^{x_0}\}_{x_0\in [N]}\) , and each \( f^{x_0}\in \mathbb{C}^N\) , we will extensively use the following discrete \( \ell_2\) -norm:

In particular, we have the following triangle inequality

The proof contains two steps. First, we show that the true solution and the fake solution differs significantly, in the sense that

Second, we prove that (define \( D_0=0\) ):

Therefore to satisfy (133), we must have \( k=\Omega(\sqrt{N})\) .

In the first step, since multiplying a phase factor \( e^{\mathrm{i} \theta}\) to \( \ket{\psi^{x_0}_k}\) does not have any physical consequences, we may choose a particular phase \( \theta\) so that (129) becomes

Therefore

This means that

On the other hand, for the “fake algorithm”, using the Cauchy-Schwarz inequality,

This violates the bound in (137), and the fake algorithm cannot solve the search problem for arbitrarily large \( k\) .

So from (137,138), and the triangle inequality, we have

This proves (133). In other words, the true solution and the fake solution must be well separated in \( \ell_2\) -norm.

In the second step, we prove (134) inductively. Clearly (134) is true for \( k=0\) . Assume this is true, then

The last inequality uses the triangle inequality of the discrete \( \ell_2\) -norm. Note that

and

we have

which finishes the induction.

Finally, combining the lower bound (133,134), we find that the necessary condition to solve the unstructured search problem is \( 4k^2=\Omega(N)\) , or \( k=\Omega(\sqrt{N})\) .

We first introduce the Hadamard test, which is a useful tool for computing the expectation value of an unitary operator with respect to a state, i.e., \( \braket{\psi|U|\psi}\) . Since \( U\) is generally not Hermitian, this does not correspond to the measurement of a physical observable. Instead the real and imaginary part of the expectation value need to be measured separately.

The (real) Hadamard test is the quantum circuit in (6) for estimating the real part of \( \braket{\psi|U|\psi}\) .

To verify this, we find that the circuit transforms \( \ket{0}\ket{\psi}\) as

The probability of measuring the qubit \( 0\) to be in state \( \ket{0}\) is

This is well defined since \( -1\le \operatorname{Re} \braket{\psi|U|\psi}\le 1\) .

To obtain the imaginary part, we can use the circuit in (7) called the (imaginary) Hadamard test, where

is called the phase gate.

Similar calculation shows the circuit transforms \( \ket{0}\ket{\psi}\) to the state

Therefore the probability of measuring the qubit \( 0\) to be in state \( \ket{0}\) is

Combining the results from the two circuits, we obtain the estimate to \( \braket{\psi|U|\psi}\) .

Example 18 (Overlap estimate using swap test)

An application of the Hadamard test is called the swap test, which is used to estimate the overlap of two quantum states \( \left\lvert \braket{\varphi|\psi}\right\rvert \) . The quantum circuit for the swap test is

Note that this is exactly the Hadamard test with \( U\) being the \( n\) -qubit swap gate. Direct calculation shows that the probability of measuring the qubit \( 0\) to be in state \( \ket{0}\) is

\[ \begin{equation} p(0)=\frac{1}{2}(1+\operatorname{Re} \braket{\varphi,\psi|\psi,\varphi})=\frac12(1+\left\lvert \braket{\varphi|\psi}\right\rvert ^2). \end{equation} \]

(150)

Example 21 (Qiskit example for phase estimation using the Hadamard test)

Here is a qiskit example of the simple phase estimation for

\[ R|1\rangle = \begin{pmatrix} 1 & 0\\ 0 & e^{i2\pi\varphi}\\ \end{pmatrix} \begin{pmatrix} 0\\ 1\\ \end{pmatrix} = e^{\I2\pi\varphi}|1\rangle, \]

with \( \varphi=0.5+\frac{1}{2^d}\equiv 0.{10… 01}\) (d bits in total). In this example, \( d=4\) and the exact value is \( \varphi=0.5625\) .

In (20), the number of measurements needed to estimate \( \varphi\) to precision \( \epsilon\) is \( \mathcal{O}(1/\epsilon^2)\) . A quadratic improvement in precision (i.e., \( \mathcal{O}(1/\epsilon)\) ) can be achieved by means of the quantum phase estimation. One such procedure is called Kitaev’s method [2, Section 13.5].

In the fixed point representation in (2.8), for simplicity we assume that the eigenvalue can be exactly represented using \( d\) bits, i.e.,

In the simplest scenario, we assume \( d=1\) , and \( \varphi=.\varphi_0,\varphi_0\in\{0,1\}\) . Then \( e^{\text{i}2 \pi \varphi}=e^{\text{i} \pi \varphi_0}\) .

Performing the real Hadamard test in (20), we have \( p(1)=0\) if \( \varphi_0=0\) , and \( p(1)=1\) if \( \varphi_0=1\) . In either case, the result is deterministic, and one measurement is sufficient to determine the value of \( \varphi_0\) .

Next, consider \( \varphi=.\underbrace{0… 0\varphi_0}_{d bits }\) . To determine the value of \( \varphi_0\) , we need to reach precision \( \epsilon<2^{-d}\) . The method in (20) requires \( \mathcal{O}(1/\epsilon^2)=\mathcal{O}(2^{2d})\) repeated measurements, or number of queries to \( U\) . The observation from Kitaev’s method is that if we can have access to \( U^j\) for a suitable power \( j\) , then the number of queries to \( U\) can be reduced. More specifically, if we can query \( U^{2^{d-1}}\) , then the circuit in (9) with \( j=d-1\) gives

The result is again deterministic. Therefore the total number of queries to \( U\) becomes \( \mathcal{O}(2^{-d})=\mathcal{O}(\epsilon^{-1})\) .

This is the basic idea behind Kitaev’s method: use a more complex quantum circuit (and in particular, with a larger circuit depth) to reduce the total number of queries. As a general strategy, instead of estimating \( \varphi\) from a single number, we assume access to \( U^{2^j}\) , and estimate \( \varphi\) bit-by-bit. In particular, changing \( U\to U^{2^j}\) in the Hadamard test allows us to estimate

One immediate difficulty of the bit-by-bit estimation is that we need to tell \( 0.0111…\) apart from \( 0.1000…\) , and the two numbers can be arbitrarily close to each other (though the two numbers can also differ at some number of digits), and some careful work is needed. We will first describe the algorithm, and then analyze its performance. The algorithm works for any \( \varphi\) , and then the goal is to estimate its \( d\) bits. For simplicity of the analysis, we assume \( \varphi\) is exactly represented by \( d\) bits. We will use extensively the distance

which is the distance on the unit circle.

First, by applying the circuit in (9) (and the corresponding circuit to determine the sign) with \( j=0,1,…,d-3\) , for each \( j\) we can estimate \( p(0)\) , so that the error in \( 2^j\varphi\) is less than \( 1/16\) for all \( j\) (this can happen with a sufficiently high success probability. For simplicity, let us assume that this happens with certainty). The measured result is denoted by \( \alpha_j\) . This means that any perturbation must be due to the \( 5\) th digit in the binary representation. For example, if \( 2^j\varphi=0.11100\) , then \( \alpha_j=0.11011\) is an acceptable result with an error \( 0.00001=1/32\) , but \( \alpha_j=0.11110\) is not acceptable since the error is \( 0.0001=1/16\) . We then round \( \alpha_j\) \( \pmod 1\) by its closest \( 3\) -bit estimate denoted by \( \beta_j\) , i.e., \( \beta_j\) is taken from the set \( \set{0.000,0.001,0.010,0.011,0.100,0.101,0.110,0.111}\) . Consider the example of \( 2^j\varphi=0.11110\) , if \( \alpha_j=0.11101\) , then \( \beta_j=0.111\) . But if \( \alpha_j=0.11111\) , then \( \beta_j=0.000\) . Another example is \( 2^j\varphi=0.11101\) , if \( \alpha_j=0.11110\) , then both \( \beta_j=0.111\) (rounded down) and \( \beta_j=0.000\) (rounded up) are acceptable. We can pick one of them at random. We will show later that the uncertainty in \( \alpha_j,\beta_j\) is not detrimental to the success of the algorithm.

Second, we perform some post-processing. Start from \( j=d-3\) , we can estimate \( .\varphi_2\varphi_1\varphi_0\) to accuracy \( 1/16\) , which recovers these three bits exactly. The values of these three bits will be taken from \( \beta_{d-3}\) directly. Then we proceed with the iteration: for \( j=d-4,…,0\) , we assign

Here \( \left\lvert \cdot\right\rvert _{\bmod 1}\) is the periodic distance on \( [0,1)\) and its value is always \( \le 1/2\) . Since the two possibilities are separated by \( 1/2\) , for each \( j\) , there will be at most one case that is satisfied. We will also show that in all circumstances, there is always one case that is satisfied, regardless of the ambiguity of the choice of \( \beta_j\) above.

After running the algorithm above, we recover \( \varphi=.\varphi_{d-1}… \varphi_0\) exactly. The total cost of Kitaev’s method measured by the number of queries to \( U\) is \( \mathcal{O}\left(\sum_{j=0}^{d-3} 2^{j}\right)=\mathcal{O}(\epsilon^{-1})\) .

If \( \varphi\) is exactly represented by \( d\) bits, we will obtain an estimate

Let us now inductively show why Kitaev’s algorithm works. Again assume \( \varphi\) is exactly represented by \( d\) bits. For \( j=d-3\) , we know that \( \varphi_2\varphi_1\varphi_0\) can be recovered exactly. Then assume \( \varphi_{d-j-2}…\varphi_0\) have all been exactly computed, at step \( j\) we would like to determine the value of \( \varphi_{d-j-1}\) . From

we know

Then

The wrong choice of \( \varphi_{d-j-1}\) denoted by \( \widetilde{\varphi}_{d-j-1}\) then satisfies

This proves the validity of (161), and hence that of Kitaev’s algorithm.

Fourier transform is used ubiquitously in scientific computing, and fast Fourier transform (FFT) is the backbone for many fast algorithms in classical computation. Similarly, the quantum Fourier transform is also an important component in many quantum algorithms, such as phase estimation, Shor’s algorithm, and inspires other fast algorithms such as fast Fermionic Fourier transform (FFFT) etc.

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

In particular

The (discrete) inverse Fourier transform is

Using the binary representation of integers

we have

Therefore the exponential can be written as

The most important step of QFT is the following direct calculation, which requires some patience with the manipulation of indices:

(173) involves a series of controlled rotations of the form

Hence before discussing the quantum circuit for QFT, let us first work out the circuit for implementing this controlled rotation. We use the relation

Example 24 (Implementation of controlled rotation)

Consider the implementation of

\[ \begin{equation} \ket{0}\ket{j}\to \frac{1}{\sqrt{2}}\left(\ket{0}+e^{\mathrm{i} 2\pi (.j_{n-1}\cdots j_0)}\ket{1}\right)\ket{j}, \end{equation} \]

(176)

Let

\[ \begin{equation} R_z(\theta)=\begin{pmatrix} 1 & 0 \\ 0 & e^{\mathrm{i} \theta} \end{pmatrix}, \end{equation} \]

(177)

and \( R_j=R_z(\pi/2^{j-1})\) . In particular, \( R_1=Z\) . The quantum circuit is

The implementation of QFT follows the same principle, but does not require the signal qubit to store the phase information. Let us see a few examples.

When \( n=1\) , we need to implement

This is the Hadamard gate:

When \( n=2\) , we need to implement

This can be implemented using the following circuit:

Comparing the result with that in (180), we find that the ordering of the qubits is reversed. To recover the desired result in QFT, we can apply a SWAP gate to the outcome, i.e.,

In order to implement the inverse Fourier transform, we only need to apply the Hermitian conjugate as

Similarly one can construct the circuit for \( U_{\mathrm{FT}}\) and its inverse for \( n=3\) .

In general, the QFT circuit is given by (10). Compare the circuit in (10) with (173), we find again that the ordering is reversed in the output. To restore the correct order as defined in QFT, we can use \( \mathcal{O}(n/2)\) swaps operations. The total gate complexity of QFT is \( \mathcal{O}(n^2)\) .

In Kitaev’s method, we use \( 1\) ancilla qubit but \( d\) different circuits of various circuit depths to perform phase estimation. In this section, we introduce the (standard) quantum phase estimation (QPE), which uses one signal quantum circuit based on QFT, but requires \( d\) -ancilla qubits to store the phase information in the quantum computer. From now, we assume \( \varphi=.\varphi_{d-1}… \varphi_0\) is exact. From the availability of \( U^j\) we can define a controlled unitary operation

When \( d=1\) , \( \mathcal{U}\) is simply the controlled \( U\) operation. For a general \( d\) , it seems that we need to implement all \( 2^d\) different \( U^j\) . However, this is not necessary. Using the binary representation of integers \( j=(j_{d-1}… j_0.)=\sum_{i=0}^{d-1} j_i 2^{i}\) , we have \( U^{j}=U^{\sum_{i=0}^{d-1} j_i 2^i}=\prod_{i=0}^{d-1}U^{j_i2^i}\) . Therefore similar to the operations in QFT,

Here the primed product \( \sideset{}{’}\prod \) is a slightly awkward notation, which means the tensor product for the first register, and the regular matrix product for the second register. It is in fact much clearer to observe the structure in the quantum circuit in (11).

Now let the initial state in the ancilla qubits be \( \ket{0^n}\) . Use QFT and \( \mathcal{U}\) , we transform the initial states according to

Since we have \( \varphi=\frac{k}{2^d}\) for some \( k\in[2^d]\) , measuring the ancilla qubits, and we will obtain the state \( \ket{k}\ket{\psi_0}\) with certainty, and we obtain the phase information. Therefore the quantum circuit for QFT based QPE is given by (12). Here we have used (169). We should note that \( U_{\mathrm{FT}}^{†}\) includes the swapping operations. (Exercise: 1. what if the swap operation is not implemented? 2. Is it possible to modify the circuit and remove the need of implementing the swap operations?)

In order to apply QPE (standard or Kitaev), we have assumed that

1. \( \ket{\psi_0}\) is an eigenstate.
2. \( \varphi_0\) has an \( d\) -bit binary representation.

In general practical calculations, neither condition can be exactly satisfied, and we need to analyze the effect on the error of the QPE. Recall the discussion in (2.9), we assume the only sources of errors are at the mathematical level (instead of noisy implementation of quantum devices). In this context, the error can be due to an inexact eigenstate \( \ket{\psi}\) , or Monte Carlo errors in the readout process due to the probabilistic nature of the measurement process. In this section, we relax these constraints and study what happens when the conditions are not exactly met. We assume \( U\) has the eigendecomposition.

Without loss of generality we assume \( 0\le \varphi_0\le \varphi_1… \le\varphi_{N-1}<1\) . We are interested in using QPE to find the value of \( \varphi_0\) .

We first relax the condition (1), i.e., assume all \( \varphi_i\) ’s have an exact \( d\) -bit binary representation, but the quantum state is given by a linear combination

Here the overlap \( p_0=\left\lvert \braket{\phi|\psi_0}\right\rvert ^2=\left\lvert c_0\right\rvert ^2<1\) .

Applying the QPE circuit in (12) to \( \ket{0^t}\ket{\phi}\) with \( t=d\) , and measure the ancilla qubits, we obtain the binary representation of \( \varphi_0\) with probability \( p_0\) . Furthermore, the system register returns the eigenstate \( \ket{\psi_0}\) with probability \( p_0\) . Of course, in order to recognize that \( \varphi_0\) is the desired phase, we need some a priori knowledge of the location of \( \varphi_0\) , e.g. \( \varphi_0\in(a,b)\) and \( \varphi_i>b\) for all \( i\ne 0\) . It would be desirable to relax both conditions (1) and (2). However, the analysis can be rather involved and additional assumptions are needed. We will discuss some of the implications in the context of estimating the ground state energy in (5.1).

For now, to simplify the analysis, we focus on the case that only the condition (2) is violated, i.e., \( \varphi_0\) cannot be exactly represented by a \( d\) -bit number, and we need to apply the QPE circuit to an initial state \( \ket{0^t}\ket{\phi}\) with \( t>d\) . The exact relation between the \( t\) and the desired accuracy \( d\) will be determined later. Similar to (183), we obtain the state

Here

Therefore if \( \varphi_0\) has an exact \( d\) -bit representation, i.e., \( \varphi_0=\widetilde{\varphi}_{k'_0}\) for some \( k'_0\) , then \( \gamma_{0,k'}=\delta_{k',k'_0}\) . We recover the previous result that one run of the QPE circuit gives the value \( \varphi_0\) deterministically.

Now assume that \( \varphi_0 \ne \widetilde{\varphi}_{k'}\) for any \( k'\) . Note that \( e^{\mathrm{i} 2\pi x}\) is a periodic function with period \( 1\) , we can only determine the value of \( x \bmod 1\) . Therefore we use the periodic distance (160). In terms of the phase, we would like to find \( k_0'\) such that

Here \( \epsilon=2^{-d}=2^{t-d}/T\) is the precision parameter. In particular, for any \( k'\) we have

Using the relation that for any \( \theta\in[-\pi,\pi]\) ,

we obtain

Let \( k'_0\) be the measurement outcome, which can be viewed as a random variable. The probability of obtaining some \( \widetilde{\varphi}_{k'_0}\) that is at least \( \epsilon\) distance away from \( \varphi_0\) is

Set \( t-d=\lceil\log_2 \delta^{-1}\rceil\) , then \( T\epsilon=2^{t-d}\ge \delta^{-1}\) . Hence for any \( 0<\delta<1\) , the failure probability

In other words, in order to obtain the phase \( \varphi_0\) to accuracy \( \epsilon=2^{-d}\) with a success probability at least \( 1-\delta\) , we need \( d+\lceil\log_2 \delta^{-1}\rceil\) ancilla qubits to store the value of the phase. On top of that, the simulation time needs to be \( T=(\epsilon\delta)^{-1}\) .

As an application, we use QPE to solve the problem of estimating the ground state energy of a Hamiltonian. Let \( H\) be a Hermitian matrix with eigendecomposition

Below are two examples of Hamiltonians commonly encountered in quantum many-body physics.

Without loss of generality we assume \( 0<\lambda_0\le\lambda_1\le … <\lambda_{N-1}<\frac12\) . Note that for the purpose of estimating the ground state energy, we do not necessarily require a positive energy gap. For simplicity of the presentation, we still assume that the ground state is non-degenerate, i.e., \( \lambda_0<\lambda_1\) . We are also provided an approximate eigenstate

of which the overlap with the ground state is \( p_0=\left\lvert \braket{\phi|\psi_0}\right\rvert ^2\) . Our goal is to estimate \( \lambda_0\) to precision \( \epsilon=2^{-d}\) . We assume \( \epsilon<\lambda_0\) . This appears in many problems in quantum many-body physics, quantum chemistry, optimization etc.

In order to use QPE (based on QFT), we assume access to the unitary evolution operator \( U=e^{\mathrm{i} 2\pi H}\) . This is called a Hamiltonian simulation problem, which will be discussed in detail in later chapters. For now we assume \( U\) can be implemented exactly. Then

This becomes a phase estimation problem, where the input vector is not an exact eigenstate. Following the discussion in (4.5), if all eigenvalues \( \lambda_j\) can be exactly represented by \( d\) -bit numbers, we obtain both the ground state and the ground state energy with probability \( p_0\) . Therefore repeating the process for \( \mathcal{O}(p_0^{-1})\) times we obtain the ground state energy.

Now we relax both conditions (1) and (2) in (4.5), and apply the QPE circuit in (12) to an initial state \( \ket{0^t}\ket{\phi}\) for some \( t>d\) . Similar to (183), we have

Here

Therefore the definition in (187) is a special case.

Our algorithm is simple: we would like to run QPE \( M\) times, and denote the output of the \( \ell\) -th run by \( \widetilde{\varphi}^{(\ell)}_{k'}\) . Then we take the minimum of the measured output \( \min_{\ell}\widetilde{\varphi}^{(\ell)}_{k'}\) as the estimate to the ground state energy. The hope is to obtain the ground state energy to accuracy \( \epsilon\) with success probability at least \( 1-\delta\) for any \( \delta>0\) . Let us now analyze this algorithm.

If \( \varphi_0\) has an exact \( d\) -bit representation, i.e., \( \lambda_0=\widetilde{\varphi}_{k'_0}\) for some integer \( k'_0\) , then \( \gamma_{k'_0,k'}=\delta_{k'_0,k'}\) . It may seem that this implies with probability \( p_0\) , we obtain the exact estimate of \( \varphi_0\) , and correspondingly the eigenstate \( \ket{\psi_0}\) is stored in the system register. This is much better than the previous assumption that all eigenvalues \( \lambda_j\) need to be represented by a \( d\) -bit number.

Unfortunately, this analysis is not correct. In fact, for any \( \lambda_k\) that does not have an exact \( t\) -bit representation (note that \( t>d\) ), we may have \( \gamma_{k,k''}\ne 0\) and \( \widetilde{\varphi}_{k''}<\lambda_0\) , i.e., we obtain an energy estimate that is lower than the ground state energy! Therefore the probability of ending up in the state \( \ket{k''}\ket{\psi_k}\) is \( \left\lvert c_k\gamma_{k,k''}\right\rvert ^2\) , i.e., it is still possible to obtain a wrong ground state energy. This is called the aliasing effect.

We demonstrate below that if \( T\) is large enough, we can control the probability of underestimating the ground state energy. Since \( \lambda_0\) is the ground state energy, and all eigenvalues are in \( (0,1/2)\) , when \( \widetilde{\varphi}_{k'}\le \lambda_0-\epsilon\) , we have

Then the probability of under estimating the ground state energy by \( \epsilon\) is

Let \( T\epsilon=\delta'^{-1}\) with \( \delta'<1/2\) , we have

Therefore after \( M\) repetitions, we have

In order to obtain \( P\left(\min_{\ell}\widetilde{\varphi}^{(\ell)}_{k'}\le \lambda_0-\epsilon\right)<\frac{\delta}{2}\) , we need to set \( \delta'=M^{-1} \delta\) .

On the other hand, we also would like to have bound \( P\left(\min_{\ell}\widetilde{\varphi}^{(\ell)}_{k'}\ge \lambda_0+\epsilon\right)\) . To this end, we first note that when \( \widetilde{\varphi}_{k'}\ge \lambda_0+\epsilon\) , we have \( \left\lvert \widetilde{\varphi}_{k'}- \lambda_0\right\rvert _1\ge\epsilon\) . Moreover,

Here we have used the normalization condition that

Therefore

This means that

We can then take \( M=\lceil \frac{2}{p_0} \log \frac{2}{\delta} \rceil\) so that

To summarize, according to the relation \( \delta'=M^{-1} \delta\) , in order to estimate the ground state energy to precision \( \epsilon=2^{-d}\) with success probability \( 1-\delta\) , we need

ancilla qubits in QPE. The circuit depth is

Taking the number of repetitions \( M\) into account, the total cost of the method is

Through the analysis above, we see that although the analysis of QPE is very clean when 1) all eigenvalues (properly scaled to be represented as phase factors) are exactly given by \( d\) -bit numbers 2) the input vector is an eigenstate, the analysis can become rather complicated and tedious when such conditions are relaxed. Such difficulty does not merely show up at the theoretical level, but can seriously impact the robust performance of QPE in practical applications. To simplify the discussion of the applications below, we will be much more cavalier about the usage of QPE and assume all eigenvalues are exactly represented by \( d\) -bit numbers whenever necessary. But we should keep such caveats in mind. Furthermore, when we move beyond QPE, the issue of having exact \( d\) -bit numbers will become much less severe in techniques based on quantum signal processing, i.e., quantum eigenvalue transformation (QET) and quantum singular value transformation (QSVT).

Let \(\ket{\psi_0}\) be prepared by an oracle \( U_{\psi_0}\) , i.e., \( U_{\psi_0}\ket{0^n}=\ket{\psi_0}\) . We have the knowledge that

and \( p_0\ll 1\) . Here \( \ket{\psi_{\mathrm{bad}}}\) is an orthogonal state to the desired state \( \ket{\psi_{\mathrm{good}}}\) , and \( \sqrt{p_0}=\sin\frac{\theta}{2}\) . The problem of amplitude estimation is to estimate \( p_0\) to \( \epsilon\) -precision. If \( p_0\) is away from \( 0\) or \( 1\) and is to be estimated directly from the Monte Carlo method, the number of samples needed is \( \mathcal{N}=\mathcal{O}(\epsilon^{-2})\) .

Let \( G=R_{\psi_0}R_{\mathrm{good}}\) be the Grover operator as in (3.3). Then in the basis \( \mathcal{B}=\{\ket{\psi_{\mathrm{good}}},\ket{\psi_{\mathrm{bad}}}\}\) , the subspace \( \mathcal{H}= \operatorname{span} B \) is an invariant subspace of \( G\) . Recall the computation of (112), the matrix representation is

Its two eigenstates are

with eigenvalues \( e^{\pm \mathrm{i} \theta}\) , respectively.

Therefore the problem of estimating \( \theta\) can be solved with phase estimation with an imperfect initial state. Note that

Consider a QPE circuit with \( t\) ancilla qubits and querying \( G\) for \( T=2^t\) times. Then each execution with the system register will be in \( \ket{\psi_+}\) or \( \ket{\psi_{-}}\) states each with probability \( 0.5\) . Let

be the number of ancilla qubits with \( \epsilon'=2^{-d}\) . Then QPE obtains an estimate denoted by \( \widetilde{\theta}\) , which approximates \( \theta\) to precision \( \epsilon'\) with success probability \( 1-\delta\) . Note that \( p_0=\sin^2\frac{\theta}{2}\) , and

Using the fact that \( \left\lvert \sin (\widetilde{\theta}-\theta)\right\rvert \le \left\lvert \widetilde{\theta}-\theta\right\rvert \le \epsilon'\) , we have

Let \( \epsilon'\) be sufficiently small. Now if \( p_0(1-p_0)=\Omega(1)\) , we can choose \( \epsilon'=\mathcal{O}(\epsilon)\) , and the total complexity of QPE is \( \mathcal{O}(\epsilon^{-1})\) .

If \( p_0\) is small, then we should estimate \( p_0\) to multiplicative accuracy \( \epsilon\) instead. Use

we have \( \epsilon'=\sqrt{p_0}\epsilon\) . Therefore the runtime of QPE is \( \mathcal{O}(p_0^{-\frac12}\epsilon^{-1})\) . If \( p_0\) is to be estimated to precision \( \epsilon'\) using the Monte Carlo method, the number of samples would be \( \mathcal{N}=\mathcal{O}(p_0^{-1}\epsilon^{-2})\) .

So, the amplitude estimation method achieves quadratic speedup in the total complexity, but the circuit depth is increased to \( \mathcal{O}(\epsilon'^{-1}\delta^{-1})\) .

Example 30 (Amplitude estimation to accelerate Hadamard test)

Consider the circuit for the Hadamard test in (6) to estimate \( \operatorname{Re}\braket{\psi|U|\psi}\) . Let the initial state \( \ket{\psi}\) be prepared by a unitary \( U_{\psi}\) , then the following combined circuit

maps \( \ket{0}\ket{0^n}\) to

\[ \begin{equation} \ket{\psi_0}=\frac{1}{2}\ket{0}(\ket{\psi}+U\ket{\psi})+ \frac{1}{2}\ket{1}(\ket{\psi}-U\ket{\psi}):=\sqrt{p_0} \ket{\psi_{\mathrm{good}}}+\sqrt{1-p_0} \ket{\psi_{\mathrm{bad}}}, \end{equation} \]

(223)

and the goal is to estimate \( p_0\) . This also gives the implementation of the reflector \( R_{\psi_0}\) .

Note that \( R_{\mathrm{good}}\) can be implemented by simply reflecting against the signal qubit, i.e.,

\[ \begin{equation} R_{\mathrm{good}}=(I-2\ket{0}\bra{0})\otimes I_n=-Z\otimes I_n. \end{equation} \]

(224)

Then we can run QPE to the Grover unitary \( G=R_{\psi_0}R_{\mathrm{good}}\) to estimate \( p(0)\) , and the circuit depth is \( \mathcal{O}(\epsilon^{-1})\) .

In this section, we consider the solution of linear systems of equations

where \( A\in \mathbb{C}^{N\times N}\) is a non-singular matrix. Without loss of generality we assume \( A\) is Hermitian. Otherwise, we can solve a dilated Hermitian matrix, which enlarges the matrix dimension by a factor of \( 2\) (i.e., use one ancilla qubit)

and solve the enlarged problem

We assume \( b\in\mathbb{C}^N\) is a normalized vector and hence can be stored in a quantum state. More specifically, we have a unitary \( U_b\) such that (may require some work registers)

On classical computers, the solution is simply \( x=A^{-1}b\) . However, the solution vector is in general not a normalized vector, and hence cannot be directly stored as a quantum state. Therefore the goal of solving the quantum linear system problem (QLSP) is to find a quantum state \( \ket{\widetilde{x}}\) so that

The normalization constant \( \left\lVert A^{-1}\ket{b}\right\rVert \) should be recovered separately.

One useful application of QLSP solvers is to evaluate the many-body Green’s function [11], based on the quantum many-body Hamiltonian in (29). We will omit the details here.

5.3.1 Algorithmic description

The HHL algorithm [12], based on QPE, is the first quantum algorithm for solving QLSP. The algorithm can be summarized as follows. Let \( A\) have the following eigendecomposition

\[ \begin{equation} A\ket{v_j}=\lambda_j\ket{v_j}. \end{equation} \]

(230)

To simplify the analysis, we assume \( 0< \lambda_0\le \lambda_1\le …\le \lambda_{N-1}< 1\) and all eigenvalues have an exact \( d\) -bit representation. The analysis can also be generalized to the case when \( A\) has negative eigenvalues, but the interpretation of the result from QPE needs to be modified accordingly.

The matrix \( A\) can be queried using Hamiltonian simulation as \( U=e^{\mathrm{i} 2\pi A}\) . Then if the input state is already one of the eigenvectors, i.e., \(\ket{b}=\ket{v_j}\), then QPE can be applied to implement the mapping

\[ \begin{equation} U_{\mathrm{QPE}}\ket{0^d}\ket{v_j} =\ket{\lambda_j}\ket{v_j}. \end{equation} \]

(231)

In general, let the input state \(|b\rangle=\sum_{j} \beta_{j}\left|v_{j}\right\rangle\) be expanded using the eigendecomposition of \( A\) . Then by linearity,

\[ \begin{equation} U_{\mathrm{QPE}}\ket{0^d}\ket{b}=\sum_{j}\beta_j\ket{\lambda_j}\ket{v_j}. \end{equation} \]

(232)

Note that the unnormalized solution satisfies

\[ \begin{equation} A^{-1}\ket{b}=\sum_{j}\frac{\beta_j}{\lambda_j}\ket{v_j}, \end{equation} \]

(233)

so all we need to do is to use the information of the eigenvalue \( \ket{\lambda_j}\) stored in the ancilla register, and perform a controlled rotation to multiply the factor \( \lambda_j^{-1}\) to each \( \beta_j\) . To this end, we see that it is crucial to store all eigenvalues in the quantum computer coherently, as achieved by QPE. We would like to implement the following controlled rotation unitary (see (5.3.2))

\[ \begin{equation} U_{\mathrm{CR}}\ket{0}\ket{\lambda_j}=\left(\sqrt{1-\frac{C^{2}}{\widetilde{\lambda}_{j}^{2}}}\ket{0}+\frac{C}{\widetilde{\lambda}_{j}}\ket{1}\right)\ket{\lambda_j}. \end{equation} \]

(234)

where each \( \widetilde{\lambda}_{j}\) approximates \( \lambda_j\) .

Finally, perform uncomputation by applying \( U^{†}_{\mathrm{QPE}}\) , we convert the information from the ancilla register \( \ket{\lambda_j}\) back to \( \ket{0^d}\) . Therefore the quantum circuit for the HHL algorithm is in (14).

Note that through the uncomputation, the \( d\) ancilla qubits for storing the eigenvalues also becomes a working register. Discarding all working registers, the resulting unitary denoted by \( U_{\mathrm{HHL}}\) satisfies

\[ \begin{equation} U_{\mathrm{HHL}}\ket{0}\ket{b}=\sum_{j}\left(\sqrt{1-\frac{C^{2}}{\widetilde{\lambda}_{j}^{2}}}\ket{0}+\frac{C}{\widetilde{\lambda}_{j}}\ket{1}\right) \beta_j\ket{v_j}. \end{equation} \]

(235)

Finally, measuring the signal qubit (the only ancilla qubit left), if the outcome is \( 1\) , we obtain the (unnormalized) vector

\[ \begin{equation} \widetilde{x}=\sum_{j}\frac{C\beta_j}{\widetilde{\lambda}_{j}}\ket{v_j} \end{equation} \]

(236)

stored as a normalized state in the system register is

\[ \begin{equation} \ket{\widetilde{x}}=\frac{\widetilde{x}}{\left\lVert \widetilde{x}\right\rVert }\approx \ket{x}, \end{equation} \]

(237)

which is the desired approximate solution to QLSP. In particular, the constant \( C\) does not appear in the solution.

Remark 9 (Recovering the norm of the solution)

The HHL algorithm returns the solution to QLSP in the form of a normalized state \( \ket{x}\) stored in the quantum computer. In order to recover the magnitude of the unnormalized solution \( \left\lVert \widetilde{x}\right\rVert \) , we note that the success probability of measuring the signal qubit in (235) is

\[ \begin{equation} p(1)=\left\lVert \sum_{j}\frac{C\beta_j}{\widetilde{\lambda}_{j}}\ket{v_j}\right\rVert ^2=\left\lVert \widetilde{x}\right\rVert ^2. \end{equation} \]

(238)

Therefore sufficient repetitions of running the circuit in (235) and estimate \( p(1)\) , we can obtain an estimate of \( \left\lVert \widetilde{x}\right\rVert \) .

More general discussion of the readout problem of the HHL algorithm will be given in (12).

5.3.2 Implementation of controlled rotation

Proposition 3 (Controlled rotation given rotation angles)

Let \(0\le\theta<1\) has exact \( d\) -bit fixed point representation \(\theta=.\theta_{d-1}…\theta_0\) be its \(d\)-bit fixed point representation. Then there is a \( (d+1)\) -qubit unitary \(U_{\theta}\) such that

\[ \begin{equation} U_{\theta}: |0\rangle|\theta\rangle \mapsto (\cos (\pi\theta)|0\rangle+\sin (\pi\theta)|1\rangle)|\theta\rangle. \end{equation} \]

(239)

Proof

First (by e.g. Taylor expansion)

\[ \begin{equation} \exp \left(-\mathrm{i} \tau \sigma_{y}\right)=\left(\begin{array}{cc}{\cos (\tau)} & {-\sin (\tau)} \\ {\sin (\tau)} & {\cos (\tau)}\end{array}\right)=: R_y(2\tau). \end{equation} \]

(240)

Here \( R_y(\cdot)\) perform a single-qubit rotation around the \( y\) -axis. For any \( j\in [2^d]\) with its binary representation \( j=j_{d-1}… j_0\) , we have

\[ \begin{equation} j/2^d=(.j_{d-1}\cdots j_0). \end{equation} \]

(241)

So choose \( \tau=\pi(.j_{d-1}… j_0)\) , and define

\[ \begin{equation} U_{\theta}=\sum_{j\in [2^d]}\exp \left(-\mathrm{i} \pi(.j_{d-1}\cdots j_0) \sigma_{y}\right)\otimes \ket{j}\bra{j}. \end{equation} \]

(242)

Applying \(U_{\theta}\) to \(\ket{0}\ket{\theta}\) gives the desired results.

The quantum circuit for the controlled rotation circuit is

This is a sequence of single-qubit rotations on the signal qubit, each controlled by a single qubit.

In order to use the controlled rotation operation, we need to store the information of \( \lambda_j\) in term of an angle \( \theta_j\) . Let \( C>0\) be a lower bound to \( \lambda_0\) , so that \( 0<C/\lambda_j<1\) for all \( j\) . Define

\[ \begin{equation} \theta_j=\frac{1}{\pi}\arcsin (C / \lambda_j), \end{equation} \]

(243)

and \( \widetilde{\theta}_j\) be its \( d'\) -bit representation. Then

\[ \begin{equation} \sin \pi \widetilde{\theta}_j\equiv \frac{C}{\widetilde{\lambda}_j}\approx \frac{C}{\lambda_j}. \end{equation} \]

(244)

Again for simplicity we assume \( d'\) is large enough so that the error of the fixed point representation is negligible in this step. The mapping

\[ \begin{equation} U_{\mathrm{angle}}\ket{0^{d'-d}}\ket{\lambda_j}=\ket{\widetilde{\theta}_j} \end{equation} \]

(245)

can be implemented using classical arithmetics circuits in (2.8), which may require \( \operatorname{poly}(d')\) gates and an additional working register of \( \operatorname{poly}(d')\) qubits, which are not displayed here. Therefore the entire controlled rotation operation needed for the HHL algorithm is given by the circuit in (15).

Therefore through the uncomputation \( U^{†}_{\mathrm{angle}}\) , the \( d'-d\) ancilla qubits also become a working register. Discard the working register, and we obtain a unitary \( U_{\mathrm{CR}}\) satisfying

\[ \begin{equation} U_{\mathrm{CR}}\ket{0}\ket{\lambda_j}= \left(\cos(\pi \widetilde{\theta}_j)\ket{0}+\sin(\pi \widetilde{\theta}_j)\ket{1}\right)\ket{\lambda_j}= \left(\sqrt{1-\frac{C^{2}}{\widetilde{\lambda}_{j}^{2}}}\ket{0}+\frac{C}{\widetilde{\lambda}_{j}}\ket{1}\right)\ket{\lambda_j}. \end{equation} \]

(246)

This is the unitary used in the HHL circuit in (14).

5.3.3 Complexity analysis of the HHL algorithm

5.3.4 Additional considerations

As an application of the HHL algorithm, let us consider a toy problem of solving the Poisson’s equation in one-dimension with Dirichlet boundary conditions

We use the central finite difference formula to discretize the Laplacian operator on a uniform grid \( r_i=(i+1)h,i\in[N]\) and \( h=1/(N+1)\) . Let \( u_i=u(r_i), b_i=b(r_i)\) , then Poisson’s equation is discretized into a linear system of equations \( Au=b\) , with a tridiagonal matrix

Our goal here is not to address the quality of the spatial discretization, but to study the cost for solving the linear system. To this end we need to compute the condition number of \( A\) . We also assume the right hand vector \( b\) is already normalized so that \( \ket{b}=b\) .

Using (4) with \( a=2/h^2,b=-1/h^2\) , the largest eigenvalue of \( A\) is \( \lambda_{\max}=\left\lVert A\right\rVert \approx 4/h^2\) , and the smallest eigenvalue \( \lambda_{\min}=\left\lVert A^{-1}\right\rVert ^{-1}\approx \pi^2\) . So

The circuit depth of the HHL algorithm is \( \mathcal{O}(N^2/\epsilon)\) , and the worst case query complexity (using AA) is \( \mathcal{O}(N^4/\epsilon)\) . So when \( N\) is large, there is little benefit in employing the quantum computer to solve this problem.

Let us now consider solving a \( d\) -dimensional Poisson’s equation with Dirichlet boundary conditions

The grid is the Cartesian product of the uniform grid in 1D with \( N\) grid points per dimension and \( h=1/(N+1)\) . The total number of grid points is \( \mathcal{N}=N^d\) . After discretization, we obtain a linear system \( \mathcal{A}u=b\) , where

where \( I\) is an identity matrix of size \( N\) . Since \( \mathcal{A}\) is Hermitian and positive definite, we have \( \left\lVert \mathcal{A}\right\rVert \approx 4d/h^2\) , and \( \left\lVert \mathcal{A}^{-1}\right\rVert ^{-1}\approx d\pi^2\) . So \( \kappa(\mathcal{A})\approx 4/(h^2\pi^2)\approx \kappa(A)\) . Therefore the condition number is independent of the spatial dimension \( d\) .

The worst case query complexity of the HHL algorithm is \( \mathcal{O}(N^2/\epsilon)\) . So when the number of grid points per dimension \( N\) is fixed and the spatial dimension \( d\) increases, and if \( U=e^{\mathrm{i} \mathcal{A}\tau}\) can be implemented efficiently for some \( \tau\left\lVert \mathcal{A}\right\rVert <1\) with \( \operatorname{poly}(d)\) cost, then the HHL algorithm will have an advantage over classical solvers, for which each matrix-vector multiplication scales linearly with respect to \( \mathcal{N}\) and is therefore exponential in \( d\) .

Consider the solution of a time-dependent linear differential equation

Here \( b(t),x(t)\in \mathbb{C}^d\) and \( A(t)\in \mathbb{C}^{d\times d}\) . (267) can be used to represent a very large class of ordinary differential equations (ODEs) and spatially discretized partial differential equations (PDEs). For instance, if \( A(t)=-\mathrm{i} H(t)\) for some Hermitian matrix \( H(t)\) and \( b(t)=0\) , this is the time-dependent Schrödinger equation

A special case when \( H(t)\equiv H\) is often called the Hamiltonian simulation problem, and the solution can be written as

which can be viewed as the problem of evaluating the matrix function \( e^{-\mathrm{i} HT}\) . This will be discussed separately in later chapters.

In this section we consider the general case of (267), and for simplicity discretize the equation in time using the forward Euler method with a uniform grid \( t_k=k\Delta t\) where \( \Delta t=T/N,k=0,…,N\) . Let \( A_k=A(t_k), b_k=b(t_k)\) . The resulting discretized system becomes

which can be rewritten as

or more compactly as a linear systems of equations

Here \( I\) is the identity matrix of size \( d\) , and \( \mathbf{x}\in\mathbb{R}^{Nd}\) encodes the entire history of the states.

To solve (271) as a QLSP, the right hand side \( \mathbf{b}\) needs to be a normalized vector. This means we need to properly normalize \( x_0,b_k\) so that

In the limit \( \Delta t \to 0\) , this can be written as

which is not difficult to satisfy as long as \( x_0,b(t)\) can be prepared efficiently using unitary circuits and \( \left\lVert x_0\right\rVert ,\left\lVert b(t)\right\rVert =\Theta(1)\) .

To solve (271) using the HHL algorithm, we need to estimate the condition number of \( \mathcal{A}\) . Note that \( \mathcal{A}\) is a block-bidiagonal matrix and in particular is not Hermitian. So we need to use the dilation method in (226) and solve the corresponding Hermitian problem.

5.5.1 Scalar case

In order to estimate the condition number, for simplicity we first assume \( d=1\) (i.e., this is a scalar ODE problem), and \( A(t)\equiv a\in\mathbb{C}\) is a constant. Then

\[ \begin{equation} \mathcal{A}=\begin{pmatrix} 1 & 0 & 0 & \cdots & 0 & 0\\ -(1+\Delta t a) & 1 & 0 & \cdots & 0 & 0\\ 0 & -(1+\Delta t a) & 1 & \cdots & 0 & 0\\ \vdots & & & & &\vdots\\ 0& 0 & 0 & \cdots & -(1+\Delta t a) & 1 \end{pmatrix}\in\mathbb{C}^{N\times N}. \end{equation} \]

(275)

Let \( \xi=1+\Delta t a\) . When \( \operatorname{Re} a\le 0\) , the absolute stability condition of the forward Euler method requires \( \left\lvert 1+\Delta t a\right\rvert =\left\lvert \xi\right\rvert <1\) . In general we are interested in the regime \( \Delta t \left\lvert a\right\rvert \ll 1\) , and in particular \( \left\lvert 1+\Delta t a\right\rvert =\left\lvert \xi\right\rvert <2\) .

Proposition 5

For any \( A\in \mathbb{C}^{N\times N}\) , \( \left\lVert A\right\rVert ^2,(1/\left\lVert A^{-1}\right\rVert )^2\) are given by the largest and smallest eigenvalue of \( A^{†}A\) , respectively.

From (5), we need to first compute

\[ \begin{equation} \mathcal{A}^{†}\mathcal{A}=\begin{pmatrix} 1+\left\lvert \xi\right\rvert ^2 & -\overline{\xi} & \cdots & 0 & 0\\ -\xi & 1+\left\lvert \xi\right\rvert ^2 & -\overline{\xi} & \cdots & 0 & 0\\ 0 & -\xi & 1+\left\lvert \xi\right\rvert ^2 & \cdots & 0 & 0\\ \vdots & & & & &\vdots\\ 0& 0 & 0 & \cdots & 1+\left\lvert \xi\right\rvert ^2 & -\overline{\xi}\\ 0& 0 & 0 & \cdots & -\xi & 1 \end{pmatrix}. \end{equation} \]

(276)

Theorem 2 (Gershgorin circle theorem, see e.g. {[14, Theorem 7.2.1]})

Let \( A\in\mathbb{C}^{N\times N}\) with entries \( a_{ij}\) . For each \( i=1,…,N\) , define

\[ \begin{equation} R_i=\sum_{j\ne i}\left\lvert a_{ij}\right\rvert . \end{equation} \]

(277)

Let \( D(a_{ii},R_i)\subseteq \mathbb{C}\) be a closed disc centered at \( a_{ii}\) with radius \( R_i\) , which is called a Gershgorin disc. Then every eigenvalue of \( A\) lies within at least one of the Gershgorin discs \( D(a_{ii},R_i)\)

Since \( \mathcal{A}^{†}\mathcal{A}\) is Hermitian, we can restrict the Gershgorin discs to the real line so that \( D(a_{ii},R_i)\subseteq \mathbb{R}\) . Then Gershgorin discs of the matrix \( \mathcal{A}^{†}\mathcal{A}\) satisfy the bound

\[ \begin{equation} \begin{split} D(a_{11},R_1)&\subseteq \left[1+\left\lvert \xi\right\rvert ^2-\left\lvert \xi\right\rvert ,1+\left\lvert \xi\right\rvert ^2+\left\lvert \xi\right\rvert \right],\\ D(a_{ii},R_i)&\subseteq \left[1+\left\lvert \xi\right\rvert ^2-2\left\lvert \xi\right\rvert ,1+\left\lvert \xi\right\rvert ^2+2\left\lvert \xi\right\rvert \right], ~ i=2,\ldots, N-1\\ D(a_{NN},R_N)&\subseteq \left[1-\left\lvert \xi\right\rvert ,1+\left\lvert \xi\right\rvert \right]. \end{split} \end{equation} \]

(278)

Applying (2) we have

\[ \begin{equation} \lambda_{\max}(\mathcal{A}^{†}\mathcal{A}) \le (1+\left\lvert \xi\right\rvert )^2< 9. \end{equation} \]

(279)

for all values of \( a\) such that \( \left\lvert \xi\right\rvert <2\) .

To obtain a meaningful lower bound of \( \lambda_{\min}(\mathcal{A}^{†}\mathcal{A})\) , we need \( \operatorname{Re} a<0\) and hence \( \left\lvert \xi\right\rvert <1\) . Use the inequality

\[ \begin{equation} \sqrt{1+x}\le 1+\frac12 x, ~ x>-1, \end{equation} \]

(280)

we have

\[ \begin{equation} 1-\left\lvert \xi\right\rvert = 1-\sqrt{(1+\Delta t\operatorname{Re} a )^2+\Delta t^2(\operatorname{Im} a )^2}\ge -\Delta t\operatorname{Re} a-\frac{(\Delta t \left\lvert a\right\rvert )^2}{2}\ge -\frac{\Delta t}{2}\operatorname{Re} a, \end{equation} \]

(281)

when

\[ \begin{equation} \Delta t\left\lvert a\right\rvert <(-\operatorname{Re} a)/\left\lvert a\right\rvert \end{equation} \]

(282)

is satisfied. Then we have

\[ \begin{equation} 1>1-\left\lvert \xi\right\rvert \ge -\frac{\Delta t}{2}\operatorname{Re} a. \end{equation} \]

(283)

Then according to (2), under the condition (282),

\[ \begin{equation} \lambda_{\min}(\mathcal{A}^{†}\mathcal{A})\ge (1-\left\lvert \xi\right\rvert )^2\ge \frac{(\Delta t \operatorname{Re} a)^2}{4}. \end{equation} \]

(284)

Therefore

\[ \begin{equation} \left\lVert \mathcal{A}\right\rVert =\sqrt{\lambda_{\max}(\mathcal{A}^{†}\mathcal{A})} \le \sqrt{1+\left\lvert \xi\right\rvert ^2+2\left\lvert \xi\right\rvert }< 3, \end{equation} \]

(285)

and

\[ \begin{equation} \left\lVert \mathcal{A}^{-1}\right\rVert ^{-1}=\sqrt{\lambda_{\min}(\mathcal{A}^{†}\mathcal{A})}\ge \frac{\Delta t (-\operatorname{Re} a)}{2}. \end{equation} \]

(286)

As a result, when \( (-\operatorname{Re} a)=\Theta(1)\) , the condition number satisfies

\[ \begin{equation} \kappa(\mathcal{A})=\left\lVert \mathcal{A}\right\rVert \left\lVert \mathcal{A}^{-1}\right\rVert =\mathcal{O}\left(1/\Delta t\right). \end{equation} \]

(287)

In summary, for the scalar problem \( d=1\) and \( \operatorname{Re} a<0\) , the query complexity of the HHL algorithm is \( \mathcal{O}((\Delta t)^{-2}\epsilon^{-2})\) .

Example 31 (Growth of the condition number when \( \operatorname{Re} a\ge0\))

The Gershgorin circle theorem does not provide a meaningful bound of the condition number when \( \operatorname{Re} a>0\) and \( 1<\left\lvert \xi\right\rvert <2\) . This is a correct behavior. To see this, just consider \( a=1,b=0\) , and the solution should grow exponentially as \( x(T)=e^{T}\) . If \( \kappa(\mathcal{A})=\mathcal{O}\left((\Delta t) ^{-1}\right)\) holds and in particular is independent of the final \( T\) , then the norm of the solution can only grow polynomially in \( T\) , which is a contradiction. See (16) for an illustration. Note that when \( a=-1.0\) , \( \Delta t=0.1\) , the condition number is less than \( 20\) , which is smaller than the upper bound above, i.e., \( 3\times\frac{2}{\Delta t (-a)}=60\) .

Remark 14

It may be tempting to modify the \( (N,N)\) -th entry of \( \mathcal{A}\) to be \( 1+\left\lvert \xi\right\rvert ^2\) to obtain

\[ \begin{equation} \mathcal{G}=\begin{pmatrix} 1+\left\lvert \xi\right\rvert ^2 & -\overline{\xi} & \cdots & 0 & 0\\ -\xi & 1+\left\lvert \xi\right\rvert ^2 & -\overline{\xi} & \cdots & 0 & 0\\ 0 & -\xi & 1+\left\lvert \xi\right\rvert ^2 & \cdots & 0 & 0\\ \vdots & & & & &\vdots\\ 0& 0 & 0 & \cdots & 1+\left\lvert \xi\right\rvert ^2 & -\overline{\xi}\\ 0& 0 & 0 & \cdots & -\xi & 1+\left\lvert \xi\right\rvert ^2 \end{pmatrix}. \end{equation} \]

(288)

Here \( \mathcal{G}\) is a Toeplitz tridiagonal matrix satisfying the requirement of (4). The eigenvalues of \( \mathcal{G}\) take the form

\[ \begin{equation} \lambda_k=1+\left\lvert \xi\right\rvert ^2+2\left\lvert \xi\right\rvert \cos \frac{k\pi}{N+1}, ~ k=1,\ldots,N. \end{equation} \]

(289)

If we take the approximation \( \lambda_{\min}(\mathcal{A}^{†}\mathcal{A})\approx \lambda_{\min}(\mathcal{G})\) , we would find that (287) holds even when \( \operatorname{Re} a>0\) . This behavior is however incorrect, despite that the matrices \( \mathcal{A}^{†}\mathcal{A}\) and \( \mathcal{G}\) only differ by a single entry!

5.5.2 Vector case

5.5.3 Computing observables

As an application of the differential equation solver in 5.5, let us consider a toy problem of solving the heat equation in one-dimension with Dirichlet boundary conditions

After spatial discretization using the central finite difference method with \( d\) grid points, this becomes a linear ODE system

where \( A\in \mathbb{R}^{N\times N}\) is a tridiagonal matrix given by (259). After applying the forward Euler method and discretize the simulation time \( T\) into \( L\) intervals with \( \Delta t=T/L\) , we obtain a linear system of size \( NL\) . The eigenvalues of \( -A\) , denoted by \( \lambda_k\) , are all negative and satisfy

The absolute stability condition requires \( \left\lvert 1+\Delta t \lambda_k\right\rvert <1\) , or \( \Delta t<h^2/4=\mathcal{O}(N^{-2})\) , which implies \( L=\mathcal{O}(N^2)\) . Since \( A\) is Hermitian, we have \( \kappa(V)=1\) . So for \( T=\mathcal{O}(1)\) , the query complexity for solving the heat equation is \( \mathcal{O}(N^2\epsilon^{-2})\) , which is the same as solving Poisson’s equation.

Again, the potential advantage of the quantum solver only appears when solving the \( d\) -dimensional heat equation

This can be written as a linear system of equations

where \( \mathcal{A}\) is given in (266). The eigenvalues of \( \mathcal{A}\) are all negative. Note that \( \left\lVert \mathcal{A}\right\rVert =\Theta(dN^2)\) , then \( h=\mathcal{O}(d^{-1}N^{-2})\) , and the query complexity of the HHL solver is \( \mathcal{O}(dN^2 \epsilon^{-2})\) . This could potentially have an exponential advantage over classical solvers.

Consider the Hamiltonian simulation problem for \( H=H_1+H_2\) , where \( e^{-\mathrm{i} H_1 \Delta t}\) and \( e^{-\mathrm{i} H_2 \Delta t}\) can be efficiently computed at least for some \( \Delta t\) . In general \( [H_1,H_2]\ne 0\) , and the splitting of the evolution of \( H_1,H_2\) needs to be implemented via the Lie product formula

When taking \( L\) to be a finite number, and let \( \Delta t=t/L\) , this gives the simplest first order Trotter method with

Therefore to perform Hamiltonian simulation to time \( t\) , the error is

So to reach precision \( \epsilon\) in the operator norm, we need

This can be improved to the second order Trotter method (also called the symmetric Trotter splitting, or Strang splitting)

Following a similar analysis to the first order method, we find that to reach precision \( \epsilon\) we need

Higher order Trotter methods are also available, such as the \( p\) -th order Suzuki formula. The local truncation error is \( (\Delta t)^{p+1}\) . Therefore to reach precision \( \epsilon\) , we need

This is often written as \( L=\mathcal{O}(t^{1+o(1)}\epsilon^{-o(1)})\) as \( p\to \infty\) .

In this section we try to refine the error bounds in (313) by evaluating the preconstant explicitly. For simplicity we only focus on the first order Trotter formula. The Trotter propagator \( \widetilde{U}(t)=e^{-\mathrm{i} t H_1}e^{-\mathrm{i} t H_2}\) satisfies the equation

with initial condition \( \widetilde{U}(0)=I\) . By Duhamel’s principle, and let \( U(t)=e^{-\mathrm{i} t H}\) , we have

So we have

Now consider \( G(t)=[e^{-\mathrm{i} t H_1},H_2]e^{\mathrm{i} t H_1}=e^{-\mathrm{i} t H_1}H_2e^{\mathrm{i} t H_1}-H_2\) , which satisfies \( G(0)=0\) and

Hence

Plugging this back to (323), we have

In the last equality, we have used the relation \( \left\lVert [H_1,H_2]\right\rVert \le 2\nu^2\) with \( \nu=\max\{\left\lVert H_1\right\rVert ,\left\lVert H_2\right\rVert \}\) . Therefore (313) can be replaced by a sharper inequality

Here the first inequality is called the commutator norm error estimate, and the second inequality the operator norm error estimate.

For the transverse field Ising model with nearest neighbor interaction, we have \( \left\lVert H_1\right\rVert ,\left\lVert H_2\right\rVert =\mathcal{O}(n)\) , and hence \( \nu^2=\mathcal{O}(n^2)\) . On the other hand, since \( [Z_iZ_j,X_k]\ne 0\) only if \( k=i\) or \( k=j\) , the commutator bound satisfies \( \left\lVert [H_1,H_2]\right\rVert =\mathcal{O}(n)\) . Therefore to reach precision \( \epsilon\) , the scaling of the total number of time steps \( L\) with respect to the system size is \( \mathcal{O}(n^2/\epsilon)\) according to the estimate based on the operator norm, but is only \( \mathcal{O}(n/\epsilon)\) according to that based on the commutator norm.

For the particle in a potential, for simplicity consider \( d=1\) and the domain \( \Omega=[0,1]\) is discretized using a uniform grid of size \( N\) . For smooth and bounded potential, we have \( \left\lVert H_1\right\rVert =\mathcal{O}(N^2)\) , and \( \left\lVert V\right\rVert =\mathcal{O}(1)\) . Therefore the operator norm bound gives \( \nu^2=\mathcal{O}(N^4)\) . This is too pessimistic. Reexamining the second inequality of (327) shows that in this case, the error bound should be \( \mathcal{O}((\Delta t)^2\nu)\) instead of \( (\Delta t)^2\nu^2\) . So according to the operator norm error estimate, we have \( L=\mathcal{O}(N^2/\epsilon)\) . On the other hand, in the continuous space, for any smooth function \( \psi(r)\) , we have

So

Here we have used that \( \left\lVert V'\right\rVert =\left\lVert V''\right\rVert =\mathcal{O}(1)\) , and \( \left\lVert \psi'\right\rVert =\mathcal{O}(N)\) in the worst case scenario. Therefore \( \left\lVert [H_1,H_2]\right\rVert =\mathcal{O}(N)\) , and we obtain a significantly improved estimate \( L=\mathcal{O}(N/\epsilon)\) according to the commutator norm.

The commutator scaling of the Trotter error is an important feature of the method. We refer readers to [17] for analysis of the second order Trotter method, and [18, 19] for the analysis of the commutator scaling of high order Trotter methods.

The query model for sparse matrices is based on certain quantum oracles. In some scenarios, these quantum oracles can be implemented all the way to the elementary gate level.

Throughout the discussion we assume \( A\) is an \( n\) -qubit, square matrix, and

If the \( \left\lVert A\right\rVert _{\max}\ge1\) , we can simply consider the rescaled matrix \( \widetilde{A}/\alpha\) for some \( \alpha>\left\lVert A\right\rVert _{\max}\) .

To query the entries of a matrix, the desired oracle takes the following general form

In other words, given \( i,j\in[N]\) and a signal qubit \( 0\) , \( O_A\) performs a controlled rotation (controlling on \( i,j\) ) of the signal qubit, which encodes the information in terms of amplitude of \( \ket{0}\) .

However, the classical information in \( A\) is usually not stored natively in terms of such an oracle \( O_A\) . Sometimes it is more natural to assume that there is an oracle

where \( \widetilde{A}_{ij}\) is a \( d'\) -bit fixed point representation of \( A_{ij}\) , and the value of \( \widetilde{A}_{ij}\) is either computed on-the-fly with a quantum computer, or obtained through an external database. In either case, the implementation of \( \widetilde{O}_A\) may be challenging, and we will only consider the query complexity with respect to this oracle.

Using classical arithmetic operations, we can convert this oracle into an oracle

where \( 0\le \widetilde{\theta}_{ij}< 1\) , and \( \widetilde{\theta}_{ij}\) is a \( d\) -bit representation of \( \theta_{ij}=\arccos(A_{ij})/\pi\)  . This step may require some additional work registers not shown here.

Now using the controlled rotation in (3), the information of \( \widetilde{A}_{ij}\) , \( \widetilde{\theta}_{ij}\) has now been transferred to the phase of the signal qubit. We should then perform uncomputation and free the work register storing such intermediate information \( \widetilde{A}_{ij}\) , \( \widetilde{\theta}_{ij}\) . The procedure is as follows

From now on, we will always assume that the matrix entries of \( A\) can be queried using the phase oracle \( O_A\) or its variants.

The simplest example of block encoding is the following: assume we can find a \( (n+1)\) -qubit unitary matrix \( U\) (i.e., \( U\in\mathbb{C}^{2N\times 2N}\) ) such that

where \( *\) means that the corresponding matrix entries are irrelevant, then for any \( n\) -qubit quantum state \( \ket{b}\) , we can consider the state

and

Here the (unnormalized) state \( \ket{\perp}\) can be written as \( \ket{1}\ket{\psi}\) for some (unnormalized) state \( \ket{\psi}\) , that is irrelevant to the computation of \( A\ket{b}\) . In particular, it satisfies the orthogonality relation.

In order to obtain \( A\ket{b}\) , we need to measure the qubit \( 0\) and only keep the state if it returns \( 0\) . This can be summarized into the following quantum circuit:

Note that the output state is normalized after the measurement takes place. The success probability of obtaining \( 0\) from the measurement can be computed as

So the missing information of norm \( \left\lVert A\ket{b}\right\rVert \) can be recovered via the success probability \( p(0)\) if needed. We find that the success probability is only determined by \( A,\ket{b}\) , and is independent of other irrelevant components of \( U_A\) .

Note that we may not need to restrict the matrix \( U_A\) to be a \( (n+1)\) -qubit matrix. If we can find any \( (n+m)\) -qubit matrix \( U_A\) so that

Here each \( *\) stands for an \( n\) -qubit matrix, and there are \( 2^m\) block rows / columns in \( U_A\) . The relation above can be written compactly using the braket notation as

A necessary condition for the existence of \( U_A\) is that \( \left\lVert A\right\rVert \le 1\) . (Note: \( \left\lVert A\right\rVert _{\max}\le 1\) does not guarantee that \( \left\lVert A\right\rVert \le 1\) , see (26)). However, if we can find sufficiently large \( \alpha\) and \( U_A\) so that

Measuring the \( m\) ancilla qubits and all \( m\) -qubits return \( 0\) , we still obtain the normalized state \( \frac{A\ket{b}}{\left\lVert A\ket{b}\right\rVert }\) . The number \( \alpha\) is hidden in the success probability:

So if \( \alpha\) is chosen to be too large, the probability of obtaining all \( 0\) ’s from the measurement can be vanishingly small.

Finally, it can be difficult to find \( U_A\) to block encode \( A\) exactly. This is not a problem, since it is sufficient if we can find \( U_A\) to block encode \( A\) up to some error \( \epsilon\) . We are now ready to give the definition of block encoding in (1).

Assume we know each matrix element of the \( n\) -qubit matrix \( A_{ij}\) , and we are given an \( (n+m)\) -qubit unitary \( U_A\) . In order to verify that \( U_A\in\operatorname{BE}_{1,m}(A)\) , we only need to verify that

and \( U_A\) applied to any vector \( \ket{0^m,b}\) can be obtained via the superposition principle.

Therefore we may first evaluate the state \( U_A\ket{0^m,j}\) , and perform inner product with \( \ket{0^m,i}\) and verify the resulting the inner product is \( A_{ij}\) . We will also use the following technique frequently. Assume \( U_A=U_B U_C\) , and then

So we can evaluate the states \( U_B^{†}\ket{0^m,i},U_C\ket{0^m,j}\) independently, and then verify the inner product is \( A_{ij}\) . Such a calculation amounts to running the circuit (18), and if the ancilla qubits are measured to be \( 0^m\) , the system qubits return the normalized state \( \sum_{i}A_{ij}\ket{i}/\left\lVert \sum_{i}A_{ij}\ket{i}\right\rVert \) .

Example 35 (Random circuit block encoded matrix)

In some scenarios, we may want to construct a pseudo-random non-unitary matrix on quantum computers. Note that it would be highly inefficient if we first generate a dense pseudo-random matrix \( A\) classically and then feed it into the quantum computer using e.g. quantum random-access memory (QRAM). Instead we would like to work with matrices that are inherently easy to generate on quantum computers. This inspires the random circuit based block encoding matrix (RACBEM) model [20]. Instead of first identifying \( A\) and then finding its block encoding \( U_A\) , we reverse this thought process: we first identify a unitary \( U_A\) that is easy to implement on a quantum computer, and then ask which matrix can be block encoded by \( U_A\) .

(34) shows that in principle, any matrix \( A\) with \( \left\lVert A\right\rVert _2\le 1\) can be accessed via a \( (1,1,0)\) -block-encoding. In other words, \( A\) can be block encoded by an \( (n+1)\) -qubit random unitary \( U_A\) , and \( U_A\) can be constructed using only basic one-qubit unitaries and CNOT gates. The layout of the two-qubit operations can be designed to be compatible with the coupling map of the hardware. A cartoon is shown in (19), and an example is given in (20).

We now give a few examples of block encodings of more general sparse matrices. We start from a \( 1\) -sparse matrix, i.e., there is only one nonzero entry in each row or column of the matrix. This means that for each \( j\in [N]\) , there is a unique \( c(j)\in[N]\) such that \( A_{c(j),j}\ne 0\) , and the mapping \( c\) is a permutation. Then there exists a unitary \( O_c\) such that

The implementation of \( O_c\) may require the usage of some work registers that are omitted here. We also have

We assume the matrix entry \( A_{c(j),j}\) can be queried via

Now we construct \( U_A=(I\otimes O_c)O_A\) , and compute

This proves that \( U_A\in\operatorname{BE}_{1,1}(A)\) .

For a more general \( s\) -sparse matrix, WLOG we assume each row and column has exactly \( s\) nonzero entries (otherwise we can always treat some zero entries as nonzeros). For each column \( j\) , the row index for the \( \ell\) -th nonzero entry is denoted by \( c(j,\ell)\equiv c_{j,\ell}\) . For simplicity, we assume that there exists a unitary \( O_c\) such that

Here we assume \( s=2^\mathfrak{s}\) and the first register is an \( \mathfrak{s}\) -qubit register. A necessary condition for this query model is that \( O_c\) is reversible, i.e., we can have \( O_c^{†}\ket{\ell}\ket{c(j,\ell)}=\ket{\ell}\ket{j}\) . This means that for each row index \( i=c(j,\ell)\) , we can recover the column index \( j\) given the value of \( \ell\) . This can be satisfied e.g.

where \( \ell_0\) is a fixed number. This corresponds to a banded matrix. This assumption is of course somewhat restrictive. We shall discuss more general query models in (7.5).

Corresponding to (356), the matrix entries can be queried via

In order to construct a unitary that encodes all row indices at the same time, we define \( D=H^{\otimes \mathfrak{s}}\) (sometimes called a diffusion operator, which is a term originated from Grover’s search) satisfying

Consider \( U_A\) given by the circuit in (21). The measurement means that to obtain \( A\ket{b}\) , the ancilla register should all return the value \( 0\) .

So far we have considered general \( s\) -sparse matrices. Note that if \( A\) is a Hermitian matrix, its \( (\alpha, m, \epsilon)\) -block-encoding \( U_A\) does not need to be Hermitian. Even if \( \epsilon=0\) , we only have that the upper-left \( n\) -qubit block of \( U_A\) is Hermitian. For instance, even the block encoding of a Hermitian, diagonal matrix in (36) may not be Hermitian (exercise). On the other hand, there are indeed cases when \( U_A=U_A^{†}\) is indeed a Hermitian matrix, and hence the definition:

The Hermitian block encoding provides the simplest scenario of the qubitization process in (8.1).

If we query the oracle (356), the assumption that for each \( \ell\) the value of \( c(j,\ell)\) is unique for all \( j\) seems unnatural for constructing general sparse matrices. So we consider an altnerative method for construct the block encoding of a general sparse matrix as below.

Again WLOG we assume that each row / column has at most \( s=2^{\mathfrak{s}}\) nonzero entries, and that we have access to the following two \( (2n)\) -qubit oracles

Here \( r(i,\ell),c(j,\ell)\) gives the \( \ell\) -th nonzero entry in the \( i\) -th row and \( j\) -th column, respectively. It should be noted that although the index \( \ell\in[s]\) , we should expand it into an \( n\) -qubit state (e.g. let \( \ell\) take the last \( \mathfrak{s}\) qubits of the \( n\) -qubit register following the binary representation of integers). The reason for such an expansion, and that we need two oracles \( O_r,O_c\) will be seen shortly.

Similar to the discussion before, we need a diffusion operator satisfying

This can be implemented using Hadamard gates as

We assume that the matrix entries are queried using the following oracle using controlled rotations

where the rotation is controlled by both row and column indices. However, if \( A_{ij}=0\) for some \( i,j\) , the rotation can be arbitrary, as there will be no contribution due to the usage of \( O_r,O_c\) .

We remark that the quantum circuit in (22) is essentially the construction in [21, Lemma 48], which gives a \( (s,n+3)\) -block-encoding. The construction above slightly simplifies the procedure and saves two extra qubits (used to mark whether \( \ell\ge s\) ).

Next we consider the Hermitian block encoding of a \( s\) -sparse Hermitian matrix. Since \( A\) is Hermitian, we only need one oracle to query the location of the nonzero entries

Here \( c(j,\ell)\) gives the \( \ell\) -th nonzero entry in the \( j\) -th column. It can also be interpreted as the \( \ell\) -th nonzero entry in the \( i\) -th column. Again the first register needs to be interpreted as an \( n\) -qubit register. The diffusion operator is the same as in (365).

Unlike all discussions before, we introduce two signal qubits, and a quantum state in the computational basis takes the form \( \ket{a}\ket{i}\ket{b}\ket{j}\) , where \( a,b\in\{0,1\},i,j\in[N]\) . In other words, we may view \( \ket{a}\ket{i}\) as the first register, and \( \ket{b}\ket{j}\) as the second register. The \( (n+1)\) -qubit SWAP gate is defined as

To query matrix entries, we need access to the square root of \( A_{ij}\) as (note that act on the second single-qubit register)

The square root operation is well defined if \( A_{ij}\ge 0\) for all entries. If \( A\) has negative (or complex) entries, we first write \( A_{ij}=\left\lvert A_{ij}\right\rvert e^{\mathrm{i} \theta_{ij}},\theta_{ij}\in[0,2\pi)\) , and the square root is uniquely defined as \( \sqrt{A_{ij}}=\sqrt{\left\lvert A_{ij}\right\rvert }e^{\mathrm{i} \theta_{ij}/2}\) .

The quantum circuit in (23) is essentially the construction in [22]. The relation with quantum walks will be further discussed in (8.2).

We first introduce some heuristic idea behind qubitization. For any \( -1< \lambda\le 1\) , we can consider a \( 2\times 2\) rotation matrix,

where we have performed the change of variable \( \lambda=\cos\theta\) with \( 0\le \theta<\pi\) .

Now direct computation shows

Using the definition of Chebyshev polynomials (of first and second kinds, respectively)

we have

Note that if we can somehow replace \( \lambda\) by \( A\) , we immediately obtain a \( (1,1)\) -block-encoding for the Chebyshev polynomial \( T_k(A)\) ! This is precisely what qubitization aims at achieving, though there are some small twists.

In the simplest scenario, we assume that \( U_A\in\operatorname{HBE}_{1,m}(A)\) . Start from the spectral decomposition

we have that for each eigenstate \( \ket{v_i}\) ,

Here \( \ket{\widetilde{\perp}_i}\) is an unnormalized state that is orthogonal to all states of the form \( \ket{0^m}\ket{\psi}\) , i.e.,

where

is a projection operator.

Since the right hand side of (384) is a normalized state, we may also write

where \( \ket{\perp_i}\) is a normalized state.

Now if \( \lambda_i=\pm 1\) , then \( \mathcal{H}_i=\operatorname{span}\{\ket{0^m}\ket{v_i}\}\) is already an invariant subspace of \( U_A\) , and \( \ket{\perp_i}\) can be any state. Otherwise, use the fact that \( U_A=U_A^{†}\) , we can apply \( U_A\) again to both sides of (384) and obtain

Therefore \( \mathcal{H}_i=\operatorname{span}\{\ket{0^m}\ket{v_i},\ket{\perp_i}\}\) is an invariant subspace of \( U_A\) . Furthermore, the matrix representation of \( U_A\) with respect to the basis \( \mathcal{B}_i=\{\ket{0^m}\ket{v_i},\ket{\perp_i}\}\) is

i.e., \( U_A\) restricted to \( \mathcal{H}_i\) is a reflection operator. This also leads to the name “qubitization”, which means that each eigenvector \( \ket{v_i}\) is “qubitized” into a two-dimensional space \( \mathcal{H}_i\) .

In order to construct a block encoding for \( T_k(A)\) , we need to turn \( U_A\) into a rotation. For this note that \( \mathcal{H}_i\) is also an invariant subspace for the projection operator \( \Pi\) :

Similarly define \( Z_{\Pi}=2\Pi-1\) , since

\( Z_{\Pi}\) acts as a reflection operator restricted to each subspace \( \mathcal{H}_i\) . Then \( \mathcal{H}_i\) is the invariant subspace for the iterate

and

is the desired rotation matrix. Therefore

Since \( \{\ket{0^m}\ket{v_i}\}\) spans the range of \( \Pi\) , we have

i.e., \( O^k=(U_A Z_{\Pi})^k\) is a \( (1,m)\) -block-encoding of the Chebyshev polynomial \( T_k(A)\) .

In order to implement \( Z_{\Pi}\) , note that if \( m=1\) , then \( Z_{\Pi}\) is just the Pauli \( Z\) gate. When \( m>1\) , the circuit

returns \( \ket{1}\ket{0^m}\) if \( b=0^m\) , and \( -\ket{1}\ket{b}\) if \( b\ne 0^m\) . So this precisely implements \( Z_{\Pi}\) where the signal qubit \( \ket{1}\) is used as a work register. We may also discard the signal qubit, and resulting unitary is denoted by \( Z_{\Pi}\) .

In other words, the circuit in (24) implements the operator \( O\) . Repeating the circuit \( k\) times gives the \( (1,m+1)\) -block-encoding of \( T_k(A)\) .

Quantum walk is one of the major topics in quantum algorithms. Roughly speaking, there are two versions of quantum walks. The continuous time quantum walk is the closely analogous to its classical counterpart, i.e., continuous time random walk. The other version, the discrete time quantum walk, or Szegedy’s quantum walk [23] is not so obviously connected to the classical random walks. We will not introduce the motivations behind the continuous and discrete time random walks, and refer readers to [8, Chapter 16,17] for detailed discussions.

8.2.1 Basics of Markov chain

8.2.2 Block encoding of the discriminant matrix

Our first goal is to construct a Hermitian block encoding of \( D\) . Assume that we have access to an oracle \( O_P\) satisfying

\[ \begin{equation} O_P\ket{0^n}\ket{j}=\sum_{k}\sqrt{P_{jk}}\ket{k}\ket{j}. \end{equation} \]

(403)

Thanks to the stochasticity of \( P\) , the right hand side is already a normalized vector, and no additional signal qubit is needed.

We also introduce the \( n\) -qubit SWAP operator Swap operator:

\[ \begin{equation} \operatorname{SWAP} \ket{i}\ket{j}=\ket{j}\ket{i}, \end{equation} \]

(404)

which swaps the value of the two registers in the computational basis, and can be directly implemented using \( n\) two-qubit SWAP gates.

We claim that the following circuit gives \( U_D\in \operatorname{HBE}_{1,n}(D)\) .

Proposition 11

(25) defines \( U_D\in\operatorname{HBE}_{1,n}(D)\) .

Proof

Clearly \( U_D\) is unitary and Hermitian. Now we compute as before

\[ \begin{equation} \ket{0^{n}}\ket{j}\xrightarrow{O_P} \sum_{k}\sqrt{P_{jk}}\ket{k}\ket{j}\xrightarrow{ \operatorname{SWAP} } \sum_{k}\sqrt{P_{jk}}\ket{j}\ket{k}. \end{equation} \]

(405)

Meanwhile

\[ \begin{equation} \ket{0^{n}}\ket{i}\xrightarrow{O_P} \sum_{k'}\sqrt{P_{ik'}}\ket{k'}\ket{i}. \end{equation} \]

(406)

So the inner product gives

\[ \begin{equation} \bra{0^{n}}\bra{i} U_D \ket{0^n}\ket{j}=\sum_{k,k'}\sqrt{P_{ik'}P_{jk}} \delta_{j,k'} \delta_{i,k}=\sqrt{P_{ij}P_{ji}}=D_{ij}. \end{equation} \]

(407)

This proves the claim.

8.2.3 Szegedy’s quantum walk

For a Markov chain defined on a graph \( G=(V,E)\) , Szegedy’s quantum walk implements a qubitization of the discriminant matrix \( D\) in (399). Let \( U_D\) be the Hermitian block encoding defined by the circuit in (25), we may readily plug it into (24), and obtain the circuit in (26) denoted by \( O_Z\) .

Let the eigendecomposition of \( D\) be denoted by

\[ \begin{equation} D\ket{v_i}=\lambda_i \ket{v_i}. \end{equation} \]

(408)

For each \( \ket{v_i}\) , the associated basis in the 2-dimensional subspace is \( \mathcal{B}_i=\{\ket{0^m}\ket{v_i},\ket{\perp_i}\}\) . Then the qubitization procedure gives

\[ \begin{equation} [O_Z]_{\mathcal{B}_i}=\begin{pmatrix} \lambda_i & -\sqrt{1-\lambda_i^2}\\ \sqrt{1-\lambda_i^2} & \lambda_i \end{pmatrix}. \end{equation} \]

(409)

The eigenvalues of \( O_Z\) in the \( 2\times 2\) matrix block are

\[ \begin{equation} e^{\pm\mathrm{i}\arccos(\lambda_i)}. \end{equation} \]

(410)

This relation is important for the following reasons. By (10), if a Markov chain is reversible and ergodic, the eigenvalues of \( D\) and \( P\) are the same. In particular, the largest eigenvalue of \( D\) is unique and is equal to \( 1\) , and the second largest eigenvalue of \( D\) is \( 1-\delta\) , where \( \delta>0\) is called the spectral gap. Since \( \arccos(1)=0\) , and \( \arccos(1-\delta)\approx \sqrt{2 \delta}\) , we find that the spectral gap of \( O_Z\) on the unit circle is in fact \( \mathcal{O}(\sqrt{\delta})\) instead of \( \mathcal{O}(\delta)\) . This is called the spectral gap amplification, which leads to e.g. the quadratic quantum speedup of the hitting time.

Example 37 (Determining whether there is a marked vertex in a complete graph)

Let \( G=(V,E)\) be a complete, graph of \( N=2^n\) vertices. We would like to distinguish the following two scenarios:

1. All vertices are the same, and the random walk is given by the transition matrix

   \[ \begin{equation} P=\frac{1}{N} e e^{\top}, ~ e=(1,\ldots,1)^{\top}. \end{equation} \]

   (411)

2. There is one marked vertex. Without loss of generality we may assume this is the \( 0\) -th vertex (of course we do not have access to this information). In this case, the transition matrix is

   \[ \begin{equation} \widetilde{P}_{ij}=\left\{ \begin{array}{ll} \delta_{ij}, & i=0,\\ P_{ij}, & i>0. \end{array}\right. \end{equation} \]

   (412)

In other words, in the case (2), the random walk will stop at the marked index. The transition matrix can also be written in the block partitioned form as

\[ \begin{equation} \widetilde{P}=\begin{pmatrix} 1 & 0\\ \frac{1}{N} \widetilde{e} & \frac{1}{N} \widetilde{e}\widetilde{e}^{\top} \end{pmatrix}. \end{equation} \]

(413)

Here \( \widetilde{e}\) is an all \( 1\) vector of length \( N-1\) .

For the random walk defined by \( P\) , the stationary state is \( \pi=\frac{1}{N}e\) , and the spectral gap is \( 1\) . For the random walk defined by \( \widetilde{P}\) , the stationary state is \( \widetilde{\pi}=(1,0,…,0)^{\top}\) , and the spectral gap of is \( \delta=N^{-1}\) . Starting from the uniform state \( \pi\) , the probability distribution after \( k\) steps of random walk is \( \pi^{\top}\widetilde{P}^{k}\) . This converges to the stationary state of \( \widetilde{P}\) , and hence reach the marked vertex after \( \mathcal{O}(N)\) steps of walks (exercise).

These properties are also inherited by the discriminant matrices, with \( D=P\) and

\[ \begin{equation} \widetilde{D}=\begin{pmatrix} 1 & 0\\ 0 & \frac{1}{N} \widetilde{e}\widetilde{e}^{\top} \end{pmatrix}. \end{equation} \]

(414)

To distinguish the two cases, we are given a Szegedy quantum walk operator called \( O\) , which can be either \( O_Z\) or \( \widetilde{O}_Z\) , which is associated with \( D,\widetilde{D}\) , respectively. The initial state is

\[ \begin{equation} \ket{\psi_0}=\ket{0^n}(H^{\otimes n} \ket{0^n}). \end{equation} \]

(415)

Our strategy is to measure the expectation

\[ \begin{equation} m_k=\braket{\psi_0|O^k |\psi_0}, \end{equation} \]

(416)

which can be obtained via Hadamard’s test.

Before determining the value of \( k\) , first notice that if \( O=O_Z\) , then \( O_{Z}\ket{\psi_0}=\ket{\psi_0}\) . Hence \( m_k=1\) for all values of \( k\) .

On the other hand, if \( O=\widetilde{O}_Z\) , we use the fact that \( \widetilde{D}\) only has two nonzero eigenvalues \( 1\) and \( (N-1)/N=1-\delta\) , with associated eigenvectors denoted by \( \ket{\widetilde{\pi}}\) and \( \ket{\widetilde{v}}=\frac{1}{\sqrt{N-1}}(0,1,1…,1)^{\top}\) , respectively. Furthermore,

\[ \begin{equation} \ket{\psi_0}=\frac{1}{\sqrt{N}}\ket{0^n}\ket{\widetilde{\pi}}+\sqrt{\frac{N-1}{N}}\ket{0^n}\ket{\widetilde{v}}. \end{equation} \]

(417)

Due to qubitization, we have

\[ \begin{equation} \widetilde{O}_Z^k\ket{\psi_0}=\frac{1}{\sqrt{N}}\ket{0^n}T_k(1)\ket{\widetilde{\pi}}+\sqrt{\frac{N-1}{N}}\ket{0^n}T_k(1-\delta)\ket{\widetilde{v}}+\ket{\perp}, \end{equation} \]

(418)

where \( \ket{\perp}\) is an unnormalized state satisfying \( (\ket{0^n}\bra{0^n})\otimes I_n \ket{\perp}=0\) . Then using \( T_k(1)=1\) for all \( k\) , we have

\[ \begin{equation} m_k=\frac{1}{N}+\left(1-\frac{1}{N}\right)T_k(1-\delta). \end{equation} \]

(419)

Use the fact that \( T_k(1-\delta)=\cos(k\arccos(1-\delta))\) , in order to have \( T_k(1-\delta)\approx 0\) , the smallest \( k\) satisfies

\[ \begin{equation} k\approx \frac{\pi}{2\arccos(1-\delta)}\approx \frac{\pi}{2\sqrt{2\delta}}=\frac{\pi\sqrt{N}}{2\sqrt{2}}. \end{equation} \]

(420)

Therefore taking \( k=\lceil\frac{\pi\sqrt{N}}{2\sqrt{2}} \rceil\) , we have \( m_k\approx 1/N\) . Running Hadamard’s test to constant accuracy allows us to distinguish the two scenarios.

Remark 17 (Without using the Hadamard test)

Alternatively, we may evaluate the success probability of obtaining \( 0^n\) in the ancilla qubits, i.e.,

\[ \begin{equation} p(0^n)=\left\lVert (\ket{0^n}\bra{0^n}\otimes I_n)O^k\ket{\psi_0}\right\rVert ^2. \end{equation} \]

(421)

When \( O=O_Z\) , we have \( p(0^n)=1\) with certainty. When \( O=\widetilde{O}_Z\) , according to (418),

\[ \begin{equation} p(0^n)=\frac{1}{N}+\left(1-\frac{1}{N}\right)T^2_k(1-\delta). \end{equation} \]

(422)

So running the problem with \( k=\lceil\frac{\pi\sqrt{N}}{2\sqrt{2}} \rceil\) , we can distinguish between the two cases.

Remark 18 (Comparison with Grover’s search)

It is natural to draw comparisons between Szegedy’s quantum walk and Grover’s search. The two algorithms make queries to different oracles, and both yield quadratic speedup compared to the classical algorithms. The quantum walk is slightly weaker, since it only tells whether there is one marked vertex or not. On the other hand, Grover’s search also finds the location of the marked vertex. Both algorithms consist of repeated usage of the product of two reflectors. The number of iterations need to be carefully controlled. Indeed, choosing a polynomial degree four times as large as (420) would result in \( m_k\approx 1\) for the case with a marked vertex.

Remark 19 (Comparison with QPE)

Another possible solution of the problem of finding the marked vertex is to perform QPE on the Szegedy walk operator \( O\) (which can be \( O_Z\) or \( \widetilde{O}_Z\) ). The effectiveness of the method rests on the spectral gap amplification discussed above. We refer to  [8, Chapter 17] for more details.

8.2.4 Comparison with the original version of Szegedy’s quantum walk

In practical applications, we are often not interested in constructing the Chebyshev polynomial of \( A\) , but a linear combination of Chebyshev polynomials. For instance, the matrix inversion problem can be solved by expanding \( f(x)=x^{-1}\) using a linear combination of Chebyshev polynomials, on the interval \( [-1,-\kappa^{-1}]\cup [\kappa^{-1},1]\) . Here we assume \( \left\lVert A\right\rVert =1\) and \( \kappa=\left\lVert A\right\rVert \left\lVert A^{-1}\right\rVert \) is the condition number. One way to implement this is via a quantum primitive called the linear combination of unitaries (LCU).

Let \(T=\sum_{i=0}^{K-1} \alpha_i U_i\) be the linear combination of unitary matrices \( U_i\) . For simplicity let \(K=2^a\), and \( \alpha_i>0\) (WLOG we can absorb the phase of \( \alpha_i\) into the unitary \( U_i\) ). Then

implements the selection of \(U_i\) conditioned on the value of the \(a\)-qubit ancilla states (also called the control register). \(U\) is called a select oracle.

Let \(V\) be a unitary operation satisfying

and \(V\) is called the prepare oracle. The 1-norm of the coefficients is given by

In the matrix form

where the first basis is \(\ket{0^m}\), and all other basis functions are orthogonal to it. Then

Then \( T\) can be implemented using the unitary given in (1) (called the LCU lemma).

The LCU Lemma is a useful quantum primitive, as it states that the number of ancilla qubits needed only depends logarithmically on \( K\) , the number of terms in the linear combination. Hence it is possible to implement the linear combination of a very large number of terms efficiently. From a practical perspective, the select and prepare oracles uses multi-qubit controls, and can be difficult to implement. If implemented directly, the number of multi-qubit controls again depends linearly on \( K\) and is not desirable. Therefore an efficient implementation using LCU (in terms of the gate complexity) also requires additional structures in the prepare and select oracles.

If we apply \(W\) to \(\ket{0^a}\ket{\psi}\) and measure the ancilla qubits, then the probability of obtaining the outcome \(0^a\) in the ancilla qubits (and therefore obtaining the state \(T\ket{\psi}/\left\lVert T\ket{\psi}\right\rVert \) in the system register) is \(\left(\left\lVert T\ket{\psi}\right\rVert /\left\lVert \alpha\right\rVert _1\right)^2\). The expected number of repetition needed to succeed is \(\left(\left\lVert \alpha\right\rVert _1/\left\lVert T\ket{\psi}\right\rVert \right)^2\). Now we demonstrate that using amplitude amplification (AA) in (3.3), this number can be reduced to \(\mathcal{O}\left(\left\lVert \alpha\right\rVert _1/\left\lVert T\ket{\psi}\right\rVert \right)\).

Example 38 (Linear combination of two matrices)

Let \( U_A,U_B\) be two \( n\) -qubit unitaries, and we would like to construct a block encoding of \( T=U_A+U_B\) .

There are two terms in total, so one ancilla qubit is needed. The prepare oracle needs to implement

\[ \begin{equation} V\ket{0}=\frac{1}{\sqrt{2}}(\ket{0}+\ket{1}), \end{equation} \]

(440)

so this is the Hadamard gate. The circuit is given by (27), which constructs \( W\in \operatorname{BE}_{\sqrt{2},1}(T)\) .

A special case is the linear combination of two block encoded matrices. Given two \( n\) -qubit matrices \( A,B\) , for simplicity let \( U_A\in\operatorname{BE}_{1,m}(A),U_B\in\operatorname{BE}_{1,m}(B)\) . We would like to construct a block encoding of \( T=A+B\) . The circuit is given by (28), which constructs \( W\in \operatorname{BE}_{\sqrt{2},1+m}(T)\) . This is also an example of a linear combination of non-unitary matrices.

In (8.1) we assume that \( U_A=U_A^{†}\) to block encode a Hermitian matrix \( A\) . For instance, \( s\) -sparse Hermitian matrices, such Hermitian block encodings can be constructed following the recipe in (7.5). However, this can come at the expense of requiring additional structures and oracles. In general, the block encoding of a Hermitian matrix may not be Hermitian itself. In this section we demonstrate that the strategy of qubitization can be modified to accommodate general block encodings.

Again start from the eigendecomposition (383), we apply \( U_A\) to \( \ket{0^m}\ket{v_i}\) and obtain

where \( \ket{\perp_i'}\) is a normalized state satisfying \( \Pi\ket{\perp_i'}=0\) .

Since \( U_A\) block-encodes a Hermitian matrix \( A\) , we have

which implies that there exists another normalized state \( \ket{\perp_i}\) satisfying \( \Pi\ket{\perp_i}=0\) and

Now apply \( U_A\) to both sides of (453), we obtain

which gives

Define

and the associated two-dimensional subspaces \( \mathcal{H}_i= \operatorname{span} B_i , \mathcal{H}'_i= \operatorname{span} B_i' \) , we find that \( U_A\) maps \( \mathcal{H}_i\) to \( \mathcal{H}_i'\) . Correspondingly \( U_A^{†}\) maps \( \mathcal{H}_i'\) to \( \mathcal{H}_i\) .

Then (451,455) give the matrix representation

Similar calculation shows that

Meanwhile both \( \mathcal{H}_i\) and \( \mathcal{H}_i'\) are the invariant subspaces of the projector \( \Pi\) , with matrix representation

Therefore

Hence \( \mathcal{H}_i\) is an invariant subspace of \( \widetilde{O}=U_A^{†}Z_{\Pi}U_A Z_{\Pi}\) , with matrix representation

Repeating \( k\) times, we have

Since any vector \( \ket{0^m}\ket{\psi}\) can be expanded in terms of the eigenvectors \( \ket{0^m}\ket{v_i}\) , we have

Therefore if we would like to construct an even order Chebyshev polynomial \( T_{2k}(A)\) , the circuit \( (U_A^{†}Z_{\Pi}U_A Z_{\Pi})^{k}\) straightforwardly gives a \( (1,m)\) -block-encoding.

In order to construct the block-encoding of an odd polynomial \( T_{2k+1}(A)\) , we note that

Using the fact that \( \mathcal{B}_i,\mathcal{B}_i'\) share the common basis \( \ket{0^m}\ket{v_i}\) , we still have the block-encoding

Therefore \( U_A Z_{\Pi}(U_A^{†}Z_{\Pi}U_A Z_{\Pi})^{k}\) is a \( (1,m)\) -block-encoding of \( T_{2k+1}(A)\) .

In summary, the block-encoding of \( T_{l}(A)\) is given by applying \( U_A Z_{\Pi}\) and \( U^{†}_A Z_{\Pi}\) alternately. If \( l=2k\) , then there are exactly \( k\) such pairs. The quantum circuit for each pair \( \widetilde{O}\) is

Otherwise if \( l=2k+1\) , then there is an extra \( U_A Z_{\Pi}\) . The effect is to map each eigenvector \( \ket{0^m}\ket{v_i}\) back and forth between the two-dimensional subspaces \( \mathcal{H}_i\) and \( \mathcal{H}_i'\) . In (8.6), we shall see that the separate treatment even/odd polynomials will play a more prominent role.

Now that we have obtained \( O_k\in\operatorname{BE}_{1,m}(T_k(A))\) for all \( k\) , we can use the LCU lemma again to construct block encodings for linear combination of Chebyshev polynomials. We omit the details here.

Let us briefly recap what we have done so far: (1) construct a block encoding of an Hermitian matrix \( A\) (the block encoding matrix itself can be non-Hermitian); (2) use qubitization to block encode \( T_k(A)\) ; (3) use LCU to block encode an arbitrary polynomial function of \( A\) (up to a subnormalization factor). This framework is conceptually appealing, but the practical implementation of the select and prepare oracles are by no means straightforward, and can come at significant costs.

8.5.1 Hermitian block encoding

Note that \( \mathrm{i} Z=e^{\mathrm{i} \frac{\pi}{2}Z}\) , if \( A\) is given by a Hermitian block encoding \( U_A\) , the block encoding of the Chebyshev polynomial in (8.1) can be written as

\[ \begin{equation} O^d=(-\mathrm{i})^d \prod_{j=1}^d (U_A e^{\mathrm{i} \frac{\pi}{2} Z_{\Pi}}). \end{equation} \]

(466)

This is a special case of the following representation. Note that \( (-\mathrm{i})^d\) is an irrelevant phase factor and can be discarded.

\[ \begin{equation} U_{\widetilde{\Phi}}=e^{\mathrm{i} \widetilde{\phi}_0 Z_{\Pi}} \prod_{j=1}^d (U_A e^{\mathrm{i} \widetilde{\phi}_j Z_{\Pi}}). \end{equation} \]

(467)

The representation (467) is called a quantum eigenvalue transformation (QET).

Due to qubitization, \( U_{\widetilde{\Phi}}\) should block encode some matrix polynomial of \( A\) . We first state the following theorem without proof.

Theorem 3 (Quantum signal processing, a simplified version)

Consider

\[ \begin{equation} U_A=\begin{pmatrix} x & \sqrt{1-x^2}\\ \sqrt{1-x^2} & -x \end{pmatrix}. \end{equation} \]

(468)

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

\[ \begin{equation} U_{\widetilde{\Phi}} := e^{\mathrm{i} \widetilde{\phi}_0 Z} \prod_{j=1}^d (U_A e^{\mathrm{i} \widetilde{\phi}_j Z}) = \begin{pmatrix} P(x) & *\\ * & * \end{pmatrix}, \end{equation} \]

(469)

where \( P\in\mathbb{C}\) satisfy

1. \( \deg(P) \leq d\) ,
2. \( P\) has parity \( (d \bmod 2)\) ,
3. \( |P(x)|\le 1, \forall x \in [-1, 1]\) .

Also define \( -\widetilde{\Phi} := (-\widetilde{\phi}_0, …, -\widetilde{\phi}_d) \in \mathbb{R}^{d+1}\) , then

\[ \begin{equation} U_{-\widetilde{\Phi}} := e^{-\mathrm{i} \widetilde{\phi}_0 Z} \prod_{j=1}^d (U_A e^{-\mathrm{i} \widetilde{\phi}_j Z}) = U^*_{\widetilde{\Phi}}= \begin{pmatrix} P^*(x) & *\\ * & * \end{pmatrix}. \end{equation} \]

(470)

Here \( P^*(x)\) is the complex conjugation of \( P(x)\) , and \( U^*_{\widetilde{\Phi}}\) is the complex conjugation (without transpose) of \( U_{\widetilde{\Phi}}\) .

Remark 22

This theorem can be proved inductively. However, this is a special case of the quantum signal processing in (6), so we will omit the proof here. In fact, (6) will also state the converse of the result, which describes precisely the class of matrix polynomials that can be described by such phase factor modulations. In (3), the condition (1) states that the polynomial degree is upper bounded by the number of \( U_A\) ’s, and the condition (3) is simply a consequence of that \( U_{\widetilde{\Phi}}\) is a unitary matrix. The condition (2) is less obvious, but should not come at a surprise, since we have seen the need of treating even and odd polynomials separately in the case of qubitization with a general block encoding. (470) can be proved directly by taking the complex conjugation of \( U_{\widetilde{\Phi}}\) .

Following the qubitization procedure, we immediately have (4).

Theorem 4 (Quantum eigenvalue transformation with Hermitian block encoding)

Let \( U_A\in\operatorname{HBE}_{1,m}(A)\) . Then for any \( \widetilde{\Phi} := (\widetilde{\phi}_0, …, \widetilde{\phi}_d) \in \mathbb{R}^{d+1}\) ,

\[ \begin{equation} U_{\widetilde{\Phi}} = e^{\mathrm{i} \widetilde{\phi}_0 Z_{\Pi}} \prod_{j=1}^d (U_A e^{\mathrm{i} \widetilde{\phi}_j Z_{\Pi}}) = \begin{pmatrix} P(A) & *\\ * & * \end{pmatrix}\in \operatorname{BE}_{1,m}(P(A)), \end{equation} \]

(471)

where \( P\in\mathbb{C}\) satisfies the requirements in (3).

Using (4), we may construct the block encoding of a matrix polynomial without invoking LCU. The cost is essentially the same as block encoding a Chebyshev polynomial.

In order to implement \( e^{\mathrm{i} \phi Z_{\Pi}}\) , we note that the quantum circuit denoted by \( \operatorname{CR} _{\phi}\) is in (30) returns \( e^{\mathrm{i}\phi}\ket{0}\ket{0^m}\) if \( b=0^m\) , and \( e^{-\mathrm{i}\phi}\ket{0}\ket{b}\) if \( b\ne 0^m\) . So omitting the signal qubit, this is precisely \( e^{\mathrm{i} \phi Z_{\Pi}}\) .

Therefore, if \( A\) is given by a Hermitian block encoding \( U_A\) , we can follow the argument in (8.1) and construct the following unitary The corresponding quantum circuit is in (31), which uses one extra ancilla qubit. When measuring the \( (m+1)\) ancilla qubits and obtain \( \ket{0}\ket{0^m}\) , the corresponding (unnormalized) state in the system register is \( P(A)\ket{\psi}\) .

The QET described by the circuit in (31) generally constructs a block encoding of \( P(A)\) for some complex polynomial \( P\) . In practical applications (such as those later in this chapter), we would like to construct a block encoding of \( P_{\operatorname{Re}}(A)\equiv (\operatorname{Re} P)(A)=\frac12(P(A)+P^*(A))\) instead. Below we demonstrate that a simple modification of (31) allows us to achieve this goal.

To this end, we use (470). Qubitization allows us to construct

\[ \begin{equation} U_{-\widetilde{\Phi}} = e^{-\mathrm{i} \widetilde{\phi}_0 Z_{\Pi}} \prod_{j=1}^d (U_A e^{-\mathrm{i} \widetilde{\phi}_j Z_{\Pi}}) = \begin{pmatrix} P^*(A) & *\\ * & * \end{pmatrix}. \end{equation} \]

(472)

So all we need is to negate all phase factors in \( \widetilde{\Phi}\) . In order to implement \( \operatorname{CR} _{-\phi}\) , we do not actually need to implement a new circuit. Instead we may simply change the signal qubit from \( \ket{0}\) to \( \ket{1}\) :

which returns \( e^{-\mathrm{i}\phi}\ket{1}\ket{0^m}\) if \( b=0^m\) , and \( e^{\mathrm{i}\phi}\ket{1}\ket{b}\) if \( b\ne 0^m\) . In other words, the circuit for \( U_{P^*(A)}\) and \( U_{P(A)}\) are exactly the same except that the input signal qubit is changed from \( \ket{0}\) to \( \ket{1}\) .

Now we claim the circuit in (32) implements a block encoding \( U_{P_{\operatorname{Re}}(A)}\in \operatorname{BE}_{1,m+1}(P_{\operatorname{Re}}(A))\) . This circuit can be viewed as an implementation of the linear combination of unitaries \( \frac12 (U_{P^*(A)}+U_{P(A)})\) .

\[ \]

To verify this, we may evaluate

\[ \begin{equation} \begin{split} \ket{0}\ket{0^m}\ket{\psi} &\xrightarrow{H\otimes I_{m+n}} \frac{1}{\sqrt{2}}(\ket{0}+\ket{1})\ket{0^m}\ket{\psi}\\ &\xrightarrow{U_{P(A)}} \frac{1}{\sqrt{2}}\ket{0}(\ket{0^m}P(A)\ket{\psi}+\ket{\perp}) + \frac{1}{\sqrt{2}}\ket{1}(\ket{0^m}P^*(A)\ket{\psi}+\ket{\perp'})\\ &\xrightarrow{H\otimes I_{m+n}} \ket{0}\left(\ket{0^m}\frac{P(A)+P^*(A)}{2}\ket{\psi}\right)+\ket{\widetilde{\perp}}\\ &=\ket{0}\ket{0^m}P_{\operatorname{Re}}(A)\ket{\psi}+\ket{\widetilde{\perp}} \end{split} \end{equation} \]

(473)

Here \( \ket{\perp},\ket{\perp'}\) are two \( (m+n)\) -qubit state orthogonal to any state \( \ket{0^m}\ket{x}\) , while \( \ket{\widetilde{\perp}}\) is a \( (m+n+1)\) -qubit state orthogonal to any state \( \ket{0}\ket{0^m}\ket{x}\) . In other words, by measuring all \( (m+1)\) ancilla qubits and obtain \( 0^{m+1}\) , the corresponding (unnormalized) state in the system register is \( P_{\operatorname{Re}}(A)\ket{\psi}\) .

8.5.2 General block encoding

If \( A\) is given by a general block encoding \( U_A\) , the quantum eigenvalue transformation should consist of an alternating sequence of \( U_A,U_A^{†}\) gates. The circuit is given by (33), and the corresponding block encoding is described in (5). Note that the Hermitian block encoding becomes a special case with \( U_A=U_A^{†}\) .

Theorem 5 (Quantum eigenvalue transformation with general block encoding)

Let \( U_A\in\operatorname{BE}_{1,m}(A)\) . Then for any \( \widetilde{\Phi} := (\widetilde{\phi}_0, …, \widetilde{\phi}_d) \in \mathbb{R}^{d+1}\) , let

\[ \begin{equation} U_{\widetilde{\Phi}}= e^{\mathrm{i} \widetilde{\phi}_0 Z_{\Pi}} \prod_{j=1}^{d/2}\left[ U^{†}_A e^{\mathrm{i} \widetilde{\phi}_{2j-1} Z_{\Pi}}U_A e^{\mathrm{i} \widetilde{\phi}_{2j} Z_{\Pi}} \right] \end{equation} \]

(474)

when \( d\) is even, and

\[ \begin{equation} U_{\widetilde{\Phi}}=(-\mathrm{i})^d e^{\mathrm{i} \widetilde{\phi}_0 Z_{\Pi}} (U_A e^{\mathrm{i} \widetilde{\phi}_{1} Z_{\Pi}}) \prod_{j=1}^{(d-1)/2}\left[ U^{†}_A e^{\mathrm{i} \widetilde{\phi}_{2j} Z_{\Pi}}U_A e^{\mathrm{i} \widetilde{\phi}_{2j+1} Z_{\Pi}} \right] \end{equation} \]

(475)

when \( d\) is odd. Then

\[ \begin{equation} U_{\widetilde{\Phi}} = \begin{pmatrix} P(A) & *\\ * & * \end{pmatrix}\in \operatorname{BE}_{1,m}(P(A)), \end{equation} \]

(476)

where \( P\in\mathbb{C}\) satisfy the conditions in (3).

Following exactly the same procedure, we find that the circuit in (32), with \( U_{P(A)}\) given by (33) implements a \( U_{P_{\operatorname{Re}}(A)}\in\operatorname{BE}_{1,m+1}(P_{\operatorname{Re}}(A))\) . This is left as an exercise.

8.5.3 General matrix polynomials

In terms of implementing matrix polynomials of Hermitian matrices, quantum eigenvalue transform provides a much simpler circuit than the method based on LCU and qubitization (i.e., linear combination of Chebyshev polynomials). The simplification is clear both in terms of the number of ancilla qubits and of the circuit architecture. However, it is not clear so far for which polynomials (either a complex polynomial \( P\in\mathbb{C}\) or a real polynomial \( P_{\operatorname{Re}}\in\mathbb{R}\) ) we can apply the QET technique, and how to obtain the phase factors. Quantum signal processing (QSP) provides a complete answer to this question.

Due to qubitization, all these questions can be answered in the context of SU(2) matrices. QSP is the theory of QET for SU(2) matrices, or the unitary representation of a scalar (real or complex) polynomial \( P(x)\) . Let \( A=x\in[-1,1]\) be a scalar with a one-qubit Hermitian block encoding

Then

is a rotation matrix.

Similar to (467), the QSP representation takes the following form

By setting \( \phi_0=…=\phi_d=0\) , we immediately obtain the block encoding of the Chebyshev polynomial \( T_d(x)\) . The representation power of this formulation is characterized by (6), which is based on slight modification of [24, Theorem 4]. In the following discussion, even functions have parity \( 0\) and odd functions have parity \( 1\) .

8.6.1 QSP for real polynomials

8.6.2 Optimization based method for finding phase factors

8.6.3 A typical workflow preparing the circuit of QET

Using QET, let us now revisit the problem of time-independent Hamiltonian simulation problem \( \mathcal{U}=e^{\mathrm{i} H t}\) . Instead of Trotter splitting, we assume that we are given access to a block encoding \( U_H\in\operatorname{BE}_{\alpha,m}(H)\) . Since \( e^{\mathrm{i} H t}=e^{\mathrm{i} (H/\alpha) (\alpha t)}\) , the subnormalization factor \( \alpha\) can be factored into the simulation time \( t\) . So WLOG we assume \( U_H\in\operatorname{BE}_{1,m}(H)\) . Since

and that \( \cos(xt),\sin(xt)\) are even and odd functions, respectively, we can construct a block encoding for the real and imaginary part directly using the circuit for real matrix polynomials in (32).

More specifically, we first use the Fourier–Chebyshev series of the trigonometric functions given by the Jacobi-Anger expansion \( [-1,1]\) :

Here \( J_{\nu}(t)\) denotes Bessel functions of the first kind.

This series converges very rapidly. With

terms, the truncated Jacobi–Anger expansion with degree up to \( d=2r+1\) can approximate \( \cos (t x),\sin(tx)\) to precision \( \epsilon/\sqrt{2}\) , respectively [24, Corollary 32]. Such a scaling also matches the complexity lower bound for Hamiltonian simulation [27, 28]. Define

where \( \beta>1\) is chosen so that \( |C_d|,|S_d|\le 1\) on \( [-1,1]\) , and \( \beta\) can be chosen to be as small as \( 1+\epsilon\) . Also let \( f_d(x)=C_d(x)+\mathrm{i} S_d(x)\) , then

(7) guarantees the existence of phase factors \( \widetilde{\Phi}_C,\widetilde{\Phi}_S\) , using which we can construct \( \mathcal{U}_C\in\operatorname{BE}_{1,m+1}(C_d(H))\) and \( \mathcal{U}_s\in\operatorname{BE}_{1,m+1}(S_d(H))\) . Finally, we can use one more ancilla qubit and LCU in (8.5.3) to construct a block encoding \( \mathcal{U}_d\in\operatorname{BE}_{2,m+2}(f_d(H))\) , or \( \mathcal{U}_d\in\operatorname{BE}_{2\beta,m+2}(\mathcal{U},\epsilon)\) . The circuit depth is \( \mathcal{O}\left(t + \log(1/\epsilon)\right)\) .

As an example, (34) shows the QSP representation of the quality of approximating \( \cos(tx)\) using \( P_{\operatorname{Re}}(x)=C_d(x)\) with \( t=4\pi,\beta=1.001,d=24\) . The quality of the approximation can be significantly improved with a larger degree \( d=50\) (see (35)). The phase factors are obtained via QSPPACK.

Given a block encoding \( U_H\in\operatorname{BE}_{1,m}(H)\) , and WLOG assume \( 0\preceq H\preceq 1\) . We assume that we are provided an initial state \( \ket{\varphi}\) so that the initial overlap \( p_0=\left\lvert \braket{\varphi|\psi_0}\right\rvert ^2\) is not small. To simplify the problem we also assume that ground and the first excited state energies \( E_0,E_1\) are known with a positive gap \( \Delta:=E_1-E_0>0\) . Our goal is to use QET to prepare an approximate quantum state \( \ket{\psi}\approx \ket{\psi_0}\) .

To this end, we can construct a threshold polynomial approximating the sign function on \( [0,1]\) . Due to the assumption that \( 0\preceq H\preceq 1\) , this function can be chosen to be an even function. Since the sign function is a discontinuous, the polynomial should only aim at approximating the sign function outside \( (\mu-\Delta/2,\mu+\Delta/2)\) , where \( \mu=(E_0+E_1)/2\) . The polynomial also needs to satisfy the conditions in (7). We need to find an even polynomial \( P_{\operatorname{Re}}(x)\) satisfying

We can achieve this by approximating the sign function, with \( \deg(P_{\operatorname{Re}})=\mathcal{O}(\log(1/\epsilon) \Delta^{-1})\) . (see e.g. [29, Corollary 7] and [24, Corollary 16]). This construction is based on an approximation to the \( \mathrm{erf}\) function. (36) gives a concrete construction of the even polynomial obtained via numerical optimization, and the phase factors are obtained via QSPPACK.

Applying the circuit \( U_{\Phi}\) to the initial state \( \ket{\varphi}\) , we have

Here \( \ket{\psi_{\perp}}\) is a state in the system register orthogonal to \( \ket{\psi_0}\) , while \( \ket{\perp}\) is orthogonal to all states \( \ket{0^m}\ket{\psi}\) . Note that

Therefore if we measure the ancilla qubits, the success probability of obtaining \( 0^m\) in the ancilla qubits, and the ground state \( \ket{\psi_0}\) in the system register is at least \( p_0 (1-\epsilon)\) . So the total number of queries to to \( U_A\) and \( U_A^{†}\) is \( \mathcal{O}(\Delta^{-1}p_0^{-1}\log(1/\epsilon) )\) .

Using amplitude amplification, the number of repetitions can be reduced to \( \mathcal{O}(\gamma^{-1})\) , and the total number of queries to to \( U_A\) and \( U_A^{†}\) becomes \( \mathcal{O}(\Delta^{-1} p_0^{-\frac12}\log(1/\epsilon) )\) . This also matches the lower bound [30].

Once the ground state is prepared, we can estimate the ground state energy by measuring the expectation value \( \braket{\psi_0|H|\psi_0}\) . The number of samples needed is \( \mathcal{O}(1/\epsilon^2)\) , which can be reduced to \( \mathcal{O}(1/\epsilon)\) using amplitude amplification. In summary, the best complexity for estimating the ground state energy \( E_0\) to accuracy \( \epsilon\) is \( \mathcal{O}(\Delta^{-1}p_0^{-\frac12}\epsilon^{-1}\log(1/\epsilon) )\) . Note that the cost of estimating the ground state energy depends on the gap \( \Delta\) . This is because the algorithm first prepares the ground state and then estimates the ground state energy. If we are only interested in estimating \( E_0\) to precision \( \epsilon\) , the gap dependence is not necessary (see (5.1) as well as  [30]).

For a square matrix \( A \in \mathbb{C}^{N \times N}\) , where for simplicity we assume \( N=2^n\) for some positive integer \( n\) , the singular value decomposition (SVD) of the normalized matrix \( A\) can be written as

or equivalently

We may apply a function \( f(\cdot)\) on its singular values and define the generalized matrix function [31, 32] as below.

Given the form in (514), we also define two other types of generalized matrix functions.

Here the left and right pointing triangles reflects that the transformation only keeps the left and right singular vectors, respectively. For a given \( A\) , somewhat confusingly, in the discussion below, the transformation \( f^{\triangleright}(A),f^{\triangleleft}(A),f^{\diamond}(A)\) will all be referred to as singular value transformations. In particular, QSVT mainly concerns \( f^{\triangleright}(A),f^{\diamond}(A)\) .

WLOG we assume access to \( U_A\in\operatorname{BE}_{1,m}(A)\) , so that the singular values of \( A\) are in \( [0,1]\) , i.e.,

In (8.4) we have observed that when \( A\) is a Hermitian matrix, the qubitization procedure introduces two different subspaces \( \mathcal{H}_i\) and \( \mathcal{H}_i'\) associated with each eigenvector \( \ket{v_i}\) . In particular, \( U_A\) maps \( \mathcal{H}_i\) to \( \mathcal{H}_i'\) , and \( U_A^{†}\) maps \( \mathcal{H}_i'\) to \( \mathcal{H}_i\) . Furthermore, both \( \mathcal{H}_i\) and \( \mathcal{H}_i'\) are the invariant subspaces of the projection operator \( \Pi\) . Therefore \( \mathcal{H}_i\) is an invariant subspace of \( U_A^{†} f(\Pi) U_A\) for any function \( f\) . Much of the same structure can be carried out to the quantum singular value transformation. The only difference is that the qubitization is now defined with respect to the singular vectors. The procedure below almost entirely parallelizes that of (8.4), except that we need to work with the singular value decomposition instead of the eigenvalue decomposition.

Start from the SVD in (514), we apply \( U_A\) to \( \ket{0^m}\ket{v_i}\) and obtain

where \( \ket{\perp_i'}\) is a normalized state satisfying \( \Pi\ket{\perp_i'}=0\) .

Since \( U_A\) block encodes a matrix \( A\) , we have

which implies that there exists another normalized state \( \ket{\perp_i}\) satisfying \( \Pi\ket{\perp_i}=0\) and

Now apply \( U_A\) to both sides of (522), we obtain

which gives

Define

and the associated two-dimensional subspaces \( \mathcal{H}_i= \operatorname{span} B_i , \mathcal{H}'_i= \operatorname{span} B_i' \) , we find that \( U_A\) maps \( \mathcal{H}_i\) to \( \mathcal{H}_i'\) . Correspondingly \( U_A^{†}\) maps \( \mathcal{H}_i'\) to \( \mathcal{H}_i\) .

Then (520,524) give the matrix representation

Similar calculation shows that

Meanwhile both \( \mathcal{H}_i\) and \( \mathcal{H}_i'\) are the invariant subspaces of the projector \( \Pi\) , with matrix representation

Therefore

Hence \( \mathcal{H}_i\) is an invariant subspace of \( \widetilde{O}=U_A^{†}Z_{\Pi}U_A Z_{\Pi}\) , with matrix representation

Repeating \( k\) times, we have

In other words,

Therefore the circuit \( (U_A^{†}Z_{\Pi}U_A Z_{\Pi})^{k}\) gives \( (1,m)\) -block-encoding of \( T_{2k}^{\triangleright}(A)\) .

Similarly,

In other words,

Therefore the circuit \( U_A Z_{\Pi}(U_A^{†}Z_{\Pi}U_A Z_{\Pi})^{k}\) gives \( (1,m)\) -block-encoding of \( T_{2k+1}^{\diamond}(A)\) .

9.3.1 Quantum circuit

Due to the close relation between eigenvalue and singular value transformation in terms of Chebyshev polynomials and qubitization in (9.2), we can obtain the QSVT circuit easily following the discussion in (8.6).

First, there are no changes to the scalar case of QSP (in terms of SU(2) matrices), and in particular (6) and (7).

For the matrix case, when \( d\) is even,

\[ \begin{equation} U_{\Phi}=(-\mathrm{i})^d e^{\mathrm{i} \widetilde{\phi}_0 Z_{\Pi}} \prod_{j=1}^{d/2}\left[ U^{†}_A e^{\mathrm{i} \widetilde{\phi}_{2j-1} Z_{\Pi}}U_A e^{\mathrm{i} \widetilde{\phi}_{2j} Z_{\Pi}} \right] \end{equation} \]

(537)

gives a \( (1,m+1)\) -block-encoding of \( P^{\triangleright}(A)\) for some even polynomial \( P\in\mathbb{C}\) .

When \( d\) is odd,

\[ \begin{equation} U_{\Phi}=(-\mathrm{i})^d e^{\mathrm{i} \widetilde{\phi}_0 Z_{\Pi}} (U_A e^{\mathrm{i} \widetilde{\phi}_{1} Z_{\Pi}}) \prod_{j=1}^{(d-1)/2}\left[ U^{†}_A e^{\mathrm{i} \widetilde{\phi}_{2j} Z_{\Pi}}U_A e^{\mathrm{i} \widetilde{\phi}_{2j+1} Z_{\Pi}} \right] \end{equation} \]

(538)

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

The quantum circuit is exactly the same as that in (33). The phase factors can be adjusted so that all polynomials \( P\) satisfying the conditions in (6) can be exactly represented. If we are only interested in some real polynomial \( P_{\operatorname{Re}}\in\mathbb{R}\) and \( P_{\operatorname{Re}}^{\diamond}(A)\) (odd) and \( P^{\triangleright}(A)\) (even), we can use (7) and the circuit in (37) (which is simply a combination of (32,33)) to implement its \( (1,m+1)\) -block-encoding. We have the following theorem. Since the conditions of QSP representation for real polynomials is simple to satisfy and is also most useful in practice, we only state the case with real polynomials.

Theorem 8 (Quantum singular value transformation with real polynomials)

Let \( A\in\mathbb{C}^{N\times N}\) be encoded by its \( (1,m)\) -block-encoding \( U_A\) . Given a polynomial \( P_{\operatorname{Re}}(x)\in\mathbb{R}\) of degree \( d\) satisfying the conditions in (7), we can find a sequence of phase factors \( \Phi\in\mathbb{R}^{d+1}\) , so that the circuit in (37) denoted by \( U_{\Phi}\) implements a \( (1,m+1)\) -block-encoding of \( P^{\diamond}_{\operatorname{Re}}(A)\) if \( d\) is odd, and of \( P^{\triangleright}_{\operatorname{Re}}(A)\) if \( d\) is even. \( U_{\Phi}\) uses \( U_A,U_A^{†}\) , m-qubit controlled NOT, and single qubit rotation gates for \( \mathcal{O}(d)\) times.

9.3.2 QSVT applied to Hermitian matrices

9.3.3 QSVT and matrix dilation

In this section, we revisit the problem of solving linear systems of equations \( Ax=b\) . With QSVT we can solve QLSP for general matrices without the need of dilating the matrix into a Hermitian matrix. Assume \( A=W\Sigma V^{†}\) is invertible, i.e., \( \forall i, \Sigma_{ii}>0\) , then

where \( f(x)=x^{-1}\) is an odd function.

WLOG assume \( A\) can be accessed by \( U_A\in\operatorname{BE}_{1,m}(A)\) . For simplicity we also assume \( \left\lVert A\right\rVert =1\) (though in general \( \left\lVert A\right\rVert \le \alpha=1\) , and we may not always be able to set \( \left\lVert A\right\rVert =1\) ). Let \( \kappa\) be the condition number of \( A\) , then the singular values of \( A\) are contained in the interval \( [\delta,1]\) , with \( \delta =\kappa^{-1}\) .

Note that \( f(\cdot)\) is not bounded by 1 and in particular singular at \( x=0\) . Therefore instead of approximating \( f\) on the whole interval \( [-1,1]\) we consider an odd polynomial \( p(x)\) such that

The \( \beta\) factor is chosen arbitrarily so that \( |p(x)|\leq 1\) for all \( x\in[-1,1]\) to satisfy the requirement of the condition (2) in (7). For instance, we may choose \( \beta=4/3\) . The precision parameter \( \epsilon'\) will be chosen later. The degree of the odd polynomial can be chosen to be \( \mathcal{O}(\frac{1}{\delta}\log(\frac{1}{\epsilon'}))\) is guaranteed by e.g. [21, Corollary 69]. This construction is not explicit (see an explicit construction in [22]). (38) gives a concrete construction of an odd polynomial obtained via numerical optimization, and the phase factors are obtained via QSPPACK.

Then (37) implements a \( (1,m+1)\) -block-encoding of \( p^{\diamond}(A^\dagger)=V p(\Sigma)W^\dagger\) denoted by \( U_{\Phi}\) . We have

The total number of queries to to \( U_A\) and \( U_A^{†}\) is

To solve QLSP, we assume the right hand side vector \( \ket{b}\) can be accessed through the oracle \( U_b\) such that

We introduce the parameter

which plays an important part in the success probability of the procedure. Let

be the unnormalized true solution, and the normalized solution state is \( \ket{x}=x/\left\lVert x\right\rVert \) . Now denote \( \widetilde{x}=p^{\diamond}(A^{\dagger})\ket{b}\) , and \( \ket{\widetilde{x}}=\widetilde{x}/\left\lVert \widetilde{x}\right\rVert \) . Then the unnormalized solution satisfies \( \|x-\widetilde{x}\|\leq\epsilon'\) . For the normalized state \( \ket{y}\) , this error is scaled accordingly. When \( \epsilon'\ll \|\ket{\widetilde{x}}\|\) , we have

Also we have

Therefore in order for the normalized output quantum state to be \( \epsilon\) -close to the normalized solution state \( \ket{x}\) , we need to set \( \epsilon'=\mathcal{O}(\epsilon \xi /\kappa)\) . This is similar to the case of the HHL algorithm in (5.3), where QPE needs to achieve \( \epsilon\) multiplicative accuracy, which means that the additive accuracy parameter \( \epsilon'\) should be set to \( \mathcal{O}(\epsilon/\kappa)\) .

The success probability of the above procedure is \( \Omega(\|\widetilde{x}\|^2)=\Omega(\xi^2/\kappa^2)\) . With amplitude amplification we can boost the success probability to be greater than \( 1/2\) with one qubit serving as a witness, i.e., if measuring this qubit we get an outcome 0 it means the procedure has succeeded, and if 1 it means the procedure has failed. It takes \( \mathcal{O}(\kappa/\xi)\) rounds of amplitude amplification, i.e., using \( U_{\Phi}^\dagger\) , \( U_{\Phi}\) , \( U_b\) , and \( U_b^\dagger\) for \( \mathcal{O}(\kappa/\xi)\) times. A single \( U_{\Phi}\) uses \( U_A\) and its inverse

times. Therefore the total number of queries to \( U_A\) and its inverse is

The number of queries to the \( U_b\) and its inverse is \( \mathcal{O}(\kappa/\xi)\) . We consider the following two cases for the magnitude of \( \xi\) .

1. In general if no further promise is given, then \( \xi\geq 1\) . The total query complexity of \( U_A\) is therefore \( \mathcal{O}(\kappa^2 \log(\kappa/\epsilon))\) . This is the typical complexity referred to in the literature.
2. If \( \ket{b}\) has a \( \Omega(1)\) overlap with the left-singular vector of \( A\) with the smallest singular value, then \( \xi=\Omega(\kappa)\) . This is the best case scenario, and the total query complexity of \( U_A\) is \( \mathcal{O}(\kappa \log(1/\epsilon))\) , and the number of queries to the right hand side vector \( \ket{b}\) is \( \mathcal{O}(1)\) , which is independent of the condition number.

So far we have assumed that we have a matrix \( A\) in mind, and the access to \( A\) is provided via the block encoding matrix \( U_A\) . The QSVT circuit takes the form

If we are further given two unitary matrices \( P,Q\) , we can equivalently rewrite the QSVT transformation as

The beginning of the equation ends with \( P\) or \( Q\) depending on the parity of \( d\) . The insertion of \( P,Q\) amounts to a basis transformation. Assume that we have access to \( \widetilde{U}_A=QU_A P^{†}\) , then \( U_A\) is the matrix representation of \( \widetilde{U}_A\) with respect to the bases given by \( P,Q\) , respectively. What is different is that the controlled rotations before \( U_A\) and \( U_A^{†}\) are now expressed with respect to two different basis sets, i.e., \( P \operatorname{CR} _{\widetilde{\phi}_d}P^{†}\) , and \( Q \operatorname{CR} _{\widetilde{\phi}_{d-1}}Q^{†}\) , respectively. This can be useful for certain applications. Let us now express these ideas more formally.

Assume that we are given an \( (n+m)\) -qubit unitary \( \widetilde{U}_A\) , and two \( (n+m)\) -qubit projectors \( \Pi,\Pi'\) . For simplicity we assume \( \operatorname{rank} (\Pi)= \operatorname{rank} (\Pi')=N\) . Define an orthonormal basis set

where the vectors \( \ket{\varphi_0},…,\ket{\varphi_{N-1}}\) span the range of \( \Pi\) , and all states \( \ket{v_i}\) are orthogonal to \( \ket{\varphi_j}\) . Similarly define an orthonormal basis set

where the vectors \( \ket{\psi_0},…,\ket{\psi_{N-1}}\) span the range of \( \Pi'\) , and all states \( \ket{w_i}\) are orthogonal to \( \ket{\psi_j}\) . We can think that the columns of \( \mathcal{B},\mathcal{B}'\) form the basis transformation matrix \( P,Q\) , respectively.

Then the matrix \( A\) is defined implicitly in terms of its matrix representation

Note that

we find that

Therefore (8) can be viewed as the singular value transformation of \( A\) , which is a submatrix of the matrix representation of \( U_A\) with respect to bases \( \mathcal{B},\mathcal{B}'\) .

The implementation of the controlled rotation \( P \operatorname{CR} _{\widetilde{\phi}}P^{†}\) relies on the implementation of \( \operatorname{C} _{\Pi} \operatorname{NOT} \) . The projectors \( \Pi,\Pi'\) can be accessed directly, and WLOG we focus on one projector \( \Pi\) . Motivated from Grover's search, we may assume access to a reflection operator

via the controlled NOT gates \( \operatorname{C} _{\Pi} \operatorname{NOT} , \operatorname{C} _{\Pi'} \operatorname{NOT} \) respectively. We can then define an \( m\) -qubit controlled NOT gate as

which can be constructed using \( R_{\Pi}\) as

Therefore assuming access to \( R_{\Pi}\) , the \( \operatorname{C} _{\Pi} \operatorname{NOT} \) gate can be implemented using the circuit in (39).

Then according to (8), we can implement the QSVT using \( \widetilde{U}_A,\widetilde{U}_A^{†}\) , \( \operatorname{C} _{\ket{0^m}\bra{0^m}} \operatorname{NOT} \) and single qubit rotation gates, where \( \ket{0^m}\bra{0^m}\otimes I_n\) is a rank-\( 2^n\) projector. In this generalized setting, \( \operatorname{C} _{\ket{0^m}\bra{0^m}\otimes I_n} \operatorname{NOT} = \operatorname{C} _{\ket{0^m}\bra{0^m}} \operatorname{NOT} \otimes I_n\) in the \( \mathcal{B},\mathcal{B}'\) basis should be implemented using \( \operatorname{C} _{\Pi} \operatorname{NOT} , \operatorname{C} _{\Pi'} \operatorname{NOT} \) , respectively. We arrive at the following theorem:

As an application, let us revisit the Grover search problem from the perspective of QSVT. Again let \( \ket{x_0}\) be the desired marked state, and \( \ket{\psi_0}\) be the uniform superposition of states. Now let us perform a basis transformation. We define an orthonormal basis set \( \mathcal{B}=\{\ket{\psi_0},\ket{v_1},…,\ket{v_{N-1}}\}\) , where all states \( \ket{v_i}\) are orthogonal to \( \ket{\psi_0}\) . Similarly define an orthonormal basis set \( \mathcal{B}'=\{\ket{x_0},\ket{w_1},…,\ket{w_{N-1}}\}\) , where all states \( \ket{w_i}\) are orthogonal to \( \ket{x_0}\) . Then the matrix of reflection operator \( R_{\psi_0}\) with respect to \( \mathcal{B},\mathcal{B'}\) is (let \( a=1/\sqrt{N}\) )

Let \( A=a\) be a \( 1\times 1\) matrix, and then \( R_{\psi_0}\in \operatorname{BE}_{1,n}(A)\) . Furthermore, the projectors \( \Pi=\ket{\psi_0}\bra{\psi_0}\) and \( \Pi'=\ket{x_0}\bra{x_0}\) can be accessed via the reflection operator \( R_{\psi_0},R_{x_0}\) , respectively. According to (570), this defines \( \operatorname{C} _{\Pi} \operatorname{NOT} , \operatorname{C} _{\Pi'} \operatorname{NOT} \) . Let \( \widetilde{U}_A=R_{\psi_0}\) , we have

and we would like to use (9) to find \( U_{\Phi}\) that block encodes

To this end, we need to find an odd, real polynomial \( P_{\operatorname{Re}}(x)\) satisfying \( P_{\operatorname{Re}}(a)\approx 1\) . More specifically, we need to find \( P_{\operatorname{Re}}(x)\) satisfying

We can achieve this by approximating the sign function, with \( \deg(P_{\operatorname{Re}})=\mathcal{O}(\log(1/\epsilon^2) a^{-1})=\mathcal{O}(\log(1/\epsilon)\sqrt{N})\) (see e.g. [29, Corollary 6]). This construction is also based on an approximation to the \( \mathrm{erf}\) function. (40) gives a concrete construction of an odd polynomial obtained via numerical optimization, and the phase factors are obtained via QSPPACK.

Then

More specifically, note that

for some unnormalized state \( \ket{\widetilde{\perp}}\) . Moreover

which gives

So

Therefore we can measure the system register to find \( x_0\) , and we achieve the same Grover type speedup.

Note that the approximation can be arbitrarily accurate, and there is no overshooting problem (though this is a small problem) as in the standard Grover search. While Grover’s search does not require the output quantum state to be exactly \( \ket{x_0}\) , this could be desirable when it is used as a quantum subroutine, such as amplitude amplification.

The immediate generalization of the procedure above is called the fixed point amplitude amplification.

Note that the ranks of \( \Pi',\Pi\) are different. This does not affect the proof.

[1] Michael A Nielsen and Isaac Chuang Quantum computation and quantum information 2000

[2] Alexei Yu Kitaev and Alexander Shen and Mikhail N Vyalyi Classical and quantum computation American Mathematical Soc. 2002 47

[3] Mikio Nakahara Quantum computing: from linear algebra to physical realizations CRC press 2008

[4] Eleanor G Rieffel and Wolfgang H Polak Quantum computing: A gentle introduction MIT Press 2011

[5] Scott Aaronson Quantum computing since Democritus Cambridge Univ. Pr. 2013

[6] John Preskill Lecture notes for Physics 219: Quantum computation Caltech Lecture Notes 1999

[7] Ronald De Wolf Quantum computing: Lecture notes arXiv:1907.09415 2019

[8] Andrew Childs Lecture Notes on Quantum Algorithms 2021

[9] Daniel Nagaj and Pawel Wocjan and Yong Zhang Fast amplification of QMA Quantum Inf. Comput. 2009 9 11 1053–1068

[10] Lin Lin and Yu Tong Optimal quantum eigenstate filtering with application to solving quantum linear systems Quantum 2020 4 361

[11] J. W. Negele and H. Orland Quantum many-particle systems Westview 1988

[12] Aram W Harrow and Avinatan Hassidim and Seth Lloyd Quantum algorithm for linear systems of equations Phys. Rev. Lett. 2009 103 150502

[13] Danial Dervovic and Mark Herbster and Peter Mountney and Simone Severini and Naïri Usher and Leonard Wossnig Quantum linear systems algorithms: a primer arXiv:1802.08227 2018

[14] G. H. Golub and C. F. Van Loan Matrix computations Johns Hopkins Univ. Press 2013 Baltimore Fourth

[15] Richard P Feynman Simulating physics with computers Int. J. Theor. Phys 1982 21 6/7

[16] Seth Lloyd Universal quantum simulators Science 1996 1073–1078

[17] Tobias Jahnke and Christian Lubich Error bounds for exponential operator splittings BIT 2000 40 4 735–744 10.1023/A:1022396519656

[18] Mechthild Thalhammer High-order exponential operator splitting methods for time-dependent Schrödinger equations SIAM J. Numer. Anal. 2008 46 4 2022–2038

[19] Andrew M Childs and Yuan Su and Minh C Tran and Nathan Wiebe and Shuchen Zhu Theory of trotter error with commutator scaling Phys. Rev. X 2021 11 1 011020

[20] Yulong Dong and Lin Lin Random circuit block-encoded matrix and a proposal of quantum LINPACK benchmark Phys. Rev. A 2021 103 6 062412

[21] András Gilyén and Yuan Su and Guang Hao Low and Nathan Wiebe Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics arXiv:1806.01838 2018

[22] Andrew M. Childs and Robin Kothari and Rolando D. Somma Quantum Algorithm for Systems of Linear Equations with Exponentially Improved Dependence on Precision SIAM J. Comput. 2017 46 1920–1950

[23] Mario Szegedy Quantum speed-up of Markov chain based algorithms 45th Annual IEEE symposium on foundations of computer science 2004 32–41

[24] András Gilyén and Yuan Su and Guang Hao Low and Nathan Wiebe Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing 2019 193–204

[25] Yulong Dong and Xiang Meng and K Birgitta Whaley and Lin Lin Efficient phase factor evaluation in quantum signal processing Phys. Rev. A 2021 103 042419

[26] J. Haah Product decomposition of periodic functions in quantum signal processing Quantum 2019 3 190

[27] Dominic W Berry and Graeme Ahokas and Richard Cleve and Barry C Sanders Efficient quantum algorithms for simulating sparse Hamiltonians Commun. Math. Phys. 2007 270 2 359–371

[28] Guang Hao Low and Isaac L. Chuang Optimal Hamiltonian Simulation by Quantum Signal Processing Phys. Rev. Lett. 2017 118 010501

[29] Guang Hao Low and Isaac L Chuang Hamiltonian simulation by uniform spectral amplification arXiv:1707.05391 2017

[30] Lin Lin and Yu Tong Near-optimal ground state preparation Quantum 2020 4 372

[31] J. B. Hawkins and A. Ben-Israel On generalized matrix functions Linear and Multilinear Algebra 1973 1 163–171

[32] F. Arrigo and M. Benzi and C. Fenu Computation of Generalized Matrix Functions SIAM J. Matrix Anal. Appl. 2016 37 836–860