Hybrid Quantum Singular Spectrum Decomposition for Time Series Analysis

A recently proposed approach to classical time series analysis is singular spectrum decomposition (SSD) Bonizzi et al. (2014). SSD aims to decompose non-linear and non-stationary time series into physically meaningful components. In stark contrast to a Fourier decomposition which employs a basis set of harmonic functions, SSD prepares an adaptive basis set of functions, whose elements depend uniquely on the time series. Applications include, but are not limited to, sleep apnea detection from an unprocessed ECG Bonizzi et al. (2015) and the processing and analysis of brain waves Lowet et al. (2017, 2016). In this work we will also cover application to gravitational wave analysis. The latter is interesting as real-time data analysis poses a challenge, and will be vital for next generations of gravitational wave detection technology. The main bottleneck for scaling this algorithm up for larger time series is the singular value decomposition (SVD) that is performed in each subroutine. For general $m\times n$ matrices, the classical SVD algorithm yields a computational complexity of $\mathcal{O}\left(mn\text{min}(m,n)\right)$ Cline and Dhillon (2006).

There exist clues that allude to quantum computing being able to provide a computational speedup for the computationally taxing SVD subroutine. Currently, quantum computing is in the Noisy Intermediate Scale Quantum (NISQ) era Preskill (2018). While some forms of error mitigation are proposed, noise still limits the amount of qubits and the fidelity of gate operations, such that their numbers are too little for full-fledged quantum computers running quantum error correction codes Kandala et al. (2018). Despite these shortcomings, hybrid quantum computers running variational quantum algorithms provide a useful platform for classically taxing calculations, one major example being the variational quantum eigensolver (VQE) for quantum chemistry O’Malley et al. (2016); Peruzzo et al. (2014); Colless et al. (2018). Hybrid algorithms complement the strengths of both processors: classically efficient tasks are performed on a classical computer, while unitary calculations such as quantum state preparation are handed to a quantum computer. Because of the iterative nature of variational algorithms, noise can be mitigated at every step, increasing robustness of the algorithm.

The SVD subroutine can be performed on a quantum computer through the variational quantum singular value decomposition (VQSVD) algorithm Wang et al. (2021). Besides its application to SSD, SVD forms a compelling application of hybrid quantum computing in general, because of its versatile utility. In a similar fashion to VQE, one can train two parameterisable quantum circuits (PQCs) to learn the left and right singular vectors of a matrix, and sample the singular values through Hadamard tests. In this work we explore both a gate-based paradigm Peruzzo et al. (2014), where the parameterisation is facilitated through rotational quantum gates, and a pulse-based paradigm de Keijzer et al. (2022), where we optimise laser pulse profiles on a neutral atom quantum computing platform. The latter is derived from quantum optimal control (QOC) theory, and is thought to be more NISQ efficient.
Our aim is to convert SSD into a variational hybrid quantum-classical algorithm called quantum singular spectrum decomposition (QSSD). Here, we pose the question whether quantum computing can provide an inherent computational speedup for SVD by surpassing the classical complexity. We employ randomised SVD to improve scalibility of the algorithm Halko et al. (2011). The classical theory behind SSD is introduced in Sec. II. The implementation of QSSD on a hybrid quantum computer is elaborated on in Sec. III, using a gate-based approach, while in Sec. IV we discuss a pulse-based approach. Simulation results for three different time series, among which the first observation of a gravitational wave, are presented in Sec. V. We conclude with a summary and discussion about the viability of QSSD in Sec. VI.

The SSD algorithm is based on the singular spectrum analysis (SSA) method, whose basic ideas are presented in Appendix A Bonizzi et al. (2014); Vautard and Ghil (1989). SSA is designed to explore the presence of specific patterns within a time series, while SSD is designed to provide a decomposition into its constituents components. Let $x(n)=\{x_{1},x_{2},\cdots,x_{N-1},x_{N}\}$ be an $N$-dimensional string of numbers, sampled at a uniform sampling frequency $F_{S}$ from some time series. We are interested in extracting (physically) meaningful components from the signal. By embedding the signal in a matrix, correlation functions between pairs of sample points can be extracted through a singular value decomposition of said matrix. Several types of components can be found, among which trends, oscillations and noise for instance. In SSD, such identification is completely automated by focusing on extracting narrow frequency band components.

The time series is first embedded in a trajectory matrix, however, in comparison to SSA, the number of columns is extended to $N$ and the time series is wrapped around, yielding:

It is to be understood that every subscript is $\text{mod }N$. The singular value decomposition (SVD) of such a matrix yields

where $\sigma_{i}$ are the singular values, and $\boldsymbol{u}_{i}$ and $\boldsymbol{v}_{i}$ the left and right singular vectors. As explained above, SSD aims to find physically meaningful components $g_{j}(n)$ through a data-driven and iterative approach. Starting from the original signal $x(n)$, at each iteration, a component $g^{(j)}(n)$ is extracted, such that a residue

remains, where $x(n)\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\scriptsize def}}}}{{=}}}v^{(0)}(n)$.

For iteration $j=1$:

The power spectral density (PSD) of the signal is calculated. The frequency $f_{\text{max}}$ associated with the dominant peak in the PSD is then estimated. If the ratio criterion $\frac{f_{\text{max}}}{F_{S}}<\delta$ is met, a trend is said to characterise the spectral contents of the residual. Ref. Bonizzi et al. (2014) sets $\delta=10^{-3}$. In such a case, the embedding dimension is set to $M=\lfloor\frac{N}{3}\rfloor$, as described in Ref. Vautard et al. (1992). The first component $g^{(1)}(n)$ is then extracted from Hankelisation of the matrix $X_{1}=\sigma_{1}\boldsymbol{u}_{1}\boldsymbol{v}_{1}^{\top}$. Since most of the energy must be contained within the trend component, only 1 singular value is required for its reconstruction. If the ratio criterion is not met, the window length is defined by $M=\lfloor 1.2\frac{F_{S}}{f_{\text{max}}}\rfloor$ Bonizzi et al. (2014).

For iterations $j\geq 2$:

For all iterations $j\geq 2$, the embedding dimension is always set equal to $M=\lfloor 1.2\frac{F_{S}}{f_{\text{max}}}\rfloor$. To extract meaningful components, only principal components with a dominant frequency in the band $[f_{\text{max}}-\delta f,f_{\text{max}}+\delta f]$ are selected. In order to estimate the value of the peak width $\delta f$, the Levenberg–Marquardt (LM) algorithm is used. First, the PSD of the residual is approximated with a sum of three Gaussian functions

where $\vec{A}$ is the vector of amplitudes, $\vec{\sigma}$ is the vector of Gaussian peak widths and $\mu_{i}$ are the peak locations, given by

where $f_{2}$ is the frequency at which the second most dominant peak is situated. Starting the algorithm from the initial state vector components

the algorithm will eventually converge towards optimal $\vec{A}^{\star},\vec{\sigma}^{\star}$. The peak width is then given by $\delta f=2.5\sigma_{1}^{\star}$. Finally, the component $g^{(j)}(n)$ is rescaled with a factor $a$, so that $\tilde{g}^{(j)}(n)=ag^{(j)}(n)$. Such rescaling makes sure the variance of the residue is minimised. Solving for

yields $a=\frac{g^{(j),\top}v^{(j-1)}}{g^{(j),\top}g^{(j)}}$ Bonizzi et al. (2014). The LM algorithm stops when the energy of the remainder has been reduced to a predefined threshold (default at 1%), ensuring proper convergence of SSD. This corresponds to the criterion

so that the reconstructed signal reads

Besides the gate paradigm, variational quantum optimal control (VQOC) de Keijzer et al. (2022) offers an alternative approach to searching through a Hilbert space $\mathfrak{H}$ (see Ref. Koch et al. (2022) for a review on QOC in general). Rather than optimising parameterised gates, we directly optimise the physical controls through which quantum operations are implemented. It is believed that for highly entangled systems, VQOC is able to outperform VQE in terms of accuracy, when resources are comparible. This is a result of the higher expressibility of optimal control, which one can prove by showing that a set of randomly initialised pulses produces unitary operations that resemble the Haar distribution more closely than a set of randomly initialised gates in the hardware efficient ansatz Hubregtsen et al. (2021).

