<figure><img src="image_1.png" /></figure>

HAL
open science

# Random Fourier Features for Operator-Valued Kernels

Romain Brault, Florence d'Alché-Buc

► To cite this version:

Romain Brault, Florence d'Alché-Buc. Random Fourier Features for Operator-Valued Kernels. 2017. [hal-01313005v3](https://hal.science/hal-01313005v3)

HAL Id: hal-01313005

https://hal.science/hal-01313005v3

Preprint submitted on 26 Oct 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Random Fourier Features for Operator-Valued Kernels

**Brault Romain**
LTCI

Télécom ParisTech
Paris, 46 rue Barrault, France
Université Paris-Saclay

ROMAIN.BRAULT@TELECOM-PARISTECH.FR

**Florence d'Alché-Buc**
LTCI

Télécom ParisTech
Paris, 46 rue Barrault, France
Université Paris-Saclay

FLORENCE.DALCHE@TELECOM-PARISTECH.FR

**Editor:** John Doe

**Abstract**

Many problems in Machine Learning can be cast into vector-valued functions approximation. Operator-Valued Kernels *Operator-Valued Kernels* and vector-valued Reproducing Kernel Hilbert Spaces provide a theoretical and versatile framework to address that issue, extending nicely the well-known setting of scalar-valued kernels. However large scale applications are usually not affordable with these tools that require an important computational power along with a large memory capacity. In this paper, we aim at providing scalable methods that enable efficient regression with *Operator-Valued Kernels*. To achieve this goal, we extend Random Fourier Features, an approximation technique originally introduced for translation-invariant scalar-valued kernels, to translation-invariant *Operator-Valued Kernels*. We develop all the machinery in the general context of Locally Compact Abelian groups, allowing for coping with *Operator-Valued Kernels*. Eventually, the provided approximated operator-valued feature map converts the nonparametric kernel-based model into a linear model in a finite-dimensional space.

**Keywords:** Random Fourier Feature, Operator-Valued Kernel

**1. Introduction**

Learning vector-valued functions is key to numerous tasks in Machine Learning such as multi-label classification, multi-task regression, vector field learning. Among the different families of tools that currently allow to address multiple output prediction, Operator-Valued Kernel methods (Micchelli and Pontil, 2005; Carmeli et al., 2010; Kadri et al., 2010; Brouard et al., 2011; Álvarez et al., 2012) benefit from a well founded theoretical background while being flexible enough to handle a large variety of problems. In a few words, Operator-Valued Kernels extend the classic scalar-valued kernels to functions with values in *output* Hilbert spaces. As in the scalar case, Operator-Valued Kernels (OVKs) are used to build Reproducing Kernel Hilbert Spaces (RKHS) in which representer theorems apply as for ridge regression or other appropriate loss functional. In these cases, learning a model in
the RKHS boils down to learning a function of the form $f(x) = \sum_{i=1}^{N} K(x, x_i)\alpha_i$ where $x_1, \dots, x_N$ are the training input data and each $\alpha_i, i = 1, \dots, N$ is a vector of the output space $\mathcal{Y}$, and each $K(x, x_i)$ is an operator on $\mathcal{Y}$.

However, OVKs suffer from the same drawbacks as classic (scalar-valued) kernel machines: they scale poorly to large datasets because they are exceedingly demanding in terms of memory and computations.

In this paper, we propose to keep the theoretical benefits of working with OVKs while providing efficient implementations of learning algorithms. To achieve this goal, we study feature map approximations $\tilde{\phi}$ of OVKs that allow for shallow architectures, namely the product of a (nonlinear) operator-valued feature $\tilde{\phi}(x)$ and a parameter vector $\theta$ such that
$$ \tilde{f}(x) = \tilde{\phi}(x)^* \theta. $$

To approximate OVKs, we extend the well-known methodology called Random Fourier Features (RFFs) (Rahimi and Recht, 2007; Le et al., 2013; Yang et al., 2015b; Sriperumbudur and Szabo, 2015; Bach, 2015; Sutherland and Schneider, 2015; Rudi et al., 2016) so far developed to speed up scalar-valued kernel machines. The RFF approach linearizes a shift-invariant kernel model by generating explicitly an approximated feature map $\tilde{\varphi}$. RFFs has been shown to be efficient on large datasets (Rudi et al., 2016). Additionally, it has been further improved by efficient matrix computations such as (Le et al., 2013, “FastFood”) and (Felix et al., 2016, “SORF”), that are considered as the best large scale implementations of kernel methods, along with Nyström approaches proposed in Drineas and Mahoney (2005). Moreover thanks to RFFs, kernel methods have been proved to be competitive with deep architectures (Lu et al., 2014; Dai et al., 2014; Yang et al., 2015a).

## 1.1 Outline and contributions

The paper is structured as follow. In Section 2 we recall briefly how to obtain RFFs for scalar-valued kernels and list the state of the art implementation of RFFs for large scale kernel learning. Then we define properly Operator-Valued Kernels, give some important theorems and properties used throughout this paper before giving a non exhaustive list of problem tackled with OVKs.

Then we move on to our contributions. In Section 3 we propose an RFF construction from $\mathcal{Y}$-Mercer shift invariant OVK that we call Operator-valued Random Fourier Feature (ORFF). Then we study the structure of a random feature corresponding to an OVK (without having to specify the target kernel). Eventually we use the framework used to construct ORFFs to study the regularization properties of OVKs in terms of Fourier Transform.

In Section 4 we assess theoretically the quality of our ORFF: we show that the stochastic ORFF estimator converges with high probability toward the target kernel and derive convergence rates. We also give a bound on the variance of the approximated OVK constructed from the corresponding ORFF.

In Section 5 we focus on Ridge regression with OVKs. First we study the relationship between finding a minimizer in the VV-RKHS induce by a given OVK and the feature
induced by the corresponding ORFF. Then we define a gradient based algorithm to tackle
Ridge regression with ORFF, show how to obtain an efficient implementation and study its
complexity.

Eventually we end this paper by some numerical experiments in Section 6 on toy and real
datasets before giving a general conclusion in Section 7.

## 2. Background

Notations used throughout this paper are summarized in Table 1.

### 2.1 Random Fourier Feature maps

The Random Fourier Features methodology introduced by Rahimi and Recht (2007) provides
a way to scale up kernel methods when kernels are Mercer and translation-invariant. We
view the input space $\mathcal{X}$ as the group $\mathbb{R}^d$ endowed with the addition law. Extensions to other
group laws such as Li et al. (2010) are described in Subsection 3.2.2 within the general
framework of operator-valued kernels.

Denote $k: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ a positive definite kernel (Aronszajn, 1950) on $\mathbb{R}^d$. A kernel $k$ is said to be *shift-invariant* or *translation-invariant* for the addition if for all $(x, z, t) \in (\mathbb{R}^d)^3$ we have $k(x + t, z + t) = k(x, z)$. Then, we define $k_0: \mathbb{R}^d \to \mathbb{R}$ the function such that $k(x, z) = k_0(x - z)$. The function $k_0$ is called the *signature* of kernel $k$. If $k_0$ is a continuous function we call the kernel “Mercer”. Then, Bochner’s theorem (Folland, 1994) is the theoretical result that leads to the Random Fourier Features.

**Theorem 1 (Bochner's theorem)**

Any continuous positive-definite function (e. g. a Mercer kernel) is the Fourier Transform of a bounded non-negative Borel measure.

It implies that any positive-definite, continuous and shift-invariant kernel $k$, have a continuous and positive-definite signature $k_0$, which is the Fourier Transform $\mathcal{F}$ of a non-negative measure $\mu$. We therefore have the $k(x, z) = k_0(x - z) = \int_{\mathbb{R}^d} \exp(-i\langle\omega, x - z\rangle)d\mu(\omega) = \mathcal{F}[k_0](\omega)$. Moreover $\mu = \mathcal{F}^{-1}[k_0]$. Without loss of generality, we assume that $\mu$ is a probability measure, i. e. $\int_{\mathbb{R}^d} d\mu(\omega) = 1$ by renormalizing the kernel since $\int_{\mathbb{R}^d} d\mu(\omega) = \int_{\mathbb{R}^d} \exp(-i\langle\omega, 0\rangle)d\mu(\omega) = k_0(0)$. and we can write the kernel as an expectation over a probability measure $\mu$. For all $x, z \in \mathbb{R}^d$

$$k_0(x - z) = \mathbf{E}_{\omega \sim \mu} [\exp(-i\langle\omega, x - z\rangle)].$$

Eventually, if $k$ is real valued we only write the real part, $k(x, z) = \mathbf{E}_{\omega \sim \mu}[\cos\langle\omega, x - z\rangle] = \mathbf{E}_{\omega \sim \mu}[\cos\langle\omega, z\rangle \cos\langle\omega, x\rangle + \sin\langle\omega, z\rangle \sin\langle\omega, x\rangle]$. Let $\bigoplus_{j=1}^D x_j$ denote the $Dd$-length column vector obtained by stacking vectors $x_j \in \mathbb{R}^d$. The feature map $\tilde{\varphi}: \mathbb{R}^d \to \mathbb{R}^{2D}$ defined as

$$\tilde{\varphi}(x) = \frac{1}{\sqrt{D}} \bigoplus_{j=1}^{D} \begin{pmatrix} \cos \langle x, \omega_j \rangle \\ \sin \langle x, \omega_j \rangle \end{pmatrix}, \quad \omega_j \sim \mathcal{F}^{-1}[k_0] \text{ i. i. d.} \quad (1)$$

is called a *Random Fourier Feature* (map). Each $\omega_j, j = 1, \dots, D$ is independently and identically sampled from the inverse Fourier transform $\mu$ of $k_0$. This Random Fourier Feature
<table><caption>Table 1: Mathematical symbols and their signification (part 1).</caption><thead><tr><th>Symbol</th><th>Meaning</th></tr></thead><tbody><tr><td>e ∈ 𝒳</td><td>The neutral element of the group 𝒳.</td></tr><tr><td>δ<sub>ij</sub></td><td>Kronecker delta function. δ<sub>ij</sub> = 0 if i ≠ j, 1 otherwise.</td></tr><tr><td>̂𝒳</td><td>The Pontryagin dual of 𝒳 when 𝒳 is a LCA group.</td></tr><tr><td>⟨·, ·⟩<sub>Y</sub></td><td>The canonical inner product of the Hilbert space Y.</td></tr><tr><td>‖·‖<sub>Y</sub></td><td>The canonical norm induced by the inner product of the Hilbert space Y.</td></tr><tr><td>F(𝒳; Y)</td><td>Topological vector space of functions from 𝒳 to Y.</td></tr><tr><td>C(𝒳; Y)</td><td>The topological vector subspace of F of continuous functions from 𝒳 to Y.</td></tr><tr><td>L(𝒳; Y)</td><td>The space bounded linear operator from a Hilbert space H to a Hilbert space Y.</td></tr><tr><td>‖·‖<sub>Y,Y'</sub></td><td>The operator norm ‖Γ‖<sub>Y,Y'</sub> = sup<sub>‖y‖<sub>Y</sub>=1</sub> ‖Γy‖<sub>Y'</sub> for all Γ ∈ L(Y, Y')</td></tr><tr><td>M<sub>m,n</sub>(ℝ)</td><td>The space of matrices of size (m, n).</td></tr><tr><td>L(Y)</td><td>The space of bounded linear operator from a Hilbert space Y to itself.</td></tr><tr><td>L<sub>+</sub>(Y)</td><td>The space of non-negative bounded linear operator from a Hilbert space H to itself.</td></tr><tr><td>B(𝒳)</td><td>Borel σ-algebra on a topological space 𝒳.</td></tr><tr><td>Leb(𝒳)</td><td>The Lebesgue measure of 𝒳.</td></tr><tr><td>Haar(𝒳)</td><td>A Haar measure of 𝒳.</td></tr><tr><td>Pr<sub>μ,ρ</sub>(𝒳)</td><td>A probability measure of 𝒳 whose Radon-Nikodym derivative (density) with respect to the measure μ is ρ.</td></tr><tr><td>F[·]</td><td>The Fourier Transform operator.</td></tr><tr><td>L<sup>p</sup>(𝒳, μ; Y)</td><td>The Banach space of ‖·‖<sub>Y</sub><sup>p</sup> (Bochner)-integrable function from (𝒳, B(𝒳), μ) to Y for p ∈ ℝ<sub>+</sub>. L<sup>p</sup>(𝒳, μ, ℝ) := L<sup>p</sup>(𝒳, μ) and L<sup>p</sup>(𝒳, μ, ℝ) = L<sup>p</sup>(𝒳, μ).</td></tr><tr><td>⊕<sup>D</sup><sub>j=1</sub> x<sub>i</sub></td><td>The direct sum of D ∈ N vectors x<sub>i</sub>'s in the Hilbert spaces H<sub>i</sub>. By definition ⟨⊕<sup>D</sup><sub>j=1</sub> x<sub>j</sub>, ⊕<sup>D</sup><sub>j=1</sub> z<sub>j</sub>⟩ = ∑<sup>D</sup><sub>j=1</sub> ⟨x<sub>j</sub>, z<sub>j</sub>⟩<sub>H<sub>i</sub></sub>.</td></tr><tr><td>‖·‖<sub>p</sub></td><td>The L<sup>p</sup>(𝒳, μ, Y) norm. ‖f‖<sup>p</sup><sub>p</sub> := ∫<sub>𝒳</sub> ‖f(x)‖<sup>p</sup><sub>Y</sub> dμ(x). When 𝒳 = N*, Y ⊆ ℝ and μ is the counting measure and p = 2 it coincide with the Euclidean norm ‖·‖<sub>2</sub> for finite dimensional vectors.</td></tr><tr><td>‖·‖<sub>∞</sub></td><td>The uniform norm ‖f‖<sub>∞</sub> = ess sup {‖f(x)‖<sub>Y</sub> | x ∈ 𝒳} = lim<sub>p→∞</sub> ‖f‖<sub>p</sub>.</td></tr><tr><td>|Γ|</td><td>The absolute value of the linear operator Γ ∈ L(Y), i. e. |Γ|<sup>2</sup> = Γ*Γ.</td></tr><tr><td>Tr[Γ]</td><td>The trace of a linear operator Γ ∈ L(Y).</td></tr><tr><td>‖·‖<sub>σ,p</sub></td><td>The Schatten p-norm, ‖Γ‖<sup>p</sup><sub>σ,p</sub> = Tr[|Γ|^p] for Γ ∈ L(Y), where Y is a Hilbert space. Note that ‖Γ‖<sub>σ,∞</sub> = ρ(Γ) ≤ ‖Γ‖<sub>Y,Y'</sub>.</td></tr><tr><td>≫</td><td>“Greater than” in the Loewner partial order of operators. Γ<sub>1</sub> ≫ Γ<sub>2</sub> if σ(Γ<sub>1</sub> − Γ<sub>2</sub>) ⊆ ℝ<sub>+</sub>.</td></tr><tr><td>≈</td><td>Given two sets 𝒳 and Y, 𝒳 ≅ Y if there exists an isomorphism ϕ : 𝒳 → Y.</td></tr></tbody></table>
map provides the following Monte-Carlo estimator of the kernel: $\tilde{k}(x, z) = \tilde{\varphi}(x) * \tilde{\varphi}(z)$. Using trigonometric identities, Rahimi and Recht (2007) showed that the same feature map can also be written

$$ \tilde{\varphi}(x) = \sqrt{\frac{2}{D}} \bigoplus_{j=1}^{D} \left( \cos(\langle x, \omega_j \rangle + b_j) \right), \qquad (2) $$

where $\omega_j \sim \mathcal{F}^{-1}[k_0]$, $b_j \sim \mathcal{U}(0, 2\pi)$ i. i. d.. The feature map defined by Equation 1 and Equation 2 have been compared in Sutherland and Schneider (2015) where they give the condition under which Equation 1 has lower variance than Equation 2. For instance for the Gaussian kernel, Equation 1 has always lower variance. In practice, Equation 2 is easier to program. In this paper we focus on random Fourier feature of the form given in Equation 1.

The dimension $D$ governs the precision of this approximation, whose uniform convergence towards the target kernel can be found in Rahimi and Recht (2007) and in more recent papers with some refinements proposed in Sutherland and Schneider (2015) and Sriperumbudur and Szabo (2015). Finally, it is important to notice that Random Fourier Feature approach only requires two steps before the application of a learning algorithm: (1) define the inverse Fourier transform of the given shift-invariant kernel, (2) compute the randomized feature map using the spectral distribution $\mu$. Rahimi and Recht (2007) show that for the Gaussian kernel $k_0(x - z) = \exp(-\gamma\|x - z\|_2)$, the spectral distribution $\mu$ is a Gaussian distribution. For the Laplacian kernel $k_0(x - z) = \exp(-\gamma\|x - z\|_1)$, the spectral distribution is a Cauchy distribution.

### 2.1.1 EXTENSIONS OF THE RFF METHOD

The seminal idea of Rahimi and Recht (2007) has opened a large literature on random features. Nowadays, many classes of kernels other than translation invariant are now proved to have an efficient random feature representation. Kar and Karnick (2012) proposed random feature maps for dot product kernels (rotation invariant) and Hamid et al. (2014) improved the rate of convergence of the approximation error for such kernels by noticing that feature maps for dot product kernels are usually low rank and may not utilize the capacity of the projected feature space efficiently. Pham and Pagh (2013) proposed fast random feature maps for polynomial kernels.

Li et al. (2010) generalized the original RFF of Rahimi and Recht (2007). Instead of computing feature maps for shift-invariant kernels on the additive group $(\mathbb{R}^d, +)$, they used the generalized Fourier transform on any locally compact abelian group to derive random features on the multiplicative group $(\mathbb{R}^d, *)$. In the same spirit Yang et al. (2014b) noticed that an theorem equivalent to Bochner's theorem exists on the semi-group $(\mathbb{R}^d, +)$. From this they derived “Random Laplace” features and used them to approximate kernels adapted to learn on histograms.

To speed-up the convergence rate of the random features approximation, Yang et al. (2014a) proposed to sample the random variable from a quasi Monte-Carlo sequence instead of i. i. d. random variables. Le et al. (2013) proposed the “Fastfood” algorithm to reduce the complexity of computing a RFF –using structured matrices and a fast Walsh-Hadarmard transform– from $O_t(Dd)$ to $O_t(D \log(d))$. More recently Felix et al. (2016) proposed also an
algorithm “SORF” to compute Gaussian RFF in $O_t(D \log(d))$ but with better convergence rates than “Fastfood” (Le et al., 2013). Mukuta and Harada (2016) proposed a data dependent feature map (comparable to the Nystroïm method) by estimating the distribution of the input data, and then finding the eigenfunction decomposition of Mercer’s integral operator associated to the kernel.

In the context of large scale learning and deep learning, Lu et al. (2014) showed that RFFs can achieve performances comparable to deep-learning methods by combining multiple kernel learning and composition of kernels along with a scalable parallel implementation. Dai et al. (2014) and Xie et al. (2015) combined RFFs and stochastic gradient descent to define an online learning algorithm called “Doubly stochastic gradient descent” adapted to large scale learning. Yang et al. (2015a) proposed and studied the idea of replacing the last fully interconnected layer of a deep convolutional neural network (LeCun et al., 1995) by the “Fastfood” implementation of RFFs.

Eventually Yang et al. (2015b) introduced the algorithm “À la Carte”, based on “Fastfood” which is able to learn the spectral distribution

## 2.2 On Operator-Valued Kernels

We now introduce the theory of Vector Valued Reproducing Kernel Hilbert Space (VV-RKHS) that provides a flexible framework to study and learn vector-valued functions. The fundations of the general theory of scalar kernels is mostly due to Aronszajn (1950) and provides a unifying point of view for the study of an important class of Hilbert spaces of real or complex valued functions. It has been first applied in the theory of partial differential equation. The theory of Operator-Valued Kernels (OVKs) which extends the scalar-valued kernel was first developped by Pedrick (1957) in his Ph. D Thesis. Since then it has been successfully applied to machine learning by many authors. In particular we introduce the notion of Operator-Valued Kernels following the propositions of Micchelli and Pontil (2005); Carmeli et al. (2006, 2010).

## 2.3 Definitions and properties

In machine learning the goal is often to find a function $f$ belonging to a class of functions $\mathcal{F}(\mathcal{X}; \mathcal{Y})$ that minimizes a criterion called the true risk. The class of functions we consider are functions living in a Hilbert space $\mathcal{H} \subset \mathcal{F}(\mathcal{X}; \mathcal{Y})$. The completeness allows to consider sequences of functions $f_n \in \mathcal{H}$ where the limit $f_n \to f$ is in $\mathcal{H}$. Moreover the existence of an inner product gives rise to a norm and also makes $\mathcal{H}$ a metric space.

Among all these functions $f \in \mathcal{H}$, we consider a subset of functions $f \in \mathcal{H}_K \subset \mathcal{H}$ such that the evaluation map $ev_x: f \mapsto f(x)$ is bounded for all $x$. i. e. such that $\|ev_x\|_K \le C_x \in \mathbb{R}$ for all $x$, where $\|\cdot\|_K^2 = \langle \cdot, \cdot \rangle_K := \langle \cdot, \cdot \rangle_{\mathcal{H}_K}$. For scalar valued kernels the evaluation map is a linear functional. Thus by Riesz's representation theorem there is an isomorphism between evaluating a function at a point and an inner product: $f(x) = ev_x f = \langle K_x, f \rangle_K$. From this we deduce the reproducing property $K(x, z) = \langle K_x, K_z \rangle_K$ which is the cornerstone of many proofs in machine learning and functional analysis. When dealing with vector-valued functions, the evaluation map $ev_x$ is no longer a linear functional, since it is vector-valued.
However, inspired by the theory of scalar valued kernel, many authors showed that if the evaluation map of functions with values in a Hilbert space $\mathcal{Y}$ is bounded, a similar reproducing property can be obtained; namely $\langle y', K(x, z)y \rangle = \langle K_x y', K_z y \rangle_K$ for all $y, y' \in \mathcal{Y}$. This motivates the following definition of a Vector Valued Reproducing Kernel Hilbert Space (VV-RKHS).

**Definition 2** (Vector Valued Reproducing Kernel Hilbert Space (Carmeli et al., 2006; Micchelli and Pontil, 2005))

Let $\mathcal{Y}$ be a (real or complex) Hilbert space. A Vector Valued Reproducing Kernel Hilbert Space on a locally compact second countable topological space $\mathcal{X}$ is a Hilbert space $\mathcal{H}$ such that

1. the elements of $\mathcal{H}$ are functions from $\mathcal{X}$ to $\mathcal{Y}$ (i. e. $\mathcal{H} \subset \mathcal{F}(\mathcal{X}, \mathcal{Y})$);

2. for all $x \in \mathcal{X}$, there exists a positive constant $C_x$ such that for all $f \in \mathcal{H}$ $\|f(x)\|_{\mathcal{Y}} \le C_x \|f\|_{\mathcal{H}}$.

Throughout this section we show that a VV-RKHS defines a unique positive-definite function called Operator-Valued Kernel (OVK) and conversely an OVK uniquely defines a VV-RKHS. The bijection between OVKs and VV-RKHSs has been first proved by Senkene and Tempel'man (1973) in 1973. In this introduction to OVKs we follow the definitions and most recent proofs of Carmeli et al. (2010).

**Definition 3** (Positive-definite Operator-Valued Kernel)

Given $\mathcal{X}$ a locally compact second countable topological space and $\mathcal{Y}$ a real Hilbert Space, a map $K: \mathcal{X} \times \mathcal{X} \to \mathcal{L}(\mathcal{Y})$ is called a positive-definite Operator-Valued Kernel kernel if $K(x, z) = K(z, x)^*$ and

$$ \sum_{i,j=1}^{N} \langle K(x_i, x_j)y_j, y_i \rangle_{\mathcal{Y}} \ge 0, \quad (3) $$

for all $N \in \mathbb{N}$, for all sequences of points $(x_i)_{i=1}^N$ in $\mathcal{X}^N$, and all sequences of points $(y_i)_{i=1}^N$ in $\mathcal{Y}^N$.

As in the scalar case any Vector Valued Reproducing Kernel Hilbert Space defines a unique positive-definite Operator-Valued Kernel and conversely a positive-definite Operator-Valued Kernel defines a unique Vector Valued Reproducing Kernel Hilbert Space.

**Proposition 4** ((Carmeli et al., 2006))

Given a Vector Valued Reproducing Kernel Hilbert Space there is a unique positive-definite Operator-Valued Kernel $K: \mathcal{X} \times \mathcal{X} \to \mathcal{L}(\mathcal{Y})$.

Given $x \in \mathcal{X}$, $K_x: \mathcal{Y} \to \mathcal{F}(\mathcal{X}; \mathcal{Y})$ denotes the linear operator whose action on a vector $y$ is the function $K_x y \in \mathcal{F}(\mathcal{X}; \mathcal{Y})$ defined for all $z \in \mathcal{X}$ by $K_x = \text{ev}_x^*$. As a consequence we have that

$$ K(x, z)y = \text{ev}_x \text{ev}_z^* y = K_x^* K_z y = (K_z y)(x). \quad (4) $$

Some direct consequences follow from the definition.
1. The kernel reproduces the value of a function $f \in \mathcal{H}$ at a point $x \in \mathcal{X}$ since for all $y \in \mathcal{Y}$ and $x \in \mathcal{X}$, $\text{ev}_x^* y = K_x y = K(\cdot, x)y$ such that $\langle f(x), y \rangle_\mathcal{Y} = \langle f, K(\cdot, x)y \rangle_\mathcal{H} = \langle K_x^* f, y \rangle_\mathcal{Y}$.

2. For all $x \in \mathcal{X}$ and all $f \in \mathcal{H}$, $\|f(x)\|_{\mathcal{Y}} \le \sqrt{\|K(x, x)\|_{\mathcal{Y}, \mathcal{Y}}} \|f\|_{\mathcal{H}}$. This comes from the fact that $\|K_x\|_{\mathcal{Y}, \mathcal{H}} = \|K_x^*\|_{\mathcal{H}, \mathcal{Y}} = \sqrt{\|K(x, x)\|_{\mathcal{Y}, \mathcal{Y}}}$ and the operator norm is sub-multiplicative.

Additionally given a positive-definite Operator-Valued Kernel, it defines a unique VV-RKHS.

**Proposition 5** ((Carmeli et al., 2006))

Given a positive-definite Operator-Valued Kernel $K: \mathcal{X} \times \mathcal{X} \to \mathcal{L}(\mathcal{Y})$, there is a unique Vector Valued Reproducing Kernel Hilbert Space $\mathcal{H}$ on $\mathcal{X}$ with reproducing kernel $K$.

Since a positive-definite Operator-Valued Kernel defines a unique Vector Valued Reproducing Kernel Hilbert Space (VV-RKHS) and conversely a VV-RKHS defines a unique Operator-Valued Kernel, we denote the Hilbert space $\mathcal{H}$ endowed with the scalar product $\langle \cdot, \cdot \rangle$ respectively $\mathcal{H}_K$ and $\langle \cdot, \cdot \rangle_K$. From now we refer to positive-definite Operator-Valued Kernels or reproducing Operator-Valued Kernels as Operator-Valued Kernels. As a consequence, given $K$ an Operator-Valued Kernel, define $K_x = K(\cdot, x)$ we have

$$K(x, z) = K_x^* K_z \quad \forall x, z \in \mathcal{X}, \qquad (5a)$$

$$\mathcal{H}_K = \overline{\text{span}} \{ K_{xy} \mid \forall x \in \mathcal{X}, \forall y \in \mathcal{Y} \}. \qquad (5b)$$

Where $\overline{\text{span}}$ is the closed span of a given set. Another way to describe functions of $\mathcal{H}_K$ consists in using a suitable feature map.

**Proposition 6** (Feature Operator (Carmeli et al., 2010))

Let $\mathcal{H}$ be any Hilbert space and $\phi: \mathcal{X} \to \mathcal{L}(\mathcal{Y}; \mathcal{H})$, with $\phi_x := \phi(x)$. Then the operator $W: \mathcal{H} \to \mathcal{F}(\mathcal{X}; \mathcal{Y})$ defined for all $g \in \mathcal{H}$, and for all $x \in \mathcal{X}$ by $(Wg)(x) = \phi_x^* g$ is a partial isometry from $\mathcal{H}$ onto the VV-RKHS $\mathcal{H}_K$ with reproducing kernel $K(x, z) = \phi_x^* \phi_z$, $\forall x, z \in \mathcal{X}$. $W^*W$ is the orthogonal projection onto $(\text{Ker } W)^\perp = \overline{\text{span}} \{ \phi_x y \mid \forall x \in \mathcal{X}, \forall y \in \mathcal{Y} \}$. Then $\|f\|_K = \inf \{ \|g\|_\mathcal{H} \mid \forall g \in \mathcal{H}, Wg = f \}$.

In this work we mainly focus on the kernel functions inducing a VV-RKHS of continuous functions. Such kernel are named $\mathcal{Y}$-Mercer kernels and generalize Mercer kernels.

**Definition 7** ($\mathcal{Y}$-Mercer kernel Carmeli et al. (2010))

A positive definite OVK $K: \mathcal{X} \times \mathcal{X} \to \mathcal{L}(\mathcal{Y})$ is called $\mathcal{Y}$-Mercer if the associated VV-RKHS $\mathcal{H}_K$ is a subspace of the space of continuous functions $\mathcal{C}(\mathcal{X}; \mathcal{Y})$.

## 2.4 Examples of Operator-Valued Kernels

In this subsection we list some Operator-Valued Kernels (OVKs) that have been used successfully in the literature. We do not recall the proof that the following kernels are well defined are refer the interested reader to the respective authors original work.
OVKs have been first introduced in Machine Learning to solve multi-task regression problems. Multi-task regression is encountered in many fields such as structured classification when classes belong to a hierarchy for instance. Instead of solving independently $p$ single output regression task, one would like to take advantage of the relationships between output variables when learning and making a decision.

**Proposition 8** (Decomposable kernel (Micheli and Glaunes, 2013))

Let $\Gamma$ be a non-negative operator of $\mathcal{L}_+(\mathcal{Y})$. $K$ is said to be a decomposable kernel<fnref n="1" /> if for all $(x, z) \in \mathcal{X}^2$, $K(x, z) := k(x, z)^\Gamma$, where $k$ is a scalar kernel.

