Fast Machine-Precision Spectral Likelihoods for Stationary Time Series (2024)

Christopher J. GeogaDept. of Statistics, University of Wisconsin-Madison

Abstract

We provide in this work an algorithm for approximating a very broad class ofsymmetric Toeplitz matrices to machine precision in $\mathcal{O}(n\log n)$ timewith applications to fitting time series models. In particular, for a symmetricToeplitz matrix $\mathbf{\Sigma}$ with values $\mathbf{\Sigma}_{j,k}=h_{|j-k|}=\int_{-1/2}^{1/2}e^{2\pi i|j-k|\omega}S(\omega%)\mathrm{d}\omega$ where $S(\omega)$ is piecewise smooth, we give an approximation $\mathbf{\mathcal{F}}\mathbf{\Sigma}\mathbf{\mathcal{F}}^{H}\approx\mathbf{D}+%\mathbf{U}\mathbf{V}^{H}$ , where $\mathbf{\mathcal{F}}$ is the DFT matrix, $\mathbf{D}$ isdiagonal, and the matrices $\mathbf{U}$ and $\mathbf{V}$ are in $\mathbb{C}^{n\times r}$ with $r\ll n$ . Studying these matrices in the context of timeseries, we offer a theoretical explanation of this structure and connect it toexisting spectral-domain approximation frameworks. We then give a completediscussion of the numerical method for assembling the approximation anddemonstrate its efficiency for improving Whittle-type likelihood approximations,including dramatic examples where a correction of rank $r=2$ to the standardWhittle approximation increases the accuracy of the log-likelihoodapproximation from $3$ to $14$ digits for a matrix $\mathbf{\Sigma}\in\mathbb{R}^{10^{5}\times 10^{5}}$ . The method and analysis of this work applieswell beyond time series analysis, providing an algorithm for extremely accuratesolutions to linear systems with a wide variety of symmetric Toeplitzmatrices whose entries are generated by a piecewise smooth $S(\omega)$ .The analysis employed here largely depends on asymptotic expansions ofoscillatory integrals, and also provides a new perspective on when existingspectral-domain approximation methods for Gaussian log-likelihoods can beparticularly problematic.

Keywords: Time series, Spectral density, Asymptoticexpansion, Oscillatory integral, Gaussian process

Acknowledgements

The author is grateful to Paul G. Beckman and the two anonymous referees fortheir thoughtful comments, discussion, and proofreading of the manuscript. Healso thanks Lydia Zoells for her careful copyediting work in the revisionstage.

1 Introduction

A mean-zero, stationary, and Gaussian time series is a stochastic process $\left\{Y_{t}\right\}_{t\in\mathbb{Z}}$ such that any contiguous collection of measurements $\bm{y}=[Y_{t_{0}},Y_{t_{0}+1},...,Y_{t_{0}+n-1}]$ with arbitrary $t_{0}\in\mathbb{Z}$ isdistributed as $\bm{y}\sim\mathcal{N}(\bm{0},\bm{\Sigma})$ , where $\bm{\Sigma}_{j,k}=h_{\left|j-k\right|}$ isthe autocovariance matrix generated by the autocovariance sequence $\left\{h_{k}\right\}_{k\in\mathbb{Z}}$ . A very special property of such covariance matrices for stationarytime series is that they are symmetric Toeplitz matrices, meaning thatthey are constant along sub- and super-diagonal bands. A fundamental andimportant problem in time series analysis is fitting a parametric model to data,which in this setting means using some parametric family of models for theautocovariance sequence $\left\{h_{k}(\bm{\theta})\right\}_{k\in\mathbb{Z}}$ that specify thedistribution $\bm{y}\sim\mathcal{N}(\bm{0},\bm{\Sigma}_{\bm{\theta}})$ (like, for example, (AR)(I)(MA)process models). The canonical estimator for parameters $\bm{\theta}\in\bm{\Theta}$ from data is the maximum likelihood estimator (MLE), which is computed byminimizing the negative log-likelihood:

\hat{\bm{\theta}}^{\text{MLE}}=\operatorname*{arg\,min}_{\bm{\theta}\in\bm{%\Theta}}\;\ell(\bm{\theta}\,|\,\bm{y})=\operatorname*{arg\,min}_{\bm{\theta}%\in\bm{\Theta}}\;\frac{1}{2}\left(\log\left|\bm{\Sigma}_{\bm{\theta}}\right|+%\bm{y}^{T}\bm{\Sigma}_{\bm{\theta}}^{-1}\bm{y}\right),

where constant terms in the negative log-likelihood $\ell(\bm{\theta}\,|\,\bm{y})$ are suppressed.

It is an elementary fact that autocovariance sequences are positive(semi-)definite, and by Herglotz’s theorem one has that

{\color[rgb]{0,0,0}\bm{\Sigma}_{j,k}=h_{|j-k|}}=\int_{[-1/2,1/2)}e^{2\pi i|j-k%|\omega}\mathop{}\!\mathrm{d}F(\omega),

where $F$ is the spectral distribution function of $\left\{Y_{t}\right\}$ [9]. If $F$ has a density with respect to the Lebesgue measureso that $\mathop{}\!\mathrm{d}F(\omega)=S(\omega)\mathop{}\!\mathrm{d}\omega$ , we call $S(\omega)$ the spectraldensity function (SDF). An elementary property of spectral densities is thatthey need only be non-negative, and symmetric about the origin in the case ofreal-valued processes $Y_{t}$ . Considering that writing parametric families offunctions that are nonnegative and symmetric is much easier than writingparametric families of sequences that are positive (semi-)definite, it comes asno surprise that modeling $\left\{h_{k}(\bm{\theta})\right\}_{k\in\mathbb{Z}}$ in terms of the spectraldensity $S_{\bm{\theta}}(\omega)$ is an attractive and liberating option forpractitioners.

Unfortunately, there is no particularly easy way to reformulate $\ell(\bm{\theta}\,|\,\bm{y})$ directly in terms of $S(\omega)$ , meaning that computing the true MLE $\hat{\bm{\theta}}^{\text{MLE}}$ requires either knowing the closed-form expression for $\bm{\Sigma}_{j,k}=h_{|j-k|}=\int_{-1/2}^{1/2}e^{2\pi i\omega|j-k|}S(\omega)%\mathop{}\!\mathrm{d}\omega$ for all lags $|j-k|$ or directly computing all the integrals numerically,and in both cases, for direct likelihood evaluation, one must in general build $\bm{\Sigma}$ entry-wise and evaluate $\ell(\bm{\theta}\,|\,\bm{y})$ using dense linear algebra. As wewill discuss, there are several popular likelihood approximations for $\bm{y}$ interms of $S_{\bm{\theta}}(\omega)$ that are convenient and computationally expedient,often running in $\mathcal{O}(n\log n)$ complexity compared to the $\mathcal{O}(n^{3})$ of adirect likelihood approximation (another fundamental challenge that motivatesmany spectral approximations). The most classical and simple of theseapproximations to $\ell(\bm{\theta}\,|\,\bm{y})$ , denoted the Whittle approximation[38], unfortunately yields very biased estimators. Improvementsto this method can mitigate the challenging new sources of bias, but still comeat the cost of additional variability compared to $\hat{\bm{\theta}}^{\text{MLE}}$ . Inthis work, we introduce a new approximation to $\ell(\bm{\theta}\,|\,\bm{y})$ that retainsthe $\mathcal{O}(n\log n)$ complexity but can be made effectively exact to doubleprecision with $14-15$ correct digits by exploiting piecewise smoothness ofspectral densities and special rank structure in the matrix $\mathbb{\mathcal{F}}\bm{\Sigma}\mathbb{\mathcal{F}}^{H}$ ,where $\mathbb{\mathcal{F}}$ is the discrete Fourier transform (DFT) matrix. As a result of thisdesign, we obtain effectively exact evaluations of the log-likelihood $\ell(\bm{\theta}\,|\,\bm{y})$ and its derivatives, as well as the ability to perform uncertaintyquantification on estimates.

