\section{Introduction}

Gaussian Processes (GPs) \citep{rasmussen2006gaussian} are flexible Bayesian nonparametric models which enable principled reasoning about distributions of functions and provide rigorous uncertainty estimates \citep{srinivas2010gaussian, deisenroth2011pilco}. Unfortunately, exact inference in GPs is impractical for large datasets because of the $\mathcal{O}(N^3)$ computational cost (for a dataset of size $N$). To overcome the computational roadblocks, sparse Gaussian processes \citep{snelson2006sparse, quinonero2005unifying} use $M$ inducing points to approximate the kernel function, reducing the computational cost to $\mathcal{O}(NM^2 + M^3)$. However,
these approaches are prone to overfitting since all inducing points are hyperparameters.
Sparse variational Gaussian Processes (SVGPs) \citep{titsias2009variational, hensman2015scalable}
offer an effective protection against overfitting by
framing a posterior approximation using the inducing points and
optimizing them
with variational inference.
Still, the $\mathcal{O}(M^3)$ complexity prevents SVGPs from scaling beyond a few thousand inducing points, creating difficulties in improving the quality of approximation.

Several approaches impose structure on the inducing points to increase the approximation capacity. Structured kernel interpolation (SKI) \citep{wilson2015kernel, wilson2015thoughts} approximates the kernel by placing inducing points over a Euclidean grid and exploiting fast structured matrix operations. SKI can use millions of inducing points, but is limited to low-dimensional problems because the grid size grows exponentially with the input dimension.  Other approaches define approximate posteriors using multiple sets of inducing points. \citet{cheng2017variational, salimbeni2018orthogonally} propose to
decouple the inducing points for modelling means and covariances,
leading to a linear complexity with respect to the number of mean inducing points. SOLVE-GP~\citep{shi2020sparse} reformulates a GP as the sum of two orthogonal processes and uses distinct groups of inducing points for each; this has the benefit of
improving the approximation at a lower cost than standard SVGPs.

\begin{figure}[t]
\centering
\includegraphics[width=\columnwidth]{figures/rotate_real.pdf}
\caption{Visualizing the harmonic kernel decomposition.
We decompose a $2$-dimensional RBF kernel as $k = \sum_{t=0}^3 k_t$ using the symmetry group of $90^\text{o}$ rotations.
We plot the real part of random functions sampled from each $\mathcal{GP}(0, k_t)$. Notice that $\mathcal{GP}(0, k_0)$ is invariant to $90^\text{o}$ rotations; $\mathcal{GP}(0, k_1)$ takes opposite values under $180^\text{o}$ rotations; $\mathcal{GP}(0, k_2)$ is invariant to $180^\text{o}$ rotations but has opposite values under $90^\text{o}$ rotations.
\label{fig:rbf-rotate-gp}}
\end{figure}

In this paper, we introduce a more scalable variational approximation for GPs via the proposed \emph{harmonic kernel decomposition (HKD)}, which decomposes the kernel as a sum of orthogonal kernels, $k(\mathbf{x}, \mathbf{x}') = \sum_{t=0}^{T-1} k_t(\mathbf{x}, \mathbf{x}')$,
using Fourier series.
The HKD reformulates the Gaussian process into an additive GP,
where each
subprocess
models a Fourier component of the function (see Figure~\ref{fig:rbf-rotate-gp} for a visualization).
We then propose the Harmonic Variational Gaussian Process~(HVGP), which uses separate sets of inducing points for the subprocesses.
Compared to a standard variational approximation, HVGPs allow us to use a large number of inducing points at a much lower computational cost.
Moreover, HVGPs have an advantage over SKI in that they allow trainable inducing points.
Finally, unlike the SOLVE-GP whose decomposition involves only two subprocesses, our HVGP is easily applicable to multiple subprocesses and can be computed efficiently with the discrete Fourier transform.

Empirically, we demonstrate the scalability and general applicability of HVGPs through a range of problems and models including using RBF kernels for modelling earth elevations and using convolutional kernels for image classification. In these experiments, HVGPs significantly outperform standard variational methods by exploiting the input-space symmetries such as translations and reflections. Our model can further exploit parallelism to achieve minimal wall-clock overhead when using 8 groups of inducing points.
In CIFAR-10 classification,
we show that our method
can be integrated with deep convolutional structures to achieve state-of-the-art results for GPs.

\section{Background}
\label{sec:bg}

\subsection{Discrete Fourier Transform}
\label{sec:bg-dft}

Fourier analysis \citep{baron1878analytical} studies the representation of functions as sums or integrals of sinusoids.
For an integrable function $f$ on $\mathbb{R}^d$, its Fourier transform is defined as
\begin{align}
    \hat{f}(\bm{\omega}) = \int_{\mathbb{R}^d} e^{-2\pi i \bm{\omega}^\top\mathbf{x}} f(\mathbf{x})\, \mathrm{d}\mathbf{x},
\end{align}
where $\hat{f}(\bm{\omega}) \in \mathbb{C}$.
In kernel theory, Bochner's Theorem \citep{bochner1959lectures} is a seminal result that uses the Fourier transform to establish a bijection between stationary kernels and positive measures in the spectral domain.

Fourier analysis can be performed over finite sequences as well, via the discrete Fourier transform (DFT) \citep{cooley1969finite}. Specifically, given a sequence $\mathbf{x} = [\mathbf{x}_0, ..., \mathbf{x}_{T-1}]^{\top}$\!, DFT computes the sequence
$\hat{\mathbf{x}} = [\hat{\mathbf{x}}_0, ..., \hat{\mathbf{x}}_{T-1}]^{\top}$\!,  with
\begin{align}
   \hat{\mathbf{x}}_{t} = \frac{1}{T}\sum_{s=0}^{T-1} \mathbf{x}_s e^{-i\frac{2\pi t s}{T}}, \; t=0,...,T-1.
\end{align}
Let $\mathbf{F}:=[\frac{1}{T} e^{-i\frac{2\pi t s}{T} } ]_{t,s=0}^{T-1} \in \mathbb{C}^{T \times T}$ denote the DFT matrix.
The DFT can be represented in vector form as
$\hat{\mathbf{x}} = \mathbf{F} \mathbf{x}$, which naturally leads to the inverse DFT: $\mathbf{x}=\mathbf{F}^{-1}\hat{\mathbf{x}}$.

More generally, if $\mathbf{X} \in \mathbb{C}^{T_1 \times \cdots \times T_k}$ is a tensor, the multidimensional DFT computes the tensor $\hat{\mathbf{X}} \in \mathbb{C}^{T_1 \times \cdots \times T_k}$,
\begin{align}
    \hat{\mathbf{X}}[t_1, ..., t_k] = \frac{1}{T} \sum_{s_1, ..., s_k} X[s_1, ..., s_k] \prod_{j=1}^k e^{-i\frac{2\pi  t_j s_j}{T_j}}, \notag
\end{align}
where $T = \prod_{j}T_j$. The DFT matrix is then a tensor product of one-dimensional DFT matrices.

\subsection{Gaussian Processes}
Given an input domain $\mathcal{X}$, a mean function $m$, and a kernel function $k: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$, the Gaussian process \citep{rasmussen2006gaussian} $\mathcal{GP}(m, k)$ is a distribution
over functions $\mathcal{X} \to \mathbb{R}$.
For any finite set $\{\mathbf{x}_1,\mathbf{x}_2, \dots, \mathbf{x}_N\} \subset \mathcal{X}$, the function values $\mathbf{f} = [f(\mathbf{x}_1), f(\mathbf{x}_2), ..., f(\mathbf{x}_N)]^\top$ have a multivariate Gaussian distribution:
\begin{equation}
\mathbf{f} \sim \mathcal{N}(m(\mathbf{X}), \mathbf{K}_{\mathbf{f} \mathbf{f}}), \notag
\end{equation}
where $m(\mathbf{X}) =[m(\mathbf{x}_1), ..., m(\mathbf{x}_N)]^\top$\!, and $\mathbf{K}_{\mathbf{f} \mathbf{f}} = [k(\mathbf{x}_i,\mathbf{x}_j)]_{i,j=1}^N$.
For simplicity we assume $m(\cdot) = 0$ throughout the paper.
The observations $y$ are modeled with a density $p(y|f(\mathbf{x}))$, often taken to be Gaussian in the regression setting: $y = f(\mathbf{x}) + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2)$.
Let  $(\mathbf{X}, \mathbf{y})$ be a training set of size $N$. The posterior distribution $p(\mathbf{f}^{\star}|\mathbf{y})$ under a Gaussian observation model is
\begin{align}
\mathcal{N}( \mathbf{K}_{\star \mathbf{f}}(\mathbf{K}_{\mathbf{f} \mathbf{f}} + \sigma^2\mathbf{I})^{-1}\mathbf{y}, \mathbf{K}_{\star \star} - \mathbf{K}_{\star \mathbf{f}}(\mathbf{K}_{\mathbf{f} \mathbf{f}} + \sigma^2\mathbf{I})^{-1}  \mathbf{K}_{\mathbf{f} \star}), \notag
\end{align}
where $\mathbf{f}^{\star}=f(\mathbf{X}^{\star})$ are function values at test locations. Unfortunately, computing the posterior mean and covariance requires inverting the kernel matrix, an $\mathcal{O}(N^3)$ computation.

Sparse variational Gaussian processes (SVGPs) \citep{titsias2009variational, hensman2013gaussian} use inducing points for scalable GP inference. Let $\mathbf{Z} = [\mathbf{z}_1, ..., \mathbf{z}_M]^{\top} \in \mathbb{R}^{M \times d}$ be $M$ inducing locations, and $\mathbf{u} = f(\mathbf{Z})$.
SVGPs consider an augmented joint likelihood, $p(f(\cdot), \mathbf{u}) = p(f(\cdot)| \mathbf{u}) p(\mathbf{u})$, and a variational approximation $q(f(\cdot), \mathbf{u}) = p(f(\cdot) | \mathbf{u}) q(\mathbf{u})$, where $ q(\mathbf{u})=\mathcal{N}(\bm{\mu}, \mathbf{S})$ is a parameterized multivariate Gaussian with mean $\bm{\mu}$ and covariance $\mathbf{S}$.
The variational approximation is optimized by maximizing the variational lower bound:
\begin{align}
    \mathcal{L} \vcentcolon= \mathbb{E}_{q(\mathbf{f}, \mathbf{u})} [\log p(\mathbf{y}|\mathbf{f}, \mathbf{X})] - \KL{q(\mathbf{u})}{p(\mathbf{u})}.
\end{align}
Since $\log p(\mathbf{y}|\mathbf{f}, \mathbf{X}) = \sum_{i=1}^N \log p(\mathbf{y}_i|f(\mathbf{x}_i))$
admits stochastic optimization,
SVGPs reduce the computational cost to $\mathcal{O}(M^3 + M^2B)$, where $B$ is the minibatch size.

\section{Orthogonal Kernel Decomposition}\label{sec:kern-decomp}
In this section we establish the framework to conduct kernel Fourier analysis, and prove that the kernel Fourier series leads to the orthogonal kernel decomposition.

\textbf{$T$-Periodic Transformation $G$.} A function $G: \mathcal{X} \to \mathcal{X}$ is $T$-periodic if $T$ is the smallest integer such that,
\begin{align}
    \forall \mathbf{x} \in \mathcal{X}, G^T (\mathbf{x}) \vcentcolon= \overbrace{G \circ \cdots \circ G}^{T} (\mathbf{x}) = \mathbf{x},
\end{align}
So $T$ consecutive applications of $G$ is the identity operator. Throughout this paper, we use $G^t$ to represent $t$ consecutive applications of $G$ for $t \in \mathbb{N}$.

\textbf{$G$-Invariant kernels.} A kernel function $k: \mathcal{X} \times \mathcal{X} \to \mathbb{C}$ is $G$-invariant if,
\begin{align}
    \forall \mathbf{x}, \mathbf{x}' \in \mathcal{X}, k(G(\mathbf{x}), G(\mathbf{x}')) = k(\mathbf{x}, \mathbf{x}'),
\end{align}
For example, a polynomial kernel $k(\mathbf{x}, \mathbf{x}')=(\mathbf{x}^{\top}\mathbf{x}'+c)^t$ is invariant to the rotation transformations.

\textbf{Kernel Fourier Series.} We define $T$ functions $k_t: \mathcal{X} \times \mathcal{X} \to \mathbb{C}, t=0,...,T-1$, using the Fourier series of the sequence $[k(\mathbf{x}, G^0(\mathbf{x}')), ..., k(\mathbf{x}, G^{T-1}(\mathbf{x}'))]$. Specifically, the $t$-th function represents the $t$-th frequency harmonics from DFT,
\begin{align}
    k_t(\mathbf{x}, \mathbf{x}') =  \sum_{s=0}^{T-1} \mathbf{F}_{t, s} k(\mathbf{x}, G^s(\mathbf{x}')),
\end{align}
where $\mathbf{F} = [\frac{1}{T} e^{-i\frac{2\pi t s}{T} } ]_{t,s=0}^{T-1}$ is the DFT matrix.

\subsection{Orthogonal Kernel Decomposition}
In this subsection we present our main result, that the kernel Fourier series form an orthogonal kernel sum decomposition of the kernel $k$.

\begin{theorem}[Kernel Sum Decomposition.]
\begin{align}
    k(\mathbf{x}, \mathbf{x}') = \sum_{t=0}^{T-1} k_t(\mathbf{x}, \mathbf{x}').
\end{align}
Moreover, for $t=0,...,T-1$, $k_t$ is a Hermitian kernel.
\end{theorem}

