A Low-Complexity Quantum Simulation Framework for Toeplitz-Structured Matrix and Its Application in Signal Processing

Toeplitz matrix reconstruction algorithms (TMRAs) are one of the central subroutines in array processing for wireless communication applications. The classical TMRAs have shown excellent accuracy in the spectral estimation for both uncorrelated and coherence sources in the recent era. However, TMRAs incorporate the classical eigenvalue decomposition technique for estimating the eigenvalues of the Toeplitz-structured covariance matrices that demand very high computational complexity for large arrays. We demonstrate a low-complexity quantum simulation framework exploiting the structured Hamiltonian of Toeplitz and circulant variants. In this framework, we consider two approaches for the estimation of the eigenvalue spectrum of a given Toeplitz-structured matrix: first, an analytical framework with Jordan form-based sparse-decomposition of a dense-Toeplitz matrix, and second, an approximation method for the conversion of a Toeplitz matrix into a circulant matrix embedding quantum subroutines. We have also compared the efficacy of the proposed algorithms with standard Hamiltonian simulation and quantum phase estimation techniques for different quantum time resolutions and gate complexities. We show quantum gate-complexity analysis for our proposed algorithms. Considering the large dimensions of the Toeplitz matrix, we have employed random matrix theory in deriving the error bounds for the estimated eigenvalues. The numerical results are obtained partly in a classical computer and in an IBM quantum simulator.

INDEX TERMS Quantum algorithms, quantum communication, quantum signal processing, quantum simulation, structured Hamiltonian, Toeplitz and circulant systems. NOMENCLATURE |ψ ∈ C n×1 Quantum state vector in column form. C, C n , C m×n Set of complex scalar, vector, and matrices. R, R n , R m×n Set of real scalar, vector, and matrices.
x , A 2-norm of vector x and matrix A, respectively. |A| Determinant of matrix A.

Diag(x)
Diagonal matrix with the vector x in diagonal.

Pr(x)
Probability of an event x.

dim(A)
Dimension of the matrix A.
x, x, X Scalar, vector, and matrix quantity. A † , a † (Complex conjugate) transpose of the (complex) matrix A, and (complex) vector a respectively. Z

Set of integers. E[x]
Expectation of a random variable x.

I. INTRODUCTION
Toeplitz matrices are ubiquitous in signal processing and communication systems. Structured matrices such as Toeplitz and circulant are used in several applications like the direction-of-arrival (DoA) estimation [1], low-rank matrix recovery [2], quadratic optimization [3], image reconstruction [4], and compressive sensing [5] for its complexity advantage over standard systems. Further, the covariance Engineering uantum Transactions on IEEE Laskar et al.: LOW-COMPLEXITY QUANTUM SIMULATION FRAMEWORK matrix obtained from the sample of a wide sense stationaryrandom process holds the Toeplitz structure, which gives practical advantages for performing certain signal processing tasks [6]. They include covariance matrix estimation for sparse array [7], and compressive covariance sampling for spectrum estimation [8]. Several properties and asymptotic behavior of Toeplitz matrices are described by Szego's theorem [9]. General practice is to transform the Toeplitz matrix T N ∈ C N×N into a circulant matrix C N ∈ C N×N using suitable procedures to get O(N log N) time-complexity for the computation of eigenvalues. The asymptotic equivalence and convergence of the circulant matrix conversion process are shown in [10] and [11].

A. TOEPLITZ MATRIX RECONSTRUCTION APPROACHES IN RECENT WIRELESS COMMUNICATION APPLICATIONS
The Toeplitz matrix reconstruction algorithm (TMRA) is widely used in the recent era for many wireless communication systems. For wireless channel estimation in a multitap intersymbol-interference scenario, the wireless channel is often modeled as a banded Toeplitz-matrix [12], [13], [14]. In recent times, for super-resolution DoA estimation problems in array processing, Toeplitz matrix-based low-complexity TMRA is proposed [1], [15], [16], [17], [18], [19]. We are interested in finding the spectrum of a Toeplitz-structured system using the quantum Hamiltonian simulation technique to augment quantum advantage for certain applications. Estimating the spectrum through the quantum phase estimation (QPE) procedure requires that the Hamiltonian system needs to be a square and a Hermitian system. We explore a Toeplitz system that arises in the DoA estimation in array processing.

B. ADVANCES IN QUANTUM SIMULATION-BASED SPECTRUM ESTIMATION
Quantum simulation of physical systems in form of operators and observables has become one of the most promising areas in quantum computation and quantum signal processing. Specifically, the quantum eigenvalue estimation (QEE) problem stands as one of the central building blocks in quantum literature [20], [21], [22], [23]. Quantum simulation-based subroutines such as Hamiltonian simulation, QPE, and inverse quantum Fourier transform (IQFT) have been used as building blocks for designing quantum circuits and systems [24], [25], [26]. Embedding quantum simulation with a QPE technique, one can solve a large-scale linear system of equations as proposed in [27] for Hermitian and sparse matrices with exponential speed-up as compared to the classical algorithm.
In the recent era, solving certain linear systems in a quantum framework possess quantum complexity advantages [27]. While solving the linear systems on a quantum computer, we need to perform a Hamiltonian simulation that takes the system (satisfying Hermitian property) as input and prepares a unitary operator [28], [29]. Often, the computational complexity advantage considers that the Hamiltonian system is sparse in nature. For a dense matrix with embedded matrix structures such as Toeplitz, Hankel, and circulant, the Hamiltonian simulation is not addressed well.
Our objective is to simulate a Toeplitz system T N through quantum simulation efficiently. First, we take a general dense Toeplitz matrix T N for finding its spectrum λ(T N ) via an approximate linear sum of Jordan forms represented with similar Hermitian sparse matrices. The sparsity augmented through the linear approximation is expected to provide a quantum complexity advantage in the quantum simulation. The Hamiltonian simulation process helps to prepare the unitary operator e −iT N τ , for the given Hamiltonian T N , and time τ . The simulation is practically performed with different methods such as product formulas [30], truncated Taylor series [31], qubitization [29], and quantum walk [28]. Here, our objective is to design a quantum algorithm to approximate a unitary operator U ≈ e −iT N τ of the Toeplitz-structured Hamiltonian for efficient eigenvalue spectrum estimation of the Toeplitz matrix T N .
Second, we show another approximation method, which is the conversion of a Toeplitz matrix to a circulant matrix C N for the estimation of spectrum λ(C N ). Integrating Szego's theorem [9] followed by IQFT, we have proposed a quantum-assisted method for spectrum estimation of the circulant system C N . We expect that the structured Hamiltonian simulation will have many practical applications due to its complexity advantages.
For a dense and non-Hermitian Toeplitz matrix, we show two methods: first, a Jordan decomposition-based sparse-Hamiltonian simulation framework for efficient quantum gate implementation of the structured matrix, and second, an approximation method for finding the spectrum of a circulant system with quantum subroutines. We show that the Toeplitzstructured systems provide a low-complexity implementation for the Hamiltonian simulation via sparse-Jordan decomposition of the underlying Toeplitz system while preparing the unitary operator for the phase estimation. Further, we have proposed a quantum spectrum estimation algorithm for circulant systems using quantum Fourier transform (QFT), which surpasses the requirement of a standard Hamiltonian simulation.

C. MOTIVATION OF THE PROPOSED WORK
From the literature, we have seen that the quantum Hamiltonian simulation is one of the central subroutines for many quantum algorithms including the quantum linear system solver [27]. However, the quantum speed-up is mostly shown for the sparse and Hermitian matrices. In engineering applications, such as signal processing and communications, the system matrices are often structured, dense, and sometimes non-Hermitian. In this work, we have considered a dense Toeplitz matrix as a Hamiltonian and look into its sparse decomposition via a proposed quantum framework. We seek to get a quantum gate-complexity advantage for the Hamiltonian simulation and the spectrum estimation of a Toeplitz-structured matrix. 5500423 VOLUME 4, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.