1.1 Whittle approximations and existing work

Let

J_{n}(\omega)=\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}e^{-2\pi ij\omega}Y_{j}

denote the discrete Fourier transform of a time series ${\color[rgb]{0,0,0}\bm{y}=}\left\{Y_{0},...,Y_{n-1}\right\}$ at frequency $\omega$ . For the entirety of this work, we willrestrict focus to $\omega$ at Fourier frequencies $\left\{\omega_{j}=\tfrac{(j-n/2-1)}{n}\right\}_{j=1}^{n}$ so that one can compute $[J_{n}(\omega_{1}),...,J_{n}(\omega_{n})]$ in $\mathcal{O}(n\log n)$ time with $\mathbb{\mathcal{F}}\bm{y}$ using an FFT, where $\mathbb{\mathcal{F}}$ is parameterized in the unitary form with the “fftshift” for convenience as $\mathbb{\mathcal{F}}=[n^{-1/2}e^{-2\pi ijk/n}]_{j\in(-n/2):(n/2-1),k\in 0:(n-1)}$ .Motivated by the elementary observations that $\mathbb{E}|J_{n}(\omega_{k})|^{2}\rightarrow_{n}S(\omega_{k})$ and $\text{Cov}(J_{n}(\omega_{k}),J_{n}(\omega_{k^{\prime}}))\rightarrow_{n}0$ (see[9] for an introduction), one can substitute thefinite-sample (co)variances with these asymptotic ones to obtain the Whittleapproximation

(1.1)

2\ell^{W}(\bm{\theta}\,|\,\bm{y})=\sum_{j=1}^{n}\log S_{\bm{\theta}}(\omega_{j%})+\frac{\left|J_{n}(\omega_{j})\right|^{2}}{S_{\bm{\theta}}(\omega_{j})},

which has this simplified form since the asymptotic covariance matrix for $[J_{n}(\omega_{1}),...,J_{n}(\omega_{n})]$ is diagonal. Optimizing this approximationinstead of $\ell(\bm{\theta}\,|\,\bm{y})$ gives the estimator $\hat{\bm{\theta}}^{W}$ . Anothercommon way to think about this approximation is to observe that it effectivelyapproximates $\bm{\Sigma}_{\bm{\theta}}$ with the circulant matrix $\bm{C}$ whose firstcolumn $\bm{c}$ has $(\mathbb{\mathcal{F}}\bm{c})_{j}=S(\omega_{j})$ . This is unfortunatelyincorrect, as $\bm{\Sigma}$ is almost never exactly circulant. Further, evenif it were, it wouldn’t in general for any finite $n$ be that circulantmatrix. As a result of these finite-sample applications of limiting identities,this estimator can exhibit severe bias for finite sample sizes and certainvarieties of spectral densities $S(\omega)$ —see [33] for aparticularly thoughtful and concrete analysis of the sources of bias of thisapproximation in terms of boundary effects. Yet another way to understand the sourceof bias with this approximation is to observe that even if $\bm{y}\sim\mathcal{N}(\bm{0},\bm{\Sigma}_{\bm{\theta}_{0}})$ with $[\bm{\Sigma}_{\bm{\theta}_{0}}]_{j,k}=\int_{-1/2}^{1/2}e^{2\pi i\omega|j-k|}S_%{\bm{\theta}_{0}}(\omega)\mathop{}\!\mathrm{d}\omega$ , $S_{\bm{\theta}_{0}}(\omega_{k})$ is not thevariance of $J_{n}(\omega_{k})$ for any finite $n$ . And while ignoring covariancestructure between values (which the Whittle approximation also does) can in some casesonly come at the cost of efficiency, using incorrect variances will almostalways come at the cost of bias. The exact (co)variances of these DFTcoefficients is a classical computation, but because they will be usedextensively in this work we provide a terse statement of their derivation here.

Proposition 1.1.

Let $\left\{Y_{t}\right\}_{t=0}^{n-1}$ be a stationary mean-zero time series with spectraldensity $S(\omega)$ and $J_{n}(\omega_{k})$ defined as above, and define $S_{n}(\omega_{k},\omega_{k^{\prime}})=\textrm{Cov}(J_{n}(\omega_{k}),J_{n}(%\omega_{k^{\prime}}))$ .Then

(1.2)

S_{n}(\omega_{k},\omega_{k^{\prime}})=\frac{e^{i\frac{n-1}{n}\pi(k-k^{\prime})%}}{n}\int_{-1/2}^{1/2}D^{s}_{n}(\omega_{k}-\omega)D^{s}_{n}(\omega_{k^{\prime}%}-\omega)S(\omega)\mathop{}\!\mathrm{d}\omega,

where $D^{s}_{n}(\omega)=\frac{\sin(\pi n\omega)}{\sin(\pi\omega)}=e^{-i(n-1)\omega/2%}\sum_{j=0}^{n-1}e^{2\pi ij\omega}$ is a “shifted” Dirichlet kernel.

Proof 1.2.

Using Herglotz’s theorem, an elementary computation shows that

	$\displaystyle\text{Cov}(J_{n}(\omega_{k}),J_{n}(\omega_{k^{\prime}}))=\int_{-1%/2}^{1/2}\left(\sum_{j=0}^{n-1}n^{-1/2}e^{2\pi ij(\omega-\omega_{k})}\right)%\left(\sum_{j=0}^{n-1}n^{-1/2}e^{2\pi ij(\omega_{k^{\prime}}-\omega)}\right)S(%\omega)\mathop{}\!\mathrm{d}\omega$
	$\displaystyle=\frac{1}{n}\int_{-1/2}^{1/2}\left(e^{-i(n-1)\pi\omega_{k}}\frac{%\sin(\pi n(\omega_{k}-\omega))}{\sin(\pi(\omega_{k}-\omega))}\right)\left(e^{i%(n-1)\pi\omega_{k^{\prime}}}\frac{\sin(\pi n(\omega_{k^{\prime}}-\omega))}{%\sin(\pi(\omega_{k^{\prime}}-\omega))}\right)S(\omega)\mathop{}\!\mathrm{d}\omega.$

Applying the definition of $D^{s}_{n}(\omega)$ , collecting terms, and bringing thecomplex exponential prefactor on the trigonometric terms outside the integralgives the results.

A particularly elegant proposition for dealing with the bias of $\hat{\bm{\theta}}^{W}$ is given in [35]: instead of using $S_{\bm{\theta}}(\omega_{j})$ in $\ell^{W}$ , use the correct finite-sample marginal variances $S_{n}(\omega_{j},\omega_{j})$ in their place. In doing so, one immediately obtains a new likelihoodapproximation whose gradient gives unbiased estimating equations (UEEs) forparameters $\bm{\theta}$ [35, 22]. This estimator, calledthe debiased Whittle estimator, is given as

\hat{\bm{\theta}}^{DW}=\operatorname*{arg\,min}_{\bm{\theta}\in\bm{\Theta}}%\ell^{DW}(\bm{\theta}\,|\,\bm{y})=\operatorname*{arg\,min}_{\bm{\theta}\in\bm{%\Theta}}\frac{1}{2}\left(\sum_{j=1}^{n}\log S_{n,\bm{\theta}}(\omega_{j},%\omega_{j})+\frac{\left|J_{n}(\omega_{j})\right|^{2}}{S_{n,\bm{\theta}}(\omega%_{j},\omega_{j})}\right),

