Compiling basic linear algebra subroutines for quantum computers
Liming Zhao Zhikuan Zhao Patrick Rebentrost Joseph Fitzsimons





Singapore University of Technology and Design,  8 Somapah Road, Singapore 487372
Centre for Quantum Technologies, National University of Singapore, 3 Science Drive 2, Singapore 117543

Singapore University of Technology and Design,  8 Somapah Road, Singapore 487372
Centre for Quantum Technologies, National University of Singapore, 3 Science Drive 2, Singapore 117543

Centre for Quantum Technologies, National University of Singapore, 3 Science Drive 2, Singapore 117543

Singapore University of Technology and Design,  8 Somapah Road, Singapore 487372
Centre for Quantum Technologies, National University of Singapore, 3 Science Drive 2, Singapore 117543
Horizon Quantum Computing, 79 Ayer Rajah Crescent, BASH #03-01, Singapore 139955







Efficiently processing basic linear algebra subroutines is of great importance for a wide range of computational problems. In this paper, we consider techniques to implement matrix functions on a quantum computer, which are composed of  basic matrix operations on a set of matrices.
These matrix operations include addition, multiplication, Kronecker sum, tensor product,  Hadamard product, and single-matrix functions.  We discuss the composed matrix functions in terms of the estimation of scalar quantities such as inner products,  trace,  determinant and Schatten p-norms. We thus provide a framework for compiling instructions for linear algebraic computations into gate sequences on actual quantum computers. 


Compiling basic linear algebra subroutines for quantum computers
    Joseph Fitzsimons
    December 30, 2023
================================================================



Quantum algorithms promise exponential speed-ups over their best known classical counterparts for certain problems. An example is the seminal Shor's algorithm that is able to find the prime factors of a given integer efficiently <cit.>. Because of such breakthroughs, substantial effort has been devoted to the study of quantum computation over the last several decades. There already exist several quantum processors with a limited number of qubits and ongoing efforts promise to increase the quantity, quality, and connectivity of the qubits. As the impressive development of quantum computers continues, it is highly desirable to study quantum algorithms for a wide range of applications.   

Matrix algebra is widely used in almost every area of science and technology. In many practical problems, the size of the input data, often in the form of vectors or matrices, is a bottleneck for efficient computations. The usual classical algorithms, for example for finding eigenvalues, can be unsuitable for large matrices. Efficient quantum algorithms for linear algebra and matrix inversion have been investigated <cit.>, with a quantum singular value transformation being a recent innovation <cit.>. These algorithms have found widespread application in quantum machine learning <cit.>. Many problems in machine learning and data processing involve a large number of matrix operations and it is useful to concatenate such operations effectively on a quantum computer. In an ideal scenario, a user of a quantum computer only declares a series of instructions on a set of matrices in a simple language, which is then compiled into a sequence of quantum gates for the quantum computer. An important intermediate step is to provide easy-to-compose techniques for matrix linear algebra.

In this paper, we present a scheme for compiling elementary linear algebraic operations on a quantum computer for a set of complex matrices, given the availability of particular unitaries generated from these matrices. The elementary matrix operations we consider here are matrix addition, multiplication, Kronecker sum, tensor product and Hadamard product, as well as arbitrary real single-matrix functions h(A) that act on the eigenvalues of a matrix A. We refer to this scheme as the Quantum Matrix Algebra Toolbox (QMAT) in the following. Specifically, given a set of matrices {A_j}, we embed these matrices into a set of Hermitian matrices and assume as input a given set of unitary operators generated by the embedding matrices. In addition, we are given a matrix function  f({A_j}) that can be divided into tiers of sub-functions, where each sub-function is one of the QMAT operations.  We are then able to compose the matrix operations and compile the tiers of sub-functions into a sequence of operations for the quantum computer. This results in an approximation of a unitary operator which encodes the computation of the matrix function f({A_j}) in a composable manner. Such a toolbox is of limited use if the desired outcome is a large vector or matrix since in this case obtaining a complete classical representation is often resource intensive <cit.>. Hence we discuss scalar outputs, such as the trace, determinant, and Schatten p-norm of f({A_j}), which often can be obtained efficiently. We  note recent efficient development of logarithmic time classical algorithms for various machine learning tasks <cit.> based on preprocessed classical data. In contrast, the scheme presented here assumes only the availability of particular unitaries and is inherently BQP complete, as it allows for the encoding of arbitrary quantum computations.



§ EMBEDDING