When $\mathcal{Y} = \mathbb{R}^p$, the operator $\Gamma$ can be represented by a matrix which can be interpreted as encoding the relationships between the outputs coordinates. If a graph coding for the proximity between tasks is known, then Evgeniou et al. (2005); Baldassarre et al. (2010); Álvarez et al. (2012) showed that $\Gamma$ can be chosen equal to the pseudo inverse $L^\dagger$ of the graph Laplacian, such that the norm in $\mathcal{H}_K$ is a graph-regularizing penalty for the outputs (tasks). When no prior knowledge is available, $\Gamma$ can be learned with one of the algorithms proposed in the literature (Dinuzzo et al., 2011; Sindhwani et al., 2013; Lim et al., 2015a). Another interesting property of the decomposable kernel is its universality (a kernel which may approximate an arbitrary continuous target function uniformly on any compact subset of the input space). A reproducing kernel $K$ is said universal if the associated VV-RKHS $\mathcal{H}_K$ is dense in the space of continuous functions $C(\mathcal{X}, \mathcal{Y})$. The conditions for a kernel to be universal have been discussed in Caponnetto et al. (2008); Carmeli et al. (2010). In particular they show that a decomposable kernel is universal provided that the scalar kernel $k$ is universal and the operator $\Gamma$ is injective. Given $(e_k)_{k=1}^p$ a basis of $\mathcal{Y}$, we recall here how the matrix $\Gamma$ act as a regularizer between the components of the outputs $f_k = \langle f(\cdot), e_k \rangle_{\mathcal{Y}}$ of a function $f \in \mathcal{H}_K$. We prove a generalized version of Proposition 9 to any Operator-Valued Kernel in Subsection 3.6.

**Proposition 9** (Kernels and Regularizers (Álvarez et al., 2012))

Let $K(x, z) := k(x, z)^\Gamma$ for all $x, z \in \mathcal{X}$ be a decomposable kernel where $\Gamma$ is a matrix of size $p \times p$. Then for all $f \in \mathcal{H}_K$, $\|f\|_K = \sum_{i,j=1}^p (\Gamma^\dagger)_{ij} \langle f_i, f_j \rangle_k$ where $f_i(\cdot) = \langle f(\cdot), e_i \rangle_{\mathcal{Y}}$ (resp $f_j = \langle f(\cdot), e_j \rangle_{\mathcal{Y}}$), denotes the i-th (resp j-th) component of $f$.

Curl-free and divergence-free kernels provide an interesting application of operator-valued kernels (Macedo and Castro, 2008; Baldassarre et al., 2012; Micheli and Glaunes, 2013) to vector field learning, for which input and output spaces have the same dimensions ($d = p$). Applications cover shape deformation analysis (Micheli and Glaunes, 2013) and magnetic fields approximations (Wahlström et al., 2013). These kernels discussed in (Fuselier, 2006) allow encoding input-dependent similarities between vector-fields. An illustration of a synthetic 2D curl-free and divergence free fields are given respectively in Figure 1 and Figure 2. To obtain the curl-free field we took the gradient of a mixture of five two dimensional Gaussians (since the gradient of a potential is always curl-free). We generated the divergence-free field by taking the orthogonal of the curl-free field.

**Proposition 10** (Curl-free and Div-free kernel (Macedo and Castro, 2008))

<footnote n="1">1. Some authors also refer to as separable kernels.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Synthetic 2D curl-free field.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Synthetic 2D divergence-free field.</figcaption></figure>

Assume $\mathcal{X} = (\mathbb{R}^d, +)$ and $\mathcal{Y} = \mathbb{R}^p$ with $d = p$. The divergence-free kernel is defined as $K^{\text{div}}(x, z) = K_0^{\text{div}}(\delta) = (\nabla\nabla^\top - \Delta I)k_0(\delta)$ and the curl-free kernel as $K^{\text{curl}}(x, z) = K_0^{\text{curl}}(\delta) = -\nabla\nabla^\top k_0(\delta)$, where $\nabla$ is the gradient operator, $\nabla\nabla^\top$ is the Hessian operator and $\Delta$ is the Laplacian operator.

## 2.5 Shift-Invariant OVK on LCA groups

The main subjects of interest of the present paper are shift-invariant Operator-Valued Kernel. When referring to a shift-invariant OVK $K: \mathcal{X} \times \mathcal{X} \to \mathcal{L}(\mathcal{Y})$ we assume that $\mathcal{X}$ is a locally compact second countable topological group with identity $e$.

**Definition 11 (Shift-invariant OVK)**

A reproducing Operator-Valued Kernel $K: \mathcal{X} \times \mathcal{X} \to \mathcal{L}(\mathcal{Y})$ is called shift-invariant if for all $x, z, t \in \mathcal{X}$, $K(x \star t, z \star t) = K(x, z)$.

A shift-invariant kernel can be characterized by a function of one variable $K_e$ called the signature of $K$. Here $e$ denotes the neutral element of the LCA group $\mathcal{X}$ endowed with the binary group operation $\star$.

**Proposition 12 (Kernel signature (Carmeli et al., 2010))**

Let $K: \mathcal{X} \times \mathcal{X} \to \mathcal{L}(\mathcal{Y})$ be a reproducing kernel. The following conditions are equivalents.

1. $K$ is a positive-definite shift-invariant Operator-Valued Kernel.

2. There is a positive-definite function $K_e: \mathcal{X} \to \mathcal{L}(\mathcal{Y})$ such that $K(x, z) = K_e(z^{-1} \star x)$.

If one of the above conditions is satisfied, then $\|K(x, x)\|_{\mathcal{Y}} = \|K_e(e)\|_{\mathcal{Y}}$, $\forall x \in \mathcal{X}$.

We recall that if $K$ is a $\mathcal{Y}$-Mercer kernel, there is a function $K_e$ such that for all $x$ and $z \in \mathcal{X}$, $K(x, z) = K_e(x \star z^{-1})$. Then an OVK $K$ is $\mathcal{Y}$-Mercer if and only if for all $y \in \mathcal{Y}$, $K_e(\cdot)y \in \mathcal{C}(\mathcal{X}; \mathcal{Y})$. In other words a $\mathcal{Y}$-Mercer kernel is nothing but a functions whose signature is continuous and positive definite (Carmeli et al., 2010), which fulfil the conditions required for the “operator-valued” Bochner theorem to apply. Note that if $K$ is a shift invariant $\mathcal{Y}$-Mercer kernel, then $\mathcal{H}_K$ contains continuous bounded functions (Carmeli et al., 2010).
## 2.6 Some applications of Operator-valued kernels

We give here a non exhaustive list of works concerning Operator-Valued Kernels. A good review of Operator-Valued Kernels has been conducted in Álvarez et al. (2012). For a theoretical introduction to OVKs the interested reader can refer to the papers Carmeli et al. (2006); Caponnetto et al. (2008); Carmeli et al. (2010). Generalization bounds for OVK have been studied in Sindhwani et al. (2013); Kadri et al. (2015); Sangnier et al. (2016); Maurer (2016). Operator-valued Kernel Regression has first been studied in the context of Ridge Regression and Multi-task learning by Micchelli and Pontil (2005). Multi-task regression (Micchelli and Pontil, 2004) and structured multi-class classification (Dinuzzo et al., 2011; Minh et al., 2013b; Mroueh et al., 2012) are undoubtedly the first target applications for working in Vector Valued Reproducing Kernel Hilbert Space. Operator-Valued Kernels have been shown useful to provide a general framework for structured output prediction (Brouard et al., 2011, 2016a) with a link to Output Kernel Regression (Kadri et al., 2013). Beyond structured classification, other various applications such as link prediction, drug activity prediction or recently metabolite identification (Brouard et al., 2016b) and image colorization (Ha Quang et al., 2010) have been developed.

Macedo and Castro (2008); Baldassarre et al. (2012) showed the interest of spectral algorithms in Ridge regression and introduced vector field learning as a new multiple output task in Machine Learning community. Wahlström et al. (2013) applied vector field learning with OVK-based Gaussian processes to the reconstruction of magnetic fields (which are curl-free). The works of Kadri et al. (2010, 2015) have been the precursors of regression with functional values, opening a new avenue of applications. Appropriate algorithms devoted to online learning have been also derived by Audiffren and Kadri (2015). Kernel learning was addressed at least in two ways: first with using Multiple Kernel Learning in Kadri et al. (2012) and second, using various penalties, smooth ones in Dinuzzo et al. (2011); Ciliberto et al. (2015) for decomposable kernels and non smooth ones in Lim et al. (2015b) using proximal methods in the case of decomposable and transformable kernels. Dynamical modeling was tackled in the context of multivariate time series modelling in Lim et al. (2013); Sindhwani et al. (2013); Lim et al. (2015b) and as a generalization of Recursive Least Square Algorithm in Amblard and Kadri (2015). Sangnier et al. (2016) recently explored the minimization of a pinball loss under regularizing constraints induced by a well chosen decomposable kernel in order to handle joint quantile regression.

## 3. Main contribution: Operator Random Fourier Features

We present in this section a construction methodology devoted to shift-invariant $\mathcal{Y}$-Mercer operator-valued kernels defined on any Locally Compact Abelian (LCA) group, noted $(\mathcal{X}, \star)$, for some operation noted $\star$. This allows us to use the general context of Pontryagin duality for Fourier Transform of functions on LCA groups. Building upon a generalization of the celebrated Bochner's theorem for operator-valued measures, an operator-valued kernel is seen as the *Fourier Transform* of an operator-valued positive measure. From that result, we extend the principle of RFF for scalar-valued kernels and derive a general methodology to build Operator-valued Random Fourier Feature (ORFF) when operator-valued kernels are
shift-invariant according to the chosen group operation. Elements of this paper have been
developped in Brault et al. (2016).

We present a construction of feature maps called Operator-valued Random Fourier Feature (ORFF), such that $f: x \mapsto \tilde{\phi}(x)^*\theta$ is a continuous function that maps an arbitrary LCA group $\mathcal{X}$ as input space to an arbitrary output Hilbert space $\mathcal{Y}$. First we define a functional Fourier feature map, and then propose a Monte-Carlo sampling from this feature map to construct an approximation of a shift-invariant $\mathcal{Y}$-Mercer kernel. Then, we prove the convergence of the kernel approximation $\tilde{K}(x, z) = \tilde{\phi}(x)^* \tilde{\phi}(z)$ with high probability on compact subsets of the LCA $\mathcal{X}$. Eventually we conclude with some numerical experiments.

## 3.1 Theoretical study

The following proposition of Zhang et al. (2012); Neeb (1998) extends Bochner's theorem to any shift-invariant $\mathcal{Y}$-Mercer kernel. We give a short introduction to LCA groups and abstract harmonic analysis in appendix A. For the sake of simplicity, the reader can take $\overline{(x, \omega)} = \exp(i\langle x, \omega \rangle_2) = \exp(-i\langle x, \omega \rangle_2)$, when $x \in \mathcal{X} = (\mathbb{R}^d, +)$.

**Proposition 13** (Operator-valued Bochner's theorem (Zhang et al., 2012; Neeb, 1998))

If a function $K$ from $\mathcal{X} \times \mathcal{X}$ to $\mathcal{Y}$ is a shift-invariant $\mathcal{Y}$-Mercer kernel on $\mathcal{X}$, then there exists a unique positive projection-valued measure $\hat{Q}: \mathcal{B}(\mathcal{X}) \to \mathcal{L}_+(\mathcal{Y})$ such that for all $x, z \in \mathcal{X}$,

$$K(x, z) = \int_{\hat{\mathcal{X}}} \overline{(x \star z^{-1}, \omega)} d\hat{Q}(\omega), \qquad (6)$$

where $\hat{Q}$ belongs to the set of all the projection-valued measures of bounded variation on the $\sigma$-algebra of Borel subsets of $\hat{\mathcal{X}}$. Conversely, from any positive operator-valued measure $\hat{Q}$, a shift-invariant kernel $K$ can be defined by Equation 6.

Although this theorem is central to the spectral decomposition of shift-invariant $\mathcal{Y}$-Mercer OVK, the following results proved by Carmeli et al. (2010) provides insights about this decomposition that are more relevant in practice. It first gives the necessary conditions to build shift-invariant $\mathcal{Y}$-Mercer kernel with a pair $(A, \hat{\mu})$ where $A$ is an operator-valued function on $\hat{\mathcal{X}}$ and $\hat{\mu}$ is a real-valued positive measure on $\hat{\mathcal{X}}$. Note that obviously such a pair is not unique and the choice of this paper may have an impact on theoretical properties as well as practical computations. Secondly it also states that any OVK have such a spectral decomposition when $\mathcal{Y}$ is finite dimensional or $\mathcal{X}$.

**Proposition 14** (Carmeli et al. (2010))

Let $\hat{\mu}$ be a positive measure on $\mathcal{B}(\hat{\mathcal{X}})$ and $A: \hat{\mathcal{X}} \to \mathcal{L}(\mathcal{Y})$ such that $\langle A(\cdot)y, y' \rangle \in L^1(\mathcal{X}, \hat{\mu})$ for all $y, y' \in \mathcal{Y}$ and $A(\omega) \ge 0$ for $\hat{\mu}$-almost all $\omega \in \hat{\mathcal{X}}$. Then, for all $\delta \in \mathcal{X}$,

$$K_e(\delta) = \int_{\hat{\mathcal{X}}} \overline{(\delta, \omega)} A(\omega) d\hat{\mu}(\omega) \qquad (7)$$

is the kernel signature of a shift-invariant $\mathcal{Y}$-Mercer kernel $K$ such that $K(x, z) = K_e(x \star z^{-1})$. The VV-RKHS $\mathcal{H}_K$ is embedded in $L^2(\hat{\mathcal{X}}, \hat{\mu}; \mathcal{Y}')$ by means of the feature operator

$$(Wg)(x) = \int_{\hat{\mathcal{X}}} \overline{(x, \omega)} B(\omega) g(\omega) d\hat{\mu}(\omega), \qquad (8)$$
Where $B(\omega)B(\omega)^* = A(\omega)$ and both integrals converge in the weak sense. If $\mathcal{Y}$ is finite dimensional or $\mathcal{X}$ is compact, any shift-invariant kernel is of the above form for some pair $(A, \hat{\mu})$.

When $p=1$ one can always assume $A$ is reduced to the scalar $1$, $\hat{\mu}$ is still a bounded positive measure and we retrieve the Bochner theorem applied to the scalar case (Theorem 1).

Proposition 14 shows that a pair $(A, \hat{\mu})$ entirely characterizes an OVK. Namely a given measure $\hat{\mu}$ and a function $A$ such that $\langle y', A(c)y \rangle_y \in L^1(\mathcal{X}, \hat{\mu})$ for all $y, y' \in \mathcal{Y}$ and $A(\omega) \ge 0$ for $\hat{\mu}$-almost all $\omega$, give rise to an OVK. Since $(A, \hat{\mu})$ determine a unique kernel we can write $\mathcal{H}_{(A, \hat{\mu})} \Rightarrow \mathcal{H}_K$ where $K$ is defined as in Equation 7. However the converse is not true: Given a $\mathcal{Y}$-Mercer shift invariant Operator-Valued Kernel, there exist infinitely many pairs $(A, \hat{\mu})$ that characterize an OVK.

The main difference between Equation 6 and Equation 7 is that the first one characterizes an OVK by a unique Positive Operator-Valued Measure (POVM), while the second one shows that the POVM that uniquely characterizes a $\mathcal{Y}$-Mercer OVK has an operator-valued density with respect to a scalar measure $\hat{\mu}$; and that this operator-valued density is not unique. Notice that obtaining a scalar-valued density from the POVM implies the addition of a (weak) $L^1$ condition on $K_e$.

Finally Proposition 14 does not provide any constructive way to obtain the pair $(A, \hat{\mu})$ that characterizes an OVK. The following Subsection 3.1.1 is based on another proposition of Carmeli et al. and shows that if the kernel signature $K_e(\delta)$ of an OVK is in $L^1$ then it is possible to construct explicitly a pair $(C, \widehat{\text{Haar}})$ from it. Additionally, we show that we can always extract a scalar-valued probability density function from $C$ such that we obtain a pair $(A, \text{Pr}_{\hat{\mu},\rho})$ where $\text{Pr}_{\hat{\mu},\rho}$ is a probability distribution absolutely continuous with respect to $\hat{\mu}$ and with associated probability density function (p. d. f) $\rho$. Thus for all $\mathcal{Z} \subset \mathcal{B}(\hat{\mathcal{X}})$, $\text{Pr}_{\hat{\mu},\rho}(\mathcal{Z}) = \int_{\mathcal{Z}} \rho(\omega) d\hat{\mu}(\omega)$. When the reference measure $\hat{\mu}$ is the Lebesgue measure, we note $\text{Pr}_{\hat{\mu},\rho} = \text{Pr}_{\rho}$. For any function $f: \mathcal{X} \times \hat{\mathcal{X}} \times \mathcal{Y} \to \mathbb{R}$, we also use the notation $\text{E}_{\widehat{\text{Haar}},\rho}[f(x, \omega, y)] = \text{E}_{\omega \sim \text{Pr}_{\widehat{\text{Haar}},\rho}}[f(x, \omega, y)] = \int_{\hat{\mathcal{X}}} f(x, \omega, y) d\text{Pr}_{\widehat{\text{Haar}},\rho}(\omega) = \int_{\hat{\mathcal{X}}} f(x, \omega, y) \rho(\omega) d\widehat{\text{Haar}}(\omega)$. where the two last equalities hold by the “law of the unconscious statistician” (change of variable formula) and the fact that $\text{Pr}_{\widehat{\text{Haar}},\rho}$ has density $\rho$. Hence depending on the context, with mild abuse of notation, we view $\omega$ either as a variable in the probability space $(\hat{\mathcal{X}}, \text{Pr}_{\widehat{\text{Haar}},\rho})$ or as a random variable with values in $\hat{\mathcal{X}}$ and distributed as $\text{Pr}_{\widehat{\text{Haar}},\rho}$.

### 3.1.1 SUFFICIENT CONDITIONS OF EXISTENCE

While Proposition 14 gives some insights on how to build an approximation of a $\mathcal{Y}$-Mercer kernel, we need a theorem that provides an explicit construction of the pair $(A, \text{Pr}_{\hat{\mu},\rho})$ from the kernel signature $K_e$. Proposition 14 in Carmeli et al. (2010) gives the solution, and also provides a sufficient condition for Proposition 14 to apply.

**Proposition 15** (Carmeli et al. (2010))

Let $K$ be a shift-invariant $\mathcal{Y}$-Mercer kernel of signature $K_e$. Suppose that for all $z \in \mathcal{X}$ and for all $y, y' \in \mathcal{Y}$, the function $\langle K_e(\cdot)y, y' \rangle_y \in L^1(\mathcal{X}, \text{Haar})$ where $\mathcal{X}$ is endowed with the
group law ⋆. Denote $C : \hat{X} \to \mathcal{L}(Y)$, the function defined for all $\omega \in \hat{X}$ that satisfies for all $y, y' \text{ in } \mathcal{Y}$:

$$ \langle y', C(\omega)y \rangle_{\mathcal{Y}} = \int_{\mathcal{X}} (\delta, \omega) \langle y', K_e(\delta)y \rangle_{\mathcal{Y}} d\hat{\text{Haar}}(\delta) = \mathcal{F}^{-1}[\langle y', K_e(\cdot)y \rangle_{\mathcal{Y}}](\omega). \quad (9) $$

Then

1. $C(\omega)$ is a bounded non-negative operator for all $\omega \in \hat{\mathcal{X}}$,

2. $\langle y, C(\cdot)y' \rangle_{\mathcal{Y}} \in L^1(\hat{\mathcal{X}}, \widehat{\text{Haar}})$ for all $y, y' \in \mathcal{X}$,

3. for all $\delta \in \mathcal{X}$ and for all $y, y' \text{ in } \mathcal{Y}$, $\langle y', K_e(\delta)y \rangle_{\mathcal{Y}} = \int_{\hat{\mathcal{X}}} \overline{(\delta, \omega)} \langle y', C(\omega)y \rangle_{\mathcal{Y}} d\widehat{\text{Haar}}(\omega) = \mathcal{F}[\langle y', C(\cdot)y \rangle_{\mathcal{Y}}](\delta)$.

We found some confusion in the literature whether a kernel is the Fourier Transform or Inverse Fourier Transform of a measure. However Lemma 16 clarifies the relation between the Fourier Transform and Inverse Fourier Transform for a translation invariant Operator-Valued Kernel. Notice that in the real scalar case the Fourier Transform and Inverse Fourier Transform of a shift-invariant kernel are the same, while the difference is significant for OVK. The following lemma is a direct consequence of the definition of $C(\omega)$ as the Fourier Transform of the adjoint of $K_e$ and also helps to simplify the definition of ORFF.

**Lemma 16**

Let $K_e$ be the signature of a shift-invariant $\mathcal{Y}$-Mercer kernel such that for all $y, y' \in \mathcal{Y}$, $\langle y', K_e(\cdot)y \rangle_{\mathcal{Y}} \in L^1(\mathcal{X}, \widehat{\text{Haar}})$ and let $\langle y', C(\cdot)y \rangle_{\mathcal{Y}} = \mathcal{F}^{-1}[\langle y', K_e(\cdot)y \rangle_{\mathcal{Y}}]$. Then

1. $C(\omega)$ is self-adjoint and $C$ is even.

2. $\mathcal{F}^{-1}[\langle y', K_e(\cdot)y \rangle_{\mathcal{Y}}] = \mathcal{F}[\langle y', K_e(\cdot)y \rangle_{\mathcal{Y}}]$.

3. $K_e(\delta)$ is self-adjoint and $K_e$ is even.

While Proposition 15 gives an explicit form of the operator $C(\omega)$ defined as the Fourier Transform of the kernel $K$, it is not really convenient to work with the Haar measure $\widehat{\text{Haar}}$ on $\mathcal{B}(\hat{\mathcal{X}})$. However it is easily possible to turn $\widehat{\text{Haar}}$ into a probability measure to allow efficient integration over an infinite domain.

The following proposition allows to build a spectral decomposition of a shift-invariant $\mathcal{Y}$-Mercer kernel on a LCA group $\mathcal{X}$ endowed with the group law ⋆ with respect to a scalar probability measure, by extracting a scalar probability density function from $C$.

**Proposition 17** (Shift-invariant $\mathcal{Y}$-Mercer kernel spectral decomposition)

Let $K_e$ be the signature of a shift-invariant $\mathcal{Y}$-Mercer kernel. If for all $y, y' \in \mathcal{Y}$, $\langle K_e(\cdot)y, y' \rangle \in L^1(\mathcal{X}, \widehat{\text{Haar}})$ then there exists a positive probability measure $\text{Pr}_{\widehat{\text{Haar}},\rho}$ and an operator-valued function $A$ such that for all $y, y' \in \mathcal{Y}$,

$$ \langle y', K_e(\delta)y \rangle_{\mathcal{Y}} = \mathbf{E}_{\widehat{\text{Haar}},\rho} \left[ \overline{(\delta, \omega)} \langle y', A(\omega)y \rangle_{\mathcal{Y}} \right], \quad (10) $$

with $\langle y', A(\omega)y \rangle_{\mathcal{Y}} \rho(\omega) = \mathcal{F}[\langle y', K_e(\cdot)y \rangle_{\mathcal{Y}}](\omega)$. Moreover
1. for all $y, y' \in \mathcal{Y}, \langle A(\cdot)y, y'\rangle_{\mathcal{Y}} \in L^1(\hat{\mathcal{X}}, \mathbf{Pr}_{\widehat{\text{Haar}},\rho})$,

2. $A(\omega)$ is non-negative for $\mathbf{Pr}_{\widehat{\text{Haar}},\rho}$-almost all $\omega \in \hat{\mathcal{X}}$,

3. $A(\cdot)$ and $\rho(\cdot)$ are even functions.

## 3.2 Examples of spectral decomposition

In this section we give examples of spectral decomposition for various $\mathcal{Y}$-Mercer kernels, based on Proposition 17.

### 3.2.1 GAUSSIAN DECOMPOSABLE KERNEL

Recall that a decomposable $\mathbb{R}^p$-Mercer introduced in Proposition 8 has the form $K(x, z) = k(x, z)\Gamma$, where $k(x, z)$ is a scalar Mercer kernel and $\Gamma \in \mathcal{L}(\mathbb{R}^p)$ is a non-negative operator. Let us focus on $K_e^{\text{dec,gauss}}(\cdot) = k_e^{\text{gauss}}(\cdot)\Gamma$, the Gaussian decomposable kernel where $K_e^{\text{dec,gauss}}$ and $k_e^{\text{gauss}}$ are respectively the signature of $K$ and $k$ on the additive group $\mathcal{X} = (\mathbb{R}^d, +)$ - i. e. $\delta = x - z$ and $e = 0$. The well known Gaussian kernel is defined for all $\delta \in \mathbb{R}^d$ as follows $k_0^{\text{gauss}}(\delta) = \exp(-\sigma^{-2}\|\delta\|_2^2)/2$ where $\sigma \in \mathbb{R}_{>0}$ is an hyperparameter corresponding to the bandwith of the kernel. The -Pontryagin- dual group of $\mathcal{X} = (\mathbb{R}^d, +)$ is $\hat{\mathcal{X}} \cong (\mathbb{R}^d, +)$ with the pairing $(\delta, \omega) = \exp(i\langle\delta, \omega\rangle_2)$ where $\delta$ and $\omega \in \mathbb{R}^d$. In this case the Haar measures on $\mathcal{X}$ and $\hat{\mathcal{X}}$ are in both cases the Lebesgue measure. However in order to have the property that $\mathcal{F}^{-1}[\mathcal{F}[f]] = f$ and $\mathcal{F}^{-1}[f] = \mathcal{R}\mathcal{F}[f]$ one must normalize both measures by $\sqrt{2\pi}^{-d}$, i. e. for all $\mathcal{Z} \in \mathcal{B}(\mathbb{R}^d)$, $\sqrt{2\pi}^d\text{Haar}(\mathcal{Z}) = \text{Leb}(\mathcal{Z})$ and $\sqrt{2\pi}^d\widehat{\text{Haar}}(\mathcal{Z}) = \widehat{\text{Leb}}(\mathcal{Z})$. Then the Fourier Transform on $(\mathbb{R}^d, +)$ is

$$ \mathcal{F}[f](\omega) = \int_{\mathbb{R}^d} \exp(-i\langle\delta, \omega\rangle_2) f(\delta) d\text{Haar}(\delta) = \int_{\mathbb{R}^d} \exp(-i\langle\delta, \omega\rangle_2) f(\delta) \frac{d\text{Leb}(\delta)}{\sqrt{2\pi}^d}. $$

Since $k_0^{\text{gauss}} \in L^1$ and $\Gamma$ is bounded, it is possible to apply Proposition 17, and obtain for all $y$ and $y' \in \mathcal{Y}$,

$$ \langle y', C^{\text{dec,gauss}}(\omega)y \rangle_2 = \mathcal{F}\left[\langle y', K_0^{\text{dec,gauss}}(\cdot)y \rangle_2\right](\omega) = \mathcal{F}[k_0^{\text{gauss}}](\omega)\langle y', \Gamma y \rangle_2. $$

Thus

$$ C^{\text{dec,gauss}}(\omega) = \int_{\mathbb{R}^d} \exp \left( -i \langle \omega, \delta \rangle_2 - \frac{\|\delta\|_2^2}{2\sigma^2} \right) \frac{d\text{Leb}(\delta)}{\sqrt{2\pi}^d} \Gamma. $$

Hence

$$ C^{\text{dec,gauss}}(\omega) = \underbrace{\frac{1}{\sqrt{2\pi^{\frac{1}{\sigma^2}}}} \exp\left(-\frac{\sigma^2}{2}\|\omega\|_2^2\right)}_{A(\cdot)=\Gamma} \sqrt{2\pi}^d \underbrace{\Gamma}_{\rho(\cdot)=\mathcal{N}(0, \sigma^{-2}I_d)\sqrt{2\pi}^d}. $$

Therefore the canonical decomposition of $C^{\text{dec,gauss}}$ is $A^{\text{dec,gauss}}(\omega) = \Gamma$ and $\rho^{\text{dec,gauss}} = \mathcal{N}(0, \sigma^{-2}I_d)\sqrt{2\pi}^d$, where $\mathcal{N}$ is the Gaussian probability distribution. Note that this decomposition is done with respect to the normalized Lebesgue measure $\widehat{\text{Haar}}$, meaning that for
all $\mathcal{Z} \in \mathcal{B}(\hat{\mathcal{X}})$,

$$ 
\begin{aligned}
\mathbf{Pr}_{\widehat{\text{Haar}}, \mathcal{N}(0, \sigma^{-2} I_d) \sqrt{2\pi}^d} (\mathcal{Z}) &= \int_{\mathcal{Z}} \mathcal{N}(0, \sigma^{-2} I_d) \sqrt{2\pi}^d d\widehat{\text{Haar}}(\omega) \\
&= \int_{\hat{\mathcal{X}}} \mathcal{N}(0, \sigma^{-2} I_d) d\mathbf{Leb}(\omega) = \mathbf{Pr}_{\mathcal{N}(0, \sigma^{-2} I_d)}(\mathcal{Z}).
\end{aligned}
$$

Thus, the same decomposition with respect to the usual –non-normalized– Lebesgue measure **Leb** yields

$$ A^{\text{dec,gauss}}(\cdot) = \Gamma \qquad (11a) $$

$$ \rho^{\text{dec,gauss}} = \mathcal{N}(0, \sigma^{-2} I_d). \qquad (11b) $$

### 3.2.2 SKEWED-$\chi^2$ DECOMPOSABLE KERNEL

The skewed-$\chi^2$ scalar kernel (Li et al., 2010), useful for image processing, is defined on the LCA group $\mathcal{X} = (-c_k; +\infty)_{k=1}^d$, with $c_k \in \mathbb{R}_{>0}$ and endowed with the group operation $\odot$. Let $(e_k)_{k=1}^d$ be the standard basis of $\mathcal{X}$. The operator $\odot: \mathcal{X} \times \mathcal{X} \to \mathcal{X}$ is defined by $x \odot z = ((x_k + c_k)(z_k + c_k) - c_k)_{k=1}^d$. The identity element $e$ is $(1 - c_k)_{k=1}^d$ since $(1 - c) \odot x = x$. Thus the inverse element $x^{-1}$ is $((x_k + c_k)^{-1} - c_k)_{k=1}^d$. The skewed-$\chi^2$ scalar kernel reads

$$ k_{1-c}^{\text{skewed}}(\delta) = \prod_{k=1}^{d} \frac{2}{\sqrt{\delta_k + c_k} + \sqrt{\frac{1}{\delta_k + c_k}}}. \qquad (12) $$