Transactions on IEEE
We have considered the Hamiltonian simulation of a Toeplitz-structured matrix for a low-complex quantum simulation. It can be both Hermitian and non-Hermitian. The non-Hermitian Toeplitz matrix has extensive usage in literature, such as the study of non-Hermitian Hamiltonian system [32], the generalized minimal residual method for numerical analysis of linear systems [33], and eigenvector localization [34]. Recently, the non-Hermitian doubly Toeplitz matrix structure has been studied in a convolution neural network to improve its performance [35]. Here, we consider a problem, namely TMRA from a signal processing and communication application point of view. It is often seen that the Toeplitz matrix is Hermitian for the TMRA problem such as the DoA estimation [36]. We hope, our work on the quantum simulation of the Toeplitz matrix will have applications for both Hermitian and non-Hermitian systems in near future to provide quantum complexity advantages.

D. CONTRIBUTIONS
With the above background, our contributions to this work are summarized as follows. 1) We have proposed a quantum framework for efficient Hamiltonian simulation considering a dense, Hermitian/non-Hermitian and structured Hamiltonian, which has several technical novelties in its fold. The proposed algorithm augments computational advantage through a Jordan-decomposition-based sparse representation and a modified Hamiltonian simulation. It further opens an area to simulate dense and structured matrices in the quantum framework. 2) We have proposed an algorithm for the spectrum estimation of a Toeplitz-structured Hamiltonian, via circulant approximation. The algorithm shows an improved complexity advantage as compared to the standard Hamiltonian simulation. The proposed one does not require the QPE method and directly computes the spectrum with the QFT circuit and measurement operation, thereby reducing the complexity. 3) We have analyzed the bound in error for the spectrum estimation via circulant approximation of a large Random matrix. In this context, we have developed several propositions and lemmas in the error bound and complexity computation, which would help in the application of our proposed framework in quantum algorithms design where the spectrum of eigenvalues is required.

E. MATHEMATICAL NOTATIONS
Some mathematical notations and symbols used in this article are described in the nomenclature.

F. ORGANIZATION OF THE ARTICLE
The rest of this article is organized as follows. Section II provides the description of Toeplitz-matrices and their representation using quantum gates. In Section III, proposed algorithms for quantum Hamiltonian simulation of Toeplitzstructured matrices are shown. The errors in the spectrum estimation for Toeplitz and circulant systems are addressed in Section IV. The computational complexity for the proposed algorithms is discussed in Section V. We discuss the application of the proposed quantum framework for array signal processing in Section VI. The numerical results for the proposed framework and application are shown in Section VII Data availability for the implementation of the algorithm on IBM quantum qiskit is discussed in Section VIII. Finally, Section IX concludes this article.

II. STRUCTURED HAMILTONIAN FRAMEWORK A. QUANTUM FORMALISM AND STANDARD HAMILTONIAN SIMULATIONS
A brief introduction to quantum states and operators in quantum computing theory is given here, considering the Dirac notations for the state vectors.

1) QUANTUM STATE VECTOR
Quantum states are complex vectors |ψ ∈ C N×1 defined on the Hilbert space H ψ ∈ C N×N . The quantum state is a superposition of the basis states written as where the coefficients β i ∈ C for i ∈ [N], and the probability amplitudes |β i | 2 follows that N i=1 |β i | 2 = ψ, ψ = 1. A qubit is a complex vector in C 2 , given by where the computational basis vectors (also called standard basis, or z basis) are denoted by |0 = [1 0] † , and |1 = [0 1] † , and |β 1 | 2 + |β 2 | 2 = 1.

2) QUANTUM OPERATOR AND OBSERVABLE
The quantum operators are related to the evolution of quantum state vectors in the Hilbert space. The physical observables are Hermitian matrices whose eigenvalues are real and nonnegative. The time-varying dynamics of the quantum state vectors in relation to the observables is described by Schrodinger's equation [37]. A Hamiltonian H is a system matrix (often a Hermitian operator in physical systems) whose time evolution for time τ is given by a unitary operator as that satisfies U † U = UU † = I. However, a quantum measurement operator projects a quantum state to measurement basis vectors to find the probability of the state to a particular basis, which is often a nonunitary operator. Hamiltonian simulation is the algorithm to prepare an approximate unitary operatorŨ from a given physical operator (or observable) H for the evolution time τ and precision as In recent literature, there are some standard Hamiltonian simulation algorithms such as Trotter-Suzuki approximation [38], [39], quantum signal processing [40], quantum walk [41], and Taylor series approximation [31].

3) COMPLETENESS RELATION
Among various possible basis vectors, the orthonormal bases are important for a quantum state to be in superposition given by as e i , e j = δ i j (satisfying the orthogonality relation).
The completeness relation for a set of basis vectors {|i } for i ∈ [N], which can represent a quantum state as shown in (5) must satisfy the following: The eigenvectors of a physical operator (or an observable) provide orthonormal basis states that are complete [42, Ch.

5.3].
A Note on Complexity of the Classical Algorithms for Toeplitz-Based Eigenvalue Spectrum Estimation: The classical eigenvalue spectrum estimation method for the Toeplitz matrix has a complexity of approximatelyÕ(N 3 ) for a fullrank matrix [43],Õ(NP 2 ) andÕ(P 3 ) for rank-deficient Toeplitz matrices (with rank P < N) [44]. In our case, we have considered the eigenvalue estimation of a largedimensional Toeplitz matrix for which the classical algorithms incur significantly very large computational complexity.
Here, we have considered a Toeplitz matrix T N ∈ C N×N as a Hamiltonian operator for estimation of its eigenvaluesspectrum λ(T N ) using a quantum formalism with lowcomputational gate-complexity. Two different approaches for estimating λ(T N ) using the quantum framework are demonstrated. The first one is an analytical framework for a generalized Toeplitz system with a proposed Jordan form-based sparse-Hamiltonian simulation, and the second one is an approximation method for Toeplitz to circulant matrix conversion with a proposed circulant-QFT framework.