Where the additional subscript of $\bm{\theta}$ in $S_{n,\bm{\theta}}$ is includedsimply to reinforce the dependence of this function on (varying) parameters $\bm{\theta}$ and is independent from the subscript $n$ , which is used to signify thefinite-sample correction. $\hat{\bm{\theta}}^{DW}$ now only comes at the cost ofadditional variance from ignoring the correlation of $J_{n}(\omega_{k})$ with $J_{n}(\omega_{k^{\prime}})$ , and serves as a useful point of comparison with the estimatorproposed here, as it better isolates the additional variance ofWhittle-type estimators compared to $\hat{\bm{\theta}}^{\text{MLE}}$ .

While the standard Whittle approximation requires a pre-computed FFT of the data,the debiased Whittle approximation requires one FFT per evaluation of $\ell^{DW}(\bm{\theta}\,|\,\bm{y})$ , since one can compute $\left\{S_{n}(\omega_{j},\omega_{j})\right\}_{j=1}^{n}$ in the time domain using the classical identity that

(1.3)

S_{n}(\omega_{j},\omega_{j})=2\text{Re}\left\{\sum_{k=0}^{n-1}(1-n^{-1}h)h_{k}%e^{-2\pi i\omega_{j}h}\right\}-h_{0}.

In [35] it is noted that this identity is useful to avoidnumerical integration. While this is true, it is also quite inconvenientbecause it requires knowing the covariance function $\left\{h_{k}\right\}$ as well asthe spectral density $S(\omega)$ . Considering that a large part of thepoint of resorting to a Whittle-type approximation is to write a parametricspectral density instead of a covariance function, this is a potentially seriousinconvenience. In Section 2, in the process of building up to our ownapproximation, we will discuss two performant and simple numerical integrationstrategies to relax this requirement and also make this method more generallyavailable to practitioners.

1.2 Our method

The approximation we introduce here is a continuation of such spectrallymotivated likelihood methods. Letting $\bm{D}$ denote thediagonal matrix with values $\left\{S(\omega_{j})\right\}_{j=1}^{n}$ , the fundamentalobservation of this work is that for the very broad class ofpiecewise-smooth spectral densities with a small number of roughpoints—including basically all standard parametric families, including(AR)(I)(MA), Matérn, rational quadratic, power law,and other models—we have that

(1.4)

\mathbb{\mathcal{F}}\bm{\Sigma}\mathbb{\mathcal{F}}^{H}\approx\bm{D}+\bm{U}\bm%{V}^{H}=\bm{D}+{\color[rgb]{0,0,0}\bm{E}},

where $\bm{U}$ and $\bm{V}$ are in $\mathbb{C}^{n\times r}$ with $r\ll n$ . Thestandard Whittle approximation views that low-rank update as being exactly zero,and in this light our method can be viewed as a “corrected” Whittleapproximation. The key observation that makes working with this approximationconvenient is that the “Whittle correction” matrix ${\color[rgb]{0,0,0}\bm{E}}=\mathbb{\mathcal{F}}\bm{\Sigma}\mathbb{\mathcal{F}}%^{H}-\bm{D}$ ,along with being severely rank-deficient in many settings, can be applied to avector in $\mathcal{O}(n\log n)$ complexity via circulant embedding of Toeplitzmatrices, and so one can utilize the methods of randomizedlow-rank approximation [21] to assemble $\bm{U}$ and $\bm{V}$ efficiently and scalably with implicit tools.Counter-examples of models where $\mathbb{\mathcal{F}}\bm{\Sigma}\mathbb{\mathcal{F}}^{H}$ does not have thisstructure are rare in practice, as most popular and natural models for spectraldensities $S(\omega)$ are smooth except at possibly a modest number oflocations. The category of models that likely will not be amenable tothis structure, as will be discussed at length in Section3, are those where $S(\omega)$ is either highly oscillatoryor has many rough points. Spectral densities like $S(\omega)=\gamma\text{sinc}(\gamma\omega)^{2}$ or $S(\omega)=\sum_{j=1}^{R}\alpha_{j}e^{-\theta|\omega-\omega_{j}|}$ for large $\gamma$ or $R$ respectively are good examples of modelsfor which the rank $r$ required for high accuracy is sufficiently large that themethod is not practical.

Outside of these pathological cases,despite having the same quasilinear runtime complexity, this representation canbe exact to computer precision for very small $r$ (sometimes as small as $r=2$ ). As one might expect, this does come at the cost of a moreexpensive prefactor. If the standard Whittle approximation requires onepre-computed FFT of the data and subsequently runs at linear complexity, and thedebiased Whittle approximation requires one additional FFT per evaluation due tocomputation of the terms in (1.3), our approximation requires $\mathcal{O}(3r)$ FFTs to assemble (although a simpler implementation using $\mathcal{O}(4r)$ FFTs isused in the computations of this work)¹¹1A software package and scriptsfor all computations done in this work is available athttps://github.com/cgeoga/SpectralEstimators.jl..For a spectral density that is neither particularly well- nor poorly-behaved, areasonable expectation for an $r$ that gives all $14$ significant digits of $\ell(\bm{\theta}\,|\,\bm{y})$ is around $r\approx 100$ —and so for full precision ofthe likelihood, this method can require several hundred FFTs. But consideringhow fast FFTs are, the prefactor cost for this method is sufficiently low thatit compares favorably with other high-accuracy alternatives.

Another material advantage of this method over other spectral approximations isthat it can be used to compute uncertainties for estimators. A basic fact aboutMLEs is that, under sufficient regularity conditions, $\bm{I}(\bm{\theta})^{-1/2}(\hat{\bm{\theta}}^{\text{MLE}}-\bm{\theta}^{\text{%true}})\rightsquigarrow\mathcal{N}(\bm{0},\bm{\mathcal{I}})$ , where $\bm{I}(\bm{\theta})$ is the expected Fisher informationmatrix that will be introduced and discussed at length later in this work,although in some settings it may actually be preferable to use the “observed”information matrix $H\ell(\bm{\theta}\,|\,\bm{y})\,|\,_{\bm{\theta}=\hat{\bm{\theta}}^{\text{MLE}}}$ [14].But when one is obtaining an estimator by optimizing some function $\tilde{\ell}(\bm{\theta})$ other than the true log-likelihood, the Hessian of $\tilde{\ell}$ does not have this interpretation. The fundamental property anapproximation $\tilde{\ell}(\bm{\theta})$ needs to satisfy to provide a statisticallyuseful estimator is that it provides unbiased estimating equations(UEEs), so that $\mathbb{E}_{\bm{\theta}_{0}}[\nabla\tilde{\ell}(\bm{\theta}\,|\,\bm{y})\,|\,_{%\bm{\theta}=\bm{\theta}_{0}}]=\bm{0}$ for any $\bm{\theta}_{0}$ [22]. But there are many such $\tilde{\ell}$ functions that satisfy this property and are not goodpointwise approximations to $\ell$ . And in those cases, one certainly cannotexpect that the Hessian of $\tilde{\ell}$ at $\hat{\bm{\theta}}^{\text{MLE}}$ resemblesthe Hessian of $\ell$ at those values. Approximations like the debiasedWhittle-induced $\tilde{\ell}$ fall into this category. The method of this work,however, can be used to obtain both estimators and approximateuncertainties via its Hessian matrix, since it provides a fully pointwise-accurateapproximation to $\ell(\bm{\theta})$ at for all $\bm{\theta}$ .