The dual of $\mathcal{X}$ is $\hat{\mathcal{X}} \cong \mathbb{R}^d$ with the pairing $(\delta, \omega) = \prod_{k=1}^d \exp(i \log(\delta_k + c_k)\omega_k)$. The Haar measure are defined for all $\mathcal{Z} \in \mathcal{B}((-c; +\infty)^d)$ and all $\hat{\mathcal{Z}} \in \mathcal{B}(\mathbb{R}^d)$ by $\sqrt{2\pi}^d \mathbf{Haar}(\mathcal{Z}) = \int_{\mathcal{Z}} \prod_{k=1}^d \frac{1}{z_k + c_k} d\mathbf{Leb}(z)$ and $\sqrt{2\pi}^d \widehat{\mathbf{Haar}}(\hat{\mathcal{Z}}) = \mathbf{Leb}(\hat{\mathcal{Z}})$. Thus the Fourier Transform is

$$ \mathcal{F}[f](\omega) = \int_{(-c; +\infty)^d} \prod_{k=1}^{d} \frac{\exp(-i \log(\delta_k + c_k)\omega_k)}{\delta_k + c_k} f(\delta) \frac{d\mathbf{Leb}(\delta)}{\sqrt{2\pi}^d}. $$

Then, applying Fubini's theorem over product space, and the fact that each dimension is independent

$$ \mathcal{F}[k_0^{\text{skewed}}](\omega) = \prod_{k=1}^{d} \int_{-c_k}^{+\infty} \frac{2 \exp(-i \log(\delta_k + c_k)\omega_k)}{(\delta_k + c_k) \left(\sqrt{\delta_k + c_k} + \sqrt{\frac{1}{\delta_k + c_k}}\right)} \frac{d\mathbf{Leb}(\delta_k)}{\sqrt{2\pi}^d}. $$

Making the change of variable $t_k = (\delta_k + c_k)^{-1}$ yields

$$ \mathcal{F}[k_0^{\text{skewed}}](\omega) = \prod_{k=1}^{d} \int_{-\infty}^{+\infty} \frac{2 \exp(-i t_k \omega_k)}{\exp\left(\frac{1}{2} t_k\right) + \exp\left(-\frac{1}{2} t_k\right)} \frac{d\mathbf{Leb}(t_k)}{\sqrt{2\pi}^d} = \sqrt{2\pi}^d \prod_{k=1}^{d} \text{sech}(\pi \omega_k). $$

Since $k_{1-c}^{\text{skewed}} \in L^1$ and $\Gamma$ is bounded, it is possible to apply Proposition 17, and obtain

$$ C^{\text{dec,skewed}}(\omega) = \mathcal{F}[k_{1-c}^{\text{skewed}}](\omega)\Gamma = \underbrace{\sqrt{2\pi}^d \prod_{k=1}^d \text{sech}(\pi\omega_k)}_{\rho(\cdot) = \mathcal{S}(0, 2^{-1})^d \sqrt{2\pi}^d} A(\cdot). $$
Hence the decomposition with respect to the usual –non-normalized– Lebesgue measure Leb yields

$$A^{\text{dec,skewed}}(\cdot) = \Gamma \qquad (13a)$$

$$\rho^{\text{dec,skewed}} = S (0, 2^{-1})^d. \qquad (13b)$$

### 3.2.3 CURL-FREE GAUSSIAN KERNEL

The curl-free Gaussian kernel is defined as $K_0^{\text{curl,gauss}} = -\nabla\nabla^T k_0^{\text{gauss}}$. Here $\mathcal{X} = (\mathbb{R}^d, +)$ so the setting is the same than Subsection 3.2.1.

$$C^{\text{curl,gauss}}(\omega)_{ij} = \mathcal{F}\left[K_{1-c}^{\text{curl,gauss}}(\cdot)_{ij}\right](\omega) = \sqrt{2\pi \frac{1}{\sigma^2}} \exp\left(-\frac{\sigma^2}{2}\|\omega\|_2^2\right) \sqrt{2\pi} \omega_i \omega_j.$$

Hence

$$C^{\text{curl,gauss}}(\omega) = \underbrace{\frac{1}{\sqrt{2\pi \frac{1}{\sigma^2}}} \exp\left(-\frac{\sigma^2}{2}\|\omega\|_2^2\right) \sqrt{2\pi}^d}_{\mu(\cdot)=\mathcal{N}(0,\sigma^{-2}I_d)\sqrt{2\pi}^d} \underbrace{\omega\omega^\top}_{A(\omega)=\omega\omega^\top}.$$

Here a canonical decomposition is $A^{\text{curl,gauss}}(\omega) = \omega\omega^\top$ for all $\omega \in \mathbb{R}^d$ and $\mu^{\text{curl,gauss}} = \mathcal{N}(0, \sigma^{-2}I_d)\sqrt{2\pi}^d$ with respect to the normalized Lebesgue measure $d\omega$. Again the decomposition with respect to the usual –non-normalized– Lebesgue measure is for all $\omega \in \mathbb{R}^d$

$$A^{\text{curl,gauss}}(\omega) = \omega\omega^\top \qquad (14a)$$

$$\mu^{\text{curl,gauss}} = \mathcal{N}(0, \sigma^{-2}I_d). \qquad (14b)$$

### 3.2.4 DIVERGENCE-FREE KERNEL

The divergence-free Gaussian kernel is defined as $K_0^{\text{div,gauss}} = (\nabla\nabla^T - \Delta)k_0^{\text{gauss}}$ on the group $\mathcal{X} = (\mathbb{R}^d, +)$. The setting is the same than Subsection 3.2.1. Hence

$$C^{\text{div,gauss}}(\omega)_{ij} = \mathcal{F}\left[K_0^{\text{div,gauss}}(\cdot)_{ij}\right](\omega) = \left(\delta_{i=j} \sum_{k=1}^{d} \omega_k^2 - \omega_i \omega_j\right) \mathcal{F}[k_0^{\text{gauss}}](\omega).$$

Hence

$$C^{\text{div,gauss}}(\omega) = \underbrace{\frac{1}{\sqrt{2\pi \frac{1}{\sigma^2}}} \exp\left(-\frac{\sigma^2}{2}\|\omega\|_2^2\right) \sqrt{2\pi}^d}_{\rho(\cdot)=\mathcal{N}(0,\sigma^{-2}I_d)\sqrt{2\pi}^d} \underbrace{\left(I_d\|\omega\|_2^2 - \omega\omega^\top\right)}_{A(\omega)=I_d\|\omega\|_2^2-\omega\omega^\top}.$$

Thus the canonical decomposition with respect to the normalized Lebesgue measure is $A^{\text{div,gauss}}(\omega) = I_d\|\omega\|_2^2 - \omega\omega^\top$ and the measure $\rho^{\text{div,gauss}} = \mathcal{N}(0, \sigma^{-2}I_d)\sqrt{2\pi}^d$. The canonical decomposition with respect to the usual Lebesgue measure is

$$A^{\text{div,gauss}}(\omega) = I_d\|\omega\|_2^2 - \omega\omega^\top \qquad (15a)$$

$$\rho^{\text{div,gauss}} = \mathcal{N}(0, \sigma^{-2}I_d). \qquad (15b)$$
### 3.3 Operator-valued Random Fourier Features (ORFF)

#### 3.3.1 BUILDING OPERATOR-VALUED RANDOM FOURIER FEATURES

As shown in Proposition 17 it is always possible to find a pair $(A, \mathbf{Pr}_{\widehat{Haar},\rho})$ from a shift invariant $\mathcal{Y}$-Mercer Operator-Valued Kernel $K_e$ such that $\mathbf{Pr}_{\widehat{Haar},\rho}$ is a probability measure, i. e. $\int_{\hat{\mathcal{X}}} \rho d\widehat{\mathbf{Haar}} = 1$ where $\rho$ is the density of $\mathbf{Pr}_{\widehat{Haar},\rho}$ and $K_e(\delta) = \mathbf{E}_{\widehat{Haar},\rho}(\delta, \omega)A(\omega)$. In order to obtain an approximation of $K$ from a decomposition $(A, \mathbf{Pr}_{\widehat{Haar},\rho})$ we turn our attention to a Monte-Carlo estimation of the expectation in Equation 10 characterizing a $\mathcal{Y}$-Mercer shift-invariant Operator-Valued Kernel.

**Proposition 18**

Let $K(x, z)$ be a shift-invariant $\mathcal{Y}$-Mercer kernel with signature $K_e$ such that for all $y, y' \in \mathcal{Y}$, $\langle y', K_e(\cdot)y \rangle_{\mathcal{Y}} \in L^1(\mathcal{X}, \mathbf{Haar})$. Then one can find a pair $(A, \mathbf{Pr}_{\widehat{Haar},\rho})$ that satisfies Proposition 17. i. e. for $\mathbf{Pr}_{\widehat{Haar},\rho}$-almost all $\omega$, and all $y, y' \in \mathcal{Y}$, $\langle y, A(\omega)y' \rangle_{\rho}(\omega) = \mathcal{F}[\langle y', K_e(\cdot)y \rangle](\omega)$. If $(\omega_j)_{j=1}^D$ be a sequence of $D \in \mathbb{N}^*$ i. i. d. random variables following the law $\mathbf{Pr}_{\widehat{Haar},\rho}$ then the operator-valued function $\tilde{K}$ defined for $(x, z) \in \mathcal{X} \times \mathcal{X}$ as

$$ \tilde{K}(x, z) = \frac{1}{D} \sum_{j=1}^{D} \overline{(x \star z^{-1}, \omega_j)} A(\omega_j) $$

is an approximation of $K$. i. e. it satisfies for all $x, z \in \mathcal{X}$, $\tilde{K}(x, z) \xrightarrow{a. s.} K(x, z)$ in the weak operator topology, where $K$ is a $\mathcal{Y}$-Mercer OVK.

Now, for efficient computations as motivated in the introduction, we are interested in finding an approximated feature map instead of a kernel approximation. Indeed, an approximated feature map will allow to build linear models in regression tasks. The following proposition deals with the feature map construction.

**Proposition 19**

Assume the same conditions as Proposition 18. Moreover, if one can define $B : \hat{\mathcal{X}} \to \mathcal{L}(\mathcal{Y}', \mathcal{Y})$ such that for $\mathbf{Pr}_{\widehat{Haar},\rho}$-almost all $\omega$, and all $y, y' \in \mathcal{Y}$, $\langle y, B(\omega)B(\omega)^*y' \rangle_{\mathcal{Y}} \rho(\omega) = \langle y, A(\omega)y' \rangle \rho(\omega) = \mathcal{F}[\langle y, K_e(\cdot)y' \rangle_{\mathcal{Y}}](\omega)$, then the function $\tilde{\phi} : \hat{\mathcal{X}} \to \mathcal{L}(\mathcal{Y}, \bigoplus_{j=1}^D \mathcal{Y}')$ defined for all $y \in \mathcal{Y}$ as follows:

$$ \tilde{\phi}(x)y = \frac{1}{\sqrt{D}} \bigoplus_{j=1}^{D} (x, \omega_j) B(\omega_j)^* y, \quad \omega_j \sim \mathbf{Pr}_{\widehat{Haar},\rho} \quad \text{i. i. d.,} $$

is an approximated feature map<fnref n="2" /> for the kernel $K$.

**Remark 20**

<footnote n="2">2. i. e. it satisfies for all $x, z \in \mathcal{X}$, $\tilde{\phi}(x)^* \tilde{\phi}(z) \xrightarrow{a. s.} K(x, z)$ in the weak operator topology, where $K$ is a $\mathcal{Y}$-Mercer OVK.</footnote>
<table><caption>**Algorithm 1:** Construction of ORFF from OVK</caption>
  <tbody>
    <tr>
      <td>Input</td>
      <td>: K(x, z) = K<sub>e</sub>(δ) a shift-invariant <i>Y</i>-Mercer kernel on (<i>X</i>, ⋆) such that<br/>&forall;y, y' &isin; <i>Y</i>, &lang;y', K<sub>e</sub>(&middot;)y&rang;<sub>y</sub> &isin; L<sup>1</sup>(<b>R</b><sup>d</sup>, <b>Haar</b>) and D the number of features.</td>
    </tr>
    <tr>
      <td>Output</td>
      <td>: A random feature <i>&phi;&#x0303;</i>(x) such that <i>&phi;&#x0303;</i>(x)<sup>*</sup><i>&phi;&#x0303;</i>(z) &asymp; K(x, z)</td>
    </tr>
    <tr>
      <td>1</td>
      <td>Define the pairing (x, ω) from the LCA group (<i>X</i>, ⋆);</td>
    </tr>
    <tr>
      <td>2</td>
      <td>Find a decomposition (<i>A</i>, <b>Pr</b><sub><b>Haar</b>,<i>ρ</i></sub>) and <i>B</i> such that</td>
    </tr>
    <tr>
      <td></td>
      <td><i>B</i>(ω)<i>B</i>(ω)<sup>*</sup><i>ρ</i>(ω) = <i>A</i>(ω)<i>ρ</i>(ω) = <i>F</i><sup>-1</sup>[<i>K</i><sub>e</sub>](ω);</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Draw <i>D</i> i. i. d. realizations (<i>ω</i><sub>j</sub>)<sup>D</sup><sub>j=1</sub> from the probability distribution <b>Pr</b><sub><b>Haar</b>,<i>ρ</i></sub>;</td>
    </tr>
    <tr>
      <td>6</td>
      <td>return
        <span style="font-size: 2.5em; vertical-align: middle; line-height: 1;">{</span>
        <table style="display: inline-table; vertical-align: middle; border: none; text-align: center;">
          <tbody>
            <tr>
              <td style="border: none; padding: 0 0.2em; line-height: 1.1em;"><i>&phi;&#x0303;</i>(x) &isin; <i>L</i>(<i>Y</i>, <b><i>&#x0303;H</i></b>)</td>
              <td style="border: none; padding: 0 0.2em; line-height: 1.1em;"><i>y</i> &mapsto; <sup>1</sup>&frasl;<sub>&radic;<i>D</i></sub> &oplus;<sup><i>D</i></sup><sub><i>j</i>=1</sub>(<i>x</i>, <i>ω</i><sub><i>j</i></sub>)<i>B</i>(<i>ω</i><sub><i>j</i></sub>)<sup>*</sup><i>y</i></td>
            </tr>
            <tr>
              <td style="border: none; padding: 0 0.2em; line-height: 1.1em;"><i>&phi;&#x0303;</i>(x)<sup>*</sup> &isin; <i>L</i>(<b><i>&#x0303;H</i></b>, <i>Y</i>)</td>
              <td style="border: none; padding: 0 0.2em; line-height: 1.1em;"><i>&theta;</i> &mapsto; <sup>1</sup>&frasl;<sub>&radic;<i>D</i></sub> &sum;<sup><i>D</i></sup><sub><i>j</i>=1</sub>(<i>x</i>, <i>ω</i><sub><i>j</i></sub>)<i>B</i>(<i>ω</i><sub><i>j</i></sub>)<i>&theta;</i><sub><i>j</i></sub></td>
            </tr>
          </tbody>
        </table>
        <span style="font-size: 2.5em; vertical-align: middle; line-height: 1;">}</span>
      </td>
    </tr>
  </tbody>
</table>

We find a decomposition such that $A(\omega_j) = B(\omega_j)B(\omega_j)^*$ for all $j \in \mathbb{N}_D^*$ either by exhibiting a closed-form or using a numerical decomposition. Such a decomposition always exists since $A(\omega)$ is positive semi-definite for all $\omega \in \tilde{\mathcal{X}}$.

Notice that an ORFF map as defined in Proposition 19 is also the Monte-Carlo sampling of the corresponding functional Fourier feature map $\phi_x : \mathcal{Y} \to L^2(\hat{\mathcal{X}}, \mathbf{Pr}_{\widehat{\text{Haar}},\rho}; \mathcal{Y}')$ as defined in Proposition 21. Indeed, for all $y \in \mathcal{Y}$ and all $x \in \mathcal{X}$, $\tilde{\phi}(x)y = \bigoplus_{j=1}^D (\phi_x y)(\omega_j)$, $\omega_j \sim \mathbf{Pr}_{\widehat{\text{Haar}},\rho}$ i. i. d. Proposition 19 allows us to define Algorithm 1 for constructing ORFF from an operator valued kernel.

We give a numerical illustration of different $\tilde{K}$ built from different i. i. d. realization $(\omega_j)_{j=1}^D$, $\omega_j \sim \mathbf{Pr}_{\widehat{\text{Haar}},\rho}$. In Figure 3, we represent the approximation of a reference function (black line) defined as $(y_1, y_2)^T = f(x_i) = \sum_{j=1}^{250} K_{ij} u_j$ where $u_j \sim \mathcal{N}(0, I_2)$ and $K$ is a Gaussian decomposable kernel. We took $\Gamma = .5I_2 + .5I_2$ such that the outputs $y_1$ and $y_2$ share some similarities. We generated 250 points equally separated on the segment $(-1; 1)$. The Gram matrix is then $K_{ij} = \exp\left(-\frac{(x_i - x_j)^2}{2(0.1)^2}\right) \Gamma$, for $i, j \in \mathbb{N}_{250}$. We took $\Gamma = .5I_2 + .5I_2$ such that the outputs $y_1$ and $y_2$ share some similarities. We generated 250 points equally separated on the segment $(-1; 1)$. Then we computed an approximate kernel matrix $\tilde{K} \approx K$ for 25 increasing values of $D$ ranging from 1 to $10^4$. The two graphs in Figure 3 on the top row shows that the more the number of features increases the closer the model $\tilde{f}(x_i) = \sum_{j=1}^{250} \tilde{K}_{ij} u_j$ is to $f$. The bottom row shows the same experiment but for a different realization of $\tilde{K}$. When $D$ is small the curves of the bottom and top rows are very dissimilar –and sine wave like– while they both converge to $f$ when $D$ increase. We introduce a *functional* feature map, we call *Fourier Feature map*, defined by the following proposition as a direct consequence of Proposition 14.

**Proposition 21 (Functional Fourier feature map)**

Let $\mathcal{Y}$ and $\mathcal{Y}'$ be two Hilbert spaces. If there exists an operator-valued function $B : \hat{\mathcal{X}} \to \mathcal{L}(\mathcal{Y}, \mathcal{Y}')$ such that for all $y, y' \in \mathcal{Y}, \langle y, B(\omega)B(\omega)^*y' \rangle_{\mathcal{Y}} = \langle y', A(\omega)y \rangle_{\mathcal{Y}} \hat{\mu}$-almost everywhere
<figure><img src="image_4.png" /><figcaption>Figure 3: Approximation of a function in a VV-RKHS using different realizations of Operator Random Fourier Feature. Top row and bottom row correspond to two different realizations of $\tilde{K}$, which are different Operator-Valued Kernel. However when $D$ tends to infinity, the different realizations of $\tilde{K}$ yield the same OVK.</figcaption></figure>

and $\langle y', A(\cdot)y \rangle \in L^1(\hat{\chi}, \hat{\mu})$ then the operator $\phi_x$ defined for all $y$ in $\mathcal{Y}$ by $(\phi_x y)(\omega) = (x, \omega)B(\omega)^*y$, is a feature map<fnref n="3" /> of some shift-invariant $\mathcal{Y}$-Mercer kernel $K$.

With this notation we have $\phi: \mathcal{X} \to \mathcal{L}(\mathcal{Y}; L^2(\hat{\chi}, \hat{\mu}; \mathcal{Y}))$ such that $\phi_x \in \mathcal{L}(\mathcal{Y}; L^2(\hat{\chi}, \hat{\mu}; \mathcal{Y}))$ where $\phi_x := \phi(x)$. Notice that an ORFF map as defined in Proposition 19 is also the Monte-Carlo sampling of the corresponding functional Fourier feature map $\phi_x: \mathcal{Y} \to L^2(\hat{\chi}, \Pr_{\widehat{\text{Haar}},\rho}; \mathcal{Y})$ as defined in Proposition 21. Indeed, for all $y \in \mathcal{Y}$ and all $x \in \mathcal{X}$,

$$ \tilde{\phi}(x)y = \bigoplus_{j=1}^{D} (\phi_x y)(\omega_j), \quad \omega_j \sim \Pr_{\widehat{\text{Haar}},\rho} \quad \text{i. i. d.} $$

## 3.4 From Operator Random Fourier Feature maps to OVKs

It is also interesting to notice that we can go the other way and define from the general form of an Operator-valued Random Fourier Feature, an operator-valued kernel.

### Proposition 22 (Operator Random Fourier Feature map)

Let $\mathcal{Y}$ and $\mathcal{Y}'$ be two Hilbert spaces. If one defines an operator-valued function on the dual of a LCA group $\mathcal{X}$, $B: \hat{\chi} \to \mathcal{L}(\mathcal{Y}, \mathcal{Y}')$, and a probability measure $\Pr_{\widehat{\text{Haar}},\rho}$ on $B(\hat{\chi})$, such

<footnote n="3">3. i. e. it satisfies for all $x, z \in \mathcal{X}$, $\phi_x^*\phi_z = K(x, z)$ where $K$ is a $\mathcal{Y}$-Mercer OVK.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4: Relationships between feature-maps. For any realization of $\omega_j \sim \text{Pr}_{\widehat{\text{Haar}},\rho}$ i. i. d., $\tilde{\mathcal{H}} = \bigoplus_{j=1}^D \mathcal{Y}'$.</figcaption></figure>

that for all $y \in \mathcal{Y}$ and all $y' \in \mathcal{Y}'$, $\langle y, B(\cdot)y' \rangle \in L^2(\hat{\mathcal{X}}, \text{Pr}_{\widehat{\text{Haar}},\rho})$, then the operator-valued function $\tilde{\phi}: \mathcal{X} \to \mathcal{L}(\mathcal{Y}, \bigoplus_{j=1}^D \mathcal{Y}')$ defined for all $x \in \mathcal{X}$ and for all $y \in \mathcal{Y}$ by

$$ \tilde{\phi}(x)y = \frac{1}{\sqrt{D}} \bigoplus_{j=1}^{D} (x, \omega_j) B(\omega_j)^* y, \quad \omega_j \sim \text{Pr}_{\widehat{\text{Haar}},\rho}, \quad \text{i. i. d.}, \qquad (16) $$

is an approximated feature map of some $\mathcal{Y}$-Mercer operator-valued kernel$^4$.

The difference between Proposition 22 and Proposition 19 is that in Proposition 22 we do not assume that $A(\omega)$ and $\text{Pr}_{\widehat{\text{Haar}},\rho}$ have been obtained from Proposition 17. We conclude by showing that any realization of an approximate feature map gives a proper operator valued kernel. Hence we can always view $\tilde{K}(x, z) = \tilde{\phi}(x)^*\tilde{\phi}(z)$ —where $\tilde{\phi}$ is defined as in Proposition 18 (construction from an OVK) or Proposition 22— as a $\mathcal{Y}$-Mercer and thus apply all the classic results of the Operator-Valued Kernel theory on $\tilde{K}$.

**Proposition 23**

Let $\omega \in \hat{\mathcal{X}}^D$. If for all $y, y' \in \mathcal{Y}$ $\langle y', \tilde{K}_e(x \star z^{-1}) y \rangle_{\mathcal{Y}} = \langle \tilde{\phi}(x)y', \tilde{\phi}(z)y \rangle_{\tilde{\mathcal{H}}} = \langle y', \frac{1}{D} \sum_{j=1}^D \overline{(x \star z^{-1}, \omega_j)} B(\omega_j) B(\omega_j)^* y \rangle_{\mathcal{Y}}$, for all $x, z \in \mathcal{X}$, then $\tilde{K}$ is a shift-invariant $\mathcal{Y}$-Mercer Operator-Valued Kernel.

Note that the above theorem does not consider the $\omega_j$'s as random variables and therefore does not shows the convergence of the kernel $\tilde{K}$ to some target kernel $K$. However is shows that any realization of $\tilde{K}$ when $\omega_j$'s are random variables yields a valid $\mathcal{Y}$-Mercer operator-valued kernel. Note that the above theorem does not considers the $\omega_j$'s as random variables and therefore does not shows the convergence of the kernel $\tilde{K}$ to some target kernel $K$. However is shows that any realization of $\tilde{K}$ when $\omega_j$'s are random variables yields a

<footnote n="4">4. i. e. it satisfies $\tilde{\phi}(x)^*\tilde{\phi}(z) \xrightarrow{a. s.} K(x, z)$ in the weak operator topology, where $K$ is a $\mathcal{Y}$-Mercer OVK</footnote>
valid $\mathcal{Y}$-Mercer operator-valued kernel. Indeed, as a result of Proposition 23, in the same way we defined an ORFF, we can define an approximate feature operator $\tilde{W}$ which maps $\tilde{\mathcal{H}}$ onto $\mathcal{H}_{\tilde{K}}$, where $\tilde{K}(x, z) = \tilde{\phi}(x)^*\tilde{\phi}(z)$, for all $x, z \in \mathcal{X}$.

**Definition 24 (Random Fourier feature operator)**

Let $\omega = (\omega_j)_{j=1}^D \in \hat{\mathcal{X}}^D$ and let $\tilde{K}_e = \frac{1}{D} \sum_{j=1}^D \overline{(\cdot, \omega_j)} B(\omega_j) B(\omega_j)^*$. We call random Fourier feature operator the linear application $\tilde{W} : \tilde{\mathcal{H}} \to \mathcal{H}_{\tilde{K}}$ defined as

$$ (\tilde{W}\theta)(x) := \tilde{\phi}(x)^*\theta = \frac{1}{\sqrt{D}} \sum_{j=1}^{D} \overline{(x, \omega_j)} B(\omega_j)\theta_j $$

where $\theta = \bigoplus_{j=1}^D \theta_j \in \tilde{\mathcal{H}}$. Then from Proposition 6, $(\mathrm{Ker} \tilde{W})^\perp = \overline{\mathrm{span}} \{\tilde{\phi}(x)y \mid \forall x \in \mathcal{X}, \forall y \in \mathcal{Y}\} \subseteq \tilde{\mathcal{H}}$.

The random Fourier feature operator is useful to show the relations between the random Fourier feature map with the functional feature map defined in Proposition 21. The relationship between the generic feature map (defined for all Operator-Valued Kernel) the functional feature map (defining a shift-invariant $\mathcal{Y}$-Mercer Operator-Valued Kernel) and the random Fourier feature map is presented in Figure 4.

**Proposition 25**

