arXiv:2108.09004v4 [quant-ph] 25 Mar 2023Date of publication xxxx 00, 0000, date of current version xxxx 00, 0000.
                                         Digital Object Identifier 10.1109/ACCESS.2022.0122113

                                         A Step-by-Step HHL Algorithm Walkthrough to
                                         Enhance Understanding of Critical Quantum
                                         Computing Concepts

                                         ANIKA ZAMAN1, HECTOR JOSE MORRELL2, HIU YUNG WONG.3, (Senior Member, IEEE)

                                         1Department of Electrical Engineering, San Jose State University, San Jose, CA 95192 USA (email: anika.zaman@sjsu.edu)
                                         2Department of Electrical Engineering, San Jose State University, San Jose, CA 95192 USA (email:
                                         hector.morrell@sjsu.edu)
                                         3Department of Electrical Engineering, San Jose State University, San Jose, CA 95192 USA (email:
                                         hiuyung.wong@sjsu.edu)
                                         Corresponding author: Hiu Yung Wong (email: hiuyung.wong@sjsu.edu). A. Zaman and H. J. Morrell have equal
                                         contributions.

                                         Part of this work was supported by the National Science Foundation under Grant 2046220.

                                           ABSTRACT After learning basic quantum computing concepts, it is desirable to reinforce the
                                           learning using an important and relatively complex algorithm through which students can observe
                                           and appreciate how qubits evolve and interact with each other. Harrow-Hassidim-Lloyd (HHL)
                                           quantum algorithm, which can solve linear system problems with exponential speed-up over the
                                           classical method and is the basis of many important quantum computing algorithms, is used to
                                           serve this purpose. The HHL algorithm is explained analytically followed by a 4-qubit numerical
                                           example in bra-ket notation. Matlab code corresponding to the numerical example is available for
                                           students to gain a deeper understanding of the HHL algorithm from a pure matrix point of view.
                                           A quantum circuit programmed using qiskit is also provided for real hardware execution in IBM
                                           quantum computers. After going through the material, students are expected to have a better
                                           appreciation of the concepts such as basis transformation, bra-ket and matrix representations,
                                           superposition, entanglement, controlled operations, measurement, quantum Fourier transformation,
                                           quantum phase estimation, and quantum programming. To help readers review these basic concepts,
                                           brief explanations augmented by the HHL numerical examples in the main text are provided in the
                                           Appendix.

                                           INDEX TERMS Harrow-Hassidim-Lloyd (HHL) quantum algorithm, Quantum Fourier transform
                                           (QFT), inverse quantum Fourier transform (IQFT), quantum phase estimation(QPE),linear system
                                           problem (LSP), Quantum Education

I. INTRODUCTION                                             technology is based on concepts very different from our
                                                            daily and classical experiences. In the early stage of
Q UANTUM Computing is promising in solving chal-            learning quantum computing, although linking to daily
       lenging engineering [1], biomedical [2] and finance  and classical experience may enhance the understanding
[3] problems. It has a tremendous advancement in the        of certain quantum concepts and such an approach
last two decades and, recently, quantum breakthrough        should not be de-emphasized, we believe a fast and
has been demonstrated using a 53-qubit system [4].How-      robust way of training a quantum workforce is to train
ever, according to the paper [5], due to less efficient     the students to be able to emulate a quantum processor
hardware implementation till date, the goal to reach        and trace the evolution of the qubits. This is partic-
superconducting quantum supremacy is yet to achieve.        ularly useful in learning quantum algorithms without
Therefore, the training of a quantum technology work-       a quantum mechanics background. Such an approach
force is an imminent goal for many countries (e.g. [23])    obviates the students from cognitive conflicts, which can
to support this fast-growing industry. However, quantum
                                                                                                                                                                   1
VOLUME 3, 2022
                                                                                                            Author et al.:

be resolved later, if possible, after they understand how   as in this paper can be found in [19].
quantum computing works. This also embraces the "Shut
up and calculate!" approach proposed by Mermin on how       A. HOW TO USE THIS PAPER
to deal with the uncomfortable feeling toward quantum       For readers who have a fresh memory of the basic con-
mechanics interpretation [6].                               cepts, they can start reading from Section II, in which
                                                            the HHL algorithm is discussed step-by-step analytically
   Besides analytical equations, matrix representation      followed by a numerical example in Section III. The basic
and computer simulations are important tools to en-         concepts mentioned in the Appendix are referred to in
hance the understanding of qubit evolution. However,        the main text and readers are encouraged to review them
available examples that include computer simulations        when needed.
are usually of simple algorithms and, very often, with-
out matrix representation. There is a lack of examples         For readers who need reviews on the basic concepts
of important and relatively complex algorithms which        first, they are encouraged to go over the Appendix first
combine some of the most important quantum comput-          before reading the main text.
ing concepts and basic algorithms. Such examples are
desirable to allow students to appreciate the roles and        For readers who have devoted substantial time to
the interplay of various basic concepts in a more realis-   learning HHL elsewhere but just need a numerical ex-
tic quantum algorithm. Harrow-Hassidim-Lloyd (HHL)          ample to reinforce the understanding, they might start
quantum algorithm [7] [8] which can be used to solve lin-   with the numerical example in Section III.
ear system problems (LSP) and can provide exponential
speedup over the classical conjugate gradient method           Equations in the Appendix begin with 'V'. If the
[9] is chosen for this purpose. HHL is the basic of many    equations are examples from the main text, the same
more advanced algorithms and is important in various        equation number is used in the Appendix.
applications such as machine learning [10] and modeling
of quantum systems [2] [11]. HHL solves system of           II. HHL ALGORITHM
linear equation which is a discretization of [12] [13].     A. DEFINITIONS AND OVERVIEW
In this paper, we detail the qubit evolution in Harrow-     We will first give an overview of the problem and the
Hassidim-Lloyd (HHL) quantum algorithm analytically         HHL algorithm. Details will be discussed in the following
with a 4-qubit circuit as a numerical example. Although     subsections with reference to the Appendix for reviewing
HHL examples are available elsewhere (e.g. [14] [15]),      basic concepts. A linear system problem (LSP) can be
this paper has certain characteristics which are not all    represented as the following
found in those examples. Firstly, the HHL algorithm is
discussed analytically step-by-step and is self-contained.                      Ax = b                          (1)
Secondly, a numerical example is given in bra-ket no-
tation mirroring the analytical equations. Thirdly, a       where A is a Nb  Nb Hermitian matrix and x and b
Matlab code corresponding exactly to the numerical          are Nb-dimensional vectors. For simplicity, it is assumed
example is available to enhance the understanding from      Nb = 2nb , where nb is the number of qubits in the
a matrix point of view. The Matlab code allows the          quantum circuit, and Nb is the total combinations due
students to trace how the wavefunction evolves instead      number of qubits,nb. In matrix representation, qubits
of just seeing the magnitudes of the coefficients as in     are represented by their total combinations(Nb  Nb),
IBM-Q. Fourthly, a qiskit code written in python [16]       or we can say for Nb number of unknowns, we need
corresponding to the numerical example is available and     nb qubits to solve unknowns. Dummy equations can
can be run in IBM simulation and hardware machines          be added otherwise to convert the system satisfy this
[17]. Finally, in the example, all the 4 qubits are traced
throughout the process without simplification.              assumption. A and b are known and x is the unknown

   The readers are assumed to have the following back-      to be solved, i.e.
ground concepts which are further enhanced through
the step-by-step walkthrough of the HHL algorithm: ba-                          x = A-1b                        (2)
sis transformation, bra-ket and matrix representations,
superposition, entanglement, encoding, controlled oper-     As an example, A =  1     -   1  ,b=      0     , and x =
ations, measurement, quantum Fourier transformation,                                      3           1
and quantum programming. To make this paper self-                                  1
contained and to help the readers better appreciate the                         -  3    1
roles of these basic concepts in the HHL, an Appendix
is devoted to briefly explaining these concepts using       3
the examples from the main text. A more detailed
explanation of these concepts using a similar approach      89 with nb = 1 and Nb = 21 = 2. Readers may refer