We assume a set of matrices {A_j}, j=1,…,J, with A_j∈ℂ^N× N. 
The matrices are not required to be Hermitian.
The restriction to square matrices is for consistency in the output of the basic matrix operations. Non-square matrices can always be padded by zeros to turn them into square matrices. The results of the QMAT operations for the original non-square matrices can be obtained by removing the corresponding zeros from the results of the extended matrices.
 In the QMAT setting, the non-Hermitian matrices {A_j} are embedded into Hermitian matrices. 
 A Hermitian matrix can be interpreted as a Hamiltonian operator in quantum mechanics and Hamiltonian simulation techniques can be used to perform matrix operations.  
 Given a matrix A∈ℂ^N× N,  we define the following embedding matrices X_1(A), X_2(A), and X_3(A)  as
 
    X_1(A)=R_1⊗ A+R_1^†⊗ A^†,	
       X_2(A)=R_2⊗ A+R_2^†⊗ A^†,	
       X_3(A)=R_3⊗ A+R_3^†⊗ A^†,

where the R_1, R_2 and R_3 are 3×3 matrices, given by

    R_1   =[[ 0 1 0; 0 0 0; 0 0 0 ]], 
    R_2=[[ 0 0 0; 0 0 1; 0 0 0 ]],   
    R_3=[[ 0 0 1; 0 0 0; 0 0 0 ]].

Analogously, we embed the vectors  x∈ℂ^N , y∈ℂ^N
 into the vectors V_1(x), V_2(y), which are defined as 

    V_1(x)=r_1 ⊗x,  
    V_2(y)=r_2⊗y,

where  r_1 and r_2 are 3-dimensional vectors

    r_1   =[[ 1 0 0 ]]^T, r_2=[[ 0 0 1 ]]^T.

It follows that

    V_1(x)^† X_3(A) V_2(y)=x^† Ay.

Hence the value of x^† Ay can be obtained by calculating the inner product of the corresponding embedded matrix and vectors. 

We now define the following three permutation operators, P_1, P_2 and P_3 

    P_1=[[ 0 0 I; I 0 0; 0 I 0 ]], 
    P_2=[[ 0 I 0; I 0 0; 0 0 I ]], 
    P_3=[[ I 0 0; 0 0 I; 0 I 0 ]],

where I denotes the N × N identity matrix and 0 denotes the N × N matrix of zeros. Then the embedding matrices X_1(A), X_2(A) and X_3(A) can be transformed into one another by applying these permutation operators

    P_i X_i(A)P_i^†   =   X_(i+1)mod  3(A),   (i=1,2,3),

which also holds for the exponentiated X_i(A) since P_1, P_2 and P_3 are unitary.
The transformations above can also be achieved in the reverse direction.



§ QUANTUM ALGORITHMS FOR MATRIX OPERATIONS


We now discuss the quantum algorithms for matrix operations. We start with the required data input model and continue with obtaining the desired output.
We then show how one can simulate the unitary operator generated by X_3(f({A_j})) when given access to a set of input unitary operators generated by X_3(A_j). We will describe the QMAT quantum subroutines for the exponentiation of matrix addition, multiplication, Kronecker sum, tensor product and Hadamard product of matrices A_1,A_2∈ℂ^N× N, as well as single matrix functions.  



 §.§ Data input

We first summarize the data input required in this work. 

Given matrices A_j for j=1,…,J, let ‖ X_3(A_j) ‖_maxτ = O(1), where τ is a  time parameter[We assume that τ can be made arbitrarily small.]. Assume access to the unitaries  e^iX_3(A_j)τ. In addition, given arbitrary ancilla qubits assume access to the controlled unitaries e^i |1⟩⟨1|⊗ X_3(A_j)τ.

Here, we use the maximum element norm ‖ A‖_max=max_ij| A_ij|.
If the matrices A_j, and hence X_3(A_j), are sparse and readily stored in a sparse matrix data structure, such unitaries are provided by using well-studied quantum walk techniques <cit.>. These techniques have been continuously improved to a nearly optimal query complexity of γlog(γ/ϵ)/loglog(γ/ϵ) <cit.>, employing the technique of linear combination of unitaries <cit.>. Here, γ=d τ‖ H‖_max, where  d is the sparsity and ϵ is the desired accuracy in operator norm.   The gate complexity of this approach scales as log N + log^5/2(d/ϵ) for each query. An optimal query complexity of  γ+ log(1/ϵ)/loglog(1/ϵ)  has been achieved by using quantum signal processing <cit.>. 
Additionally, if the matrix is provided as quantum density matrix we can use the quantum state exponentiation technique <cit.>.