B. TOEPLITZ-STRUCTURED HAMILTONIAN SIMULATION
We take a Toeplitz matrix T N with entries denoted by t[h, k] = t[h − k] for h, k ∈ {0, 1, . . . , N − 1} given as follows: The Hamiltonian simulation of matrix T N needs to be performed to prepare a unitary operator for further processing in the QPE technique to find the spectrum. However, the direct implementation of matrix T N requires it to be in symmetric or Hermitian form, i.e., T † N = T N . In this work, we consider the quantum architecture design for both Hermitian and non-Hermitian matrices. However, for the largedimensional-matrix simulation, we consider the Hermitian case to analyze the error. We show that the Hamiltonian matrix T N will possess an efficient low-complex quantum simulation while it is expressed with a sparse decomposition technique.
Note on Non-Hermitian Toeplitz Matrix for Preparing the Hamiltonian: For the non-Hermitian Toeplitz matrix, there is the classical procedure to prepare a unitary similar to a Hermitian matrix (as discussed in [45] where J N is the N × N backward identity operator (see [45 Th. 3.3]). The elements of the Hermitian matrix S N is given as where t h−k := t[h − k], and i = √ −1 is the complex imaginary number. It can be further noted that T N can be viewed as a linear combination of the N × N Jordan form J N (0), the transpose of J N (0), and its power. The Jordan form J N (0) can be written as For the quantum simulation, we may require a possible Hermitian Hamiltonian representation for every Jordan form J N (0). The Jordan form J N (0) is unitarily similar to the Hermitian matrix S N (0) via the operation of U N = 1/ √ 2(I N + 5500423 VOLUME 4, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. Hence, we may express that a Toeplitz matrix is unitarily similar to a Hermitian matrix (see [46,Th. 4]), with a possible Jordan form representation. A Note on Sparse-Decomposition of the Toeplitz Matrix: Based on the literature [45], [46], it is evident that a Hermitian Toeplitz matrix can be expressed in a Jordan normal form, which is a sparse decomposition of the Toeplitz matrix. Further, we see that a non-Hermitian Toeplitz matrix possesses a similar decomposition with a Jordan form and its Hermitian form, which we discuss next.

1) SPARSE AND JORDAN-CANONICAL REPRESENTATION FOR THE HAMILTONIAN T N
For the matrix T N , we may augment the sparse representation of the Jordan form for the efficient Hamiltonian simulation. Matrix T N ∈ C N×N can be represented in terms of linear combinations of sparse matrices (Jordan forms, and identity matrices) and their powers as follows: where J N = J N (0) ∈ C N×N is the Jordan operator corresponding to a zero-eigenvalue of some matrix [45]. We will exploit this representation under the quantum framework for an efficient quantum circuit realization for the Toeplitz matrix T N .

2) PROPOSED QUANTUM REALIZATION CIRCUIT FOR A TOEPLITZ HAMILTONIAN OPERATOR
We show the circuit realization of a Toeplitz operator of dimension 4 × 4 using a proposed composite quantum Jordan gate, and the extension to the generalized case. For a 4 × 4, we propose a quantum Jordan operator as given by the following lemma. Lemma 1: A composite quantum gate can be realized for the Jordan operator J N for N = 4 using the elementary quantum gates as where σ 0 , σ x , l u , and l L denote identity, Pauli-X, raising ladder, and lower ladder gates, respectively. Here, the elementary quantum gates (2 × 2 operators) are discussed as Proof: The expression in the right-hand side of (13) can be realized as follows: Similarly, it holds true for J † 4 using gates l u , σ 0 , and σ x . In Fig. 1(a A Toeplitz matrix can be obtained using the proposed Jordan operators as The first term t[0] × I 4 generates the diagonal terms of the Toeplitz matrix. Matrix J 4 and its different powers are given as follows: Engineering uantum Similarly, we can get the different powers of matrix J † 4 which will give us realization for the off-diagonals below the main diagonal of T 4 . Following (12) for N = 4, we show a proposed quantum architecture for the representation of Toeplitz operator T 4 in Fig. 1. In our circuit representation, the dark dots represent connections, the circle ⊗ with notation "T " denotes tensor operation, and circle ⊕ denotes adder circuit (addition or subtraction operation are denoted by the sign of the inputs). For the composite J 4 gate, we need two lower ladder operators, and one Pauli-X and one Pauli-Y gate, respectively. For the implementation of J † 4 , we need additionally two raising ladder gates. For a 4 × 4 Toeplitz matrix with symbols, [t 3 , . . . , t 0 , . . . , t 3 ] can be implemented using the proposed Jordan gates, as shown in Fig. 1(b). Here, I 4 is prepared first using the tensor product of two 2 × 2 identity gates.
Note: One can prepare both Hermitian and Non-Hermitian Toeplitz matrices using the same circuit as shown in Fig. 1(b), by changing the input symbols only (for Hermitian

b) Realization of a Toeplitz matrix of dimension N × N with the proposed quantum Jordan gates
We can extend the circuit realization for the generalized case with dimension N × N.
Lemma 2: For a Hermitian Toeplitz matrix T N = T † N ∈ C N×N for N > 4 and N = 2 n q with input qubit size of n q , there exists a recursive representation for the Hermitian Toeplitz operator using quantum-Jordan gates given by where σ ⊗N 0 is the identity operator of dimension N × N, and the Jordan operators J N and J † N are the matrices of the form Proof: A Toeplitz matrix can be represented with Jordan operators as shown in (12). For the first diagonal, we can prepare an identity operator of size N × N, where N = 2 n q for input qubits n q . We will show the recursive relation for the Jordan operators using the principle of mathematical induction.
For N = 8 with input size n q = 3, we will use the elementary quantum gates as discussed in Section II-B1 to prepare the operators for the higher dimension. Operator J 4 , and J † 4 can be implemented using Lemma 1 [as shown in Fig. 1(a)]. Operator l u N 2 := l u ⊗ l u is a 4 × 4 raising ladder operator, and 0 4 is a zero-matrix of dimension 4 × 4. The right-hand 5500423 VOLUME 4, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.

Engineering uantum
Transactions on IEEE side of (19) can be expanded for N = 8 as Similarly, we get gate J † 8 using J † 4 , 0 4 , and lower ladder gate l L4 = l L ⊗ l L of size 4 × 4. Using the principle of mathematical induction for n q = n q + 1 input qubits, one can implement the recursive operator for J N and J † N . Corollary 1: For a family of rank-deficient Toeplitz matrices, T N of size N × N and rank r < N, the decomposition given by (12) exists, and there is an efficient quantum algorithm that can simulate every T N ∈ C N×N .
Proof: For every T N ∈ C N×N with rank r < N, there is at least an eigenvalue 0 ∈ λ(T N ) for which we can find a Jordan canonical block in the form of J N (0), as given in (10). Hence, there will be an efficient sparse decomposition of T N in the form of J N (0), as given by (12), which can be implemented by a quantum simulation algorithm efficiently.
Note on Circuit Depth and Gate Counts: In most of the recent quantum simulators, such as IBM-QISKIT [47] and Google Cirq [48], the elementary gates such as Pauli gates, and Hadamard gates are available. The ladder gates (which are often used in quantum photonic circuits) can be implemented using combinations of Pauli gates on the quantum simulator. We discuss the circuit depth and the number of elementary gates required for the implementation of the Toeplitz matrix given as follows. Note that the depth of a quantum circuit is a metric that calculates the longest path between the input and the output of the circuit. 1) A 2 × 2 Toeplitz matrix is a simple Hermitian matrix, and hence, it can be implemented with one Pauli-0 and one Pauli-x gate. Here, the circuit depth for the Toeplitz matrix is 2. 2) A 4 × 4 Toeplitz matrix needs two Pauli-0 gates for the preparation of its principal diagonal. Further, we need to prepare six off-diagonals using Jordan gates as given in (16). To prepare a 4 × 4 Jordan gate J 4 , four elementary gates are required (one Pauli-0, two lower ladder gates, and one Pauli-x gate, respectively). Similarly, we need four elementary gates for preparing J † 4 (one Pauli-0, two upper ladder gates, and one Pauli-x gate, respectively). From Fig. 1(a), we can see that the circuit depth for every J 4 gate is 4. 3) For a Toeplitz matrix T 8 , we need three Pauli-0 gates to prepare its diagonal structure, i.e., I 8 = σ ⊗3 0 . We need to prepare 14 Jordan gates (including seven symmetric Jordan gates) J 8 , and J † 8 ∈ R 8×8 . Here, for the preparation of 8 × 8 Jordan gate using the recursive algorithm, we need two J 4 gates and two elementary l u gates. 4) The quantum architecture of a Toeplitz matrix with dimension N × N (considering N as the power of 2) can be prepared with elementary quantum gates as follows. a) To prepare the principal diagonal, we need log 2 N number of σ 0 gates, b) We have (2N − 2) off-diagonals that can be prepared with combinations of J N gates including (N − 1) number of J † N gates. Considering each Jordan gate J N as a unit, the depth of the circuit is N. Using the recursive implementation of the Jordan gate, we need two Jordan subcircuits J N 2 and an additional upper ladder operator of Here, to implement a ladder operator l u N 2 , we need log N 2 elementary l u of dimension 2 × 2 gates. 5) The number of elementary gates required for the simulation of a unitary matrix U(2 n ) is given by (n 3 4 n ) (see [49,Sec. VIII]). Considering, N = 2 n , the gate complexity is approximately (N 2 (log N) 3 ). In our case, the depth of the circuit for a Toeplitz matrix of dimension N × N is given by N, and the number of elementary gates is in polynomial with log N approximately. Hence, the approximate gate complexity in simulating a Toeplitz matrix with elementary gates is given by˜ (NPoly(log N)).