In this paper, we consider application to a neutral atom quantum computer that is realised on a lattice of trapped Rydberg atoms, whose electronic state manifold offers an embedding of qubit states $|0\rangle$ and $|1\rangle$. Consider 2 linear arrays of respectively $q_{m}$ and $q_{n}$ atoms, which are treated as 2-level systems where the $|0\rangle$-state is mapped to a hyperfine electronic ground state, and the $|1\rangle$-state is mapped to a Rydberg state. Qubit states are addressed by monochromatic laser pulses, which we aim to optimise in order to find an optimal pulse that approximates the target state through propagators $U(t)\in\mathfrak{L}(\mathfrak{H}_{q_{m}})$ and $V(t)\in\mathfrak{L}(\mathfrak{H}_{q_{n}})$. For a total evolution time $T$ and $t\in(0,T)$, the circuit unitaries are subject to the propagator Schrödinger equation

Here, $H_{i}$, $i\in\{U,V\}$ are the total Hamiltonians that drive the evolution of the qubit states. On a neutral atom quantum computing platform, the native qubit Hamiltonian reads

where $H_{d}$ is the drift Hamiltonian, which represents ’always-on’ interactions in the system, and $H_{c}[Z(t)]$ is the control Hamiltonian, parameterised by the controls $Z(t)$ we impose on our system. For a linear array of $N$ qubits and $L\in\mathbb{N}$ controls, we find the van der Waals interaction

and the control Hamiltonian

Here, $C_{6}$ characterises the Rydberg-Rydberg interaction strength, $R$ is the interatomic distance between neighbouring qubits, and $Q_{l}$ are control operators associated with our control parameters. In this work, we control the Rabi frequency and laser phase $|\Omega(t)|e^{i\varphi(t)}$, and the detuning $\Delta$ of the $|1\rangle$-state, so that their controls are given by

These controls parameterise both circuits by modeling pulses as a set of $T/\delta t\in\mathbb{N}$ time-discretised control pulses with duration $\delta t$. The physical implementation is a result from smoothing out the piecewise constant equidistant pulses at a minor loss of fidelity. We aim to optimise a new loss function

by adjusting (29) to penalise strong controls for $\lambda>0$. We enforce these constraints by adding time-dependent Lagrange multiplier terms to $\mathcal{L}_{\text{QOC}}$:

where the multipliers are called adjoints, the inner product is given by

and $\Omega_{N}=L^{2}([0,T],\mathfrak{H}_{N})$. To ensure a realistic implementation, we restrict our controls to the space $\mathcal{Z}_{\text{ad}}$ of admissible pulses

In Appendix D, we prove that it is possible to optimise this Lagrangian under these constraints by sampling the necessary quantities on a quantum computer, where we have omitted the proof of existence and uniqueness of the solutions. Our update scheme at iteration $k$ becomes

for the $U$-controls, and idem dito for $V$. For optimisation, we adopt the same Armijo optimiser protocol where we check for the first Wolfe condition at every step. In Appendix D we also explicitly solve for the adjoints. In conclusion, the VQOC paradigm requires to solve the propagator Schrödinger equations (42), (43), the adjoint equations (80) and (81), and the analytical gradient (53).

To benchmark the performance of QSSD, several simulations have been run, using a classical emulator of a noise-free quantum circuit. For efficiency purposes we assume full access to all relevant matrix elements, without resorting to sampling amplitudes through measurements. We compare gate-based and pulse-based results to the retrieved classical SSD components for a variety of applications, though we stress that because of the inherent inefficiency of the algorithm, we have not attempted to make a fair comparison between the two implementations. In order to draw a fair contrast, either implementation needs to be run on equivalent evolution circuits, where both the gate-based and pulse-based have equivalent resources. This alludes to a similar circuit evolution time $T$, equivalent controls such as Bloch sphere rotational control only, and the number of quantum evaluations $\#$QE${}_{\text{VQOC}}$ should scale linearly in the quantum resources compared to the number of quantum evaluations $\#$QE${}_{\text{VQE}}$.

In this work, we adopted the hardware efficient ansatz for the gate-based circuits, with depth $D=10$ and an iteration number it = 1000, and we assumed full control for the pulse-based circuits, in natural units where $\hbar=1$ and $C_{6}=1$, for a total evolution time $T=10$, a step size discretisation $T/\delta t=100$ and an iteration number it = 150. We aim to show that with general resources, one can approximate the function space of the SVD reasonably well on a quantum computer, without employing equivalent evolution circuits. Note that the VQOC parameters are not fine-tuned for any implementation in particular. For a more realistic implementation on a neutral atom quantum computer, the parameters should lie in experimental bounds, and the space $\mathcal{Z}_{\text{ad}}$ should be adjusted accordingly.