If we are interested in inner products such as x^† f({ A_j})y, we require the preparation of quantum states corresponding to x and y.  This input assumption is not required if we are interested in the inherent scalar quantities of the output matrices such as the trace or determinant.

Assume routines that prepare the quantum states |V_i(v)⟩=V_i(v)/|v| for the classical vectors v = x,y and i=1,2.

Such data access can be provided if the data are stored in a quantum random access memory (QRAM) as discussed in <cit.>.



 §.§ Output



We start with a description of a method for calculating an inner product ⟨V_1(x)|X_3(f({A_j}))|V_2(y)⟩.  
A method for estimating inner products of the form ⟨V_1(x)|V_2(y)⟩ was presented in <cit.>, which we describe in Appendix <ref>. This procedure can be extended to compute ⟨V_1(x)|X_3(f({A_j}))|V_2(y)⟩ by performing the matrix-vector multiplication X_3(f({A_j}))|V_2(y)⟩, 
analogous to the algorithms presented in <cit.>. If given access to the controlled unitary e^i |1⟩⟨1|⊗ X_3(f({A_j}))t, we can perform the quantum phase estimation algorithm <cit.> using |V_2(y)⟩ as the input state. This phase estimation results in 

    ∑_k=0^N-1γ_k|μ_k(f)⟩|λ̃_k(f)⟩,

where λ̃_k(f) approximates the eigenvalues λ_k(f) ≡λ_k(X_3(f({A_j}))) of the matrix 
X_3(f({A_j})) and |μ_k(f)⟩ are the associated eigenvectors. We have defined  γ_k=⟨μ_k(f)|V_2(y)⟩.
The phase estimation requires a runtime 1/ϵ_λ for an ϵ_λ approximation of the eigenvalues <cit.>. 

Subsequently, we perform a controlled rotation of an ancilla register initialized in |0⟩ conditioned on the eigenvalue register, which acts on 
|λ_k(f)⟩|0⟩ as 

    |λ_k(f)⟩(√(1-c^2 λ_k^2(f))|0⟩+ cλ_k(f)|1⟩).

The constant c is chosen such that cλ_k(f)≤1 for all k. We then uncompute the eigenvalue register and measure the ancilla qubit. A measurement result of |1⟩ leaves the remaining system in the desired output proportional to

    ∑_k=0^N-1γ_k λ_k(f) |μ_k(f)⟩=X_3(f({A_j}))|V_2(y)⟩.

The success probability of the measurement is given by ∑_k γ^2_k c^2 λ_k^2.

We continue with the QMAT subroutines.
We present the subroutines for convenience using the unitaries without the control register, see Input Assumption
<ref>.  We note that the subroutines can be made controlled by simply using the corresponding controlled unitaries. 












 §.§ Addition


The addition of matrices via Trotter-based methods is well-known and we discuss it here for completeness in the context of the embedding matrices.
Given access to e^iX_3(A_1)τ and e^iX_3(A_2)τ,  we are able to approximate e^i(X_3(A_1)+X_3(A_2))t with bounded error by using the Trotter formula, see Subroutine <ref>. Hence given access to e^iX_3(A_1)t/n and e^iX_3(A_2)t/n, a constant ϵ_1-error approximation of e^i(X_3(A_1)+X_3(A_2))t can be achieved in n=O(t^2/ϵ_1) steps.
The spectral norm ‖·‖ is used for quantifying the error.
It is worth mentioning that 
the constant error can be reduced with the higher-order Suzuki-Trotter formula <cit.>.
For example, given further access to e^iX_3(A_i)t/2n, the same error can be achieved in n=O(t^3/2/√(ϵ_1)) steps using longer sequences. 
Note that our Input Assumption <ref> allows for arbitrary-τ unitaries and thus for all higher-order refinements.


  









 §.§ Multiplication


One way of performing matrix multiplication applied to a quantum state is via two consecutive phase estimations and controlled rotations. Here, we present a version in the QMAT setting by using the commutator Lie formula.
We show how to carry out a matrix multiplication by using the commutator 

    [X_1(A_1), X_2(A_2)]    = [[    0    0    M;    0    0    0; -M^†    0    0 ]],

where M=A_1A_2. Note that the result of the commutator shown in the above equation is not a Hermitian matrix, so we  construct a Hermitian matrix via an imaginary factor

    i[X_1(A_1), X_2(A_2) ]=X_3(iM).

It follows that

    e^iX_3(iM)=e^-[X_1(A_1), X_2(A_2) ]=e^[iX_1(A_1), iX_2(A_2) ].