2                                                              8

                                                            to Appendix V-13 and Appendix V-14 to review how

                                                            LSP is solved classically using Gaussian Elimination and

                                                            Conjugate Gradient Method, respectively.

                                                            A is assumed to be Hermitian (See Appendix V-1). If

                                                            it is not Hermitian, the A can be converted to   0  A  ,
                                                                                                            A   0

                                                            which is Hermitian. Readers may refer to [20] for the

                                                            more advanced treatment when A is not Hermitian.

                                                                                                            VOLUME 3, 2022
Author et al.:

FIGURE 1. Schematic of the HHL quantum circuit flowing from left to right. The circuit is decomposed into top and bottom portions for clarity. Note that
the lowest qubit in the diagram is the most significant bit (MSB) while the top one is the least significant bit (LSB).

   Figure 1 shows the schematic of the HHL algorithm        convention used in qiskit [16] and the IBM-Q platform
and the corresponding circuit to solve LSP. In the HHL      [17]
quantum algorithm, the Nb components of b and x
are encoded as the amplitudes/coefficients (amplitude          As shown in Figure 1, besides the b-register, which
encoding) of basis states of the nb-qubits, | b, which      belongs to the more significant bits, there are two
form a CNb Hilbert space. These nb qubits are called        more sets of inputs to the algorithm. The first set is
b-register. Qubit nb is chosen to be large enough to        sometimes called the c-register because it is related
encode b, i.e. 2nb needs to be the same as the length of    to the time (clock) in the controlled rotation in the
the vectors b and x. The matrix A is simulated through      QPE part. Therefore, they are also called the clock
Hamiltonian encoding by encoding it as the Hamiltonian      qubits.The c-register stores the values of the phase of
of a unitary gate. Appendix V-9 reviews examples of         the eigenvalues of the A matrix after the QPE. There are
various encoding schemes.                                   n qubits in the c-register. Since basis encoding is used
                                                            (i.e. the phase value is encoded as the basis number (See
   The HHL algorithm has 5 main components, namely          Appendix V-9), the value of n determines how accurately
state preparation, quantum phase estimation (QPE),          the phase can be stored. A larger n results in higher
ancilla bit rotation, inverse quantum phase estimation      accuracy when the encoding is not exact. We set N = 2n.
(IQPE), and measurement. In this paper, the little-
endian convention is used. In a little-endian convention,      The last set of qubits is the ancilla qubit | a which
the rightmost (ending) qubit represents the least signifi-  is the LSB. The ancilla qubit, as its name implies, is
cant bit (LSB). For example, in a 4-qubit system, |0001     important to help achieve the goal although it will be
in binary is |1 in decimal because the 1 in the basis       discarded at the end, as will be detailed later.
state |0001 is the LSB, representing 20 instead of 23 (if
it were the most significant bit, MSB). Moreover, in the       The matrix A, which is a Hamiltonian, may be written
circuit diagrams, the lowest qubit represents the MSB       as a linear combination of the outer products of its
and the topmost qubit represents the LSB, which is a        eigenvectors, |ui ui| weighted by its eigenvalues, i,in
                                                            Eq.(3), (See Appendix V-8).
VOLUME 3, 2022
                                                                                                                                                                   3
                                                                                                              Author et al.:

                  2nb -1                                         0 

           A=               i |ui ui|          (3)               1 

                       i=0                                      b=     ...      0  |0  +1  |1    +  +Nb-1  |Nb  -  1  = |b
                                                                            
                                                                     
Since A is diagonal in its eigenvector basis, its inverse                   

is simply, A-1 =  2nb -1    i-1   |ui     ui|. b can be also    Nb-1
                  i=0                                                                                                        (9)

expressed in the basis formed by the eigenvectors of A,         The vector b is represented in a column form on the

such that                 2nb -1                                left with coefficients  s. This is also a valid represen-

                  |b =            bj |uj       (4)              tation of |b . On the right, the corresponding basis of

                          j=0                                   the Hilbert space formed by the nb qubits is written
                                                                explicitly. Therefore,
Therefore, Eq. (2) can be encoded as,

                            2nb -1                                              |1 = |b b |0    0 c |0 a              (10)

           |x = A-1 |b =            i-1bi |ui  (5)                 From now on, some of the subscripts of the kets will
                                                                be omitted when there is no ambiguity. Since the state
                               i=0                              preparation depends on the actual value of b, it will be
                                                                discussed in more detail in the numerical example.
by using the fact that ui|uj = ij. The goal of the HHL
algorithm is to find the solution in this form and |x is        C. QUANTUM PHASE ESTIMATION
stored in the b-register.Storing the right hand side of
(5) in the b-register is equivalent to storing |x in the b-     Quantum phase estimation (QPE) is also an eigen-
register. But the solution is encoded as the amplitudes         value estimation algorithm. QPE has three components,
of the basis vectors |0 / |1 (the measurement basis).           namely the superposition of the clock qubits through
Therefore, the solutions are not -i 1bi, which are the          Hadamard gates, controlled rotation, and inverse quan-
amplitudes of the eigenvector basis vectors. However,           tum Fourier transform (IQFT). The goal of QPE is to
one will naturally get the correct amplitudes if it is          estimate the phase of the eigenvalues of the unitary
measured in the |0 / |1 basis and this is only possible         rotation matrix, U = eiAt, in the controlled gate, C - U,
if the qubits are not entangled with other qubits. This         (Fig. 1) used in the QPE. Again, this gate encodes the
can be checked mathematically replacing |ui by |0 / |1          matrix A as its Hamiltonian. It is also instructive to
based on their relationship. If |ui are entangled with          note that the eigenvalues of U must be roots of unity
other qubits, one cannot obtain the desired solution.           (i.e. in the form of ei) as U is unitary. Therefore, the
This will be clear after the ancilla bit rotation to be         phase of the eigenvalue of the gate is proportional to the
detailed later.                                                 eigenvalue of A. As a result, by using QPE in the HHL
                                                                algorithm, it is expected the eigenvalues of A will be
   The equation needs also to be prepared so that the           encoded in the c-register after the QPE, i.e. at |4 . As
eigenvectors, |ui , and |b are normalized so they can           it will be clear later, the eigenvalues are only encoded
be properly represented as a unit vector in quantum             through basis encoding. The c-register does not store the
computing. Therefore, (4) and Eq. (5) require                   exact eigenvalues.

                  2nb -1                                           Here we assume the readers are already familiar with
                                                                IQFT and it will not be explained in detail. Readers may
                        |bj|2 = 1                               review the basic concepts in Appendices V-11 and V-12.

                    j=0                                            In the first step of QPE, Hadamard gates are applied
                                                                to the clock qubits to create a superposition of the clock
                                                           (6)  qubits,

                  2nb -1

                          |i-1bi|2 = 1         (7)

                  i=0

                                                                       |2 =                Inb  Hn  I |1                 (11)
                                                                             =
B. STATE PREPARATION                                                                   |b  1     (|0  + |1 )n |0         (12)

There are total nb + n + 1 qubits and they are initialized                                 2  n
as                                                                                            2

    |0 = |0    0 b |0    0 c |0 a = |0 nb |0 n |0 (8)        In the controlled rotation part, controlled gates are
                                                                applied to |b with the clock qubits as the controlling
   In the state preparation, |0    0 b in the b-register     qubits (Figure 2). The number of the qubit, n, of the c-
needs to be rotated to have the amplitudes correspond           register determines the number of the controlled gates.
to the coefficients of b. That is                               The gates are in the form of U2r , where r is the index of
                                                                the clock qubit. Also, U = eiAt. For the most significant
4                                                               bit in the c-register, |cn-1 , it controls the gate U2n-1

                                                                                                              VOLUME 3, 2022
Author et al.:

                                                                                       due to destructive interference. By ignoring the states
                                                                                       with zero amplitude, we may rewrite |4 as

                                                                                             |4      =  1       2n -1            |0  a
                                                                                                        2n |b
                                                                                                                      e2ik0 |N

                                                                                                                k=0

                                                                                                     = |b |N |0 a                       (16)

FIGURE 2. The controlled-rotation part of QPE. U is replaced by eiAt in the               Therefore, in QPE, the clock qubits are used to
HHL algorithm.                                                                         represent the phase information of U, which is , and
                                                                                       the accuracy depends on the number of qubits, n.

                                                                                          Since in Hamiltonian encoding, U is related to A
                                                                                       through

on the b-register while the least significant one, |c0 ,                                                       U = eiAt                 (17)
controls the gate U20 = U on the b-register.
                                                                                       where t is the evolution time for that Hamiltonian. U is
   To understand how it works, we begin by assuming                                    also diagonal in A s eigenvector, |ui , basis. If |b = |uj ,
that |b is an eigenvector of U with eigenvalue e2i.
The eigenvalue is written in this form so that the phase,                                               U |b = eijt |uj                 (18)
, will be encoded as the basis state in the c-register
(See Eq. (16)). Therefore, based on the definition of                                     By equating ijt to 2i in Eq. (13), we get  =
eigenvalues and eigenvectors (See Section V-8),                                        jt/2 and Eq. (16) becomes

                            U |b = e2i |b                               (13)

   When the control clock qubit is |0 , |b will not be                                               |4 = |uj |Njt/2 |0 a               (19)

affected. If the clock bit is |1 , U will be applied to |b .                              Thus the eigenvalues of A have been encodeded in the
This is equivalent to multiplying e2i2j in front of the                                clock qubits (basis encoding). However, in general, given
                                                                                       in (4), by superposition,
|1 of the jth clock qubit, |cj , as one can assign the
prefactor to the controlling qubit. Therefore, after the                                                2nb -1

controlled-U operation, we have

                                                                                             |4 =               bj |uj |Njt/2 |0 a      (20)

                      1             + e2i2n-1 |1 )  (|0                                                 j=0

|3  = |b              2n    (|0                              +                            The j are usually not integers. We will choose t so
                         2                                                             that ~j = Njt/2 are integers. Therefore, ~j are usually
                                                                                       scaled version of j.
    e2i2n-2 |1 )      (|0 + e2i20 |1 )  |0 a
                                                                                          4 can be rewritten as
    = |b         1    2n -1              |0 a                           (14)
                2n
                            e2ik |k
                   2
                      k=0
                                                                                                           2nb -1

   In the IQFT part,(15), only the clock qubits are                                                  |4 =          bj |uj |~j |0 a      (21)
affected. Note that in certain literature, this is called
Quantum Fourier Transform (QFT) (Appendix V-11).                                                             j=0

|4  = |b        IQFT(            1  2n -1          ) |0                                D. CONTROLLED ROTATION AND MEASUREMENT OF THE
                                                                                       ANCILLA QUBIT
                            2    n        e2ik |k        a
                                 2                                                     The next step is to rotate the ancilla qubit, |0 a, based
                                                                                       on the encoded eigenvalues in the c-register, such that,
                                    k=0

    = |b         1    2n -1                        ) |0  a
                2n
                            e2ik(IQFT |k
                   2
                         k=0                                                                 2nb -1
                                                                                                                         C2          C
                      2n -1              2n -1                                         |5 =          bj |uj |~j (  1 - ~j2 |0 a + ~j |1 a) (22)
                1
    = |b        2n               e2ik(          e-2iyk/N |y    ) |0  a                       j=0

                      k=0                y=0                                           where C is a constant. The goal is to show why this is
                                                                                       useful.
        1             2n-1 2n-1                          |0 a           (15)
    = 2n |b                                                                               When the ancilla qubit is measured, the ancilla qubit
                                  e2ik(-y/N) |y                                        wavefunction will collapse to either |0 or |1 . If it is |0 ,
                                                                                       the result will be discarded and the computation will
                      y=0 k=0                                                          be repeated until the measurement is |1 . Therefore, the
                                                                                       final wavefunction of interest is
Due to interference, only |y satisfying the condition
                                         ias mp2klin=t-u01dee2oifk(-2kyn=/-0N1)e=0 =0
 - y/N = 0 will have a finite
2n. Otherwise, the amplitude

VOLUME 3, 2022                                                                                                                          5
                                                                                                                                                       Author et al.:

                                      2nb -1                               |8 =                            1                      2nb -1        bjC    |uj
                                                                                                                                    j=0          ~j
|6 =                 1                     bj |uj  |~j      C        (23)                    2n/2          2nb -1    |  bj C  |2
                                                            ~j |1 a                                        j=0           ~j
           2nb -1           bj C      j=0
           j=0           |   ~j   |2                                                                   2n -1

                                                                                                   (         e-ijtye2iy~j/N |y ) |1 a                         (25)

                                                                                                       y=0

where the prefactor is due to normalization after mea-                        Since we earlier chose to set ~j = Njt/2, therefore,
                                                                           the two exponential terms cancel each other and
surement.  Since  |  C   |2  is   the  probabily        of  obtaining  |1
                     ~j
when the ancilla bit is measured, C should be chosen to

be as large as possible. Compared to (5), the result re-

sembles the answer |x that we are looking for. However,                                             1                    2nb -1      bjC             2n -1
                                                                                                                           j=0        ~j
we can only obtain the correct result if the b-register is                 |8 =                    2nb -1     bj C   |2                         |uj        |y  |1 a
                                                                                       2n/2        j=0         ~j
measured in the eigenvector basis (i.e. |uj instead of                                                     |                                         y=0
                                                                                 =
|0 / |1 ). However, the b-register is entangled with the                                                                 C                           2n -1
clock qubits, |~j . This means that we cannot factorize
                                                                                                                        2nb -1       bj C       |x          |y |1 a
                                                                                                                        j=0           j
the result into a tensor product of the c-register and b-                                              2n/2                       |         |2         y=0

register (See the discussion after (5) and Appendix V-6).                                                                                                      (26)

As a result, we cannot convert the b-register into the                        The clock qubits and the b-register are now unen-
                                                                           tangled and the b-register stores |x . By applying the
|0 / |1 measurement basis with the desired amplitudes.                     Hadamard gate on the clock qubits, finally, we have

We will need to uncompute the state so that it gives the

right results in the |0 / |1 measurement during which

the b-register and c-register will be unentangled.

   The measurement of the ancilla qubit can be and                         |9 =                     1                2nb -1       bjC      |uj     |0 n |1 a
is usually performed after uncomputation. However,                               =                                     j=0         j
since the ancilla bit is not involved in any operations                                         2nb -1  |  bj C  |2
after the controlled rotation, measuring the ancilla bit                                        j=0         j
before the uncomputation gives the same result. For
simplicity in the derivation, it is thus performed before                                           C                |x     b |0     n     |1                  (27)
the uncomputation.                                                                                                                   c
                                                                                                2nb -1  |  bj C  |2                             a
                                                                                                j=0         j

                                                                           If C is real and by using (7),

                                                                                 |9 =                      1                  |x        |0  n      |1
                                                                                                                                            c
E. UNCOMPUTATION - INVERSE QPE                                                                          2nb -1   |   bj  |2          b                 a
                                                                                                        j=0          j
The uncomputation is done by using inverse QPE.
Firstly, QFT is applied to the clock qubits as,                                              =  |x  b  |0  n     |1      a                                     (28)
                                                                                                           c

                                      2nb -1                                  Here, the solution |x (Eq. 5) is stored in the b-register
                                        j=0                                successfully.

|7 =              1                           bjC  |uj  QFT |~j      |1 a
                                               ~j
           2nb -1        |  bj C  |2                                       III. NUMERICAL EXAMPLE
           j=0               ~j                                            We will present a numerical example and apply HHL
                                                                           to it step-by-step. The implementation is shown in
=                 1                   2nb -1  bjC  |uj                     Figure 3. Firstly, we will discuss how to implement the
                                        j=0    ~j                          controlled-U and ancilla qubit rotations.
           2nb -1        |  bj C  |2
           j=0               ~j

     1     2n -1                       ) |1   a                      (24)  A. ENCODING SCHEME
 ( 2n/2                                                                   In this example, the matrix A and vector b are set to be
                 e2iy~j/N |y

           y=0

                                                                                                    A=           1          -  1                               (29)
                                                                                                                               3
   Then inverse controlled-rotations of the b-register by                                                            1
the clock qubits are applied with U-1 = e-iAt. Similar to                                                        -   3        1
the forward process, when the controlling rth clock qubit
is |0 , |uj will not be affected. If the rth clock qubit is                                                b=            0                                     (30)
|1 , (U-1)2r will be applied to |uj . This is equivalent                                                                 1
to multiplying e-ijty if the c-register is |y . This is due
to the similar argument in Eq. (14) and the fact that                                                                                   -1                     -1
2i = ijt. Therefore,                                                                                                                      2                      2
                                                                           The eigenvectors of A are u0 =                               -1         , u1 =      1
6
                                                                                                                                            2                  2

                                                                           with  eigenvalues    0      =   2     and        1     =     4   respectively.      We
                                                                                                           3                            3

                                                                                                                                                       VOLUME 3, 2022
Author et al.:

need to using basis encoding to encode the eigenvalues

in the basis formed by the clock qubit and 2 qubits are                                               -1      0
                                                                                                       0      1
needed by encoding 0 as |01 and 1 as |10 so that it                           Ud2iag = UdiagUdiag =                 (35)
maintains the ratio of 1/0 = 2. This means ~0 = 1 and
~1 = 2 or in other words, |~0 = |01 and |~1 = |10 .                          It is worth noting that both are naturally unitary
                                                                          which is a requirement for a quantum operation.
This gives a perfect encoding with n = 2 (i.e. N = 4).
                                                                             To obtain U and U2 in the original basis, we need
Therefore,      t  is  chosen  to  be     3   to  achieve  the  encoding  to apply similarity transformation again in the reverse
                                           4                              direction,
scheme since ~j = Njt/2.

Since b is a 2-dimensional complex vector, it can be

encoded using 1 qubit and, thus, nb = 1.

The solution to the LSP is found to be

                                          3                                   U = VUdiagV

                               x=         8                     (31)                -1 -1       i0    -1 -1
                                          9                                           2    2   0 -1       2      2
                                                                              =     -1   1
                                          8                                                    1+i    -1 1
                                                                                              -1 + i
whereby, the ratio of |x0|2 to |x1|2 is 1 : 9.                                      2    2                2      2

                                                                                 1  -1 + i                          (36)
                                                                              =      1+i
B. CONTROLLED-U IMPLEMENTATION
                                                                                 2
In reality, we expect the controlled-U operation to be
implemented by Hamiltonian simulation [18]. However,                      U2 = VU2diagV
to understand the algorithm, we will derive the matrix
for U and then map this to the Controlled - U gate used                             -1 -1     -1 0    -1 -1
in IBM-Q directly. Since n = 2, there are two operations                              2    2   01         2      2
needed, namely U21 = U2 and U20 = U, controlled by                            =     -1   1            -1 1
c1 and c0, respectively.
                                                                                    2    2                2      2
   In order to find the corresponding matrix for U2 =
ei2At and U = eiAt, we need to perform similarity                             =      0   -1                         (37)
transformation on i2At and iAt, exponentiate then, and                              -1    0
transforms back to the original basis.
                                                                             To implement U and U2, a 4-parameter arbitrary
   The transformation matrix from the original basis to                   unitary gate with global phase [19],
the eigenvector basis is
                                                                                ei cos(/2)    -ei(+) sin(/2)
                                                                          U=  ei(+) sin(/2)   ei(++) cos(/2)        (38)

                       V=          u0 u1
                          =
                                   -1 -1                                     By choosing  = ,  = ,  = 0,  = 0, U2 is
                                       2      2                           implemented.
                                   -1 1                         (32)
                                                                             By choosing  = /2,  = -/2,  = /2,  = 3/4,
                                       2      2                           U is implemented.

   Since V is real and symmetric,its conjugate, V equals                     For the IQPE part, we also need to implement U-1
itself.                                                                   and U-2. Since in this example, (U2)-1 = U2, one can
                                                                          use the same set of parameters to implement (U2)-1.
   The diagonalized A, i.e. expressed in the basis formed
by u0 and u0, is                                                             However,

                       Adiag = VAV

                                   =      2   0                 (33)             U-1 = 1      -1 - i   1-i          (39)
                                                                                          2    1-i    -1 - i
                                          3   4
                                              3
                                          0

As it is diagonal, U can be obtained by exponentiation                       We need to choose  = /2,  = /2,  = -/2,  =
                                                                          -3/4 to implement U-1.
of the elements accordingly.
                                                                             The controlled version of matrix U can then be
                       Udiag =     ei0t 0                                 constructed using
                              =      0 ei1t
                              =                                               C - U = I  |0 0| + U  |1 1|           (40)
                                   ei/2 0
                                     0 ei                       (34)         Note that in this equation, only the controlling clock
                                                                          bit and the b-register are included for simplicity. For
                                   i0                                     example,
                                   0 -1
                                                                                                                                                                                 7
where t = 3/4 as mentioned earlier is used. Also,

VOLUME 3, 2022
                                                                                                                                                 Author et al.:

                                                                                       The following function can achieve the goal,

C - U-1 =     1  0          1    0        1     -1 - i           1-i                                                            
              0  1          0    0     +         1-i            -1 - i
                                                                                                        (c) = (c1c0) = 3 c1 + c0                         (47)
                                          2

              0  0
              0  1
                                                                                      Therefore, the controlled rotation can be implemented
           1 0 0 0 0 0                                             0 0             as

   =       0      0      0  0    +  1  0         -1 - i            0 1-i 
           0      0      1  0       2  0            0
                                                                   0      0                                          
                                                                                               |1 1|  I  RY( ) + |0 0|  I  I+

              0000                     0 1-i                       0 -1 - i                                          3
                                                                                                  I  |1 1|  RY() + I  |0 0|  I (48)
              2 0 0 0 
                                                                                      where the operators operate on qubits |c1 , |c0 , and
   =       1  0      -1 - i      0     1-i                                (41)     |a from left to right, respectively.
           2  0         0        2
                                       0         
                                                 

                 0 1 - i 0 -1 - i

C. IMPLEMENTATAION OF THE CONTROLLED-ROTATION                                      D. QUANTUM CIRCUIT
                                                                                   An HHL circuit for the numerical example is then built
OF ANCILLA QUBIT                                                                   and shown in Figure 3. We will then walk through the
                                                                                   circuit using numerical substitution.
The coefficients of |0 and |1 of the ancilla bit after ro-

tation in Eq. (22) are      1-   C2    and   C   ,    respectively.       The
                                 ~j 2        ~j

sum of the square of the magnitudes of the coefficients is

1 as required. This means also C  ~j. Since the minimal                            E. NUMERICAL SUBSTITUTION
~j is 1, we will set C = 1 to maximize the probability of                          The algorithm begins with

measuring |1 during the ancilla bit measurement.

   The transformation of |0 a to             1  -      1    |0  a  +  1   |1  a                   |0 = |0 b  |00 c  |0 a = |0000                         (49)
                                                      ~j 2            ~j

is known to be equivalent to RY() rotation,

                                                                                   X-gate is then applied to convert |0 b to |1 b with

           RY() =           cos        - sin                              (42)                       |1 = X  I  I |0 = |1000                             (50)
                                 2                    2
                                                   
                            sin  2     cos         2

with  = 2 arcsin     1      . One can check this by multiply-                         After applying the Hadamard gates to create a super-
                     ~j                                                            position among the clock qubits,

ing RY() to |0 a.

                                                                                       |2 = I  Hn  I |1
                                 2                  2
   RY() |0 a =              cos        - sin                    1                                       1          + |1 )2 |0
                                                                0
                            sin        cos                                                     = |1     22    (|0
                                 2               2                                                         2

                                                                                                  1
                    = cos 2 |0 a + sin 2 |1 a (43)                                             = (|1000 + |1010 + |1100 + |1110 ) (51)

                                                                                                  2

   Therefore, we will establish a function to implement                                Before applying the CU3(controlled rotation of an-
this rotation and this function only need to be valid
when the input are the encoded eigenvalues because                                 cillary qubit) gates in the bra-ket notation, it will be
only encoded eigenvalues have zero magnitudes in the
c-register as shown in (21). The function is defined as                            convenient to perform a basis change to the eigenvector

                                                                                   basis of A. Since |1             =      1    (-  |u0  + |u1 ), we have
                                                                                                                             2
                                                                                          -1                   1 .
                                                                                   b0  =    2  and   b1    =        Therefore,
                                                                                                                 2

                                                       1                                                1
           (c) = (c1c0) = 2 arcsin c                                      (44)
                                                                                       |2      = |1      (|000      + |010      + |100   + |110 )
                                                                                                        2

   where c is the value of the clock qubits and c1c0 is its                                       1                         1
binary form.                                                                                   =        (-    |u0   + |u1  ) (|000       + |010  + |100
                                                                                                    2                       2
   Since only |~j has non-zero amplitude in (21), we only
need to set up (44) such that it is correct for |c =|01                                        + |110 )

and |10 , namely                                                                                     1
                                                                                               =           (-  |u0  |000    - |u0   |010  - |u0    |100
                                                                                                  22

                                                1                                              - |u0 |110 + |u1 |000 + |u1 |010

           (1) = (01) = 2 arcsin =                                        (45)                 + |u1 |100 + |u1 |110 )                                   (52)
                                                1

                                             1                            (46)        In the controlled rotation operations, when the cor-
           (2) = (10) = 2 arcsin =                                                 responding c-register is |k c, a phase change of j =
                                             23
                                                                                                                                                                                      VOLUME 3, 2022
8
Author et al.:

FIGURE 3. The HHL circuit corresponding to the numerical example built to run in IBM-Q. This circuit is partitioned in the same way as the general HHL
schematic in Figure 1.

kjt/2 is added (i.e. multiplied by e2ij ) for |uj . Since              Now apply IQFT to the clock qubits, e.g.

t  =  3   ,  0  =  2  and  1  =  4  ,  we  have                        IQFT |10 = IQFT |2
       4           3             3

|3    =      1   (- |u0    |000  - e2i0 |u0      |010                                   =        1  22 -1
             
             22                                                                             22/2         e-2i2y/4 |y

      - e2i20 |u0 |100 - e2i30 |u0 |110 +                                                           y=0

      |u1 |000 + e2i1 |u1 |010                                                             1
                                                                                        = (|0 - |1 + |2 - |3
      + e2i21 |u1 |100 + e2i31 |u1 |110 )
                                                                                           2
             1                                                                             1
                                                                                        = (|00 - |01 + |10 - |11 ) (55)
                                                                                           2

      =      22  (- |u0    |000  - ei0t |u0  |010                      Similarly,

      - ei20t |u0 |100 - ei30t |u0 |110 +                                               1
                                                                       IQFT |00 = (|00 + |01 + |10 + |11 ) (56)
      |u1 |000 + ei1t |u1 |010
                                                                                        2

      + ei21t |u1 |100 + ei31t |u1 |110 )                                               1
                                                                       IQFT |01 = (|00 - i |01 - |10 + i |11 ) (57)
             1                   - ei/2 |u0            - ei |u0
      =          (-   |u0  |000                  |010            |100                   2
             22
                                                                                              1
      - ei3/2 |u0 |110 + |u1 |000 + ei |u1 |010                              IQFT |11 = (|00 + i |01 - |10 - i |11 ) (58)

      + ei2 |u1 |100 + ei3 |u1 |110 )                                                         2

             1                                                         Therefore, applying IQFT to |3 and substituting
      =          (-   |u0  |000  - i |u0   |010  + |u0  |100           Eq. (55) to (58),
             22

      + i |u0 |110 + |u1 |000 - |u1 |010 +                             |4 = IQFT |3
                                                                             = 1
      |u1 |100 - |u1 |110 )                                      (53)            42

   Before applying IQFT, the terms are regrouped for                   ((- |u0 + |u1 )(|00 + |01 + |10 + |11 )+
simplicity.
                                                                       (-i |u0 - |u1 )(|00 - i |01 - |10 + i |11 )+

                                                                       (|u0 + |u1 )(|00 - |01 + |10 - |11 )+

                1                                                      (i |u0 - |u1 )(|00 + i |01 - |10 - i |11 )) |0
   |3     =         ((-  |u0  + |u1 ) |00  + (-i |u0    - |u1 ) |01
                22                                                                 1
                                                                       =                (-  |u0  |01  + |u1  |10 ) |0  (59)
          + (|u0 + |u1 ) |10 + (i |u0 - |u1 ) |11 ) |0 (54)                          2

VOLUME 3, 2022                                                                                                                                          9
                                                                                                                                           Author et al.:

It can be seen that after IQFT, the eigenvalues are                          e-i0t = -i, e-i20t = -1, e-i1t = -1, e-i21t = 1, and
                                                                             Nt/2=3/2
encoded in the clock qubits as |01 and |11 with non-

zero  amplitudes     due    constructive     interference.       b0  =  -1                81                   1
                                                                          2
             1 .
and   b1  =       We    clearly    see  the  entanglement        between     |8 =         5  (-        |u0      (|00   + |01       + |10   + |11 ) |1
               2                                                                                   2           2

the b-register and the c-register that |u0 goes with |01                                  1         1

and |u1 goes with |11 .                                                             +  |u1           (|00      + |01      + |10       + |11 )) |1
                                                                                       22           2
After performing the ancilla qubit rotation,
                                                                                    18 1                               1
                                                                                    =           (-       |u0      +          |u1   )(|00   + |01
                                                                                       2     5        2              22
                21 -1
                                                C2            C                     + |10 + |11 ) |1
      |5 =             bj |uj |~j (        1 - ~j2 |0 + ~j |1 )

                j=0                                                                 12 8 1                                      1
                                                                                    = ()            (-            |u0  +              |u1  )(|00   + |01
                                                                                       23        5       2    2              4     2
                  1                             1          1                                             3                   3
             =  -       |u0   |01  (    1 - 12 |0      + |1 )+
                     2                                    1                         + |10 + |11 ) |1                                               (65)

             1                          1           1                           Finally, by applying Hadamard gate to the clock
              |u1 |10 (          1 - 22 |0      + |1 )               (60)    qubits,
               2                                   2

If the measurement of the ancilla bit is |1 ,                                                               28 1
                                                                                                 |9    =             (-            |u0
                                                                                                          3       5       2    2
                                                                                                                          3

                                                                                                               1
                              81                                                                       +             |u1  ) |00    |1              (66)
                |6 =             (-        |u0  |01    |1                                                   4     2
                              5         2                                                                   3

                             1                                                  It can be verified that |9 is a normalized vector as
                        +  |u1 |10 |1 )                              (61)    it should be because every operation in the HHL circuit
                            22

                                                                             is unitary and preserves the norm.

Applying QFT to the encoded eigenvalues, we have                             Equation (66) can be simplified by substituting |u0 =

                                                                             -1 |0  + -1 |1      and |u1       = -1 |0       + 1 |1 . We obtain,

                                                                               2           2                          2             2

      QFT |10 = QFT |2                                                                                 12                                          (67)
                                                                                             |9  =           (|0 + 3 |1 ) |00 |1
                              22 -1                                                                 2       5
                        1
                  =                e2i2y/4 |y
                        22/2                                                    The probability ratio of obtaining |0 and |1 when
                                y=0                                          b-register is measured is thus 1 : 9 as expected.

                     1                                               (62)    F. SIMULATION RESULTS
                  = (|00 - |01 + |10 - |11 )

                     2

                                                                             Matlab code implementing the numerical example using
                                                                             matrix approach is created and available at [21]. In the

      QFT |01 = QFT |1                                                       Matlab code, measurement is not performed (i.e. not
                       1
                                                                             partical tracing of the matrix). 9 is found to be,
                   = (|00 + i |01 - |10 - i |11 ) (63)
                       2                                                                                          -0.4330

   Therefore, applying QFT to |6 and substituting                                                                  0.2500 
Eq. (62) to (63), we obtain
                                                                                                                     0.0000     
                                                                                                                                

                                                                                                                  -0.0000

                                                                                                                     0.0000     
                                                                                                                                

                                                                                                                  -0.0000

                     81                 1                                                                            0.0000     
                                                                                                                                

      |7 =           5  (-       |u0     (|00   + i |01    - |10                                    |9      =        0.0000                        (68)
                              2         2                                                                            0.4330     
                                                                                                                                
                                                                                                                                

                                   1            1                                                                    0.7500     
             - i |11 ) |1     +  |u1             (|00  - |01                                                                    
                                 22             2
                                                                                                                  -0.0000

             + |10 - |11 )) |1                                       (64)                                            0.0000     
                                                                                                                                

                                                                                                                  -0.0000

                                                                                                                     0.0000     
                                                                                                                                

   For the controlled rotation, the state is multiplied by                                                        -0.0000
e-ijt and e-i2jt if c0 = 1 and c1 = 1, respectively. Since
                                                                                                                     0.0000
10
                                                                                                                                           VOLUME 3, 2022
Author et al.:

Since |0 c are discarded during the measurement step,                       example is also presented in the bra-ket notation. The
only |0001 and |1001 are left. Their amplitude ratio is                     numerical example echos the analytical derivation to
0.252 : 0.752 = 1 : 9 as expected.                                          help students understand how qubits evolve in this
                                                                            important and relatively complex algorithm. A Matlab
   The circuit in Fig. 3 is also simulated in the IBM-Q                     code corresponding to the numerical example is con-
system (code available at [21]). Since only the b-register                  structed to help understand the algorithm from the ma-
and the ancilla qubit are measured, there are only four                     trix point of view. Qiskit circuit of the corresponding cir-
possible outputs as shown in Figure 4. Again, only |1 a                     cuit which can be simulated in IBM-Q and run on their
should be considered. The ratio of the measurement                          quantum computing hardware is also available. Through
probability of |0 b |1 a to |1 b |1 a is 0.063 : 0.564 = 1 :                this self-contained and step-by-step walkthrough, the
8.95, which is close to the expected value.                                 basic concepts in quantum computing are reinforced.

FIGURE 4. Simulation result of the circuit in Figure 3 using IBM - Q. Only  V. APPENDIX
the MSB | b, and the LSB | a are measured.                                  1) Hermitian matrix
                                                                            A Hermitian matrix is a matrix that equals to its adjoint
   On the other hand, due to the imperfection and noise                     matrix (transpose followed by complex conjugation).
in a real quantum computer, the hardware execution                          That is, if A is a Hermitian matrix, then it is defined
of the same circuit does not give a satisfactory result                     as,
(Figure 5). The ratio of the measurement probability of
|0 b |1 a to |1 b |1 a is only 0.1422 : 0.3612 = 1 : 2.54.                                       A = A = (AT)        (V.69)

FIGURE 5. Hardware result of the circuit in Figure 3 run in machine         where AT is the transpose of A.
ibmq_santiago. Only the MSB, | b and the LSB, | a are measured.                In this paper, the matrix, A, in the LPS to be solved

IV. CONCLUSION                                                              is assumed to be Hermitian.
In this paper, we presented the HHL algorithm through a                        Another example is in (32), where V is Hermitian.
step-by-step walkthrough of the derivation. A numerical
                                                                                                    V = u0 u1
VOLUME 3, 2022
                                                                                                           -1 -1
                                                                                                             2    2
                                                                                                    =      -1   1    (32)

                                                                                                           2    2

                                                                            2) Bra-ket Notation

                                                                            Bra-ket notation is commonly used in quantum mechan-

                                                                            ics. A vector v is represented as |v in its ket form. The

                                                                            bra form of the vectors forms a dual space to the space

                                                                            of the kets. The bra form of v is v|.

                                                                            In matrix representation, ket is the complex conjugate

                                                                            transpose of bra and vice versa. For example, if |v =

                                                                            1   , then  v| =     1     i.
                                                                            -i

                                                                            3) Superposition

                                                                            Superposition or Quantum Superposition is a quantum
                                                                            state which is the linear combination of two or more
                                                                            basis states. For example, a superposition state can be
                                                                            |v = c1 |1 +c2 |0 , where c1 and c2 are complex number
                                                                            and |1 and |0 are basis states. A Hadamard gate is
                                                                            a gate commonly used to create a superposition state
                                                                            (Appendix V-5).

                                                                            4) Basis Transformation and Quantum Gate
                                                                            In quantum computing, we only care about the basis
                                                                            transformation due to rotation in the hyperspace. The
                                                                            transformation is equivalent to the multiplication of
                                                                            the basis vectors by a unitary matrix, U, which is
                                                                            the transformation matrix. All quantum gates can be
                                                                            defined as how the basis vectors are transformed from

                                                                                                                                                                                 11
                                                                                                                               Author et al.:

the initial basis vector to the final basis vectors. Usually,            In this paper, after the ancilla bit rotation, we have
a quantum gate rotates a basis state into another basis
state (e.g. the NOT gate) and has its classical counter-                       81                                 1
part. But there are some gates that rotate a basis state                 |6 =     (-       |u0    |01  |1   +        |u1  |10  |1 )  (61)
to a superposition of two or more basis states. Such gates                     5        2                      22
have no classical counterparts. For example, a Hadamard
gate defines how an initial basis vector is rotated to an                where the b-register and the c-register are entangled
equal superposition of two basis vectors (Appendix V-5).                 and |u0 (|u1 ) always appears with |01 (|10 ) after the
                                                                         measurement.
5) Hadamard Gate
                                                                            If the b-register were not entangled with the c-register,
                                                                         we have

The Hadamard gate is a quantum gate that does not

have a classical counterpart. It rotates the basis state to                          81                        1

create an equal superposition of the basis states. For a                       |6 =  5     (-       |u0    +  |u1 )             (V.75)
                                                                                                 2            22
1-qubit  case,  this  means    it  has  equal  probability  (i.e.  1  )
                                                                   2

of measuring |0 and |1 .                                                 By substituting |u0      = -1 |0 + -1 |1         and |u1    =

    The matrix form of the Hadamard gate is,                                                             2           2
                                                                         -1 |0 + 1 |1 and after simplification, we have
                                                                         2     2
                       1           1   1
                                   1  -1                    (V.70)

                        2                                                      8 1 -1                      -1             1 -1

                                                                   0     |6 =     (-  (  |0 +  |1 ) +  (  |0 +
                                                                   1           5 22                         2           22 2
When it is applied on the basis state |1 , which is
                                                                            1 |1 ))