There is a significant body of literature on the more general problem ofapproximating covariance matrices for the purpose of estimation problems. Amongthe most accurate varieties of approximations in this space is the use ofhierarchical matrices, which employ rank-structure of matrixsub-blocks in a way that provides high accuracy but retains good run-timecomplexity. These methods have been applied to the Gaussian process problem in avariety of ways[8, 1, 15, 28, 26, 17],and unlike the tools proposed in this work are not specialized toone-dimensional gridded measurements. Similarly, tools specific to Toeplitzmatrices that exploit rank-structure in off-diagonal blocks of $\mathbb{\mathcal{F}}\bm{\Sigma}\mathbb{\mathcal{F}}^{H}$ have been proposed for fast direct linear solvers[27, 10].The downside to these methods is that thealgorithms are often quite complicated, and more domain-specific quantities instatistical computing—such as log-determinants, derivatives with respect tomodel parameters, and certain matrix-matrix products—are more involved tocompute, and are often computed to lower accuracy[3, 32, 28, 17]. Other popular methods inthis space that provide less accurate likelihood approximations in the pointwisesense but are often easier to implement and analyze include sparseapproximations to $\bm{\Sigma}^{-1}$ [30, 25, 29, 34, 24], direct low-rankapproximation methods [11], and matrix tapering [16].The method presented here fits into this body of work as a specialized methodfor stationary time series data. In exchange for that specialization, however,it achieves the crushing accuracy of hierarchical matrix methods in a muchsimpler and more tractable approximation format, and permits the fast andaccurate computation of domain-specific quantities like log-determinantsand derivatives of the Gaussian log-likelihood.

1.3 Outline

In the next several sections, we will introduce tools for implicitly assemblingand working with matrices of the form (1.4). We start with adiscussion of using quadrature and asymptotic expansions to efficiently andaccurately work with $\bm{\Sigma}$ using only the spectral density $S(\omega)$ , and wethen provide a discussion of (i) why the perturbation term in (1.4)should be low-rank, which also offers a new theoretical perspective for whenWhittle approximations are particularly problematic, and (ii) the details ofassembling the low-rank $\bm{U}\bm{V}^{H}$ entirely from fast matrix-vectorproducts. We then close by providing several tests and benchmarks to demonstratethe speed and accuracy of the method.

2 Oscillatory integration of spectral densities

As discussed in the previous section, creating the low-rank approximation $\mathbb{\mathcal{F}}\bm{\Sigma}\mathbb{\mathcal{F}}^{H}\approx\bm{D}+\bm{U}\bm%{V}^{H}$ requires working with $\bm{\Sigma}$ , theToeplitz matrix with values $\bm{\Sigma}_{j,k}=h_{|j-k|}$ where $\left\{h_{k}\right\}$ is theautocovariance sequence of the time series. Unlike in [35] whereboth kernel and spectral density values are used to avoid numerical integration,we will now discuss strategies for efficiently and accurately computingintegrals of the form $h_{k}=\int_{-1/2}^{1/2}S(\omega)e^{2\pi ik\omega}\mathop{}\!\mathrm{d}\omega$ for $k\in 0,...,n-1$ , so that one can create this factorization given onlythe spectral density. Let us briefly review the challenges of such a task.First, the standard method for obtaining the values $\left\{h_{k}\right\}$ using thetrapezoidal rule via the FFT may require a large number of nodes for evenmoderate accuracy if, for example, $S(\omega)$ has rough points (like how $S(\omega)=e^{-|\omega|}$ has no derivatives at the origin). A second common issue arisesif $S(\omega)$ is not periodic at its endpoints, which limits the convergence rateof the trapezoidal rule to its baseline rate of $\mathcal{O}(m^{-2})$ for $m$ nodes forgeneral non-smooth and/or non-periodic integrands. As a final complication, bythe Shannon-Nyquist sampling theorem, one requires at least $\mathcal{O}(k)$ quadraturenodes to resolve the oscillations of $e^{2\pi ik\omega}$ on $[-1/2,1/2]$ , andso any direct summation method to obtain $\left\{h_{k}\right\}_{k=0}^{n-1}$ would scale atleast as $\mathcal{O}(n^{2})$ .

To bypass these issues, we propose the use of adaptive quadraturefor a fixed number of lags $k$ , but then a transition to the use of asymptoticexpansions to $h_{k}$ for lags greater than some $k_{0}$ . For integrands that arenot too oscillatory, standard adapative integration tools can be veryefficient and capable for any piecewise-smooth integrand. In this work, theJulia lanuage package QuadGK.jl [7, 23] was used,although many other tools would likely have worked equally well. Briefly,Gaussian quadrature is a tool for achieving high-order accuracy in the numericalintegration of smooth functions. It gives approximations to integrals as $\int_{a}^{b}f(x)\mathop{}\!\mathrm{d}x\approx\sum_{j=1}^{M}\alpha_{j}f(x_{j})$ that are exact for polynomials up to order $2M-1$ , where $\left\{\alpha_{j}\right\}_{j=1}^{M}$ are weights and $\left\{x_{j}\right\}_{j=1}^{M}$ arenodes. In the case of $[a,b]=[-1,1]$ , the Legendre polynomials can beused to obtain the weights and nodes and achieve an error rate of $\mathcal{O}(M^{-2m-1})$ for functions $f\in\mathcal{C}^{(m)}([a,b])$ [18, 36].

More interesting, however, is the discussion of methods for accurately andefficiently computing $h_{k}$ when $k$ is large. Unlike traditional integrationtools, asymptotic expansion-type methods get more accurate as $k$ increases. As in the introduction, the following result is not new (see[12] for a comprehensive discussion), but we state it in thespecific form that is useful for this work since it will be referred tofrequently. Since the idea of the proof is also used repeatedly, we againprovide it here.

Proposition 2.1.

Let $S(\omega)\in\mathcal{C}^{(m)}([-1/2,1/2])$ and $\int_{-1/2}^{1/2}|S^{(m)}|(\omega)\mathop{}\!\mathrm{d}\omega<C<\infty$ . Then

Proof 2.2.

Simply apply integration by parts as many times as possible:

	$\displaystyle\int_{-1/2}^{1/2}S(\omega)e^{2\pi ik\omega}\mathop{}\!\mathrm{d}\omega$	$\displaystyle=\frac{S(1/2)e^{\pi ik}-S(-1/2)e^{-\pi ik}}{-2\pi ik}+(-2\pi ik)^%{-1}\int_{-1/2}^{1/2}S^{\prime}(\omega)e^{2\pi ik\omega}\mathop{}\!\mathrm{d}\omega$
		$\displaystyle\vdots$
		$\displaystyle=\sum_{j=0}^{m-1}(-2\pi ik)^{-(j+1)}(S^{(j)}(1/2)e^{\pi ik}-S^{(j%)}(-1/2)e^{-\pi ik})$
		$\displaystyle+(-2\pi ik)^{-(m+1)}\int_{-1/2}^{1/2}S^{(m)}(\omega)e^{2\pi i%\omega k}\mathop{}\!\mathrm{d}\omega.$

Remarkably, evaluating this expansion for any lag $k$ only requires a fewderivatives of $S(\omega)$ at its endpoints. For this reason, with this toolevaluating the tail of $\left\{h_{k}\right\}$ for high lags $k$ actually becomes thefastest part of the domain to handle. And while the supposition of the abovetheorem that $S\in\mathcal{C}^{(m)}([-1/2,1/2])$ is obviously restrictive, asimple observation extends this result to a much broader class of functions.

Corollary 2.3.