However, our goal is to embed the matrix M  instead of iM. Fortunately, the term  e^iX_3(M) can be obtained from e^iX_3(iM) by using an unitary operator

    U_1=[[ √(-i)I      0      0;      0      I      0;      0      0  √(i)I ]],

it follows that,

    e^iX_3(M)   =U_1e^iX_3(iM)U_1^†.


Based on  Eq. (<ref>), we present a quantum subroutine to find the exponentiation of the product of matrices A_1 and A_2 in subroutine <ref>.  For simulating the commutator, 
 we require the sequence

    l̃(x_1,x_2)   :=    e^x_1e^x_2e^-x_1e^-x_2 e^-x_1e^-x_2e^x_1e^x_2

where x_1 and x_2 are placeholders for arbitrary matrices. 
As we show in Appendix <ref>, the error of l̃(X_1(A_1)t/n,X_2(A_2)t/n)^n^2/2t for simulating the commutator is bounded by  ϵ_2=O(t^3/n^2)  which requires n'=O(t^2/ϵ_2) applications of e^iX_1(A_1)t/n and e^iX_2(A_2)t/n.  These unitaries are obtained from Input Assumption <ref> by using the permutation relations in Eq. (<ref>).









 §.§ Kronecker sum and tensor product 


Building on the above results, we can also find an approximation of the exponentiation of the Kronecker sum and tensor product of matrices A_1 and A_2. 
Simulating Kronecker sums on quantum computers is the topic of the initial works on quantum simulation <cit.>, where the Hamiltonian is a sum of local Hamiltonians. For completeness, we present the simulation of Kronecker sums here in the context of the embedding matrices.
The Kronecker sum of these two matrices denoted by ⊕ is defined as

    A_1⊕ A_2:=A_1⊗ I +I⊗ A_2,

where I is the identity matrix.  Since A_1⊗ I  and I⊗ A_2 commute, one can obtain the exponentiation of the Kronecker sum by 

    e^A_1⊕ A_2= (e^A_1⊗ I) (I⊗ e^ A_2).

Now we consider the exponentiation of A_1⊕  A_2  with the embedding matrices by using subroutine <ref>. The exponentiation of embedded A_1⊕ A_2 is 

    e^iX_3(A_1⊕ A_2)t=e^i(X_3(A_1⊗ I) +X_3(I⊗ A_2))t.

Using the definition of X_3(A), we have

    X_3(A_1⊗ I)   =   
    R_3⊗ A_1⊗ I+R_3^†⊗ A_1^†⊗ I,
    
    X_3(I⊗ A_2)   =    R_3⊗ I⊗ A_2+R_3^†⊗ I⊗ A_2^†.

Considering the three registers,  note that X_3(A_1) acts on the first and second register and X_3(A_2) acts on first and third register.   
Thus, we can approximate e^iX_3(A_1⊕ A_2)t
using subroutine <ref> with the inputs e^iX_3(A_1)τ and e^iX_3(A_2)τ on the appropriate registers.


We can also derive the exponentiation of the tensor product of matrices A_1 and A_2, using the permutation relations in Eq. (<ref>) and subroutine <ref>. Since

    A_1⊗ A_2=(A_1⊗ I )(I⊗ A_2),

the Hamiltonian simulation of A_1⊗ A_2 can be performed with the embedding e^iX_3(A_1⊗ A_2)t using subroutine <ref> with the inputs e^± iX_3(A_1⊗ I)τ and e^± iX_3(I ⊗ A_2)τ, noting that X_1(A) and X_2(A) can be obtained from X_3(A) via the permutations defined in Eq. (<ref>).








 §.§ Hadamard product


We denote the Hadamard product of matrices A_1 and A_2  as A_1∘ A_2, which is a matrix given by

    (A_1∘ A_2)_ij= (A_1)_ij (A_2)_ij.

We now discuss a method to obtain such a Hadamard product in the QMAT framework. Define the following (non-Hermitian) matrix 

    S=∑_i=0^N-1|i⟩⟨i|⊗|0̅⟩⟨i|,

where |0̅⟩:=|0… 0⟩. Since the matrix S is sparse, there exists an efficient quantum algorithm to simulate the embedded sparse matrix X_3(S) <cit.>.
We can obtain the Hadamard product of A_1 and A_2 from the tensor product in the following way

    S(A_1⊗ A_2)S^†
    =(A_1 ∘ A_2)⊗|0̅⟩⟨0̅|.