in matrix form, we have,                                                      2

                          1    1       1  0                 (V.71)             81              3                                (V.76)
                      =        1      -1  1                 (V.72)          = ( |0 + |1 )
                                                                               54              4
                            2

                                   1       1                                This is the same as (67). The probability of measuring
                               =          -1                             |0 and |1 has the ratio of 1:9 as expected.

                                     2                                      However, when there is entanglement, the probability
                                                                         of measuring |0 and |1 would not be 1:9 because the
    which can also be represented in bra-ket form as,                    previous grouping is impossible.

                               |0 - |1                      (V.73)
                                    2

   In this paper, Hadamard gates are applied in clock                          8 1 -1                    -1               1 -1
qubit to create superposition from |1 to |2 . For                        |6 =     (-  (  |0 +  |1 ) |01 +  (  |0
example, in (11),                                                              5 22                         2             22 2

                |2 = Inb  Hn  I |1                          (11)            + 1 |1 )) |10                                       (V.77)
                                                                                 2
where |2 is obtained by applying tensor product of
identity gates and an n-qubit Hadamard gate to |1 .                      7) Controlled Operation
The identity gates are applied to the b-register and the
ancilla qubit while the Hadamard gate is applied to the                  Controlled operation requires more than one qubit. For
clock qubits. In this equation, the n-qubit Hadamard                     a 2-qubit controlled gate, an operation is applied to a
gate is represented as Hn, i.e. the tensor product of n                  qubit (the target qubit), if the value of the controlling
1-qubit Hadamard gates.                                                  qubit is 1 in the basis vector.