Let $S(\omega)\in\mathcal{C}^{(m)}([-1/2,1/2]\setminus\left\{\omega_{1}^{r},...,%\omega_{L}^{r}\right\})$ ,so that it is smooth except at locations $\left\{\omega_{1}^{r},...,\omega_{L}^{r}\right\}$ ,further assume that at each $\omega_{l}^{r}$ it has at least $m$ directionalderivatives $S^{(m\pm)}(\omega_{l}^{r})$ from both the left ( $m-$ ) and the right( $m+$ ) for $l\in 1,...,L$ . Then

	$\displaystyle\int_{-1/2}^{1/2}S(\omega)e^{2\pi ik\omega}\mathop{}\!\mathrm{d}\omega$	$\displaystyle=\sum_{l=1}^{L-1}\sum_{j=0}^{m-1}(-2\pi ik)^{-(j+1)}(S^{(j-)}(%\omega_{l+1}^{r})e^{2\pi ik\omega_{l+1}^{r}}-S^{(j+)}(\omega_{l}^{r})e^{2\pi ik%\omega_{l}^{r}})$
		$\displaystyle+\mathcal{O}(k^{-m-1}).$

This corollary, which is proven by simply breaking up the domain $[-1/2,1/2]$ into segments with endpoints at each $\omega_{l}^{r}$ and applying the above argumentto each segment, now means that this asymptotic expansion approach can be usedto accurately calculate $h_{k}$ for very large $k$ for any spectral density thatis just piecewise smooth. The prefactor on the error term in this settingis natually increased compared to the setting of Proposition 2.1.Nonetheless, however, the convergence rates are on our side. For example, if $k=3\,000$ and one uses five derivatives, then the error term in the simple caseis given by $C\cdot(6,000\pi)^{-6}\approx C\cdot 10^{-26}$ for someconstant $C$ . For a spectral density like $S(\omega)=10e^{-10|\omega|}$ , anexample that will be studied extensively in later sections due to beingwell-behaved except at the origin, we see that $h_{3000}\approx 5\cdot 10^{-7}$ . For that particular function, then, this bound indicates that unless $C$ is quite large, one can reasonably expect a full $14-15$ digits of accuracy.For functions with few rough points, this evaluation method onlyrequires derivatives of $S(\omega)$ at a few locations, and can be evaluated injust a few nanoseconds on a modern CPU once those have been pre-computed.Combined with adaptive quadrature for sufficiently low lags $k$ that theexpansion is not yet accurate, this gives $\left\{h_{k}\right\}_{k=0}^{n-1}$ for all lagsquickly and accurately.

An alternative way to compute $\left\{h_{k}\right\}_{k=1}^{n}$ for potentially large $n$ in $\mathcal{O}(n\log n)$ time is to use nonuniform fast Fourier transforms (NUFFTs) andGauss-Legendre quadrature to evaluate $h_{k}$ at all lags[13, 5]. While the traditional FFT algorithm criticallydepends on measurements being given on a regular grid and exploits thesymmetries that that creates, NUFFTs are fast algorithms to evaluate sums ofcomplex exponentials at potentially irregular locations and target frequenciesin quasilinear runtime. Per the discussion above, we see that computing $\left\{h_{k}\right\}_{k=1}^{n}$ would require $M=\mathcal{O}(n)$ quadrature nodes (to resolve thehighest frequencies) and naturally would need to be evaluated at $n$ target values,which means that naive computation of those integrals would require $\mathcal{O}(n^{2})$ work. With the NUFFT, however, the sums can be computed as a matrix-vectorproduct

h_{k}=\sum_{j=1}^{M}\alpha_{j}e^{2\pi ik\omega_{j}^{q}}S(\omega_{j}^{q}),\quadk%=0,...,n-1

in $\mathcal{O}(n\log n)$ time, where $\omega_{j}^{q}$ denotes the $j$ -th quadrature node.This method, while much more expensive than using asymptotic expansions, has theadvantage that it does not require computing any derivatives to evaluate thesequence $\left\{h_{k}\right\}_{k=0}^{n-1}$ . Because of this, it can be more accurate for SDFs $S(\omega)$ with higher-order derivatives that potentially get very large or arefor some other reason numerically unstable to compute. With that said, however,as discussed above, the accuracy of Gauss-Legendre quadrature depends heavily onhow smooth the integrand is, so we still recommend splitting the domain $[-1/2,1/2)$ into panels such that the rough points of $S(\omega)$ are endpoints as withthe splitting of the asymptotic expansions. If $S(\omega)$ is amenable to thesplit asymptotic expansion method and $n$ is large, then the asymptoticexpansion approach will be significantly faster. But if $S(\omega)$ has some properties that make working with higher-order derivatives infeasible,this NUFFT-based method is an attractive alternative that still enablespractitioners to work exclusively with $S(\omega)$ for evaluating thelog-likelihood and retain $\mathcal{O}(n\log n)$ complexity. For a much more detaileddiscussion of numerical Fourier integration of spectral densities and thequadrature-based tools introduced here for doing it quickly, we refer readers to[6].

3 Low rank structure of the Whittle correction

We now turn to the question of why one should expect the Whittlecorrection matrix ${\color[rgb]{0,0,0}\bm{E}}=\mathbb{\mathcal{F}}\bm{\Sigma}\mathbb{\mathcal{F}}%^{H}-\bm{D}$ where $\bm{D}$ is diagonal with $\bm{D}_{j,j}=S(\omega_{j})$ , to be of lownumerical rank. The primary tool for this investigation will again be asymptoticexpansions, and the crucial observation is that for every $k$ and $k^{\prime}$ , ${\color[rgb]{0,0,0}\bm{E}}_{k,k^{\prime}}$ is given by an oscillatory integral with the same high frequencyof $2\pi n$ . While the below analysis is not actually how we choose to computeand assemble this low-rank approximation in this work, the following derivationprovides some idea for why the Whittle correction matrix ${\color[rgb]{0,0,0}\bm{E}}$ has low-rankstructure, even when $S$ has non-smooth points or $\left\{h_{k}\right\}$ decays slowly. For theduration of this section, we will assume that $S(\omega)\in\mathcal{C}^{(m)}([-1/2,1/2])$ for convenience, although the results here canagain be extended using directional derivatives if $S$ is not differentiable butis smooth from the left and the right at rough points.

We begin by adding and subtracting $S(\omega_{k})$ and $S(\omega_{k^{\prime}})$ to the inner integrandin (1.2), which with minor simplification steps can be expressed as

(3.1)		$\displaystyle\frac{1}{n}$	$\displaystyle\int_{-1/2}^{1/2}D_{n}^{s}(\omega_{k}-\omega)D_{n}^{s}(\omega_{k^%{\prime}}-\omega)S(\omega)\mathop{}\!\mathrm{d}\omega$
	$\displaystyle=\frac{1}{2n}$	$\displaystyle\int_{-1/2}^{1/2}D_{n}^{s}(\omega_{k}-\omega)D_{n}^{s}(\omega_{k^%{\prime}}-\omega)(S(\omega)-S(\omega_{k}))\mathop{}\!\mathrm{d}\omega$
	$\displaystyle+\frac{1}{2n}$	$\displaystyle\int_{-1/2}^{1/2}D_{n}^{s}(\omega_{k}-\omega)D_{n}^{s}(\omega_{k^%{\prime}}-\omega)(S(\omega)-S(\omega_{k^{\prime}}))\mathop{}\!\mathrm{d}\omega$
	$\displaystyle+\frac{S(\omega_{k})+S(\omega_{k^{\prime}})}{2n}$	$\displaystyle\int_{-1/2}^{1/2}D_{n}^{s}(\omega_{k}-\omega)D_{n}^{s}(\omega_{k^%{\prime}}-\omega)\mathop{}\!\mathrm{d}\omega.$