The size of the resulting matrix is extended by the operator |0̅⟩⟨0̅|. In the embedding matrix form, we note the relation

    e^iX_3(S(A_1⊗ A_2)S^†)t 
    =e^iX_3(A_1 ∘  A_2)⊗|0̅⟩⟨0̅| t
           = e^iX_3(A_1 ∘  A_2)t⊗|0̅⟩⟨0̅|+ I⊗(I-|0̅⟩⟨0̅|),

where the last equality follows from the series expansion of the exponential. 
Thus, we approximate  e^iX_3(A_1 ∘ A_2)t  by using e^iX_3(S(A_1⊗ A_2)S^†)t with an ancillary register in state |0̅⟩. 
The operator e^iX_3(S(A_1⊗ A_2)S^†)t can be constructed by  
combining the methods for the matrix tensor product and the multiplication subroutine <ref> with access to e^± iX_3(S)τ, e^± iX_3(S^†)τ, e^± iX_3(A_1)τ and e^± iX_3(A_2)τ. 







 §.§ Real operator functions


We define an arbitrary real operator function h(A) of a Hermitian matrix A as
 
    h(A)=∑_λh(λ)|μ_λ⟩⟨μ_λ|,

 where λ and |μ_λ⟩ are the eigenvalues and the corresponding eigenvectors respectively.
Such a function can be applied in the QMAT framework using results presented in Ref. <cit.>, where a procedure was given to construct a quantum circuit that performs W'=∑_λe^ih(λ)|u_λ⟩⟨u_λ| for a real odd and periodic function h(λ):(-π,π]→ (-π,π], given a controlled unitary W=|0⟩⟨0|⊗ I + |1⟩⟨1|⊗∑_λe^iλ|u_λ⟩⟨u_λ|. One can apply the formalism to even functions h(λ) by simply dividing the even function into two odd functions  h(λ)/λ and λ, then multiplying them with using the functionality of  subroutine <ref>. More generally, if h(λ) is neither even nor odd, we first divide it into an even function 1/2(h(λ)+h(-λ)) and an odd function 1/2(h(λ)-h(-λ)), then recombine them using the functionality of subroutine <ref>. 

Now we consider an odd  function h(A) of a Hermitian matrix A in the embedding formula. According to the following equation

    [[ 0 0 A; 0 0 0; A 0 0 ]]
    [[  |u_λ⟩;      0; ±|u_λ⟩ ]]
    =±λ[[  |u_λ⟩;      0; ±|u_λ⟩ ]],

we see that the eigenvalues of X_3(h(A)) are ± h(λ) whereas the eigenvalues of h(X_3(A)) are h(±λ). Thus we have X_3(h(A))=h(X_3(A)) since ± h(λ)=h(±λ) for an odd function. 

As such we have shown that given an input of the form e^iX_3(A)τ, it is possible to obtain e^iX_3(h(A))t for general real operator functions. The outcome of this subroutine is in the desired matrix exponent form, and hence can be concatenated with all other subroutines discussed in this paper.

To summarize, these subroutines allow performing a mix of operations of matrices in a concatenated fashion. Given a set of  unitaries generated by {X_3(A_j)}, j=1, …, J, and a large class of functions f({A_j}) of these matrices, 
we are able to construct the unitary operator e^iX_3(f({A_j}))t. We note that each procedure produces an output satisfying the Input Assumption <ref>, which allows for concatenation of the procedures. The value of ⟨x| f({A_j})|y⟩ can be estimated efficiently, as was shown in Sections <ref>. Other scalar quantities can also be estimated, as will be shown in the next section.









 §.§ Estimating norms, traces and determinants

Let A=f({A_j}), we now show a method to approximate the Schatten p-norm of A using the QMAT embedding matrices. The Schatten p-norm of matrix A is defined as ‖ A‖_p=(∑_k=1^Nσ_k^p)^1/p for  σ_1 ≥σ_2...≥σ_N ≥ 0 the singular values of A. Since the eigenvalues of X_3(A) are {±σ_k}, we can estimate the Schatten p-norm by sampling the absolute value of the p-th power of the eigenvalues of X_3(A) then calculating the p-th root of the result.  It can be achieved by performing e^iX_3(A) on a maximally mixed state then using phase estimation and measurement to get the eigenvalues. A discussion on the Schatten p-norm via the DQC1 protocol can be found in <cit.>.

