Fast Machine-Precision Spectral Likelihoods for Stationary Time Series
Christopher J. Geoga


[
    [
    31/12/2021
==============




  We provide in this work an algorithm for approximating a very broad class of
symmetric Toeplitz matrices to machine precision in 𝒪(n log n)
time. In particular, for a Toeplitz matrix  with values _j,k =
h_j-k = ∫_-1/2^1/2 e^2 π i j-k S()
where S() is piecewise smooth, we give an approximation ^H
≈D + UV^H, where  is the DFT matrix, D is
diagonal, and the matrices U and V are in 𝒞^n ×
r with r ≪ n.  Studying these matrices in the context of time series, we
offer a theoretical explanation of this structure and connect it to existing
spectral-domain approximation frameworks. We then give a complete discussion
of the numerical method for assembling the approximation and demonstrate its
efficiency for improving Whittle-type likelihood approximations, including
dramatic examples where a correction of rank r = 2 to the standard Whittle
approximation increases the accuracy from 3 to 14 digits for a matrix ∈^10^5 × 10^5. The method and analysis of this work applies well
beyond time series analysis, providing an algorithm for extremely accurate
direct solves with a wide variety of symmetric Toeplitz matrices.  The analysis
employed here largely depends on asymptotic expansions of oscillatory integrals,
and also provides a new perspective on when existing spectral-domain
approximation methods for Gaussian log-likelihoods can be particularly
problematic. 





§ INTRODUCTION
 

Let ∈^n × n be a symmetric Toeplitz matrix with entries
_j,k = h_j-k = ∫_-1/2^1/2 e^2 π i j-k
S(). In this work we discuss approximations ^H
≈D + UV^H where D is diagonal, U, V∈ℂ^n × r form a matrix of rank r ≪ n, and  is the
discrete Fourier transform (DFT) matrix, parameterized in the unitary form with
the “fftshift" for convenience as = [n^-1/2 e^-2 π i j k/n ]_j ∈
(-n/2):(n/2-1), k ∈ 0:(n-1).  These matrices are of fundamenal importance in
time series analysis, and several motivating aspects of the derivation will be
done in that language.  But while time series analysis provides a convenient and
elegant way to study the problem, all of the tools discussed and applied below
are applicable to any symmetric Topelitz matrix whose “generating" function
S() (called a spectral density in the time series literature when
S() ≥ 0 and S() = S(-)) is known and meets some basic
requirements like being piecewise smooth. As will be demonstrated, however, the
method still applies when the function S() takes negative values, and so
this structure can potentially be exploited for a wide variety of symmetric
Toeplitz matrices and applications beyond just time series analysis.

To establish statistical notation and introduce relevant concepts for the
remainder of this work, consider the length-n mean-zero time series
Y_t_t=0^n-1.  If Y_t is stationary, which we will assume it
is for the duration of this work, then it has an autocovariance sequence given
by h_k = Cov(Y_t+k, Y_t) for any t. The Cramér spectral
representation theorem <cit.> states that there exists an
orthogonal increment process Z(ν) defined on (-1/2, 1/2] such that

    Y_t = ∫_(-1/2, 1/2] e^2 π i t ν dZ(ν).

Corresponding to Z there is a spectral distribution function F such that
F(τ) = 0 for τ≤ -1/2, F(τ) = h_0 for τ≥ 1/2, and
F(μ) - F(τ) = Var(Z(μ) - Z(τ)) for -1/2 < τ <= μ < 1/2
<cit.>. Herglotz's theorem gives that 

    h_k = ∫_(-1/2, 1/2] e^2 π i k ν dF(ν),

and if F has a density with respect to the Lebesgue measure, which we will
denote S, that density is called the spectral density. A basic fact
about covariance functions is that they must be positive (semi-)definite, and
that spectral densities, when they exist, must be non-negative (and symmetric
about the origin for real-valued processes). In light of this fact, it is
clearly easier to specify a parametric family of valid spectral densities than
it is a family of valid covariance functions. Unsurprisingly, then, in many
settings modeling dependence structure in the spectral domain is an appealing
choice for practitioners.

A second reason that makes modeling covariance structure in the spectral domain
appealing is that the DFT of the vector  (with _j = Y_j) often has a
better-behaved covariance matrix than  itself does. If Y_t is
measured with no gaps and a constant sampling rate (which will be assumed to be
unit sampling for the duration of this work without loss of generality), the
covariance matrix of , denoted  with _j,k = h_|j-k|, is
Toeplitz, meaning that it is constant along super- and sub-diagonal bands.
Toeplitz matrices made with absolutely summable sequences h_n_n
∈ are asymptotically connected to circulant matrices, a special sub-class of
Toeplitz matrices that are diagonalized by the DFT <cit.> and
correspond to periodic processes. This observation motivates the Whittle
approximation to the negative log-likelihood, which (omitting constants) is
given by

    2 ℓ^W() = ∑_j=1^nlog S_(_j) +
    J(_j)^2/S_(_j),

where S_() is a parametric spectral density, _j = 2 π
(j-n/2-1)n is a Fourier frequency with 1 ≤ j < n (note the offset here
to accomodate the “fftshift" in  but make the indices of _j and
(^H)_j,k in subsequent sections easier), and J() =
1√(n)∑_j=0^n-1 e^-2 π i j  Y_j is the DFT
coefficient for frequency  (which will only be computed for Fourier
frequencies via the FFT algorithm). Here and throughout, ℓ will be used to
denote the negative log-likelihood. Comparing this approximation to the
exact negative log-likelihood, given by (again omitting constants)

    2 ℓ() = log() + ^T ()^-1,

we can see that (<ref>) is an approximation that is motivated by
acting as if  were circulant—and specifically, the “right" circulant
matrix, where the DFT of the first column is precisely the vector of
S_(_j)_j=1^n. After pre-computing the DFT of ,
ℓ^W can be evaluated in linear complexity. This compares very well with the
cubic complexity of evaluating the true ℓ, which involves a very costly
Cholesky factorization. This approximation method is exceptionally popular for
its speed and ease of implementation.

Unfortunately,  is almost never exactly circulant for most popular classes
of models. And even if the circulant approximation is reasonably good, the
implied diagonals of ^H are unfortunately not the spectral density
at Fourier frequencies. While many textbooks and existing works in the
literature discuss the rates of convergence of Var(J(_j)) to
S(_j), it is useful and easy to compute the exact covariance of
finite-sample DFT coefficients, which also gives a convenient form for the
entries of ^H that doesn't involve a double sum. This is not a new
observation, but because it is referenced later in this work we present the
result as a proposition. And while again no aspect of this proof is new, the
argument will be used several times later in the work so we provide a terse
version of the proof here.

  Let Y_t_t=0^n-1 be a stationary mean-zero time series with spectral
  density S() and J(_k) defined as above, and define
  S_n(_k, _k') = Cov(J(_k), J(_k')) 
  Then
  
    S_n(_k, _k')
        =
        e^i n-1/nπ (k - k')/n∫_-1/2^1/2 
        D^s_n(k/n - )
        D^s_n(k'/n - )
        S() ,

  where D^s_n() = sin (π n )/sin (π) = e^-i (n-1)
  /2∑_j=0^n-1 e^2 π i j is a “shifted" Dirichlet kernel.


Using the spectral representation theorem stated above, the Fubini-Tonelli
theorem, and properties of orthogonal increment processes (see
<cit.> for an introduction), we see that

    Cov(J(_k), J(_k')) 
       = 
    ∫_-1/2^1/2( ∑_j=0^n-1 n^-1/2 e^2 π  i j ( - _k))
    ( ∑_j=0^n-1 n^-1/2 e^2 π i j (_k' - ))
     Z() Z() 
    
       =
    1/n∫_-1/2^1/2( 
    e^-i (n-1) π_ksin(π n (_k - ))/sin(π (_k - )))
    ( 
    e^i (n-1) π_k'sin(π n (_k' - ))/sin(π (_k' - )))
    S() ,

Applying the definition of D^s_n(), collecting terms, and bringing the
complex exponential prefactor on the sin terms outside the integral gives the
results.


There are many excellent discussions and investigations into the rate at which
this expression converges to S(_k) (see <cit.> for a
standard introduction), so we will not provide much discussion here. The
unfortunate result is that (<ref>) converges to S(_k)
very slowly when _k = _k' and zero otherwise, and while many
works give the rate of (n^-1) under ideal conditions, a more practical
issue is that the prefactor for this rate depends greatly on properties of S
and can be arbitrarily large <cit.>.  One example of a property that
affects |S(_k) - S_n(_k, _k)| is the dynamic range of S: when
_k = _k', (<ref>) is a convolution of the spectral density
with the Fejer kernel (denoted F_n() = |D_n^s()^2| here), which
has a dynamic range of about five orders of magnitude. If the ratio of the
largest and smallest values of S is much greater than 10^5, then the
convolution with F_n will dramatically affect the function at frequencies for
which the power is comparatively low. This observation was the fundamental
motivation for tapering; see <cit.> for more information.

Turning back to the Whittle estimator, it is well-known that, due to supplying
the incorrect variances for J(_j)_k=1^n, the minimizer of
ℓ^W can be quite biased, especially for small sample sizes
<cit.>. A particularly elegant proposition for dealing with this
issue is given in <cit.>: instead of using S_(_j) in
ℓ^W, directly use S_n(_j, _j) in its place.  By providing the
correct finite-sample variances for the DFT coefficients, one immediately
obtains a new likelihood approximation whose gradient gives unbiased estimating
equations for parameters θ <cit.>. In
<cit.>, they compute S_n(_j, _j)_j=1^n in the
time domain, using the classical identity that

    S_n(_j, _j) = 2 Re∑_k=0^n-1 (1 - n^-1 h) h_k e^- 2 π_j h - h_0,

and comment on how this identity is useful to avoid numerical integration.
While this is true, it is also quite inconvenient in that it
requires knowing the covariance function as well as the spectral density.
Considering that a large part of the point of resorting to a Whittle-type
approximation is to write a parametric spectral density instead of a covariance,
if one is unwilling to numerically integrate spectral densities then this
method is only available for the even more narrow class of models for which the
spectral density and covariance function can be written in closed form
and conveniently evaluated.

The approximation we introduce here is a continuation of such spectrally
motivated likelihood methods. Letting  denote the Toeplitz covariance
matrix for a time series with spectral density S() and D denote the
diagonal matrix with values S(_j)_j=1^n, the fundamental
observation of this work is that for a very broad class of spectral densities we
have that

    ^H ≈D + UV^H,

where U and V are in ℂ^n × r with r ≪ n. The
standard Whittle approximation views that low-rank update as being exactly zero,
and in this light our method can be viewed as a “corrected" Whittle
approximation. The key observation that makes working with this approximation
convenient is that the “Whittle correction" matrix = ^H - D,
along with being severely rank-deficient in many settings, can be applied to a
vector in (n log n) complexity via circulant embedding of Toeplitz
matrices, and so one can utilize the very powerful methods of randomized
low-rank approximation <cit.> to assemble U and V
efficiently and scalably.

Despite having the same quasilinear runtime complexity, for very small r
(sometimes as small as r=2), this representation can be exact to
computer precision. As one might expect, this does come at the cost of a more
expensive prefactor than a standard Whittle-type approximation.  If the standard
Whittle approximation requires one pre-computed FFT of the data and subsequently
runs at linear complexity, and the debiased Whittle approximation requires one
additional FFT per evaluation due to computation of the terms in
(<ref>), our approximation requires (3 r) FFTs to assemble
(although a simpler implementation using (4 r) FFTs is used in the
computations of this work[A software package and scripts for all
computations done in this work is available at
.]). For a spectral
density that is neither particularly well- nor poorly-behaved, a reasonable
expectation for an r that gives all 14 significant digits of ℓ() is around r ≈ 100—and so for full precision of the likelihood ,
this method can require several hundred FFTs. This is an undeniably more
expensive likelihood approximation, and there will of course be circumstances in
which a practitioner may choose a cheaper option at the cost of efficiency or
bias. But considering how fast FFTs are, the prefactor cost for this method is
sufficiently low that it compares favorably with other high-accuracy
alternatives.

More so then other spectral-domain approximation methods, this method is
fundamentally motivated by a matrix approximation. But unlike very technical
hierarchical matrix approximations that are often used to approximate covariance
matrices <cit.>,
this very simple representation means that it is trivial to compute exact
gradients of the log-likelihood, for example, which is a challenge that
hierarchical matrix-based methods often have to circumvent by using stochastic
trace approximations <cit.>. Further, it is
also trivial to implement symmetric factorizations, perform matrix-matrix
products, fully compute and instantiate inverses, and many other operations that
are often unavailable or problematically expensive to do with more sophisticated
matrix compression structures.

In the next several sections, we will introduce tools for implicitly assembling
and working with matrices of the form (<ref>). We start with a
discussion of using quadrature and asymptotic expansions to efficiently and
accurately work with  using only the spectral density S(), and we
then provide a discussion of (i) why the perturbation term in (<ref>)
should be low-rank, which also offers a new theoretical perspective for when
Whittle approximations are particularly problematic, and (ii) the details of
assembling the low-rank UV^H entirely from fast matrix-vector
products. We then close by providing several tests and benchmarks to demonstrate
the speed and accuracy of the method.





§ OSCILLATORY INTEGRATION OF SPECTRAL DENSITIES
 

As discussed in the previous section, creating the low-rank approximation ^H ≈D + UV^H requires working with , the
Toeplitz matrix with values _j,k = h_|j-k| where h_k is the
autocovariance sequence of the time series. Unlike in <cit.> where
both kernel and spectral density values are used to avoid numerical integration,
we will now discuss strategies for efficiently and accurately computing
integrals of the form

    h_k = ∫_-1/2^1/2 S() e^2 π i k

for k ∈ 0, ..., n so that one can create this factorization using
only the spectral density. Let us briefly review the challenges of such a task.
First, the standard method for obtaining the values h_k using the
trapezoidal rule via the FFT may require a large number of nodes for even
moderate accuracy if, for example, S() has rough points (like how S()
= e^-|| has no derivatives at the origin). A second common issue arises
if S() is not periodic at its endpoints, which limits the convergence rate
of the trapezoidal rule to its baseline rate of (m^-2) for m nodes
for general non-smooth and/or non-periodic integrands.  Even if one instead opts
to use a higher-order quadrature rule like Gauss-Legendre, for example, the
Shannon-Nyquist sampling theorem still requires at least (k) many nodes to
resolve the oscillations of e^2 π i k. The difficulty that this poses
is that if you need to evaluate h_k out to, say, k = 10 000, those last
values in the sequence will pose a significant computational burden.

To bypass these issues, we propose the use of adaptive Gauss-Legendre quadrature
for a fixed number of lags k, but then a transition to the use of asymptotic
expansions to h_k for lags greater than k. For integrands that are not too
oscillatory, say for k ≤ 2 000 (the cutoff used in all computations of
this work), standard adapative integration tools can be very efficient and
capable for any piecewise-smooth integrand. In this work, the Julia lanuage
package  <cit.> 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 numerical integration of
smooth functions. It gives approximations to integrals as

    ∫_a^b f(x)  x ≈∑_j=1^M w_j f(x_j)

that are exact for polynomials up to order 2M - 1, where w_j_j=1^M
are weights and x_j_j=1^M are nodes. In the case of [a,
b] = [-1, 1] the Legendre polynomials can be used to obtain the weights and
nodes and achieve an error rate of (M^-2m - 1)
<cit.>. These integrals can easily be made adaptive
by computing their discretized approximation for orders 2M-1 and 3M+1, as is
done in Gauss-Kronrod quadrature, for example, and concluding that the
approximation has converged when these two quantities agree to a prescribed
tolerance <cit.>.

More interesting, however, is the discussion of methods for accurately and
efficiently computing h_k when k is large. Unlike traditional integration
tools, asymptotic expansion-type methods get more accurate as k
increases. As in the introduction, the following result is not new (see, for
example, <cit.> for a comprehensive discussion), but we state it in
the specific form that is useful for this work since it will be referred to
frequently. Since the idea of the proof is also used repeatedly, we again
provide a terse proof.
 
  Let S() ∈𝒞^m([-1/2, 1/2]) and ∫_-1/2^1/2
S^(m)()  < C < ∞. Then

    ∫_-1/2^1/2 S() e^2 π i k 
      =
      - ∑_j=0^m-1 (-i 2 π k)^-(j+1) ( S^(j)(1/2) e^π i k - S^(j)(-1/2) e^- π i k)
    +
    (k^-m - 1).



  Simply apply integration by parts as many times as possible:

    ∫_-1/2^1/2 S() e^2 π i k    = 
      S(1/2) e^π i k - S(-1/2) e^-π i k/-2 π i k
      +
      (-2 π i k)^-1∫_-1/2^1/2 S'() e^2 π i k 
       ⋮
       =
      ∑_j=0^m-1 (-2 π i k)^-(j+1) ( S^(j)(1/2) e^π i k - S^(j)(-1/2) e^- π i k)
      
       +
      (- 2 π i k)^-(m+1)∫_-1/2^1/2 S^(m)() e^2 π i  k.


Remarkably, evaluating this expansion for any lag k only requires a few
derivatives of S() at its endpoints. For this reason, with this tool
evaluating the tail of h_k for high frequencies k actually becomes the
fastest part of the domain to handle. And while the supposition of the above
theorem that S ∈𝒞^m([-1/2, 1/2]) is obviously restrictive, a
simple observation extends this result to a much broader class of functions.


Let S() ∈𝒞^m([-1/2, 1/2] _1, ..., _L),
so that it is smooth except at locations _1, ..., _L, and
further assume that at each _j it has at least m directional
derivatives S^(m±)(_l) from both the left (m-) and the right
(m+) for l ∈ 1, ..., L.  Then

    ∫_-1/2^1/2 S() e^2 π i k    =
      ∑_l=1^L-1∑_j=0^m-1 (-2 π i k)^-(j+1) 
      (
      S^(j-)(_l+1) e^2 π i k _l+1 
      - 
      S^(j+)(_l) e^2 π i k _l 
      )
    +
    (k^-m - 1).


This corollary, which is proven by simply breaking up the domain [-1/2, 1/2]
into segments with endpoints at each _l and applying the above argument to
each segment, now means that this asymptotic expansion approach can be used to
accurately calculate h_k for very large k for any spectral density that is
just piecewise smooth. The prefactor on the error term in this setting is
natually increased compared to the setting of Proposition <ref>.
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 case
is given by C · (6,000 π)^-6≈ C · 10^-26. For a spectral
density like S() = 10 e^-10||, an example that will be studied
extensively in later sections due to being well-behaved except at the origin, we
see that h_3000≈ 5 · 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. Considering that this evaluation method uses
only derivatives of S() at the endpoint and takes just a few nanoseconds
on a modern CPU once those have been pre-computed, this combination of adaptive
Gauss-Legendre quadrature for sufficiently low lags k and asymptotic
expansions for the rest gives h_k_k=0^n-1 for all lags quickly and
accurately.





§ LOW RANK STRUCTURE OF THE WHITTLE CORRECTION
 

We now turn to the question of why one should expect the matrix

    C
      =
      ^H - D,

where D is diagonal with D_j,j = S(_j), to be of low
numerical rank. The primary tool for this investigation will again be asymptotic
expansions, and perhaps the crucial observation is that for every k and k',
C_k,k' is given by an oscillatory integral with the same high frequency
of 2 π n.  While the below analysis is not actually how we choose to compute
and assemble this low-rank approximation in this work, the following derivation
provides some idea for why the Whittle correction matrix C has low-rank
structure, even when S is not smooth and/or h_k decays slowly. For the
duration of this section, we will assume that S() ∈𝒞^(m)([-1/2, 1/2]) for convenience, although the results here can
again be extended using directional derivatives if S is not differentiable but
is smooth from the left and the right at rough points.

We begin by adding and subtracting S(k/n) and S(k'/n) to the inner integrand
in (<ref>), which with minor simplification steps can be expressed as

    1/n   ∫_-1/2^1/2 
      D_n^s(k/n - ) D_n^s(k'/n - )
      S()
      
    
      =
      1/2n   ∫_-1/2^1/2 
      D_n^s(k/n - ) D_n^s(k'/n - )
      (S() - S(k/n))
      
    
      +
      1/2n   ∫_-1/2^1/2 
      D_n^s(k/n - ) D_n^s(k'/n - )
      (S() - S(k'/n))
      
    
      +
      S(k/n) + S(k'/n)/2n   ∫_-1/2^1/2
      D_n^s(k/n - ) D_n^s(k'/n - )
      .

While unwieldy, this representation of the integral already provides a
reasonably direct explanation about several features of ^H.
Recalling that

    ∫_-1/2^1/2 D_n^s(k/n - ) D_n^s(k'/n - ) 
      =
      
        n     k = k' 
        
    
        0     k ≠ k' 
      ,

we see that the third term is precisely the diagonal contribution of
S(_k). Now we will now argue that the first two terms in the above sum
correspond to severely rank-deficient matrices.  Since they can of course be
analyzed in the exact same way, we study only the first in detail here. 

To begin, we slightly rewrite the first term in (<ref>) as

    ∫_-1/2^1/2sin(π n (k/n - )) sin(π n (k'/n - ))
      (π(k/n - )) (π(k'/n - ))
      (S() - S(k/n)) .

From here, we partition the domain into 

    [-1/2, 1/2] 
      =
      [-1/2, k/n - γ]_Type I∪B_γ(k/n)_Type II∪[k/n + γ, k'/n - γ]_Type I∪B_γ(k'/n)_Type II∪[k'/n + γ, 1/2]_Type I,

where B_γ(x) = [x - γ, x + γ] and γ is some small
number chosen to keep distance from the singularities of the cosecant
terms—for example, a value of γ = 0.001. By design, then, in Type I
intervals the cosecant terms are simple analytic functions. Defining 

    S̃_n, k, k'() = 
      (π(k/n - )) (π(k'/n - )) (S() - S(k/n)),

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 will
be used many times in this section.


If g() ∈𝒞^m([a,b]) and g^(m) is integrable on [a,b],
then
  
    ∫_a^b 
        sin(π n (k/n - )) sin(π n (k'/n - ))
        g() 
           =
        1/2[
          ∫_a^b g() ] 
        (-1)^k - k'
       -
        1/2
        e^i π n k+k'/2[
          e^2 π i b∑_j=0^m-1g^(j)(b)/(-i 2 π n)^j+1
          -
          e^2 π i a∑_j=0^m-1g^(j)(a)/(-i 2 π n)^j+1]
        
       +
        (n^-m-1).



  Using the product-to-sum angle formula that sin(θ) sin(ϕ) =
  12(cos(θ - ϕ) - cos(θ + ϕ)) and standard
  manipulations of the complex exponential, the left hand side equation can be
  re-written as
  
    1/2cos(π n k - k'/n) ∫_a^b g() 
      +
      1/2
        e^i π n k + k'/2∫_a^b e^-i 2 π n  g() .

  But cos(π n (k - k')) = (-1)^k-k', which provides the simplification of
  the first term. For the second term, just as in the proof of Proposition
  <ref>, we simply do integration by parts as many times as possible
  on the now standard-form oscillatory integral in the second term to get
  
    ∫_a^b e^-2 π i n  g() 
        =
        - ∑_j=0^m-1 (-2 π i n)^-j-1
        g^(j)(b) e^-2 π i n b
        -
        g^(j)(a) e^-2 π i n a
        +
        (n^-m-1).

  An elementary rearrangement gives the result.

While this proposition is presented in generality, there is a reasonable amount
of additional simplification that one can do with the right-hand side based on
the value of k + k' 4. In particular, we see that e^i π n
(k+k')/2 = i^k-k' 4, which means that every single complex
exponential in the right-hand side of (<ref>) simplifies nicely. If
we assume that k = k' 4 so that e^i π n (k+k')/2 = 1, the
asymptotic expansion term can be more concretely expanded to

    (cos(2 π b) + i sin(2 π b))
      ∑_j=0^m-1g^(j)(b)/(-i 2 π n)^j+1
      -
      (cos(2 π a) + i sin(2 π a))
      ∑_j=0^m-1g^(j)(a)/(-i 2 π n)^j+1
    
    =
       cos(2 π b) 
    ∑_j=0,   j  odd^m-1 (-1)^1+⌊ j/2 ⌋g^(j)(b)/(2 π n)^j+1
    -
    cos(2 π a) 
    ∑_j=0,   j  odd^m-1 (-1)^1+⌊ j/2 ⌋g^(j)(a)/(2 π n)^j+1
    
    +
       sin(2 π b) 
    ∑_j=0,   j  even^m-1 (-1)^2+⌊ j/2 ⌋g^(j)(b)/(2 π n)^j+1
    -
    sin(2 π a) 
    ∑_j=0,   j  even^m-1 (-1)^2+⌊ j/2 ⌋g^(j)(a)/(2 π n)^j+1.

And while we don't enumerate the other cases for the three values of e^i π n
(k+k')/2, the only difference is in the odd/even indexing and the alternating
sign inside the sum. The key takeaway here is that the entire integral
(<ref>) can actually be written as a very smooth integral term
(∫_a^b g()), four trigonometric functions that are
non-oscillatory since they have unit wavelength and a, b ∈ [-1/2, 1/2], and
a remainder term depending on high-order derivatives of g() on [a, b].
As it pertains to the Whittle correction matrix , for each of the Type-I
regions we have endpoints (a = -1/2, b = _k - γ), (a = _k +
γ, b = _k' - γ), and (a = _k'+γ, b = 1/2)
respectively, and so the effective contribution of the Type-I integrals to
_k,k' is given by 
∫_[-1/2, 1/2]  (B_γ(_k) ∪ B_γ(_k'))S̃_n,k,k'()
and a total of 12 trigonometric functions at the three sets of given endpoints
a and b. Given that the trigonometric functions are analytic and the
integral function is at least m+1-times differentiable everywhere by
assumption, by the observation that very smooth kernel matrices often exhibit
rapid spectral decay (an observation dating back to at least
<cit.> with the Fast Multipole Method (FMM)),  we expect the
contribution of the Type-I intervals to  to have exceptionally fast
spectral decay and be severely rank-deficient so long as the remainder from the
asymptotic expansion is small. With slightly more effort, we may
repeat this analysis on the Type-II regions with singularities.

For the first Type II region of [k/n - γ, k/n + γ], the story is
only slightly more complicated. This time we introduce the condensed notation of

    S̃_n, k'() = (π(k'/n - ))(S() - S(k/n))

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

    ∫_k/n + γ^k/n - γsin(π n (k/n - )) sin(π n (k'/n - ))
      (π(k/n - )) S̃_n, k'() .

The important observation to make here is that the singularity presented by
(π (k/n - )) is simple. Recalling the Laurent expansion (t) =
t^-1 + t6 + 7t^2360 + ..., we see that one can simply
“subtract off" the singularity and obtain a standard power-series type
representation of (t) - t^-1≈ P_l(t) for t ∈ [-γ,
γ] and some low-order polynomial P_l based on the Laurent series.  This
motivates the decomposition of (<ref>) into

    ∫_k/n + γ^k/n - γsin(π n (k/n - )) sin(π n (k'/n
      - )) P_l(k/n - ) S̃_n, k'() 
    
      +
         ∫_k/n + γ^k/n - γsin(π n (k/n - )) sin(π n (k'/n
      - )) S̃_n, k'()/k/n - .

The first term above can again be expanded as a small number of unit-wavelength
trigonometric functions via Proposition <ref> with g() =
P_l(k/n - ) S̃_n,k'(). Because S̃_n,k'() = 0 at
= k/n and has been assumed to be smooth, the second term is also not
singular and is a standard oscillatory integral. Another application of
Proposition <ref> with g() = (k/n - )^-1S̃_n,k'() can be applied, making the entire contribution of
(<ref>) expressible as a linear combination of a small number
of unit-wavelength trigonometric functions.

The final segment of [-1/2, 1/2] to study is the Type II region [k'/n - ,
k/n - ], where the singularity caused by (π (k'/n - )) has not
already been factored out. Introducing one final condensed integrand notation of

    S̃_n, k() = (π(k/n - ))(S() - S(k/n)),

we again decompose the contribution of that segment as

    ∫_k'/n + γ^k'/n - γsin(π n (k/n - )) sin(π n (k'/n - ))
      (π(k'/n - )) S̃_n, k() 
    
      =
         ∫_k'/n + γ^k'/n - γsin(π n (k/n - )) sin(π n
      (k'/n - )) P_l(k'/n - ) S̃_n, k() 
    
      +
         ∫_k'/n + γ^k'/n - γsin(π n (k/n - )) sin(π n
      (k'/n - )) S̃_n, k()/k'/n - .
 
The first term in the divided integral can be handled yet again by Proposition
<ref> just as above. The only new wrinkle is that the second term
in the right-hand side now still has a singularity. To handle this last term,
recall the elementary trick that, for simple singular integrals of smooth
functions, one can again subtract off the singularity in the sense of

    ∫_a^b x^-1 g(x)  x
      =
      ∫_a^b x^-1 (g(x) - g(0))  x + g(0) logb/a,

where 0 ∈ (a, b), g(x) is differentiable at zero, and the logarithmic term
is the result of a Cauchy principal-value interpretation of the integral. With
this trick in mind, we just split that last term one more time to obtain

    ∫_k'/n + γ^k'/n - γsin(π n (k/n - )) sin(π n
      (k'/n - )) S̃_n, k()/k'/n - 
    
      =
         ∫_k'/n + γ^k'/n - γsin(π n (k/n - )) sin(π n
      (k'/n - )) S̃_n, k() - S̃_n,k(k'/n)/k'/n - 
    
      +
      S̃_n,k(k'/n) 
         ∫_k'/n + γ^k'/n - γsin(π n (k/n - )) sin(π n (k'/n - )) 1/k'/n - .

Now that the singularity has been removed, the first term on the right-hand side
can once again be expressed as a small sum of unit-length trigonometric
functions via Proposition <ref>. With only a slightly more
technical application of the singularity subtraction trick, the final term can
be written in the Cauchy principle value sense as a sum of a small number of
logarithmic terms including k/n, k'/n, and γ, although we omit the
full expression due to its length.

Combining all of these segmented contributions, we finally conclude in
representing _k,k' as a small number of very smooth functions. And while
ostensibly breaking the integral into five pieces, two of which need to be
broken up again into two or three additional terms to be suitable for expansion,
should lead to so many sines and cosines that the matrix is not particularly
rank-deficient, we note that all of these trigonometric functions have the exact
same wavelength, and  many or all are nearly or exactly linearly dependent with
each other. As we show in the numerical demonstrations section, for sufficiently
well-behaved spectral densities S one can go from three correct digits with a
diagonal approximation to 14 correct digits with only a rank r=2
approximation to , even when ∈^10^5 × 10^5. And even for
less well-behaved S(), such as one with no derivatives at the origin, a
rank of r=128 is sufficient for 14 digits even for a  of the same size.

This derivation is also illuminating about the sources of error in Whittle
approximations and the degree of rank-structure in the Whittle correction
matrices. For example, if S has rough points where it is continuous but not
smooth, then the above domain partitioning will have to be refined so that
Corollary 1 can be applied to fix the accuracy of the asymptotic
approximation. This will then mean that _k,k' will need more asymptotic
expansion-sourced terms for the expansion to be accurate. Moreover, we see that
if the remainder term in the asymptotic expansion is large, there is potentially
a great deal of structure in that matrix that this examination does not account for.
One such setting where that can easily happen is when S() itself is highly
oscillatory. As an example, A frustrating challenge for many fast matrix methods is the
triangular covariance function with a large bandwidth, such as h_k = max(0,
1- |k|/r) with large r. Picking r = 3n/4 will provide a 
that is likely to be problematic for hierarchical matrix-type methods because it
is zero on the outer triangles of the largest off-diagonal n/2 × n/2
blocks, making low-rank approximations for them impossible (or at least
inaccurate). In this setting, that problem is more directly observed in the
spectral domain. The corresponding spectral density to that sequence h_k
is S() = sin(2 π r /2)^2/r sin(2 π/2)^2, which
produces a significant remainder term ∫_-1/2^1/2 S^(m)() e^2 π
i  (k-k') that is not accounted for in the asymptotic expansion
terms. This makes any order of asymptotic expansion inaccurate, and leads to
Whittle correction matrices  that are nearly full rank.

As a final note on the structure being exploited here, it is interesting to
consider the Whittle correction matrix  in the context of other fast
algorithms for special matrix-vector products.  Particularly after subtracting
off singularities, one may view  as a kernel matrix made with a bandlimited
function, which is precisely the type of structure that is exploited in existing
methods like the fast sinc transform and fast Gauss transform
<cit.>. But unlike in those cases, this compactly
supported Fourier transform is highly oscillatory for every entry. This
is exactly what we exploit here, however, to understand the low-rank structure
of : every entry of C is the result of a high frequency
oscillatory integral of a compactly supported function, and so asymptotic
expansion-type methods can be employed equally accurately for every component. 




§ NUMERICAL ASSEMBLY AND ACTION OF THE WHITTLE CORRECTION
 
With this theoretical justification for the low-rank structure of the Whittle
correction matrix, we now turn to discussing how to actually compute and
assemble the matrices U and V such that ^H - D≈UV^H. First, we note that the action of  onto vectors
can be done in (n log n) with the familiar circulant embedding technique
<cit.> once the values h_k_k=0^n-1 have been computed.
While we refer readers to the references for a more complete discussion of
circulant embedding, we will provide a brief discussion here. The symmetric
Toeplitz matrix  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
c = [h_0, h_1, ..., h_n-1, h_n-2, ..., h_1], and if one assembles a
Toeplitz matrix with this length 2n - 1 column the matrix is in fact
circulant, which means it is diagonalized by the FFT <cit.>. Since the
FFT can be applied in (n log n) complexity, so can the action of this
augmented circulant matrix. Thus, one can compute  for any vector 
by extracting the first n components of C [, 0_n-1], where
C is the circulant matrix made with c and 0_n-1 is n-1
many zeros appended to .

Letting _n be the FFT of size n (briefly now defined without the
“fftshift"), we can thus summarize the action of _n _n^H - D
on a vector  with

    (_n _n^H - D) 
      =
      _n {_2n-1^H
        (
        _2n-1c∘_2n-1_n^H 
    0_n-1)
      }_1:n
      -
      D.

Assuming that _2n-1c is pre-computed, this means that the action of
^H - D on  requires four FFTs, two of size n and two of
size 2n-1. 

With this established, we turn to the process of approximating  from the
position of only being able to apply it to vectors efficiently.  The field of
randomized low-rank approximations has become an important part of modern
numerical linear algebra, and we refer readers to review papers like
<cit.> and the more recent <cit.> for broad introductions
and context. As it pertains to this work, it is sufficient to discuss the
problem of using randomized algorithms to estimate the range of a matrix.
Letting A∈^n × n denote an arbitrary matrix of rank r,
<cit.> shows that the orthogonal matrix Q such that

    AΩ = QR,

where Ω∈^n × (r+p) is a “sketching"/test matrix, which
in this work will be made of i.i.d. standard normal random variables (although
other and potentially faster choices are available), is an approximate basis for
the column space of A. Here p is an oversampling parameter which is
often picked to be some small number like p=5. From there, one obtains a
simple low-rank representation of A as

    A≈Q_UQ^H
      A_V^H.

This low-rank representation can easily be converted to other truncated
factorization types like a partial SVD or eigendecomposition (the representation
used in the software implementation of this paper) <cit.>. Connecting
this to Section <ref>, then, we see that obtaining a rank r
approximation ^H - D requires applying it to r + p many
random vectors. While there is sufficient structure in this specific
matrix-vector application that one could reduce the number of FFTs from
4(r+p), the implementation used in this work does not employ any particular
optimization of that form and is nonetheless satisfyingly fast for an
effectively exact method.

Armed with this low-rank approximation, one now has all of the pieces for
working with 

    ^H ≈D + UV^H.

Evaluating the likelihood for a vector  (Equation <ref>) in (n) time
is now trivial with this representation, first by using the matrix determinant
lemma to see that

    logD + UV^H
      =
      log_r + V^H D^-1U + logD

and then using the Sherman-Morrison-Woodbury formula to compute the quadratic
form as 

    ^T (D + UV^H)^-1
    =
      ^T D^-1
      -
      ^T D^-1U(_r + V^H D^-1U)^-1V^H
      D^-1.

Since all of the matrix-matrix operations in the right-hand sides of those two
equations are for small ^r × r matrices, we see that both the
log-determinant and the quadratic form can be evaluted in (n) complexity.
And since the assembly of D + UV^H is done in (n log n),
we conclude that the full job of assembling the approximation and then
evaluating the approximated log-likelihood runs in quasilinear time complexity and
linear storage complexity.

As will be highlighted in the next section, one particularly attractive aspect
of this specific approximation structure is that it is very simple. Due to this
simplicity, two approximated DFT-conjugated Toeplitz matrices can be exactly
multiplied in (n) complexity, for example, and even the inverse of one can
be applied to the other at the same cost. Let us now discuss how to
continue exploiting this structure of ^H in obtaining similarly
accurate gradients and information matrices for a wide range of models.





§ GRADIENTS, FACTORIZATION, AND INFORMATION MATRICES
 

If one assumes that a parametrically indexed spectral density S_()
has partials derivatives ∂_θ_j S_() that share the
same smoothness structure as S_ with respect to  (such as being
piecewise 𝒞^m([-1/2, 1/2])), then the above arguments about the low-rank
structure of ^H apply equally well to ∂_θ_j^H. For notational clarity, introduce the
subscript j and let D_j, U_j, and V_j be matrices such
that

    ∂/∂θ_j() ^H
      ≈D_j + U_j V_j^H.

The j-th term of the gradient of the Gaussian negative log-likeklihood
ℓ(; ) is given by

    2 [∇ℓ()]_j
      =
      tr(
      ()^-1∂/∂θ_j()
      )
      -
      ^H
      ()^-1∂/∂θ_j()
      ()^-1.

Observing this equivalent structure in (<ref>), we note that the
quadratic form term in (<ref>) can easily be evaluated in (n)
complexity once the two low-rank approximations have been assembled by again
making use of the Sherman-Morrison-Woodbury formula and the speed of simple
matrix-vector products. More interesting, however, and what sets this
approximation method apart from more complex general-purpose matrix compression
methods is that the trace can also be computed exactly without resorting to
stochastic trace estimation (see <cit.> and references therein for
discussion). In particular, note that the Sherman-Morrison-Woodbury formula
given in the above section implies that (D + UV^H)^-1 can be
again represented as a low rank perturbation of a diagonal matrix, which we will
denote as D̃ + ŨC̃Ṽ^H.
Applying this, we see that

    (D + UV^H)^-1(D_j + U_j V_j^H)
         =
      (D̃ + ŨC̃Ṽ^H)
      (D_j + U_j V_j^H)
      
       =
      D̃D_j
      +
      D̃U_j V^H
      +
      ŨC̃Ṽ_H D_j
      +
      ŨC̃Ṽ^H
      D_j + U_j V_j^H
      
       =
      D̃D_j
      +
      M.

While this looks problematic, we note that the matrix M defined as the
last three terms in the above equation is a sum of three matrices whose rank is
≤ r, and so the rank of M is at most 3 r. Second, because each of
those terms is already a low-rank representation that can be applied to vectors
in (n) time, we can simply repeat the randomized low-rank approximation
strategy from earlier, computing first a randomized basis for the column space
with 

    (D̃U_j V^H
      +
      ŨC̃Ṽ_H D_j
      +
      ŨC̃Ṽ^H
      D_j + U_j V_j^H)
      Ω
      =
      QR

and then re-compressing to a single low-rank representation with QQ^H M. From there, we observe the simple fact that tr(A
+ BC^T) = tr(A) + tr(C^T B) and
conclude that the trace in (<ref>) can be computed exactly and in (n)
time. As will be demonstrated in the next section, these gradients are computed
to similarly high accuracy as the log-likelihood evaluation itself, providing
effectively the same number of digits.

A second pleasant property of this approximation structure for  is that it
admits a simple and direct symmetric factorization. In the above discussions for
notational simplicity we have worked with the low rank form D + UV^H. But as previously mentioned, these low-rank representations can be
easily and quickly converted to other structures like a partial SVD or
eigendecomposition. For this factorization, it is easiest to work with a
representation like D + UΛU^H, where Λ∈^r × r is a diagonal matrix of eigenvalues of the low-rank Whittle
correction. In this form, a slight manipulation and then applying Theorem 3.1
of <cit.> gives a very simple symmetric factorization as

    D + UΛU^H
      =
      D^1/2( + ŨΛŨ^H)D^1/2
      =
      D^1/2
      ( + ŨXŨ^H)
      ( + ŨXŨ^H)^H
      D^1/2
      =
      WW^H,

where Ũ = D^-1/2U and X = L^-H (M
- ) L^-1, U^T U = LL^T, and + L^T
ΛL = MM^T. While that is many equations, the actual
computation of the matrix X used in assembling W is is done in
just a few ^r × r matrix-matrix operations (Algorithms 1 and 2 in
<cit.>).

This symmetric factor is useful for several purposes. For one, it gives a means
of obtaining exact simulations of any time series in (n log n) time,
although of course more standard periodic embedding would be a more simple and
efficient way of achieving the same thing. In the setting of derivative
information, it offers a more novel benefit in that it enables very accurate
stochastic estimation of Fisher information matrices.

The Fisher information matrix for parameters  in a Gaussian model like the
one used in this work is given by, again in terms of just (),

    I()_j,k = 1/2tr(
        ()^-1∂/∂θ_j()
        ()^-1∂/∂θ_k()
        ).

This matrix is the asymptotic precision of the MLE under sufficient regularity
conditions, and is often substituted in place of the Hessian matrix of a
log-likelihood ℓ() due to the fact that it is easier to compute and is
always positive definite. As discussed above, it is fully possible to scalably
evaluate this sequence of matrix-matrix products to obtain an exact Fisher
information matrix, and in the following section we will provide a verification
that this computation still runs in quasilinear complexity despite the large
number of matrix-matrix operations. But <cit.> introduced a fast “symmetrized"
stochastic estimator the matrix I() that is very accurate and can be
computed in a fraction of the runtime cost of the exact I().
Stochastic trace estimation is a rich field with a broad literature (we refer
readers to <cit.> for a useful introduction and overview), but for this
work the fundamental observation to make is that for any random vector u
with = 0 and =, we have that ^T A =
tr(A) for any matrix A.  The sample average
approximation (SAA) trace estimator is then based on drawing several i.i.d.
vectors , denoted _j_j=1^M, and instead of evaluating the
trace directly using 

    M^-1∑_j=1^M _j^T A≈tr(A).

The variance of this estimator depends on several properties of A, but a
specific choice of  having i.i.d. random signs is particular popular
because ^T A = 2([F]A^2 - ∑_j A_j,j^2), so that
for diagonally concentrated matrices this estimator can be quite accurate. 

Stochastic trace estimation has been applied to the Gaussian process computation
problem in many works, notably first in <cit.>. In <cit.>, a
theorem was provided indicating that if = WW^T, then using the
simple property that tr(AB) = tr(BA)
combined with this factorization one can instead compute I_j,k with 

    I_j,k = 1/2tr(
        W^-T∂/∂θ_j()
        ()^-1∂/∂θ_k()
        W^-1),

and the variance of the SAA estimator for (<ref>) is bounded above by
the variance of the SAA estimator for (<ref>). This theoretical
observation was verified in <cit.> and subsequent works, and
<cit.> provided an even further “symmetrized" trace to estimate for
off-diagonal elements. Introducing the notation _j = ∂_θ_j(), one may compute a fully “symmetrized" trace with 

    I_j,k = 
      (4 M)^-1∑_l=1^M
      _l^T W^-1 (_j + _k) ^-1 (_j + _k) W^-T_l
      - 1/2I_j,j - 1/2I_k,k,

where diagonal elements I_j,j are trivial to fully symmetrize and thus
can be computed in advance of the off-diagonal ones.

This estimator I for I enjoys both very high accuracy
and great computational benefits. The primary computational benefit is that this
can be computed in a single pass over the derivative matrices _j, since all
one needs to evaluate (<ref>) is ^-1 and _j
W^-T_l_l=1^M. Since one commonly uses M ≈ 70 or some
similarly small number of SAA vectors, for example, it is very easy to pre-solve
the SAA vectors _l with W^-T and then in a single for-loop
assemble and apply _j to each of them. Whether those pre-applied vectors
are saved to disk or kept in RAM, there is no more need for derivative matrices
_j after that point, and so one need not ever even have _j and _k
instantiated at the same time to fully evaluate I. For models
with many parameters, the speedup that this affords can be substantial.

The Hessian of ℓ(), while again absolutely computable in (n
log n) complexity, is a similar situation to the expected Fisher information
matrix. It has terms given by (again using just  instead of ^H)

    [H ℓ()]_j,k
      =
      I_j,k
      + 
      tr(
        ^-1∂_θ_k∂_θ_j)
      -
      ^H ( ∂_θ_k^-1_j ^-1)
      .

From this expression one sees that the same things are possible as with the
gradient, but the prefactor on all operations will be larger: one must compute
and re-compress matrix products involving second derivatives of D + UV^H with respect to model paramaeters, and one must again do a full
matrix-matrix multiply for the additional trace (unless one opts to again use
SAA). All of these things are perfectly computable using the tools already
introduced in this work. But because of the high prefactor cost, we do not
explore them further here.
  




§ NUMERICAL DEMONSTRATIONS
 

In this section we will demonstrate the efficiency and accuracy of this method
for approximating . In particular, we will study two spectral densities:

    S_1()    = θ_1 (1 - 2 θ_2 cos(2 π) + θ_2^2)^-1
    
      S_2()    = θ_1 e^-θ_2||,

which have corresponding autocovariance sequences

    h_1,k   = θ_1 (1 - θ_2^2)^-1θ_2^k  
    
      h_2,k   = 2 θ_1 e^-θ_2 /2(-θ_2 cos(π k) + e^θ_2
      /2θ_2 + 2 π k sin(π k))/θ_2^2 + (2 π k)^2.

The model S_1 was chosen because it is not only very smooth and with a limited
dynamic range, but its periodic extension is even continuous at the endpoints
(which the above section demonstrates are often the dominating source of
structure in the Whittle correction for otherwise smooth SDFs).
This spectral density is in some sense maximally well-behaved, and its
corresponding autocovariance sequence decays predictably quickly. From the
perspective of this approximation framework, this is an ideal setting.
S_2(), on the other hand, is more challenging because of its roughness at
the origin. This model was chosen to demonstrate that the directional derivative
correction introduced above can fully restore the accuracy of asymptotic
expansions even for spectral densities that aren't even always once
differentiable, and that the rank structure of  is still very much exploitable.

As a first investigation, we verify the quasilinear runtime complexity of
assembling the approximation (<ref>). In particular, we approximate
ℓ() with 

    2 ℓ()
      ≈ 
      2 ℓ̃()
      =
      logD() + () ()^H)
      +
      ()^H D() + () ()^H
      ^-1
      (),

which follows from the basic fact that =  () ^H.

Figure <ref> shows the runtime cost of assembling D + ^H for a variety of data sizes ranging from n=1 000 to n=96 000 along
with a theoretical (n log n) line. As one would expect for a fixed-rank
approximation assembled using a (1)-sized number of FFTs, the empirical
agreement with the theoretical complexity is good. While there is interesting
variation in runtime cost when the entire assembly runs in significantly under one second,
by the time the program runtime cost approaches a second the stability of the
scaling is clear.

As a second investigation, we look at the relative error in the negative
log-likelihood computed with the implied approximation to (<ref>).  This
study, summarized in Figure <ref>, clearly demonstrates the way
that this approximation depends on the specific spectral density being modeled.
For S_1, we see that even a rank of r=2 provides an effectively exact
log-likelihood with 14+ digits of precision. The second two sub-plots in
Figure <ref> show two different applications of S_2. The center
figure shows the relative accuracy of the negative log-likelihood when one does
not correct the asymptotic expansion for the rough point at the origin
when computing the autocovariance sequence h_k_k=0^n-1, instead
simply pretending that the function has many derivatives everywhere.  In this
work the asymptotic expansions were used for k > 2 000, and it is clear in
the figure that as soon as those expansions start being used that one goes from
14 digits in the rank r=128 case to 8 digits, and the quality of the
approximation does not meaningfully improve even as the rank varies by an order
of magnitude.


The final panel in Figure <ref> shows the relative error of the
likelihood approximation (<ref>) with the corrected asymptotic expansion
that applies Corollary <ref> at the splitting point _1 = 0.
Here we now see that there is no visible loss in accuracy once the expansions
start being used for tails of the autocovariance sequence, indicating that the
domain-splitting correction restores the accuracy of the expansions. As the
theory would imply and  this figure empirically demonstrates, the rank of the
“Whittle correction" depends on properties of the spectral density.  In this
particular case, a rank of r ≈ 100 is necessary to get almost every
significant digit of the log-likelihood. This low rank makes sense in light of
Section <ref>'s asymptotic expansion argument, as S_2()
is non-smooth at the origin and at the periodic endpoints.


Finally, we validate the runtime cost and accuracy of the gradient, the exact
expected Fisher matrix, and the symmetrized stochastic expected Fisher matrix of
the negative log-likelihood, although due to the exact values requiring
matrix-matrix products we now only test to matrices of size 48 000.  Figure
<ref> provides a summary of the runtime cost of computing
the three quantities. As can be seen, the agreement with the theoretically
expected (n log n) is clear. And even for n ≈ 50 000 data points,
we see that the gradient of the likelihood—which requires a matrix-matrix
product for each term—can be computed to full precision in under ten seconds.
By virtue of requiring only matrix-vector products, which are exceptionally fast
with this simple low-rank structure, the symmetrized stochastic expected Fisher
information matrix is actually slightly faster to compute than the gradient.
Each term of the exact expected Fisher information matrix requires a product of
four matrices, and so the prefactor of the (n log n) complexity is
significantly higher, with the rank r=128 value for n=48 000 taking about
40 seconds compared to the ≈ 5 s that its stochastic
counterpart required. This is of course more expensive, but it is also
sufficiently manageable that one oculd certainly compute it once to
high accuracy to obtain the true asymptotic precision of the MLE once the
estimate _MLE has been computed, for example.

To study the accuracy of these approximations, Figure <ref> shows
two quantities. Letting g() denote the true gradient ∇ℓ() and g̃() denote the approximated gradient
∇ℓ̃(), the top row shows a standard max-norm type
metric of 

    max_j=1, ..., pg_j() -
      g̃_j()/g_j()

for models S_1() and S_2().  Letting I denote the true
expected Fisher information matrix and Ĩ denote the approximated
one (with just the approximated  or additionally by using the stochastic
trace estimator), the second row of Figure <ref> shows the
relative operator norm error
I^-1I - Ĩ,
which is arguably a more relevant metric for a Hessian approximation than a
pointwise max-norm. 

Many of the conclusions from examining the relative error in the log-likelihood
above still apply here. As before, for spectral density S_1 the matrix 
is so well-behaved that even a rank r=4 approximation to ^H -
D and its derivatives gives 12+ digits of accuracy in the gradient and
exact expected Fisher matrix. For the likelihood of S_2, due to the roughness
at the origin the rank of the Whittle correction again needs to be increased.
But as one would hope, since a rank of r=128 gave effectively every digit in
the log-likelihood, it also gives effectively every digit in the gradient and
exact expected Fisher matrix. Considering that the derivatives of S_1 and
S_2 with respect to θ_2 are not strictly non-negative, this also
provides a verification that the structure this method exploits works equally
well for symmetric Toeplitz matrices that are not positive-definite.

For the stochastic expected Fisher matrix, every trace approximation was done
using M = 72 SAA vectors and the SAA-induced error is clearly the dominant
source of disagreement with the true expected Fisher matrix, even for r=2.
While it is tempting conclude that one could simply set r=2 and obtain a good
estimator for the expected Fisher information matrix in under a second even for
n = 50 000 data points, which may well be true in plenty of circumstances,
the situation can be slightly more complicated than that. When using exact
matrices , _j, and _k, the symmetrized stochastic estimator is an
unbiased estimator. And while the SAA-induced variability in the approximated
expected Fisher using the matrix approximation here may be much larger
than the bias induced by setting the rank r too small, the estimator now
having some amount of bias may be a potential source of issues, so we advise
practitioners to reduce r from what they deemed necessary for an accurate
log-likelihood with caution.




§ DISCUSSION
 

This work introduces a pleasingly simple approximation for the covariance matrix
of the DFT of a stationary time series, which incidentally may provide useful
methods to a much broader range of settings where one encounters symmetric
Toeplitz matrices whose diagonal sequences have well-behaved Fourier transforms.
In many settings and more general Gaussian process applications, for example
with irregular locations or in multiple dimensions, low rank-type methods are
known to have severe limitations <cit.>. And for purely approximating
the Toeplitz matrix , many of those issues would still be present, so it
is interesting to observe that conjugating with the DFT matrix changes the
algebraic structure of  to the degree where a low-rank perturbation of a
diagonal matrix represents ^H to almost every achievable
significant digit at double precision.

The theory in Section <ref> provides a new and complementary
perspective to existing analyses of the accuracy of Whittle approximations.  In
<cit.>, for example, a theoretical analysis of spectral-domain
approximations is provided where error is characterized by an integer K and
the assumption that the corresponding covariance sequence has a finite K-th
moment ∑_j j^K h_j < ∞. This work, on the other hand, discusses
non-smooth points of spectral densities and large high-order derivatives as the
principle source of error.  These are of course intimately connected, and the
correspondence of a function's decay rate with the smoothness of its Fourier
transform is well-known. If one is designing parametric models via the spectral
density, the considerations that this work motivatives for designing models that
can be approximated well—namely, selecting and designing the number of
non-smooth points, including at endpoints, and deciding on the magnitude of
quantities like ∫_-1/2^1/2 S^(m)() e^2 π i  k—may be more direct in some circumstances.

Many methods already exist for exploiting structure in covariance matrices like
 and achieving similar or even slightly better runtime complexity.
<cit.> provides a hierarchically semi-separable-type (HSS)
approximation to covariance matrices  that runs in linear time, and
<cit.> use hierachically off-diagonal low-rank structure
(HODLR) to achieve quasilinear complexity that is comparable to that of this
work. Many other hierarchical matrix formats exist <cit.>, and in some
cases can be pushed to a similar degree of accuracy as the method presented
here. What makes this structure representation interesting, in our opinion, is
its simplicity. This approximation is easy to implement in software (and the
released code can be inspected to prove that point), and the resulting D +
UV^H has a very transparent structure. Moreover, its runtime
complexity is very simple and predictable since it is only controlled by one
tuning parameter (the rank r). Despite this, it achieves high accuracy for a
very broad class of spectral densities. And for users who are so inclined,
making its rank adaptive would be a very simple modification to the code,
following, for example, guidelines from <cit.>. Moreover, error control
and analysis for this structure is much more direct than in a more general
hierarchical matrix format, which must control the error of individual low rank
approximations of off-diagonal blocks.

With that said, there are of course circumstances in which the cost of using
this log-likelihood evaluation strategy is not worth the higher prefactor than,
say, the debiased Whittle approximation <cit.>. If one has a
dataset with 50 000 measurements that needs to be fitted once, for example,
it is probably worth spending an entire two to three minutes to obtain an
effectively exact MLE instead of the 15 seconds it might take to obtain an
approximate estimator from, say, the debiased Whittle method. Or at the least,
it is likely worth refining that more computationally expedient estimator using
the more accurate methods introduced here. But for an online application, on the
other hand, or a setting in which one must fit tens or hundreds of datasets, a
much cheaper method that runs in seconds instead of minutes may be worth the
potential efficiency cost.






  The author is grateful to Paul G. Beckman for his thoughtful comments, 
  discussion, and proofreading of the manuscript.