3) EFFICIENT HAMILTONIAN SIMULATIONS WITH JORDAN FORM AND IMPLEMENTATION OF THE UNITARY OPERATOR
For the estimation of eigenspectrum λ(T N ) of a generalized Toeplitz matrix T N , we will perform the Hamiltonian simulation-based QPE method. For the QPE, we need to prepare a unitary operator from the Toeplitz matrix given as follows: The above unitary matrix U N is expressed in an exact form with a matrix exponential function, which is often very difficult to simulate with realizable quantum gates. The Hamiltonian simulation algorithm helps to find an approximate unitary matrixŨ N closer to the exact operator U N (with the minimum difference in the spectral norm). For the practical implementation of this approximate unitary operator via the Hamiltonian simulation, we consider a precision T given by The problem in (23) is addressed by suitable quantum frameworks, such as product formula (e.g., Trotter-Suzuki approximation) [38], [39], quantum walk [41], quantum signal processing [40], and Taylor series approximation [31].
Here, we will use the truncated Taylor series approximation technique for the quantum simulation of the Toeplitz matrix. The approximation of the Hamiltonian up to orders VOLUME 4, 2023 Engineering uantum With a bigger size of L, the Taylor series approximation error T can be minimized.

4) QPE FOR SPECTRUM ESTIMATION OF T N
QPE can be implemented in several ways [50], [51] with an objective to estimate phases θ j of eigenvalue λ j for j ∈ [N] for the unitary operatorŨ N , given as Here, |v j ∈ C N is an eigenstate ofŨ N ∈ C N×N (known in prior from oracle) and the corresponding eigenvalue is λ j . The eigenvalue estimation for Hamiltonian T N is now mapped to phase estimation for the unitary operatorŨ N . The efficient implementation of QPE needs unitaryŨ N to be conditioned on the state of the ancillary qubits in a controlled fashion given by where I represents identity operator, and |0 and |1 are computational bases for the ancillary qubits. The superposition states of the qubits are generated by Hadamard gates, while the controlled operation is performed by U-rotation gates.

5) QUANTUM TIME RESOLUTION
The time-evolution of a quantum state through Hamiltonian simulation requires a time resolution for the preparation of unitary operatorŨ N ∈ C N×N . Here, we define a terminology called, quantum time resolution (QTR) for Hamiltonian simulation.
Definition: The minimum amount of time τ required to prepare the unitary operatorŨ N ∈ C N×N for the Hamiltonian simulation of a Hermitian matrix A QTR is a sufficient amount of time, such that the expected error in the simulation converges asymptotically.

C. CIRCULANT APPROXIMATION-BASED QFT FRAMEWORK
A Toeplitz matrix is asymptotically (for N → ∞) equivalent to a circulant matrix C N with bounded error C N − T N . Here, C N and T N are asymptotically equivalent and both are bounded as follows [9]: where A is the approximation error in generating the circulant matrix from original T N . Note that matrix C N may become sparse if there are sufficient zero (0) entries in c. The generating functiont( f ) for a Toeplitz matrix T N following Szego's theorem [9] is given bỹ where t[k] is the Fourier series expansion for the generating functiont( f ), given as follows: Szego's method for the Toeplitz matrix in an asymptotic case follows: where λ(T N ) is the spectrum of T N . For practical estimation of the spectrum λ(T N ) via quantum algorithm,t is approximated at first by partial Fourier sum [10] given by for some finite N. Circulant matrix C N is constructed in a way that it possesses eigenvalues λ i for i ∈ [N], which are samples of P N−1 ( f ), and are denoted as P N−1 ( n N ). We are interested to compute the first row c = [c 0 , c 1 , . . . , c N−1 ] of the circulant matrix C N , as spectrum λ(C N ) for the circulant matrix C N can be easily computed from c N . Vector c N can be generated from P N−1 ( n N ) as a sequence given by Note: In (32), we have considered the fact that , C N becomes Hermitian. For finite N, the error analysis for the approximation of C N from a general Toeplitz matrix T N is discussed in Section IV.

1) QUANTUM STATE PREPARATION
After generating vector c N ∈ C N as an approximate alternate of T N , the first step is to prepare a quantum state vector |c 5500423 VOLUME 4, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
A specialized memory model called quantum random access memory (QRAM) helps to encode state |c in quantum superposition. The oracle QRAM that is the standard model for qubit encoding requires time complexity of O( √ N nz (c)), where N nz (c) denotes the number of nonzero elements in the sparse-vector c. Further reduction of time-complexity for vector and matrix encoding can be done with the augmented-QRAM method [52], where an ancillary qubit with superposition state is prepared with conditional rotation, amplitude amplification, and quantum key-value mapping techniques. With augmented QRAM, the below transformation can be implemented as a unitary operator as follows:   [9]. Suppose, U = [u (1) , u (2) , . . . , u (N ) ] be the unitary matrix, and = Diag(λ m ) is a diagonal matrix. Hence, it follows that C N = U U † is unitarily similar to a diagonal matrix, and hence, it is normal. It is observed that only the first row of C N is required in the eigenvalue expression, instead of the whole matrix. As we have the quantum state vector |c prepared via some QRAM method, we propose the below statement.
Statement: For a circulant matrix C N ∈ C N×N , eigenvalues λ m for m ∈ [N] can be computed by the application of QFT on the state vector |c ∈ C N without implementation of the quantum eigenvalue decomposition (EVD) method. The below expressions show the QFT conversion of probability amplitude c j to λ as follows: Algorithm 1: Quantum Jordan Gate-Based Spectrum Estimation of a Toeplitz Matrix.
break; 8: End If 9: End For 10: For j = 1 : N, do: Here, the QFT works as a linear operator for the conversion of probability amplitudes from basis | j to |l . For the efficient implementation of QFT for |c , we seek N = 2 n for n ∈ Z.