Similarly, the trace of a Hermitian matrix can be estimated by sampling the eigenvalues. Notice that, in the embedding formula, the trace of X_3(A) is zero, thus we cannot estimate the trace of A from the eigenvalues of X_3(A) directly. If A is non-Hermitian, we construct the Hermitian matrices A+A^† and i(A-A^†). Then the real part of (A) is equal to 1/2(A+A^†) and the imaginary part is equal to i/2(A-A^†).
Now we show the method to estimate (A) where A is Hermitian. 
The eigenvalues of A are {λ_i}, such that the eigenvalues of X_3(A) are {±λ_i}. Construct a matrix A'=A+cI, with c such that all eigenvalues of A are shifted to positive values. 
The sampling method proceeds as follows: apply e^iX_3(A') on a maximally mixed state, then perform phase estimation and measurement which yield one of the eigenvalues of X_3(A'), say ±λ'_i = ± (λ_i + c). Using the absolute value of the measurement outcome, the corresponding eigenvalue of A can be extracted by λ_i=|λ'_i| -c, which is then used to evaluate the desired trace. This procedure works, because by operating on the maximally mixed state we are effectively choosing an eigenvector uniformly at random, and hence sampling the eigenvalues uniformly at random. An alternate approach to trace estimation would be to sample diagonal elements of the matrix in a basis chosen uniformly at random from a maximal unbiased set. In the classical problem of stochastic trace estimation, such a procedure has been shown to yield the same expectation value with greatly reduced variance when compared to sampling in any fixed basis <cit.>.

 If A is a Hermitian positive definite matrix, one can estimate the determinant [A] as well. This can be achieved by leveraging the equality log([A])=(log(A)), and hence [A] = exp((log(A))). Hence by combining a trace estimation procedure with a procedure to generate evolution under log(A), and post processing the result, one obtains a procedure to estimate [A] <cit.>. The error analysis of the estimation methods is given in Appendix <ref>. It shows that the Schatten-p norms, the trace and determinant of matrix A can be estimated with small relative error with high success probability. 





 §.§ Example


Given a matrix function f({A_j}), we express the function f(·) as an abstract syntax tree, composed of basic single-matrix and two-matrix operations. Consider an example in which we have six matrices, A_j for j=1,…,6, and a matrix function f({A_j})=((A_1A_2)⊗ h(A_3))∘ (A_6 ⊕ (A_4+A_5)). Figure <ref> shows the corresponding tree. The unitary  e^iX_3(f({A_j})) is generated by concatenating the different subroutines using the embedding matrices.
 
  







§ DISCUSSION AND CONCLUSION


The Quantum Matrix Algebra Toolbox that we have presented here gives a consistent and composable set of procedures for manipulating matrices, represented as unitary evolution of a quantum system, and for estimating scalar valued functions of these matrices. This toolbox allows for the estimation of scalar properties of matrix functions corresponding to arbitrary abstract syntax trees containing real eigenvalue functions of a single matrix as well as the most common two-matrix products and sums. Matrix embeddings allow us to simulate complex unitaries with sequences of simple unitaries, in the spirit of the usual Trotter-based Hamiltonian simulation techniques. For a set of matrices {A_j}, the methods can obtain ⟨x| f({A_j}) |y⟩  and the trace of f({A_j}) with bounded errors, and also individual elements of f({A_j}) if we choose ⟨x| and |y⟩ as the row and column basis vectors.  As the underlying sparse Hamiltonian techniques are logarithmic in the dimension of the matrices, we may achieve the estimation of such quantities much faster than classical algorithms which often scale polynomial in the dimension.

In our approach, we have built up the full set of matrix operations by using Trotter formulae to approximate evolution under the sum or commutator of matrices. Alternative approaches exist for calculating matrix functions with exponentially improved error scaling. In Ref. <cit.> it was shown that one can compute addition and multiplication of two matrices if the block-encodings of these matrices are given. The techniques in the present work can potentially be combined with the block-encoding framework to construct alternative compiling schemes with improved error scaling and hence improved efficiency.

In conclusion, the techniques we present here represent a systematic method for the compilation of user-specified instruction sets of desired matrix operations into gate sequences for actual quantum computers. We believe that these techniques will prove useful in constructing further quantum algorithms based on fast matrix algebra.




This material is based on research supported by the Singapore National Research Foundation under grant NRF-NRFF2013-01 and NRF2017-NRF-ANR004, and by the US Airforce Office of Scientific Research under grant FA2386-18-1-4003. The authors also acknowledge support from Singapore's Ministry of Education.





§ COMMUTATOR SIMULATION




Given access to e^ H_1t/n and e^H_2t/n  for two  matrices H_1, H_2 ∈ℂ^N× N for small t, we can approximate e^[H_1,H_2]t with bounded error by using the second order of the Baker-Campbell-Hausdorff formula <cit.> as follows

    e^[H_1,H_2]t+O(t^2/n^2)
    = (e^H_1t/ne^H_2 t/ne^-H_1 t/ne^-H_2 t/n)^n^2/t.
 