For any $g \in \mathcal{H} = L^2(\hat{\mathcal{X}}, \mathrm{Pr}_{\overline{\mathrm{Haar}},\rho}; \mathcal{Y}')$, let $\theta := \frac{1}{\sqrt{D}} \bigoplus_{j=1}^D g(\omega_j)$, $\omega_j \sim \mathrm{Pr}_{\overline{\mathrm{Haar}},\rho}$ i. i. d. Then

$$ 1. \quad (\tilde{W}\theta)(x) = \tilde{\phi}(x)^*\theta \xrightarrow{a. s.} \phi_x^*g = (Wg)(x), $$

$$ 2. \quad \|\theta\|_{\tilde{\mathcal{H}}}^2 \xrightarrow{a. s.} \|g\|_{\mathcal{H}}^2, $$

We write $\tilde{\phi}(x)^*\tilde{\phi}(x) \approx K(x, z)$ when $\tilde{\phi}(x)^*\tilde{\phi}(x) \xrightarrow{a. s.} K(x, z)$ in the weak operator topology when $D$ tends to infinity. With mild abuse of notation we say that $\tilde{\phi}(x)$ is an approximate feature map of the functional feature map $\phi_x$ i. e. $\tilde{\phi}(x) \approx \phi_x$, when for all $y', y \in \mathcal{Y}$,

$$ \langle y, K(x, z)y' \rangle_y = \langle \phi_x y, \phi_z y' \rangle_{L^2(\hat{\mathcal{X}}, \mathrm{Pr}_{\overline{\mathrm{Haar}},\rho}; \mathcal{Y}')} \approx \langle \tilde{\phi}(x)y, \tilde{\phi}(x)y' \rangle_{\tilde{\mathcal{H}}} := \langle y, \tilde{K}(x, z)y' \rangle_y $$

where $\phi_x$ is defined in the sense of Proposition 21.

## 3.5 Examples of Operator Random Fourier Feature maps

We now give two examples of operator-valued random Fourier feature map. First we introduce the general form of an approximated feature map for a matrix-valued kernel on the additive group $(\mathbb{R}^d, +)$.

**Example 1 (Matrix-valued kernel on the additive group)**
In the following let $K(x, z) = K_0(x - z)$ be a $\mathcal{Y}$-Mercer matrix-valued kernel on $\mathcal{X} = \mathbb{R}^d$, invariant w. r. t. the group operation $+$. Then the function $\tilde{\phi}$ defined as follow is an Operator-valued Random Fourier Feature of $K_0$. For all $y \in \mathcal{Y}$,

$$ \tilde{\phi}(x)y = \frac{1}{\sqrt{D}} \bigoplus_{j=1}^{D} \begin{pmatrix} \cos \langle x, \omega_j \rangle_2 B(\omega_j)^* y \\ \sin \langle x, \omega_j \rangle_2 B(\omega_j)^* y \end{pmatrix}, \quad \omega_j \sim \mathbf{Pr}_{\widehat{\text{Haar}}, \rho} \text{ i. i. d..} $$

In particular we deduce the following features maps for the kernels proposed in Subsection 3.2.

* For the decomposable Gaussian kernel $K_0^{\text{dec,gauss}}(\delta) = k_0^{\text{gauss}}(\delta)\Gamma$ for all $\delta \in \mathbb{R}^d$, let $BB^* = \Gamma$. A bounded-and unbounded-ORFF map is

$$ \tilde{\phi}(x)y = \frac{1}{\sqrt{D}} \bigoplus_{j=1}^{D} \begin{pmatrix} \cos \langle x, \omega_j \rangle_2 B^* y \\ \sin \langle x, \omega_j \rangle_2 B^* y \end{pmatrix} = (\tilde{\phi}(x) \otimes B^*)y, $$

where $\omega_j \sim \mathbf{Pr}_{\mathcal{N}(0, \sigma^{-2}I_d)}$ i. i. d. and $\tilde{\phi}(x) = \frac{1}{\sqrt{D}} \bigoplus_{j=1}^{D} \begin{pmatrix} \cos \langle x, \omega_j \rangle_2 \\ \sin \langle x, \omega_j \rangle_2 \end{pmatrix}$ is a scalar RFF map (Rahimi and Recht, 2007).

* For the curl-free Gaussian kernel, $K_0^{\text{curl,gauss}} = -\nabla \nabla^\top k_0^{\text{gauss}}$ an unbounded ORFF map is

$$ \tilde{\phi}(x)y = \frac{1}{\sqrt{D}} \bigoplus_{j=1}^{D} \begin{pmatrix} g \cos \langle x, \omega_j \rangle_2 \omega_j^\top y \\ \sin \langle x, \omega_j \rangle_2 \omega_j^\top y \end{pmatrix}, \quad (17) $$

$\omega_j \sim \mathbf{Pr}_{\mathcal{N}(0, \sigma^{-2}I_d)}$ i. i. d. and a bounded ORFF map is

$$ \tilde{\phi}(x)y = \frac{1}{\sqrt{D}} \bigoplus_{j=1}^{D} \begin{pmatrix} \cos \langle x, \omega_j \rangle_2 \frac{\omega_j^\top}{\|\omega_j\|} y \\ \sin \langle x, \omega_j \rangle_2 \frac{\omega_j^\top}{\|\omega_j\|} y \end{pmatrix}, \quad \omega_j \sim \mathbf{Pr}_\rho \text{ i. i. d..} $$

where $\rho(\omega) = \frac{\sigma^2 \|\omega\|^2}{d} \mathcal{N}(0, \sigma^{-2}I_d)(\omega)$ for all $\omega \in \mathbb{R}^d$.

* For the divergence-free Gaussian kernel $K_0^{\text{div,gauss}}(x, z) = (\nabla \nabla^\top - \Delta I_d)k_0^{\text{gauss}}(x, z)$ an unbounded ORFF map is

$$ \tilde{\phi}(x)y = \frac{1}{\sqrt{D}} \bigoplus_{j=1}^{D} \begin{pmatrix} \cos \langle x, \omega_j \rangle_2 B(\omega_j)^\top y \\ \sin \langle x, \omega_j \rangle_2 B(\omega_j)^\top y \end{pmatrix} \quad (18) $$

where $\omega_j \sim \mathbf{Pr}_\rho$ i. i. d. and $B(\omega) = (\|\omega\|I_d - \omega\omega^\top)$ and $\rho = \mathcal{N}(0, \sigma^{-2}I_d)$ for all $\omega \in \mathbb{R}^d$. A bounded ORFF map is

$$ \tilde{\phi}(x)y = \frac{1}{\sqrt{D}} \bigoplus_{j=1}^{D} g \begin{pmatrix} \cos \langle x, \omega_j \rangle_2 B(\omega_j)^\top y \\ \sin \langle x, \omega_j \rangle_2 B(\omega_j)^\top y \end{pmatrix}, \quad \omega_j \sim \mathbf{Pr}_\rho \text{ i. i. d.,} $$

where $B(\omega) = (I_d - \frac{\omega\omega^\top}{\|\omega\|^2})$ and $\rho(\omega) = \frac{\sigma^2 \|\omega\|^2}{d} \mathcal{N}(0, \sigma^{-2}I_d)$ for all $\omega \in \mathbb{R}^d$.
The second example extends scalar-valued Random Fourier Features on the skewed multiplicative group –described in appendix A.3 and Subsection 3.2.2– to the operator-valued case.

**Example 2** (Matrix-valued kernel on the skewed multiplicative group)

In the following, $K(x, z) = K_{1-c}(x \odot z^{-1})$ is a $\mathcal{Y}$-Mercer matrix-valued kernel on $\mathcal{X} = (-c; +\infty)^d$ invariant w. r. t. the group operation$^5 \odot$. Then the function $\tilde{\phi}$ defined as follow is an Operator-valued Random Fourier Feature of $K_{1-c}$. For all $y \in \mathcal{Y}$,

$$ \tilde{\phi}(x)y = \frac{1}{\sqrt{D}} \bigoplus_{j=1}^{D} \left( \begin{array}{l} \cos \langle \log(x+c), \omega_j \rangle_2 B(\omega_j)^* y \\ \sin \langle \log(x+c), \omega_j \rangle_2 B(\omega_j)^* y \end{array} \right), \quad \omega_j \sim \text{Pr}_{\widehat{\text{Haar}}, \rho} \quad \text{i. i. d..} $$

## 3.6 Regularization property

We have shown so far that it is always possible to construct a feature map that allows to approximate a shift-invariant $\mathcal{Y}$-Mercer kernel. However we could also propose a construction of such map by studying the regularization induced with respect to the Fourier Transform of a target function $f \in \mathcal{H}_K$. In other words, what is the norm in $L^2(\hat{\mathcal{X}}, \widehat{\text{Haar}}; \mathcal{Y}')$ induced by $\|\cdot\|_K$?

**Proposition 26**

Let $K$ be a shift-invariant $\mathcal{Y}$-Mercer Kernel such that for all $y, y'$ in $\mathcal{Y}$, $\langle y', K_e(\cdot)y \rangle_y \in L^1(\mathcal{X}, \widehat{\text{Haar}})$. Then for all $f \in \mathcal{H}_K$

$$ \|f\|_K^2 = \int_{\hat{\mathcal{X}}} \frac{\langle \mathcal{F}[f](\omega), A(\omega)^{\dagger} \mathcal{F}[f](\omega) \rangle_y}{\rho(\omega)} d\widehat{\text{Haar}}(\omega). \quad (19) $$

where $\langle y', A(\omega)y \rangle_y \rho(\omega) := \mathcal{F}[\langle y', K_e(\cdot)y \rangle](\omega)$.

Note that if $K(x, z) = k(x, z)$ is a scalar kernel then for all $\omega$ in $\hat{\mathcal{X}}$, $A(\omega) = 1$. Therefore we recover the well known result for kernels that is for any $f \in \mathcal{H}_K$ we have $\|f\|_K = \int_{\hat{\mathcal{X}}} \mathcal{F}[k_e](\omega)^{-1} \mathcal{F}[f](\omega)^2 d\widehat{\text{Haar}}(\omega)$ (Yang et al., 2012; Vert; Smola et al., 1998). Eventually from this last equation we also recover Proposition 9 for decomposable kernels. If $A(\omega) = \Gamma \in \mathcal{L}_+(\mathbb{R}^p)$, $\|f\|_K = \sum_{i,j=1}^p (\Gamma^\dagger)_{ij} \langle f_i, f_j \rangle_k$. We also note that the regularization property in $\mathcal{H}_K$ does not depends (as expected) on the decomposition of $A(\omega)$ into $B(\omega)B(\omega)^*$. Therefore the decomposition should be chosen such that it optimizes the computation cost. For instance if $A(\omega) \in \mathcal{L}(\mathbb{R}^p)$ has rank $r$, one could find an operator $B(\omega) \in \mathcal{L}(\mathbb{R}^p, \mathbb{R}^r)$ such that $A(\omega) = B(\omega)B(\omega)^*$. Moreover, in light of Equation 19 the regularization property of the kernel with respect to the Fourier Transform, it is also possible to define an approximate feature map of an Operator-Valued Kernel from its regularization properties in the VV-RKHS as proposed in Algorithm 2.

<footnote n="5">5. The group operation $\odot$ is defined in Subsection 3.2.2.</footnote>
**Algorithm 2:** Construction of ORFF

**Input :**

* The pairing $(x, \omega)$ of the LCA group $(\mathcal{X}, \star)$.

* A probability measure $\mathrm{Pr}_{\widehat{\mathrm{Haar}}, \rho}$ with density $\rho$ w.r.t. the haar measure $\widehat{\mathrm{Haar}}$ on $\widehat{\mathcal{X}}$.

* An operator-valued function $B: \widehat{\mathcal{X}} \to \mathcal{L}(\mathcal{Y}, \mathcal{Y}')$ such that for all $y, y' \in \mathcal{Y}$, $\langle y', B(\cdot)B(\cdot)^*y \rangle \in L^1(\widehat{\mathcal{X}}, \mathrm{Pr}_{\widehat{\mathrm{Haar}}, \rho})$.

* $D$ the number of features.

**Output:** A random feature $\tilde{\phi}(x)$ such that $\tilde{\phi}(x)^*\tilde{\phi}(z) \approx K(x, z)$.

1. Draw $D$ random vectors $(\omega_j)_{j=1}^D$ i. i. d. from the probability law $\mathrm{Pr}_{\widehat{\mathrm{Haar}}, \rho}$;

2. return $\begin{cases} \tilde{\phi}(x) \in \mathcal{L}(\mathcal{Y}, \tilde{\mathcal{H}}) & : y \mapsto \frac{1}{\sqrt{D}} \bigoplus_{j=1}^D (x, \omega_j) B(\omega_j)^* y; \\ \tilde{\phi}(x)^* \in \mathcal{L}(\tilde{\mathcal{H}}, \mathcal{Y}) & : \theta \mapsto \frac{1}{\sqrt{D}} \sum_{j=1}^D (x, \omega_j) B(\omega_j) \theta_j \end{cases}$

4. Main contribution: convergence with high probability of the ORFF estimator

We are now interested in a non-asymptotic analysis of the ORFF approximation of shift-invariant $\mathcal{Y}$-Mercer kernels on LCA group $\mathcal{X}$ endowed with the operation group $\star$ where $\mathcal{X}$ is a Banach space (The more general case where $\mathcal{X}$ is a Polish space is discussed in the appendix appendix B.2). For a given $D$, we study how close is the approximation $\tilde{K}(x, z) = \tilde{\phi}(x)^*\tilde{\phi}(z)$ to the target kernel $K(x, z)$ for any $x, z$ in $\mathcal{X}$.

If $A \in \mathcal{L}_+(\mathcal{Y})$ we denote $\|A\|_{\mathcal{Y},\mathcal{Y}}$ its operator norm (the induced norm). For $x$ and $z$ in some non-empty compact $\mathcal{C} \subset \mathbb{R}^d$, we consider: $F(x \star z^{-1}) = \tilde{K}(x, z) - K(x, z)$ and study how the uniform norm $\|\tilde{K} - K\|_{\mathcal{C} \times \mathcal{C}} := \sup_{(x,z) \in \mathcal{C} \times \mathcal{C}} \|\tilde{K}(x, z) - K(x, z)\|_{\mathcal{Y},\mathcal{Y}}$ behaves according to $D$. All along this document we denote $\delta = x \star z^{-1}$ for all $x$ and $z \in \mathcal{X}$. Figure 5 empirically shows convergence of three different OVK approximations for $x, z$ sampled from the compact $[-1, 1]^4$ and using an increasing number of sample points $D$. The log-log plot shows that all three kernels have the same convergence rate, up to a multiplicative factor.

A typical application is the study of the deviation of the empirical mean of independent identically distributed random variables to their expectation. This means that given an error $\epsilon$ between the kernel approximation $\tilde{K}$ and the true kernel $K$, if we are given enough samples to construct $\tilde{K}$, the probability of measuring an error greater than $\epsilon$ is essentially zero (it drops at an exponential rate with respect to the number of samples $D$). To measure the error between the kernel approximation and the true kernel at a given point many metrics are possible. e. g. any matrix norm such as the Hilbert-Schmidt norm, trace norm, the operator norm or Schatten norms. In this work we focus on measuring the error in terms of
operator norm. For all $x, z \in \mathcal{X}$ we look for a bound on

$$ \mathbf{Pr}_{\rho} \left\{ (\omega_j)_{j=1}^D \left| \left\| \tilde{K}(x, z) - K(x, z) \right\|_{\mathcal{Y}, \mathcal{Y}} \ge \epsilon \right. \right\} \\ = \mathbf{Pr}_{\rho} \left\{ (\omega_j)_{j=1}^D \left| \sup_{0 \ne y \in \mathcal{Y}} \frac{\left\| (\tilde{K}(x, z) - K(x, z))y \right\|_{\mathcal{Y}}}{\|y\|_{\mathcal{Y}}} \ge \epsilon \right. \right\} $$

In other words, given any vector $y \in \mathcal{Y}$ we study how the residual operator $\tilde{K} - K$ is able to send $y$ to zero. We believe that this way of measuring the “error” to be more intuitive. Moreover, on contrary to an error measure with the Hilbert-Schmidt norm, the operator norm error does not grows linearly with the dimension of the output space as the Hilbert-Schmidt norm does. On the other hand the Hilbert-schmidt norm makes the studied random variables Hilbert space valued, for which it is much easier to derive concentration inequalities (Smale and Zhou, 2007; Pinelis, 1994; Naor, 2012). Note that in the scalar case ($A(\omega) = 1$) the Hilbert-Schmidt norm error and the operator norm are the same and measure the deviation between $\tilde{K}$ and $K$ as the absolute value of their difference.

A raw concentration inequality of the kernel estimator gives the error on one point. If one is interesting in bounding the maximum error over $N$ points, applying a union bound on all the point would yield a bound that grows linearly with $N$. This would suggest that when the number of points increase, even if all of them are concentrated in a small subset of $\mathcal{X}$, we should draw increasingly more features to have an error below $\epsilon$ with high probability. However if we restrict ourselves to study the error on a compact subset of $\mathcal{X}$ (and in practice data points lies often in a closed bounded subset of $\mathbb{R}^d$), we can cover this compact subset by a finite number of closed balls and apply the concentration inequality and the union bound only on the center of each ball. Then if the function $\|\tilde{K}_e - K_e\|$ is smooth enough on each ball (i. e. Lipschitz) we can guarantee with high probability that the error between the centers of the balls will not be too high. Eventually we obtain a bound in the worst case scenario on all the points in a subset $\mathcal{C}$ of $\mathcal{X}$. This bound depends on the covering number $\mathcal{N}(\mathcal{C}, r)$ of $\mathcal{X}$ with ball of radius $r$. When $\mathcal{X}$ is a Banach space, the covering number is proportional to the diameter of the diameter of $\mathcal{C} \subseteq \mathcal{X}$.

Prior to the presentation of general results, we briefly recall the uniform convergence of RFF approximation for a scalar shift invariant kernel on the additive LCA group $\mathbb{R}^d$ and introduce a direct corollary about decomposable shift-invariant OVK on the LCA group $(\mathbb{R}^d, +)$.

## 4.1 Random Fourier Features in the scalar case and decomposable OVK

Rahimi and Recht (2007) proved the uniform convergence of Random Fourier Feature (RFF) approximation for a scalar shift-invariant kernel on the LCA group $\mathbb{R}^d$ endowed with the group operation $\star = +$. In the case of the shift-invariant decomposable OVK, an upper bound on the error can be obtained as a direct consequence of the result in the scalar case obtained by Rahimi and Recht (2007) and other authors (Sutherland and Schneider, 2015; Sriperumbudur and Szabo, 2015).
<figure><img src="image_6.png" /><figcaption>Figure 5: Error reconstructing the target operator-valued kernel $K$ with ORFF approximation $\tilde{K}$ for the decomposable, curl-free and divergence-free kernel.</figcaption></figure>

**Theorem 27** (Uniform error bound for RFF, Rahimi and Recht (2007))

Let $C$ be a compact of subset of $\mathbb{R}^d$ of diameter $|C|$. Let $k$ be a shift invariant kernel, differentiable with a bounded second derivative and $\mathbf{Pr}_\rho$ its normalized Inverse Fourier Transform such that it defines a probability measure. Let $\tilde{k} = \sum_{j=1}^D \cos \langle \cdot, \omega_j \rangle \approx k(x, z)$ and $\sigma^2 = \mathbf{E}_\rho \|\omega\|_2^2$. Then we have

$$ \mathbf{Pr}_\rho \left\{ (\omega_j)_{j=1}^D \left| \left\| \tilde{k} - k \right\|_{C \times C} \ge \epsilon \right. \right\} \le 2^8 \left( \frac{\sigma |C|}{\epsilon} \right)^2 \exp \left( - \frac{\epsilon^2 D}{4(d+2)} \right) $$

From Theorem 27, we can deduce the following corollary about the uniform convergence of the ORFF approximation of the decomposable kernel. We recall that for a given pair $x, z$ in $C$, $\tilde{K}(x, z) = \tilde{\phi}(x)^* \tilde{\phi}(z) = \Gamma \tilde{k}(x, z)$ and $K_0(x - z) = \Gamma \mathbf{E}_{\widehat{\text{Haar}}, \rho}[\tilde{k}(x, z)]$.

**Corollary 28** (Uniform error bound for decomposable ORFF)

Let $C$ be a compact of subset of $\mathbb{R}^d$ of diameter $|C|$. Let $K$ be a decomposable kernel built from a positive operator self-adjoint $\Gamma$, and $k$ a shift invariant kernel with bounded second derivative such that $\tilde{K} = \sum_{j=1}^D \cos \langle \cdot, \omega_j \rangle \Gamma \approx K$ and $\sigma^2 = \mathbf{E}_\rho \|\omega\|_2^2$. Then

$$ \mathbf{Pr}_\rho \left\{ (\omega_j)_{j=1}^D \left| \left\| \tilde{K} - K \right\|_{C \times C} \ge \epsilon \right. \right\} \le 2^8 \left( \frac{\sigma \|\Gamma\|_{y,y} |C|}{\epsilon} \right)^2 \exp \left( - \frac{\epsilon^2 D}{4 \|\Gamma\|_2^2 (d+2)} \right) $$
Note that a similar corollary could have been obtained for the recent result of Sutherland and Schneider (2015) who refined the bound proposed by Rahimi and Recht by using a Bernstein concentration inequality instead of the Hoeffding inequality. More recently Sriperumbudur and Szabo (2015) showed an optimal bound for Random Fourier Feature. The improvement of Sriperumbudur and Szabo (2015) is mainly in the constant factors where the bound does not depend linearly on the diameter $|C|$ of $C$ but exhibit a logarithmic dependency $\log(|C|)$, hence requiring significantly less random features to reach a desired uniform error with high probability. Moreover, Sutherland and Schneider (2015) also considered a bound on the expected max error $\mathbf{E}_{\text{Haar},p} \left\| \tilde{K} - K \right\|_{\infty}$, which is obtained using Dudley's entropy integral (Dudley, 1967; Boucheron et al., 2013) as a bound on the supremum of an empirical process by the covering number of the indexing set. This useful theorem is also part of the proof of Sriperumbudur and Szabo (2015).

## 4.2 Uniform convergence of ORFF approximation on LCA groups

Before introducing the new theorem, we give the definition of the Orlicz norm which gives a proxy-bound on the norm of subexponential random variables.

**Definition 29** (Orlicz norm (Van Der Vaart and Wellner, 1996))

Let $\psi: \mathbb{R}_+ \to \mathbb{R}_+$ be a non-decreasing convex function with $\psi(0) = 0$. For a random variable $X$ on a measured space $(\Omega, \mathcal{T}(\Omega), \mu)$, the quantity $\|X\|_{\psi} = \inf\{C > 0 \mid \mathbf{E}_{\mu}[\psi(|X|/C)] \le 1\}$. is called the Orlicz norm of $X$.

Here, the function $\psi$ is chosen as $\psi(u) = \psi_{\alpha}(u)$ where $\psi_{\alpha}(u) := e^{u^{\alpha}} - 1$. When $\alpha = 1$, a random variable with finite Orlicz norm is called a *subexponential variable* because its tails decrease at an exponential rate. Let $X$ be a self-adjoint random operator. Given a scalar-valued measure $\mu$, we call *variance* of an operator $X$ the quantity $\mathrm{Var}_{\mu}[X] = \mathbf{E}_{\mu}[X - \mathbf{E}_{\mu}[X]]^2$. Among the possible concentration inequalities adapted to random operators (Tropp et al., 2015; Minsker, 2011; Ledoux and Talagrand, 2013; Pinelis, 1994; Koltchinskii et al., 2013), we focus on the results of Tropp et al. (2015); Minsker (2011), for their robustness to high or potentially infinite dimension of the output space $\mathcal{Y}$. To guarantee a good scaling with the dimension of $\mathcal{Y}$ we introduce the notion of intrinsic dimension (or effective rank) of an operator.

**Definition 30**

Let $A$ be a trace class operator acting on a Hilbert space $\mathcal{Y}$. We call intrinsic dimension the quantity: $\mathrm{IntDim}(A) = \|A\|_{\mathcal{Y},\mathcal{Y}}^{-1} \mathrm{Tr}[A]$.

Indeed the bound proposed in our first publication at ACML (Brault et al., 2016) based on Koltchinskii et al. (2013) depends on $p$ while the present bound depends on the intrinsic dimension of the variance of $A(\omega)$ which is always smaller than $p$ when the operator $A(\omega)$ is Hilbert-Schmidt ($p \le \infty$).

**Corollary 31**

Let $K: \mathcal{X} \times \mathcal{X} \to \mathcal{L}(\mathcal{Y})$ be a shift-invariant $\mathcal{Y}$-Mercer kernel, where $\mathcal{Y}$ is a finite dimensional Hilbert space of dimension $p$ and $\mathcal{X}$ a finite dimensional Banach space of dimension $d$.
Moreover, let $\mathcal{C}$ be a closed ball of $\mathcal{X}$ centred at the origin of diameter $|\mathcal{C}|$, $A: \hat{\mathcal{X}} \to \mathcal{L}(\mathcal{Y})$ and $\mathbf{Pr}_{\widehat{\text{Haar}},\rho}$ a pair such that

$$ \tilde{K}_e = \sum_{j=1}^{D} \cos(\cdot, \omega_j) A(\omega_j) \approx K_e, \quad \omega_j \sim \mathbf{Pr}_{\widehat{\text{Haar}},\rho} \quad \text{i. i. d...} $$

Let $\mathcal{D}_C = \mathcal{C} \star \mathcal{C}^{-1}$ and $V(\delta) \asymp \mathbf{Var}_{\widehat{\text{Haar}},\rho} \tilde{K}_e(\delta)$ for all $\delta \in \mathcal{D}_C$ and $H_\omega$ be the Lipschitz constant of the function $h: x \mapsto (x, \omega)$. If the three following constants exist

$$ m \geq \int_{\hat{\mathcal{X}}} H_{\omega} \|A(\omega)\|_{\mathcal{Y},\mathcal{Y}} d\mathbf{Pr}_{\widehat{\text{Haar}},\rho}(\omega) < \infty $$

and

$$ u \geq 4 \left( \left\| \|A(\omega)\|_{\mathcal{Y},\mathcal{Y}} \right\|_{\psi_1} + \sup_{\delta \in \mathcal{D}_C} \|K_e(\delta)\|_{\mathcal{Y},\mathcal{Y}} \right) < \infty $$

and

$$ v \geq \sup_{\delta \in \mathcal{D}_C} D \|V(\delta)\|_{\mathcal{Y},\mathcal{Y}} < \infty. $$

Define $p_{\text{int}} \geq \sup_{\delta \in \mathcal{D}_C} \text{IntDim}(V(\delta))$, then for all $0 < \epsilon \leq m|\mathcal{C}|$,

$$ \mathbf{Pr}_{\widehat{\text{Haar}},\rho} \left\{ (\omega_j)_{j=1}^D \left| \left\| \tilde{K} - K \right\|_{\mathcal{C} \times \mathcal{C}} \geq \epsilon \right. \right\} \\ \leq 8\sqrt{2} \left( \frac{m|\mathcal{C}|}{\epsilon} \right) (p_{\text{int}} r_{v/D}(\epsilon))^{\frac{1}{d+1}} \begin{cases} \exp \left( -D \frac{\epsilon^2}{8v(d+1)(1+\frac{1}{p})} \right), & \epsilon \leq \frac{v}{u} \frac{1+1/p}{K(v,p)} \\ \exp \left( -D \frac{\epsilon}{8u(d+1)K(v,p)} \right), & \text{otherwise,} \end{cases} $$

$$ \text{where } K(v,p) = \log(16\sqrt{2}p) + \log\left(\frac{u^2}{v}\right) \text{ and } r_{v/D}(\epsilon) = 1 + \frac{3}{\epsilon^2 \log^2(1+D\epsilon/v)}. $$

We give a comprehensive full proof of the theorem in appendix B.2. It follows the usual scheme derived in Rahimi and Recht (2007) and Sutherland and Schneider (2015) and involves Bernstein concentration inequality for unbounded symmetric matrices (Theorem 45).

## 4.3 Dealing with infinite dimensional operators

We studied the concentration of ORFFs under the assumption that $\mathcal{Y}$ is finite dimensional. Indeed a $d$ term characterizing the dimension of the input space $\mathcal{X}$ appears in the bound proposed in Corollary 31, and when $d$ tends to infinity, the exponential part goes to zero so that the probability is bounded by a constant greater than one. Unfortunately, considering unbounded random operators Minsker (2011) doesn't give any tighter solution.

In our first bound presented at ACML, we presented a bound based on a matrix concentration inequality for unbounded random variable. Compared to this previous bound, Corollary 31 does not depend on the dimensionality $p$ of the output space $\mathcal{Y}$ but on the intrinsic dimension of the operator $A(\omega)$. However to remove the dependency in $p$ in the exponential part, we must turn our attention to operator concentration inequalities for
bounded random variable. To the best of our knowledge we are not aware of concentration inequalities working for “unbounded” operator- valued random variables. Following the same proof than Corollary 31 we obtain Corollary 32.

**Corollary 32**

Let $K: \mathcal{X} \times \mathcal{X} \to \mathcal{L}(\mathcal{Y})$ be a shift-invariant $\mathcal{Y}$-Mercer kernel, where $\mathcal{Y}$ is a Hilbert space and $\mathcal{X}$ a finite dimensional Banach space of dimension $D$. Moreover, let $C$ be a closed ball of $\mathcal{X}$ centered at the origin of diameter $|C|$, subset of $\mathcal{X}$, $A: \hat{\mathcal{X}} \to \mathcal{L}(\mathcal{Y})$ and $\mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho}$ a pair such that

$$ \tilde{K}_e = \sum_{j=1}^{D} \cos(\cdot, \omega_j) A(\omega_j) \approx K_e, \quad \omega_j \sim \mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho} \quad \text{i. i. d..} $$

where $A(\omega_j)$ is a Hilbert-Schmidt operator for all $j \in \mathbb{N}_D^*$. Let $\mathcal{D}_C = C \star C^{-1}$ and $V(\delta) \asymp \mathbf{Var}_{\widehat{\mathcal{Haar}},\rho} \tilde{K}_e(\delta)$ for all $\delta \in \mathcal{D}_C$ and $H_\omega$ be the Lipschitz constant of the function $h: x \mapsto (x, \omega)$. If the three following constants exists

$$ m \geq \int_{\hat{\mathcal{X}}} H_{\omega} \|A(\omega)\|_{\mathcal{Y},\mathcal{Y}} d\mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho}(\omega) < \infty $$

and

$$ u \geq \sup_{\omega \in \hat{\mathcal{X}}} \|A(\omega)\|_{\mathcal{Y},\mathcal{Y}} + \sup_{\delta \in \mathcal{D}_C} \|K_e(\delta)\|_{\mathcal{Y},\mathcal{Y}} < \infty $$

and

$$ v \geq \sup_{\delta \in \mathcal{D}_C} D \|V(\delta)\|_{\mathcal{Y},\mathcal{Y}} < \infty. $$

define $p_{int} \geq \sup_{\delta \in \mathcal{D}_C} \mathrm{IntDim}(V(\delta))$ then for all $\sqrt{\frac{v}{D}} + \frac{u}{3D} < \epsilon < m|C|,$

$$ \mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho} \left\{ \left( \omega_j \right)_{j=1}^{D} \left| \sup_{\delta \in \mathcal{D}_C} \|F(\delta)\|_{\mathcal{Y},\mathcal{Y}} \geq \epsilon \right. \right\} \leq 8\sqrt{2} \left( \frac{m|C|}{\epsilon} \right) p_{int}^{\frac{1}{d+1}} \exp(-D\psi_{v,d,u}(\epsilon)) $$

$$ \text{where } \psi_{v,d,u}(\epsilon) = \frac{\epsilon^2}{2(d+1)(v+u\epsilon/3)}. $$

Again a full comprehensive proof is given in appendix B.2 of the appendix. Notice that in this result, The dimension $p = \mathrm{dim} \mathcal{Y}$ does not appear. Only the intrinsic dimension of the variance of the estimator. Moreover when $d$ is large, the term $p_{int}^{\frac{1}{d+1}}$ goes to one, so that the impact of the intrinsic dimension on the bound vanish when the dimension of the input space is large. subsectionVariance of the ORFF approximation We now provide a bound on the norm of the variance of $\tilde{K}$, required to apply Corollaries 31 and 32. This is an extension of the proof of Sutherland and Schneider (2015) to the operator-valued case, and we recover their results in the scalar case when $A(\omega) = 1$. An illustration of the bound is provided in Figure 6 for the decomposable and the curl-free OVK.

**Proposition 33 (Bounding the variance of $\tilde{K}_e$)**

Let $K$ be a shift invariant $\mathcal{Y}$-Mercer kernel on a second countable LCA topological space $\mathcal{X}$. Let $A: \hat{\mathcal{X}} \to \mathcal{L}(\mathcal{Y})$ and $\mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho}$ a pair such that $\tilde{K}_e = \sum_{j=1}^{D} \cos(\cdot, \omega_j)A(\omega_j) \approx K_e$, $\omega_j \sim \mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho}$ i. i. d. Then,

$$ \mathbf{Var}_{\widehat{\mathcal{Haar}},\rho} \left[ \tilde{K}_e(\delta) \right] \asymp \frac{1}{2D} \left( (K_e(2\delta) + K_e(e)) \mathbf{E}_{\widehat{\mathcal{Haar}},\rho} [A(\omega)] - 2K_e(\delta)^2 + \mathbf{Var}_{\widehat{\mathcal{Haar}},\rho} [A(\omega)] \right) $$
<figure><img src="image_7.png" /><figcaption>Figure 6: Comparison between an empirical bound on the norm of the variance of the decomposable (left) and curl-free (right) ORFF obtained and the theoretical bound proposed in Proposition 33 versus D.</figcaption></figure>

## 4.4 Application on decomposable, curl-free and divergence-free OVK

First, the two following examples discuss the form of $H_\omega$ for the additive group and the skewed-multiplicative group. Here we view $\mathcal{X} = \mathbb{R}^d$ as a Banach space endowed with the Euclidean norm. Thus the Lipschitz constant $H_\omega$ is bounded by the supremum of the norm of the gradient of $h_\omega$.

### Example 3 (Additive group)

On the additive group, $h_\omega(\delta) = \langle \omega, \delta \rangle$. Hence $H_\omega = \|\omega\|_2$.

### Example 4 (Skewed-multiplicative group)

On the skewed multiplicative group, $h_\omega(\delta) = \langle \omega, \log(\delta + c) \rangle$. Therefore $\sup_{\delta \in \mathcal{C}} \|\nabla h_\omega(\delta)\|_2 = \sup_{\delta \in \mathcal{C}} \|\omega/(\delta + c)\|_2$. Eventually $\mathcal{C}$ is compact subset of $\mathcal{X}$ and finite dimensional thus $\mathcal{C}$ is closed and bounded. Thus $H_\omega = \|\omega\|_2/(\min_{\delta \in \mathcal{C}} \|\delta\|_2 + c)$.