6) Entanglement                                                             For example, in Figure 2, |b is the target qubit and
Entanglement refers to the quantum state of a 2- or                      |cn-1 is the controlling qubit. The operation of U2n-1
more-qubit system that cannot be expressed as a tensor                   is applied to |b only if |cn-1 is 1 in the basis state (e.g.
product of the individual qubit. This is an important                    |bcn-1    = |01    .
feature that quantum computing uses often. As an
example,                                                                    In general, the controlled version of a unitary gate,
                                                                         U , can be implemented using the following equation if
                                                                         the LSB is the controlling qubit.

                | = 1 (|00 + |11 )                                             C - U = I  |0 0| + U  |1 1|                           (40)
                           2
                                                            (V.74)

is an entangled state. It cannot be expressed as a tensor                which literally means that if the controlling qubit is
product of two individual qubit states.                                  0, Identity gate is applied to the target qubit (MSB).
                                                                         Otherise, U is applied to the target qubit.
12
                                                                                                                                                                            VOLUME 3, 2022
Author et al.:

8) Eigenvalue and Eigenvector                                    Hamiltonian in a unitary gate. For example, in this
When a non-zero n  n matrix A is applied to an n-               paper, Eq. (17) shows that
dimensional vector V and has the following relationship,
                                                                                               U = eiAt                 (17)

                          AV = V                   (V.78)        where it encodes matrix A as the Hamiltonian of the
                                                                 unitary gate U. Matrix A needs to be Hermitian as