III. PROPOSED ALGORITHM A. QUANTUM SPECTRUM ESTIMATION OF A TOEPLITZ SYSTEM USING JORDAN-FORM-BASED REPRESENTATION
In Algorithm 1, we have presented a pseudocode for the spectrum estimation of a Toeplitz matrix T N with Jordonform-based sparse-representation.
Note: With a slight abuse of notation, we consider "spectrum estimation" as a quantum eigenspectrum estimation. In our subsequent algorithms, we have shown QPE based on eigenspectrum estimation as follows.  Description of Algorithm 1: Among the input variables, T N denotes the Toeplitz matrix, N = 2 n represents the size of the spectrum λ(T N ); U N , τ th , T , |ψ , and L are exact unitary matrix (classically prepared), time of evolution, precision, quantum state vector, and the number of terms for the Taylor approximation, respectively. τ th can be obtained based on the experimental setup as detailed in the numerical section. Note that v j is an eigenstate ofŨ N prepared by an oracle, with eigenvalue e 2πiθ j (where θ j is unknown) for j ∈ [N]. In the proposed algorithm, we have shown a modified Hamiltonian simulation for the preparation of unitary operator U N . Here, we have incorporated Jordan's decomposition-based representation of the Toeplitz matrix T N for the preparation of U N . The Taylor series expansion gives an approximate unitary evolution U N of the Hamiltonian T N , with L-terms (with the tradeoff between complexity and precision). The standard QPE method is performed on oracle-prepared eigenstate v j to estimate the eigenvalues of the unitary operator U N . Here, an N-length state vector is initialized with tensor representation. Next, we perform Hadamard operations for the superposition of the ancillary qubit φ 0 . The superposition state gives the quantum advantage of parallelism to augment quantum speed in the computation. Next, we create quantum interference through successive controlled-U gate operations on the superposition states to control phase rotation on |1 -computational bases. Further, IQFT is executed to the controlled-state vectors via the auxiliary resisters |h . The measurement is performed at the output of the IQFT circuit on a computational basis to detect the phase changes of the qubits acting on the quantum eigenstates. The eigenvalues are obtained from the measurement that is stored in variablẽ λ j for j ∈ [N]. Here, we have shown steps 23-26 as steps of the quantum measurement. Here, n q = log N is the required number of qubits for the quantum experiment. The quantum shots are referred to the number of experiments performed to get the histogram. We have specified the number of quantum shots by M 1 . Function Measure(m, n q ) performs the measurement of the quantum state (at the output of IQFT) on computational bases for M 1 shots and n q qubits. The computation of the probability from the measurement is discussed in Section IV-B.
A Note on Fig. 2: We have shown a block-level circuit description of Algorithm 1 for the spectrum estimation of a Toeplitz-structured matrix in Fig. 2. Here, we have shown the sparse-decomposition-based representation of the Toeplitz matrix using identity and Jordan operators for N × N dimension. A detailed circuit of this part (implementation of the Jordan operators and the Toeplitz operator) for a smaller dimension of 4 × 4 is shown in Fig. 1. In the Hamiltonian simulator subcircuit, the approximate unitary operatorŨ N is prepared. In the QPE circuit, the approximate unitary operator (in our caseŨ N ) is conditionally applied to the phase registers to raise its power with 2 0 to 2 t for t qubits for the initialization of state |0 following the approach as shown in [37,Ch. 5.2]. Through the QPE circuit, the phases of the approximated unitary operator are estimated. Here, H denotes Hadamard gate, c denotes the classical register, |0 is the controlled qubit, and |v represents the eigenstate. Measurements are performed at the output of IQFT on a computational basis and stored in the classical registers.
A Note on Input Eigenstates: For the computation of the eigenvalue spectrum of the Toeplitz matrix using Algorithm 1, we have assumed that the input eigenstates |v j are known. For a circulant matrix, these states are known a priori. However, the eigenvectors are not known for the generalized Toeplitz matrix in advance. Hence, input eigenstates need to be prepared for the QPE algorithm, which must be sufficiently close to an actual eigenvector [20]. In [53], the authors have shown quantum subroutines for finding the eigenvectors of a Hamiltonian operator and discussed its application to atomic physics. A further refinement of eigenvector computation with a small number of quantum gates is shown in [54], where a classical eigenvector is taken initially to prepare it on a quantum state using fine grid approximation. Recently, an iterative QPE method is proposed to estimate phase θ using polarization qubits described by a higher dimensional Hilbert 5500423 VOLUME 4, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.  [55]. In [20], the authors have shown theoretically that for a sparse and Hermitian operator, the generalized eigenvalue problem can be solved.

B. QUANTUM SPECTRUM ESTIMATION OF A CIRCULANT APPROXIMATE TOEPLITZ SYSTEM
We have shown a pseudocode for the spectrum estimation of the Toeplitz matrix via circulant matrix approximation in Algorithm 2.
Description of Algorithm 2: For a given Toeplitz matrix with structure generating symbols t[k] for k = N − 1, . . . , 1, 0, −1, . . . , −(N − 1), vector c N containing the symbols of the circulant matrix can be generated using (33). We will see that error C N − T N is bounded by desired precision for the Toeplitz to circulant approximation. For a circulant system C N , term c N is its first row. We take symbols t[k] as input of the algorithm, and we are interested to find the spectrum λ(C N ) of the circulant system, generated using the symbol conversion as shown. Before further processing, we need to prepare the quantum state vector |c with normalization followed by quantum superposition, as shown in (34). The augmented QRAM technique can be used to store the quantum state vector |c in the memory efficiently. After preparing the state vector |c , we are ready to perform the IQFT operation. We apply the Fourier bases | j and apply on state |c to get the amplitudes λ l on bases |l for l ∈ [N]. Here, we have given steps 17-20 for the quantum measurement. Here, we take input qubits of size n q = log N for the quantum experiment. M 2 number of quantum shots are performed to get the histogram. Here, Function Measure(m, n q ) performs the measurement of the quantum state (at the output of IQFT) on computational bases for M 2 shots and n q qubits.
The computation of the probability from the measurement is discussed in Section IV-B and the the result in Section VII-E.

IV. ERROR BOUND IN ESTIMATED SPECTRUM FOR CIRCULANT APPROXIMATION
We incur an approximation error A for approximating the circulant variant C N from the Toeplitz matrix T N given by In the asymptotic case (N → ∞), the approximation error A → 0. However, for finite N, we get an estimated spectrum λ(C N ) = {λ c l for l ∈ [N]} for the circulant matrix C N obtained from the Toeplitz matrix T N . In this spectrum estimation process, the true eigenvalues of the Toeplitz matrix can be written as whereλ c ∈ λ(C N ) denotes the estimated eigenvalue, λ t ∈ λ(T N ) represents the true eigenvalue, QFT is the quantization error in the QFT process, and Q is the estimation error in the eigenvalues for the propagation of the approximation error A in generating the circulant matrix. Our objective is to quantify the estimation error based on (38) given by Note that we have used . to denote the 2-norm of the corresponding vector. Computing the error bound for QFT , and Q separately is often numerically intractable. Hence, we define a variable for the combined errors as