In order to reduce the error, we rearrange the order of  every term in the right side of  Eq. (<ref>). 
Now we define   
 
    l(H_1t/n,H_2t/n):=e^H_1t/ne^H_2t/ne^-H_1t/ne^-H_2 t/n,
 
 and we have

    l(-H_1t/n,-H_2t/n)   =   e^-H_1t/ne^-H_2t/ne^H_1t/ne^H_2t/n.

 Then, combining the above two equations, we obtain

    l̃(H_1 t/n, H_2 t/n)   :=   l(H_1 t/n, H_2 t/n)l(-H_1 t/n, -H_2 t/n)
       =    e^2[H_1,H_2]t^2/n^2+O((t/n)^4).

Hence, the term  e^[H_1, H_2]t can be approximated with bounded error as
  
    l̃(H_1t/n,H_2t/n)^n^2/2t=e^[H_1,H_2]t+O(t^3/n^2).

Let n'=n^2/2t, 
we see that e^[H_1,H_2]t can be approximated to a constant ϵ_2-error by using n'=O(t^2/ϵ_2) copies of e^H_1t/n and e^H_2t/n.




§ QUANTUM INNER PRODUCT ESTIMATION



We now review a method for obtaining ⟨x|y⟩ for given vectors |x⟩,|y⟩∈ℂ^N. Assume we are give a state 

    |φ⟩=1/√(2)(|0⟩|x⟩+|1⟩|y⟩).

Applying a Hadamard operator on the first qubit, this yields

    |φ⟩=1/2(|0⟩(|x⟩+|y⟩)+|1⟩(|x⟩-|y⟩)).

Measuring the first qubit in computational basis, the probability to obtain |0⟩ is given by

    p=1/2(1+ Re(⟨x|y⟩)).

By repeating this procedure a constant number of times, we can estimate the real part of ⟨x|y⟩ to fixed precision. Likewise, the imaginary part of ⟨x|y⟩ can be obtained by applying a phase rotation on the first register in Eq. (<ref>), to produce

    |φ⟩=1/√(2)(|0⟩|x⟩-i|1⟩|y⟩),
 
followed by applying a Hadamard matrix and measuring the first qubit in Pauli Z basis, the probability of obtaining |0⟩ is given by 

    p=1/2(1+ Im(⟨x|y⟩)).

This can be repeated a constant number of times to get an estimate of the imaginary part of ⟨x|y⟩ to fixed precision. In this way, we can approximate the inner product of vector x and y.
 
As the variance of binomial distribution specified by the probability p and the number of trials m 
is given by mp(1-p), the variance of the estimate for p is equal to p(1-p)/m. The variance of the real part is 4 times the variance of the estimate for p. Thus the error of the real part of the inner product is then given by 

    (1/m(1- Re(⟨x|y⟩)^2) )^1/2,

 which is then bounded by 1/√(m), and similarly for the imaginary part.



§ ERROR ANALYSIS OF THE SCHATTEN-P NORM ESTIMATION



We would like to estimate the Schatten-p norm of a matrix A by uniform sampling its singular values. 
Assume a relative error of the singular values σ_i from quantum phase estimation be ϵ_σ>0.
Concretely, we are interested in estimating ‖ A ‖ ^p_p to fixed accuracy, independent of the dimension of A, with high probability.

Formally, we have random variables σ̃_i with expectation value E[σ̃_i] = σ_i and variance Var[σ̃_i] = ϵ_σ^2 σ_i^2.
In addition, we sample the index i uniformly, thus we define a random variable X̃ to be the random variable that has as outcomes each σ̃_i^p with probability 1/N. Let K_p the Lipshitz constant for the function f_p(x) = x^p in the interval [0,σ_max]. Note that

    |σ_i^p - σ̃_i^p |≤ K_p |σ_i - σ̃_i |,
 
and also the special case σ̃_i^2p≤ K_p^2  σ̃_i^2.
Let Ỹ be the random variable N/T∑_ j = 1^T X̃_j, where X̃_j denote independent instances of the random variable X̃.
This random variable (Ỹ) has the expectation value, 

    E[Ỹ]    =   N/T E[∑_j=1^T X̃_j]
       =   N E[X̃] = N ∑_j=1^N 1/N E[σ̃_j^p] 
       =   ‖ A‖_p^p+ ∑_j=1^N E[σ̃_j^p]-‖ A‖_p^p
       ≤   ‖ A‖_p^p+ ∑_j=1^N E[|σ̃_j^p - σ^p_j| ]  
       ≤   ‖ A‖_p^p+  K_p ∑_j=1^N  E[|σ̃_j -σ_j| ] 
       ≤   ‖ A‖_p^p+ K_p ∑_j=1^N  √(E[|σ̃_j -σ_j|^2 ])
       =    ‖ A‖_p^p+ K_p ∑_j=1^N  √(Var[σ̃_j])
       ≤   ‖ A‖_p^p+ K_p ∑_j=1^N  ϵ_σσ_i 
       =   ‖ A‖_p^p+K_p ϵ_σ‖ A‖_1.