where  is a scalar, then, by definition, V and  are              it is used to represent the Hamiltonian (the energy)
the eigenvector and eigenvalue of A, respectively. This is       of the system. However, A does not need to be
similar to (3), where the matrix A is expressed as a linear      unitary and U will be unitary due to its definition
combination of the outer products of its eigenvectors,           in (17).
|ui ui|.

                          2nb -1                             10) Discrete Fourier Transform (DFT)
                                                             The discrete Fourier Transform (DFT) transforms an N-
                   A=             i |ui ui|            (3)   dimensional vector X to another N-dimensional vector
                                                             Y. The transformation matrix  contains the powers of
                          i=0                                the N-th root of unity,  = ei2/N. The transformation
                                                             is represented as,
   This can be checked by applying A to its eigenvector
|uj ,

                          2nb -1

                A |uj  =          i |ui   ui| |uj  (V.79)              Y=                                                 X
                                         ij
                            i=0                               y0                    -00              -0(N-1) 
                          2nb -1
                                                              y1            1  -10                   -1(N-1) 
                       =          i |ui
                                                               ...     =                  ...       ...          ...    
                                                                              N                                         
                          i=0                                                                                           

                       = j |uj                               yN-1                  -(N-1)0    -(N-1)(N-1)

which meets the definition in Eq. (V.78).                                 x0 

                                                                          x1                                            (V.83)
                                                                         