Now we compute upper bounds on the norm of the variance and Orlicz norm of the three ORFFs we took as examples.

### 4.4.1 DECOMPOSABLE KERNEL

notice that in the case of the Gaussian decomposable kernel, i. e. $A(\omega) = A$, $e = 0$, $K_0(\delta) = Ak_0(\delta)$, $k_0(\delta) \ge 0$ and $k_0(\delta) = 1$, then we have

$$D \| \mathbf{Var}_{\mu} [\tilde{K}_0(\delta)] \|_{y,y} \le (1 + k_0(2\delta)) \|A\|_{y,y} / 2 + k_0(\delta)^2.$$

### 4.4.2 CURL-FREE AND DIVERGENCE-FREE KERNELS:

recall that in this case $p = d$. For the (Gaussian) curl-free kernel, $A(\omega) = \omega\omega^*$ where $\omega \in \mathbb{R}^d \sim \mathcal{N}(0, \sigma^{-2}I_d)$ thus $\mathbf{E}_\mu[A(\omega)] = I_d/\sigma^2$ and $\mathbf{Var}_\mu[A(\omega)] = (d+1)I_d/\sigma^4$. Hence,

$$D \| \mathbf{Var}_{\mu} [\tilde{K}_0(\delta)] \|_{y,y} \le \frac{1}{2} \left\| \frac{1}{\sigma^2} K_0(2\delta) - 2K_0(\delta)^2 \right\|_{y,y} + \frac{(d+1)}{\sigma^4}.$$
This bound is illustrated by Figure 5 B, for a given datapoint. Eventually for the Gaussian divergence-free kernel, $A(\omega) = I\|\omega\|_2^2 - \omega\omega^*$, thus $E_\mu[A(\omega)] = I_d(d-1)/\sigma^2$ and $Var_\mu[A(\omega)] = d(4d-3)I_d/\sigma^4$. Hence,

$$D \left\| \mathrm{Var}_\mu \left[ \tilde{K}_0(\delta) \right] \right\|_{y,y} \le \frac{1}{2} \left\| \frac{(d-1)}{\sigma^2} K_0(2\delta) - 2K_0(\delta)^2 \right\|_{y,y} + \frac{d(4d-3)}{\sigma^4}.$$

To conclude, we ensure that the random variable $\|A(\omega)\|_{y,y}$ has a finite Orlicz norm with $\psi = \psi_1$ in these three cases.

### 4.4.3 COMPUTING THE ORLICZ NORM

for a random variable with strictly monotonic moment generating function (MGF), one can characterize its inverse $\psi_1$ Orlicz norm by taking the functional inverse of the MGF evaluated at 2 (see Lemma 43 of the appendix). In other words $\|X\|_{\psi_1}^{-1} = \mathrm{MGF}(x)_X^{-1}(2)$. For the Gaussian curl-free and divergence-free kernel,

$$\left\| A^{\mathrm{div}}(\omega) \right\|_{y,y} = \left\| A^{\mathrm{curl}}(\omega) \right\|_{y,y} = \|\omega\|_2^2,$$

where $\omega \sim \mathcal{N}(0, I_d/\sigma^2)$, hence $\|A(\omega)\|_2 \sim \Gamma(p/2, 2/\sigma^2)$. The MGF of this gamma distribution is $\mathrm{MGF}(x)(t) = (1 - 2t/\sigma^2)^{-(p/2)}$. Eventually

$$\left\| A^{\mathrm{div}}(\omega) \right\|_{y,y}^{-1}_{\psi_1} = \left\| A^{\mathrm{curl}}(\omega) \right\|_{y,y}^{-1}_{\psi_1} = \frac{\sigma^2}{2} \left( 1 - 4^{-\frac{1}{p}} \right).$$

## 5. Learning with ORFF

Before focusing on learning function with an ORFF model, we briefly review the context of supervised learning in VV-RKHS. model.

### 5.1 Supervised learning within VV-RKHS

Let $s = (x_i, y_i)_{i=1}^N \in (\mathcal{X} \times \mathcal{Y})^N$ be a sequence of training samples. Given a local loss function $L: \mathcal{X} \times \mathcal{F} \times \mathcal{Y} \to \mathbb{R}$ such that $L$ is proper, convex and lower semi-continuous in $\mathcal{F}$, we are interested in finding a vector-valued function $f_s: \mathcal{X} \to \mathcal{Y}$, that lives in a VV-RKHS and minimize a tradeoff between a data fitting term $L$ and a regularization term to prevent from overfitting. Namely finding $f_s \in \mathcal{H}_K$ such that

$$f_s = \arg \min_{f \in \mathcal{H}_K} \frac{1}{N} \sum_{i=1}^{N} L(x_i, f, y_i) + \frac{\lambda}{2} \|f\|_K^2 \quad (20)$$

where $\lambda \in \mathbb{R}_+$ is a (Tychonov) regularization hyperparameter. We call the quantity

$$\mathcal{R}_{\lambda}(f, s) = \frac{1}{N} \sum_{i=1}^{N} L(x_i, f, y_i) + \frac{\lambda}{2} \|f\|_{K}^{2}, \quad \forall f \in \mathcal{H}_{K}, \forall s \in (\mathcal{X} \times \mathcal{Y})^{N}.$$

the (Tychonov) regularized risk of the model $f \in \mathcal{H}_K$ according the local loss $L$. A common choice for $L$ is the squared error loss $L: (x, f, y) \mapsto \|f(x) - y\|_y^2$ which yields the vector-valued ridge regression problem.
### 5.1.1 REPRESENTER THEOREM AND FEATURE EQUIVALENCE

Regression in Vector Valued Reproducing Kernel Hilbert Space has been well studied (Álvarez et al., 2012; Argyriou et al., 2009; Minh et al., 2013a, 2016; Sangnier et al., 2016; Kadri et al., 2015; Micchelli and Pontil, 2005; Brouard et al., 2016a), and a cornerstone of learning in VV-RKHS is the representer theorem<sup>6</sup>, which allows to replace the search of a minimizer in a infinite dimensional VV-RKHS by a finite number of parameters $(u_i)_{i=1}^N, u_i \in \mathcal{Y}$.

In the following we suppose we are given a cost function $c: \mathcal{Y} \times \mathcal{Y} \to \mathbb{R}$, such that $c(f(x), y)$ returns the error of the prediction $f(x)$ w.r.t. the ground truth $y$. A loss function of a model $f$ with respect to an example $(x, y) \in \mathcal{X} \times \mathcal{Y}$ can be naturally defined from a cost function as $L(x, f, y) = c(f(x), y)$. Conceptually the function $c$ evaluates the quality of the prediction versus its ground truth $y \in \mathcal{Y}$ while the loss function $L$ evaluates the quality of the model $f$ at a training point $(x, y) \in \mathcal{X} \times \mathcal{Y}$.

**Theorem 34 (Representer theorem)**

Let $K$ be a $\mathcal{Y}$-Mercer Operator-Valued Kernel and $\mathcal{H}_K$ its corresponding $\mathcal{Y}$-Reproducing Kernel Hilbert space. Let $c: \mathcal{Y} \times \mathcal{Y} \to \mathbb{R}$ be a cost function such that $L(x, f, y) = c(Vf(x), y)$ is a proper convex lower semi-continuous function in $f$ for all $x \in \mathcal{X}$ and all $y \in \mathcal{Y}$. Eventually let $\lambda \in \mathbb{R}_{>0}$ be the Tychonov regularization hyperparameters The solution $f_s \in \mathcal{H}_K$ of the regularized optimization problem

$$f_s = \arg\min_{f \in \mathcal{H}_K} \frac{1}{N} \sum_{i=1}^{N} c(f(x_i), y_i) + \frac{\lambda}{2} \|f\|_{K}^{2} \quad (21)$$

has the form $f_s = \sum_{j=1}^{N} K(\cdot, x_j)u_{s,j}$ where $u_{s,j} \in \mathcal{Y}$ and

$$u_s = \arg\min_{u \in \bigoplus_{i=1}^{N} \mathcal{Y}} \frac{1}{N} \sum_{i=1}^{N} c \left( \sum_{j=1}^{N} K(x_i, x_j)u_j, y_i \right) + \frac{\lambda}{2} \sum_{j=1}^{N} u_j^* K(x_i, x_j)u_j. \quad (22)$$

The first representer theorem was introduced by Wahba (1990) in the case where $\mathcal{Y} = \mathbb{R}$. The extension to an arbitrary Hilbert space $\mathcal{Y}$ has been proved by many authors in different forms (Brouard et al., 2011; Kadri et al., 2015; Micchelli and Pontil, 2005). The idea behind the representer theorem is that even though we minimize over the whole space $\mathcal{H}_K$, when $\lambda > 0$, the solution of Equation 21 falls inevitably into the set $\mathcal{H}_{K,s} = \{\sum_{j=1}^{N} K_{x_j} u_j \mid \forall (u_i)_{i=1}^N \in \mathcal{Y}^N\}$. Therefore the result can be expressed as a finite linear combination of basis functions of the form $K(\cdot, x_k)$. Notice that we can perform the kernel expansion of $f_s = \sum_{j=1}^{N} K(\cdot, x_j)u_{s,j}$ even though $\lambda = 0$. However $f_s$ is no longer the solution of Equation 21 over the whole space $\mathcal{H}_K$ but a projection on the subspace $\mathcal{H}_{K,s}$. The representer theorem show that minimizing a functional in a VV-RKHS yields a solution which depends on all the points in the training set. Assuming that for all $x_i$ and $x \in \mathcal{X}$ and for all $u_i \in \mathcal{Y}$ it takes time $O(P)$ to compute $K(x_i, x)u_i$, making a prediction using the representer theorem takes $O(NP)$. Obviously If $\mathcal{Y} = \mathbb{R}^p$, Then $P = O(p^2)$ thus making a prediction cost $O(Np^2)$ operations.

<footnote n="6">6. Sometimes referred to as minimal norm interpolation theorem.</footnote>
## 5.2 Learning with Operator Random Fourier Feature maps

Instead of learning a model $f$ that depends on all the points of the training set, we would like to learn a parametric model of the form $\tilde{f}(x) = \tilde{\phi}(x)^*\theta$, where $\theta$ lives in some space $\tilde{\mathcal{H}}$. We are interested in finding a parameter vector $\theta_s$ such that

$$ \theta_s = \arg \min_{\theta \in \tilde{\mathcal{H}}} \mathfrak{R}_{\lambda}(\theta, s) = \arg \min_{\theta \in \tilde{\mathcal{H}}} \frac{1}{N} \sum_{i=1}^{N} c(\tilde{\phi}(x_i)^* \theta, y_i) + \frac{\lambda}{2} \|\theta\|_{\tilde{\mathcal{H}}}^2 \quad (23) $$

The following theorem states that when $\lambda > 0$ then learning with a feature map is equivalent to learn with a kernel. Moreover if $f_s \in \mathcal{H}_K$ is a solution of Equation 21 and $\theta_s \in \mathcal{H}$ is the solution of Equation 24, then $f_s = \phi(\cdot)^*\theta_s$. This equivalence could have been obtained by means of Lagrange duality. However in this proof we do not use such tool: we only focus on the representer theorem and the fact that there exists a partial isometry $W$ between the VV-RKHS and a feature space $\mathcal{H}$. We show that if $\theta_s$ is a solution of Equation 23, then $\theta_s$ belongs to $(\text{Ker } W)^{\perp}$, thus there is an isometry between $\theta_s \in \tilde{\mathcal{H}}$ and $\mathcal{H}_{\tilde{K}}$: namely $W$.

**Theorem 35 (Feature equivalence)**

Let $\tilde{K}$ be an Operator-Valued Kernel such that for all $x, z \in \mathcal{X}$, $\tilde{\phi}(x)^*\tilde{\phi}(z) = \tilde{K}(x, z)$ where $\tilde{K}$ is a $\mathcal{Y}$-Mercer OVK and $\mathcal{H}_{\tilde{K}}$ its corresponding $\mathcal{Y}$-Reproducing kernel Hilbert space. Let $c: \mathcal{Y} \times \mathcal{Y} \to \mathbb{R}$ be a cost function such that $L(x, \tilde{f}, y) = c(\tilde{f}(x), y)$ is a proper convex lower semi-continuous function in $\tilde{f} \in \mathcal{H}_{\tilde{K}}$ for all $x \in \mathcal{X}$ and all $y \in \mathcal{Y}$. Eventually let $\lambda \in \mathbb{R}_{>0}$ be the Tychonov regularization hyperparameter. The solution $f_s \in \mathcal{H}_{\tilde{K}}$ of the regularized optimization problem

$$ \tilde{f}_s = \arg \min_{\tilde{f} \in \mathcal{H}_{\tilde{K}}} \frac{1}{N} \sum_{i=1}^{N} c(\tilde{f}(x_i), y_i) + \frac{\lambda}{2} \|\tilde{f}\|_{\tilde{K}}^2 \quad (24) $$

has the form $\tilde{f}_s = \tilde{\phi}(\cdot)^*\theta_s$, where $\theta_s \in (\text{Ker } \tilde{W})^{\perp}$ and

$$ \theta_s = \arg \min_{\theta \in \tilde{\mathcal{H}}} \frac{1}{N} \sum_{i=1}^{N} c(\tilde{\phi}(x_i)^* \theta, y_i) + \frac{\lambda}{2} \|\theta\|_{\tilde{\mathcal{H}}}^2 \quad (25) $$

In the aforementioned theorem, we use the notation $\tilde{K}$ and $\tilde{\phi}$ because our main subject of interest is the ORFF map. However this theorem works for any feature maps $\phi(x) \in \mathcal{L}(\mathcal{Y}, \mathcal{H})$ even when $\mathcal{H}$ is infinite dimensional.<fnref n="7" /> This shows that when $\lambda > 0$ the solution of Equation 22 with the approximated kernel $K(x, z) \approx \tilde{K}(x, z) = \tilde{\phi}(x)^*\tilde{\phi}(z)$ is the same than the solution of Equation 25 up to an isometric isomorphism (see appendix B.3.2). Namely, if $u_s$ is the solution of Equation 22, $\theta_s$ is the solution of Equation 25 and $\lambda > 0$ we have

$$ \theta_s = \sum_{i=1}^{N} \tilde{\phi}(x_i)(u_s)_i \in (\text{Ker } W)^{\perp} \subseteq \tilde{\mathcal{H}}. $$

<footnote n="7">7. If $\phi(x) : \mathcal{L}(\mathcal{Y}, \mathcal{H})$ and $\dim(\mathcal{H}) = \infty$, the decomposition $\mathcal{H} = (\text{Ker } W) \oplus (\text{Ker } W)^{\perp}$ holds since $\mathcal{H}$ is a Hilbert space and $W$ is a bounded operator.</footnote>
If $\lambda_K = 0$ we can still find a solution $u_s$ of Equation 22. By construction of the kernel expansion, we have $u_s \in (\text{Ker } W)^\perp$. However looking at the proof of Theorem 35 we see that $\theta_s$ might not belong to $(\text{Ker } W)^\perp$. We can compute a residual vector $r_s = \sum_{i=1}^N \tilde{\phi}(x_i)(u_s)_i - \theta_s$. Since $\sum_{j=1}^N \tilde{\phi}(x_j) \in (\text{Ker } W)^\perp$ by construction, if $r_s = 0$, it means that $\lambda_K$ is large enough for both representer theorem and ORFF representer theorem to apply. If $r_s \neq 0$ but $\tilde{\phi}(\cdot)^*r_s = 0$ it means that both $\theta_s$ and $\sum_{j=1}^N \tilde{\phi}(x_j)u_s$ are in $(\text{Ker } W)^\perp$, thus the representer theorem fails to find the “true” solution over the whole space $\mathcal{H}_{\tilde{K}}$ but returns a projection onto $\mathcal{H}_{\tilde{K},s}$ of the solution. If $r_s \neq 0$ and $\tilde{\phi}(\cdot)^*r_s \neq 0$ means that $\theta_s$ is not in $(\text{Ker } W)^\perp$, thus the feature equivalence theorem fails to apply. Since $r_s = \sum_{i=1}^N \tilde{\phi}(x_i)(u_s)_i - \theta_s^\perp - \theta_s^\| \text{ and } \sum_{i=1}^N \tilde{\phi}(x_i)(u_s)_i \text{ is in } (\text{Ker } W)^\perp$, with mild abuse of notation we write $r_s = \theta^\|$. This remark is illustrated in Figure 7.

In Figure 7, we generated the data from a since wave to which we add some Gaussian noise. We learned a Gaussian kernel based RFF model (blue curve) and a kernel model (yellow curve) where the kernel is obtained from the RFF map. The left column represents the fit of the model to the points for four different valued of $\lambda$ (top to bottom: $10^{-2}$, $10^{-5}$, $10e^{10}$, 0). The middle column shows if the RFF solution $\theta_s$ is in $(\text{Ker } \tilde{W})^\perp$. This is true for all values of $\lambda$. The right column shows that even though $\theta_s$ is in $(\text{Ker } \tilde{W})^\perp$, when $\lambda \to 0$ learning with RFF is different from learning with the kernel constructed from the RFF maps since the coefficients of $\theta^\|$ are all different from 0.

## 5.3 Solving ORFF-based regression