Here, we have used Jensen's inequality twice and Lipshitz continuity.
For the quadratic term, we have, 

    E[Ỹ^2]   =   N^2/T^2 E[ ∑_j=1^T X̃_j^2  ]  
       =   N^2/T E[ X̃^2  ] 
       =   N/T∑_j=1^N E[ σ̃_j^2p ] 
       ≤   N K_p^2/T∑_j=1^N  E[ σ̃_j^2 ] 
       =   N K_p^2 (1+ϵ_σ^2)/T∑_j=1^N σ_j^2 
       =   N K_p^2 (1+ϵ_σ^2)/T‖ A ‖_2^2


According to Chebyshev's inequality, we obtain

    P(|Ỹ- E[Ỹ]|≥ε‖ A‖_p^p )    ≤   Var(Ỹ)/ε^2 ‖ A‖_p^2p≤E(Ỹ^2)/ε^2 ‖ A‖_p^2p
       ≤   N/T ε^2K_p^2 (1+ϵ_σ^2)/‖ A‖_p^2p‖ A ‖_2^2.

This failure probability is employed to bound the final error. 
Note that 

    |Ỹ - ‖ A‖_p^p|≤|Ỹ - E[Ỹ]  | + K_p ϵ_σ‖ A‖_1.

Thus, with probability 1-P_fail the error is bounded as

    |Ỹ - ‖ A‖_p^p|≤ε‖ A‖_p^p + K_p ϵ_σ‖ A‖_1,

where

    P_ fail≤N/T ε^2K_p^2 (1+ϵ_σ^2)/‖ A‖_p^2p‖ A ‖_2^2.

Assume that σ_i=Θ(1) for all i. Then the above expression can be simplified as

    P_fail≤C K_p^2 (1+ϵ_σ^2)  /T ε^2.
 
where C is a constant. Let C K_p^2 (1+ϵ_σ^2)  /T ε^2=a, then
the number of samples required is T = ⌈C K_p^2 (1+ϵ_σ^2)  /ε^2 a ⌉. We can take a=0.01, for example, to achieve a 99% confidence.
In summary, the outcome of the estimation procedure is bounded by

    ‖ A‖_p(1±ε± K_p ϵ_σ‖ A‖_1/‖ A‖_p^p)^1/p,
 
with high probability. Since x^1/p is a concave function, and since we have already assumed that the singular values satisfy σ_i=Θ(1), this implies that the relative error in the estimation of ‖ A‖_p is bounded by 

    ε_R    ≤    |ε|^1/p + (K_p |ϵ_σ| ‖ A‖_1/‖ A‖_p^p)^1/p
       =    |ε|^1/p + κ_p |ϵ_σ|^1/p,
 
where κ_p is some constant that depends only on p and the allowed (constant) range of the singular values of A. Thus the relative error does not depend on the dimensions of A, provided that the assumption on the relative size of the non-zero singular values of A holds. 

The error analysis of trace estimation is similar to the Schatten 1-norm. Assume the relative error from phase estimation of the eigenvalues be ε_λ. The sampling outcome is then bounded by (A)(1±ε±ϵ_λ) with a number of samples T= C(1+ϵ_λ^2)/ε^2 a, which is obtained from Eq. (<ref>) taking K_p=1 and p=1. 

Finally, for the determinant, the error analysis is as follows. Let ϵ_log(λ) be the relative error of the eigenvalues of log(A) from the phase estimation step. Since (A)=e^(log(A)), the determinant of a matrix A can be estimated with a number of samples T=C(1+ϵ_log(λ)^2)/ε^2 a. The error is then given by

    | e^(log(A))-e^(log(A))|   ≤    K_e |(log(A))-(log(A))|
        =     K_e |(log(A))(ϵ_log(λ)+ϵ)|,

with high probability, where K_e is the Lipshitz constant for the function f(x)=e^x. The Lipshitz constant is 1 if x ≤ 0, which corresponds to the case when the spectral norm is ‖ A‖≤ 1.

apsrev