9) Different Types of Encoding                                                ...  
The three common types of encodings are explained                                  
here.                                                                              

    Basis Encoding- Basis encoding converts classical                      xN-1
      information such as numbers or matrix to quantum
      information in the form of basis states. For example,  11) Inverse Quantum Fourier Transform (IQFT) and
                                                             Quantum Fourier Transform(QFT)
x = 2 binary 10 quantum state |10
                   ----           ----------                 Mathematically, IQFT is similar to DFT (See Ap-
                                                             pendix V-10). The transformation matrix, UI, is N  N
                                                             for an N-dimensional Hilbert space. Therefore, N = 2n
                                                             for an n-qubit system. Eq. (V.83) becomes

                2          10                                                             |Y = UI |X                    (V.84)
                3          11
x=                 binary         quantum state |1011        and UI has the same expression as  in DFT.
                   ----           ----------

                                                   (V.80)                    -00                     -0(N-1) 

 Amplitude Encoding- Amplitude encoding encodes                       1  -10                        -1(N-1) 
   the information as the coefficients of the basis
   vectors. For example, for                                 UI     =              ...         ...       ...            (V.85)
                                                                         N                                            
                                                                                                                      

                                                                            -(N-1)0    -(N-1)(N-1)

                           v=     v0               (V.81)       Note that in some literature, e.g. [19], this form of
                                  v1                         IQFT is called QFT. |X and |Y are the quantum
                                                             states/vectors in the N-dimensional Hilbert space. IQFT
which is assumed to be normalized (|v| = 1), it can          can be treated as the rotation of |X to |Y .
be encoded as in the following quantum state,
                                                                If |X is a basis vector |k , applying IQFT to |k using
                       |v = v0 |0 + v1 |1          (V.82)    Eq. (V.85), we have

   where v0 and v1 become the coefficients of the basis                            UI |k       1    N-1                 (V.86)
   states, |0 and |0 , respectively. In the main text,                                    =
   Eq. (9) shows how the values of the components of                                                     -jk |j
   vector |b are encoded using amplitude encoding.                                             N j=0
 Hamiltonian Encoding- One type of the Hamil-
   tonian encoding is to encode the matrix as the               This is the equation used often in this paper. It tells
                                                             us that by applying IQFT to a basis vector, the basis

VOLUME 3, 2022                                                                                                            13
                                                                                                                                             Author et al.:

vector is rotated to a superposition of all basis vectors
weighted by the powers of the N-th root of unity.

   For example, in (55) in the main text,

        IQFT |10 = IQFT |2

                        =  1       22 -1                                      FIGURE 6. Implementation of a 2qubit inverse quantum Fourier
                                                                              transformation.
                           22/2         e-2i2y/4 |y

                                    y=0

                           1                                                      In  general,   the  phase      shift  angle  is      =     -2   ,  r-1
                        = (|0 - |1 + |2 - |3                                                                                                  2r

                           2                                                  is the distance between the controlling qubit and the
                           1
                        = (|00 - |01 + |10 - |11 ) (55)                       target qubit. For the 2-qubit IQFT case, there is only
                           2

where N = 22 = 4, k = 2, j = y in (V.86) is used.                             one controlled phase shift gate and r = 2 and this results
The basis state |10 becomes a superposition of all other
basis states after IQFT.                                                      in  the  phase     =    -   .
                                                                                                       2
   Another more complex example is the general equa-
tion,in (15), for the IQFT in Figure 1.                                           For QFT, the circuit is the same as the IQFT, but

                                                                              the  phase   shift  is    negated,  i.e.        =    2   .  This    can   be
                                                                                                                                   2r

                                                                              appreciated by the fact that the elements in the IQFT

                                                                              and QFT have opposition signs in (V.85) and (V.87),

                                                                              respectively.

    |4     = |b  IQFT(          1   2n -1          ) |0

                              2  n        e2ik |k        a                    13) Gaussian elimination method
                                 2                                            Here, Gaussian elimination is demonstrated by solving
                                                                              (1) using the numerical example in Section III.
                                    k=0

           = |b   1        2n -1                   ) |0  a
                 2n
                                 e2ik(IQFT |k
                    2
                           k=0                                                                               Ax = b                                     (1)

                 1      2n -1            2n -1
                 2n
           = |b                 e2ik(           e-2iyk/N |y    ) |0  a                              1        -1   x0       =     0                (V.89)
                                                                                                                  x1             1
                           k=0           y=0                                                        -1        3
                                                                                                     3
                                                                                                             1

               1        2n-1 2n-1                        |0 a           (15)  which is rewritten as an augmented matrix followed by
           = 2n |b
                                    e2ik(-y/N) |y                             Gaussian method of Elimination to solve for x0 and x1,

                           y=0 k=0

   Here, IQFT is applied to the c-register which is a                             1    -1     0  -R-o-w-2--=--3---R-o-w-2--+-R--o-w1     1  -1         0
superposition of basis states, |k . Using the distribution                                    1                                           0   3         3
law of matrix operations, IQFT is applied to individual                           -1    3                                                     8
|k and (V.86) is used with y = j and N = 2n.                                       3                                                          3
                                                                                       1
   QFT is the inverse of IQFT and can be treated as the
rotation of the basis. The rotation matrix is given by                                                3                 3     -1       0
                                                                                           Row2 =  Row2                0      1
                                                                                           --------8-----                              9
                                                                                                                                       8

                1 1                             1                               3   -1     0             1                             1  0       3
                                                                                   0    1        Row1 =  Row1 + Row2                     0  1       8
                                                                                              9  --------3------------                               9
                                                                                              8                                                      8

                1 1                          (N-1)                         This the solution x is

    UQ  =  1    1          2                 2(N-1)         (V.87)                                        x0         3
             N              ...          ...                                                                 x1         8
                                                                                                                        9
                   ...                             ...                                                           =      8
                                                            
                                                            

                  1 (N-1)    (N-1)(N-1)                                       The complexity of Gaussian Elimination is O(N3).
                                                                              This is much slower than the classical conjugate gradient
    Equivalently,

                                         N-1                                  method (Appendix V-14), to which HHL is compared.

                   UQ |k              1       jk |j            (V.88)         14) Conjugate Gradient Method
                                 =                                            The Conjugate Gradient Method (CGM) solves the LSP
                                         N j=0                                with a complexity of O(N) and is the fastest known
                                                                              classical solver. Therefore, the speed of HHL, which has
    It can be shown that UI = UQ-1 or UIUQ = I.                               a complexity of O(log(N)), is often compared to the
                                                                              speed of CGM [22]. Thus, HHL provides an exponential