While unwieldy, this representation of the integral already provides areasonably direct explanation about several features of $\mathbb{\mathcal{F}}\bm{\Sigma}\mathbb{\mathcal{F}}^{H}$ .Recalling that

\int_{-1/2}^{1/2}D_{n}^{s}(\omega_{k}-\omega)D_{n}^{s}(\omega_{k^{\prime}}-%\omega)\mathop{}\!\mathrm{d}\omega=\begin{cases}n&k=k^{\prime}\\0&k\neq k^{\prime}\end{cases},

we see that the third term is precisely the diagonal contribution of $S(\omega_{k})$ . We will now argue that the first two terms in the above sumcorrespond to severely rank-deficient matrices. Since they can of course beanalyzed in the exact same way, we study only the first in detail here.

To begin, we slightly rewrite the first term in (3.1) as

(3.2)

\int_{-1/2}^{1/2}\sin(\pi n(\omega_{k}-\omega))\sin(\pi n(\omega_{k^{\prime}}-%\omega))\csc(\pi(\omega_{k}-\omega))\csc(\pi(\omega_{k^{\prime}}-\omega))(S(%\omega)-S(\omega_{k}))\mathop{}\!\mathrm{d}\omega.

From here, we partition the domain into

(3.3)

[-1/2,1/2]=\underbrace{[-1/2,\omega_{k}-\gamma]}_{\text{Type I}}\cup%\underbrace{B_{\gamma}(\omega_{k})}_{\text{Type II}}\cup\underbrace{[\omega_{k%}+\gamma,\omega_{k^{\prime}}-\gamma]}_{\text{Type I}}\cup\underbrace{B_{\gamma%}(\omega_{k^{\prime}})}_{\text{Type II}}\cup\underbrace{[\omega_{k^{\prime}}+%\gamma,1/2]}_{\text{Type I}},

where $B_{\gamma}(x)=[x-\gamma,x+\gamma]$ and $\gamma$ is some smallnumber chosen to keep distance from the singularities of the cosecant terms. Bydesign, then, in Type I intervals the cosecant terms are simple analyticfunctions. Defining

(3.4)

\tilde{S}_{n,k,k^{\prime}}(\omega)=\csc(\pi(\omega_{k}-\omega))\csc(\pi(\omega%_{k^{\prime}}-\omega))(S(\omega)-S(\omega_{k})),

we see that in Type I regions this function is bounded above (assuming that $S$ itself is) and as smooth as $S$ . This motivates the following result that willbe used many times in this section.

Proposition 3.1.

If $g(\omega)\in\mathcal{C}^{(m)}([a,b])$ and $g^{(m)}$ is integrable on $[a,b]$ ,then

(3.5)		$\displaystyle\int_{a}^{b}\sin(\pi n(\omega_{k}-\omega))\sin(\pi n(\omega_{k^{%\prime}}-\omega))g(\omega)\mathop{}\!\mathrm{d}\omega=$
(3.6)		$\displaystyle\frac{(-1)^{k-k^{\prime}}}{2}\int_{a}^{b}g(\omega)\mathop{}\!%\mathrm{d}\omega-\frac{1}{2}\Re\left\{e^{i\pi n\frac{k+k^{\prime}}{2}}\left[e^%{2\pi ib}\sum_{j=0}^{m-1}\frac{g^{(j)}(b)}{(-i2\pi n)^{j+1}}-e^{2\pi ia}\sum_{%j=0}^{m-1}\frac{g^{(j)}(a)}{(-i2\pi n)^{j+1}}\right]\right\}$
	$\displaystyle+\mathcal{O}(n^{-m-1}),$

where $\Re(a+ib)=a$ denotes the real part of a complex number.

Proof 3.2.

Using the product-to-sum angle formula that $\sin(\theta)\sin(\phi)=\tfrac{1}{2}(\cos(\theta-\phi)-\cos(\theta+\phi))$ and standardmanipulations of the complex exponential, the left hand side equation can bere-written as

\frac{1}{2}\cos\left(\pi n\frac{k-k^{\prime}}{n}\right)\int_{a}^{b}g(\omega)%\mathop{}\!\mathrm{d}\omega+\frac{1}{2}\Re\left\{e^{i\pi n\frac{k+k^{\prime}}{%2}}\int_{a}^{b}e^{-i2\pi n\omega}g(\omega)\mathop{}\!\mathrm{d}\omega\right\}.

But $\cos(\pi n(k-k^{\prime}))=(-1)^{k-k^{\prime}}$ , which provides the simplification ofthe first term. For the second term, just as in the proof of Proposition2.1, we simply do integration by parts as many times as possibleon the now standard-form oscillatory integral in the second term to get

\int_{a}^{b}e^{-2\pi in\omega}g(\omega)\mathop{}\!\mathrm{d}\omega=-\sum_{j=0}%^{m-1}(-2\pi in)^{-j-1}\left\{g^{(j)}(b)e^{-2\pi inb}-g^{(j)}(a)e^{-2\pi ina}%\right\}+\mathcal{O}(n^{-m-1}).

An elementary rearrangement gives the result.

While this proposition is presented in generality, there is a reasonable amountof additional simplification that one can do with the right-hand side based onthe value of $k+k^{\prime}\pmod{4}$ . In particular, we see that $e^{i\pi n(k+k^{\prime})/2}=i^{k-k^{\prime}\pmod{4}}$ , which means that every single complexexponential on the right-hand side of (3.5) simplifies nicely. Ifwe assume that $k=k^{\prime}\pmod{4}$ so that $e^{i\pi n(k+k^{\prime})/2}=1$ , theasymptotic expansion term can be more concretely expanded to

		$\displaystyle\Re\left\{(\cos(2\pi b)+i\sin(2\pi b))\sum_{j=0}^{m-1}\frac{g^{(j%)}(b)}{(-i2\pi n)^{j+1}}-(\cos(2\pi a)+i\sin(2\pi a))\sum_{j=0}^{m-1}\frac{g^{%(j)}(a)}{(-i2\pi n)^{j+1}}\right\}$
(3.7)		$\displaystyle=$	$\displaystyle\cos(2\pi b)\sum_{j=0,\,j\text{ odd}}^{m-1}(-1)^{1+\lfloor j/2%\rfloor}\frac{g^{(j)}(b)}{(2\pi n)^{j+1}}-\cos(2\pi a)\sum_{j=0,\,j\text{ odd}%}^{m-1}(-1)^{1+\lfloor j/2\rfloor}\frac{g^{(j)}(a)}{(2\pi n)^{j+1}}$
	$\displaystyle+$	$\displaystyle\sin(2\pi b)\sum_{j=0,\,j\text{ even}}^{m-1}(-1)^{2+\lfloor j/2%\rfloor}\frac{g^{(j)}(b)}{(2\pi n)^{j+1}}-\sin(2\pi a)\sum_{j=0,\,j\text{ even%}}^{m-1}(-1)^{2+\lfloor j/2\rfloor}\frac{g^{(j)}(a)}{(2\pi n)^{j+1}}.$