A. ERROR BOUND ON QQ
For finding the bound QQ , we consider that the Toeplitz matrix is Hermitian (i.e., T † N = T) and it is a large dimensional (N → ∞) random matrix whose elements are independent and identical samples taken from a circularly symmetric Gaussian distribution with zero mean and variance of real and imaginary part be σ 2 T 2 each, as follows: Following (33), the elements c[k] that are drawn from T N also follow the circularly symmetric Gaussian distribution with zero mean and modified variance for real and imaginary parts to be σ 2 C 2 each, given as follows: The estimation error QQ can be written (avoiding the error for implementation of the QFT circuit) as In general, the joint probability density function p(λ t ,λ c ) is not tractable for a large matrix. Hence, we consider the We assume that the joint probability density function p(δ, δ 1 ) for the variables δ and δ 1 are independent, which can be represented as where J denotes the Jacobian matrix. The Jacobian matrix for (44) is with |J | = 1. Now, the error probability function p(δ) can be written as First, we will find the probability function p(λ c ), i.e., p(δ 1 − δ) with the following lemma. Note: In many signal processing and communication applications, the system matrix can be Hermitian. However, the elements of the matrix need not be Gaussian, and independent identically distributed. For the large dimensional Hermitian matrix, as in our case, the probability density function of eigenvalues follows that of the Wigner matrix [56, Sec. II-B1], irrespective of the distribution of individual elements in the matrix.
For a large-dimensional matrix, the probability density of λ c [n] will follow Gaussian distribution (applying the central limit theorem). Here, c[k] follows the independent and identical, complex circularly symmetric Gaussian distribution CN (0, σ 2 C ) from (42). Hence, the sum of the elements c[k] ∀k ∈ [N] for every eigenvalue will follow the circularly symmetric Gaussian distribution with zero mean and variance of Nσ 2 C . We can write the density function p(λ c ) as Replacing λ c by δ 1 − δ in (49), we get the distribution p(δ 1 − δ) as follows: Next, we will find the probability density function p(λ t ), i.e., p(δ) for the Hermitian Toeplitz matrix T N = T † N for a large dimension (N → ∞). The eigenvalues λ t [ j] ∈ λ(T N ) ∀ j ∈ [N] for the Toeplitz matrix T N ∈ C N×N for large N will follow the distribution given by [56] p(λ t ) = As we have taken T N = T † N , the eigenvalues will lie in the range 0 ≤ λ t ≤ 2σ T , considering that the elements t[k] for k ∈ [N] has zero mean and variance σ 2 T . Proposition 1: The probability density function of the total estimation error, p( QQ ) for a Hermitian and large random Toeplitz matrix T N of size N × N is bounded by the following expression (taking Taylor series approximation): Proof: Using expressions (50) and (51) in (47), we get the integral as We take the Taylor series approximation for the exponential function to avoid intractability of the higher order terms as follows: Now, p(δ) can be written as follows: where the integrals where a = 2σ T , and x = δ 1 . Hence, using result (58), integral (56) becomes We compute integral (57) as Integral I 4 := 2σ T 0 sin −1 δ 1 2σ T dδ 1 has form a 0 sin −1 x a with a = 2σ T and x = δ 1 , which can be computed as follows: Hence, (62) and (57) become Hence, probability p( QQ ) is approximated as with the constraint that √ 2π Nσ C > 0. Variance of QQ : For computing the moments (mean and variance) of QQ , we consider a truncated probability density for range (− B , B ). Value B depends on the system bit resolution and parameter estimation error using quantum circuits. The mean of the total error QQ is given by The variance of the total error QQ can be computed as follows: Engineering uantum where the constants are c 1 := μ , c 2 := 1 √ 2π Nσ C , and , respectively. Hence, the bound of the total error can be given by (69), shown at the bottom of this page. Note on QFT Precision: The quantization error QFT depends on the number of bits to represent an element from array c N . If the distance between two consecutive quantization levels is Q , the total variance of the error in representing N number of elements with a precision of N q -bits is given by where the quantization error Q = 1 2 (Nq−2) . By increasing the number of qubits, the quantization error QFT can be reduced.

B. MEASUREMENT NOISE AND PROBABILITY OF CORRECT MEASUREMENT OUTCOME
Each eigenvalue is estimated from a quantum algorithm using approximated QFT and QPE. The experiment needs to be repeated multiple times to maximize the probability of the outcome of a measurement basis, mitigating the measurement noise. Here, we consider that the quantum circuit is simulated for m ∈ [M] times, with one or multiple rounds given by r = 1, . . . , R n . The ancillary qubits are prepared in a superposition state through the application of Hadamard gates controls U τ r (26) with integer variable τ r per round. The ancillary qubits are read out in the X-basis and give measurement outcome m r ∈ {0, 1}. The total number of controlled-U rotation over M experiments can be given by T = M n=1 R n r=1 τ r . However, for each experiment, the number of controlled-U rotations is = R n r=1 τ r . The coherence length coh is defined in [57] as Note that coh limits the maximum of the controlled-U operations per experiment with ≤ coh . In (71), τ err N q represents time-to-failure of N q qubits, τ U denotes the time required to implement single controlled-U gate, and τ err is the time-toerror of a single qubit. Hence, the value of is a crucial parameter that is related to the QTR of the quantum hardware as discussed in Section II-B4.
In our simulation, we have considered that the ancilla qubit is rotated by R z (γ 1 ) = e −iγ 1 /2 . Here, γ 1 is the phase for rotation around the Z-axis direction for the controlled-rotation U τ 1 . Hence, the eigenstate output from the operator U τ 1 is where β j is the normalizing coefficient for eigenstates |ψ j corresponding to eigenvalues λ j . There is uncertainty in getting the result on a particular X-basis as the experiment is repeated for M-times. The probability of measurement of the first ancillary qubit (say c 1 ) on m r ∈ {0, 1} basis can be written as where D j ≡ |β j | 2 , and the outcome-state of the ancillary qubit is given by Thus, the probability of correct measurement outcome for the first ancillary qubit is given by where λ k ∈ λ(T N ) denotes the kth eigenvalue.
We have already assumed that the quantum system has N a -bits resolution and the probability of correctly measuring every single bit is determined as Pr b in (75). Let y i be the decimal value of the ith binary string of length N q . Let y k i,l be an N q length binary string with l bits reversed with respect to y i and k indicates the kth realization of such string.
Proposition 2: The variance of the quantum measurement noise w m for a N q -bits quantum resolution system after the 5500423 VOLUME 4, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.

Engineering uantum
Transactions on IEEE quantum simulation task can be computed as follows: Proof: The probability of j-bits mismatch for y i is Pr It is considered that the sign bit is equally prone to flip incurring an error. Therefore, the mean of the error will be zero. There will be N q j number of binary strings having j-bits in mismatch with respect to y i and l = 1, 2, . . . , N q . The index i varies as i = 1, 2, . . . , 2 N q . Therefore, the average value of the error can be given by (76).

V. COMPUTATIONAL COMPLEXITY
We will consider that eigenvalues λ l ∈ λ(T N ) of the Hamiltonian T N (or equivalently the phases θ j of the eigenvalues of unitaryŨ N ) has the binary b-bits representation. Here, the overall complexity of estimating N eigenvalues is the sum of N × (3b + b 2 /2) and the complexity of the Hamiltonian simulation [37]. With Taylor series approximation [31], the complexity of the Hamiltonian simulation is given by where n q denotes the number of qubits required for simulating a d-sparse Hamiltonian with simulation error bounded by T , and the term κ = d 2 T N max τ . For a dense matrix, the Hamiltonian simulation complexity is significantly large and the quantum speed-up may be compromised. Next, we will discuss the modified Hamiltonian simulation for a possible quantum complexity advantage. We will further see that the circulant approximation technique will produce guaranteed complexity improvement by avoiding the QPE approach for the estimation of the spectrum for bounded error.

A. COMPLEXITY FOR HAMILTONIAN SIMULATION OF A TOEPLITZ MATRIX REPRESENTED IN SPARSE AND JORDAN-CANONICAL FORM
The form in (18) is a generalized decomposition of a Toeplitz matrix in Jordan canonical form-based Hermitian representation that can be used in many quantum simulations with required normalization or modification. Here, its computational complexity for the Hamiltonian simulation framework is discussed. Proposition 3: The complexity for the Hamiltonian simulation of a Hermitian Toeplitz-structured system T N ∈ C N×N in form (18) to prepare a unitary operatorŨ N for time τ with the error in spectral norm given by Ũ N − exp (−iT N τ ) ≤ T can be approximated with Taylor series method as where ω denotes a constant term, and n q is the required number of qubits. Proof: Here, the simulation complexity for t[0]I N is of orderÕ(1) as it needs to simulate the first element of the Toeplitz matrix multiplied with identity operator. As T N = T † N , we need onlyÕ(1) complexity for each t[i], where i = 1, . . . , N. Hence, the simulation complexity for N n=1 t[n] J is (N − 1) ×Õ (1). Given that the number of qubits n q is taken as N = 2 n q , the Hamiltonian complexity for n q -qubits resolution following Taylor series approximation can be computed as follows. In