12) Implementation of QFT and IQFT                                            speedup over the fastest known classical method.
QFT and IQFT are constructed using Hadamard gates,
controlled phase shift gates, and SWAP gates. Readers                            When we solve a system of linear equation, according
may refer to other sources for the details (e.g. [19]). Here,                 to (Eq.(1)),where A is a matrix, b is a vector and x is to
we show the circuit of a 2-qubit IQFT gate (Fig. 6).
                                                                                                                                                                                 VOLUME 3, 2022
14
Author et al.:

be solved. If A is a 22 matrix and b is 21 ,then x can              [10] S. Dmitry et al., "The Potential of Quantum Computing and
be solved easily. But is A is a long matrix, for example
1000, 000, 000  1000, 000, 000 and b is 1000, 000, 000  1           Machine Learning to Advance Clinical Research and Change
vector, and N in this case is 1000, 000, 000. To solve
x in Classical Gaussian Elimination method we need                    the Practice of Medicine." Missouri medicine 115 (5), 463467
O(N3) speed, where O is omega. In Classical Conjugate
Gradient Method with sparse matrix that contains lots                 (2018).
of zeros, it will take O(N) speed. And for HHL Quantum
Algorithm, with sparse matrix, it takes O(log(N)) speed.              [11] Bojia Duan, Jiabin Yuan, Chao-Hua Yu, Jianbang Huang,

   However, according to the paper [24], the speed of                 and Chang-Yu Hsieh, "A survey on HHL algorithm: From
inner product in HHL Quantum Algorithm is only
log(mn)/ steps when certain amplitudes is obtained                    theory to application in quantum machine learning", Physics
and distinguished among other amplitudes, otherwise
the speed is only quadratically faster than classical                 Letters A 384, 126595 (2020).
algorithm.
                                                                      [12] Shengbin Wang, Zhimin Wang, Wendong Li, Lixin Fan,
   To solve (1) in CGM method,
                                                                      Zhiqiang Wei, and Yongjian Gu, "Quantum fast Poisson

                                                                      solver: the algorithm and complete and modular circuit de-

                                                                      sign," Quantum Information Processing 19, Article number:

                                                                      170 (2020).

                                                                      [13] H. J. Morrell and H. Y. Wong, "Study of using Quantum

                                                                      Computer to Solve Poisson Equation in Gate Insulators,"

                                                                      2021 International Conference on Simulation of Semiconduc-

                                                                      tor Processes and Devices (SISPAD), 2021, pp. 69-72, doi:

                                                                      10.1109/SISPAD54002.2021.9592604.

                                                                      [14] Schleich, P., 2019. How to solve a linear system of

                                                                      equations using a quantum computer. [Online]. Avail-

                                                                      able: http://www.acom.rwth-aachen.de/_media/3teaching/

                                                                      00projects/schleich.pdf.

                Ax = b        (1)                                     [15] HHL Example using Qiskit. [Online]. Available:

                                                                      https://qiskit.org/textbook/ch-applications/hhl_tutorial.html.

initial guess of x is used as the starting point. The resid-          [16] Gadi Aleksandrowicz, et al., (2019). Qiskit: An Open-
ual is then found and the search direction is determined
by finding the steepest descent. This is repeated until a             source Framework for Quantum Computing (0.7.2). Zenodo.
stable condition is met.
                                                                      doi.org/10.5281/zenodo.2562111.
   The residual in the ksearch is given as,
                                                                      [17] IBM  Quantum         Site.    [Online].  Available:

                                                                      https://quantumcomputing.ibm.com/.

                                                                      [18] Dominic W. Berry, Graeme Ahokas, Richard Cleve, and

                                                                      Barry C. Sanders, "Efficient quantum algorithms for simu-

                                                                      lating sparse Hamiltonians," arXiv:quant-ph/0508139, 2007.

                Rk = b - Axk  (V.90)                                  [19] Hiu Yung Wong, Introduction to Quantum Computing: From

                                                                      a Layperson to a Programmer in 30 Steps. Switzerland:

   The readers do not need to understand CGM to                       Springer Nature, 2022, pp. 170. doi.org/10.1007/978-3-030-
understand HHL. Interested readers may refer to the
literature (e.g. [9]) for more details.                               98339-0. ISBN-10: 3030983382.

                                                                      [20] Danial Dervovic, Mark Herbster, Peter Mountney, Simone

                                                                      Severini, Nari Usher, and Leonard Wossnig, "Quantum lin-

                                                                      ear systems algorithms: a primer," arXiv:1802.08227v1.

                                                                      [21] Matlab code and Jupyter Notebook, [Online]. Available:

REFERENCES                                                            https://github.com/hywong2/HHL_Example.

 [1] P. W. Shor, "Algorithms for quantum computation: discrete        [22] Vandenbrocque, Adrien. "On Quantum Algorithms
      logarithms and factoring," Proceedings 35th Annual Sympo-
      sium on Foundations of Computer Science, 1994, pp. 124134,     for Solving Linear Systems of Equations."Master's
      doi: 10.1109/SFCS.1994.365700.
                                                                      semester Project I", 2019. [Online]. Available:
 [2] C. Outeiral, M. Strahm, J. Shi, G. M. Morris, S. C. Benjamin,
      and C. M. Deane, "The prospects of quantum computing            https://adrienvdb.com/projects-and-reports/.
      in computational molecular biology," WIREs Comput. Mol.
      Sci., 11, e1481 (2021).                                         [23] National Strategic Overview for Quantum In-

 [3] D. J. Egger et al., "Quantum Computing for Finance: State-       formation Science Report (National Science and
      of-the-Art and Future Prospects," in IEEE Transactions on
      Quantum Engineering, 1, pp. 124, 2020, Art no. 3101724,        Technology Council, 2018). [Online]. Available:
      doi: 10.1109/TQE.2020.3030314.
                                                                      https://www.quantum.gov/wp-content/uploads/2020/10/
 [4] Arute, F., Arya, K., Babbush, R. et al. Quantum supremacy
      using a programmable superconducting processor. Nature          2018_NSTC_National_Strategic _Overview_QIS.pdf.
      574, 505510 (2019). doi.org/10.1038/s41586-019-1666-5
                                                                      [24] Aaronson, S. Read the fine print. Nature
 [5] Madsen, L.S., Laudenbach, F., Askarani, M.F. et al. Quan-
      tum computational advantage with a programmable photonic        Phys 11, 291293 (2015). [Online]. Available:
      processor. Nature 606, 7581 (2022). doi.org/10.1038/s41586-
      022-04725-x                                                     https://doi.org/10.1038/nphys3272

 [6] N. David Mermin, "Could Feynman Have Said This?,"                                                                        15
      Physics Today 57 (5), 10 (2004); doi: 10.1063/1.1768652

 [7] A. Harrow, A. Hassidim, and S. Lloyd, "Quantum algorithm
      for linear systems of equations," Phys. Rev. Lett. 103, 150502
      (2009).

 [8] Yudong Cao, Anmer Daskin, Steven Frankel, and Sabre
      Kais, "Quantum circuit design for solving linear systems of
      equations," Molecular Physics 110, 1516 (2011).

 [9] R. Chandra, S.C. Eisenstat, and M.H. Schultz, "Conjugate
      Gradient Methods for Partial Differential Equations," in the
      Proceedings of the AICA International Symposium on Com-
      puter Methods for Partial Differential Equations, Bethlehem,
      Pennsylvania, June 1975.

VOLUME 3, 2022