In order to find a solution to Equation 23, we turn our attention to gradient descent methods. We define an algorithm (Algorithm 3) to find efficiently a solution to Equation 23 when $c(y, y') = \|y = y'\|^2_y$ and study its complexity.

### 5.3.1 GRADIENT METHODS

Since the solution of Equation 23 is unique when $\lambda > 0$, a sufficient and necessary condition is that the gradient of $\mathfrak{R}_\lambda$ at the minimizer $\theta_s$ is zero. We use the Frechet derivative, the strongest notion of derivative in Banach spaces (Conway, 2013; Kurdila and Zabarankin, 2006) which directly generalizes the notion of gradient to Banach spaces. The chain rule is valid in this context (Kurdila and Zabarankin, 2006, theorem 4.1.1 page 140). Hence

$$ \nabla_{\theta} c(\tilde{\phi}(x_i)^* \theta, y_i) = \tilde{\phi}(x_i) \left( \left. \frac{\partial}{\partial y} c(y, y_i) \right|_{y=\tilde{\phi}(x_i)^* \theta} \right)^*, \text{ and } \nabla_{\theta} \|\theta\|^2_{\tilde{\mathcal{H}}} = 2\theta. $$

Provided that $c(y, y_i)$ is Frechet differentiable w.r.t. $y$, for all $y$ and $y_i \in \mathcal{Y}$ we have $\nabla_{\theta}\mathfrak{R}_\lambda(\theta, s) \in \tilde{\mathcal{H}}$ and

$$ \nabla_{\theta} \mathfrak{R}_{\lambda}(\theta, s) = \frac{1}{N} \sum_{i=1}^{N} \tilde{\phi}(x_i) \left( \left. \frac{\partial}{\partial y} c(y, y_i) \right|_{y=\tilde{\phi}(x_i)^* \theta} \right)^* + \lambda \theta \quad (26) $$

**Example 5** (Naive closed form for the squared error cost)
<figure><img src="image_8.png" /><figcaption>Figure 7: ORFF equivalence theorem.</figcaption></figure>
**Algorithm 3:** Naive closed form for the squared error cost.

**Input:**

* $s = (x_i, y_i)_{i=1}^N \in (\mathcal{X} \times \mathbb{R}^p)^N$ a sequence of supervised training points,

* $\tilde{\phi}(x_i) \in \mathcal{L}(\mathbb{R}^p, \mathbb{R}^r)$ a feature map defined for all $x_i \in \mathcal{X}$,

* $\lambda \in \mathbb{R}_{>0}$ the Tychonov regularization term,

**Output:** A model $h: \mathcal{X} \to \mathbb{R}^p$, $h(x) = \tilde{\phi}(x)^\top \theta_s$. such that $\theta_s$ minimize Equation 23, where $c(y, y') = \|y - y'\|_2^2$ and $\mathbb{R}^r$ and $\mathbb{R}^p$

1. $\mathbf{P} \leftarrow \frac{1}{N} \sum_{i=1}^{N} \tilde{\phi}(x_i) \tilde{\phi}(x_i)^\top \in \mathcal{L}(\mathbb{R}^r, \mathbb{R}^r)$;

2. $\mathbf{Y} \leftarrow \frac{1}{N} \sum_{i=1}^{N} \tilde{\phi}(x_i) y_i \in \mathbb{R}^r$;

3. $\theta_s \leftarrow \text{solve}_\theta((\mathbf{P} + \lambda I_r)\theta = \mathbf{Y})$;

4. $\text{return } h: x \mapsto \tilde{\phi}(x)^\top \theta_s$;

Consider the cost function defined for all $y, y' \in \mathcal{Y}$ by $c(y, y') = \frac{1}{2}\|y - y'\|_{\mathcal{Y}}^2$. Then
$$ \left( \left. \frac{\partial c(y, y_i)}{\partial y} \right|_{y=\tilde{\phi}(x_i)^*\theta} \right)^* = (\tilde{\phi}(x_i)^*\theta - y_i). $$
Thus, since the optimal solution $\theta_s$ verifies
$$ \nabla_{\theta_s} \mathfrak{R}_\lambda(\theta_s, s) = 0 \quad \text{we have} \quad \frac{1}{N} \sum_{i=1}^N \tilde{\phi}(x_i) (\tilde{\phi}(x_i)^* \theta_s - y_i) + \lambda \theta_s = 0. $$
Therefore,

$$ \left( \frac{1}{N} \sum_{i=1}^{N} \tilde{\phi}(x_i) \tilde{\phi}(x_i)^* + \lambda I_{\tilde{\mathcal{H}}} \right) \theta_s = \frac{1}{N} \sum_{i=1}^{N} \tilde{\phi}(x_i) y_i. \quad (27) $$

Suppose that $\mathcal{Y} \subseteq \mathbb{R}^p$, and for all $x \in \mathcal{X}$, $\tilde{\phi}(x) : \mathbb{R}^r \to \mathbb{R}^p$ where all spaces are endowed with the Euclidean inner product. From this we can derive Algorithm 3 which returns the closed form solution of Equation 23 for $c(y, y') = \frac{1}{2}\|y - y'\|_2^2$.

### 5.3.2 COMPLEXITY ANALYSIS

Algorithm 3 constitutes our first step toward large-scale learning with Operator-Valued Kernels. We can easily compute the time complexity of Algorithm 3 when all the operators act on finite dimensional Hilbert spaces. Suppose that $p = \text{dim}(\mathcal{Y}) < \infty$ and for all $x \in \mathcal{X}$, $\tilde{\phi}(x) : \mathcal{Y} \to \tilde{\mathcal{H}}$ where $r = \text{dim}(\tilde{\mathcal{H}}) < \infty$ is the dimension of the redescription space $\tilde{\mathcal{H}} = \mathbb{R}^r$. Since $p$ and $r < \infty$, we view the operators $\tilde{\phi}(x)$ and $I_{\tilde{\mathcal{H}}}$ as matrices. Step 1 costs $O_t(Nr^2p)$. Steps 2 costs $O_t(Nrp)$. For step 3, the naive inversion of the operator costs $O_t(r^3)$. Eventually the overall complexity of Algorithm 3 is $O_t(r^2(Np + r))$, while the space complexity is $O_s(r^2)$.

This complexity is to compare with the kernelized solution. Let

$$ \mathbf{K}: \begin{cases} \mathcal{Y}^N \to \mathcal{Y}^N \\ u \mapsto \bigoplus_{i=1}^{N+U} \sum_{j=1}^{N+U} K(x_i, x_j) u_j \end{cases} $$
When $\mathcal{Y} = \mathbb{R}$,

$$ \mathbf{K} = \begin{pmatrix} K(x_1, x_1) & \dots & K(x_1, x_{N+U}) \\ \vdots & \ddots & \vdots \\ K(x_{N+U}, x_1) & \dots & K(x_{N+U}, x_{N+U}) \end{pmatrix} $$

is called the Gram matrix of $\mathbf{K}$. When $\mathcal{Y} = \mathbb{R}^p$, $\mathbf{K}$ is a matrix-valued Gram matrix of size $pN \times pN$ where each entry $\mathbf{K}_{ij} \in \mathcal{M}_{p,p}(\mathbb{R})$. Then the equivalent kernelized solution $u_s$ of Theorem 34 is

$$ \left( \frac{1}{N} \mathbf{K} + \lambda I_{\bigoplus_{i=1}^N \mathcal{Y}} \right) u_s = \bigoplus_{i=1}^N y_i. $$

which has time complexity $O_t(N^3p^3)$ and space complexity $O_s(N^2p^2)$. Suppose we are given a generic ORFF map (see Subsection 3.5). Then $r = 2Dp$, where $D$ is the number of samples. Hence Algorithm 3 is better that its kernelized counterpart when $r = 2Dp$ is small compared to $Np$. Thus, roughly speaking it is better to use Algorithm 3 when the number of features, $r$, required is small compared to the number of training points. Notice that Algorithm 3 has a linear complexity with respect to the number of supervised training points $N$ so it is better suited to large scale learning provided that $D$ does not grows linearly with $N$. Yet naive learning with Algorithm 3 by viewing all the operators as matrices is still problematic. Indeed learning $p$ independent models with scalar Random Fourier Features would cost $O_t(D^2p^3(N+D))$ since $r = 2Dp$. This Means that learning vector-valued function has increased the (expected) complexity from $p$ to $p^3$. However in some cases we can drastically reduce the complexity by viewing the feature-maps as linear operators rather than matrices.

## 5.4 Efficient learning with ORFF

When developing Algorithm 3 we considered that the feature map $\tilde{\phi}(x)$ was a matrix from $\mathbb{R}^p$ to $\mathbb{R}^r$ for all $x \in \mathcal{X}$, and therefore that computing $\tilde{\phi}(x)\tilde{\phi}(z)^\mathrm{T}$ has a time complexity of $O(r^2p)$. While this holds true in the most generic scenario, in many cases the feature maps present some structure or sparsity allowing to reduce the computational cost of evaluating the feature map. We focus on the Operator-valued Random Fourier Feature given by Algorithm 1, developed in Subsection 3.3 and Subsection 3.5 and treat the decomposable kernel, the curl-free kernel and the divergence-free kernel as an example. We recall that if $\mathcal{Y}' = \mathbb{R}^{p'}$ and $\mathcal{Y} = \mathbb{R}^p$, then $\tilde{\mathcal{H}} = \mathbb{R}^{2Dp'}$ thus the Operator-valued Random Fourier Features given in Section 3 have the form

$$ \left\{ \begin{array}{ll} \tilde{\phi}(x) \in \mathcal{L} \left( \mathbb{R}^p, \mathbb{R}^{2Dp'} \right) & : y \mapsto \frac{1}{\sqrt{D}} \bigoplus_{j=1}^D (x, \omega_j) B(\omega_j)^\mathrm{T} y \\ \tilde{\phi}(x)^\mathrm{T} \in \mathcal{L} \left( \mathbb{R}^{2Dp'}, \mathbb{R}^p \right) & : \theta \mapsto \frac{1}{\sqrt{D}} \sum_{j=1}^D (x, \omega_j) B(\omega_j) \theta_j \end{array} \right. $$

where $\omega_j \sim \mathbf{Pr}_{\tilde{\mathcal{H}}a\mathcal{r},p}$ i. i. d. and $B(\omega_j) \in \mathcal{L}(\mathbb{R}^p, \mathbb{R}^{p'})$ for all $\omega_j \in \tilde{\mathcal{X}}$. Hence the Operator-valued Random Fourier Feature can be seen as the block matrix $\in \mathcal{M}_{2Dp',p}(\mathbb{R})$

$$ \tilde{\phi}(x) = \left( \cos\langle x, \omega_1 \rangle B(\omega_1) \quad \sin\langle x, \omega_1 \rangle B(\omega_1) \quad \dots \quad \cos\langle x, \omega_D \rangle B(\omega_D) \quad \sin\langle x, \omega_D \rangle B(\omega_D) \right)^\mathrm{T} \quad (28) $$
### 5.4.1 CASE OF STUDY: THE DECOMPOSABLE KERNEL

Throughout this section we show how the mathematical formulation relates to a concrete (Python) implementation. We propose a Python implementation based on NumPy (Oliphant, 2006), SciPy (Jones et al., 2014) and Scikit-learn (Pedregosa et al., 2011). Following Equation 28, the feature map associated to the decomposable kernel would be

$$ \tilde{\phi}(x) = \frac{1}{\sqrt{D}} \left( \begin{array}{cc} \cos\langle x, \omega_1 \rangle B & \sin\langle x, \omega_1 \rangle B \\ \vdots & \vdots \\ \cos\langle x, \omega_D \rangle B & \sin\langle x, \omega_D \rangle B \end{array} \right)^{\mathrm{T}} \\ = \frac{1}{\sqrt{D}} \left( \begin{array}{cc} \cos\langle x, \omega_1 \rangle & \sin\langle x, \omega_1 \rangle \\ \vdots & \vdots \\ \cos\langle x, \omega_D \rangle & \sin\langle x, \omega_D \rangle \end{array} \right)^{\mathrm{T}} \otimes B^{\mathrm{T}}, $$

$\omega_j \sim \mathbf{Pr}_{\widehat{\text{Haar}},\rho}$ i. i. d., which would lead to the following naive python implementation for the Gaussian (RBF) kernel of parameter $\gamma$, whose associated spectral distribution is $\mathbf{Pr}_{\rho} = \mathcal{N}(0, 2\gamma)$. Let $\theta \in \mathbb{R}^{2Dp'}$ and $y \in \mathbb{R}^1$. With such implementation evaluating a matrix vector product such as $\tilde{\phi}(x)^{\mathrm{T}}\theta$ or $\tilde{\phi}(x)y$ have $O_t(2Dp'p)$ time complexity and $O_s(2Dp'p)$ of space complexity, which is utterly inefficient. Indeed, recall that if $B \in \mathcal{M}_{p,p'}(\mathbb{R}^{p'})$ is matrix, the operator $\tilde{\phi}(x)$ corresponding to the decomposable kernel is

$$ \tilde{\phi}(x)y = \frac{1}{\sqrt{D}} \bigoplus_{j=1}^{D} \left( \begin{array}{c} \cos\langle x, \omega_j \rangle B^{\mathrm{T}} y \\ \sin\langle x, \omega_j \rangle B^{\mathrm{T}} y \end{array} \right) = \left( \frac{1}{\sqrt{D}} \bigoplus_{j=1}^{D} \left( \begin{array}{c} \cos\langle x, \omega_j \rangle \\ \sin\langle x, \omega_j \rangle \end{array} \right) \right) \otimes (B^{\mathrm{T}} y) \quad (29) $$

and

$$ \tilde{\phi}(x)^{\mathrm{T}}\theta = \frac{1}{\sqrt{D}} \sum_{j=1}^{D} \cos\langle x, \omega_j \rangle B\theta_j + \sin\langle x, \omega_j \rangle B\theta_j = B \left( \frac{1}{\sqrt{D}} \sum_{j=1}^{D} (\cos\langle x, \omega_j \rangle + \sin\langle x, \omega_j \rangle) \theta_j \right). \quad (30) $$

Which requires only evaluation of $B$ on $y$ and can be implemented easily in Python thanks to SciPy's LinearOperator. Note that the computation of these expressions can be fully vectorized<fnref n="8" /> using the vectorization property of the Kronecker product. In the following we consider $\Theta \in \mathcal{M}_{2D,u'}(\mathbb{R})$ and the operator $\mathrm{vec}: \mathcal{M}_{p',2D}(\mathbb{R}) \to \mathbb{R}^{2Dp'}$ which turns a matrix into a vector (i. e. $\theta_{p'i+j} = \mathrm{vec}(\Theta_{ij})$, $i \in \mathbb{N}_{(2D-1)}$ and $j \in \mathbb{N}_{p'}^*$). Then $(\tilde{\phi}(x) \otimes B^{\mathrm{T}})^{\mathrm{T}}\theta = (\tilde{\phi}(x)^{\mathrm{T}} \otimes B) \mathrm{vec}(\Theta) = \mathrm{vec}(B\Theta\tilde{\phi}(x))$. with this trick, many authors (Sindhwani et al., 2013; Brault et al., 2016; Rosasco et al., 2010; Carmeli et al., 2010) notice that the decomposable kernel usually yields a Stein equation (Penzl, 1998). Indeed rewriting step 3 of Algorithm 3 gives a system to solve of the form

$$ \tilde{\phi}(X)\tilde{\phi}(X)^{\mathrm{T}}\Theta B^{\mathrm{T}}B + \lambda\Theta - Y = 0 \Leftrightarrow (\tilde{\phi}(X)\tilde{\phi}(X)^{\mathrm{T}} \otimes B^{\mathrm{T}}B + \lambda I_{2Dp'})\Theta - Y = 0 $$

Many solvers exists to solve efficiently this kind of systems<fnref n="9" />, but most of them share the particularity that they are not just restricted to handle Stein equations. Broadly speaking,

<footnote n="8">8. See Walt et al. (2011).</footnote>

<footnote n="9">9. For instance Sleijpen et al. (2010).</footnote>
iterative solvers (or matrix free solvers) are designed to solve any systems of equation of
ther form $PX = C$, where $P$ is a linear operator (not a matrix). This is exacly our case
where $\tilde{\varphi}(x) \otimes B^T$ is the matrix form of the operator $\Theta \mapsto \text{vec}(B\Theta\tilde{\varphi}X)$.

This leads us to the following (more efficient) Python implementation of the Decomposable
ORFF "operator" to be feed to a matrix-free solvers.

```python
def EfficientDecomposableGaussianORFF(X, A, gamma=1.,
                                     D=100, eps=1e-5, random_state=0):
    r"""Return the efficient ORFF map associated with the data X.

    Parameters
    -----------
    X : {array-like}, shape = [n samples, n features}
        Samples.
    A : {array-like}, shape = [n targets, n targets]
        Operator of the Decomposable kernel (positive semi-definite)
    gamma : {float},
        Gamma parameter of the RBF kernel.
    D : {integer}
        Number of random features.
    eps : {float}
        Cutoff threshold for the singular values of A.
    random_state : {integer}
        Seed of the generator.

    Returns
    -------
    \tilde{t} \tilde{d} \{\Phi}(X) : Linear Operator, callable
    """
    # Decompose A=BB^T\ transpose
    u, s, v = svd(A, full matrices=False, compute_uv=True)
    B = dot(diag(sqrts(s[s > eps]))), v[s > eps, :])

    # Sample a RFF from the scalar Gaussian kernel
    phi_s = RBFSampler(gamma=gamma, n_components=D, random_state(random_state)
    phiX = phi_s.fit_transform(X)

    # Create the ORFF linear operator
    return LinearOperator((phiX.shape[0] * B.shape[1], D * B.shape[0]),
                          matvec=lambda b: dot phiX, dot(b reshape((D, B.shape[0附加]), B)),
                          rmatvec=lambda r: dot phiX.T, dot(r reshape((X.shape[0], B.shape[1附加]), B.T)))
```

It is worth mentioning that the same strategy can be applied in many different language.
For instance in C++, the library Eigen (Guennebaud et al., 2010) allows to wrap a sparse
matrix with a custom type, where the user overloads the transpose and dot product
operator (as in Python). Then the custom user operator behaves as a (sparse) matrix –
see https://eigen.tuxfamily.org/dox/group__MatrixfreeSolverExample.html. With
this implementation the time complexity of $\tilde{\phi}(x)^T\theta$ and $\tilde{\phi}(x)y$ falls down to $O_t((D + p)p')$
and the same holds for space complexity.

A quick experiment shows the advantage of seeing the decomposable kernel as a linear
operator rather than a matrix. We draw $N = 100$ points $(x_i)_{i=1}^N$ in the interval $(0, 1)^{20}$ and
use a decomposable kernel with matrix $\Gamma = BB^T \in \mathcal{M}_{p,p}(\mathbb{R})$ where $B \in \mathcal{M}_{p,p}(\mathbb{R})$ is a random
matrix with coefficients drawn uniformly in $(0, 1)$. We compute $\tilde{\phi}(x)^T\theta$ for all $x_i$'s, where
$\theta \in \mathcal{M}_{2D,1}(\mathbb{R})$, $D = 100$, with the implementation EfficientDecomposableGaussianORFF,
Equation 30. The coefficients of $\theta$ were drawn at random uniformly in $(0, 1)$. We report the
<figure><img src="image_9.png" /><figcaption>Figure 8: Efficient decomposable Gaussian ORFF (lower is better).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Efficient curl-free Gaussian ORFF (lower is better).</figcaption></figure>

execution time in Figure 8 for different values of $p$, $1 \le p \le 100$. The left plot reports the execution time in seconds of the construction of the feature. The middle plot reports the execution time of $\tilde{\phi}(x)^{\top}\theta$, and the right plot the memory used in bytes to store $\tilde{\phi}(x)$ for all $x_i$'s. We averaged the results over ten runs.

**Curl-free kernel.** We use the unbounded ORFF map presented in Equation 17. We draw $N = 1000$ points $(x_i)_{i=1}^N$ in the interval $(0, 1)^p$ and use a curl-free kernel. We compute $\tilde{\phi}(x)^{\top}\theta$ for all $x_i$'s, where $\theta \in \mathcal{M}_{2D,1}(\mathbb{R})$, $D = 500$, with the matrix implementation and the LinearOperator implementation. The coefficients of $\theta$ were drawn at random uniformly in $(0, 1)$. We report the execution time in Figure 9 for different values of $p$, $1 \le p \le 100$. The left plot reports the execution time in seconds of the construction of the features. The middle plot reports the execution time of $\tilde{\phi}(x)^{\top}\theta$, and the right plot the memory used in bytes to store $\tilde{\phi}(x)$ for all $x_i$'s. We averaged the results over fifty runs. As we can see the
<figure><img src="image_11.png" /><figcaption>Figure 10: Efficient divergence-free Gaussian ORFF (lower is better).</figcaption></figure>

linear-operator implementation is one order of magnitude slower than its matrix counterpart. However it uses considerably less memory.

**Divergence-free kernel.** We use the unbounded ORFF map presented in Equation 18. We draw $N = 100$ points $(x_i)_{i=1}^N$ in the interval $(0, 1)^p$ and use a curl-free kernel. We compute $\tilde{\phi}(x)^T\theta$ for all $x_i$'s, where $\theta \in \mathcal{M}_{2D_p,1}(\mathbb{R})$, $D = 100$, with the matrix implementation and the LinearOperator implementation. The coefficients of $\theta$ were drawn at random uniformly in $(0,1)$. We report the execution time in Figure 9 for different values of $p$, $1 \le p \le 100$. The left plot reports the execution time in seconds of the construction of the feature. The middle plot reports the execution time of $\tilde{\phi}(x)^T\theta$, and the right plot the memory used in bytes to store $\tilde{\phi}(x)$ for all $x_i$'s. We averaged the results over ten runs. We draw the same conclusions as the curl-free kernel.

# 6. Numerical experiments

We present a set of experiments to complete the theoretical contribution and illustrate the behavior of ORFF-regression. First we study how well the ORFF regression recover the result of operator-valued kernel regression. Second we show the advantages of ORFF regression over independent RFF regression. A code implementing ORFF is available at https://github.com/operalib/operalib a framework for OVK Learning.

## 6.1 Learning with ORFF vs learning with OVK

### 6.1.1 DATASETS

The first dataset considered is the handwritten digits recognition dataset MNIST<fnref n="10" />. We select a training set of 12,000 images and a test set of 10,000 images. The inputs are images represented as a vector $x_i \in [0, 255]^{784}$ and the targets $y_i \in \mathbb{N}_9$ are integers between 0 and 9. First we scaled the inputs such that they take values in $[-1, 1]^{784}$. Then we binarize

<footnote n="10">10. available at http://yann.lecun.com/exdb/mnist</footnote>
the targets such that each number is represented by a unique binary vector of dimension 10. The vector $y_i$ is zero everywhere except on the dimension corresponding to the class where it is one. For instance the class 4 is encoded $(0 \ 0 \ 0 \ 0 \ 1 \ 0 \ 0 \ 0 \ 0 \ 0)^T$. To predict classes, we use the simplex coding method presented in Mroueh et al. (2012). The intuition behind simplex coding is to project the binarized labels of dimension $p$ onto the most separated vectors on the hypersphere of dimension $p-1$. For ORFF we can encode directly this projection in the $B$ matrix of the decomposable kernel $K_0(\delta) = BB^*k_0(\delta)$ where $k_0$ is a Gaussian kernel. The matrix $B$ is computed via the recursion

$$B_{p+1} = \begin{pmatrix} 1 & u^T \\ 0_{p-1} & \sqrt{1-p^{-2}}B_p \end{pmatrix}, \quad B_2 = \begin{pmatrix} 1 & -1 \end{pmatrix},$$

where $u = (-p^{-2} \ \dots \ -p^{-2})^T \in \mathbb{R}^{p-1}$ and $0_{p-1} = (0 \ \dots \ 0)^T \in \mathbb{R}^{p-1}$. For Operator-Valued Kernels we project the binarized targets on the simplex as a preprocessing step, before learning with the decomposable $K_0(\delta) = I_p k_0(\delta)$, where $k_0$ is a scalar Gaussian kernel.

The second dataset is a simulated five dimensional (5D) vector field with structure. We generate a scalar field as a random function $f: [-1, 1]^5 \to \mathbb{R}$, where $\tilde{f}(x) = \tilde{\varphi}(x)^*\theta$ where $\theta$ is a random matrix with each entry following a standard normal distribution, $\tilde{\varphi}$ is a scalar Gaussian RFF with bandwidth $\sigma = 0.4$. The input data $x$ are generated from a uniform probability distribution. We take the gradient of $\tilde{f}$ to generate the curl-free 5D vector field.

The third dataset is a synthetic of data from $\mathbb{R}^{20} \to \mathbb{R}^4$ as described in Audiffren and Kadri (2015). In this dataset, inputs $(x_1, \dots, x_{20})$ are generated independently and uniformly over $[0, 1]$ and the different outputs are computed as follows. Let $\varphi(x) = (x_1^2, x_4^2, x_1x_2, x_3x_5, x_2, x_4, 1)$ and $(w_i)$ denotes the i. i. d. copies of a seven dimensional Gaussian distribution with zero mean and covariance $\Sigma \in \mathcal{M}_{7,7}(\mathbb{R})$ such that $(0.5 \ 0.25 \ 0.1 \ 0.05 \ 0.15 \ 0.1 \ 0.15)$. Then, the outputs of the different tasks are generated as $y_i = w_i\varphi(x)$. We use this dataset with $p = 4, 10^5$ instances and for the train set and also $10^5$ instances for the test set.

### 6.1.2 RESULTS

**Performance of ORFF regression on the first dataset.** We trained both ORFF and OVK models on MNIST dataset with a decomposable Gaussian kernel with signature $K_0(\delta) = \exp(-\|\delta\|/(2\sigma^2))\Gamma$. To apply Algorithm 3 after noticing that in the case of the decomposable kernel with $\lambda_M = 0$, it boilds down to a Stein equation (Brault et al., 2016, section 5.1), we use an off-the-shelf solver<fnref n="11" /> able to handle Stein's equation. For both methods we choose $\sigma = 20$ and use a 2-fold cross validation on the training set to select the optimal $\lambda$. First, Figure 11 compares the running time between OVK and ORFF models using $D = 1000$ Fourier features against the number of datapoints $N$. The log-log plot shows ORFF scaling better than the OVK w.r.t. the number of points. Second, Figure 11 shows the test prediction error versus the number of ORFFs $D$, when using $N = 1000$ training

<footnote n="11">11. Available at http://ta.twi.tudelft.nl/nw/users/gijzen/IDR.html</footnote>
<figure><img src="image_12.png" /><figcaption>Figure 11: Empirical comparison of ORFF and OVK regression on MNIST dataset and empirical behavior of ORFF regression versus *D* and *N*.</figcaption></figure>

points. As expected, the ORFF model converges toward the OVK model when the number of features increases.

**Performance of ORFF regression on the second dataset.** We perform a similar experiment on the second dataset (5D-vector field with structure). We use a Gaussian curl-free kernel with bandwidth equal to the median of the pairwise distances and tune the hyperparameter $\lambda$ on a grid. Here we optimize Equation 23, where $c$ is the squared error cost, using Scipy's L-BFGS-B (Byrd et al., 1995) solver<fnref n="12" /> with the gradients given in Equation 26 and the efficient linear operator described in Subsection 5.4 (e. g. Equations 29 and 30). Figure 12 (bottom row) reports the $R^2$ (coefficient of determination) score on the test set versus the number of curl-ORFF $D$ with a comparison with curl-OVK. In this experiment, we see that curl-ORFF can even be better than curl-OVK, suggesting that ORFF might play an additional regularizing role. It also shows the computation time of curl-ORFF and curl-OVK. We see that OVK regression does not scale with large datasets, while ORFF regression does. When $N > 10^4$, OVK regression exceeds memory capacity.

**Structured prediction vs Independent (RFF) prediction.** On the second dataset, Figure 12 (top row) compares $R^2$ score and time of ORFF regression using the trivial identity decomposable kernel, e. g. independent RFFs, to curl-free ORFF regression. Curl-free ORFF outperforms independent RFFs, as expected, since the dataset involves structured outputs.

**Impact of the number of random features ($D$).** In this setting we solved the optimization problem for both ORFF and OVK using a L-BFGS-B. Figure 13 top row shows that for a fixed number of instance in the train set, OVK performs better than ORFF in terms of accuracy ($R^2$). However ORFF scales better than OVK w.r.t. the number of data. ORFF is able to process more data than OVK in the same time and thus reach a better accuracy for a given amount of time. Bottom row shows that ORFF tends to reach OVK's accuracy for a fixed number of data when the number of features increase.

**Multitask learning.** In this experiment we are interested in multitask learning with operator-valued random Fourier features, and see whether the approximation of a joint OVK performs better than an independent OVK. In this setting we assume that for each entry $x_i \in \mathbb{R}^d$ we only have access to one observation $y_i \in \mathbb{R}$ corresponding to a task $t_i$. We used

<footnote n="12">12. Available at http://docs.scipy.org/doc/scipy/reference/optimize.html</footnote>
<figure><img src="image_13.png" /><figcaption>Figure 12: Empirical comparison between curl-free ORFF, curl-free OVK, independent ORFF, independent OVK on a synthetic vector field regression task.</figcaption></figure>

the SARCOS dataset, taken from http://www.gaussianprocess.org/gpml/data/ website. This is an inverse dynamics problem, i. e. we have to predict the 7 joint torques given the joint positions, velocities and accelerations. Hence, we have to solve a regression problem with 21 inputs and 7 outputs which is a very nonlinear function. It has 45K inputs data. Suppose that we are given a collection of inputs data $x_1, ..., x_N \in \mathbb{R}^{21}$ and a collection of output data $((y_1, t_1), ..., (y_N, t_N)) \in (\mathbb{R} \times \mathbb{N}_T)^N$ where T is the number of tasks. We consider the following multitask loss function $L(h(x), (y, t)) = \frac{1}{2}(\langle h(x), e_t \rangle_2 - y)^2$, This loss function is adapted to datasets where the number of data per tasks is unbalanced (i. e. for one input data we observe the value of only one task and not all the tasks.). We optimise the regularized risk

$$ \frac{1}{N} \sum_{i=1}^{N} L(h(x_i), (y_i, t_i)) + \frac{\lambda}{2N} ||h||_{\mathcal{H}}^2 = \frac{1}{2N} \sum_{i=1}^{N} (\langle h(x_i), e_{t_i} \rangle - y_i)^2 + \frac{\lambda}{2N} ||h||_{\mathcal{H}}^2 $$

We used a model h based on the decomposable kernel $h(x) = (\varphi(x)^T \otimes B)\theta$ we chose B such that $BB^T = A$, where A is the inverse graph Laplacian L of the similarities between the tasks, parametrized by an hyperparameter $\gamma \in \mathbb{R}_+$. $L_{kl} = \exp(-\gamma \sqrt{\sum_{i=1}^{N} (y_i^k - y_i^l)^2})$. We draw N data randomly for each task, hence creating a dataset of $N \times 7$ data and computed the nMSE on the proposed test set (4.5K points). We repeated the experiments 80 times to avoid randomness. We choose $D = \frac{\max(N, 500)}{2}$ features, and optimized the problem with a second order batch gradient. Table 3 shows that using the ORFF approximation of an operator-valued kernel with a good prior on the data improves the performances w.r.t.
<figure><img src="image_14.png" /><figcaption>Figure 13: Decomposable kernel on the third dataset: $R^2$ score vs number of data in the train set (N)</figcaption></figure>

<table><caption>Table 3: Error (% of nMSE) on SARCOS dataset.</caption><thead><tr><th>N</th><th>Independent (%)</th><th>Laplacian (%)</th><th>p-value</th><th>T</th></tr></thead><tbody><tr><td>50 × 7</td><td>23.138 ± 0.577</td><td>22.254 ± 0.536</td><td>2.68%</td><td>4(s)</td></tr><tr><td>100 × 7</td><td>16.191 ± 0.221</td><td>15.568 ± 0.187</td><td>&lt; 0.1%</td><td>16(s)</td></tr><tr><td>200 × 7</td><td>12.713 ± 0.0978</td><td>12.554 ± 0.0838</td><td>1.52%</td><td>12(s)</td></tr><tr><td>400 × 7</td><td>10.785 ± 0.0579</td><td>10.651 ± 0.0466</td><td>&lt; 0.1%</td><td>10(s)</td></tr><tr><td>800 × 7</td><td>7.512 ± 0.0344</td><td>7.512 ± 0.0344</td><td>100%</td><td>15(s)</td></tr><tr><td>3200 × 7</td><td>5.658 ± 0.0187</td><td>5.658 ± 0.0187</td><td>100%</td><td>20(s)</td></tr></tbody></table>

the independent ORFF. However the advantage seems to be less important more data are available.

## 7. Conclusion

OVKs naturally extend the celebrated kernel method used to learn scalar-valued functions, to the case of learning vector-valued functions. Although OVKs are appealing from a theoretical aspect, these methods scale poorly in terms of computation time when the number of data is high. Indeed, to evaluate the value of function with an Operator-Valued Kernel, it requires to evaluate an Operator-Valued Kernel on all the point in the given dataset. Hence naive learning with kernels usually scales cubically in time with the number
<figure><img src="image_15.png" /><figcaption>Figure 14: Decomposable kernel on the third dataset: $R^2$ score vs number of data in the train set (N) for different number of random samples (D).</figcaption></figure>

of data. In the context of large-scale learning such scaling is not acceptable. Through this work we propose a methodology to tackle this difficulty.

Enlightened by the literature on large-scale learning with *scalar*-valued kernel, in particular the work of Rahimi and Recht (Rahimi and Recht, 2007), we propose to replace an OVK by a random feature map that we called Operator-valued Random Fourier Feature. Our contribution starts with the formal mathematical construction of this feature from an OVK. Then we show that it is also possible to obtain a kernel from an ORFF. Eventually we analyse the regularization properties in terms of Fourier Transform of $\mathcal{Y}$-Mercer kernels. Then we moved on giving a bound on the error due to the random approximation of the OVK with high probability. We showed that it is possible to bound the error even though the ORFF estimator of an OVK is not a bounded random variable. Moreover we also give a bound when the dimension of the output data infinite.

After ensuring that an ORFF is a good approximation of a kernel, we moved on giving a framework for supervised learning with Operator-Valued Kernels. We showed that learning with a feature map is equivalent to learn with the reconstructed OVK under some mild conditions. Then we focused on an efficient implementation of ORFF by viewing them as linear operators rather than matrices and using matrix-free (iterative) solvers and concluded with some numerical experiments.

Following Rahimi and Recht a generalization bound for ORFF kernel ridge would probably suggest that the number of feature to draw is proportional to the number of data. However new results of Rudi et al. (2016) suggest that the number of feature should be proportional
to the square root of the number of data. In a future work, we shall investigate this results and extend it to ORFF.

Since the construction of ORFF is valid for infinite dimensional Hilbert spaces such as function spaces, we would also like to investigate learning function valued functions in an efficient manner.

## Acknowledgments

The authors are grateful to Maxime Sangnier (UPMC, France) for the insightful discussions and Markus Heinonen (Aalto University, Sweden) for the preliminary experiments.

# Appendix A. Reminder on Abstract Harmonic Analysis

## A.1 Locally compact Abelian groups

**Definition 36** (Locally Compact Abelian (LCA) group.)

A group $\mathcal{X}$ endowed with a binary operation $\star$ is said to be a Locally Compact Abelian group if $\mathcal{X}$ is a topological commutative group w.r.t. $\star$ for which every point has a compact neighborhood and is Hausdorff (T2).

Moreover given a element $z$ of a LCA group $\mathcal{X}$, we define the set $z \star \mathcal{X} = \mathcal{X} \star z = \{z \star x | \forall x \in \mathcal{X}\}$ and the set $\mathcal{X}^{-1} = \{x^{-1} | \forall x \in \mathcal{X}\}$. We also note $e$ the neutral element of $\mathcal{X}$ such that $x \star e = e \star x = e$ for all $x \in \mathcal{X}$. Throughout this paper we focus on positive definite function. Let $\mathcal{Y}$ be a complex separable Hilbert space. A function $f: \mathcal{X} \to \mathcal{Y}$ is positive definite if for all $N \in \mathbb{N}$ and all $y \in \mathcal{Y}$,

$$ \sum_{i,j=1}^{N} \langle y_i, f(x_j^{-1} \star x_i) y_j \rangle_{\mathcal{Y}} \geq 0 \qquad (31) $$

for all sequences $(y_i)_{i \in \mathbb{N}_N^*} \in \mathcal{Y}^N$ and all sequences $(x_i)_{i \in \mathbb{N}_N^*} \in \mathcal{X}^N$. If $\mathcal{Y}$ is real we add the assumption that $f(x^{-1}) = f(x)^*$ for all $x \in \mathcal{X}$.

## A.2 Even and odd functions

Let $\mathcal{X}$ be a LCA group and $\mathbb{K}$ be a field viewed as an additive group. We say that a function $f: \mathcal{X} \to \mathbb{K}$ is even if for all $x \in \mathcal{X}$, $f(x) = f(x^{-1})$ and odd if $f(x) = -f(x^{-1})$. The definition can be extended to operator-valued functions.

**Definition 37** (Even and odd operator-valued function on a LCA group)

Let $\mathcal{X}$ be a measured LCA group and $\mathcal{Y}$ be a Hilbert space, and $\mathcal{L}(\mathcal{Y})$ the space of bounded linear operators from $\mathcal{Y}$ to itself viewed as an additive group. A function $f: \mathcal{X} \to \mathcal{L}(\mathcal{Y})$ is (weakly) even if for all $x \in \mathcal{X}$ and all $y, y' \in \mathcal{Y}$, $\langle y, f(x^{-1}) y' \rangle_{\mathcal{Y}} = \langle y, f(x)y' \rangle_{\mathcal{Y}}$ and (weakly) odd if $\langle y, f(x^{-1}) y' \rangle_{\mathcal{Y}} = -\langle y, f(x)y' \rangle_{\mathcal{Y}}$.
<table><caption>Table 4: Classification of Fourier Transforms in terms of their domain and transform domain.</caption><thead><tr><th>X =</th><th>Х̃ ≅</th><th>Operation</th><th>Pairing</th></tr></thead><tbody><tr><td>R<sup>d</sup></td><td>R<sup>d</sup></td><td>+</td><td>(x, ω) = exp (i⟨x, ω⟩<sub>2</sub>)</td></tr><tr><td>R<sup>d</sup><sub>*,+</sub></td><td>R<sup>d</sup></td><td>·</td><td>(x, ω) = exp (i⟨log(x), ω⟩<sub>2</sub>)</td></tr><tr><td>(-c; +∞)<sup>d</sup></td><td>R<sup>d</sup></td><td>⊙</td><td>(x, ω) = exp (i⟨log(x + c), ω⟩<sub>2</sub>)</td></tr></tbody></table>

It is easy to check that if $f$ is odd then $\int_X \langle y, f(x)y' \rangle_Y d\text{Haar}(x) = 0$. Besides the product of an even and an odd function is odd. Indeed for all $f, g \in \mathcal{F}(X; \mathcal{L}(Y))$, where $f$ is even and $g$ odd. Define $h(x) = \langle y, f(x)g(x)y' \rangle$. Then we have $h(x^{-1}) = \langle y, f(x^{-1})g(x^{-1})y' \rangle = \langle y, f(x)(-g(x))y' \rangle = -h(x)$.

### A.3 Characters

LCA groups are central to the general definition of Fourier Transform which is related to the concept of Pontryagin duality (Folland, 1994). Let $(X, \star)$ be a LCA group with $e$ its neutral element and the notation, $x^{-1}$, for the inverse of $x \in X$. A character is a complex continuous homomorphism $\omega: X \to U$ from $X$ to the set of complex numbers of unit module $U$. The set of all characters of $X$ forms the Pontryagin dual group $\hat{X}$. The dual group of an LCA group is an LCA group so that we can endow $\hat{X}$ with a “dual” Haar measure noted $\overline{\text{Haar}}$. Then the dual group operation is defined by $(\omega_1 \star' \omega_2)(x) = \omega_1(x)\omega_2(x) \in U$. The Pontryagin duality theorem states that $\hat{\hat{X}} \cong X$. i.e. there is a canonical isomorphism between any LCA group and its double dual. To emphasize this duality the following notation is usually adopted: $\omega(x) = (x, \omega) = (\omega, x) = x(\omega)$, where $x \in X \cong \hat{\hat{X}}$ and $\omega \in \hat{\hat{X}}$. The form $(\cdot, \cdot)$ defined in appendix A.3 is called (duality) pairing. Another important property involves the complex conjugate of the pairing which is defined as $(\overline{x}, \omega) = (x^{-1}, \omega) = (x, \omega^{-1})$. We notice that for any pairing depending of $\omega$, there exists a function $h_\omega: X \to \mathbb{R}$ such that $(x, \omega) = \exp(ih_\omega(x))$ since any pairing maps into $U$. Moreover, $(x \star z^{-1}, \omega) = \omega(x)\omega(z^{-1}) = \exp(+ih_\omega(x))\exp(+ih_\omega(z^{-1})) = \exp(+ih_\omega(x))\exp(-ih_\omega(z))$. Table 4 provides an explicit list of pairings for various groups based on $\mathbb{R}^d$ or its subsets. The interested reader can refer to Folland (1994) for a more detailed construction of LCA, Pontryagin duality and Fourier Transforms on LCA.

### A.4 The Fourier Transform

For a function with values in a separable Hilbert space, $f \in L^1(X, \text{Haar}; Y)$, we denote $\mathcal{F}[f]$ its Fourier Transform (FT) which is defined by

$$ \forall \omega \in \hat{\mathcal{X}}, \quad \mathcal{F}[f](\omega) = \int_{\mathcal{X}} \overline{(x, \omega)} f(x) d\text{Haar}(x). $$
The Inverse Fourier Transform (IFT) of a function $g \in L^1(\hat{\mathcal{X}}, \widehat{\text{Haar}}; \mathcal{Y})$ is noted $\mathcal{F}^{-1}[g]$ defined by $\forall x \in \mathcal{X}, \mathcal{F}^{-1}[g](x) = \int_{\hat{\mathcal{X}}} (x, \omega)g(\omega)d\widehat{\text{Haar}}(\omega)$, We also define the flip operator $\mathcal{R}$ by $(\mathcal{R}f)(x) := f(x^{-1})$.

**Theorem 38 (Fourier inversion)**

Given a measure **Haar** defined on $\mathcal{X}$, there exists a unique suitably normalized dual measure $\widehat{\text{Haar}}$ on $\hat{\mathcal{X}}$ such that for all $f \in L^1(\mathcal{X}, \text{Haar}; \mathcal{Y})$ and if $\mathcal{F}[f] \in L^1(\hat{\mathcal{X}}, \widehat{\text{Haar}}; \mathcal{Y})$ we have

$$f(x) = \int_{\hat{\mathcal{X}}} (x, \omega) \mathcal{F}[f](\omega) d\widehat{\text{Haar}}(\omega), \quad \text{for Haar-almost all } x \in \mathcal{X}. \qquad (32)$$

i. e. such that $(\mathcal{R}\mathcal{F}\mathcal{F}[f])(x) = \mathcal{F}^{-1}\mathcal{F}[f](x) = f(x)$ for Haar-almost all $x \in \mathcal{X}$. If $f$ is continuous this relation holds for all $x \in \mathcal{X}$.

Thus when a Haar measure **Haar** on $\mathcal{X}$ is given, the measure on $\hat{\mathcal{X}}$ that makes Theorem 38 true is called the dual measure of **Haar**, noted $\widehat{\text{Haar}}$. Let $c \in \mathbb{R}_*$ If $c\text{Haar}$ is the measure on $\mathcal{X}$, then $c^{-1}\widehat{\text{Haar}}$ is the dual measure on $\hat{\mathcal{X}}$. Hence one must replace $\widehat{\text{Haar}}$ by $c^{-1}\widehat{\text{Haar}}$ in the inversion formula to compensate. Whenever $\widehat{\text{Haar}} = \text{Haar}$ we say that the Haar measure is self-dual. For the familiar case of a scalar-valued function $f$ on the LCA group $(\mathbb{R}^d, +)$, we have for all $\omega \in \hat{\mathcal{X}} = \mathbb{R}^d$

$$\mathcal{F}[f](\omega) = \int_{\mathcal{X}} \overline{(x, \omega)} f(x) d\widehat{\text{Haar}}(x) = \int_{\mathbb{R}^d} \exp(-i\langle x, \omega \rangle_2) f(x) d\text{Leb}(x), \qquad (33)$$

the Haar measure being here the Lebesgue measure. Notice that the normalization factor of $\widehat{\text{Haar}}$ on $\hat{\mathcal{X}}$ depends on the measure **Haar** on $\mathcal{X}$ and the duality pairing. For instance let $\mathcal{X} = (\mathbb{R}^d, +)$. If one endow $\mathcal{X}$ with the Lebesgue measure as the Haar measure, the Haar measure on the dual is defined for all $\mathcal{Z} \in \mathcal{B}(\mathbb{R}^d)$ by

$$\text{Haar}(\mathcal{Z}) = \text{Leb}(\mathcal{Z}), \quad \text{and} \quad \widehat{\text{Haar}}(\mathcal{Z}) = \frac{1}{(2\pi)^d} \text{Leb}(\mathcal{Z}),$$

in order to have $\mathcal{F}^{-1}\mathcal{F}[f] = f$. If one use the cleaner equivalent pairing $(x, \omega) = \exp(2i\pi\langle x, \omega \rangle_2)$ rather than $(x, \omega) = \exp(i\langle x, \omega \rangle_2)$, then $\widehat{\text{Haar}}(\mathcal{Z}) = \text{Leb}(\mathcal{Z})$. The pairing $(x, \omega) = \exp(2i\pi\langle x, \omega \rangle_2)$ looks more attractive in theory since it limits the messy factor outside the integral sign and make the Haar measure self-dual. However it is of lesser use in practice since it yields additional unnecessary computation when evaluating the pairing. Hence for symmetry reason on $(\mathbb{R}^d, +)$ and reduce computations we settle with the Haar measure on $\mathbb{R}^d$ groups (additive and multiplicative) defined as $\widehat{\text{Haar}}(\mathcal{Z}) = \text{Haar}(\mathcal{Z}) = \sqrt{2\pi}^{-d}\text{Leb}(\mathcal{Z})$. We conclude this subsection by recalling the injectivity property of the Fourier Transform.

**Corollary 39 (Fourier Transform injectivity)**

Given $\mu$ and $\nu$ two measures, if $\mathcal{F}[\mu] = \mathcal{F}[\nu]$ then $\mu = \nu$. Moreover given two functions $f$ and $g \in L^1(\mathcal{X}, \text{Haar}; \mathcal{Y})$ if $\mathcal{F}[f] = \mathcal{F}[g]$ then $f = g$

**Appendix B. Proofs**

In this section we give the proofs of our contributions stated in the main body of the paper.
## B.1 Construction

### B.1.1 PROOF OF LEMMA 16

**Proof** For any function $f$ on $(\mathcal{X}, \star)$ define the flip operator $\mathcal{R}$ by $(\mathcal{R}f)(x) := f(x^{-1})$. For any shift invariant $\mathcal{Y}$-Mercer kernel and for all $\delta \in \mathcal{X}$, $K_e(\delta) = K_e(\delta^{-1})^*$. Indeed from the definition of a shift-invariant kernel, $K_e(\delta^{-1}) = K(\delta^{-1}, e) = K(e, \delta) = K(\delta, e)^* = K_e(\delta)^*$.

Item 1: taking the Fourier Transform yields, $\langle y', C(\omega)y \rangle_{\mathcal{Y}} = \mathcal{F}^{-1}[\langle y', K_e(\cdot)y \rangle_{\mathcal{Y}}](\omega) = \mathcal{R}\mathcal{F}^{-1}[\langle K_e(\cdot)y', y \rangle_{\mathcal{Y}}](\omega) = \mathcal{R}\langle C(\cdot)y', y \rangle_{\mathcal{Y}}(\omega) = \langle y', C(\omega^{-1})^*y \rangle_{\mathcal{Y}}$. Hence $C(\omega) = C(\omega^{-1})^*$. Suppose that $\mathcal{Y}$ is a complex Hilbert space. Since for all $\omega \in \hat{\mathcal{X}}$, $C(\omega)$ is bounded and non-negative so $C(\omega)$ is self-adjoint. Besides we have $C(\omega) = C(\omega^{-1})^*$ so $C$ must be even. Suppose that $\mathcal{Y}$ is a real Hilbert space. The Fourier Transform of a real valued function obeys $\mathcal{F}[f](\omega) = \overline{\mathcal{F}[f](\omega^{-1})}$. Therefore since $C(\omega)$ is non-negative for all $\omega \in \hat{\mathcal{X}}$, $\langle y', C(\omega)y \rangle = \overline{\langle y', C(\omega^{-1})y \rangle} = \langle y, C(\omega^{-1})^*y' \rangle = \langle y, C(\omega)y' \rangle$. Hence $C(\omega)$ is self-adjoint and thus $C$ is even.

Item 2: simply, for all $y, y' \in \mathcal{Y}$, $\langle y, C(\omega^{-1})y' \rangle = \langle y', C(\omega)y \rangle$ thus $\mathcal{F}^{-1}[\langle y', K_e(\cdot)y \rangle_{\mathcal{Y}}](\omega) = \langle y', C(\omega)y \rangle = \mathcal{R}\langle y', C(\cdot)y \rangle(\omega) = \mathcal{R}\mathcal{F}^{-1}[\langle y', K_e(\cdot)y \rangle_{\mathcal{Y}}](\omega) = \mathcal{F}[\langle y', K_e(\cdot)y \rangle_{\mathcal{Y}}](\omega)$.

Item 3: from Item 2 we have $\mathcal{F}^{-1}[\langle y', K_e(\cdot)y \rangle] = \mathcal{F}^{-1}\mathcal{R}\langle y', K_e(\cdot)y \rangle$. By injectivity of the Fourier Transform, $K_e$ is even. Since $K_e(\delta) = K_e(\delta^{-1})^*$, we must have $K_e(\delta) = K_e(\delta)^*$. ■

### B.1.2 PROOF OF PROPOSITION 17

**Proof** This is a simple consequence of Proposition 15 and Lemma 16. By taking $\langle y', C(\omega)y \rangle = \mathcal{F}^{-1}[\langle y', K_e(\cdot)y \rangle](\omega) = \mathcal{F}[\langle y', K_e(\cdot)y \rangle](\omega)$ we can write the following equality concerning the OVK signature $K_e$: $\langle y', K_e(\delta)y \rangle(\omega) = \int_{\hat{\mathcal{X}}} \overline{(\delta, \omega)}\langle y', C(\omega)y \rangle d\widehat{\text{Haar}}(\omega) = \int_{\hat{\mathcal{X}}} \overline{(\delta, \omega)}\langle y', \frac{1}{\rho(\omega)}C(\omega)y \rangle \rho(\omega)d\widehat{\text{Haar}}(\omega)$. It is always possible to choose $\rho(\omega)$ such that $\int_{\hat{\mathcal{X}}} \rho(\omega)d\widehat{\text{Haar}}(\omega) = 1$. For instance choose

$$ \rho(\omega) = \frac{\|C(\omega)\|_{\mathcal{Y},\mathcal{Y}}}{\int_{\hat{\mathcal{X}}} \|C(\omega)\|_{\mathcal{Y},\mathcal{Y}} d\widehat{\text{Haar}}(\omega)} $$

Since for all $y, y' \in \mathcal{Y}$, $\langle y', C(\cdot)y \rangle \in L^1(\hat{\mathcal{X}}, \widehat{\text{Haar}})$ and $\mathcal{Y}$ is a separable Hilbert space, by Pettis measurability theorem, $\int_{\hat{\mathcal{X}}} \|C(\omega)\|_{\mathcal{Y},\mathcal{Y}} d\widehat{\text{Haar}}(\omega)$ is finite and so is $\|C(\omega)\|_{\mathcal{Y},\mathcal{Y}}$ for all $\omega \in \hat{\mathcal{X}}$. Therefore $\rho(\omega)$ is the density of a probability measure $\text{Pr}_{\widehat{\text{Haar}},\rho}$, i. e. conclude by taking $\text{Pr}_{\widehat{\text{Haar}},\rho}(\mathcal{Z}) = \int_{\mathcal{Z}} \rho(\omega) d\widehat{\text{Haar}}(\omega)$, for all $\mathcal{Z} \in \mathcal{B}(\hat{\mathcal{X}})$. ■

### B.1.3 PROOF OF PROPOSITION 18

**Proof** Suppose that for all $y, y' \in \mathcal{Y}$, $\langle y', A(\omega)y \rangle \rho(\omega) = \mathcal{F}^{-1}[\langle y', K_e(\cdot)y \rangle](\omega)$ where $\rho$ is a probability distribution (see Proposition 17). From the strong law of large numbers $\frac{1}{D} \sum_{j=1}^D \overline{(x \star z^{-1}, \omega_j)} A(\omega_j) \xrightarrow{\text{a. s.}} \mathbf{E}_{\widehat{\text{Haar}},\rho}[(x \star z^{-1}, \omega_j) A(\omega)]$ where the integral converges
in the weak operator topology. Then by Proposition 17 we recover $K_e$ when $D \to \infty$ since, $\mathbf{E}_{\widehat{\text{Haar}},\rho}[(x \star z^{-1}, \omega_j)A(\omega)] = K_e(x \star z^{-1})$. ■

### B.1.4 PROOF OF PROPOSITION 19

**Proof** Let $(\omega_j)_{j=1}^D$ be a sequence of $D \in \mathbb{N}^*$ i.i.d. random variables following the law $\text{Pr}_{\widehat{\text{Haar}},\rho}$. For all $x, z \in \mathcal{X}$ and all $y, y' \in \mathcal{Y}$,

$$ \langle \tilde{\phi}(x)y, \tilde{\phi}(z)y' \rangle_{\bigoplus_{j=1}^D y'} = \frac{1}{D} \left\langle \bigoplus_{j=1}^D ((x, \omega_j)B(\omega_j)^*y), \bigoplus_{j=1}^D ((z, \omega_j)B(\omega_j)^*y') \right\rangle $$

By definition of the inner product in direct sum of Hilbert spaces,

$$ \frac{1}{D} \left\langle \bigoplus_{j=1}^{D} ((x, \omega_j) B(\omega_j)^* y), \bigoplus_{j=1}^{D} ((z, \omega_j) B(\omega_j)^* y') \right\rangle \\ = \frac{1}{D} \sum_{j=1}^{D} \left\langle y, \overline{(x, \omega_j) B(\omega_j)(z, \omega_j) B(\omega_j)^* y'} \right\rangle_y = \left\langle y, \left( \frac{1}{D} \sum_{j=1}^{D} \overline{(x \star z^{-1}, \omega_j) A(\omega_j)} \right) y' \right\rangle_y, $$

Eventually apply Proposition 18 to obtain the convergence of the Monte-Carlo plug-in estimator to the true kernel $K$. ■

### B.1.5 PROOF OF PROPOSITION 21

**Proof** For all $y, y' \in \mathcal{Y}$ and $x, z \in \mathcal{X}$,

$$ \langle y, \phi_x^* \phi_z y' \rangle_y = \langle \phi_x y, \phi_z y' \rangle_{L^2(\hat{\mathcal{X}}, \hat{\mu}; \mathcal{Y})} = \int_{\hat{\mathcal{X}}} \overline{(x, \omega)} \langle y, B(\omega)(z, \omega) B(\omega)^* y' \rangle d\hat{\mu}(\omega) \\ = \int_{\hat{\mathcal{X}}} \overline{(x \star z^{-1}, \omega)} \langle y, B(\omega) B(\omega)^* y' \rangle d\hat{\mu}(\omega) = \int_{\hat{\mathcal{X}}} \overline{(x \star z^{-1}, \omega)} \langle y, A(\omega) y' \rangle d\hat{\mu}(\omega), $$

which defines a $\mathcal{Y}$-Mercer according to Proposition 14 of Carmeli et al. (2010). ■

### B.1.6 PROOF OF PROPOSITION 22

**Proof** From the strong law of large numbers $\frac{1}{D} \sum_{j=1}^D \overline{(x \star z^{-1}, \omega_j) A(\omega_j)} \xrightarrow{a.s.} \mathbf{E}_{\widehat{\text{Haar}},\rho}[(x \star z^{-1}, \omega_j)A(\omega)]$ where the integral converges in the weak operator topology. Then by Proposition 14, $\mathbf{E}_{\widehat{\text{Haar}},\rho}[(x \star z^{-1}, \omega_j)A(\omega)] = K_e(x \star z^{-1})$. ■
### B.1.7 PROOF OF PROPOSITION 23

**Proof** Apply Proposition 6 to $\tilde{\phi}$ considering the Hilbert space $\tilde{\mathcal{H}}$ to show that $\tilde{K}$ is an OVK. Then Proposition 12 shows that $\tilde{K}$ is shift-invariant since $\tilde{K}(x, z) = \tilde{K}_e (x \star z^{-1})$. Since $B(\omega)$ is a bounded operator, $\tilde{K}$ is $\mathcal{Y}$-Mercer because all the functions in the sum are continuous. ■

### B.1.8 PROOF OF PROPOSITION 25

**Proof** [of item 1] Since $(\omega_j)_{j=1}^D$ are i. i. d. random vectors, for all $y \in \mathcal{Y}$ and for all $y' \in \mathcal{Y}'$, $\langle y, B(\cdot)y' \rangle \in L^2(\hat{\mathcal{X}}, \mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho})$ and $g \in L^2(\hat{\mathcal{X}}, \mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho}; \mathcal{Y}')$,

$$ (\tilde{W}\theta)(x) = \tilde{\phi}(x)^*\theta = \frac{1}{D}\sum_{j=1}^D \overline{(x, \omega_j)} B(\omega_j)g(\omega_j), \quad \omega_j \sim \mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho} \quad \text{i. i. d.} \\ \xrightarrow{D \to \infty} \int_{\hat{\mathcal{X}}} \overline{(x, \omega)} B(\omega)g(\omega) d\mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho}(\omega) = (Wg)(x) := \phi_x^*g. $$

from the strong law of large numbers. ■

**Proof** [of item 2] Again, since $(\omega_j)_{j=1}^D$ are i. i. d. random vectors and $g \in L^2(\hat{\mathcal{X}}, \mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho}; \mathcal{Y}')$,

$$ \|\theta\|_{\tilde{\mathcal{H}}}^2 = \frac{1}{D} \sum_{j=1}^{D} \|g(\omega_j)\|_{\mathcal{Y}'}^2, \quad \omega_j \sim \mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho} \quad \text{i. i. d.} \\ \xrightarrow{D \to \infty} \int_{\hat{\mathcal{X}}} \|g(\omega)\|_{\mathcal{Y}'}^2 d\mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho}(\omega) = \|g\|_{L^2(\hat{\mathcal{X}}, \mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho}; \mathcal{Y}')^2. $$

from the strong law of large numbers. ■

### B.1.9 PROOF OF PROPOSITION 26

**Proof** We first show how the Fourier Transform relates to the feature operator. Since $\mathcal{H}_K$ is embedded into $\mathcal{H} = L^2(\hat{\mathcal{X}}, \mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho}; \mathcal{Y}')$ by means of the feature operator $W$, we have for all $f \in \mathcal{H}_k$, for all $f \in \mathcal{H}$ and for all $x \in \mathcal{X}$

$$ \mathcal{F}[\mathcal{F}^{-1}[f]](x) = \int_{\hat{\mathcal{X}}} \overline{(x, \omega)} \mathcal{F}^{-1}[f](\omega) d\widehat{\mathcal{Haar}}(\omega) = f(x) \\ (Wg)(x) = \int_{\hat{\mathcal{X}}} \overline{(x, \omega)} \rho(\omega) B(\omega) g(\omega) d\widehat{\mathcal{Haar}}(\omega) = f(x). $$

By injectivity of the Fourier Transform, $\mathcal{F}^{-1}[f](\omega) = \rho(\omega)B(\omega)g(\omega)$. From Proposition 6 we have

$$ \|f\|_K^2 = \inf \left\{ \|g\|_\mathcal{H}^2 \mid \forall g \in \mathcal{H}, Wg = f \right\} \\ = \inf \left\{ \int_{\hat{\mathcal{X}}} \|g(\omega)\|_{\mathcal{Y}'}^2 d\mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho}(\omega) \mid \forall g \in \mathcal{H}, \mathcal{F}^{-1}[f] = \rho(\cdot)B(\cdot)g(\cdot) \right\}. $$
The pseudo inverse of the operator $B(\omega)$ – noted $B(\omega)^{\dagger}$ – is the unique solution of the system $\mathcal{F}^{-1}[f](\omega) = \rho(\omega)B(\omega)g(\omega)$ w.r.t. $g(\omega)$ with minimal norm<fnref n="13" />. Eventually, $\|f\|_K^2 = \int_{\hat{\chi}} \frac{\|B(\omega)^{\dagger}\mathcal{F}^{-1}[f](\omega)\|_Y^2}{\rho(\omega)^2} d\mathbf{Pr}_{\widehat{\text{Haar}},\rho}(\omega)$ Using the fact that $\mathcal{F}^{-1}[\cdot] = \mathcal{F}\mathcal{R}[\cdot]$ and $\mathcal{F}^2[\cdot] = \mathcal{R}[\cdot]$,

$$ \begin{aligned} \|f\|_K^2 &= \int_{\hat{\chi}} \frac{\|\mathcal{R}[B(\cdot)^{\dagger}\rho(\cdot)](\omega)\mathcal{F}[f](\omega)\|_Y^2}{\rho(\omega)^2} d\widehat{\text{Haar}}(\omega) \\ &= \int_{\hat{\chi}} \frac{\langle B(\omega)^{\dagger}\mathcal{F}[f](\omega), B(\omega)^{\dagger}\mathcal{F}[f](\omega)\rangle_Y}{\rho(\omega)} d\widehat{\text{Haar}}(\omega) \\ &= \int_{\hat{\chi}} \frac{\langle \mathcal{F}[f](\omega), A(\omega)^{\dagger}\mathcal{F}[f](\omega)\rangle_Y}{\rho(\omega)} d\widehat{\text{Haar}}(\omega). \end{aligned} $$

<figure><img src="image_16.png" /></figure>

## B.2 Convergence with high probability of the ORFF estimator

We recall the notations $\delta = x \star z^{-1}$, for all $x, z \in \mathcal{X}$, $\tilde{K}(x, z) = \tilde{\phi}(x)^* \tilde{\phi}(z)$, $\tilde{K}^j(x, z) = \phi_x(\omega_j)^* \phi_z(\omega_j)$, where $\omega_j \sim \mathbf{Pr}_{\widehat{\text{Haar}},\rho}$ and $K_e(\delta) = K(x, z)$ and $\tilde{K}_e(\delta) = \tilde{K}(x, z)$. For the sake of readability, we use throughout the proof the quantities: $F(\delta) := \tilde{K}(x, z) - K(x, z)$ and $F^j(\delta) := \frac{1}{D} (\tilde{K}^j(x, z) - K(x, z))$. We also view $\mathcal{X}$ as a metric space endowed with the distance $d_{\mathcal{X}}: \mathcal{X} \times \mathcal{X} \to \mathbb{R}_+$. Compared to the scalar case, the proof follows the same scheme as the one described in (Rahimi and Recht, 2007; Sutherland and Schneider, 2015), but we consider an operator norm as measure of the error and therefore concentration inequality dealing with these operator norm. The main feature of Proposition 46 is that it covers the case of bounded ORFF as well as unbounded ORFF. In the case of bounded ORFF, a Bernstein inequality for matrix concentration such that the one proved in Mackey et al. (2014, Corollary 5.2) or the formulation of Tropp (2012) recalled in Koltchinskii et al. (2013) is suitable. However some kernels like the curl and the divergence-free kernels do not have obvious bounded $\|F^j\|_{Y,Y}$ but exhibit $F^j$ with subexponential tails. Therefore, we use an operator Bernstein concentration inequality adapted for random matrices with subexponential norms.

### B.2.1 EPSILON-NET

Let $C \subseteq \mathcal{X}$ be a compact subset of $\mathcal{X}$. Let $D_C = \{ x \star z^{-1} | x, z \in C \}$ with diameter at most $2|C|$ where $|C|$ is the diameter of $C$. Since $C$ is supposed compact, so is $D_C$. Since $D_C$ is also a metric space it is well known that a compact metric space is totally bounded. Thus it is possible to find a finite $\epsilon$-net covering $D_C$. We call $T = N(D_C, r)$ the number of closed balls of radius $r$ required to cover $D_C$. For instance if $D_C$ is a subspace finite dimensional Banach space with diameter at most $2|C|$ it is possible to cover the space with at most $T = (4|C|/r)^d$ balls of radius $r$ (see Cucker and Smale (2001, proposition 5)). Let us call $\delta_i, i = 1, \dots, T$ the center of the $i$-th ball, also called anchor of the $\epsilon$-net. Denote $L_F$ the Lipschitz constant

<footnote n="13">13. Note that since $B(\omega)$ is bounded the pseudo inverse of $B(\omega)$ is well defined for $\widehat{\text{Haar}}$-almost all $\omega$.</footnote>
of F. Let $\|\cdot\|_{y,y}$ be the operator norm on $\mathcal{L}(\mathcal{Y})$ (largest eigenvalue). We introduce the following technical lemma.

**Lemma 40**

$\forall \delta \in \mathcal{D}_c$, if

$$L_F \le \frac{\epsilon}{2r} \qquad (34)$$

and

$$\|F(\delta_i)\|_{y,y} \le \frac{\epsilon}{2}, \quad \text{for all } i \in \mathbb{N}_T^* \qquad (35)$$

then $\|F(\delta)\|_{y,y} \le \epsilon$.

**Proof** $\|F(\delta)\|_{y,y} = \|F(\delta) - F(\delta_i) + F(\delta_i)\|_{y,y} \le \|F(\delta) - F(\delta_i)\|_{y,y} + \|F(\delta_i)\|_{y,y}$ for all $0 < i < T$. Using the Lipschitz continuity of F we have $\|F(\delta) - F(\delta_i)\|_{y,y} \le d_X(\delta, \delta_i)L_F \le rL_F$ hence $\|F(\delta)\|_{y,y} \le rL_F + \|F(\delta_i)\|_{y,y} = \frac{r\epsilon}{2r} + \frac{\epsilon}{2} = \epsilon$. ■

To apply the lemma, we must bound the Lipschitz constant of the operator-valued function F (Equation 34) and $\|F(\delta_i)\|_{y,y}$, for all $i = 1, \dots, T$ as well (Equation 35).

## B.2.2 BOUNDING THE LIPSCHITZ CONSTANT

This proof is a slight generalization of Minh (2016) to arbitrary metric spaces. It differs from our first approach (Brault et al., 2016), based on the proof of Sutherland and Schneider (2015) which was only valid for a finite dimensional input space $\mathcal{X}$ and imposed a twice differentiability condition on the considered kernel.

**Lemma 41**

Let $H_\omega \in \mathbb{R}_+$ be the Lipschitz constant of $h_\omega(\cdot)$ and assume that $\int_{\hat{X}} H_\omega \|A(\omega)\|_{y,y} d\mathbf{Pr}_{\widehat{Haar},\rho}(\omega) < \infty$. Then the operator-valued function $K_e: \mathcal{X} \to \mathcal{L}(\mathcal{Y})$ is Lipschitz with

$$\|K_e(x) - K_e(z)\|_{y,y} \le d_X(x, z) \int_{\hat{X}} H_\omega \|A(\omega)\|_{y,y} d\mathbf{Pr}_{\widehat{Haar},\rho}(\omega). \qquad (36)$$

**Proof** We use the fact that the cosine function is Lipschitz with constant 1 and $h_\omega$ Lipschitz with constant $H_\omega$. For all $x, z \in \mathcal{X}$ we have

$$\begin{aligned} \|\tilde{K}_e(x) - K_e(z)\|_{y,y} &= \left\| \int_{\hat{X}} (\cos h_\omega(x) - \cos h_\omega(z)) A(\omega) d\mathbf{Pr}_{\widehat{Haar},\rho} \right\|_{y,y} \\ &\le \int_{\hat{X}} |\cos h_\omega(x) - \cos h_\omega(z)| \|A(\omega)\|_{y,y} d\mathbf{Pr}_{\widehat{Haar},\rho} \\ &\le \int_{\hat{X}} |h_\omega(x) - h_\omega(z)| \|A(\omega)\|_{y,y} d\mathbf{Pr}_{\widehat{Haar},\rho} \\ &\le d_X(x, z) \int_{\hat{X}} H_\omega \|A(\omega)\|_{y,y} d\mathbf{Pr}_{\widehat{Haar},\rho} \end{aligned}$$

■

In the same way, considering $\tilde{K}_e(\delta) = \frac{1}{D} \sum_{j=1}^D \cos h_{\omega_j}(\delta) A(\omega_j)$, where $\omega_j \sim \mathbf{Pr}_{\widehat{Haar},\rho}$, we
can show that $\tilde{K}_e$ is Lipschitz with $\|\tilde{K}_e(x) - \tilde{K}_e(z)\|_{y,y} \le d_X(x, z) \frac{1}{D} \sum_{j=1}^D H_{\omega_j} \|A(\omega_j)\|_{y,y}$.
Combining the Lipschitz continuity of $\tilde{K}_e$ and $\tilde{K}$ (Lemma 41) we obtain

$$ \begin{aligned} \|F(x) - F(z)\|_{y,y} &= \|\tilde{K}_e(x) - \tilde{K}_e(x) - \tilde{K}_e(z) + K_e(z)\|_{y,y} \\ &\le \|\tilde{K}_e(x) - \tilde{K}_e(z)\|_{y,y} + \|K_e(x) - K_e(z)\|_{y,y} \\ &\le d_X(x, z) \left( \int_{\hat{X}} H_\omega \|A(\omega)\|_{y,y} d\mathbf{Pr}_{\widehat{\text{Haar}},\rho} + \frac{1}{D} \sum_{j=1}^D H_{\omega_j} \|A(\omega_j)\|_{y,y} \right) \end{aligned} $$

Taking the expectation yields $\mathbf{E}_{\widehat{\text{Haar}},\rho}[L_F] = 2 \int_{\hat{X}} H_\omega \|A(\omega)\|_{y,y} d\mathbf{Pr}_{\widehat{\text{Haar}},\rho}(\omega)$ Thus by Markov's inequality,

$$ \mathbf{Pr}_{\widehat{\text{Haar}},\rho} \{ (\omega_j)_{j=1}^D | L_F \ge \epsilon \} \le \frac{\mathbf{E}_{\widehat{\text{Haar}},\rho}[L_F]}{\epsilon} \le \frac{2}{\epsilon} \int_{\hat{X}} H_\omega \|A(\omega)\|_{y,y} d\mathbf{Pr}_{\widehat{\text{Haar}},\rho}. \quad (37) $$

### B.2.3 BOUNDING $F$ ON A GIVEN ANCHOR POINT $\delta_i$

To bound $\|F(\delta_i)\|_{y,y}$, Hoeffding inequality devoted to matrix concentration (Mackey et al., 2014) can be applied. We prefer here to turn to tighter and refined inequalities such as Matrix Bernstein inequalities (Sutherland and Schneider (2015) also pointed that for the scalar case). The first non-commutative (matrix) concentration inequalities are due to the pioneer work of Ahlswede and Winter (2002), using bound on the moment generating function. This gave rise to many applications Tropp (2012); Oliveira (2009); Koltchinskii et al. (2013) ranging from analysis of randomized optimization algorithm to analysis of random graphs and generalization bounds usefull in machine learning. The concentration inequatilty of Koltchinskii et al. (2013) we used in our original paper (Brault et al., 2016) has the default to grow linearly with the dimension $p$ of the output space $\mathcal{Y}$. However if the evaluation of the operator-valued kernel at two points yields a low-rank matrix, this bound could be improved since only a few principal dimensions are relevant. Moreover this bound cannot be used when dealing with operator-valued kernel acting on infinite dimensional Hilbert spaces. Recent results of Minsker (2011) consider the notion of intrinsic dimension to avoid this “curse of dimensionality” (see definition 30 for the definition). When $A$ is approximately low-rank (i. e. many eigenvalues are small), or go quickly to zero, the intrinsic dimension can be much lower than the dimensionality. Indeed, $1 \le \text{IntDim}(A) \le \text{Rank}(A) \le \text{dim}(A)$.

**Theorem 42** (Bounded non-commutative Bernstein with intrinsic dimension (Minsker, 2011; Tropp et al., 2015))

Consider a sequence $(X_j)_{j=1}^D$ of $D$ independent Hilbert-Schmidt self-adjoint random operators acting on a separable Hilbert $\mathcal{Y}$ space that satisfy $\mathbf{E}X_j = 0$ for all $j \in \mathbb{N}_D^*$. Suppose that there exist some constant $U \ge 2\|X_j\|_{y,y}$ almost surely for all $j \in \mathbb{N}_D^*$. Define a semi-definite upper bound for the operator-valued variance $V \ge \sum_{j=1}^D \mathbf{E}X_j^2$. Then for all $\epsilon \ge \sqrt{\|V\|_{y,y}} + U/3$,

$$ \mathbf{Pr} \left\{ \left\| \sum_{j=1}^{D} X_j \right\|_{y,y} \ge \epsilon \right\} \le 4 \text{IntDim}(V) \exp(-\psi_{V,U}(\epsilon)) $$
where $\psi_{V,U}(\epsilon) = \frac{\epsilon^2}{2\|V\|_{y,y} + 2U\epsilon/3}$

the concentration inequality is restricted to the case where $\epsilon \ge \sqrt{\|V\|_{y,y}} + U/3$ since the probability is vacuous on the contrary. The assumption that $X_j$'s are Hilbert-Schmidt operators comes from the fact that the product of two such operator yields a trace-class operator, for which the intrinsic dimension is well defined.

However, to cover the general case including unbounded ORFFs like curl and divergence-free ORFFs, we choose a version of Bernstein matrix concentration inequality proposed in Koltchinskii et al. (2013) that allows to consider matrices that are not uniformly bounded but have subexponential tails. In the following we use the notion of Orlicz norm to bound random variable by their tail behavior rather than their value (see definition 29). For the sake of simplicity, we now fix $\psi(t) = \psi_1(t) = \exp(t) - 1$. Although the Orlicz norm should be adapted to the tail of the distribution of the random operator we want to quantify to obtain the sharpest bounds. We also introduce two technical lemmas related to Orlicz norm. The first one relates the $\psi_1$-Orlicz norm to the moment generating function (MGF).

**Lemma 43**

Let $X$ be a random variable with a strictly monotonic moment-generating function. We have $\|X\|_{\psi_1}^{-1} = MGF_{|X|}^{-1}(2)$.

**Proof** We have

$$ \|X\|_{\psi_1} = \inf \{ C > 0 \mid \mathbf{E}[\exp(|X|/C)] \le 2 \} = \frac{1}{\sup \{ C > 0 \mid \text{MGF}_{|X|}(C) \le 2 \}} $$

$X$ has strictly monotonic moment-generating thus $C^{-1} = \text{MGF}_{|X|}^{-1}(2)$. Hence $\|X\|_{\psi_1}^{-1} = \text{MGF}_{|X|}^{-1}(2)$. ■

The second lemma gives the Orlicz norm of a positive constant.

**Lemma 44**

If $a \in \mathbb{R}_+$ then $\|a\|_{\psi_1} = \frac{a}{\ln(2)} < 2a$.

**Proof** We consider $a$ as a positive constant random variable, whose Moment Generating Function (MGF) is $\text{MGF}_a(t) = \exp(at)$. From Lemma 43, $\|a\|_{\psi_1} = \frac{1}{\text{MGF}_X^{-1}(2)}$. Then $\text{MGF}_{|a|}^{-1}(2) = \frac{\ln(2)}{|a|}$, $a \ne 0$. If $a = 0$ then $\|a\|_{\psi_1} = 0$ by definition of a norm. Thus $\|a\|_{\psi_1} = \frac{a}{\ln(2)}$. ■

We now turn our attention to Minsker (2011)'s theorem to for unbounded random variables.

**Theorem 45** (Unbounded non-commutative Bernstein with intrinsic dimension)

Consider a sequence $(X_j)_{j=1}^D$ of $D$ independent self-adjoint random operators acting on a finite dimensional Hilbert space $\mathcal{Y}$ of dimension $p$ that satisfy $\mathbf{EX}_j = 0$ for all $j \in \mathbb{N}_D^*$. Suppose that there exist some constant $U \ge \|\|X_j\|_{\mathcal{Y},\mathcal{Y}}\|$ for all $j \in \mathbb{N}_D^*$. Define a semi-definite upper bound for the operator-valued variance $V \succ \sum_{j=1}^D \mathbf{EX}_j^2$. Then for all
$\epsilon > 0,$

$$
\mathrm{Pr} \left\{ \left\| \sum_{j=1}^{D} X_j \right\|_{y,y} \ge \epsilon \right\} \le \begin{cases} 2 \mathrm{IntDim}(V) \exp \left( -\frac{\epsilon^2}{2\|V\|_{y,y} (1+\frac{1}{p})} \right) r_V(\epsilon), & \epsilon \le \frac{\|V\|_{y,y} (1+1/p)}{2U K(V,p)} \\ 2 \mathrm{IntDim}(V) \exp \left( -\frac{\epsilon}{4UK(V,p)} \right) r_V(\epsilon), & \text{otherwise.} \end{cases}
$$

$$
\text{where } K(V,p) = \log(16\sqrt{2}p) + \log\left(\frac{DU^2}{\|V\|_{y,y}}\right) \text{ and } r_V(\epsilon) = 1 + \frac{3}{\epsilon^2 \log^2(1+\epsilon/\|V\|_{y,y})}
$$

Let $\psi = \psi_1$. To use Theorem 45, we set $X_j = F^j(\delta_i)$. We have indeed $\mathbf{E}_{\widehat{\text{Haar}},\rho}[F^j(\delta_i)] = 0$
since $\tilde{K}(\delta_i)$ is the Monte-Carlo approximation of $K_e(\delta_i)$ and the matrices $F^j(\delta_i)$ are self-
adjoint. We assume we can bound all the Orlicz norms of the $F^j(\delta_i) = \frac{1}{D}(\tilde{K}^j(\delta_i) - K_e(\delta_i))$. In
the following we use constants $u_i$ such that $u_i = DU$. Using Lemma 44 and the sub-additivity
of the $\|\cdot\|_{y,y}$ and $\|\cdot\|_{\psi_1}$ norm,

$$
\begin{align*}
u_i = 2D \max_{1 \le j \le D} \left\| \|F^j(\delta_i)\|_{y,y} \right\|_{\psi_1} &\le 2 \max_{1 \le j \le D} \left\| \|\tilde{K}^j(\delta_i)\|_{y,y} \right\|_{\psi_1} + 2 \left\| \|K_e(\delta_i)\|_{y,y} \right\|_{\psi_1} \\
&< 4 \max_{1 \le j \le D} \left\| \|A(\omega_j)\|_{y,y} \right\|_{\psi_1} + 4 \|K_e(\delta_i)\|_{y,y} = 4 \left( \left\| \|A(\omega)\|_{y,y} \right\|_{\psi_1} + \|K_e(\delta_i)\|_{y,y} \right)
\end{align*}
$$

In the same way we defined the constants $v_i = DV$, $v_i = D \sum_{j=1}^{D} \mathbf{E}_{\widehat{\text{Haar}},\rho} F^j(\delta_i)^2 = D \text{Var}_{\widehat{\text{Haar}},\rho} [\tilde{K}(\delta_i)]$. Then applying Theorem 45, we get for all $i \in \mathbb{N}^*_{\mathcal{N}(\mathcal{D}_C, r)}$ ($i$ is the index of each anchor)

$$
\mathbf{Pr}_{\widehat{\text{Haar}},\rho} \left\{ (\omega_j)_{j=1}^D \middle| \|F(\delta_i)\|_{y,y} \ge \epsilon \right\} \\
\le \begin{cases}
4 \operatorname{IntDim}(v_i) \exp \left( -D \frac{\epsilon^2}{2\|v_i\|_{y,y} (1+\frac{1}{p})} \right) r_{v_i/D}(\epsilon), & \epsilon \le \frac{\|v_i\|_{y,y} (1+1/p)}{2u_i K(v_i,p)} \\
4 \operatorname{IntDim}(v_i) \exp \left( -D \frac{\epsilon}{4u_i K(v_i,p)} \right) r_{v_i/D}(\epsilon), & \text{otherwise.}
\end{cases}
$$

with

$$
K(v_i, p) = \log(16\sqrt{2}p) + \log\left(\frac{u_i^2}{\|v_i\|_{y,y}}\right)
$$

and

$$
r_{v_i/D} = 1 + \frac{3}{\epsilon^2 \log^2(1 + D\epsilon / \|v_i\|_{y,y})}.
$$

To unify the bound on each anchor we define two constant

$$
u = 4 \left( \left\| \|A(\omega)\|_{y,y} \right\|_{\psi_1} + \sup_{\delta \in \mathcal{D}_C} \|K_e(\delta)\|_{y,y} \right) \ge \max_{i=1,...T} u_i
$$

and

$$
v = \sup_{\delta \in \mathcal{D}_C} D \mathrm{Var}_{\widehat{\text{Haar}}, \rho} [\tilde{K}_e(\delta)] \geq \max_{i=1, \dots, T} v_i.
$$
### B.2.4 UNION BOUND AND EXAMPLES

Taking the union bound over the anchors yields

$$ \mathbf{Pr}_{\widehat{\text{Haar}},\rho} \left\{ (\omega_j)_{j=1}^D \left| \bigcup_{i=1}^{\mathcal{N}(\mathcal{D}_C, r)} \|F(\delta_i)\|_{\mathcal{Y},\mathcal{Y}} \ge \epsilon \right. \right\} \\ \le 4\mathcal{N}(\mathcal{D}_C, r) r_{v/D}(\epsilon) \text{IntDim}(v) \left\{ \begin{array}{ll} \exp\left(-D\frac{\epsilon^2}{2\|v\|_{\mathcal{Y},\mathcal{Y}}(1+\frac{1}{p})}\right), & \epsilon \le \frac{\|v\|_{\mathcal{Y},\mathcal{Y}}}{2u} \frac{1+1/p}{K(v,p)} \\ \exp\left(-D\frac{\epsilon}{4uK(v,p)}\right), & \text{otherwise.} \end{array} \right. \quad (38) $$

Hence combining Equation 37 and Equation 38 gives and summing up the hypothesis yields the following proposition

**Proposition 46**

Let $K: \mathcal{X} \times \mathcal{X} \to \mathcal{L}(\mathcal{Y})$ be a shift-invariant $\mathcal{Y}$-Mercer kernel, where $\mathcal{Y}$ is a finite dimensional Hilbert space of dimension $p$ and $\mathcal{X}$ a metric space. Moreover, let $C$ be a compact subset of $\mathcal{X}$, $A: \hat{\mathcal{X}} \to \mathcal{L}(\mathcal{Y})$ and $\mathbf{Pr}_{\widehat{\text{Haar}},\rho}$ a pair such that $\tilde{K}_e = \sum_{j=1}^D \cos(\cdot, \omega_j)A(\omega_j) \approx K_e$ $\omega_j \sim \mathbf{Pr}_{\widehat{\text{Haar}},\rho}$ i. i. d.. Let $V(\delta) \succcurlyeq \mathbf{Var}_{\widehat{\text{Haar}},\rho} \tilde{K}_e(\delta)$, for all $\delta \in \mathcal{D}_C$ and $H_\omega$ be the Lipschitz constant of the function $h: x \mapsto (x, \omega)$. If the three following constant exists

$$ m \ge \int_{\hat{\mathcal{X}}} H_{\omega} \|A(\omega)\|_{\mathcal{Y},\mathcal{Y}} d\mathbf{Pr}_{\widehat{\text{Haar}},\rho} < \infty $$

and

$$ u \ge 4 \left( \left\| \|A(\omega)\|_{\mathcal{Y},\mathcal{Y}} \right\|_{\psi_1} + \sup_{\delta \in \mathcal{D}_C} \|K_e(\delta)\|_{\mathcal{Y},\mathcal{Y}} \right) < \infty $$

and

$$ v \ge \sup_{\delta \in \mathcal{D}_C} D \|V(\delta)\|_{\mathcal{Y},\mathcal{Y}} < \infty. $$

Define $p_{\text{int}} \ge \sup_{\delta \in \mathcal{D}_C} \text{IntDim}(V(\delta))$ then for all $r \in \mathbb{R}_+^*$ and all $\epsilon \in \mathbb{R}_+^*$,

$$ \mathbf{Pr}_{\widehat{\text{Haar}},\rho} \left\{ (\omega_j)_{j=1}^D \left| \left\| \tilde{K} - K \right\|_{\mathcal{C} \times \mathcal{C}} \ge \epsilon \right. \right\} \le 4 \left( \begin{array}{ll} \frac{rm}{\epsilon} + p_{\text{int}} \mathcal{N}(\mathcal{D}_C, r) r_{v/D}(\epsilon) & \\ \left\{ \begin{array}{ll} \exp\left(-D\frac{\epsilon^2}{8v(1+\frac{1}{p})}\right), & \epsilon \le \frac{v}{u} \frac{1+1/p}{K(v,p)} \\ \exp\left(-D\frac{\epsilon}{8uK(v,p)}\right), & \text{otherwise.} \end{array} \right. & \end{array} \right) $$

where

$$ K(v,p) = \log(16\sqrt{2p}) + \log\left(\frac{u^2}{\|v\|_{\mathcal{Y},\mathcal{Y}}}\right) $$

and

$$ r_{v/D}(\epsilon) = 1 + \frac{3}{\epsilon^2 \log^2(1 + D\epsilon/\|v\|_{\mathcal{Y},\mathcal{Y}})} $$
**Proof** Let $m = \int_{\hat{\chi}} H_{\omega} \|A(\omega)\|_{y,y} d\mathbf{Pr}_{\widehat{\text{Haar}},\rho}$. From Lemma 41, $\mathbf{Pr}_{\widehat{\text{Haar}},\rho} \{ (\omega_j)_{j=1}^D \mid L_F \ge \frac{\epsilon}{2r} \} \le \frac{4rm}{\epsilon}$. Thus from Lemma 40, for all $r \in \mathbb{R}_+^*$,

$$ \mathbf{Pr}_{\widehat{\text{Haar}},\rho} \left\{ (\omega_j)_{j=1}^D \left| \sup_{\delta \in \mathcal{D}_c} \|F(\delta)\|_{y,y} \ge \epsilon \right. \right\} \le \\ \mathbf{Pr}_{\widehat{\text{Haar}},\rho} \left\{ (\omega_j)_{j=1}^D \left| L_F \ge \frac{\epsilon}{2r} \right. \right\} + \mathbf{Pr}_{\widehat{\text{Haar}},\rho} \left\{ (\omega_j)_{j=1}^D \left| \bigcup_{i=1}^{\mathcal{N}(\mathcal{D}_c, r)} \|F(\delta_i)\|_{y,y} \ge \epsilon \right. \right\} = 4 \frac{rm}{\epsilon} \\ + 4\mathcal{N}(\mathcal{D}_c, r) r_{v/D}(\epsilon) \text{IntDim}(v) \begin{cases} \exp\left(-D \frac{\epsilon^2}{8\|v\|_{y,y} \left(1+\frac{1}{p}\right)}\right), & \epsilon \le \frac{\|v\|_{y,y} \left(1+\frac{1}{p}\right)}{u K(v,p)} \\ \exp\left(-D \frac{\epsilon}{8uK(v,p)}\right), & \text{otherwise.} \end{cases} $$

<figure><img src="image_17.png" /></figure>

With minor modifications we can obtain a second inequality for the case where the random operators $A(\omega_j)$ are bounded almost surely. This second bound with more restrictions on $A$ has the advantage of working in infinite dimension as long as $A(\omega_j)$ is a Hilbert-Schmidt operator.

**Proposition 47**

Let $K: \mathcal{X} \times \mathcal{X} \to \mathcal{L}(\mathcal{Y})$ be a shift-invariant $\mathcal{Y}$-Mercer kernel, where $\mathcal{Y}$ is a Hilbert space and $\mathcal{X}$ a metric space. Moreover, let $C$ be a compact subset of $\mathcal{X}$, $A: \hat{\chi} \to \mathcal{L}(\mathcal{Y})$ and $\mathbf{Pr}_{\widehat{\text{Haar}},\rho}$ a pair such that $\tilde{K}_e = \sum_{j=1}^D \cos(\cdot, \omega_j) A(\omega_j) \approx K_e$, $\omega_j \sim \mathbf{Pr}_{\widehat{\text{Haar}},\rho}$ i. i. d.. where $A(\omega_j)$ is a Hilbert-Schmidt operator for all $j \in \mathbb{N}_D^*$. Let $\mathcal{D}_c = \mathcal{C} \star \mathcal{C}^{-1}$ and $V(\delta) \succ \mathbf{Var}_{\widehat{\text{Haar}},\rho} \tilde{K}_e(\delta)$, for all $\delta \in \mathcal{D}_c$ and $H_\omega$ be the Lipschitz constant of the function $h: x \mapsto (x, \omega)$. If the three following constant exists

$$ m \ge \int_{\hat{\chi}} H_{\omega} \|A(\omega)\|_{y,y} d\mathbf{Pr}_{\widehat{\text{Haar}},\rho} < \infty $$

and

$$ u \ge \underset{\omega \in \hat{\chi}}{\operatorname{ess\,sup}} \|A(\omega)\|_{y,y} + \sup_{\delta \in \mathcal{D}_c} \|K_e(\delta)\|_{y,y} < \infty $$

and

$$ v \ge \sup_{\delta \in \mathcal{D}_c} D \|V(\delta)\|_{y,y} < \infty. $$

define $p_{\text{int}} \ge \sup_{\delta \in \mathcal{D}_c} \text{IntDim}(V(\delta))$ then for all $r \in \mathbb{R}_+^*$ and all $\epsilon > \sqrt{\frac{v}{D}} + \frac{1}{3D}u$,

$$ \mathbf{Pr}_{\widehat{\text{Haar}},\rho} \left\{ (\omega_j)_{j=1}^D \left| \sup_{\delta \in \mathcal{D}_c} \|F(\delta)\|_{y,y} \ge \epsilon \right. \right\} \le 4 \left( \frac{rm}{\epsilon} + p_{\text{int}} \mathcal{N}(\mathcal{D}_c, r) \exp(-D\psi_{v,u}(\epsilon)) \right) $$

where $\psi_{v,u}(\epsilon) = \frac{\epsilon^2}{2(v+u\epsilon/3)}$.
When the covering number $\mathcal{N}(\mathcal{D}_C, r)$ of the metric space $\mathcal{D}_C$ has an analytical form, it is possible to optimize the bound over the radius $r$ of the covering balls. As an example, we refine Proposition 46 and Proposition 47 in the case where $\mathcal{C}$ is a finite dimensional Banach space.

**Corollary 48**

Let $K: \mathcal{X} \times \mathcal{X} \to \mathcal{L}(\mathcal{Y})$ be a shift-invariant $\mathcal{Y}$-Mercer kernel, where $\mathcal{Y}$ is a finite dimensional Hilbert space of dimension $p$ and $\mathcal{X}$ a finite dimensional Banach space of dimension $d$. Moreover, let $\mathcal{C}$ be a closed ball of $\mathcal{X}$ centered at the origin of diameter $|\mathcal{C}|$, $A: \tilde{\mathcal{X}} \to \mathcal{L}(\mathcal{Y})$ and $\mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho}$ a pair such that $\tilde{K}_e = \sum_{j=1}^D \cos(\cdot, \omega_j)A(\omega_j) \approx K_e$, $\omega_j \sim \mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho}$ i. i. d.. Let $\mathcal{D}_C = \mathcal{C} \star \mathcal{C}^{-1}$ and $V(\delta) \ge \mathbf{Var}_{\widehat{\mathcal{Haar}},\rho} \tilde{K}_e(\delta)$, for all $\delta \in \mathcal{D}_C$ Let $H_\omega$ be the Lipschitz constant of $h_\omega: x \mapsto (x, \omega)$. If the three following constant exists

$$m \ge \int_{\tilde{\mathcal{X}}} H_{\omega} \|A(\omega)\|_{\mathcal{Y},\mathcal{Y}} d\mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho} < \infty$$

and

$$u \ge 4 \left( \left\| \|A(\omega)\|_{\mathcal{Y},\mathcal{Y}} \right\|_{\psi_1} + \sup_{\delta \in \mathcal{D}_C} \|K_e(\delta)\|_{\mathcal{Y},\mathcal{Y}} \right) < \infty$$

and

$$v \ge \sup_{\delta \in \mathcal{D}_C} D \|V(\delta)\|_{\mathcal{Y},\mathcal{Y}} < \infty.$$

Define $p_{int} \ge \sup_{\delta \in \mathcal{D}_C} \mathrm{IntDim}(V(\delta))$, then for all $0 < \epsilon \le m|\mathcal{C}|$,

$$\mathbf{Pr}_{\widehat{\mathcal{Haar}},\rho} \left\{ (\omega_j)_{j=1}^D \left| \left\| \tilde{K} - K \right\|_{\mathcal{C} \times \mathcal{C}} \ge \epsilon \right. \right\} \\ \le 8\sqrt{2} \left( \frac{m|\mathcal{C}|}{\epsilon} \right) (p_{int} r_{v/D}(\epsilon))^{\frac{1}{d+1}} \left\{ \begin{array}{ll} \exp \left( -D \frac{\epsilon^2}{8v(d+1)(1+\frac{1}{p})} \right), & \epsilon \le \frac{v}{u} \frac{1+1/p}{K(v,p)} \\ \exp \left( -D \frac{\epsilon}{8u(d+1)K(v,p)} \right), & \text{otherwise,} \end{array} \right.$$

$$\text{where } K(v,p) = \log(16\sqrt{2}p) + \log\left(\frac{u^2}{v}\right) \text{ and } r_{v/D}(\epsilon) = 1 + \frac{3}{\epsilon^2 \log^2(1+D\epsilon/v)}.$$

**Proof** As we have seen in appendix B.2.1, suppose that $\mathcal{X}$ is a finite dimensional Banach space. Let $\mathcal{C} \subset \mathcal{X}$ be a closed ball centered at the origin of diameter $|\mathcal{C}| = \mathcal{C}$ then the difference ball centered at the origin

$$\mathcal{D}_C = \mathcal{C} \star \mathcal{C}^{-1} = \{ x \star z^{-1} \mid \|x\|_{\mathcal{X}} \le C/2, \|z\|_{\mathcal{X}} \le C/2, (x, z) \in \mathcal{X}^2 \} \subset \mathcal{X}$$
is closed and bounded, so compact and has diameter $|C| = 2C$. It is possible to cover it with $\log(\mathcal{N}(\mathcal{D}_C, r)) = d \log(\frac{2|C|}{r})$ closed balls of radius $r$. Plugging back into Equation yields

$$ \mathbf{Pr}_{\widehat{\text{Haar}},\rho} \left\{ (\omega_j)_{j=1}^D \left| \left\| \tilde{K} - K \right\|_{\mathcal{C} \times \mathcal{C}} \ge \epsilon \right. \right\} \le 4 \left( \frac{rm}{\epsilon} + p_{\text{int}} \left( \frac{2|C|}{r} \right)^d r_{v/D}(\epsilon) \right. \\ \left. \left\{ \begin{array}{ll} \exp \left( -D \frac{\epsilon^2}{8v \left(1+\frac{1}{p}\right)} \right), & \epsilon \le \frac{v}{u} \frac{1+1/p}{K(v,p)} \\ \exp \left( -D \frac{\epsilon}{8uK(v,p)} \right), & \text{otherwise.} \end{array} \right. \right) $$

The right hand side of the equation has the form $ar + br^{-d}$ with $a = \frac{m}{\epsilon}$ and

$$ b = p_{\text{int}}(2|C|)^d r_{v/D}(\epsilon) \left\{ \begin{array}{ll} \exp\left(-D\frac{\epsilon^2}{8v\left(1+\frac{1}{p}\right)}\right), & \epsilon \le \frac{v}{u}\frac{1+1/p}{K(v,p)} \\ \exp\left(-D\frac{\epsilon}{8uK(v,p)}\right), & \text{otherwise.} \end{array} \right. $$

Following Rahimi and Recht (2007); Sutherland and Schneider (2015); Minh (2016), we optimize over $r$. It is a convex continuous function on $\mathbb{R}_+$ and achieve minimum at $r = \left(\frac{bd}{a}\right)^{\frac{1}{d+1}}$ and the minimum value is $r_* = a^{\frac{d}{d+1}} b^{\frac{1}{d+1}} \left(d^{\frac{1}{d+1}} + d^{-\frac{d}{d+1}}\right)$, hence

$$ \mathbf{Pr}_{\widehat{\text{Haar}},\rho} \left\{ (\omega_j)_{j=1}^D \left| \left\| \tilde{K} - K \right\|_{\mathcal{C} \times \mathcal{C}} \ge \epsilon \right. \right\} \\ \le C_d \left( \frac{2m|C|}{\epsilon} \right)^{\frac{d}{d+1}} (p_{\text{int}} r_{v/D}(\epsilon))^{\frac{1}{d+1}} \left\{ \begin{array}{ll} \exp\left(-D\frac{\epsilon^2}{8v(d+1)\left(1+\frac{1}{p}\right)}\right), & \epsilon \le \frac{v}{u}\frac{1+1/p}{K(v,p)} \\ \exp\left(-D\frac{\epsilon}{8u(d+1)K(v,p)}\right), & \text{otherwise,} \end{array} \right. \\ \le 8\sqrt{2} \left( \frac{m|C|}{\epsilon} \right) (p_{\text{int}} r_{v/D}(\epsilon))^{\frac{1}{d+1}} \left\{ \begin{array}{ll} \exp\left(-D\frac{\epsilon^2}{8v(d+1)\left(1+\frac{1}{p}\right)}\right), & \epsilon \le \frac{v}{u}\frac{1+1/p}{K(v,p)} \\ \exp\left(-D\frac{\epsilon}{8u(d+1)K(v,p)}\right), & \text{otherwise,} \end{array} \right. $$

where $C_d = 4 \left(d^{\frac{1}{d+1}} + d^{-\frac{d}{d+1}}\right)$. Eventually when $\mathcal{X}$ is a Banach space, the Lipschitz constant of $h_\omega$ is the supremum of the gradient $H_\omega = \sup_{\delta \in \mathcal{D}_C} \|(\nabla h_\omega)(\delta)\|_{\hat{\mathcal{X}}}$. ■

Following the same proof technique we obtain the second bound for bounded ORFF.

### Corollary 49

Let $K: \mathcal{X} \times \mathcal{X} \to \mathcal{L}(\mathcal{Y})$ be a shift-invariant $\mathcal{Y}$-Mercer kernel, where $\mathcal{Y}$ is a Hilbert space and $\mathcal{X}$ a finite dimensional Banach space of dimension $D$. Moreover, let $C$ be a closed ball of $\mathcal{X}$ centered at the origin of diameter $|C|$, subset of $\mathcal{X}$, $A: \hat{\mathcal{X}} \to \mathcal{L}(\mathcal{Y})$ and $\mathbf{Pr}_{\widehat{\text{Haar}},\rho}$ a pair such that $\tilde{K}_e = \sum_{j=1}^D \cos(\cdot, \omega_j)A(\omega_j) \approx K_e$, $\omega_j \sim \mathbf{Pr}_{\widehat{\text{Haar}},\rho}$ i. i. d.. where $A(\omega_j)$ is a Hilbert-Schmidt operator for all $j \in \mathbb{N}_D^*$. Let $\mathcal{D}_C = C \star C^{-1}$ and $V(\delta) \ge \mathbf{Var}_{\widehat{\text{Haar}},\rho} \tilde{K}_e(\delta)$ for
all $\delta \in \mathcal{D}_C$ and $H_\omega$ be the Lipschitz constant of the function $h : x \mapsto (x, \omega)$. If the three following constant exists

$$m \geq \int_{\hat{\chi}} H_{\omega} \|A(\omega)\|_{y,y} d\mathbf{Pr}_{\widehat{\text{Haar}},\rho} < \infty$$

and

$$u \geq \sup_{\omega \in \hat{\chi}} \|A(\omega)\|_{y,y} + \sup_{\delta \in \mathcal{D}_C} \|K_e(\delta)\|_{y,y} < \infty$$

and

$$v \geq \sup_{\delta \in \mathcal{D}_C} D \|V(\delta)\|_{y,y} < \infty.$$

define $p_{\text{int}} \geq \sup_{\delta \in \mathcal{D}_C} \text{IntDim}(V(\delta))$ then for all $\sqrt{\frac{v}{D}} + \frac{u}{3D} < \epsilon < m|C|$,

$$\mathbf{Pr}_{\widehat{\text{Haar}},\rho} \left\{ (\omega_j)_{j=1}^D \left| \sup_{\delta \in \mathcal{D}_C} \|F(\delta)\|_{y,y} \ge \epsilon \right. \right\} \le 8\sqrt{2} \left( \frac{m|C|}{\epsilon} \right) p_{\text{int}}^{\frac{1}{d+1}} \exp(-D\psi_{v,d,u}(\epsilon))$$

$$\text{where } \psi_{v,d,u}(\epsilon) = \frac{\epsilon^2}{2(d+1)(v+u\epsilon/3)}.$$

## B.2.5 PROOF OF THE ORFF ESTIMATOR VARIANCE BOUND (PROPOSITION 33).

We use the notations $\delta = x \star z^{-1}$ for all $x, z \in \mathcal{X}$, $\tilde{K}(x, z) = \tilde{\phi}(x)^* \tilde{\phi}(z)$, $\tilde{K}^j(x, z) = \phi_x(\omega_j)^* \phi_z(\omega_j)$ and $K_e(\delta) = K_e(x, z)$.

**Proof** Let $\delta \in \mathcal{D}_C$ be a constant. From the definition of the variance of a random variable and using the fact that the $(\omega_j)_{j=1}^D$ are i.i.d. random variables,

$$\begin{aligned} \mathrm{Var}_{\widehat{\text{Haar}},\rho} [\tilde{K}_e(\delta)] &= \mathbf{E}_{\widehat{\text{Haar}},\rho} \left[ \frac{1}{D} \sum_{j=1}^D \tilde{K}_e^j(\delta) - K_e(\delta) \right]^2 = \frac{1}{D^2} \mathbf{E}_{\widehat{\text{Haar}},\rho} \left[ \sum_{j=1}^D \tilde{K}_e^j(\delta) - K_e(\delta) \right]^2 \\ &= \frac{1}{D} \mathbf{E}_{\widehat{\text{Haar}},\rho} \left[ \tilde{K}_e^j(\delta)^2 - \tilde{K}_e^j(\delta) K_e(\delta) - K_e(\delta) \tilde{K}_e^j(\delta) + K_e(\delta)^2 \right] \end{aligned}$$

From the definition of $\tilde{K}_e^j$, $\mathbf{E}_{\widehat{\text{Haar}},\rho} \tilde{K}_e^j(\delta) = K_e(\delta)$, which leads to

$$\mathrm{Var}_{\widehat{\text{Haar}},\rho} [\tilde{K}_e(\delta)] = \frac{1}{D} \mathbf{E}_{\widehat{\text{Haar}},\rho} [\tilde{K}_e^j(\delta)^2 - K_e(\delta)^2]$$

A trigonometric identity gives us $(\cos(\delta, \omega))^2 = \frac{1}{2}(\cos(2\delta, \omega) + \cos(e, \omega))$. Thus

$$\mathrm{Var}_{\widehat{\text{Haar}},\rho} [\tilde{K}_e(\delta)] = \frac{1}{2D} \mathbf{E}_{\widehat{\text{Haar}},\rho} [ (\cos(2\delta, \omega) + \cos(e, \omega)) A(\omega)^2 - 2K_e(\delta)^2 ].$$

Also,

$$\begin{aligned} \mathbf{E}_{\widehat{\text{Haar}},\rho} [\cos(2\delta, \omega) A(\omega)^2] &= \mathbf{E}_{\widehat{\text{Haar}},\rho} [\cos(2\delta, \omega) A(\omega)] \mathbf{E}_{\widehat{\text{Haar}},\rho} [A(\omega)] \\ &\quad + \mathrm{Cov}_{\widehat{\text{Haar}},\rho} [\cos(2\delta, \omega) A(\omega), A(\omega)] \\ &= K_e(2\delta) \mathbf{E}_{\widehat{\text{Haar}},\rho} [A(\omega)] + \mathrm{Cov}_{\widehat{\text{Haar}},\rho} [\cos(2\delta, \omega) A(\omega), A(\omega)] \end{aligned}$$