And while we don’t enumerate the other cases for the three values of $e^{i\pi n(k+k^{\prime})/2}$ , the only difference is in the odd/even indexing and the alternatingsign inside the sum. The key takeaway here is that the entire integral(3.2) can actually be written as a very smooth integral term( $\int_{a}^{b}g(\omega)\mathop{}\!\mathrm{d}\omega$ ), four trigonometric functions that arenon-oscillatory since they have unit wavelength and $a,b\in[-1/2,1/2]$ withcoefficients taken from functions that are themselves smooth (since $g\in\mathcal{C}^{(m)}([a,b])$ ), and a remainder term depending on high-orderderivatives of $g(\omega)$ on $[a,b]$ . As it pertains to the Whittle correctionmatrix ${\color[rgb]{0,0,0}\bm{E}}$ , for each of the Type-I regions we have endpoints $(a=-1/2,b=\omega_{k}-\gamma)$ , $(a=\omega_{k}+\gamma,b=\omega_{k^{\prime}}-\gamma)$ , and $(a=\omega_{k^{\prime}}+\gamma,b=1/2)$ respectively, and so the effective contribution ofthe Type-I integrals to ${\color[rgb]{0,0,0}\bm{E}}_{k,k^{\prime}}$ is given by $\int_{[-1/2,1/2]\setminus(B_{\gamma}(\omega_{k})\cup B_{\gamma}(\omega_{k^{%\prime}}))}\tilde{S}_{n,k,k^{\prime}}(\omega)\mathop{}\!\mathrm{d}\omega$ and a total of $12$ trigonometric functions at the three sets of given endpoints $a$ and $b$ with smoothly varying coefficients. Given that the trigonometricfunctions are analytic and that $g\in\mathcal{C}^{(m)}([-1/2,1/2])$ , sincevery smooth kernel matrices often exhibit rapid spectral decay (an observationdating back to at least [20] with the Fast Multipole Method(FMM)), we expect the contribution of the Type-I intervals to ${\color[rgb]{0,0,0}\bm{E}}$ to haveexceptionally fast spectral decay and be severely rank-deficient so long as theremainder from the asymptotic expansion is small. With slightly more effort, wemay repeat this analysis on the Type-II regions with singularities.

For the first Type II region of $[\omega_{k}-\gamma,\omega_{k}+\gamma]$ , the story isonly slightly more complicated. This time we introduce the condensed notation of

\tilde{S}_{n,k^{\prime}}(\omega)=\csc(\pi(\omega_{k^{\prime}}-\omega))(S(%\omega)-S(\omega_{k}))

for the bounded and non-oscillatory part of the integrand, and we study

(3.8)

\int_{\omega_{k}+\gamma}^{\omega_{k}-\gamma}\sin(\pi n(\omega_{k}-\omega))\sin%(\pi n(\omega_{k^{\prime}}-\omega))\csc(\pi(\omega_{k}-\omega))\tilde{S}_{n,k^%{\prime}}(\omega)\mathop{}\!\mathrm{d}\omega.

The important observation to make here is that the singularity presented by $\csc(\pi(\omega_{k}-\omega))$ is simple. Recalling the Laurent expansion $\csc(t)=t^{-1}+\tfrac{t}{6}+\tfrac{7t^{2}}{360}+...$ , we see that one can simply“subtract off” the singularity and obtain a standard power-series typerepresentation of $\csc(t)-t^{-1}\approx P_{l}(t)$ for $t\in[-\gamma,\gamma]$ and some low-order polynomial $P_{l}$ based on the Laurent series. Thismotivates the decomposition of (3.8) into

		$\displaystyle\int_{\omega_{k}+\gamma}^{\omega_{k}-\gamma}\sin(\pi n(\omega_{k}%-\omega))\sin(\pi n(\omega_{k^{\prime}}-\omega))P_{l}(\omega_{k}-\omega)\tilde%{S}_{n,k^{\prime}}(\omega)\mathop{}\!\mathrm{d}\omega$
	$\displaystyle+$	$\displaystyle\int_{\omega_{k}+\gamma}^{\omega_{k}-\gamma}\sin(\pi n(\omega_{k}%-\omega))\sin(\pi n(\omega_{k^{\prime}}-\omega))\frac{\tilde{S}_{n,k^{\prime}}%(\omega)}{\omega_{k}-\omega}\mathop{}\!\mathrm{d}\omega.$

The first term above can again be expanded as a small number of unit-wavelengthtrigonometric functions via Proposition 3.1 with $g(\omega)=P_{l}(\omega_{k}-\omega)\tilde{S}_{n,k^{\prime}}(\omega)$ . Because $\tilde{S}_{n,k^{\prime}}(\omega)=0$ at $\omega=\omega_{k}$ and has been assumed to be smooth, the second term is also notsingular and is a standard oscillatory integral. Another application ofProposition 3.1 with $g(\omega)=(\omega_{k}-\omega)^{-1}\tilde{S}_{n,k^{\prime}}(\omega)$ can be applied, making the entire contribution of(3.8) expressible as a linear combination of a small numberof unit-wavelength trigonometric functions with smoothly varying coefficients.

The final segment of $[-1/2,1/2]$ to study is the Type II region $[\omega_{k^{\prime}}-\gamma,\omega_{k^{\prime}}+\gamma]$ , where the singularity caused by $\csc(\pi(\omega_{k^{\prime}}-\omega))$ has not already been factored out. Introducing one finalcondensed integrand notation of

\tilde{S}_{n,k}(\omega)=\csc(\pi(\omega_{k}-\omega))(S(\omega)-S(\omega_{k})),

we again decompose the contribution of that segment as

		$\displaystyle\int_{\omega_{k^{\prime}}+\gamma}^{\omega_{k^{\prime}}-\gamma}%\sin(\pi n(\omega_{k}-\omega))\sin(\pi n(\omega_{k^{\prime}}-\omega))\csc(\pi(%\omega_{k^{\prime}}-\omega))\tilde{S}_{n,k}(\omega)\mathop{}\!\mathrm{d}\omega$
	$\displaystyle=$	$\displaystyle\int_{\omega_{k^{\prime}}+\gamma}^{\omega_{k^{\prime}}-\gamma}%\sin(\pi n(\omega_{k}-\omega))\sin(\pi n(\omega_{k^{\prime}}-\omega))P_{l}(%\omega_{k^{\prime}}-\omega)\tilde{S}_{n,k}(\omega)\mathop{}\!\mathrm{d}\omega$
	$\displaystyle+$	$\displaystyle\int_{\omega_{k^{\prime}}+\gamma}^{\omega_{k^{\prime}}-\gamma}%\sin(\pi n(\omega_{k}-\omega))\sin(\pi n(\omega_{k^{\prime}}-\omega))\frac{%\tilde{S}_{n,k}(\omega)}{\omega_{k^{\prime}}-\omega}\mathop{}\!\mathrm{d}\omega.$

The first term in the divided integral can be handled yet again by Proposition3.1 just as above. The only new wrinkle is that the second termin the right-hand side now still has a singularity. Thankfully, with one moreintegral splitting we can re-express this last term with

		$\displaystyle\int_{\omega_{k^{\prime}}+\gamma}^{\omega_{k^{\prime}}-\gamma}%\sin(\pi n(\omega_{k}-\omega))\sin(\pi n(\omega_{k^{\prime}}-\omega))\frac{%\tilde{S}_{n,k}(\omega)}{\omega_{k^{\prime}}-\omega}\mathop{}\!\mathrm{d}\omega$
	$\displaystyle=$	$\displaystyle\int_{\omega_{k^{\prime}}+\gamma}^{\omega_{k^{\prime}}-\gamma}%\sin(\pi n(\omega_{k}-\omega))\sin(\pi n(\omega_{k^{\prime}}-\omega))\frac{%\tilde{S}_{n,k}(\omega)-\tilde{S}_{n,k}(\omega_{k^{\prime}})}{\omega_{k^{%\prime}}-\omega}\mathop{}\!\mathrm{d}\omega$
	$\displaystyle+\tilde{S}_{n,k}(\omega_{k^{\prime}})$	$\displaystyle\int_{\omega_{k^{\prime}}+\gamma}^{\omega_{k^{\prime}}-\gamma}%\sin(\pi n(\omega_{k}-\omega))\sin(\pi n(\omega_{k^{\prime}}-\omega))\frac{1}{%\omega_{k^{\prime}}-\omega}\mathop{}\!\mathrm{d}\omega.$