B. COMPUTATIONAL COMPLEXITY FOR SIMULATING A TOEPLITZ MATRIX WITH CIRCULANT APPROXIMATION
The circulant matrix C N approximated from a Toeplitzstructured Hamiltonian T N can be generated using a single array c N = {c[k] for k ∈ [N]} with approximation error bounded by A as shown in (33). The array vector c N represents the first row of the circulant matrix C N , which generated all the below consecutive rows of the matrix by a right cyclic shift of the array elements. The following lemma shows the complexity in obtaining the estimated spectrum λ(C N ) for matrix C N . Lemma 4: For an n q -qubit representation of the input array c N of size N = 2 n q , which generates the circulant matrix, the eigenvalues λ c l ∀ l ∈ [N] can be obtained using an ideal n q -qubit unitary QFT circuit with gate-complexity of O(n q ( n q +1 2 )), while an approximate QFT circuit comprising Clifford and T -gates may incur a gate-complexity of O(n q log(  [37, p. 220]. It shows that the QFT circuit-based spectrum estimation method can provide exponentially faster operation while running on a quantum computer than the classical FFT-based eigenvalue estimation technique. In approximate quantum Fourier transform (AQFT) techniques, the rotation gates with rotation angles smaller than a certain threshold value are removed. Recently, an n q -qubit AQFT circuit shows that circuit built on T -gates can reduce the QFT complexity from O(n 2 q ) to O(n q log(n q )) operations [58]. While we incur an approximation error bounded by A , the complexity of eigenvalue computation by the AQFT method becomes A complexity comparison for different Hamiltonian simulation algorithms for a dense Toeplitz matrix is shown in Table 1.

VI. APPLICATION IN ARRAY SIGNAL PROCESSING
Spectral estimation for a Hermitian Toeplitz matrix is pervasive in many array processing applications [1], [17], [19]. We have considered a DoA estimation problem here to demonstrate the application of our proposed quantum framework in modern array processing.
We assume that P narrow far-field signals impinge on a uniform linear array (ULA) with M (M > P) antennas from the directions of θ = {θ [P] } simultaneously at time t, as shown in Fig. 3. The received signal vector at the array output can be written as Here, x(t ) = [x 1 (t ), x 2 (t ), . . . , x P (t )] † denotes the incident signal vector and the received signal vector y(t ) = [y 1 (t ), y 2 (t ), . . . , y P (t )] † , v(t ) represents the additive noise vector (assumed to be circular complex Gaussian white noise), and A(t ) = [a(θ 1 ), a(θ 2 ), . . . , a(θ P )] † is the array manifold matrix with manifold vectors a(θ i ) = [1, e −iπ sin(θ i ) , e −i2π sin(θ i ) , . . . , e −iπ (M−1) sin(θ i ) ] † (assuming the distance between antenna elements be d x = λ x 2 for wavelength of λ x ). The source signal x(t ), and noise vector v(t ) are assumed to be uncorrelated given by where σ 2 x = [σ 2 1 , σ 2 2 , . . . , σ 2 P ] represents the power parameter of x(t ), σ 2 v denotes the noise variance, and δ t 1 ,t 2 = 1 for t 1 = t 2 or 0 elsewhere. The covariance matrix of the received signal vector can be obtained as where T x = P i=1 σ 2 i a(θ i )a(θ i ) † = A Diag(σ 2 x )A † is a Hermitian Toeplitz matrix of the source signals. The expression 5500423 VOLUME 4, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.   (83) shows that the covariance matrix T y is also a Hermitian Toeplitz matrix. We take an estimated Toeplitz matrix in practical applications for L number of snapshots given bŷ For the spectrum estimation of the estimated Hermitian Toeplitz matrix given in (84), the classical TMRA approaches perform EVD techniques. We can apply our proposed Algorithm 1 here for spectrum estimation with quantum complexity advantage. The performance of the classical TMRA and the proposed quantum algorithm is shown in Section VII. Note: The standard DoA estimation method (as shown in [1]) can be applied after the quantum spectrum estimation for finding DoAs. Here, we restrict our result to spectrum estimation problems to augment the complexity advantage (as discussed in Section V).

VII. NUMERICAL RESULTS
In this section, we have obtained numerical results related to the proposed algorithm for the Toeplitz-structured Hamiltonian simulation and its spectrum estimation. Here, we have shown the simulation result of our proposed algorithm in sections, partly performed on a classical computer and partly on a real-time IBM quantum simulator (IBM "Statevector" quantum simulator). We have shown the numerical results of estimation error and algorithmic complexity (in comparison to the standard Hamiltonian simulation) in Sections VII-A-VII-C, which are performed on a classical simulator (MATLAB-2021b) following the parameters, as shown in Table 2. Further, the application of the proposed algorithms on a TMRA problem is discussed in Section VII-D. Finally, we have shown the implementation of the algorithm of finding eigenvalues of a Toeplitz matrix on an IBM quantum machine (IBM "Statevector" simulator) in Section VII-E. We have limited the application to find the second principal eigenvalue and the minimum eigenvalue with the highest eigenvalue being 1. As we have limited quantum resources presently for academic research, we have shown the analysis part (which requires simulation of large-dimensional matrices) on the classical computer and a small-scale application of the proposed algorithm on the real-time quantum simulator.

A. ESTIMATED SPECTRUM AND ERROR-PROBABILITY FOR THE CIRCULANT APPROXIMATION
Here, we have shown a numerical simulation framework for the circulant approximation of a Toeplitz matrix in Wigner form, whose elements follow independent and identically distributed Gaussian probability distribution, with zero mean and unit variance. We have given the spectrum estimated by Algorithm 2 and computed the error probability. Fig. 4 shows the plot of eigenvalues (normalized) of a Toeplitz matrix T N in Wigner form, where the matrix dimension is taken to be 32 × 32. Here, λ t i for i ∈ [N] denotes the true eigenvalues (computed using the classical eigenvalue decomposition technique), and λ c i for i ∈ [N] represents the estimated eigenvalues using proposed Algorithm 2. The eigenvalues are plotted here in ascending order. It can be observed that the circulant approximations for larger N will produce almost the same eigenvalues with quantum simulation for a Toeplitz matrix.
We increase the matrix dimension to ensure the convergence of the error variance for the approximation method used in Algorithm 2. We have shown the probability of approximation error Pr(δ) with varying matrix dimensions in Fig. 5. Based on the result in (66), an empirical result is shown here, taking the variance σ C = 1, and the matrix dimension is varied from N = 10 to 1000 with increment 30. For a dimension of 1000 × 1000, the error probability is Pr(δ) < 12 × 10 −3 approximately. Further, the figure shows that for the increase in matrix dimension, the error probability decreases.

B. EFFECT OF QTR ON ESTIMATED SPECTRUM FOR THE TOEPLITZ-STRUCTURED HAMILTONIAN SIMULATION
We have performed the eigenvalue spectrum estimation problem for a given Hermitian Toeplitz matrix using the proposed Algorithm 1 and compared it with the classical EVD approach. The proposed simulation technique depends on the   evolution time parameter or QTR (τ ). Here, we have taken an expression for the QTR τ given by where η is a constant (controlling parameter for the QTR). The different choices of η are considered in the numerical simulations, as the ideal value of the τ is often unknown. For improper choice of QTR, we have seen the estimated spectrum. deviates from the true spectrum estimated classically. In Fig. 6, we have shown the normalized eigenvalues for a Toeplitz matrix of dimension 8 × 8 with elements taken from a normal distribution. For τ 1 , τ 2 , and τ 3 , the values of η are taken as 2, 1, and 0.1 respectively, and T N 2 = 1. Here, the estimated spectrum (using the proposed Algorithm 1) has deviated significantly from the true eigenvalues for different values of the QTR. Here, we have seen the estimated eigenvalues are closest to the true eigenvalues of the Toeplitz matrix for the parameter τ 1 as compared to τ 2 and τ 3 .

C. COMPARISON OF COMPUTATIONAL GATE-OPERATION COMPLEXITY
An empirical simulation of gate-operation complexity of the proposed structured Hamiltonian simulation as given in (78) in comparison with standard Hamiltonian simulation is shown in Fig. 7. Here, the term ω is taken as 4 approximately, τ = 0.48, T max = 1 and the simulation error is taken to be T = 0.01. In this plot, the computational complexity between the proposed method and standard Hamiltonian simulation is computed for different sizes of inputs n q = log 2 (N), where N × N is the dimension of the Hamiltonian matrix. As N grows, the structured Hamiltonian simulation shows significant complexity improvement as compared to standard Hamiltonian simulation. For a matrix dimension of 2 10 × 2 10 , the proposed structured Hamiltonian simulation-based spectrum estimation algorithm takes 468 gate operations approximately, whereas the standard Hamiltonian simulation takes 1400 gate operations approximately. The proposed framework enhances about 66.57% operation complexity advantage for a qubit size of 10. Hence, embedding a Toeplitz-structured Hamiltonian in the simulation shows efficient quantum operations and reduced quantum circuit complexity.
The gate-operation complexity in the circulant approximation method as shown in Algorithm 2 is given in Fig. 8 and compared with the complexity of the standard Hamiltonian method, as well as the Toeplitz-structured Hamiltonian simulation method proposed in Algorithm 1. The AQFT-based approximation method takes the least quantum gate operations in comparison with the other methods, as it surpasses the requirement of Hamiltonian simulation and directly uses  the QFT (which takes lesser quantum gates than the Hamiltonian simulation).
For an input qubit size of n q = 10, and matrix dimension of 1024 × 1024, the required quantum gate operations in the circulant approximation method is 199 approximately, whereas the Toeplitz-structured method and the standard Hamiltonian simulation method take 468 and 1400 operations, respectively. Hence, the circulant approximation method achieves more complexity efficiency of approximately 57.48% than the Toeplitz-structured Hamiltonian simulation, and about 85.79% than the standard Hamiltonian simulation, at the cost of approximation error.

D. COMPARISON OF ESTIMATED SPECTRUM IN ARRAY PROCESSING
In Fig. 9, we show the estimated spectrum for a Hermitian Toeplitz matrix using the classical TMRA approach and proposed quantum Algorithm 1. The parameters taken in the simulation are given in Table 2. We take a simulation framework with sensor (antenna) size of M = 16 uniformly distributed with the distance between elements be λ 2 with λ = 3.7 mm considering a transmission frequency of 80 GHz. We have kept the signal power to be 0.5 W for a sinusoidal signal, and the standard deviation of the Gaussian additive noise to be 0.01. We have taken ten-objects in front of the receiver whose angular positions θ can be uniformly distributed within (− π 2 , π 2 ). We took 5000 simulations to take an estimate of the Hermitian Toeplitz matrixT y for the TMRA problem in our setting. We have found that l 2 and Frobenius norm of the estimated matrix is approximately T y 2 = 1.0087, and T y F = 1.99, respectively.
In the quantum simulation of matrixT y , we take QTR = 0.48 with preparing the unitary operator exp(−iT y τ ). For computing the QTR, we have used expression (85) with T y 2 = 1.0087. We have simulated the quantum algorithm on a classical computer due to quantum resource limitations.
In Fig. 9, the (normalized) eigenvalues for matrixT y computed by classical method and the proposed quantum Algorithm 1 is shown. These two methods have an almost similar eigenvalue spectrum. We have seen that the mean absolute error in spectrum estimation between these methods is approximately 3 × 10 −3 .
In Fig. 10, we have shown a comparison plot for the probability error in detecting objects using the classical and quantum algorithms. For the source detection in array processing, we have followed the QA-SDA algorithm as shown in [59] using the estimated Toeplitz matrixT y . We have varied the SNR and kept all other parameters the same as Table 2. The probability of error (P e ) at an SNR of −11 dB for the classical TMRA approach is 0.23 approximately, whereas the P e = 0.08 for the proposed quantum approach. The error probability decreases at a higher SNR for both algorithms.

TABLE 3. Parameters on IBM Simulator
In comparison with the classical TMRA approach, we have seen that the error probability P e due to missed detecting objects for the proposed algorithm is improved for the DoA estimation of multiple objects.

E. SIMULATION OF THE PROPOSED ALGORITHM ON IBM QUANTUM SIMULATOR
We have simulated a small-scale experiment for the eigenvalue estimation of a Toeplitz matrix (T N with N = 16 generated using Table 2) on an IBM real-time quantum simulator platform with QISKIT script. Here, we are using a normalized Hermitian Toeplitz matrix and the highest eigenvalue is 1 (see Fig. 9). Hence, we are computing the second principal eigenvalue and the minimum eigenvalue (which are double numbers) to see the precision in the estimated eigenvalues. In Table 3, we have given the parameters for the quantum simulation on an IBM quantum computer as follows. The circuit implementation for both applications is given in Fig. 12. Here, we have used six control qubits, one target quantum register, and six classical register, and we have performed the simulation 10 5 times to get the histogram. We kept our measurement bases with 6-bit binary representation and obtained our result up to 4-b precision. The histogram of eigenstate probability for the second principal eigenvalue is shown in Fig. 11(a), and the histogram corresponding to the minimum eigenvalue state is shown in Fig. 11(b). We have run the experiment on an IBM real-time Statevector quantum machine for 10 5 times and counted the number of repetitions of every computational basis to get the probability of occurrence. The estimated eigenvalue corresponds to the measurement basis is the one that has the highest probability of occurrence, given bŷ where N dec denotes the decimal value of the binary string of the measurement basis (having the highest probability) and n q is the number of ancillary qubits (here, n q = 6). With the above (see Table 3) parameters and simulation environment, we have obtained that the second estimated eigenvalue is 0.5781 and the estimated minimum eigenvalue is 0.0469. The estimated eigenvalues are relatively closer to the true eigenvalues (which are 0.5748, and 0.0187 respectively). However, with more qubit and quantum resources, the error can be further minimized and the full experiment (computation of the eigenvalue spectrum) may be performed.

VIII. DATA AVAILABILITY FOR THE IMPLEMENTATION OF THE ALGORITHM ON IBM QUANTUM QISKIT
The QEE circuit is implemented on an IBM quantum machine, as shown in Fig. 12. Here, we have used six quantum registers as control registers, one quantum register for eigenstate initialization, and six classical registers to store the result. We have made an open-access QISKIT script for this simulation as available on [60].

IX. CONCLUSION
In this article, we have proposed a Toeplitz-structured quantum Hamiltonian simulation framework having a significant gate-operation complexity advantage for spectrum estimation problems. We have proposed two algorithms, viz., one using a proposed quantum Jordan-gatebased Hamiltonian simulation, and another using a circulant approximation-based QFT approach. The proposed algorithms can help simulate many classical signal processing applications faster where Toeplitz and circulant matrices are relevant. We have devised the bound on the error considering large-dimensional systems with practical assumptions. An application of the proposed framework is shown for TMRA in array signal processing. Looking forward, the proposed low-complex algorithms for structured Hamiltonian are promising avenues to explore shortly for many large-scale applications.