We first apply QSSD to a simple toy example to showcase that properties of the SSD algorithm, such as pinpointing non-stationarity of a signal, carry over to the quantum algorithm. Then we apply QSSD to simulated cortical local field potentials (LFPs), electric potentials originating from the extracellular space in brain tissue, as well as GW150914 data, resulting from the first observation of a gravitational wave event. As a figure of merit we introduce the error measure

where $p(\omega)$ is the PSD of the simulated individual band components, and where $p_{\text{SSD}}(\omega)$ is the PSD of the SSD-recovered component. $\Lambda_{\omega}$ is the relevant frequency band of the PSD over which we integrate, and $|\Lambda_{\omega}|$ is its measure. This error focuses more closely on the frequency content of the retrieved components than an $l_{2}$-norm would, and is not cumulative over time due to small displacements.

In this work we have established a unified framework that combines the promising method of singular spectrum decomposition for classical time series analysis with the variational quantum singular value decomposition algorithm into quantum singular spectrum decomposition (QSSD). The singular value decomposition subroutine is designed to be performed on quantum hardware, while the retrieval of the SSD components is done classically. Using analytical gradient optimisers, we circumvent finite difference errors, and we guarantee optimisation at every iteration. QSSD could find interesting applications to near-term gravitational wave analysis because of the variational formulation. Here, we demonstrated this algorithm for providing a more crisp waveform for matched filtering.

QSSD translates a classical algorithm into a hybrid quantum-classical framework. In contrast to the classical algorithm, the quantum formulation does not have efficient access to all matrix elements of $U$ and $V$. Unpacking the dense trajectory matrix into its unitary basis elements and reading out all the columns of $U$ and $V$ are two processes that required us to perform an exponential number of quantum circuit runs, in the number of qubits $q_{m}$. An efficient quantum algorithm would rely on knowing only few parameters that are readily measured on a quantum device, while employing few resources. To mitigate the effects of the exponential scaling laws, we have adopted randomised SVD. While rendering QSSD more economic, the same trick can be applied to the classical algorithm, providing no inherent quantum speedup. An efficient rendition of SSD would circumvent measuring $|u_{i}\rangle$ and $|v_{i}\rangle$, and instead use a different pipeline to read out SSD components. For further research, we pose the question whether this is possible within this variational framework.

The quality of SSD is related to the quality of the singular vectors, not the singular values, as dictated by Eq. (8), since a rescaling is performed after every iteration to ensure optimal energy extraction. It is therefore crucial for the convergence of the quantum routine to select a circuit ansatz that can approximate the solution space well. For the gate-based approach we have adopted the hardware efficient ansatz (HEA), and for the pulse-based approach we have assumed full control. Despite the HEA being often invoked as a standard method for quantum state preparation, there is a lack of understanding on its capabilities. Since data analysis does not know an effective design or model, we pose the question if an effective model can be established for specific time series applications, and if certain ansätze, either gate-based and pulse-based, can improve convergence for finding the orthonormal matrices.

In conclusion, it remains an open problem whether a true quantum speedup can be achieved through the method of QSSD. The variational formulation of this quantum algorithm does not allow for polynomial scaling laws in the number of qubits, and requires measuring exponentially many quantum amplitudes. At the same time, as mentioned above, several aspects of the theory can be improved upon, which are not necessarily restricted to QSSD applications only. For quantum state preparation purposes, it is interesting to investigate if one could apply quantum natural gradient (QNG) theory to the optimiser routine to improve convergence Stokes et al. (2020), or quantify the geometry of the search space to find better initial states Katabarwa et al. (2022). Overall, results from this preliminary attempt to translate SSD to a quantum computing framework show that this should be likely achievable with proper improvements in the theory and methods employed. In turn, this may pave the way to a true quantum speedup of QSSD over classical SSD. We leave the adaptations of all these ideas to future studies of QSSD and alternative approaches to time series analysis on a quantum computer.

We thank Robert de Keijzer, Madhav Mohan and Menica Dibenedetto for discussions. This research is financially supported by the Dutch Ministry of Economic Affairs and Climate Policy (EZK), as part of the Quantum Delta NL programme, and by the Netherlands Organisation for Scientific Research (NWO) under Grant No. 680.92.18.05.

The data that support the findings of this study are available from the corresponding author upon request. The GW150914 gravitational wave event data is publicly available LIGO-Virgo collaboration (2019). Simulations were performed using QuTiP Johansson et al. (2012).