\begin{proof}
The equality can be directly proven,
\begin{align}
   \sum_{t=0}^{T-1} k_t(\mathbf{x}, \mathbf{x}')
    &= \sum_{t=0}^{T-1} \sum_{s=0}^{T-1} \frac{1}{T}e^{-i\frac{2\pi st}{T}} k(\mathbf{x}, G^s(\mathbf{x}')) \notag \\
    &=\frac{1}{T}  \sum_{s=0}^{T-1} k(\mathbf{x}, G^s(\mathbf{x}')) \sum_{t=0}^{T-1} e^{-i\frac{2\pi st}{T} } \notag \\
    &= \frac{1}{T}  \sum_{s=0}^{T-1} k(\mathbf{x}, G^s(\mathbf{x}')) T\delta_{s} = k(\mathbf{x}, \mathbf{x}'). \notag
\end{align}
To prove that $k_t$ is a kernel, we claim that $k_t(\mathbf{x}, \mathbf{x}') = \mathbf{F}_{t,:}^H \mathbf{K} \mathbf{F}_{t,:}$, where $\mathbf{K} = [k(G^{s_1}(\mathbf{x}), G^{s_2}(\mathbf{x}'))]_{s_1, s_2=0}^{T-1}$. Expanding the matrix vector products,
\begin{align}
    &\frac{1}{T^2}\sum_{s_1=0}^{T-1}\sum_{s_2=0}^{T-1} e^{-i\frac{2\pi t}{T}(s_2-s_1) }k(G^{s_1}(\mathbf{x}), G^{s_1 + (s_2-s_1)}(\mathbf{x}')) \notag \\
    &=\frac{1}{T}\sum_{s_0=0}^{T-1} e^{-i\frac{2\pi t }{T}s_0 }k(\mathbf{x}, G^{s_0}(\mathbf{x}')) = k_t(\mathbf{x}, \mathbf{x}'), \notag
\end{align}
where we changed the variable $s_2-s_1$ to $s_0$. Furthermore, since $k$ is a kernel and $\mathbf{K}$ is positive semidefinite,  we conclude that $k_t(\mathbf{x}, \mathbf{x}') = \mathbf{F}_{t,:}^H \mathbf{K} \mathbf{F}_{t,:}$ is a kernel as well.
\end{proof}

Besides the kernel sum decomposition, we further show that the kernels $k_t, t=0,...,T-1$ are orthogonal to each other, as identified by the following lemma.

\begin{lemma}[Orthogonality]\label{lem:ortho} For any $0\leq t_1 \neq t_2 \leq T-1$, let $\mathcal{H}_k, \mathcal{H}_1, \mathcal{H}_2$ be the RKHSs corresponding to the kernel $k, k_{t_1}, k_{t_2}$, respectively. Then for any  $f \in \mathcal{H}_{1}$ and $g \in \mathcal{H}_{2}$, $\innerprod{f}{g}_{\mathcal{H}_k} = 0$.
\end{lemma}
The proof of this lemma can be found in Subsection~\ref{app:proof:ortho} in the appendix. Because the $\mathcal{H}_k$ inner product of $f \in \mathcal{H}_{1}, g \in \mathcal{H}_{2}$ is always zero, we immediately obtain that the RKHSs for $k_t$ are disjoint expect the function $f \equiv 0$.
\begin{proposition}[Disjoint] For $t=0,...,T-1$, let $\mathcal{H}_t$ be the RKHS associated with $k_t$, then for any $0\leq t_1 \neq t_2 \leq T-1$,
\begin{align}
    \mathcal{H}_{t_1} \cap \mathcal{H}_{t_2} = \{0\},
\end{align}
\end{proposition}

The kernel summation and the orthogonality of RKHSs, translate to the RKHS orthogonal sum decomposition as follows,
\begin{theorem}[RKHS Orthogonal Sum Decomposition.] The RKHS $\mathcal{H}_k$ admits an orthogonal sum decomposition,
\begin{align}
    \mathcal{H}_k = \bigoplus_{t=0}^{T-1} \mathcal{H}_{k_t},
\end{align}
Specifically, for any function $f \in \mathcal{H}_k $, $f$ has the unique decomposition $f = \sum_{t=0}^{T-1} f_t, f_t \in \mathcal{H}_{k_t}$,
\begin{align}
    f_t(\mathbf{x}) = \int k_t(\mathbf{x}, \mathbf{y}) f(\mathbf{y}) d\mathbf{x},
\end{align}
The RKHS norm of $f$ equals to,
\begin{align}
    \|f\|^2_{\mathcal{H}_k} = \sum_{t=0}^{T-1} \|f_t\|^2_{\mathcal{H}_{k_t}}.
\end{align}
\end{theorem}

\subsection{Examples of Kernel Fourier Series}\label{subsec:example}
\textbf{A Teaser Example.} We start with a toy problem for illustrating the kernel decomposition. We use $\theta \in [0, 2\pi)$ as the input, then $G(\theta) = (\theta + \frac{2\pi}{T}) \mod 2\pi$ is  $T$-periodic.
the kernel $k(\theta, \theta') = e^{-i(\theta - \theta')} + e^{-2i(\theta - \theta')}$ is then $G$-invariant. The kernel Fourier series can be explicitly obtained as, $k_{1}(\theta, \theta') = e^{-i(\theta - \theta')}$, $k_{2}(\theta, \theta') = e^{-2i(\theta - \theta')}$, and $k_t = 0$ otherwise. The RKHSs of $k_1$ and $k_2$ contain periodic functions of period $2\pi$ and $\pi$, respectively. Therefore, the kernel Fourier series decompose the original RKHS $\mathcal{H}_k$ into orthogonal RKHSs based on the frequencies.

\textbf{Inner-product kernels} are kernels of the form,
\begin{align}
    k(\mathbf{x}, \mathbf{x}') = \kappa(\mathbf{x}^H\mathbf{x}, \mathbf{x}^H\mathbf{x}', \mathbf{x}'^H\mathbf{x}'),
\end{align}
where the function $\kappa$ assures $k$ being positive semi-definite. For any matrix $\mathbf{R} \in \mathbb{C}^{d \times d}$, which is unitary $\mathbf{R} \mathbf{R}^{*} = \mathbf{I}$ and $\mathbf{R}^{T} = \mathbf{I}$, the mapping $G(\mathbf{x}) = \mathbf{R} \mathbf{x}$ is a $T$-periodic transformation. Because $\mathbf{R}$ is unitary, the kernel $k$ is $G$-invariant,
\begin{align}
    k(G(\mathbf{x}), G(\mathbf{x}')) &= \kappa(\mathbf{x}^H \mathbf{R}^H \mathbf{R} \mathbf{x}, \mathbf{x}^H \mathbf{R}^H \mathbf{R} \mathbf{x}', \mathbf{x}'^H \mathbf{R}^H \mathbf{R}  \mathbf{x}') \notag \\
    &= \kappa(\mathbf{x}^H\mathbf{x}, \mathbf{x}^H\mathbf{x}', \mathbf{x}'^H\mathbf{x}') \notag,
\end{align}
Therefore, we can compute the kernel Fourier series $k_t$. $k_t$ captures the $t$-th harmonic of the kernel $k$ along the trajectory of rotations $\mathbf{R}$.

\textbf{Stationary kernels} are kernels of the form,
\begin{align}
    k(\mathbf{x}, \mathbf{x}') = \kappa(\mathbf{x} - \mathbf{x}'),
\end{align}
where $\kappa$ is a positive-type function. Let $T=2$ and $G(\mathbf{x}) = -\mathbf{x}$, then $G$ is $T$-periodic. For real kernels whose $k(\mathbf{x}, \mathbf{x}') \in \mathbb{R}$, the kernel is symmetric $k(\mathbf{x}, \mathbf{x}')=k(\mathbf{x}', \mathbf{x})$. Then $k$ is $G$-invariant,
\begin{align}
    k(G(\mathbf{x}), G(\mathbf{x}')) = \kappa(\mathbf{x}' - \mathbf{x}) = k(\mathbf{x}', \mathbf{x}) = k(\mathbf{x}, \mathbf{x}'). \notag
\end{align}
Notably, the inner-product kernels are invariant to negations as well, thus $G(\mathbf{x})=-\mathbf{x}$ is applicable to both stationary kernels and inner-product kernels.

\subsection{Multi-dimensional Transformations}\label{subsec:multi-dim}
In the previous subsections we considered the Fourier series along one transformation trajectory: $k(\mathbf{x}, G^0(\mathbf{x}')), k(\mathbf{x}, G^1(\mathbf{x}')), ..., k(\mathbf{x}, G^{T-1}(\mathbf{x}'))$. Furthermore, we can extend the single transformation sequence to a multidimensional transformation tensor map, which would provide a tensor sum decomposition of the kernel. Specifically, let $T_1, ..., T_M \in \mathbb{N}$ and $G_m$ is a $T_m$-periodic transformation for $m=1,...,M$, respectively. We further assume that all transformations commute, i.e., $\forall 1\leq m_1, m_2 \leq M$,
\begin{align}
\forall \mathbf{x} \in \mathcal{X}, G_{m_1} (G_{m_2} (\mathbf{x})) =  G_{m_2} (G_{m_1} (\mathbf{x})) ,
\end{align}
Due to the commuting property, we can use indices $(t_1, ..., t_M)$ and the character $G$ to represent applying each $G_m$ by $t_m$ times,
\begin{align}
    G^{(t_1, ..., t_M)}(\mathbf{x}) \vcentcolon= G_1^{t_1}\cdots G_M^{t_M}(\mathbf{x}),
\end{align}
Now that we have defined the $(T_1, ..., T_M)$-periodic tensor transformation $G$. If a kernel $k$ is $G_m$-invariant for all $m=1,...,M$, then $k$ is $G$-invariant.

We define the kernel Fourier series correspondingly. Let $(t_1, ..., t_M)$ be the tensor index, we compute the kernel harmonics from multidimensional DFT,
\begin{align}
    &k_{(t_1, ..., t_M)}(\mathbf{x}, \mathbf{x}') \notag \\
    &= \sum_{s_1=0}^{T_1 - 1} \cdots \sum_{s_M=0}^{T_M - 1} \prod_{m=1}^M \mathbf{F}^{(m)}_{t_m, s_m} k(\mathbf{x}, G^{(s_1, ..., s_M)}(\mathbf{x}')), \notag
\end{align}
where $\mathbf{F}^{(m)} = [\frac{1}{T_m} e^{-i \frac{2\pi t_m s_m}{T_m}}]_{t_m, s_m=0}^{T_m-1}$ is the DFT matrix for the $m$-th transformation.

\subsection{Kernel Fourier Transform}
We have shown that the Fourier series of length $T$ decomposes the kernel into $T$ orthogonal kernels based on the frequencies. Intuitively, if we let $T \to \infty$, then we obtain a continuous frequency representation of the kernel, which would be akin to a Fourier transform. In this subsection, we present rigorous evidences to connect the kernel Fourier series with kernel Fourier transform.

Consider the transformation $G^\mathbf{s}: \mathcal{X} \to \mathcal{X}, \mathbf{s} \in \mathbf{R}^M$, where $M\geq 1$ corresponds to the multi-dimensional transformations. $G$ is periodic with period $\mathbf{1}$: $G^\mathbf{0}(\mathbf{x})=G^\mathbf{1}(\mathbf{x})=\mathbf{x}, G^{\mathbf{s}_1+\mathbf{s}_2}(\mathbf{x})=G^{\mathbf{s}_1}(G^{\mathbf{s}_2}(\mathbf{x})), \forall \mathbf{s}_1, \mathbf{s}_2 \in \mathbb{R}^M$. A kernel is $G$-invariant if for any $\mathbf{s} \in \mathbf{R}^M$, $k(G^\mathbf{s}(\mathbf{x}), G^\mathbf{s}(\mathbf{x}')) = k(\mathbf{x}, \mathbf{x}')$. Then we define the \emph{kernel Fourier transform} as the following: for $\mathbf{t} \in \mathbb{R}^M$,
\begin{align}
    k_{\mathbf{t}}(\mathbf{x}, \mathbf{x}') = \int_{[0,1]^M} e^{-2\pi i \mathbf{s}^{\top}\mathbf{t}} k(\mathbf{x}, G^\mathbf{s}(\mathbf{x}')) \mathrm{d}\mathbf{s},
\end{align}
Therefore, $k_\mathbf{t}(\mathbf{x}, \mathbf{x}')$ captures the frequency $\mathbf{t}$  of the function $\mathbf{s} \to k(\mathbf{x}, G^\mathbf{s}(\mathbf{x}'))$. Similar to kernel series, the kernel Fourier transform is an alternative representation of $k$, as shown in the following theorem.
\begin{theorem}[Kernel Frequency Representation]\label{thm:kern-ft}
\begin{align}
    k(\mathbf{x}, \mathbf{x}') = \int_{\mathbb{R}^M} k_\mathbf{t}(\mathbf{x}, \mathbf{x}') \mathrm{d} \mathbf{t}.
\end{align}
Moreover, $k_\mathbf{t}$ is a kernel for all $\mathbf{t} \in \mathbb{R}^M$.
\end{theorem}
The proof can be found in Subsection~\ref{app:subsec:proof-kern-ft} in the appendix.

\section{Kernel Spectral Analysis}
\section{Background}

Throughout the paper, we use lowercase characters $a$ to represent scalars, bold lowercase characters $\mathbf{a}$ to represent vectors, and bold uppercase characters $\mathbf{A}$ and to represent matrices. We use $\mathbf{A}^{\top}$ to represent transposes, and use $\mathbf{A}^H$ to represent the conjugate transposes.

\subsection{Fourier Transform and Fourier Series}
Fourier analysis \citep{baron1878analytical} is concerned in representing the functions as sums of harmonic series. For an integrable function $f$ with domain $\mathbb{R}^d$, its Fourier transform returns a function $\hat{f}$ of the frequence $\bm{\omega}$,
\begin{align}
    \hat{f}(\bm{\omega}) = \int_{\mathbb{R}^d} \exp(-2\pi i \bm{\omega}^\top\mathbf{x}) f(\mathbf{x}) \mathrm{d}\mathbf{x},
\end{align}
where $\hat{f}(\bm{\omega}) \in \mathbb{C}$, whose magnitude characterizes the amount of the frequency $\bm{\omega}$ present in $f$. The inverse Fourier transform maps $\hat{f}$ back to $f$, $f(\mathbf{x})  = \int_{\mathbb{R}^d} \exp(2\pi i \bm{\omega}^\top\mathbf{x}) \hat{f}(\bm{\omega}) \mathrm{d}\bm{\omega}$.

Fourier analysis can be performed over finite sequences as well, via the discrete Fourier transform (DFT). Specifically, given a finite sequence $\mathbf{x} = [\mathbf{x}_0, ..., \mathbf{x}_{N-1}]^{\top}$, DFT computes the Fourier series $\hat{\mathbf{x}} = [\hat{\mathbf{x}}_0, ..., \hat{\mathbf{x}}_{N-1}]^{\top}$,
\begin{align}
   \hat{\mathbf{x}}_{m} = \frac{1}{N}\sum_{n=0}^{N-1} \mathbf{x}_n e^{-i\frac{2\pi n m}{N}}, \; m=0,...,N-1,
\end{align}
Let $\mathbf{F} = [\frac{1}{N} e^{-i\frac{2\pi n m}{N} } ]_{n,m=0}^{N-1} \in \mathbb{C}^{N \times N}$, DFT can be written compactly as : $\hat{\mathbf{x}} = \mathbf{F} \mathbf{x}$. We can also compute the inverse discrete Fourier transform: $\mathbf{x} = \mathbf{F}^{-1}\hat{\mathbf{x}}$, where $\mathbf{F}^{-1} = [e^{i\frac{2\pi nm}{N}}]_{n,m=0}^{N-1}$. More generally, if $\mathbf{X} \in \mathbb{C}^{N_1 \times \cdots \times N_k}$ is a tensor, DFT computes $\hat{\mathbf{X}} \in \mathbb{C}^{N_1 \times \cdots \times N_k}$,
\begin{align}
    \hat{\mathbf{X}}[m_1, ..., m_k] = \frac{1}{N} \sum_{n_1, ..., n_k} X[n_1, ..., n_k] \prod_{j=1}^k e^{-i\frac{2\pi n_j m_j}{N_j}}, \notag
\end{align}
where $N = \prod_{j}N_j$ and $\hat{\mathbf{X}}$ is the frequency tensor of the same shape as $\mathbf{X}$.

\subsection{Gaussian Processes}
Given the input domain $\mathcal{X}$ and a kernel function $k: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$, a zero-mean Gaussian process \citep{rasmussen2006gaussian} $\mathcal{GP}(0, k)$ is a distribution $p(f)$ over functions $\mathcal{X} \to \mathbb{R}$. For any finite set $\{\mathbf{x}_1, ..., \mathbf{x}_n\} \subset \mathcal{X}$, the function values $\mathbf{f} = (f(\mathbf{x}_1), f(\mathbf{x}_2), ..., f(\mathbf{x}_n))$ have a multivariate Gaussian distribution.
\begin{equation}
\mathbf{f} \sim \mathcal{N}(\mathbf{0}, \mathbf{K}_{\mathbf{f} \mathbf{f}}), \notag
\end{equation}
where $\mathbf{K}_{\mathbf{f} \mathbf{f}} = [k(\mathbf{x}_i,\mathbf{x}_j)]_{i,j=1}^n$. The observation $y$ is generated by the likelihood $p(y|f)$. For regression problems, the likelihood usually models independent Gaussian noise: $y = f(\mathbf{x}) + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2)$. Let  $(\mathbf{X}, \mathbf{y})$ be the training set of size $N$, the posterior $p(\mathbf{f}^{\star}|\mathbf{y})$ is an explicit Gaussian under Gaussian likelihoods,
\begin{align}
\mathcal{N}( \mathbf{K}_{\star \mathbf{f}}(\mathbf{K}_{\mathbf{f} \mathbf{f}} + \sigma^2\mathbf{I}_N)^{-1}\mathbf{y}, \mathbf{K}_{\star \star} - \mathbf{K}_{\star \mathbf{f}}(\mathbf{K}_{\mathbf{f} \mathbf{f}} + \sigma^2\mathbf{I}_N)^{-1}  \mathbf{K}_{\mathbf{f} \star}), \notag
\end{align}
where $\mathbf{f}^{\star}=f(\mathbf{X}^{\star})$ are function values at test locations. However, the posterior involves inverting the kernel matrix, which requires $\mathcal{O}(N^3)$ computational costs.

Sparse variational Gaussian Processes (SVGP) \citep{titsias2009variational, hensman2013gaussian} use inducing point for scalable GP inference. Let $\mathbf{Z} = [\mathbf{z}_1, ..., \mathbf{z}_M]^{\top} \in \mathbb{R}^{M \times d}$ be $M$ inducing locations, and $\mathbf{u} = f(\mathbf{Z})$. SVGP considers an augmented joint likelihood, $p(f, \mathbf{u}) = p(f | \mathbf{u}) p(\mathbf{u})$, and a variational posterior $q(f, \mathbf{u}) = p(f | \mathbf{u}) q(\mathbf{u}), q(\mathbf{u})=\mathcal{N}(\bm{\mu}, \mathbf{S})$. The variational posterior is optimized by maximizing the Evidence Lower Bound (ELBO):
\begin{align}
    \mathcal{L} \vcentcolon= \mathbb{E}_{q(f, \mathbf{u})} [\log p(\mathbf{y}|f, \mathbf{X})] - \KL{q(\mathbf{u})}{p(\mathbf{u})}.
\end{align}
Since $\log p(\mathbf{y}|f, \mathbf{X}) = \sum_{i=1}^N \log p(\mathbf{y}_i|f, \mathbf{x}_i)$, the ELBO admits stochastic optimizations. Therefore, the SVGP reduces the computational costs from $\mathcal{O}(N^3)$ to $\mathcal{O}(M^3)$.

\section{Harmonic Kernel Decomposition}\label{sec:kern-decomp}
In this section, we introduce kernel Fourier series and use them
to form the harmonic kernel decomposition.
All proofs can be found in \Cref{app:proof:ortho}.

\subsection{Kernel Fourier Series}

We first propose a general method for representing a kernel as a sum of functions.
The idea is based on the DFT (see Sec.~\ref{sec:bg-dft}). Let $k: \mathcal{X}\times\mathcal{X} \to \mathbb{C}$ be a positive definite kernel.
To apply the DFT, we fix the first input $\mathbf{x}$, and construct a finite sequence of kernel values using a transformation $G: \mathcal{X}\to\mathcal{X}$ that applies to the second input:
\begin{equation} \label{eq:dft-seq}
	[k(\mathbf{x}, G^0(\mathbf{x}')), k(\mathbf{x}, G^1(\mathbf{x}')), ..., k(\mathbf{x}, G^{T-1}(\mathbf{x}'))],
\end{equation}
where $G^0(\mathbf{x}) := \mathbf{x}$ and $G^t := G \circ G^{t - 1}$.
Note that in signal processing, the sequence under DFT usually contains equally-spaced samples along the time domain. Here, we adopt a more general form by using $G$ to exploit symmetries in the input domain.

\begin{definition}[Kernel Fourier Series]
We define $T$ complex-valued functions $k_t: \mathcal{X} \times \mathcal{X} \to \mathbb{C}$, $t=0, ~\dots,~T-1$ using the DFT of the sequence \eqref{eq:dft-seq}:
\begin{align} \label{eq:def-kt}
	k_t(\mathbf{x}, \mathbf{x}') =  \sum_{s=0}^{T-1} \mathbf{F}_{t, s} k(\mathbf{x}, G^s(\mathbf{x}')),
\end{align}
where $\mathbf{F}$ is the DFT matrix (see Sec.~\ref{sec:bg-dft}).
The \emph{kernel Fourier series} of $k(\mathbf{x}, G^s(\mathbf{x}'))$ is given by the inverse DFT:
\begin{align} \label{eq:kernel-fourier-series}
	k(\mathbf{x}, G^s(\mathbf{x}')) = \sum_{t=0}^{T-1} \mathbf{F}_{s,t}^{-1} k_t(\mathbf{x}, \mathbf{x}').
\end{align}
The inverse DFT matrix is $\mathbf{F}^{-1} = T\mathbf{F}^H$, where $\cdot^H$ is the conjugate transpose.
\end{definition}

Given the positive definiteness of $k$, it is tempting to ask whether $k_t$ is also a (complex-valued) kernel.
In the next section, we will study the conditions when this holds and use it to form an orthogonal kernel decomposition.

\subsection{Harmonic Kernel Decomposition}\label{subsec:kernel-decom}
We first introduce the following definitions.

\begin{definition}[$T$-Cyclic Transformation]
A function $G: \mathcal{X} \to \mathcal{X}$ is $T$-cyclic if $T$ is the smallest integer such that,
\begin{align}
	\forall \mathbf{x} \in \mathcal{X},\;G^T (\mathbf{x}) \vcentcolon= \overbrace{G \circ \cdots \circ G}^{T} (\mathbf{x}) = \mathbf{x}.
\end{align}
\end{definition}
In group theory, $\{G^0, G^1, \dots, G^{T - 1}\}$ forms a cyclic group of order $T$, and $G$ is the generator of this group.
Interestingly, given a $T$-cyclic $G$, multiplying $k_t$ with $e^{i\frac{2\pi t}{T}}$ corresponds to a shift by $G$ in the second input.
\begin{proposition}[Shift]\label{prop:invariance} For any $0\leq t \leq T-1$,
\begin{align}
    k_t(\mathbf{x}, G(\mathbf{x}')) = e^{i\frac{2\pi t}{T}} k_t(\mathbf{x}, \mathbf{x}').
\end{align}
\end{proposition}
From \Cref{prop:invariance} we have $k_0(\mathbf{x}, G(\mathbf{x}')) = k_0(\mathbf{x}, \mathbf{x}')$, and when $T$ is even, $k_{T/2}(\mathbf{x}, G^2(\mathbf{x}')) = k_{T/2}(\mathbf{x}, \mathbf{x}')$.
This property is illustrated in \Cref{fig:rbf-rotate-gp}.

\begin{definition}[$G$-Invariant kernels]
A kernel function $k: \mathcal{X} \times \mathcal{X} \to \mathbb{C}$ is $G$-invariant if,
\begin{align}
	\forall \mathbf{x}, \mathbf{x}' \in \mathcal{X},\; k(G(\mathbf{x}), G(\mathbf{x}')) = k(\mathbf{x}, \mathbf{x}').
\end{align}
\end{definition}

For example, a polynomial kernel $k(\mathbf{x}, \mathbf{x}')=(\mathbf{x}^{\top}\mathbf{x}'+c)^t$ is invariant to the rotation transformations.

\begin{theorem}[Harmonic Kernel Decomposition] \label{thm:kern-decomp}
Let $G$ be a $T$-cyclic transformation, and $k$ be a $G$-invariant kernel. Then, the following decomposition holds:
\begin{align}
    k(\mathbf{x}, \mathbf{x}') = \sum_{t=0}^{T-1} k_t(\mathbf{x}, \mathbf{x}'),
\end{align}
where $k_t,\;t=0, \dots,T-1$ are defined as in Eq.~\eqref{eq:def-kt}. Moreover, for any $0\leq t \leq T-1$, $k_t$ is a Hermitian kernel.
\end{theorem}
The equation follows from the kernel Fourier series of $k(\mathbf{x}, G^0(\mathbf{x}'))$ by noticing that $\mathbf{F}^{-1}_{0,:} = \mathbf{1}$. To prove that $k_t$ is a kernel, we show that $k_t(\mathbf{x}, \mathbf{x}') = \mathbf{F}_{:, t}^H \mathbf{K}(\mathbf{x}, \mathbf{x}') \mathbf{F}_{:, t}$, where $\mathbf{K}(\mathbf{x}, \mathbf{x}') = [k(G^{s_1}(\mathbf{x}), G^{s_2}(\mathbf{x}'))]_{s_1, s_2=0}^{T-1}$,

Besides the kernel sum decomposition, we further show that the kernels $k_t,\;t=0,...,T-1$ are orthogonal to each other, as identified by the following lemma.

\begin{lemma}[Orthogonality]\label{lem:ortho} For any $0\leq t_1 \neq t_2 \leq T-1$, let $\mathcal{H}_k, \mathcal{H}_{k_{t_1}}, \mathcal{H}_{k_{t_2}}$ be the RKHSs corresponding to the kernel $k, k_{t_1}, k_{t_2}$, respectively. Then for any  $f \in \mathcal{H}_{k_{t_1}}$ and $g \in \mathcal{H}_{k_{t_2}}$, $\innerprod{f}{g}_{\mathcal{H}_k} = 0$.
\end{lemma}
Because the $\mathcal{H}_k$ inner product of $f \in \mathcal{H}_{1}, g \in \mathcal{H}_{2}$ is always zero, we immediately obtain that the RKHSs for $k_t$ are disjoint except for the function $f \equiv 0$.
\begin{proposition}[Disjoint]\label{prop:disjoint} For any $0\leq t_1 \neq t_2 \leq T-1$,
\begin{align}
    \mathcal{H}_{k_{t_1}} \cap \mathcal{H}_{k_{t_2}} = \{0\},
\end{align}
\end{proposition}

The kernel decomposition and orthogonality translate to the RKHS orthogonal sum decomposition as follows:
\begin{theorem}[Orthogonal Sum Decomposition of RKHS]\label{thm:rkhs-decom} The RKHS $\mathcal{H}_k$ admits an orthogonal sum decomposition,
\begin{align}
    \mathcal{H}_k = \bigoplus_{t=0}^{T-1} \mathcal{H}_{k_t}.
\end{align}
Specifically, for any function $f \in \mathcal{H}_k $, $f$ has the unique decomposition $f = \sum_{t=0}^{T-1} f_t,\; f_t \in \mathcal{H}_{k_t}$, and
$f_t(\mathbf{x}) = \innerprod{f}{k_t(\mathbf{x}, \cdot)}_{\mathcal{H}_k}$.
The RKHS norm of $f$ is equal to
\begin{align}
    \|f\|^2_{\mathcal{H}_k} = \sum_{t=0}^{T-1} \|f_t\|^2_{\mathcal{H}_{k_t}}.
\end{align}
\end{theorem}

\subsection{Examples of Harmonic Kernel Decomposition}\label{subsec:example}
The HKD relies on the $(G, k)$ pair where $G$ is a $T$-cyclic transformation and $k$ is a kernel invariant to $G$.
In this section we provide examples of such kernels and transformations.
Notably, all inner-product kernels and stationary kernels\footnote{This includes, e.g., polynomial, Gaussian, Mat\'ern, periodic, arccosine, and rational quadratic kernels.} can be decomposed with the HKD when paired with an appropriately chosen $G$.

\paragraph{An opening example.} We start with a toy example to illustrate the kernel decomposition. Let $k(\theta, \theta') = e^{-i(\theta - \theta')} + e^{-2i(\theta - \theta')}$ for $\theta \in [0, 2\pi)$. The transformation $G(\theta) = (\theta + \frac{2\pi}{T}) \mod 2\pi$ is  $T$-cyclic. Based on the kernel Fourier series, we obtain $k_{1}(\theta, \theta') = e^{-i(\theta - \theta')}$, $k_{2}(\theta, \theta') = e^{-2i(\theta - \theta')}$, and $k_t = 0$ otherwise. We observe that the RKHS of $k_2$ contains periodic functions with basic period $\pi$, while the RKHS of $k_1$ contains functions with basic period $2\pi$.
In this way, our method decomposes $\mathcal{H}_k$ into orthogonal RKHSs.

\textbf{Inner-product kernels} are kernels of the form,
\begin{align}
    k(\mathbf{x}, \mathbf{x}') = h(\mathbf{x}^H\mathbf{x}, \mathbf{x}^H\mathbf{x}', \mathbf{x}'^H\mathbf{x}'),
\end{align}
where the function $h$ ensures that $k$ is positive semi-definite \citep{hofmann2008kernel}. For a matrix $\mathbf{R} \in \mathbb{C}^{d \times d}$, which is unitary (i.e.~$\mathbf{R} \mathbf{R}^{H} = \mathbf{I}$),
the kernel $k$ is $G$-invariant:
\begin{align}
    \!k(G(\mathbf{x}), G(\mathbf{x}'))\!&=\!h(\mathbf{x}^H \mathbf{R}^H \mathbf{R} \mathbf{x}, \mathbf{x}^H \mathbf{R}^H \mathbf{R} \mathbf{x}', \mathbf{x}'^H \mathbf{R}^H \mathbf{R}  \mathbf{x}') \notag \\
    \!&=\! h(\mathbf{x}^H\mathbf{x}, \mathbf{x}^H\mathbf{x}', \mathbf{x}'^H\mathbf{x}') = k(\mathbf{x}, \mathbf{x}') \notag,
\end{align}
Examples include reflections, rotations, and permutations.
Moreover, if $\underbrace{\mathbf{R} \cdots \mathbf{R}}_T = \mathbf{I}$,
the mapping $G(\mathbf{x}) = \mathbf{R} \mathbf{x}$ is  $T$-cyclic.

\textbf{Stationary kernels} are kernels of the form,
\begin{align}
    k(\mathbf{x}, \mathbf{x}') = \kappa(\mathbf{x} - \mathbf{x}'),
\end{align}
where $\kappa$ is a positive-type function \citep{berlinet2011reproducing}. Let $T=2$ and $G(\mathbf{x}) = -\mathbf{x}$; then $G$ is $T$-cyclic. For real kernels whose $k(\mathbf{x}, \mathbf{x}') \in \mathbb{R}$, the kernel is symmetric (i.e.~$k(\mathbf{x}, \mathbf{x}')=k(\mathbf{x}', \mathbf{x})$). Then $k$ is $G$-invariant:
\begin{align}
    k(G(\mathbf{x}), G(\mathbf{x}')) = \kappa(\mathbf{x}' - \mathbf{x}) = k(\mathbf{x}', \mathbf{x}) = k(\mathbf{x}, \mathbf{x}'). \notag
\end{align}
Similarly, we can prove that inner-product kernels are negation-invariant.
Stationary kernels are also invariant to the $T$-cyclic transformation $G_i(\mathbf{x}) = \mathbf{x} + \frac{2\pi}{T} \mathbf{e}_i$ on a multidimensional torus
$\mathbb{T}^d = \underbrace{\mathbb{S}^1 \times \cdots \times \mathbb{S}^1}_d$, where $\mathbb{S}^1$ represents a one-dimensional circle.

\subsection{Resolving Complex-Valued Kernels}\label{subsec:resolve-hermitian}

From $\mathbf{F} = [\frac{1}{T} e^{-i\frac{2\pi t s}{T} } ]_{t,s=0}^{T-1}$, we know that the DFT introduces complex values whenever $T > 2$.
Therefore, $k_t$ is Hermitian but not necessarily real-valued.
For example, the decomposition in Fig.~\ref{fig:rbf-rotate-gp} introduces imaginary values when $t = 1, 3$.
Since $k$ is real-valued, we can obtain a real-valued kernel decomposition by pairing up $k_t$s.
Specifically, for a $T$-cyclic transformation $G$, we have
\begin{align}
	(k_t + k_{T-t})(\mathbf{x}, \mathbf{x}') = \frac{2}{T}\sum_{s=0}^{T-1} \cos(\frac{2\pi t s}{T}) k(\mathbf{x}, G^s(\mathbf{x}')), \notag
\end{align}
In this way we obtain a real-valued decomposition with $\lfloor T/2 \rfloor+1$ kernels.

\subsection{Multi-Way Transformations}\label{subsec:multi-dim}

Previously we considered the Fourier series along one transformation orbit: $k(\mathbf{x}, G^0(\mathbf{x}')), ..., k(\mathbf{x}, G^{T-1}(\mathbf{x}'))$. We can extend it to multi-way transformations, akin to a multidimensional DFT. Let $T_1, ..., T_J \in \mathbb{N}$ and $G_j$ be a $T_j$-cyclic transformation for $j=1,...,J$, respectively. We further assume that all transformations commute, i.e., $\forall 1\leq j_1, j_2 \leq J$,
\begin{align}
\forall \mathbf{x} \in \mathcal{X}, G_{j_1} (G_{j_2} (\mathbf{x})) =  G_{j_2} (G_{j_1} (\mathbf{x})).
\end{align}
Due to commutativity, we can use the indices $(t_1, ..., t_J)$ to represent applying each $G_j$ for $t_j$ times,
\begin{align}
    G^{(t_1, ..., t_J)}(\mathbf{x}) \vcentcolon= G_1^{t_1}\cdots G_J^{t_J}(\mathbf{x}),
\end{align}
where $G:= G_1 \otimes \cdots \otimes G_J$. Moreover, if a kernel $k$ is $G_j$-invariant for all $j=1,...,J$, then $k$ is $G$-invariant.

Letting $t=(t_1, ..., t_J)$ be a multi-index, we compute the $J$-way kernel Fourier series from a multidimensional DFT:
\begin{align}
    k_{t}(\mathbf{x}, \mathbf{x}')
    = \sum_{s=(0, ..., 0)}^{(T_1 - 1, ..., T_J-1)} \prod_{j=1}^J \mathbf{F}^{(j)}_{t_j, s_j} k(\mathbf{x}, G^{s}(\mathbf{x}')), \notag
\end{align}
where $\mathbf{F}^{(j)} \in \mathbb{C}^{T_j \times T_j}$ is the DFT matrix of
the $j$-th transformation.
Similar to \Cref{thm:kern-decomp}, these $k_t$s also form an HKD of $k$.

Taking a two-dimensional RBF kernel as an example, we can check that it is invariant to negation along either dimension: $G_1([x_1, x_2]^{\top}) = [-x_1, x_2]^{\top}, G_2([x_1, x_2]^{\top}) = [x_1, -x_2]^{\top} $. Because $G_1$ and $G_2$ commute, this forms a 2-way transformation $G = G_1 \otimes G_2$, which corresponds to an HKD with $2 \times 2 = 4$ sub-kernels.

\section{Spectrally Decomposed Variational Gaussian Processes}
In this section, we explore the implications of the kernel orthogonal decomposition, and propose the \emph{Spectrally Decomposed Variational Gaussian Processes (SDVGP)}, an inference strategy for stochastic variational Gaussian processes.

\subsection{Decoupled Variational Inference}
Given the kernel decomposition $k = \sum_{t=0}^{T-1} k_t$ \footnote{$t$ can be a multi-index for multi-dimensional transformations.}, the Gaussian Processes $\mathcal{GP}(0, k)$ can be reformulated as,
\begin{align}\label{eq:additive-gp}
    f = \sum_{t=0}^{T-1} f_t, \; f_t \sim \mathcal{GP}(0, k_t), t=0,...,T-1
\end{align}
Given the GP decomposition, SDVGP uses a variational posterior separately for each GP, similarly to \citet{shi2020sparse}. Specifically, for $t=0,...,T-1$, let $\mathbf{u}_t$ be inducing variables on $\mathbf{Z}_t$ and $q_t(\mathbf{u}_t) = \mathcal{N}(\bm{\mu}_t, \mathbf{S}_t)$ be the variational posterior of $\mathbf{u}_t$, then the overall variational posterior can be formulated as,
\begin{align} \label{eq:q-k}
    f &= \sum_{t=0}^{T-1} f_t, \\
    q_t(f_t, \mathbf{u}_t) &= p_t(f_t | \mathbf{u}_t) q_t(\mathbf{u}_t), t=0,...,T-1, \notag
\end{align}
where $p_t$ represents the Gaussian Process $\mathcal{GP}(0, k_t)$. Specifically, if each $q_t$ uses $m$ inducing points, we term the model as a $T\times m$ model. The variational posterior can be optimized by maximizing the Evidence Lower Bound (ELBO),
\begin{align}
    \mathcal{L} \vcentcolon= \mathbb{E}_{q(f, \mathbf{u})} [\log p(\mathbf{y}|f, \mathbf{X})] & - \sum_{t=0}^{T-1} \KL{q_t(\mathbf{u}_t)}{p(\mathbf{u}_t)}. \notag
\end{align}

\subsection{Approximation Capacity}\label{subsec:appr-capa}
How well does the decoupled variational inference approximate the true posterior, compared to the standard SVGP? This question boils down to two questions,
\begin{itemize}
\item \textbf{Kernel Approximation.} How well do decoupled inducing points $\{\mathbf{Z}_t\}_{t=0}^{T-1}$ approximate the kernel $k$ ?
\item \textbf{Posterior Approximation.} How well does the decoupled variational posterior $\prod_{t=0}^{T-1} q_t(\mathbf{u}_t)$ approximate the optimial posterior $q^{\star}(\{\mathbf{u}_t\}_{t=0}^{T-1})$?
\end{itemize}

\textbf{Kernel Approximation.} Intuitively, how well does the Nystrom approximation using $\mathbf{Z}$ approximates the kernel, depends on how large is the space spanned by $k(\mathbf{x}, \mathbf{Z})$. For decoupled inducing points, the space is then $\text{span}\{k_t(\mathbf{x}, \mathbf{Z}_t)\}_{t=0}^{T-1}$. We let $\mathbf{Z}_0 =  \cdots \mathbf{Z}_{T-1} =: \tilde{\mathbf{Z}}$, since the Fourier series is a bijective linear operation,
\begin{align}
    \text{span}\{k_t(\mathbf{x}, \mathbf{Z}_t)\}_{t=0}^{T-1} = \text{span}\{k(\mathbf{x}, G^{t}(\tilde{\mathbf{Z}}))\}_{t=0}^{T-1}, \notag
\end{align}
Therefore, the decoupled inducing points $\{\mathbf{Z}_t\}_{t=0}^{T-1}$ approximate the kernel at least as well as the standard inducing points $\{G^t(\tilde{\mathbf{Z}})\}_{t=0}^{T-1}$. If the transformation $G$ is set so that $\{G^t(\tilde{\mathbf{Z}})\}_{t=0}^{T-1}$ covers the input domain well, then the decoupled inducing points approximate the kernel accurately.

\textbf{Posterior Approximation.} Given the decoupled inducing points, the optimal posterior covariance is $\mathbf{S}^{\star} \in \mathbb{R}^{TM \times TM}$, while we induce a block diagonal structure in SDVGP. However, we can show that $\mathbf{S}^{\star}$ is approximately block diagonal indeed if the input distribution is invariant to $G$.

\begin{lemma}\label{lem:decoupled} If the data distribution is invariant to $G$, i.e.,
\begin{align}
    \forall \mathbf{x}, p(G(\mathbf{x})) = p(\mathbf{x}),
\end{align}
Then $\mathbf{S}^{\star}$ becomes block diagonal when the training size $N \to \infty$.
\end{lemma}
The proof can be found in Subsection~\ref{app:prof:lemma-dvi} in the appendix. This lemma highlights that, if the density of the data distribution is approximately invariant to the transformation $G$, then the decoupled variational posterior induces little error compared to the optimal posterior.

Summing up the analyses in kernel approximation and posterior approximation, we conclude that the SDVGP approximates the true GP posterior accurately when \emph{the trajectories under the transformation $G$ cover the input domain well}. We illustrate this gist by a \textit{flip-mnist} dataset, where each digit in the MNIST dataset is randomly flipped up-and-down or left-and-right. For a RBF kernel, we consider two variations of SDVGP in terms of the transformation. 1) Negation. We split input dimensions into two groups and negate them separately, resulting in a $4\times 100$ model. 2) Flipping the image up-and-down or left-and-right, resulting in a $4\times 100 $ model. We compare them with the standard SVGP using $100, 400$ inducing points, shown in Figure~\ref{fig:flipmnist}.

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{figures/flipmnist_M100_reordered_two.pdf}
\caption{Flip-mnist.  From left to right, we plot how \textit{train loss}, \textit{test acc} and \textit{test lld} evolve with training. We observe SDVGP using \textit{FLIP: $4 \times 100$} performs on par with standard SVGP using \textit{$400$} inducing points while SDVGP using \textit{NEG: $4 \times 100$} performs substantially worse. This experiment highlights the importance of matching the transformation $G$ with the data distribution. \label{fig:flipmnist}}
\end{figure}

\subsection{Computational Costs}

\textbf{Computational Costs.} The computational costs boil down to the cost of
computing $k_t(\mathbf{Z}_t, \mathbf{Z}_t)$ and the cost of decoupled variational inference. Let $T\times m$ represents using $m$ inducing points for each $t=0, ..., T-1$, and assume the mini-batch size $n = \mathcal{O}(m)$. The computational cost of decoupled variational inference is then $\mathcal{O}(T m^3)$; To compute $k_t(\mathbf{Z}_t, \mathbf{Z}_t)$, we need the kernel values $k(\mathbf{Z}_t, G^s(\mathbf{Z}_t))$ for $s=0,...,T-1$. If $c_G$ represents the cost of applying $G$, then computing $\mathbf{K}_{\mathbf{u}, \mathbf{u}}$ requires $\mathcal{O}(Tm \times T c_G + Tm^2 \times T)$ operations. Therefore, the overall computational cost is
\begin{align}
    \mathcal{O}(T m^3  + T^2 m^2 + T^2 m c_G),
\end{align}
When $T \ll m$, the cost scales approximately linearly with $T$. Furthermore, the proposed SDVGP supports parallel computations straightforwardly: we can locate all computations of $k_t$ on the $t$-th device and enjoy the parallelism.

\textbf{Memory Costs.}  The memory cost overhead stems from computing the kernel matrix $\mathbf{K}_{\mathbf{u}\bu}$. For computing $\mathbf{K}_{\mathbf{u}\bu}$, we need the kernel values $k(\mathbf{Z}_t, G^s(\mathbf{Z}_t))$ for $s=0,...,T-1$, which implies the memory cost $\mathcal{O}(T \times m^2 \times T)$. Since the memory cost of decoupled variational inference is $\mathcal{O}(T \times m^2)$, the overall memory cost is,
\begin{align}
    \mathcal{O}(T^2m^2),
\end{align}

\textbf{Storage Costs.} Storing a $T\times m$ model needs to keep $T$ variational posteriors, requiring $\mathcal{O}(Tm^2 + Tmd)$ storage.

\subsection{Resolving Hermitian Kernels}\label{subsec:resolve-hermitian}
Our orthogonal kernel decomposition arises from the discrete Fourier transform, thus $k_t$ is Hermitian but not necessarily real-valued. This elicits the difficulty for $k_t$ to be used in Gaussian processes as a GP is commonly defined as a real process \citep{rasmussen2006gaussian}. So when is $k_t$ not real-valued and how can we resolve this issue ?

For a single dimensional transformation, the DFT matrix $\mathbf{F} = [\frac{1}{T} e^{-i\frac{2\pi t s}{T} } ]_{t,s=0}^{T-1}$. Then the  DFT introduces imaginary values whenever $T > 2$. For $M$-dimensional transformations where $t = (t_1, ..., t_M)$ is a multi-index, subsection~\ref{subsec:multi-dim} shows that the DFT matrix is the tensor product of all single-dimensional DFT matrices. Thus the DFT introduces imaginary values whenever $T_m > 2$ for some $1\leq m\leq M$. Therefore, rotations by $90^{\text{o}}$ in Figure~\ref{fig:rbf-rotate-gp} does introduce imaginary values because its period is $4$, but flipping left-right and up-down in Figure~\ref{fig:flipmnist} does not because it involves only period $2$ transformations.

Fortunately we can resolve Hermitian kernels to be all real-valued. Consider a single dimensional transformation with period $T$, we sum up $k_t$ and $k_{T-t}$ together,
\begin{align}
    (k_t + k_{T-t})(\mathbf{x}, \mathbf{x}') = \sum_{s=0}^{T-1} \left(\mathbf{F}_{t, s}+\mathbf{F}_{T-t, s} \right) k(\mathbf{x}, G^s(\mathbf{x}')), \notag
\end{align}
Because $\mathbf{F}_{t, s}+\mathbf{F}_{T-t, s} = \frac{2}{T} \cos(\frac{2\pi t s}{T})$, $k_t + k_{T-t}$ is a real-valued kernel. By pairing up all Hermitian kernels, we obtain a real-valued orthogonal kernel decomposition,
\begin{align}
    k(\mathbf{x}, \mathbf{x}') = \sum_{t=0}^{\lfloor T/2 \rfloor+1} \left(k_t + k_{T-t}\right)(\mathbf{x}, \mathbf{x}') .
\end{align}
Similarly, for $M$-dimensional transformations with periods $(T_1, ..., T_M)$, we can obtain a real-valued decomposition with $(\lfloor T_1/2 \rfloor+1) \times \cdots \times (\lfloor T_M/2 \rfloor+1)$ kernels.

\section{Related Works}

The idea of applying Fourier analysis to kernels goes back at least to \citet{bochner1959lectures}.
In machine learning, this led to a flowering of large-scale kernel methods based on random features~\citep{rahimi2008random,yu2016orthogonal,dao2017gaussian}.
Bochner's theorem also allows designing stationary kernels by modeling a spectral density~\citep{wilson2013gaussian,samo2015generalized,parra2017spectral,benton2019function}.
On hyperspheres, zonal kernels are the counterpart of stationary kernels.
Their spectral decomposition is given by spherical harmonics~\citep{thomson1888treatise, morimoto1998analytic}.
Although closely related, none of these works have considered the discrete Fourier transform adopted in our method.

HVGPs share many similarities with the works that propose decoupled~\citep{cheng2017variational, salimbeni2018orthogonally} and orthogonal~\citep{shi2020sparse} inducing points.
In particular, \citet{shi2020sparse} is also based on an orthogonal decomposition of the kernel and uses distinct groups of inducing points for them.
However, their decomposition
involves matrix inversion while ours can be computed using fast Fourier transforms.

Structured Kernel Interpolation (SKI) \citep{wilson2015kernel,wilson2015thoughts, evans2018scalable, izmailov2018scalable} places inducing points on a grid, leading to a structured $\mathbf{K}_{\mathbf{u}\bu}$ that allows fast matrix-vector multiplications.
For one-dimensional data, SKI exploits the Toeplitz structure of $\mathbf{K}_{\mathbf{u}\bu}$ generated by stationary kernels.
They first embed the Toeplitz matrix into a circulant matrix $\mathbf{C}$, and use the fact that circulant matrices can be diagonalized by the DFT~\citep{tee2007eigenvectors} to enable fast computations:
\begin{align}
	\mathbf{C} = \mathbf{F}^{-1} \mathrm{diag}(\mathbf{F} \mathbf{c}) \mathbf{F}.
\end{align}
Here $\mathbf{F}$ is the DFT matrix, and $\mathbf{c}$ is the first column of $\mathbf{C}$.
This equation highlights a connection with our HKD:
If we let $\mathbf{C} = [k(G^{t_1}(\mathbf{x}'), G^{t_2}(\mathbf{x}))]_{t_1, t_2=0}^{T-1}$, then $\mathbf{c} = [k(\mathbf{x}, G^0(\mathbf{x}')), ..., k(\mathbf{x}, G^{T-1}(\mathbf{x}'))]$ is the sequence we constructed in Eq.~\eqref{eq:dft-seq},
and the eigenvalues $\mathbf{F} \mathbf{c}$
recover our decomposition $[k_t(\mathbf{x}, \mathbf{x}')]_{t=0}^{T-1}$ by the discrete Fourier transform.
In other words, our approach generalizes the structure of one-dimensional equally-spaced grids in SKI into arbitrary cyclic groups.
Moreover, our method allows trainable inducing locations, which plays an important role in combating the curse of dimensionality.

Besides inducing points in the data space, a number of works have investigated inducing features in the frequency domain.
However, these inducing features are either limited to specific kernels \citep{lazaro2009inter,hensman2017variational} or involve numerical approximations \citep{dutordoir2020sparse,burt2020variational}.
The implementation of \citet{dutordoir2020sparse} only supports data up to 8 dimensions.

Incorporating invariances with respect to input-space transformations into Gaussian processes is also investigated in a stream of works \citep{ginsbourger2016degeneracy, van2019learning}.
Our work is orthogonal to them since we are not designing invariant models.
Instead, we proposed a general inference method for GPs that can benefit from invariances in the data distribution.
Relatedly, \citet{solin2020hilbert, borovitskiy2020} studied Gaussian processes on Riemannian manifolds.

\section{Harmonic Variational Gaussian Processes}\label{sec:hvgp}
In this section, we explore the implications of the HKD, and propose a scalable inference strategy for variational Gaussian processes. All proofs can be found in Sec~\ref{app:prof:lemma-dvi} in the appendix.

\subsection{Variational Inference for Decomposed GPs }\label{subsec:hf-gp}
Given the kernel decomposition\footnote{$t$ can be a multi-index for multi-way transformations.} $k = \sum_{t=0}^{T-1} k_t$, the Gaussian process can be represented in an additive formulation,
\begin{align}\label{eq:additive-gp}
    f = \sum_{t=0}^{T-1} f_t, \quad f_t \sim \mathcal{GP}(0, k_t).
\end{align}
For $t=0,...,T-1$, we introduce inducing points $\mathbf{Z}_t$ and denote by $\mathbf{u}_t := f_t(\mathbf{Z}_t)$ the inducing variables. Let $p_t$ represent $\mathcal{GP}(0, k_t)$. We consider an augmented model,
\begin{align}\label{eq:joint-ll}
    f &= \sum_{t=0}^{T-1} f_t, \quad
    p_t(f_t(\cdot), \mathbf{u}_t) = p_t(f_t(\cdot) | \mathbf{u}_t) p_t(\mathbf{u}_t),
\end{align}
and define the variational posterior approximation as
\begin{equation}
\begin{aligned} \label{eq:full-q-h}
    f &= \sum_{t=0}^{T-1} f_t, \\
    f_t(\cdot) &\sim p_t(f_t(\cdot) | \mathbf{u}_t), \; \mathbf{u}_{0:T-1} \sim q(\mathbf{u}_{0:T-1}).
\end{aligned}
\end{equation}
To understand how well this variational distribution approximates the true GP posterior,
we compare it with a standard SVGP, for which
the quality of approximation has been studied extensively by \citet{burt2019rates}.
For simplicity, we focus our analysis on the case where inducing points are shared across all subprocesses: $\mathbf{Z}_0 = ... = \mathbf{Z}_{T-1} := \mathbf{Z}$, and we assume a complex-valued kernel decomposition without using the techniques in \Cref{subsec:resolve-hermitian}.
Then, we demonstrate that Eq.~\eqref{eq:full-q-h}
is equivalent to an SVGP with inducing points $\{G^t(\mathbf{Z})\}_{t=0}^{T-1}$.
\begin{theorem}\label{lem:q-equi} Consider an SVGP with inducing points $\{G^t(\mathbf{Z})\}_{t=0}^{T-1}$.
Let $\mathbf{v}_t := f(G^t(\mathbf{Z}))$ be the inducing variables and $\mathbf{V}:=[\mathbf{v}_0, ..., \mathbf{v}_{T-1}]^{\top} \in \mathbb{C}^{T \times m}$.
Suppose its variational distribution is
\begin{align}
    q_{\mathrm{svgp}}(f(\cdot), \mathbf{V}) &= p(f(\cdot)|\mathbf{V}) \mathcal{N}(\vectorize{\mathbf{V}}| \vectorize{\mathbf{M}_{v}}, \mathbf{S}_v),  \notag
\end{align}
where $\mathbf{M}_v \in \mathbb{C}^{T \times m}$, $\mathbf{S}_v \in \mathbb{C}^{Tm \times Tm}$ are the mean and covariance, respectively.
Let $\mathbf{U}:=[\mathbf{u}_0, ..., \mathbf{u}_{T-1}]^{\top} \in \mathbb{C}^{T \times m}$.
Then, Eq.~\eqref{eq:full-q-h}
and
$q_{\mathrm{svgp}}$ have the same marginal distribution of $f(\cdot)$ if $q(\mathbf{u}_{0:T-1})$ is defined as
\begin{align}
    q(\vectorize{\mathbf{U}}) = \mathcal{N}(\vectorize{\mathbf{F}^H \mathbf{M}_v}, (\mathbf{I} \otimes \mathbf{F}^H) \mathbf{S}_v (\mathbf{I} \otimes \mathbf{F})). \notag
\end{align}

\end{theorem}
The proof is based on showing the bijective linearity $\mathbf{U} = \mathbf{F}^H \mathbf{V}$.
Since the theorem assumes shared inducing points, our variational approximation in Eq.~\eqref{eq:full-q-h}
has a larger capacity than
SVGPs with inducing points $\{G^t(\mathbf{Z})\}_{t=0}^{T-1}$.
Therefore, if the inducing points $\{G^t(\mathbf{Z})\}_{t=0}^{T-1}$ match the input distribution well, our variational posterior can approximate the true posterior accurately.

\subsection{Harmonic Variational Gaussian Processes}\label{subsec:hvgp}

\begin{figure}[t]
\centering
\includegraphics[width=\columnwidth]{figures/flipmnist_M100_reordered_two_long.pdf}
\caption{Flip-MNIST.  We plot how \textit{left}: train loss and \textit{right}: test accuracy evolve with training. We observe the HVGP using \textit{FLIP: $4 \times 100$} perform similarly with the SVGP using \textit{$400$} inducing points while the HVGP using \textit{NEG: $4 \times 100$} perform substantially worse. \label{fig:flipmnist}}
\end{figure}

The additive GP reformulation in Eq.~\eqref{eq:additive-gp} ensures the independence between $\mathbf{u}_{0:T-1}$ in the prior, and we demonstrated the orthogonality of the decomposed RKHSs in Theorem~\ref{thm:rkhs-decom}. Thus it is tempting to modelling the variational posterior separately within each RKHS. Now we introduce the Harmonic Variational Gaussian Process (HVGPs), which enforces independence between $\mathbf{u}_t$ by letting $q(\mathbf{u}_{0:T-1}) = \prod_{t=0}^{T-1} q_t(\mathbf{u}_t)$. Then the variational posterior becomes
\begin{align} \label{eq:diag-q-h}
    f &= \sum_{t=0}^{T-1} f_t, \;
    q_t(f_t(\cdot), \mathbf{u}_t) = p_t(f_t(\cdot) | \mathbf{u}_t) q_t(\mathbf{u}_t).
\end{align}
In other words, HVGPs use a variational posterior independently for each GP. We set $q_t(\mathbf{u}_t) = \mathcal{N}(\bm{\mu}_t, \mathbf{S}_t)$ as Gaussians. In particular, if each $q_t$ uses $m$ inducing points, we term the model a $T\times m$ model. The variational posterior can be optimized by maximizing the ELBO,
\begin{align*}

    \resizebox{\columnwidth}{!}{$\displaystyle
    \mathbb{E}_{q(\{f_t\}^{T-1}_{t=0})} [\log p(\mathbf{y}|\sum_{t=0}^{T-1}f_t, \mathbf{X})]  - \sum_{t=0}^{T-1} \KL{q_t(\mathbf{u}_t)}{p_t(\mathbf{u}_t)}.$}
\end{align*}
where $q(\{f_t(\cdot)\}^{T-1}_{t=0}) := \prod_{t=0}^{T - 1} \int p_t(f_t(\cdot)|\mathbf{u}_t)q_t(\mathbf{u}_t) d\mathbf{u}_t$.

How well does $\prod_{t=0}^{T-1} q_t(\mathbf{u}_t)$  approximate the optimal Gaussian variational posterior $q^{\star}(\mathbf{u}_{0:T-1})$?
$q^{\star}$ has a covariance $\mathbf{S}^{\star} \in \mathbb{R}^{Tm \times Tm}$, while HVGPs induce block diagonal structures.
Fortunately, we can show that $\mathbf{S}^{\star}$ is approximately block diagonal if the input distribution is invariant to $G$.

\begin{theorem}\label{lem:decoupled} If the input distribution $p$ is invariant to $G$, i.e., the random variable $\mathbf{x} \sim p$ and the random variable $G(\mathbf{x}), \mathbf{x} \sim p$ are identically distributed, then
$\mathbf{S}^{\star}$ becomes block diagonal when the training size $N \to \infty$.
\end{theorem}

Theorem~\ref{lem:decoupled} indicates that the independent variational distributions $\prod_{t=0}^{T-1} q_t(\mathbf{u}_t)$ in HVGPs accurately approximate $q^{\star}(\mathbf{u}_{0:T-1})$
when the input distribution
is symmetric under the transformation $G$. The symmetry further makes it easy for $\{G^t(\mathbf{Z})\}_{t=0}^{T-1}$ to match the input distribution, which by Theorem~\ref{lem:q-equi} renders that the variational approximation in Eq.~\eqref{eq:full-q-h} with the optimal $q^{\star}(\mathbf{u}_{0:T-1})$ would be close to the true posterior.

We illustrate the gist with a \textit{flip-mnist} problem, where each digit in the MNIST dataset is randomly flipped up-and-down or left-and-right.
For a RBF kernel, we consider two variations of HVGPs in terms of the transformation: 1) Negation. We split input dimensions into two groups and negate them separately, resulting in a $4\times 100$ model. 2) Flipping the image up-and-down or left-and-right, resulting in a $4\times 100 $ model. We compare them with SVGPs using $100, 400$ inducing points, shown in Figure~\ref{fig:flipmnist}.
This experiment highlights the importance of matching the transformation $G$ with the data distribution.

\subsection{Computational Cost}

Assume that we have a $J$-way transformation, and each way is $\tilde{T}_j$-cyclic.
After decomposition this results in $\tilde{T} = \prod_{j=1}^J \tilde{T}_j$ complex-valued kernels and subsequently, $T = \prod_{j=1}^J (\lfloor  \tilde{T}_j / 2 \rfloor + 1) \geq \tilde{T} / 2^J$ real-valued kernels.
Let $T\times m$ represent using $m$ inducing points for each $t \in \{0, \dots, T-1\}$, and assume the mini-batch size $B = \mathcal{O}(m)$.

\paragraph{Time Complexity.}
The computational cost boils down to the cost of
computing $k_t(\mathbf{Z}_t, \mathbf{Z}_t)$ and the cost of variational inference.
To compute $k_t(\mathbf{Z}_t, \mathbf{Z}_t)$, we need the kernel values $k(\mathbf{Z}_t, G^s(\mathbf{Z}_t))$ for $s=0,...,\tilde{T}-1$.
If we assume the cost of applying $G$ is $c_G$, then computing $\mathbf{K}_{\mathbf{u}, \mathbf{u}}$ requires $\mathcal{O}(Tm \times \tilde{T} c_G + Tm^2 \times \tilde{T})$ operations.
Variational inference costs $\mathcal{O}(T m^3)$ time.
Therefore, the overall complexity is
$\mathcal{O}(T m^3  + 2^J T^2 m^2 + 2^J T^2 m c_G)$.
We note $2^J \leq T$, and for a single-way transformation, the cost simplifies to $\mathcal{O}(T m^3  + T^2 m^2 + T^2 m c_G)$.
In contrast, a SVGP with $Tm$ inducing points has the time complexity $\mathcal{O}(T^3m^3)$.
Furthermore, HVGPs support straightforward parallelisms
by locating computations of $k_t$
on separate devices.

\paragraph{Space Complexity.} For computing $k_t(\mathbf{Z}_t, \mathbf{Z}_t)$, we need the kernel values $k(\mathbf{Z}_t, G^s(\mathbf{Z}_t))$ for $s=0,\dots,\tilde{T}-1$, which implies the memory cost $\mathcal{O}(T m^2 \times \tilde{T} )$.
Adding the $\mathcal{O}(T m^2)$ memory for keeping variational approximations, the overall space complexity is, $\mathcal{O}(2^J T^2m^2)$.

\section{Experiments}
We present empirical evaluations in this section. All results were obtained using NVIDIA Tesla P100 GPUs, except in Sec~\ref{subsec:reg-ben} we used NVIDIA Tesla T4.
Code is available at \url{https://github.com/ssydasheng/Harmonic-Kernel-Decomposition}.

\subsection{Earth Elevation}

We adopt GPs to fit the ETOPO1 elevation data of the earth \citep{amante2009etopo1}. ETOPO1 bedrock models the Earth's elevations from the bedrock surface underneath the ice sheets. A location is represented by the (longitude, latitude) pair, where longitude $\in [-180, 180]$ and latitude $\in [-90, 90]$. We build the dataset by choosing a location every $0.1$ degrees of longitude and latitude, resulting in $6,480,000$ data points. The dataset is randomly split for $72\%$ training,  $8\%$ validating, and $20\%$ testing. We use a three dimensional RBF kernel between the Euclidean coordinates of any two (longitude, latitude) locations. Because moving two locations eastwards by the same amount of longitudes preserves their Euclidean distance, the kernel is invariant to the $T$-cyclic longitude translation:
\begin{align}
G([\mathrm{lon}, \mathrm{lat}]^{\top}) = [\mathrm{lon} + \frac{360}{T}, \mathrm{lat}]^{\top},
\end{align}
We set the period $T=12$ and $T=24$, so that $G$ moves a point eastwards by $30$ and $15$ degrees, respectively. Then we resolve Hermitian kernels to obtain $\lfloor 12 / 2 \rfloor + 1 = 7, \lfloor 24 / 2 \rfloor + 1 = 13$ real-valued kernels following Sec~\ref{subsec:resolve-hermitian}.

\begin{table}[t]
\centering
\small
\begin{tabular}{ccccc}
\toprule
 Model        & Test RMSE   & Test NLL    & Time \\
\midrule
SKI      & \textbf{0.145} & 1.313 & 1.18h \\
\midrule
1k       & 0.252 & 0.040 & 0.38h  \\
3k       & 0.208 & -0.146 & 2.47h \\
5k       & 0.196 & -0.203 & 8.70h \\
7x1k     & 0.189 & -0.246 & 1.55h \\
13x1k    & \textbf{0.177} & \textbf{-0.314} & 4.25h  \\
\bottomrule
\end{tabular}
\caption{Test performances on Earth elevations. \label{tab:elevation}}
\end{table}

We compare SVGPs with $1k, 3k, 5k$ inducing points and the HVGPs with $7 \times 1k, 13 \times 1k$ inducing points. We parallelize HVGPs using 4 GPUs, while SVGPs use only 1 GPU since it cannot be easily parallelized. All models are optimized using the Adam optimizer with learning rate 0.01 for 100K iterations. We also compare with SKI \citep{wilson2015kernel}. SKI runs into an out-of-memory error because of the large dataset, so we train it using a random $600,000$ subset of the training data. The performances are shown in Table~\ref{fig:elevation} and the predictive means are visualized in Figure~\ref{fig:elevation}. From both the table and the figure, we observe using more inducing points in variational GP models fits the dataset substantially better. Moreover, because of the decomposed structures and the parallelisms, HVGPs use more inducing points but run faster. In comparison, SKI uses $1M$ inducing points and achieves the best RMSE, but its NLL is much worse compared to variational GPs.

\begin{figure}[t]
    \centering
    \includegraphics[width=\columnwidth]{figures/translatemnist_share_loss_acc_M50.pdf}
    \caption{Translate-MNIST.  We plot how \textit{left: train loss} and \textit{right: test acc} evolve with training. We compare $16 \times M$ HVGPs with 1) \textit{NEG}: negations; 2) \textit{TRAN}: translations; 3) \textit{TRAN-S}: translations with shared inducing points. We observe the HVGPs using translations even outperform the SVGP with \textit{$800$} inducing points while the HVGP using negations performs substantially worse. Moreover, though the HVGP with \textit{TRAN-S} has only $50$ trainable inducing points, it performs similarly with the \textit{TRAN} model. \label{fig:trans-mnist}}
    \vspace{-0.2cm}
\end{figure}

\begin{figure*}[t]
\centering
\includegraphics[width=0.98\textwidth]{figures/bar_M1K_inducing_test_nll-matern32.pdf}
\caption{Test negative log-likelihoods on regression benchmarks. We compare using $M, 2M, 8\times M, 8\times 2M$ for $M=1000$. We observe that the $8\times M, 8 \times 2M$ outperform the standard $M$ and $2M$ inducing points for the most datasets. \label{fig:reg}}
\end{figure*}

\begin{figure*}[t]
\centering
\includegraphics[width=0.98\textwidth]{figures/nystrom_opt.pdf}
\caption
Nystr\"om approximation errors measured by $\tr(\mathbf{K}_{\mathbf{f}\bbf} - \mathbf{K}_{\mathbf{f} \mathbf{u}} \mathbf{K}_{\mathbf{u} \mathbf{u}}^{-1} \mathbf{K}_{\mathbf{u} \mathbf{f}})$.
We compare SVGPs using $M, 2M$ inducing points, and HVGPs with $8 \times M$ inducing points. For the transformation in HVGPs, we include both the negation along the standard axes and along the principal components.
We observe that negations along PCA directions usually outperform negations along standard axes, by demonstrating smaller approximation errors. And both variations of HVGPs perform better than the SVGPs. We also observe a consistency between the Nystr\"om approximation errors and the regression performances. For example, the SVGP ($2M$) performs better than the HVGP ($8\times M$) for the \textit{buzz} dataset in Figure~\ref{fig:reg}, and this is similarly reflected by the Nystr\"om approximation errors.
\label{fig:nystrom-appr}}
\end{figure*}

\subsection{Translate-MNIST}
The Elevation experiment uses one-way translations in HVGPs. In this section, we consider two-way translations for the \textit{Translate-MNIST} dataset. The dataset is obtained by translating every MNIST image leftwards and downwards by random numbers of pixels. We use an RBF kernel with shared lengthscales. The HVGP uses a 2-way transformation $G$ by translating the image leftwards or downwards by $4$ pixels. Since the MNIST images are of size $28\times 28$, $G$ is $(7, 7)$ cyclic. After resolving Hermitian kernels, we arrive at $(1 + \lfloor 7 /2 \rfloor) \times (1 + \lfloor 7 /2 \rfloor) = 16$ groups.

We compare the $16 \times 50$ translation HVGP with SVGPs using $50$ and $800$ inducing points. We further consider a variation of the HVGP by sharing the inducing points $\mathbf{Z}$ as in Theorem~\ref{lem:q-equi}. We also include a $16 \times 50$ HVGP with 4-way negations whose transformation does not match the input-space distribution. We optimize all models using the Adam optimizer with learning rate 0.001 for 100K iterations. The results are shown in Figure~\ref{fig:trans-mnist}.

\subsection{Regression Benchmarks}\label{subsec:reg-ben}

We also evaluate our method on standard regression benchmarks, whose training data sizes range from $200$ thousand to $1$ million. Following \citet{wang2019exact}, we use the Mat\'ern 3/2 kernel with shared lengthscales. We consider the $J$-way composition of negations.
Specifically, we conduct negations over PCA directions. We split the PCA directions into $J$ subsets, and applying negations over these subsets results in a $2^J \times m$ model.  A visual comparison between the negation along axes and the negation along principal directions is shown in Figure~\ref{fig:pca}.

\begin{figure}[h]
\centering
\includegraphics[width=\columnwidth]{figures/pca1.pdf}
\caption{Negation along axes (\textit{left}) and along principled directions (\textit{right}). The shaded area represents the input distribution. We observe that $G^{1,0}(\mathbf{x}), G^{0,1}(\mathbf{x})$ are out of the data distribution when transforming along axes. In comparison, when transforming along PCA directions, the whole orbit is in-distribution.\label{fig:pca}}
\end{figure}

We compare SVGPs using $M$ and $2M$ inducing points with HVGPs using $8\times M$ and $8 \times 2M$ inducing points for $M=1000$. For HVGPs, we use $3$-way negations over PCA directions, and we use $8$ GPUs to place the computations of each GP in parallel.  The results for negative log likelihoods (NLLs) are reported in Figure~\ref{fig:reg}. We also report the root mean squared error (RMSE) performances in Figure~\ref{app:fig:reg-rmse} in the Appendix.

\begin{figure}[t]
\centering
\includegraphics[width=\columnwidth]{figures/all_M1K_training-wilson_3droad-matern32.pdf}
\caption{Test negative log-likelihoods during training and the training time per iteration for the \textit{3droad} dataset.\label{fig:reg-training}}
\end{figure}

In Figure~\ref{fig:reg-training} we plot the evolution of the test negative log-likelihoods during training, and the training time per iteration, for the \textit{3droad} dataset. In Figure~\ref{fig:reg-training}, we observe that using more inducing points enables learning the dataset more quickly and converging to a better minima. Furthermore, due to the benefit of parallelism, the $8\times M$ HVGP has comparable running time compared to the standard SVGP with $M$ inducing points. And it is much faster compared to the SVGP with $2M$ inducing points in spite of the improved performance. Moreover, the computational bottleneck of SVGPs lies in the Cholesky decomposition, which does not support easy parallelism.

The performance of variational GPs relies largely on how well can the inducing points summarize the dataset, which can be measured by the accuracy of Nystr\"om approximation $\mathbf{K}_{\mathbf{f} \mathbf{f}} \approx \mathbf{K}_{\mathbf{f} \mathbf{u}} \mathbf{K}_{\mathbf{u} \mathbf{u}}^{-1} \mathbf{K}_{\mathbf{u} \mathbf{f}}$ \citep{drineas2005nystrom, titsias2009variational, burt2019rates}.
We compare the Nystr\"om approximation errors with all methods, using the trace norm defined as $\tr(\mathbf{K}_{\mathbf{f} \mathbf{f}} - \mathbf{K}_{\mathbf{f} \mathbf{u}} \mathbf{K}_{\mathbf{u} \mathbf{u}}^{-1} \mathbf{K}_{\mathbf{u} \mathbf{f}})$.
For HVGPs,
the trace norm is computed as
\begin{align}
    \sum_{t=0}^{T-1}\tr(\mathbf{K}_{t, \mathbf{f} \mathbf{f}} - \mathbf{K}_{t, \mathbf{f} \mathbf{u}} \mathbf{K}_{t, \mathbf{u} \mathbf{u}}^{-1} \mathbf{K}_{t, \mathbf{u} \mathbf{f}}).
\end{align}
where we use $\mathbf{K}_{t, \cdot}$ to represent the kernel $k_t$. For each dataset, we randomly sample $3000$ points as $\mathbf{X}$ and use a Mat\'ern 3/2 kernel whose lengthscales are set based on the median
heuristic. We initialize the inducing points $\mathbf{Z}$ using K-means and optimize them to minimizing the trace error. We compare SVGPs, HVGPs with negations along axes, and HVGPs with negations along principled directions. The results are shown in Figure~\ref{fig:nystrom-appr}.

\subsection{CIFAR-10 Classification}
In this subsection we  conduct experiments on the CIFAR10 classification problem using deep convolutional Gaussian processes \citep{blomqvist2019deep, dutordoir2019translation}, which combine the deep GP with the convolutional inducing features \citep{van2017convolutional}. Following the settings in \citet{shi2020sparse}, we compare HVGP with SVGP on both one-layer and multi-layer convolutional GPs.

For HVGPs we use the negation transformations on the inducing filters $G(\mathbf{z})=-\mathbf{z}$. We compare HVGPs using \textcolor[HTML]{3182bd}{2xM}, \textcolor[HTML]{3182bd}{4xM} inducing points with SVGPs using \textcolor[HTML]{31a354}{M}, \textcolor[HTML]{31a354}{2M} inducing points. For the HVGP (\textcolor[HTML]{3182bd}{4xM}), we use $4$ GPUs to achieve parallelism. We also compare with the 2-way decomposed model in \citet{shi2020sparse} termed as \textcolor[HTML]{d95f0e}{M+M}. The results are summarized in Table~\ref{tab:cifar}. We observe that using more inducing filters results in better performances. In particular, the HVGP (\textcolor[HTML]{3182bd}{4xM}) achieves the best NLLs. Because of the parallelism, the HVGP (\textcolor[HTML]{3182bd}{4xM}) also has comparable running time with the HVGP (\textcolor[HTML]{3182bd}{2xM}), and both are faster than \textcolor[HTML]{31a354}{2M} and \textcolor[HTML]{d95f0e}{M+M} for deep models.

\begin{table}[t] \vskip -0.05in
\resizebox{\columnwidth}{!}{%
\begin{tabular}{ccccc}
\toprule
M                          & Model        & ACC   & NLL    & sec/iter \\
\midrule
\multirow{5}{*}{384x0, 1K}      & \textcolor[HTML]{31a354}{M}   & 65.70$\pm$0.06 &  1.65$\pm$0.00 & 0.21  \\
                                & \textcolor[HTML]{31a354}{2M}   & \textbf{67.84$\pm$0.07} &  1.52$\pm$0.00 & 0.39  \\
                                & \textcolor[HTML]{d95f0e}{M+M}  & 67.67$\pm$ 0.07 & \textbf{1.50 $\pm$0.01} & 0.39 \\
                                & \textcolor[HTML]{3182bd}{2xM} & 66.26$\pm$1.11 &  1.76$\pm$0.17 & 0.45  \\
                                & \textcolor[HTML]{3182bd}{4xM}  & 67.76$\pm$0.05 &  \textbf{1.51$\pm$0.01} & 0.52  \\
\midrule
\multirow{5}{*}{384x1, 1K}      & \textcolor[HTML]{31a354}{M}   & 76.40$\pm$0.02 &  1.03$\pm$0.00 & 0.16  \\
                                &\textcolor[HTML]{31a354}{2M}   & 77.11$\pm$0.10 &  1.00$\pm$0.00 & 0.47  \\
                                & \textcolor[HTML]{d95f0e}{M+M}  & \textbf{77.48$\pm$0.10} & 0.98$\pm$0.01 & 0.41 \\
                                & \textcolor[HTML]{3182bd}{2xM}  & 77.09$\pm$0.18 &  1.00$\pm$0.00 & 0.37  \\
                                & \textcolor[HTML]{3182bd}{4xM}  & 77.30$\pm$0.17 &  \textbf{0.95$\pm$0.00} & 0.36  \\
\midrule
\multirow{5}{*}{384x2, 1K}      & \textcolor[HTML]{31a354}{M}    & 79.01$\pm$0.11 &  0.86$\pm$0.00 & 0.17  \\
                                &\textcolor[HTML]{31a354}{2M}   & 80.27$\pm$0.04 &  0.81$\pm$0.00 & 0.52  \\
                                & \textcolor[HTML]{d95f0e}{M+M}  & 79.98 $\pm$0.21 & 0.80$\pm$0.01  & 0.46 \\
                                & \textcolor[HTML]{3182bd}{2xM}   & 80.04$\pm$0.04 &  0.80$\pm$0.00 & 0.37  \\
                                & \textcolor[HTML]{3182bd}{4xM}   & \textbf{80.52$\pm$0.20} &  \textbf{0.75$\pm$0.01} & 0.37  \\
\midrule
\multirow{5}{*}{384x3, 1K}      & \textcolor[HTML]{31a354}{M}  & 82.41$\pm$0.08 &  0.73$\pm$0.01 & 0.40  \\
                                &\textcolor[HTML]{31a354}{2M}   & - & - & -  \\
                                & \textcolor[HTML]{d95f0e}{M+M}  & 83.26$\pm$0.19 & 0.69$\pm$0.01 & 1.24 \\
                                & \textcolor[HTML]{3182bd}{2xM}   & \textbf{84.97$\pm$0.08} &  0.60$\pm$0.00 & 0.90 \\
                                &  \textcolor[HTML]{3182bd}{4xM}  & \textbf{84.85$\pm$0.11} &  \textbf{0.58$\pm$0.00} & 0.90  \\
\bottomrule
\end{tabular}}
\caption{Deep Convolutional GPs for CIFAR-10 classification. Previous SOTA \citep{shi2020sparse} achieves ACC=$80.33$, NLL=0.82, and $1.25$ sec/iter.
We use $384\text{x}\ell, 1K$ to represent a $(\ell+1)$-layer model with a respective number of inducing points in each layer.
We compare \textcolor[HTML]{31a354}{M}, \textcolor[HTML]{31a354}{2M}, \textcolor[HTML]{d95f0e}{M+M}, \textcolor[HTML]{3182bd}{2xM},  \textcolor[HTML]{3182bd}{4xM}. We used $4$ GPUs for the \textcolor[HTML]{3182bd}{4xM} model to achieve parallelism. For the four-layer model, using \textcolor[HTML]{31a354}{2M} inducing points did not fit in memory. Instead we used a model with~$(700\text{x}3, 1600)$ inducing points and achieved ACC=$82.89\pm 0.05$, NLL=$0.73 \pm 0.00$, and $1.10$ sec/iter. \label{tab:cifar}}
\end{table}

\section{Kernel Fourier Transform}
We have shown that the Fourier series of length $T$ form the harmonic kernel decomposition with $T$ kernels. Intuitively, if $T \to \infty$, we obtain a “continuous” frequency representation of the kernel, which would be akin to a Fourier transform.

Consider the transformation $G^\mathbf{s}: \mathcal{X} \to \mathcal{X}, \mathbf{s} \in \mathbb{R}^J$, corresponding to $J$-way transformations. We assume the transformation $G$ is $\mathbf{1}$-periodic: $G^\mathbf{0}(\mathbf{x})=G^\mathbf{1}(\mathbf{x})=\mathbf{x}, G^{\mathbf{s}_1+\mathbf{s}_2}(\mathbf{x})=G^{\mathbf{s}_1}(G^{\mathbf{s}_2}(\mathbf{x})), \forall \mathbf{s}_1, \mathbf{s}_2 \in \mathbb{R}^J$. A kernel is $G$-invariant if for any $\mathbf{s} \in \mathbb{R}^J$, $k(G^\mathbf{s}(\mathbf{x}), G^\mathbf{s}(\mathbf{x}')) = k(\mathbf{x}, \mathbf{x}')$.

Given the inputs $\mathbf{x}, \mathbf{x}'$, we consider the space of kernel values: $k(\mathbf{x}, G^{\mathbf{s}}(\mathbf{x}')), \mathbf{s} \in \mathbb{R}^J$. For $\mathbf{t} \in \mathbb{R}^J$, we define the complex-valued function $k_{\mathbf{t}}: \mathcal{X} \times \mathcal{X} \to \mathbb{C}$ using the Fourier transform,
\begin{align}
    k_{\mathbf{t}}(\mathbf{x}, \mathbf{x}') = \int_{\mathbb{R}^J} e^{-2\pi i \mathbf{s}^{\top}\mathbf{t}} k(\mathbf{x}, G^\mathbf{s}(\mathbf{x}')) \mathrm{d}\mathbf{s},
\end{align}
In this way, $k_\mathbf{t}(\mathbf{x}, \mathbf{x}')$ captures the frequency of $\mathbf{t}$  in the function $\mathbf{s} \to k(\mathbf{x}, G^\mathbf{s}(\mathbf{x}'))$. Similar to the harmonic kernel decomposition, we show an alternative representation of the kernel using $k_t$.
\begin{theorem}[Harmonic Kernel Representation]\label{thm:kern-ft}
\begin{align}
    k(\mathbf{x}, \mathbf{x}') = \int_{\mathbb{R}^J} k_\mathbf{t}(\mathbf{x}, \mathbf{x}') \mathrm{d} \mathbf{t}.
\end{align}
Moreover, $k_\mathbf{t}$ is a kernel for all $\mathbf{t} \in \mathbb{R}^J$.
\end{theorem}

\begin{proof}[\textbf{Proof of Theorem~\ref{thm:kern-ft}}]
We prove this theorem by the following derivation,
\begin{align}
    \int_{\mathbb{R}^J} k_{\mathbf{t}}(\mathbf{x}, \mathbf{x}') \mathrm{d}\mathbf{t}     &=  \int_{\mathbb{R}^J} \int_{\mathbb{R}^J} e^{-2\pi i \mathbf{s}^{\top}\mathbf{t}} k(\mathbf{x}, G^\mathbf{s}(\mathbf{x}')) \mathrm{d}\mathbf{s} \mathrm{d}\mathbf{t}
    = \int_{\mathbb{R}^J}  k(\mathbf{x}, G^\mathbf{s}(\mathbf{x}'))  \int_{\mathbb{R}^J}e^{-2\pi i \mathbf{s}^{\top}\mathbf{t}} \mathrm{d}\mathbf{t} \mathrm{d}\mathbf{s} \notag \\
    &= \int_{\mathbb{R}^J}  k(\mathbf{x}, G^\mathbf{s}(\mathbf{x}'))  \delta_{\mathbf{s}} \mathrm{d}\mathbf{s}
    = k(\mathbf{x}, \mathbf{x}'). \notag
\end{align}
where we used the property that the Fourier transform of the constant function is the delta function.

To show that $k_t$ is a kernel, we prove the following equality,
\begin{align}
    \int_{\mathbb{R}^J}  \int_{\mathbb{R}^J} e^{-2\pi i \mathbf{t}^{\top}(\mathbf{s}_2 - \mathbf{s}_1)} k(G^{\mathbf{s}_1}(\mathbf{x}), G^{\mathbf{s}_2}(\mathbf{x}')) \mathrm{d}\mathbf{s}_1 \mathrm{d}\mathbf{s}_2
    &=\int_{\mathbb{R}^J}  \int_{\mathbb{R}^J} e^{-2\pi i \mathbf{t}^{\top}\mathbf{s}_2} k(\mathbf{x}, G^{\mathbf{s}_2}(\mathbf{x}')) \mathrm{d}\mathbf{s}_1 \mathrm{d}\mathbf{s}_2 \notag \\
    &= \int_{\mathbb{R}^J} e^{-2\pi i \mathbf{t}^{\top}\mathbf{s}_2} k(\mathbf{x}, G^{\mathbf{s}_2}(\mathbf{x}')) \mathrm{d}\mathbf{s}_2
    = k_\mathbf{t}(\mathbf{x}, \mathbf{x}'). \notag
\end{align}
\end{proof}

We demonstrate the Kernel Fourier Transform by considering a stationary kernel on the unit circle. We denote the input $x$ as the angle, then the kernel admits the form $k(x, x') = \kappa(x - x')$, where $\kappa$ is a periodic function of period $2\pi$. Let $\kappa_0(t) = \kappa(t) \mathbb{I}[0 \leq t < 2\pi]$, then
\begin{align}
    k(x, x') = \kappa(x -x') = \sum_{n \in \mathbb{Z}} \kappa_0(x - x' - 2 \pi n), \notag
\end{align}
Let $G^{s}(x) = x + 2\pi s$, we obtain,
\begin{align}
    k_t(x, x')
    &= \sum_{n \in \mathbb{Z}}  \int_{\mathbb{R}} e^{-2\pi i s t} \kappa_0(x - x' - 2\pi s - 2 \pi n) \mathrm{d}s
    = \frac{1}{2\pi}\sum_{n \in \mathbb{Z}}  \int_{\mathbb{R}} e^{-i (x-x'-2\pi n - w) t} \kappa_0(w) \mathrm{d}w \notag \\
    &= e^{-it(x - x')}  \left[\frac{1}{2\pi} \int_{\mathbb{R}} e^{i w t} \kappa_0(w)  \mathrm{d}w \right] \sum_{n \in \mathbb{Z}} e^{2\pi i nt }
    = e^{-it(x - x')} \hat{\kappa}_0(t) \sum_{n \in \mathbb{Z}} \delta(t - n). \notag
\end{align}
where $\hat{\kappa}_0$ is the inverse Fourier transform of $\kappa_0$. Then we have the Fourier series,
\begin{align}
    k(x, x') = \int_\mathbb{R} k_t(x, x') dt = \sum_{n \in \mathbb{Z}} \hat{\kappa}_0(n) e^{-i n(x - x')}. \notag
\end{align}

\section{Inter-domain Inducing Points Formulation}

We present an inter-domain inducing points interpretation of the harmonic kernel decomposition. An inter-domain inducing point is a function $w: \mathcal{X} \to \mathbb{C}$ whose inducing variable is defined as,
\begin{align}
    u_w = \int f(\mathbf{x}) w(\mathbf{x}) d\mathbf{x},
\end{align}
We introduce $T$ kinds of inter-domain inducing points. For $t=0, ..., T-1$, given $\mathbf{z}_t \in \mathcal{X}$, the inter-domain inducing point of the $t$-th kind is,
\begin{align}
    w_t &= \sum_{s=0}^{T-1} \mathbf{F}^H_{t, s} \delta_{G^s(\mathbf{z}_t)}, \\
    u_{w_t} &= \int f(\mathbf{x}) w_t(\mathbf{x}) d\mathbf{x} = \sum_{s=0}^{T-1} \mathbf{F}^H_{t, s} f(G^s(\mathbf{z}_t)),
\end{align}
Therefore, we can generalize the kernel function to include inter-domain inputs,
\begin{align}
    k(\mathbf{x}, w_t) &= \mathbb{E}[f(\mathbf{x}) u^H_{w_t}] = \sum_{s=0}^{T-1} \mathbf{F}_{t, s} k(\mathbf{x}, G^s(\mathbf{z}_t)), \\
    k(w_t, w'_t) &= \mathbb{E}[u_{w_t}u^H_{w'_t} ] = \sum_{s=0}^{T-1}\sum_{s'=0}^{T-1} \mathbf{F}^H_{t, s} \mathbf{F}_{t, s'} k(G^s(\mathbf{z}_t), G^{s'}(\mathbf{z}'_t))= \sum_{s=0}^{T-1} \mathbf{F}_{t, s} k(\mathbf{z}_t, G^s(\mathbf{z}'_t)),
\end{align}
where the last equality is based on Lemma~\ref{app:lem:kfs-1}. Furthermore, for $0 \leq t \neq t' \neq T-1$,
\begin{align}
    k(w_{t}, w_{t'}) &= \sum_{s=0}^{T-1}\sum_{s'=0}^{T-1} \mathbf{F}_{t, s} \mathbf{F}_{t', s'}^H k(G^s(\mathbf{z}_t), G^{s'}(\mathbf{z}'_{t'})) = 0. \notag
\end{align}
where the last equality is based on Lemma~\ref{app:lem:kfs-orth}. We observe,
\begin{align}
    k(\mathbf{x}, w_t) = k_t(\mathbf{x}, \mathbf{z}_t); \; k(w_t, w'_t) = k_t(\mathbf{z}_t, \mathbf{z}'_t).
\end{align}
Now we find that \emph{the proposed inter-domain inducing points formulation is equivalent to the kernel Fourier series.} Furthermore, the equivalence also reinterprets HVGPs as standard SVGPs using inter-domain inducing points while enforcing block diagonal posterior covariances.

\section{More Experiments and Details}

\subsection{Toy Visualization}

HVGPs are based on the decomposed GP formulation and assume independent variational posteriors. Therefore, the predictions on a target location $\mathbf{x}^{\star}$ can be decomposed as the combination of independent elements,
\begin{align}
    &\mathcal{N}(\mathbf{0}, \mathbf{K}_{\star \star} - \sum_{t=0}^{T-1} \mathbf{K}_{t, \star \mathbf{u}} \mathbf{K}_{t, \mathbf{u} \mathbf{u}}^{-1} \mathbf{K}_{t, \mathbf{u} \star})  ) + \sum_{t=0}^{T-1} \mathcal{N}(\mathbf{K}_{t, \star \mathbf{u}_t} \mathbf{K}_{t, \mathbf{u}_t \mathbf{u}_t}^{-1} \bm{\mu}_t, \mathbf{K}_{t, \star \mathbf{u}_t} \mathbf{K}_{t, \mathbf{u}_t \mathbf{u}_t}^{-1}  \mathbf{S}_t \mathbf{K}_{t, \mathbf{u}_t \mathbf{u}_t}^{-1} \mathbf{K}_{t, \star \mathbf{u}_t} ) \notag
\end{align}
where we use $\mathbf{K}_{t, \cdot}$ to represent the kernel $k_t$. The first term in the prediction represents the error of the Nystr\"om approximation, and the remaining terms contain the predictions from all subprocesses.

In this section we conduct a Snelson's 1D toy experiment to visualize the posterior predictions and each term. We set $T=2, G(x)=-x, m=5$, which results in HVGP ($2\times 5$). Because the original training inputs are positive, we preprocess it by subtracting the inputs by the mean.
The results are shown in Figure~\ref{fig:toy}. We find that using $2 \times 5$ inducing points fit the training data well, and generate reasonable predictive uncertainty as well. The predictions for the two GPs correspond to the symmetric and the antisymmetric fraction, respectively.

\begin{figure}[t]
\centering
\includegraphics[width=0.9\textwidth]{figures/toy.png}
\caption{Posterior predictions on the Snelson dataset, where shaded bands correspond to intervals of $\pm 3$  standard deviation. The leftmost \textit{posterior prediction} plot is the combination of the right three plots, plus the observation variance. We also visualize the associated inducing points for each plot. We observe that the HVGP predictions are separated as the symmetric fraction and the antisymmetric fraction.\label{fig:toy}}
\end{figure}

\subsection{Regression Benchmarks}

We use the Mat\'ern 3/2 kernel with shared lengthscales across input dimensions.
For HVGPs, the transformation $G$ is by negating over PCA directions. We split the PCA directions into $J$ subsets, then apply negations over which results in a $2^J \times M$ model. We let the $j$-th subset contain the directions with the $j \text{th large}, (J+j) \text{th large}, ...$ eigenvalues, so that the principal subspace is covered well. Except for the \textit{year} dataset which has a standard train/test split, each dataset is randomly split into $64\%$ training, $16\%$ validating, $20\%$ testing sets and is averaged over $3$ random splits. We initialize the inducing points using K-means and initialize the kernel lengthscale using the median heuristic. The Gaussian likelihood variance is initialized at $0.1$. For all experiments, we optimize for 30k iterations with the Adam optimizer using learning rate 0.003 and batch size 256. We visualize the results for test RMSEs in Figure~\ref{app:fig:reg-rmse}, and how each criterion evolves along training in Figure~\ref{app:fig:reg-training-all}.

\begin{figure}[h]
\centering
\includegraphics[width=0.98\textwidth]{figures/bar_M1K_inducing_test_rmse_sq-matern32.pdf}
\caption{Test RMSEs on regression benchmarks. We compare SVGPs using $M, 2M$ inducing points and HVGPs using $8\times M, 8\times 2M$ inducing points, for $M=1000$. \label{app:fig:reg-rmse}}
\end{figure}

\begin{figure*}[t]
\centering
\includegraphics[width=0.98\textwidth]{figures/all_M1K_training-matern32-v2.pdf}
\caption{How \textit{train loss}, \textit{test rmse}, and \textit{test nll} evolve during training. We compare SVGPs using $M, 2M$ inducing points and HVGPs using $8\times M, 8\times 2M$ inducing points, for $M=1000$.\label{app:fig:reg-training-all}}
\end{figure*}

\subsection{CIFAR-10 Classification}

\begin{minipage}{\textwidth}
  \hspace{2em}
  \begin{minipage}[b]{0.2\textwidth}
    \centering
\includegraphics[width=\textwidth]{figures/convgp-G.pdf}\captionof{figure}{Group Splitting for the HVGP ($4\times M$). \label{app:fig-conv-g}}
  \end{minipage}
  \hspace{-2em}
  \begin{minipage}[b]{0.8\textwidth}
    \centering
    \begin{tabular}{ccccc}
\toprule
             & 1-layer   & 2-layer   & 3-layer   & 4-layer   \\
\midrule
filter size  & 5         & 5, 4      & 5,4,5     & 5,4,5,4   \\
stride size  & 1         & 1,2       & 1,2,1     & 1,1,1,1   \\
channel num  & -         & 10        & 10,10     & 16,16,16  \\
pooling      & -         & -         & -         & mean      \\
pooling size & -         & -         & -         & 1,2,1     \\
padding      & SAME      & SAME      & SAME      & SAME      \\
M            & 384x0, 1K & 384x1, 1K & 384x2, 1K & 384x3, 1K \\
\bottomrule
\end{tabular}
      \captionof{table}{Model Configurations for Deep Convolutional Gaussian processes. \label{app:tab:dcgp-conf}}
    \end{minipage}
  \end{minipage}

For deep Gaussian processes, we let $M_l, h_{l}$ be the number of inducing points and the number of input units in the $l$-th layer, respectively. The variational posterior for the inducing points $\mathbf{U}^l, \mathbf{U}^l \in \mathbb{R}^{M_l \times h_{l+1}}$ in the $l$-th layer is usually a multivariate Gaussian,
\begin{align}
    q(\mathrm{vec}(\mathbf{U}^l)) = \mathcal{N}(\mathrm{vec}(\mathbf{M}^l), \bm{\Sigma}^l),
\end{align}
where $\mathbf{M}^l \in \mathbb{R}^{M_l \times h_{l+1}}, \bm{\Sigma}^l \in \mathbb{R}^{(M_l  h_{l+1})\times (M_l  h_{l+1})}$  are the mean and the covariance, respectively. A commonly-used structure for $\bm{\Sigma}^l$ is the block-diagonal covariance \citep{salimbeni2017doubly} , i.e., assuming independence between output channels. However, the true posterior is not independent. Moreover, such covariance involves $h_{l+1}$  covariances of shape $M_l \times M_l$, which are both memory intensive and computation intensive. Therefore, following \citet{park2018deep}, we use the Kronecker-factored structure for the covariance, i.e., $\bm{\Sigma}^l = \bm{\Sigma}_o^{l} \otimes \bm{\Sigma}_i^{l}$, where $\bm{\Sigma}_o \in \mathbb{R}^{h_{l+1} \times h_{l+1}}, \bm{\Sigma}_i \in \mathbb{R}^{M_l \times M_l}$ correspond to the output covariance and the input covariance, respectively.

Following \citet{shi2020sparse}, all models were optimized using 270k iterations with the Adam optimizer using a learning rate 0.003 and a batch size 64. We anneal the learning rate by 0.25 every 50k iterations to ensure convergence. Unlike \citet{shi2020sparse} which used a zero mean function, we used a convolution mean function whose filter is 1 for the center pixel and 0 everywhere else, since we observe it with a better performance. We used the robust multi-class classification likelihood. For lower layers in the deep convolutional GP, we used multi-output GPs for each input patch \citep{blomqvist2019deep}; for the output layer, we used the TICK kernel \citep{dutordoir2019translation}.  The patch kernels are RBF kernels with shared lengthscales, whose lengthscales and variances are initialized at $5$. The TICK location kernel is a Mat\'ern 3/2 kernel whose lengthscales and variances are initialized at $1$ and $3$, respectively. To initialize the inducing filters, we use K-means samples from $\min(100*M, 10000)$ random input patches, while the inputs in all layers are obtained by forwarding the image through a random Xavier convnet.

For the HVGP ($2\times M$) we use the negation transformation on the inducing points $G(\mathbf{z})=-\mathbf{z}$. For the HVGP ($4 \times M$), we also use the negation transformation over two groups that are determined by pixel locations, as shown in Figure~\ref{app:fig-conv-g}.

\section{Proofs}

\subsection{Lemmas}

\begin{lemma}\label{app:lem:kfs-1} For any $t=0,...,T-1$, $\mathbf{x}, \mathbf{x}' \in \mathcal{X}$,
\begin{align}
    \sum_{s=0}^{T-1}\sum_{s'=0}^{T-1} \mathbf{F}^{H}_{t, s} \mathbf{F}_{t, s'} k(G^s(\mathbf{x}), G^{s'}(\mathbf{x}'))
    = \sum_{s=0}^{T-1} \mathbf{F}_{t, s} k(\mathbf{x}, G^s(\mathbf{x}')). \notag
\end{align}
\end{lemma}
\begin{proof}[Proof] We prove the equality by the expression of $\mathbf{F}$,
\begin{align}
    &\frac{1}{T^2}\sum_{s=0}^{T-1}\sum_{s'=0}^{T-1} e^{-i\frac{2\pi t}{T}(s'-s) }k(G^{s}(\mathbf{x}), G^{s + (s'-s)}(\mathbf{x}')) =\frac{1}{T}\sum_{s_0=0}^{T-1} e^{-i\frac{2\pi t }{T}s_0 }k(\mathbf{x}, G^{s_0}(\mathbf{x}')), \notag
\end{align}
where we used the kernel invariance to $G$. Also, since $G$ is $T$-cyclic, we changed the variable $s'-s$ to $s_0$ and $s_0$ still ranges from $0$ to $T-1$.
\end{proof}

\begin{lemma}\label{app:lem:kfs-orth}For any $0 \leq t_1 \neq t_2 \leq T-1$, $\mathbf{x}, \mathbf{x}' \in \mathcal{X}$,
\begin{align}
    \sum_{s_1=0}^{T-1} \sum_{s_2=0}^{T-1} \mathbf{F}^H_{t_1,s_1} \mathbf{F}_{t_2,s_2} k(G^{s_1}(\mathbf{x}), G^{s_2}(\mathbf{x}'))  = 0, \label{eq:kfs-orth-1} \\
     \sum_{s_1=0}^{T-1} \sum_{s_2=0}^{T-1} \mathbf{F}_{t_1,s_1} \mathbf{F}^H_{t_2,s_2} k(G^{s_1}(\mathbf{x}), G^{s_2}(\mathbf{x}'))  = 0, \label{eq:kfs-orth-2}
\end{align}
\end{lemma}
\begin{proof}[Proof]
Below we prove \eqref{eq:kfs-orth-1}.
The proof of \eqref{eq:kfs-orth-2} follows similarly.
\begin{align}
& \sum_{s_1=0}^{T-1} \sum_{s_2=0}^{T-1} \mathbf{F}^H_{t_1,s_1} \mathbf{F}_{t_2,s_2} k(G^{s_1}(\mathbf{x}), G^{s_2}(\mathbf{x}')) = \sum_{s_1=0}^{T-1} \sum_{s_2=0}^{T-1}  \mathbf{F}^H_{t_1,s_1} \mathbf{F}_{t_2,s_1 + s_2} k(G^{s_1}(\mathbf{x}), G^{s_1+s_2}(\mathbf{x}')) \notag \\
&= \sum_{s_2=0}^{T-1} k( \mathbf{x}, G^{s_2}(\mathbf{x}')) \sum_{s_1=0}^{T-1}\mathbf{F}^H_{t_1,s_1} \mathbf{F}_{t_2,s_1 + s_2}
= \sum_{s_2=0}^{T-1} k( \mathbf{x}, G^{s_2}(\mathbf{x}')) e^{-i \frac{2\pi t s_2}{T}}\sum_{s_1=0}^{T-1}\mathbf{F}^H_{t_1,s_1} \mathbf{F}_{t_2,s_1}
= 0. \notag
\end{align}
In the last step, $\sum_{s_1=0}^{T-1}\mathbf{F}^H_{t_1,s_1} \mathbf{F}_{t_2,s_1} = 0$ whenever $t_1 \neq t_2$.
This is because the columns of $\mathbf{F}$ form an orthogonal basis over the set of $T$-dimensional complex vectors.
\end{proof}

\begin{lemma}\label{app:lem:gp-decomp-value} Under the harmonic formulation, for $\mathbf{x} \in \mathcal{X}$,
\begin{align}
    f_t(\mathbf{x}) = \sum_{s=0}^{T-1} \mathbf{F}^H_{t,s} f(G^s(\mathbf{x})).
\end{align}
\end{lemma}
\begin{proof}[Proof] We consider a marginal distribution on a subset of function values,
\begin{align}
    p(\{f(G^s(\mathbf{x}))\}_{s=0}^{T-1}, \{f_0(G^s(\mathbf{x}))\}_{s=0}^{T-1}, ..., \{f_{T-1}(G^s(\mathbf{x}))\}_{s=0}^{T-1}), \notag
\end{align}
The distribution can be represented as,
\begin{align}
    f(G^s(\mathbf{x})) &= \sum_{t=0}^{T-1} f_t(G^s(\mathbf{x})) , s=0,...,T-1,\notag \\
    \mathbf{f}_t(G^{0:T-1}(\mathbf{x})) &\sim \mathcal{N}(0, \mathbf{K}_t(G^{0:T-1}(\mathbf{x}), G^{0:T-1}(\mathbf{x}))), \notag
\end{align}
We first investigate the structure of the kernel matrix $\mathbf{K}_t$,
\begin{align}
    k_t(G^j(\mathbf{x}), G^{j'}(\mathbf{x}))
    &= \sum_{s=0}^{T-1} \mathbf{F}_{t,s} k(G^j(\mathbf{x}), G^{s+j'}(\mathbf{x}))
    =\sum_{s=0}^{T-1} \mathbf{F}_{t,s} k(\mathbf{x}, G^{s+j'-j }(\mathbf{x})) = \sum_{s=0}^{T-1} \mathbf{F}_{t,s+j-j'} k(\mathbf{x}, G^{s}(\mathbf{x})) \notag \\
    &= e^{-\frac{2\pi i t(j-j')}{T}}\sum_{s=0}^{T-1} \mathbf{F}_{t,s} k(\mathbf{x}, G^{s}(\mathbf{x})) = e^{-\frac{2\pi it(j-j')}{T}} k_t(\mathbf{x}, \mathbf{x}),\notag
\end{align}
Therefore, the matrix $\mathbf{K}_t = [k_t(\mathbf{x}, \mathbf{x}) e^{-\frac{2\pi it(j-j')}{T}}]_{j,j'=0}^{T-1}$. Let $\epsilon_t \in \mathbb{R}$ be a random Gaussian noise, then the random vector of $f_t$ can be written as,
\begin{align}
    \mathbf{f}_t(G^{0:T-1}(\mathbf{x})) = [\sqrt{k_t(\mathbf{x}, \mathbf{x})} e^{-\frac{2\pi it j}{T}}  \epsilon_t]_{j=0}^{T-1},
\end{align}

Now we can compute the RHS in the lemma,
\begin{align}
    \sum_{s=0}^{T-1} \mathbf{F}^H_{t,s} f(G^s(\mathbf{x}))
    &= \sum_{s=0}^{T-1} \mathbf{F}^H_{t,s} \sum_{t'=0}^{T-1} f_{t'}(G^s(\mathbf{x})) = \sum_{t'=0}^{T-1}  \sum_{s=0}^{T-1} \mathbf{F}^H_{t,s} f_{t'}(G^s(\mathbf{x})) , \notag
\end{align}
If $t'=t$,
\begin{align}
    \sum_{s=0}^{T-1} \mathbf{F}^H_{t,s} f_{t}(G^s(\mathbf{x}))
    &=  \sum_{s=0}^{T-1} \mathbf{F}^H_{t,s} \sqrt{k_t(\mathbf{x}, \mathbf{x})} e^{-\frac{2\pi it s}{T}}  \epsilon_t
    = \sqrt{k_t(\mathbf{x}, \mathbf{x})}\epsilon_t = f_t(\mathbf{x}).
\end{align}
If $t' \neq t$,
\begin{align}
    \sum_{s=0}^{T-1} \mathbf{F}^H_{t,s} f_{t'}(G^s(\mathbf{x}))
    &=  \sum_{s=0}^{T-1} \mathbf{F}^H_{t,s} \sqrt{k_{t'}(\mathbf{x}, \mathbf{x})} e^{-\frac{2\pi it' s}{T}}  \epsilon_t = \sqrt{k_{t'}(\mathbf{x}, \mathbf{x})} \epsilon_t \sum_{s=0}^{T-1} e^{\frac{2\pi i(t-t')}{T}} = 0.
\end{align}
Therefore,
\begin{align}
    \sum_{s=0}^{T-1} \mathbf{F}^H_{t,s} f(G^s(\mathbf{x})) =f_t(\mathbf{x}).
\end{align}
Because this holds for all marginal distributions, it holds as well for the function samples from Gaussian processes.\end{proof}

\subsection{Proofs for Sec~\ref{sec:kern-decomp}}\label{app:proof:ortho}
\begin{proof}[\textbf{Proof of Proposition~\ref{prop:invariance}}]
The equality can be directly proven,
\begin{align}
    k_t(\mathbf{x}, G(\mathbf{x}'))
    &= \sum_{s=0}^{T-1} \mathbf{F}_{t, s} k(\mathbf{x}, G^{s+1}(\mathbf{x}')) = \sum_{s=0}^{T-1} \mathbf{F}_{t, s-1} k(\mathbf{x}, G^{s}(\mathbf{x}')) \notag \\
    &= \sum_{s=0}^{T-1} \frac{1}{T} e^{-i\frac{2\pi t (s-1)}{T} } k(\mathbf{x}, G^{s}(\mathbf{x}')) = e^{i\frac{2\pi t}{T}} k_t(\mathbf{x}, \mathbf{x}').
\end{align}
\end{proof}

\begin{proof}[\textbf{Proof of Theorem~\ref{thm:kern-decomp}}]\label{proof:app:kern-decomp}
The equality can be directly proven,
\begin{align}
   \sum_{t=0}^{T-1} k_t(\mathbf{x}, \mathbf{x}')
    &= \sum_{t=0}^{T-1} \sum_{s=0}^{T-1} \frac{1}{T}e^{-i\frac{2\pi st}{T}} k(\mathbf{x}, G^s(\mathbf{x}'))
    =\frac{1}{T}  \sum_{s=0}^{T-1} k(\mathbf{x}, G^s(\mathbf{x}')) \sum_{t=0}^{T-1} e^{-i\frac{2\pi st}{T} }
    = \frac{1}{T}  \sum_{s=0}^{T-1} k(\mathbf{x}, G^s(\mathbf{x}')) T\delta_{s} = k(\mathbf{x}, \mathbf{x}'). \notag
\end{align}
To prove that $k_t$ is a kernel, we observe from Lemma~\ref{app:lem:kfs-1} that $k_t(\mathbf{x}, \mathbf{x}') = \mathbf{F}_{t,:}^H \mathbf{K} \mathbf{F}_{t,:}$, where $\mathbf{K} = [k(G^{s_1}(\mathbf{x}), G^{s_2}(\mathbf{x}'))]_{s_1, s_2=0}^{T-1}$. Since $k$ is a kernel,  we conclude that $k_t(\mathbf{x}, \mathbf{x}') = \mathbf{F}_{t,:}^H \mathbf{K} \mathbf{F}_{t,:}$ is a kernel as well.

\end{proof}

\begin{proof}[\textbf{Proof of Lemma~\ref{lem:ortho}}]\label{proof:kern-ortho}
We firstly prove that, for all $t_1 \neq t_2$ and $\mathbf{x}, \mathbf{x}'$, $\innerprod{k_{t_1}(\cdot, \mathbf{x})}{k_{t_2}(\cdot, \mathbf{x}')}_{\mathcal{H}_k} = 0$. The RKHS inner product can be computed as,
\begin{align}
&\innerprod{k_{t_1}(\cdot, \mathbf{x})}{k_{t_2}(\cdot, \mathbf{x}')}_{\mathcal{H}_k} = \sum_{s_1=0}^{T-1} \sum_{s_2=0}^{T-1} \mathbf{F}_{t_1,s_1} \mathbf{F}^H_{t_2,s_2} k(G^{s_1}(\mathbf{x}), G^{s_2}(\mathbf{x}')) =0,
\end{align}
where the last equality is due to Lemma~\ref{app:lem:kfs-orth}.

Moreover, if the functions $f, g$ can be written as linear combinations of the corresponding kernels,
\begin{align}
 f(\mathbf{x}) = \sum_s a_s k_{t_1}(\mathbf{x}, \mathbf{x}^s_{t_1}) ; \;
 g(\mathbf{x}) = \sum_s b_s k_{t_2}(\mathbf{x}, \mathbf{x}^s_{t_2}) \notag
\end{align}
Following that $\innerprod{k_{t_1}(\cdot, \mathbf{x})}{k_{t_2}(\cdot, \mathbf{x}')}_{\mathcal{H}_k} $ for all $\mathbf{x}, \mathbf{x}'$, $\innerprod{f}{g}_{\mathcal{H}_k} = 0$ as well.

Based on Moore-Aronszajn Theorem \citep{aronszajn1950theory, berlinet2011reproducing}, the RKHS spaces $\mathcal{H}_{k_{t_1}}$ and $\mathcal{H}_{k_{t_1}}$ are the set of functions which are pointwise limits of Cauchy sequences in the form $f_n(\mathbf{x}) = \sum_s a_s k_{t_1}(\mathbf{x}, \mathbf{x}^s_{t_1}) $ and $g_n(\mathbf{x}) = \sum_s b_s k_{t_2}(\mathbf{x}, \mathbf{x}^s_{t_2})$, respectively. Moreover, based on the \citet[Lemma 5]{berlinet2011reproducing}, the inner product of two pointwisely convergent Cauchy sequences also converges. We conclude that for any $f \in \mathcal{H}_{k_{t_1}}, g \in \mathcal{H}_{k_{t_2}}$, $\innerprod{f}{g}_{\mathcal{H}_k} = 0$.
\end{proof}

\begin{proof}[\textbf{Proof of Proposition~\ref{prop:disjoint}}]
Without loss of generality, we only need to prove that,
\begin{align}
\mathcal{H}_1 \cap \mathcal{H}_2 = \{0\}, \notag
\end{align}
Firstly, $0 \in \mathcal{H}_1, 0\in \mathcal{H}_2$ because $\mathcal{H}_1$ and $\mathcal{H}_2$ are Hilbert spaces. Then we assume another function $f \neq 0$ and $f \in \mathcal{H}_1 \cap \mathcal{H}_2$. By Lemma~\ref{lem:ortho}, $\|f\|_{\mathcal{H}_k} = \innerprod{f}{f}_{\mathcal{H}_k} = 0$, which is contradictory to $f \neq 0$ and $\mathcal{H}_k$ being a Hilbert space.
\end{proof}

\begin{proof}[\textbf{Proof of Theorem~\ref{thm:rkhs-decom}}]
We use $\mathcal{H}_t$ to represent the RKHS corresponding to the kernel $k_t$. Given a function $f \in \mathcal{H}_k$, we firstly assume $f$ can be written as a linear combination of the kernel functions,
\begin{align}
f(\mathbf{x}) = \sum_s a_s k(\mathbf{x}, \mathbf{x}^s) , \notag
\end{align}
Based on the kernel sum decomposition, we can rewrite $f$,
\begin{align}
    f(\mathbf{x}) = \sum_s a_s \sum_{t=0}^{T-1} k_t(\mathbf{x}, \mathbf{x}^s) = \sum_{t=0}^{T-1} \underbrace{\sum_s a_s k_t(\mathbf{x}, \mathbf{x}^s) }_{:=f_t(\mathbf{x})}, \notag
\end{align}
Because $f_t$ is a linear combination of $k_t$, $f_t \in \mathcal{H}_{t}$, for $t=0,...,T-1$. Proposition~\ref{prop:disjoint} states that the RKHSs $\mathcal{H}_{t_1}, \mathcal{H}_{t_2}$ are disjoint except the zero function, thus $f = \sum_{t=0}^{T-1} f_t$ is a unique expansion of $f$ to these RKHSs. Moreover, we can represent the function $f$ alternatively,
\begin{align}
    f(\mathbf{x}) &= \innerprod{f}{k(\mathbf{x}, \cdot)}_{\mathcal{H}_k} = \sum_{t=0}^{T-1}\innerprod{f}{k_t(\mathbf{x}, \cdot)}_{\mathcal{H}_k}
    = \sum_{t=0}^{T-1}\innerprod{\sum_{t'=0}^{T-1} f_{t'}}{k_t(\mathbf{x}, \cdot)}_{\mathcal{H}_k}  = \sum_{t=0}^{T-1}\innerprod{f_t}{k_t(\mathbf{x}, \cdot)}_{\mathcal{H}_k}, \notag
\end{align}
where the last equality uses the orthogonality between RKHSs. By using the orthogonality again, we also show that,
\begin{align}
    \innerprod{f_t}{k_t(\mathbf{x}, \cdot)}_{\mathcal{H}_k} &= \innerprod{f_t}{\sum_{t'=0}^{T_1}k_{t'}(\mathbf{x}, \cdot)}_{\mathcal{H}_k}
     = \innerprod{f_t}{k(\mathbf{x}, \cdot)}_{\mathcal{H}_k} = f_t(\mathbf{x}),
\end{align}
Therefore, $f(\mathbf{x}) = \sum_{t=0}^{T-1}\innerprod{f}{k_t(\mathbf{x}, \cdot)}_{\mathcal{H}_k}$ uniquely separates $f$ into these RKHSs. More generally, if $f$ is the pointwise limits of Cauchy sequences of functions in the form of linear combinations of the kernel function. Based on the \citet[Lemma 5]{berlinet2011reproducing}, the inner product of two pointwisely convergent Cauchy sequences also converges. We conclude that for any $f \in \mathcal{H}_{k}$,
\begin{align}
    f(\mathbf{x}) = \sum_{t=0}^{T-1}\innerprod{f}{k_t(\mathbf{x}, \cdot)}_{\mathcal{H}_k},
\end{align}
uniquely decomposes the function $f$ into the RKHSs $\mathcal{H}_t, t=0,...,T-1$.

Based on \citet[Theorem 5]{berlinet2011reproducing}, the squared RKHS norm  of $f$ can be written as the sum of squared RKHS norms,
\begin{align}
    \|f\|_{\mathcal{H}_k}^2 = \sum_{t=0}^{T-1} \|f_t\|_{\mathcal{H}_t}^2 .\notag
\end{align}

\end{proof}

\subsection{Proof of Sec~\ref{sec:hvgp}}\label{app:prof:lemma-dvi}

\begin{proof}[\textbf{Proof of Theorem~\ref{lem:q-equi}}]

Under the HVGP formulation, the inducing variable $u_t = f_t(\mathbf{z})$. From Lemma~\ref{app:lem:gp-decomp-value}, we have $f_t(\mathbf{z}) = \sum_{s=0}^{T-1} \mathbf{F}^H_{t, s} f(G^s(\mathbf{z}))$.

Under the inter-domain formulation, let $w_t$ be the inter-domain inducing point corresponding to $\mathbf{z}$ in $k_t$,
\begin{align}
     w_t = \sum_{s=0}^{T-1} \mathbf{F}^H_{t, s} \delta_{G^s(\mathbf{z})}, \notag
\end{align}
Then the inducing variable corresponding to $w_t$ is,
\begin{align}\label{app:eq:ut}
    u_{w_t} = \int f(\mathbf{x}) w_t(\mathbf{x}) d\mathbf{x} = \sum_{s=0}^{T-1} \mathbf{F}^H_{t, s} f(G^s(\mathbf{z})),
\end{align}
Therefore, the two inducing variables are the same,
\begin{align}
    u_{w_t} = u_t.
\end{align}

Therefore, the variational posterior under the harmonic formulation can be rewritten in an inter-domain SVGP form,
\begin{align}
    q^{inter}(f, \mathbf{U}) = p(f|\mathbf{U}; \{\mathbf{w}_t\}_{t=0}^{T-1}) q(\vectorize{\mathbf{U}}),
\end{align}
where $\mathbf{U} :=[\mathbf{u}_0, ..., \mathbf{u}_{T-1}]^{\top}$ and $p$ is the inter-domain Gaussian process.

Now we connect the inter-domain SVGP to the standard SVGP using $\{G^t(\mathbf{Z})\}_{t=0}^{T-1}$. For the standard SVGP, the inducing variables are $\mathbf{v}_t = f(G^t(\mathbf{Z}))$, and $\mathbf{V}:=[\mathbf{v}_0, ..., \mathbf{v}_{T-1}]^{\top} \in \mathbb{C}^{T \times m}$. For the inter-domain SVGP, the inducing variables are $\mathbf{u}_t$. As shown in Eq.~\eqref{app:eq:ut}, $\mathbf{u}_t = \sum_{s=0}^{T-1} \mathbf{F}^H_{t, s} \mathbf{v}_s$, then we have the equality,
\begin{align}
    \mathbf{U} = \mathbf{F}^H \mathbf{V},
\end{align}
Because of the bijective linearity,
\begin{align}
p(f|\mathbf{U}; \{\mathbf{w}_t\}_{t=0}^{T-1}) = p(f|\mathbf{V}; \{G^t(\mathbf{Z})\}_{t=0}^{T-1}),
\end{align}
Furthermore, the variational posterior for $\mathbf{V}$ is $\mathcal{N}(\vectorize{\mathbf{V}}| \vectorize{\mathbf{M}_{v}}, \mathbf{S}_v)$, which is equivalent to the variational posterior for $\mathbf{U}$,
\begin{align}
    q(\vectorize{\mathbf{U}}) = \mathcal{N}(\vectorize{\mathbf{U}}| \vectorize{\mathbf{F}^H \mathbf{M}_{v}}, (\mathbf{I} \otimes \mathbf{F}^H) \mathbf{S}_v  (\mathbf{I} \otimes \mathbf{F})). \notag
\end{align}
So the argument has been proved.

\end{proof}

\begin{lemma}[Variational Gaussian Approximations]\label{lem:vari-gau-apprx} Let $\mathcal{N}(\bm{\mu}, \mathbf{S})$ be a Gaussian variational posterior for a SVGP, then the optimal $\mathbf{S}^{\star}$ is in the form of,
\begin{align}
    \mathbf{S}^{\star} =  \mathbf{K}_{\mathbf{u} \mathbf{u}} \left( \mathbf{K}_{\mathbf{u} \mathbf{u}} + \mathbf{K}_{\mathbf{u} \mathbf{f}} \bm{\Lambda} \mathbf{K}_{\mathbf{f} \mathbf{u}} \right)^{-1} \mathbf{K}_{\mathbf{u} \mathbf{u}}.
\end{align}
where $\bm{\Lambda} = \mathrm{diag}([\lambda_n]_{n=1}^N)$ is diagonal,
\begin{align}
    \lambda_n = -2\nabla_{\sigma^2_n}\mathbb{E}_{q(f_n)} [\log p(y_n | f_n)],
\end{align}
where $\sigma^2_n$ is the predictive variance of $f_n$ under the variational posterior.
\end{lemma}
\begin{proof}[Proof]
Given the variational posterior, the predictive distribution of $f_n$ can be computed as,
\begin{align}
    \mathcal{N}(\mathbf{k}_{f \mathbf{u}} \mathbf{K}_{\mathbf{u}\bu}^{-1}\bm{\mu}, k_{ff} +  \mathbf{k}_{f \mathbf{u}}\mathbf{K}_{\mathbf{u}\bu}^{-1} (\mathbf{S} - \mathbf{K}_{\mathbf{u}\bu})\mathbf{K}_{\mathbf{u}\bu}^{-1} \mathbf{k}_{\mathbf{u} f}), \notag
\end{align}
where we denote the predictive variance as $\sigma_n^2$. The variational posterior is optimized by maximizing the ELBO, which can be computed as,
\begin{align}
    \mathcal{L} = \sum_{n=1}^N \mathbb{E}_{q(f_n)} [\log p(y_n | f_n)] - \KL{\mathcal{N}(\bm{\mu}, \mathbf{S})}{\mathcal{N}(\mathbf{0}, \mathbf{K}_{\mathbf{u}\bu})}, \notag
\end{align}
We compute the derivatives of $\mathcal{L}$ towards $\mathbf{S}$,
\begin{align}
    \nabla_{\mathbf{S}} \mathcal{L} &= -\frac{1}{2}\sum_{n=1}^N \lambda_n \nabla_{\mathbf{S}} \sigma_n^2 - \frac{1}{2}(\mathbf{K}_{\mathbf{u}\bu}^{-1} - \mathbf{S}^{-1})
    = -\frac{1}{2}\mathbf{K}_{\mathbf{u}\bu}^{-1} (\sum_{n=1}^N \lambda_n \mathbf{k}_{ \mathbf{u} f_n} \mathbf{k}_{f_n \mathbf{u}})\mathbf{K}_{\mathbf{u}\bu}^{-1} - \frac{1}{2}(\mathbf{K}_{\mathbf{u}\bu}^{-1} - \mathbf{S}^{-1}) \notag \\
    &=  -\frac{1}{2}\mathbf{K}_{\mathbf{u}\bu}^{-1} \mathbf{K}_{\mathbf{u}\mathbf{f}} \bm{\Lambda} \mathbf{K}_{\mathbf{f} \mathbf{u}} \mathbf{K}_{\mathbf{u}\bu}^{-1} - \frac{1}{2} \mathbf{K}_{\mathbf{u}\bu}^{-1} + \frac{1}{2}\mathbf{S}^{-1},
\end{align}
Let the derivative be zero, we obtain the optimal $\mathbf{S}^{\star}$,
\begin{align}
    \mathbf{S}^{\star} = \mathbf{K}_{\mathbf{u}\bu}(\mathbf{K}_{\mathbf{u}\bu} + \mathbf{K}_{\mathbf{u}\mathbf{f}} \bm{\Lambda} \mathbf{K}_{\mathbf{f} \mathbf{u}})^{-1}\mathbf{K}_{\mathbf{u}\bu}. \notag
\end{align}
\end{proof}

\begin{proof}[\textbf{Proof of Theorem~\ref{lem:decoupled}}] Based on Lemma~\ref{lem:vari-gau-apprx}, the optimal posterior covariance is
\begin{align}
    \mathbf{S}^{\star} = \mathbf{K}_{\mathbf{u} \mathbf{u}} \left( \mathbf{K}_{\mathbf{u} \mathbf{u}} +\mathbf{K}_{\mathbf{u} \mathbf{f}}\bm{\Lambda} \mathbf{K}_{\mathbf{f} \mathbf{u}} \right)^{-1} \mathbf{K}_{\mathbf{u} \mathbf{u}},
\end{align}
Given that $\mathbf{K}_{\mathbf{u} \mathbf{u}}$ is block diagonal, by the continuous mapping theorem, it remains to prove that $\mathbf{K}_{\mathbf{u} \mathbf{f}}\bm{\Lambda}  \mathbf{K}_{\mathbf{f} \mathbf{u}}$ approaches block diagonal.

Firstly we assume that Hermitian kernels are not resolved, thus $\mathbf{K}_{\mathbf{f}\mathbf{u}} = \mathbf{K}_{\mathbf{u}\mathbf{f}}^H$.
Because $\lambda_n$ only depends on $(\mathbf{x}_n, y_n)$, for the $(\mathbf{z}_t, \mathbf{z}_{t'})$ off-diagonal element in $\mathbf{K}_{\mathbf{u} \mathbf{f}}\bm{\Lambda} \mathbf{K}_{\mathbf{f} \mathbf{u}} $,
\begin{align}
    &\frac{1}{N}\sum_{n=1}^N \lambda_n k_t(\mathbf{z}_t, \mathbf{x}_n) k_{t'}^H(\mathbf{z}_{t'}, \mathbf{x}_n) \to \mathbb{E}_{p(\mathbf{x})p(y|\mathbf{x})}[ \lambda(\mathbf{x}, y) k_t(\mathbf{z}_t, \mathbf{x})k^H_{t'}(\mathbf{z}_{t'}, \mathbf{x}) ]
    = \mathbb{E}_{p(\mathbf{x})}[ \mathbb{E}_{p(y|\mathbf{x})}[\lambda(\mathbf{x}, y)] k_t(\mathbf{z}_t, \mathbf{x})k^H_{t'}(\mathbf{z}_{t'}, \mathbf{x}) ],  \notag
\end{align}
We let $\hat{\lambda}(\mathbf{x}) := \mathbb{E}_{p(y|\mathbf{x})}[\lambda(\mathbf{x}, y)]$, then the formula can be further computed as,
\begin{align}
    & \mathbb{E}_{p(\mathbf{x})}[ \hat{\lambda}(\mathbf{x}) k_t(\mathbf{z}_t, \mathbf{x})k^H_{t'}(\mathbf{z}_{t'}, \mathbf{x}) ] \notag \\
    &=\mathbb{E}_{p(\mathbf{x})}[ \hat{\lambda}(\mathbf{x}) \sum_{s=0}^{T-1} \sum_{s'=0}^{T-1} \mathbf{F}_{t,s} \mathbf{F}^H_{t',s'} k(\mathbf{x}, G^s (\mathbf{z}_t)) k^H(\mathbf{x}, G^{s'} (\mathbf{z}_{t'}))  ] \notag \\
    &= \mathbb{E}_{p(\mathbf{x})}[ \hat{\lambda}(\mathbf{x}) \sum_{s=0}^{T-1} \sum_{s'=0}^{T-1} \mathbf{F}_{t,s} \mathbf{F}^H_{t',s+s'} k(\mathbf{x}, G^s(\mathbf{z}_t)) k^H(\mathbf{x}, G^{s+s'}(\mathbf{z}_{t'}))  ] \notag \\
    &= \mathbb{E}_{p(\mathbf{x})}[\hat{\lambda}(G^s(\mathbf{x}))  \sum_{s=0}^{T-1} \sum_{s'=0}^{T-1}\mathbf{F}_{t,s} \mathbf{F}^H_{t',s+s'}k(G^{s} (\mathbf{x}), G^{s}(\mathbf{z}_t)) k^H(G^{s} (\mathbf{x}), G^{s+s'} (\mathbf{z}_{t'}))  ] \notag \\
    &=\mathbb{E}_{p(\mathbf{x})}[\hat{\lambda}(G^s(\mathbf{x})) k( \mathbf{x}, \mathbf{z}_t) \sum_{s'=0}^{T-1} k^H(\mathbf{x}, G^{s'} (\mathbf{z}_{t'}))  \sum_{s=0}^{T-1} \mathbf{F}_{t,s} \mathbf{F}^H_{t',s+s'} ] = 0, \notag
\end{align}
In the second equality we used the periodicity of $G$; In the third equality we used the assumption that $G^s(\mathbf{x})$ has the same distribution as $\mathbf{x}$; In the last equality we used the property that $ \sum_{s=0}^{T-1} \mathbf{F}_{t,s} \mathbf{F}^H_{t',s+s'}  = 0$ for all $t \neq t'$.

Furthermore, if the Hermitian kernels are resolved in HVGP, let $T$ be the period, then $\mathbf{K}_{\mathbf{u} \mathbf{f}}\bm{\Lambda} \mathbf{K}_{\mathbf{f} \mathbf{u}}$ is a matrix of $(1+\lfloor T/2 \rfloor) \times (1+\lfloor T/2 \rfloor)$. For any off-diagonal element at $(t, t')$, $[\mathbf{K}_{\mathbf{u} \mathbf{f}}\bm{\Lambda} \mathbf{K}_{\mathbf{f} \mathbf{u}}]_{t,t'}$ equals to,
\begin{align}
    \frac{1}{N}\sum_{n=1}^N \lambda_n \left(k_t(\mathbf{z}_t, \mathbf{x}_n) + k_{T-t}(\mathbf{z}_t, \mathbf{x}_n)\right)\left(k^H_{t'}(\mathbf{z}_{t'}, \mathbf{x}_n) + k^H_{T-t'}(\mathbf{z}_{t'}, \mathbf{x}_n)\right) \notag
\end{align}
Given previous results, because $t, T-t$ are both different with $t', T-t'$, the formula becomes $0$ as well, as $N \to \infty$.
\end{proof}

\section{Discussion: Connections to Group Theory}

\section{Conclusion}
We presented the harmonic kernel decomposition which exploited input-space symmetries to obtain an orthogonal kernel sum decomposition, based on which we introduced a scalable variational GP model and analyzed how well the model approximates the true posterior of the GP. We validated its superior performances in terms of scalability and accuracy through a range of empirical evaluations.

\section*{Acknowledgements}
We thank Wesley Maddox, Greg Benton, Sanyam Kapoor, Michalis Titsias, Radford Neal, and anonymous reviewers for their insightful comments and discussions on this project. We also thank the Vector Institute for providing the scientific computing resources. SS was supported by the Connaught Fellowship. RG acknowledges support from the CIFAR Canadian AI Chairs program.