The first term on the right-hand side is no longer singular because $\tilde{S}_{n,k}(\omega)-\tilde{S}_{n,k}(\omega_{k^{\prime}})$ is smooth and goes to zeroat $\omega=\omega_{k^{\prime}}$ faster than $(\omega_{k^{\prime}}-\omega)^{-1}$ blows up, and so onemay apply Proposition 3.1. The second integral is the integralof $\sin$ against a bounded function since $\sin(\pi n(\omega_{k^{\prime}}-\omega_{k}))$ isalso smooth and goes to zero as $\omega\rightarrow\omega_{k^{\prime}}$ , and thus canbe handled with the even more simple asymptotic expansion analagous to the oneused in Proposition 2.1 that again represents it as a small numberof smooth functions.

Combining all of these segmented contributions, we finally conclude inrepresenting ${\color[rgb]{0,0,0}\bm{E}}_{k,k^{\prime}}$ as a small number of very smooth functions. And whileostensibly breaking the integral into five pieces, two of which need to bebroken up again into two or three additional terms to be suitable for expansion,should lead to so many smooth terms that the matrix is not particularlyrank-deficient, we note that all of these trigonometric functions have unitwavelength and are potentially linearly dependent with each other, and thatderivatives of the corresponding function $g$ in each application of Proposition3.1 are by supposition also well-behaved.

This derivation is also illuminating about the sources of error in Whittleapproximations and the degree of rank-structure in the Whittle correctionmatrices. For example, if $S$ has rough points where it is continuous but notsmooth, then the above domain partitioning will have to be refined so thatCorollary $1$ can be applied to fix the accuracy of the asymptoticapproximation. This will then mean that ${\color[rgb]{0,0,0}\bm{E}}_{k,k^{\prime}}$ will need more asymptoticexpansion-sourced terms for the expansion to be accurate. For spectraldensities $S(\omega)$ that have many such rough points, the number of unique termsin the split expansion can certainly add up to a degree that ruins this rankstructure. And if $S(\omega)$ and its derivatives are highly oscillatory, thenterms like (3) may not actually yield particularlylow-rank matrices despite being smooth.

4 Numerical assembly and action of the Whittle correction

With this theoretical justification for the low-rank structure of the Whittlecorrection matrix, we now turn to discussing how to actually compute andassemble the matrices $\bm{U}$ and $\bm{V}$ such that $\mathbb{\mathcal{F}}\bm{\Sigma}\mathbb{\mathcal{F}}^{H}-\bm{D}\approx\bm{U}\bm%{V}^{H}$ . First, we note that the action of $\bm{\Sigma}$ onto vectorscan be done in $\mathcal{O}(n\log n)$ with the familiar circulant embedding technique[39] once the values $\left\{h_{k}\right\}_{k=0}^{n-1}$ have been computed.While we refer readers to the references for a more complete discussion ofcirculant embedding, we will provide a brief discussion here. The symmetricToeplitz matrix $\bm{\Sigma}$ is specified entirely by its first column $[h_{0},h_{1},...,h_{n-1}]$ . This column can be extended to one of length $2n-1$ given by $\bm{c}=[h_{0},h_{1},...,h_{n-1},h_{n-1},h_{n-2},...,h_{1}]$ , and if oneassembles a Toeplitz matrix with this length $2n-1$ column, the matrix is infact circulant and diagonalized by the FFT [19]. Sincethe FFT can be applied in $\mathcal{O}(n\log n)$ complexity, so can thisaugmented circulant matrix. Thus, one can compute $\bm{\Sigma}\bm{v}$ for any vector $\bm{v}$ by extracting the first $n$ components of $\bm{C}[\bm{v},\bm{0}_{n-1}]$ , where $\bm{C}$ is the circulant matrix made with $\bm{c}$ , and $\bm{0}_{n-1}$ is $n-1$ many zeros appended to $\bm{v}$ .

Letting $\mathbb{\mathcal{F}}_{n}$ be the FFT of size $n$ (briefly now defined without the“fftshift”), we can thus summarize the action of $\mathbb{\mathcal{F}}_{n}\bm{\Sigma}\mathbb{\mathcal{F}}_{n}^{H}-\bm{D}$ on a vector $\bm{v}$ with

(4.1)

(\mathbb{\mathcal{F}}_{n}\bm{\Sigma}\mathbb{\mathcal{F}}_{n}^{H}-\bm{D})\bm{v}%=\mathbb{\mathcal{F}}_{n}\left\{\mathbb{\mathcal{F}}_{2n-1}^{H}\left(\left\{%\mathbb{\mathcal{F}}_{2n-1}\bm{c}\right\}\circ\mathbb{\mathcal{F}}_{2n-1}%\begin{bmatrix}\mathbb{\mathcal{F}}_{n}^{H}\bm{v}\\\bm{0}_{n-1}\end{bmatrix}\right)\right\}_{1:n}-\bm{D}\bm{v}.

Assuming that $\mathbb{\mathcal{F}}_{2n-1}\bm{c}$ is pre-computed, this means that the action of $\mathbb{\mathcal{F}}\bm{\Sigma}\mathbb{\mathcal{F}}^{H}-\bm{D}$ on $\bm{v}$ requires four FFTs, two of size $n$ and two ofsize $2n-1$ .

With this established, we turn to the process of approximating $\bm{C}$ from theposition of only being able to apply it to vectors efficiently. The field ofrandomized low-rank approximations has become an important part of modernnumerical linear algebra, and we refer readers to review papers like[21] and the more recent [37] for broad introductionsand context. As it pertains to this work, it is sufficient to discuss theproblem of using randomized algorithms to estimate the range of a matrix.Letting $\bm{A}\in\mathbb{R}^{n\times n}$ denote an arbitrary matrix of rank $r$ ,[21] shows that the orthogonal matrix $\bm{Q}$ such that

\bm{A}\bm{\Omega}=\bm{Q}\bm{R},

where $\bm{\Omega}\in\mathbb{R}^{n\times(r+p)}$ is a “sketching”/test matrix, whichin this work will be made of i.i.d. standard normal random variables (althoughother and potentially faster choices are available), is an approximate basis forthe column space of $\bm{A}$ . The matrix $\bm{Q}$ itself may be obtained byperforming a standard QR factorization on the product $\bm{A}\bm{\Omega}$ , andthe parameter $p$ is an oversampling parameter which is often picked to be somesmall number like $p=5$ [21].From there, one obtains a simple low-rank representation of $\bm{A}$ as

\bm{A}\approx\underbrace{\bm{Q}}_{\bm{U}}\underbrace{\bm{Q}^{H}\bm{A}}_{\bm{V}%^{H}}.

This entire process is summarized in Algorithm 1.

1. Using a method from Section 2, obtain theautocovariance terms $\left\{h_{j}\right\}_{j=0}^{n-1}$ .

2. Create diagonal matrix $\bm{D}_{j,j}=S(\omega_{j})$ at Fourierfrequencies $\left\{\omega_{j}\right\}_{j=1}^{n}$ .

3. Draw sketching matrix $\bm{\Omega}\in\mathbb{R}^{n\times(p+r)}$ with $\bm{\Omega}_{j,k}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,1)$ .

4. Using Equation 4.1, compute $(\mathbb{\mathcal{F}}_{n}\bm{\Sigma}\mathbb{\mathcal{F}}_{n}^{H}-\bm{D})\bm{\Omega}$