Year: 2019

Abstract:
Distributional robust risk (DRR) minimisation has arisen as a flexible and effective framework for machine learning. Approximate solutions based on dualisation have become particularly favorable in addressing the semi-infinite optimisation, and they also provide a certificate of the robustness for the worst-case population loss. However existing methods are restricted to either linear models or very small perturbations, and cannot find the globally optimal solution for restricted nonlinear models such as kernel machines. In this paper we resolve these limitations for a general class of kernel space, and our approach is based on a new upper bound of DRRs using an empirical risk regularised by the Lipschitz constant of the model, e.g., deep neural networks and kernel methods. As an application, we showed that it also provides a certificate for adversarial training, and global solutions can be achieved on product kernel machines in polynomial time.

heading: 1 INTRODUCTION
text: Regularised risk minimisation has been the workhorse of learning nonlinear hypotheses such as deep neural networks and kernel machines. Recently, distributional robust risk (DRR) minimization has emerged as a promising instance with marked efficacy and flexibility. Instead of perturbing the observed data points, DRRs consider perturbations to the empirical distribution, constituting an ambiguity set P that lives in the space of data distributions. Let Ω be an outcome space with (true) distribution µ, e.g., the joint space of input and output. Given a loss function `, a model f suffers a loss value `f (ω) over an outcome ω, and the risk of f under µ is risk`(f, µ) := Eµ[`f ]. In DRR minimisation, a model f is sought that minimises the expectation of loss ` over an ambiguity set P , i.e., that minimises supν∈P risk`(f, ν) (Delage & Ye, 2010; Goh & Sim, 2010; Wiesemann et al., 2014). The ambiguity sets can be constructed by moment matching (Bhattacharyya et al., 2005; Farnia & Tse, 2016), divergence balls (Ben-Tal et al., 2013; Duchi et al., 2016; Hu & Hong, 2016), or Wasserstein distance balls (Kantorovitch, 1958). In this work we focus on the last due to its favorable properties in statistics and computation, along with extensive applications in DRR (Esfahani & Kuhn, 2018; Gao & Kleywegt, 2016; Zhao & Guan, 2018).
Despite the generality of DRR, its computational efficiency remains a challenge, since the supremum is over a (typically) uncountably infinite dimension space. Tractable equivalent convex programs can be derived only for a limited range of loss functions along with linear hypothesis spaces (Blanchet et al., 2016; El Ghaoui & Lebret, 1997; Shafieezadeh-Abadeh et al., 2015; 2017; Xu et al., 2009a;b). Although Shafieezadeh-Abadeh et al. (2017) developed lifted variants for reproducing kernel Hilbert spaces (RKHS) to accommodate nonlinear hypotheses, the perturbation was applied to Φ(ω), where Φ is the implicit feature map. This still falls short of robustness with respect to distributions over Ω.
A more promising technique for optimizing DRRs over nonlinear hypothesis spaces—including deep neural networks and kernel machines—is by dualising it to a form that is amenable to (approximate) optimisation (Esfahani & Kuhn, 2018). The fundamental strong duality result was established independently by Blanchet & Murthy (2019) and Gao & Kleywegt (2016), and has been applied to various tasks such as specification of regularisation parameter (Blanchet et al., 2016), design of transport cost (Blanchet et al., 2017), and selection of ambiguity region size for optimal confidence interval (Duchi et al., 2016). In particular, Sinha et al. (2018) used it to construct an efficiently computable certificate on the level of robustness for the worst-case population loss. However, these methods are still subject to marked restrictions when applied to smooth nonlinear models. First, Sinha et al. (2018) restricts the perturbation to be small, which despite the common interest of
imperceptible perturbations, leaves unaddressed the equally interesting regimes of medium to large perturbations; see discussions in Openreview (2018). Moreover, although the global solvability of the inner Lagrangian penalty problem (robust surrogate loss) can be ensured by small enough perturbation, there is no practical procedure to compute the threshold. Finally, the overall optimisation of the nonlinear model is still subject to nonconvexity, precluding tractable global solutions for restricted but still general classes of nonlinear models such as kernel methods.
The first goal of this work, therefore, is to develop a novel certificate on distributional robustness that dispenses with these restrictions (§3). Specifically, we will leverage the McShane-Whitney extension theorem (McShane, 1934; Whitney, 1934) to upper bound DRRs by the empirical risk regularized with the Lipschitz constant of the model f , while additionally accounting for the underlying transport cost and the loss `. The result vastly generalises the vector norm regularisation in linear binary classification (Shafieezadeh-Abadeh et al., 2017, Thm. 3.11) to nonlinear models and extended real-valued cost functions that encode constraints, along with an arbitrary metric space of labels that is general enough for multiclass problems. Appealing to any magnitude of perturbation, it also enjoys improved computational efficiency compared with the robust surrogate loss in Sinha et al. (2018).
A particularly effective domain to apply this new certificate is adversarial learning (Szegedy et al., 2014), where models are trained to be resilient to malicious distortions on the data. Although Lipschitz regularisation has been a popular recipe for robustness (e.g., Anil et al., 2019; Cisse et al., 2017; Farnia et al., 2019; Gouk et al., 2018; Huster et al., 2018; Scaman & Virmaux, 2018) and generalisation accuracy (Miyato et al., 2018; Yoshida & Miyato, 2017), it remains a heuristic and therefore our second major contribution is to reveal in §3.1 that adversarial risks (Goodfellow et al., 2015; Madry et al., 2018; Shaham et al., 2018) can be bounded by a DRR. Such a rigorous justification has hitherto been restricted to logistic loss (Suggala et al., 2019, Thm. 9), and a similar tightness result has been established only for linear models (Shafieezadeh-Abadeh et al., 2017, Thm. 3.20). As a result, our new certificate amounts to a new bound on the worst-case risk under attacks, complementing the existing certificates (Raghunathan et al., 2018; Tsuzuku et al., 2018; Weng et al., 2018; Wong & Kolter, 2018; Wong et al., 2018) with a more computationally efficient approach. It further achieved state-of-the-art accuracy under a range of attacks on standard benchmark datasets (§5).
In practice, however, the evaluation of Lipschitz constant L is NP-hard for neural networks (Scaman & Virmaux, 2018), compelling approximations of it, or explicit engineering of layers to respect Lipschitz while analyzing the expressiveness in specific cases (e.g., `∞ norm in Anil et al. (2019)). We, instead, pursue a new path and explore the following question: does there exist a hypothesis space which: a) is expressive enough in modeling; b) allows the exact value of L to be computed efficiently; c) enforcing the Lipschitz constant leads to a convex constraint that renders efficient optimisation.
Interestingly, kernel machines satisfy all these requirements for some kernels. For example, Gaussian kernels are universal, whose RKHS can approximate any continuous function on a compact set in a uniform sense (Micchelli et al., 2006). The RKHS of multi-layer inverse kernels compactly encompasses `1-regularized neural networks (Shalev-Shwartz et al., 2011), degrading the generalisation performance by only a polynomial constant (Zhang et al., 2016; 2017). Similar results have been conjectured for Gaussian kernels (Shalev-Shwartz et al., 2011). Our third contribution proves that b) can be achieved for product kernels such as Gaussian kernels with high probability by using the Nyström approximation (Drineas & Mahoney, 2005; Williams & Seeger, 2000), and ε approximation error of L requires only O(1/ε2) samples (§4). Empirically this approximation is also effective for non-product kernels like inverse kernels. Such a sampling based approach also leads to a single convex constraint, making it scalable to 60k examples with even an interior-point solver (§5). The convenience in evaluating L renders our certificate of DRR even more favorable than those based on robust surrogate losses (Blanchet & Murthy, 2019; Gao & Kleywegt, 2016; Sinha et al., 2018).

heading: 2 PRELIMINARIES
text: The vast majority of our technical results and proofs are deferred to Appendix A, for which theoremlike statements are numbered to be consistent. The extended real line is R̄ := [−∞,+∞], R̄≥0 := [0,∞], and [n] := {1, 2, . . . , n}. For topological spaces X,Y , the Borel subsets are B(X), and the Borel probability measures are P(X). The universal sigma algebra is U(X) := ⋂ µ∈P(X) Bµ(X) where Bµ(X) is the completion of the Borel sets with respect to µ ∈ P(X). Let L0(X,Y ) denote the Borel measurable mappings X → Y , and L1(X,µ) denote the Borel functions f ∈ L0(X,R) with ∫ |f |dµ <∞ for µ ∈ P(X).
For two measures µ, ν ∈ P(Ω) the set of (µ, ν)-couplings is Π(µ, ν) := {π ∈ P(Ω ×Ω) | µ = ∫ π( · ,dω), ν = ∫ π(dω, · )}.
Let c ∈ L0(Ω ×Ω, R̄). The c-transportation cost of µ, ν ∈ P(Ω), and c-transportation cost ball of radius r ≥ 0, centred at µ ∈ P(Ω) are respectively
costc(µ, ν) := inf { ∫ cdπ : π ∈ Π(µ, ν) } and Bc(µ, r) := {ν ∈ P(Ω) | costc(µ, ν) ≤ r}. (1)
A function f : Ω → R̄ is c-Lipschitz if there exists L ≥ 0 such that ∀ω1, ω2 ∈ dom f : |f(ω1)− f(ω2)| ≤ Lc(ω1, ω2). (2)
The least c-Lipschitz constant of f (cf. Cranko et al., 2019) is the infimum over L ≥ 0 satisfying (2), and is denoted by lipc(f), so that when (X, d) is a metric space lipd(f) agrees with the usual Lipschitz notion. When c : X → R̄ (e.g., when c is a norm), we take c(x, y) := c(x − y) for all x, y ∈ X in (1) and (2).

heading: 3 CERTIFICATE FOR DISTRIBUTIONAL ROBUSTNESS
text: While an elegant concept, the DRR suffers from a lack of tractability. That is, in order to effectively minimise it, we first need to be able to compute or estimate it. When the loss function is convex with respect to the input space this is straight-forward, however in general approximations are necessary. Our first contribution is an upper bound for it.
For a function f : X → R̄ there is another function co f : X → R̄, called the convex envelope of f . It is the greatest closed convex function that minorises f . The quantity ρ(f) := supx∈X(f(x)− co f(x)) was first suggested by Aubin & Ekeland (1976) to quantify the lack of convexity of a function, and has since shown to be of considerable interest for a variety of nonconvex applications (Askari et al., 2019; Kerdreux et al., 2019; Lemaréchal & Renaud, 2001; Udell & Boyd, 2016). When X = Rn there are well-known ways to compute both co f and ρ(f), and a brief discussion on these appears in the appendix (Remark 2 on p. 18). ρ(f) = 0 when f is closed convex. Theorem 1. Assume X is a separable Fréchet space and fix µ ∈ P(X). Assume c : X → R̄≥0 is sublinear and continuous, and f ∈ L1(X,µ) is upper semicontinuous. Then for all r ≥ 0,
DRR := sup ν∈Bc(µ,r) risk`(f, ν) ≤ r lipc(`f ) + risk`(f, µ). (3)
The tightness of the bound can be quantified as follows. Let ∆(µ) := r lipc(`f ) + risk`(f, µ) − supν∈Bc(µ,r) risk`(f, ν). If lipc(f) <∞ then
∆(µ) ≤ r ( lipc(`f )− [ lipc(co `f )− 1
r
∫ (`f − co `f ) dµ ] + ) , (4)
where [ · ]+ := max{ · , 0} and 1/0 :=∞, so that when `f is closed convex there is equality in (3).
Clearly (4) is tight for convex `f . Furthermore, Proposition 1 shows that (4) is also tight for a large family of nonconvex functions and distributions — particularly the upper-semicontinuous loss functions on a compact set X0 ⊆ X , with the collection of probability distributions supported on X0. Proposition 1. Assume X is a separable Fréchet space with X0 ⊆ X . Assume c : X → R̄≥0 is sublinear and continuous, and `f ∈ ⋂ µ∈P(X0) L1(X,µ) is upper semicontinuous, has lipc(`f ) <∞, and attains its maximum on X0. Then for all r ≥ 0 with 1/0 :=∞,
sup µ∈P(X0)
∆(µ) = r ( lipc(`f )− [ lipc(co `f )− 1
r ρ(`f ) ] + ) .
Theorem 1 subsumes many existing results (viz. Gao & Kleywegt, 2016, Cor. 2 (iv), Cisse et al., 2017, §3.2, Sinha et al., 2018, Shafieezadeh-Abadeh et al., 2017, Thm. 3.20) with a great deal more generality, applying to a very broad family of models, loss functions, and outcome spaces. It is the first time to our knowledge that the slackness in (3) has been characterised tightly.
The extension of Theorem 1 for robust classification in the absence of label noise is straight-forward. Corollary 1. Assume X is a separable Fréchet space and Y is a topological space. Fix µ ∈ P(X × Y ). Assume c : (X × Y )× (X × Y )→ R̄ satisfies c(x, y, x′, y′) = cX(x− x′) whenever
y = y′ and c(x, y, x′, y′) =∞ whenever y 6= y′, where cX : X → R̄ is symmetric, sublinear, and continuous, and f ∈ L1(X × Y, µ) is upper semicontinuous. Then for all r ≥ 0 there is (3). To see the tightness of the bound, if lipc(`f ) < ∞ there is (4), where the closed convex hull is interpreted as co(`f )(x, y) := co(`f ( · , y))(x). If additionally `f ( · , y) is closed convex for all y ∈ Y , there is equality in (3).

heading: 3.1 DISTRIBUTIONAL ROBUSTNESS AS ADVERSARIAL ROBUSTNESS
text: We next show how Theorem 1 can be useful for adversarial learning. The following objective function has been proposed to build a robust classifier. Let X and Y be topological spaces, fix µ ∈ P(X × Y ) and let d be a metric on X . The following objective has been proposed (viz Goodfellow et al., 2015; Madry et al., 2018; Shaham et al., 2018) as a means of learning models that are robust to adversarial perturbations
adversarial risk := ∫
sup x̃∈Bd(x,r)
`f (x̃, y)µ(dx× dy) = ∫
sup ω̃∈Bd̃(ω,r) `f (ω̃)µ(dω), (5)
where in the equality we extend d to a metric on Ω := X × Y with d̃((x, y), (x′, y′)) := d(x, x′) +∞ Jy 6= y′K .
We refer to (5) as the adversarial risk. Theorem 2. Assume (X, c) is a separable Banach space. Fix µ ∈ P(X) and let Rµ(r) := {g ∈ L0(X,R≥0) | ∫ g dµ ≤ r}. Then for f ∈ L0(Ω, R̄), r > 0 there is
variable-radius risk := sup g∈Rµ(r)
∫ µ(dω) sup
ω′∈Bc(ω,g(ω)) `f (ω ′) ≤ sup ν∈Bc(µ,r) risk`(f, ν) = DRR. (6)
The equality holds in (6) if µ is non-atomically concentrated on a compact subset of X , on which f is continuous with the subspace topology.
We refer to the left-hand side (LHS) of (6) as the variable-radius risk. The variable-radius risk has appeared in various forms in similar results, usually formulated using empirical distributions, that is, an average of Dirac masses, (viz. Gao & Kleywegt, 2016; Shafieezadeh-Abadeh et al., 2017). Of course any finite set is compact, and so any empirical distribution satisfies the concentration assumption. Likewise the subspace topology on a finite set is the discrete topology, which makes the continuity assumption trivial.
Both the adversarial risk and the variable-radius risk imply an uncertainty set over a collection of adversaries that may perturb the data. Figure 5 in the appendix (on p. 20) shows the practical difference between the kinds of adversaries in these uncertainty sets. Immediately there is a corollary similar to Corollary 1 for Theorem 2.
It is easy to see that the variable-radius risk upper bounds the adversarial risk (5) by observing that the constant function gr ≡ r is included in the supremum over Rµ(r) in (6). As a result,
adversarial risk ≤ variable-radius risk (a) ≤ DRR (b)
≤ Lipschitz regularised risk (RHS of (3)), (7) where (a) is by Theorem 2 and (b) is by Theorem 1.
In general, it is difficult to characterise the tightness of the upper bounds in Theorem 1 and 2. So we resorted to an empirical demonstration that the sum of all the three gaps in (7) is relatively low. We randomly generated 100 Gaussian kernel classifiers f = ∑100 i=1 γik(x
i, ·), with xi sampled from the MNIST dataset and γi sampled uniformly from [−2, 2]. The bandwidth was set to the median of pairwise distances. In Figure 1, the x-axis is the adversarial risk in (5) where the perturbation δ is bounded in `p ball and computed by PGD. The y-axis is the Lipschitz regularised empirical risk. The scattered dots lie closely to the diagonal, demonstrating that the above bounds are tight in practice.

heading: 4 PROVABLE LIPSCHITZ REGULARISATION FOR KERNEL METHODS
text: Theorems 1 and 2 open up a new path to optimising the adversarial risk (5) by Lipschitz regularisation (RHS of (4)), where the upper bounding relationship is established through DRR. In general, however, it is still hard to compute the Lipschitz constant for a nonlinear model. However, we will show that
0 5 10 0
5
10
(a) ‖δ‖2 ≤ 3
0 10 20 0
10
20
(b) ‖δ‖∞ ≤ 0.3
Figure 1: Empirical evaluation of the sum of the gaps from Theorems 1 and 2. The Lipschitz constants supx∈X ‖∇f(x)‖q (left: p = 2, right: p=∞, 1/p+1/q=1) were estimated by BFGS.
0 20 40 60 80 0
20
40
60
80
(a) 5-layer inverse kernel
0 10 20 30 0
10
20
30
(b) Gaussian kernel (σ=3)
Figure 2: Comparison of λmax(G>G) and the RHS of (8), as upper bounds for the Lipschitz constant. Smaller values are tighter. 100 functions sampled in the same way as in Figure 1.
for some types of kernels, this can be done efficiently on functions in its RKHS. Thanks to the known connections between kernel method and deep learning, this technique will also potentially benefit the latter. For example, `1-regularised neural networks are compactly contained in the RKHS of multi-layer inverse kernels k(x, y) = (2− x>y)−1 with ‖x‖2 ≤ 1 and ‖y‖2 ≤ 1 (Zhang et al., 2016, Lemma 1 and Theorem 1) and (Shalev-Shwartz et al., 2011; Zhang et al., 2017), and even possibly Gaussian kernels k(x, y) = exp(−‖x− y‖2 /(2σ2)) (Shalev-Shwartz et al., 2011, §5). Let us consider a Mercer’s kernel k on a convex domain X ⊆ Rd, with the corresponding RKHS denoted as H. The standard kernel method seeks a discriminant function f from H with the conventional form of finite kernel expansion f(x) = 1l ∑l a=1 γak(x
a, ·), such that the regularised empirical risk can be minimised with the standard (hinge) loss and RKHS norm. We start with real-valued f for univariate output such as binary classification, and later extend it to multiclass.
Our goal here is to additionally enforce, while retaining a convex optimisation in γ := {γa}, that the Lipschitz constant of f falls below a prescribed threshold L > 0, which is equivalent to supx∈X ‖∇f(x)‖2 ≤ L thanks to the convexity of X . A quick but primitive solution is to piggyback on the standard RKHS norm constraint ‖f‖H ≤ C, in view that it already induces an upper bound on ‖∇f(x)‖2 as shown in Example 3.23 of Shafieezadeh-Abadeh et al. (2017),
sup x∈X ‖∇f(x)‖2 ≤ ‖f‖H sup z>0 z−1g(z), where g(z) ≥ sup x,x′∈X:‖x−x′‖2=z ‖k(x, ·)− k(x′, ·)‖H . (8) For Gaussian kernels, g(z) = max{σ−1, 1}z. For exponential and inverse kernels, g(z) = z (Bietti & Mairal, 2019). Bietti et al. (2019) justified that the RKHS norm of a neural network may serve as a surrogate for Lipschitz regularisation. But the quality of such an approximation, i.e., the gap in (8), can be loose as we will see later in Figure 2. Besides, C and L are independent parameters.
How can we tighten the approximation? A natural idea is to directly bound the gradient norm at n random locations {ws}ns=1 sampled i.i.d. from X . These are obviously convex constraints on γ. But how many samples are needed in order to ensure ‖∇f(x)‖2 ≤ L+ ε for all x ∈ X? Unfortunately, as shown in Appendix A.1, n may have to grow exponentially by 1/εd for a d-dimensional space. Therefore we seek a more efficient approach by first slightly relaxing ‖∇f(x)‖2. Let gj(x) := ∂jf(x) be the partial derivative with respect to the j-th coordinate of x, and ∂i,jk(x, y) be the partial derivative to xi and yj . i or j being 0 means no derivative. Assuming supx∈X k(x, x) = 1 and gj ∈ H (true for various kernels considered by Assumptions 1 and 2 below), we get a new bound
sup x∈X ‖∇f(x)‖22 = sup x∈X ∑d j=1 〈gj , k(x, ·)〉2H ≤ sup
ϕ:‖ϕ‖H=1
∑d j=1 〈gj , ϕ〉2H = λmax(G >G), (9)
where λmax evaluates the maximum eigenvalue, and G := (g1, . . . , gd). The “matrix” is only a notation because each column is a function in H, and obviously the (i, j)-th entry of G>G is 〈gi, gj〉H. Interestingly, λmax(G
>G) delivers significantly lower (i.e., tighter) value in approximating the Lipschitz constant supx∈X ‖∇f(x)‖2, compared with ‖f‖Hmaxz>0 g(z) z from (8). Figure 2 compared these two approximants, where λmax(G>G) was computed from (11) derived below, and the landmarks {ws} consisted of all training examples; drawing more samples led to little difference. Such a positive result motivated us to develop refined algorithms to address the only remaining obstacle to leveraging λmax(G>G): no analytic form for computation. Interestingly, it is readily approximable in both theory and practice. Indeed, the role of gj can be approximated by g̃j , where g̃j ∈ Rn is the Nyström approximation (Drineas & Mahoney, 2005; Williams & Seeger, 2000):
a g̃j := K−1/2(gj(w1), . . . , gj(wn))>=(Z>Z)−1/2Z>gj (noting gj(w1) = 〈 gj , k(w 1, ·) 〉 H) (10)
where K := [k(wi, wi ′ )]i,i′ , Z := (k(w
1, ·), k(w2, ·), . . . , k(wn, ·)), G̃ := (g̃1, . . . , g̃d). So to ensure λmax(G>G) ≤ L2 + ε, intuitively we can resort to enforcing λmax(G̃>G̃) ≤ L2, which also retains the convexity in the constraint in γ. However, to guarantee ε error, the number of samples (n) required is generally exponential (Barron, 1994). Fortunately, we will next show that n can be reduced to polynomial for quite a general class of kernels that possess some decomposed structure.

heading: 4.1 A COORDINATE-WISE NYSTRÖM APPROXIMATION FOR PRODUCT KERNELS
text: A number of kernels factor multiplicatively over the coordinates, such as periodic kernels (MacKay, 1998), Gaussian kernels, and Laplacian kernels. We will consider k(x, y) = ∏d j=1 k0(xj , yj) where X = Xd0 and k0 is a base kernel on X0. Let the RKHS of k0 be H0, and let µ0 be a finite Borel measure with supp[µ0] = X0. Periodic kernels have k0(xj , yj) = exp ( −sin ( π v (xj − yj) )2 /(2σ2) ) .
The key benefit of this decomposition is that the derivative ∂0,1k(x, y) can be written as ∂0,1k0(x1, y1) ∏d j=2 k0(xj , yj). Since k0(xj , yj) can be easily dealt with, approximation will be
needed only for ∂0,1k0(x1, y1). Applying this idea to g1 = 1l ∑l a=1 γa∂ 0,1k(xa, ·), we can derive
‖g1‖2H = l −2 ∑l
a,b=1 γaγb
〈 ∂0,1k0(x a 1 , ·), ∂0,1k0(xb1, ·) 〉 H0 ∏d j=2 k0(x a j , x b j), (11)
〈g1, g2〉H = l −2 ∑l
a,b=1 γaγb∂
0,1k0(x a 1 , x b 1)∂ 0,1k0(x b 2, x a 2) ∏d
j=3 k0(x
a j , x b j).
So the off-diagonal entries ofG>G can be computed exactly. To approximate the diagonal, we sample {w11, . . . , wn1 } from µ0, set Z1 = (k0(w11, ·), . . . , k0(wn1 , ·)), and apply Nyström approximation:〈
∂0,1k0(x a 1 , ·), ∂0,1k0(xb1, ·) 〉 H0 ≈ ∂0,1k0(xa1 , ·)>Z1 · (Z>1 Z1)−1 · Z>1 ∂0,1k0(xb1, ·) (12)
where Z>1 ∂ 0,1k0(x a 1 , ·) = (∂0,1k0(xa1 , w11), . . . , ∂0,1k0(xa1 , wn1 ))>, (13)
and analogously for Z>1 ∂ 0,1k0(x b 1, ·). We will denote this approximation of G>G as P̃G. Clearly, λmax(P̃G) ≤ L2 is a convex constraint on γ, based on i.i.d. samples {wsj : s ∈ [n], j ∈ [d]} from µ0. It is now important to analyse how many samples wsj are needed, such that
λmax(P̃G) ≤ L2 =⇒ λmax(G>G) ≤ L2 + ε with high probability.

heading: 4.2 GENERAL SAMPLE COMPLEXITY AND ASSUMPTIONS ON THE PRODUCT KERNEL
text: Fortunately, product kernels only require approximation bounds for each coordinate, making the sample complexity immune to the exponential growth in the dimensionality d. Specifically, we first consider base kernels k0 with a scalar input, i.e., X0 ⊆ R. Recall from Steinwart & Christmann (2008, Chapter 4) that the integral operator for k0 and µ0 is defined by
Tk0 = I ◦ S : L2(X0, µ0)→ L2(X0, µ0) where S : L2(X0, µ0)→ C(X0), (Sf)(x) = ∫ k0(x, y)f(y)dµ0(y), f ∈ L2(X0, µ0),
and I: C(X0) ↪→ L2(X0;µ0) is the inclusion operator. By the spectral theorem, if Tk0 is compact, then there is an at most countable orthonormal set {ẽj}j∈J of L2(X0, µ0) and {λj}j∈J with λ1 ≥ λ2 ≥ . . . > 0 such that Tk0f = ∑ j∈J λj 〈f, ẽj〉L2(X0;µ0) ẽj for all f ∈ L2(X0, µ0). It is easy to
see that ϕj := √ λjej is an orthonormal basis ofH0 (Steinwart & Christmann, 2008).
Our proof is built upon the following two assumptions on the base kernel. The first one asserts that fixing x, the energy of k0(x, ·) and ∂0,1k0(x, ·) “concentrates” on the leading eigenfunctions. Assumption 1. Suppose k0(x, x) = 1 and ∂0,1k0(x, ·) ∈ H0 for all x ∈ X0. For all ε > 0, there exists Nε ∈ N such that the tail energy of ∂0,1k0(x, ·) beyond the Nε-th eigenpair is less than ε, uniformly for all x ∈ X0. That is, denoting Φm := (ϕ1, . . . , ϕm),
Nε := inf m {∥∥∂0,1k0(x, ·)− ΦmΦ>m∂0,1k0(x, ·)∥∥H0 < ε for all x ∈ X0 and∥∥k0(x, ·)− ΦmΦ>mk0(x, ·)∥∥H0 < ε for all x ∈ X0} <∞.
The second assumption asserts the smoothness and range of eigenfunctions in a uniform sense.
Assumption 2. Under Assumption 1, {ej(x) : j ∈ Nε} is uniformed bounded over x ∈ X0, and the RKHS inner product of ∂0,1k0(x, ·) with {ej : j ∈ Nε} is also uniformly bounded over x ∈ X0:
Mε := sup x∈X0 max j∈[Nε] ∣∣∣〈∂0,1k0(x, ·), ej〉H0 ∣∣∣ <∞, and Qε := supx∈X0 maxj∈[Nε] |ej(x)| <∞. Theorem 3. Suppose k0, X0, and µ0 satisfy Assumptions 1 and 2. Let {wsj : s ∈ [n], j ∈ [d]} be sampled i.i.d. from µ0. Then for any f whose coordinate-wise Nyström approximation (11) and (12) satisfy λmax(P̃G) ≤ L2, the Lipschitz condition λmax(G>G) ≤ L2 + ε is met with probability 1− δ, as long as n ≥ Θ̃ ( 1 ε2N 2 εM 2 εQ 2 ε log dNε δ ) , almost independent of d. Here Θ̃ hides all poly-log terms.
Satisfaction of Assumptions. In Appendix A.4 and A.5, we will show that for periodic kernel and Gaussian kernel, Assumptions 1 and 2 hold true with Õ(1) values ofNε,Mε, andQε. It remains open whether non-product kernels such as inverse kernel also enjoy this polynomial sample complexity. Appendix A.6 suggests that the complexity is quasi-polynomial for inverse kernels.

heading: 5 EXPERIMENTAL RESULTS
text: We studied the empirical robustness and accuracy of the proposed Lipschitz regularisation technique for adversarial training of kernel methods, under both Gaussian kernel and inverse kernel. Comparison will be made with state-of-the-art defense algorithms under effective attacks.
Datasets. We tested on three datasets: MNIST, Fashion-MNIST, and CIFAR10. The number of training/validation/test examples for the three datasets are 54k/6k/10k, 54k/6k/10k, 45k/5k/10k, respectively. Each image in MNIST and Fashion-MNIST is represented as a 784-dimensional feature vector, with each feature/pixel normalised to [0, 1]. For CIFAR10, we trained it on a residual network to obtain a 512-dimensional feature embedding, which were subsequently normalised to [0, 1]. They were used as the input for training all the competing algorithms and were subject to attack.
Attacks. To evaluate the robustness of the trained model, we attacked them on test examples using the random initialized Projected Gradient Descent method with 100 steps (PGD, Madry et al., 2018) under two losses: cross-entropy and C&W loss (Carlini & Wagner, 2017). The perturbation δ was constrained in an `2 or `∞ ball. To evaluate robustness, we scaled the perturbation bound δ from 0.1 to 0.6 for `∞ norm, and from 1 to 6 for `2 norm (when δ = 6, the average magnitude per coordinate is 0.214).
Algorithms. We compared four training algorithms. The Parseval network orthonormalises the weight matrices to enforce the Lipschitz constant (Cisse et al., 2017). We used three hidden layers of 1024 units and ReLU activation (Par-ReLU). Also considered is the Parseval network with MaxMin activations (Par-MaxMin), which enjoys much improved robustness (Anil et al., 2019). Both algorithms can be customised for `2 or `∞ attacks, and were trained under the corresponding norms. Using multi-class hinge loss, they constitute strong baselines for adversarial learning.
Both Gaussian and inverse kernel machines applied Lipschitz regularisation by randomly and greedily selecting {ws}, and they will be referred to as Gauss-Lip and Inverse-Lip, respectively. In practice, Gauss-Lip with the coordinate-wise Nyström approximation (λmax(P̃G) from Eq (12)) can approximate λmax(G>G) with a much smaller number of sample than if using the holistic approximation as in (10). Furthermore, we found an even more efficient approach. Inside the iterative training algorithm, we used L-BFGS to find the input that yields the steepest gradient under the current solution, and then added it to the set {ws} (which was initialized with 15 random points). Although L-BFGS is only a local solver, this greedy approach empirically reduces the number of samples by an order of magnitude. See the empirical convergence results in Appendix A.9. Its theoretical analysis is left for future investigation. We also applied this greedy approach to Inverse-Lip.
Extending binary kernel machines to multiclass. The standard kernel methods learn a discriminant function f c := ∑ a γ c ak(x
a, ·) for each class c ∈ [10], based on which a large supply of multiclass classification losses can be applied, e.g., CS (Crammer & Singer, 2001) which was used in our experiment. Since the Lipschitz constant of the mapping from {f c} to a real-valued loss is typically at most 1, it suffices to bound the Lipschitz constant of x 7→ (f1(x), . . . , f10(x))> by
aa max x λmax(G(x)G(x) >) ≤ max‖ϕ‖H=1 λmax
(∑10
c=1 G>c ϕϕ >Gc
)
≤ L2, (14)
where G(x) := [∇f1(x), · · · ,∇f10(x)] = [G>1 k(x, ·), . . . , G>10k(x, ·)] with Gc := (gc1, . . . , gcd). The last term in (14) can be approximated using the same technique as in the binary case. Furthermore, the principle can be extended to `∞ attacks, whose details are relegated to Appendix A.10.
Parameter selection. We used the same parameters as in Anil et al. (2019) for training Par-ReLU and Par-MaxMin. To defend against `2 attacks, we set L = 100 for all algorithms. GaussLip achieved high accuracy and robustness on the validation set with bandwidth σ = 1.5 for FashionMNIST and CIFAR-10, and σ = 2 for MNIST. To defend against `∞ attacks, we set L = 1000 for all the four methods as in Anil et al. (2019). The best σ for Gauss-Lip is 1 for all datasets. Inverse-Lip used 5 stacked layers.
Results. Figures 3 and 4 show how the test accuracy decays as an increasing amount of perturbation (δ) in `2 and `∞ norm is added to the test images, respectively. Clearly Gauss-Lip achieves higher accuracy and robustness than Par-ReLU and Par-MaxMin on the three datasets, under both `2 and `∞ bounded PGD attacks with C&W loss. In contrast, Inverse-Lip only performs similarly to Par-ReLU. Interestingly, we noticed that `2 based Par-MaxMin are only slightly better than Par-ReLU under `2 attacks, although the former does perform significantly better under `∞ attacks.
For the sake of space, the results for cross-entropy PGD attacks are deferred to Figures 8 and 9 in Appendix A.11. Here cross-entropy PGD attackers find stronger attacks to Parseval networks but not to our kernel models. Our Gauss-Lip again significantly outperforms Par-MaxMin on all the three datasets and under both `2 and `∞ norms. The improved robustness of Gauss-Lip does not seem to be attributed to the obfuscated gradient (Athalye et al., 2018), because as shown Figures 3, 4, 8, 9, increased distortion bound does increase attack success, and unbounded attacks drive the success rate to very low. In practice, we also observed that random sampling finds much weaker attacks, and taking 10 steps of PGD is much stronger than just one step.
Visualization. The gradient with respect to inputs is plotted in Figure 10 (in the appendix on p. 31) for `2 trained Par-MaxMin and Gauss-Lip. The i-th row and j-th column corresponds to the targeted attack of turning the original class j into a new class i, hence the gradient is on the cross-entropy loss with class i as the ground truth. These two figures also explained why Gauss-Lip is more robust than Par-MaxMin: the attacker can easily reduce the targeted cross-entropy loss by following the gradient as shown in Figure 10a, and hence successfully attack Par-MaxMin. In contrast, the gradient shown in Figure 10b does not provide much information on how to flip the class.
Conclusion. In this paper, we derived a new certificate for distributional robust risk minimization by using Lipschitz regularization. Application to adversarial learning based on kernel methods exhibited superior robustness, with provably polynomial sample complexity for product kernels. We will apply this function space to GANs to witness the difference between probability distributions, leading to a more stable training scheme as the inner level optimization becomes convex.

heading: Appendix
text: The pseudo-code of training binary SVMs by enforcing Lipschitz constant is given in Algorithm 1.
Algorithm 1: Training binary SVMs by enforcing Lipschitz constant L 1 Initialise the constraint set S by some random samples from X . 2 for i = 1, 2, . . . do 3 Train SVM using one of the following constraints:
1© Brute-force: ‖∇f(w)‖22 ≤ L2, ∀ w ∈ S
2© Nyström holistic: λmax(G̃>G̃) ≤ L2 using S as the set {w1, . . . , wn} in Eq (10)
3© Nyström coordinate wise: λmax(P̃G) ≤ L2 using S as the set {w1, . . . , wn} in Eq (12)
4 Let the trained SVM be f (i). 5 Find a new w to add to S by one of the following methods:
a© Random: randomly sample w from X . b© Greedy: find argmaxx∈X ∥∥∇f (i)(x)∥∥ (local optimisation) by L-BFGS with 10 random initialisations. Add the distinct results upon convergence to S.
6 Return if L(i) := maxx∈X ∥∥∇f (i)(x)∥∥ falls below L.
Finding the exact argmaxx∈X ∥∥∇f (i)(x)∥∥ is intractable, so we used a local maximum found by L-BFGS with 10 random initialisations as the Lipschitz constant of the current solution f (i) (L(i) in step 6). The solution found by L-BFGS is also used as the new greedy point added in step 5b.
Furthermore, the kernel expansion f(x) = 1l ∑l a=1 γak(x
a, ·) can lead to high cost in optimisation (our experiment used l = 54000), and therefore we used another Nyström approximation for the kernels. We randomly sampled 1000 landmark points, and based on them we computed the Nyström approximation for each k(xa, ·), denoted as ϕ̃(xa) ∈ R1000. Then f(x) can be written as 1 l ∑l a=1 γaϕ̃(x a)>ϕ̃(x). Defining w = 1l ∑l a=1 γaϕ̃(x
a), we can equivalently optimise over w, and the RKHS norm bound on f can be equivalently imposed as the `2-norm bound on w.
To summarise, Nyström approximation is used in two different places: one for approximating the kernel function, and one for computing ‖gj‖H either holistically or coordinate wise. For the former, we randomly sampled 1000 landmark points; for the latter, we used greedy selection as option b in step 5 of Algorithm 1.
Detailed algorithm for multiclass classification. It is easy to extend Algorithm 1 to multiclass. For example, with MNIST dataset, we solve the following optimisation problem to defend `2 attacks:
min γ1,...,γ10 n∑ i=1 `(F (x),y), where F = [ n∑ i=1 γ1i k(xi, ·); . . . ; n∑ i=1 α10i k(xi, ·) ]
s.t. sup ‖ϕ‖H≤1 λmax ( 10∑ c=1 G>c ϕϕ >Gc ) ≈ sup ‖v‖2≤1 λmax ( 10∑ c=1 G̃>c vv >G̃c ) ≤ L2,
where `(F (x),y) is the Crammer & Singer loss, and the constraint is derived from (14) by using its Nyström approximation G̃c = [g̃c1, . . . , g̃ c d], which depends on {γ1, . . . ,γ10} linearly. Note that the constraint itself is a supremum problem:
sup ‖v‖2≤1 λmax ( 10∑ c=1 G̃>c vv >G̃c ) = sup ‖v‖2≤1,‖u‖2≤1 u> ( 10∑ c=1 G̃>c vv >G̃c ) u.
Since there is only one constraint, interior point algorithm is efficient. It requires the gradient of the constraint, which can be computed by Danskin’s theorem. In particular, we alternates between
updating v and u, until they converge to the optimal v∗ and u∗. Finally, the derivative of the constraint with respect to {γc} can be calculated from ∑10 c=1(u > ∗ G̃ > c v∗)
2, as a function of {γc}. To defend `∞ attacks, we need to enforce the `∞ norm of the Jacobian matrix:
sup x∈X ∥∥∥[g1(x), . . . , g10(x)]>∥∥∥ ∞ = sup x∈X max 1≤c≤10 ‖gc(x)‖1
= max 1≤c≤10 sup x∈X ‖gc(x)‖1
≤ max 1≤c≤10 sup ‖ϕ‖2≤1,‖u‖∞≤1 u>G̃>c ϕ,
where the last inequality is due to
sup x∈X ‖g(x)‖1 = sup x∈X sup ‖u‖∞≤1 u>g(x) ≤ sup ‖v‖2≤1,‖u‖∞≤1 u>G̃>v.
Therefore, the overall optimisation problem to defense `∞ attacks is
min γ1,...,γ10 n∑ i=1 `(F (x),y), where F = [ n∑ i=1 γ1i k(xi, ·); . . . ; n∑ i=1 γ10i k(xi, ·) ] s.t. sup ‖v‖2≤1,‖u‖∞≤1 u>G̃>c v ≤ L, ∀c ∈ {1, . . . , 10} (15)
For each c, we alternatively update v and u in (15), converging to the optimal v∗ and u∗. Finally, the derivative of sup‖v‖2≤1,‖u‖∞≤1 u >G̃>c v with respect to γ c can be calculated from u>∗ G̃ > c v∗, as a function of γc.

heading: A PROOFS OF RESULTS
text: The following appendix contains the complete set of proofs and auxiliary results.

heading: PROOFS FOR §3: CERTIFICATE FOR DISTRIBUTIONAL ROBUSTNESS
text: Duality results like Lemma 1 have been the basis of a number of recent theoretical efforts in the theory of adversarial learning (Blanchet et al., 2016; Gao & Kleywegt, 2016; Shafieezadeh-Abadeh et al., 2017; Sinha et al., 2018), the results of Blanchet & Murthy (2019) being the most general to date. Lemma 1 (Blanchet & Murthy (2019, Thm. 1)). Assume Ω is a Polish space and fix µ ∈ P(Ω). Let c : Ω ×Ω → R̄≥0 be lower semicontinuous with c(ω, ω) = 0 for all ω ∈ Ω, and f ∈ L1(Ω,µ) is upper semicontinuous. Then for all r ≥ 0 there is
sup ν∈Bc(µ,r)
∫ f dν = inf
λ≥0
( λr + ∫ fλc dµ ) . (16)
The necessity for such duality results like Lemma 1 is because while the supremum on the left hand side of (16) is over a (usually) infinite dimensional space, the right hand side only involves only a finite dimensional optimisation. The generalised conjugate in (16) also hides an optimisation, but when the outcome space Ω is finite dimensional, this too is a finite dimensional problem.
The following is sometimes stated a consequence of or in the proof of the McShane–Whitney extension theorem, but it is immediate to observe. Lemma 2 (McShane–Whitney). Let X be a set. Assume c : X ×X → R̄≥0 satisfies c(x, x) = 0 for all x ∈ X , and f : X → R. Then
∀x, y ∈ X : f(x)− f(y) ≤ λc(x, y) =⇒ ∀y ∈ X : f(y) = sup x∈X
( f(x)− λc(x, y) ) .
Lemma 3. Assume X is a locally convex Hausdorff topological vector space. Let c : X → R̄ be lower semicontinuous, sublinear, and continuous at 0, let f : X → R̄ be closed convex. Then for λ > 0 there is
∀y ∈ X : sup x∈X
( f(x)− λc(x− y) ) = { f(y) ∂f(X) ⊆ ∂λc(0) ∞ ∂f(X) 6⊆ ∂λc(0).
Proof. Because f is closed convex, it is equal to its biconjugate (Zălinescu, 2002, Thm. 2.3.3), because c is sublinear and lower semicontinuous λc(x) = supx∗∈∂λc(0) 〈x, x∗〉 for all x ∈ X (Zălinescu, 2002, Thm. 2.4.14 (iv)). It follows that
sup x∈X
( f(x)− λc(x− y) ) = sup x∈X sup x∗∈∂f(X) inf g∗∈∂λc(0) ( 〈x, x∗〉 − f∗(x∗)− 〈g∗, x− y〉 ) = sup x∈X sup x∗∈∂f(X) inf g∗∈∂λc(0) ( 〈x, x∗ − g∗〉+ 〈y, g∗〉 − f∗(x∗) ) .
Because c is continuous at 0, ∂c(0) is weak∗-compact and convex (Zălinescu, 2002, Thm. 2.4.9), and so we can apply a minimax theorem (Zălinescu, 2002, Thm. 2.10.2) to produce
sup x∈X sup x∗∈∂f(X) inf g∗∈∂λc(0)
( 〈x, x∗ − g∗〉+ 〈y, g∗〉 − f∗(x∗) ) = sup x∗∈∂f(X) inf g∗∈∂λc(0) sup x∈X ( 〈x, x∗ − g∗〉+ 〈y, g∗〉 − f∗(x∗)
) = sup x∗∈∂f(X) inf g∗∈∂λc(0) { 〈y, x∗〉 − f∗(x∗) g∗ = x∗ ∞ g∗ 6= x∗
= sup x∗∈∂f(X) { 〈y, x∗〉 − f∗(x∗) x∗ ∈ ∂λc(0) ∞ x∗ /∈ ∂λc(0)
= { f(y) ∂f(X) ⊆ ∂λc(0) ∞ ∂f(X) 6⊆ ∂λc(0),
as claimed.
Remark 1. The minimisation of g−h, where g and h are convex functions, is called difference convex (DC) programming (Hiriart-Urruty, 1989). The condition in Lemma 3 bears a striking resemblance to the common necessary condition (e.g. Hiriart-Urruty, 1989; Penot, 1998) for such problems
x ∈ arginf x′∈X
f(x′) =⇒ ∂h(x) ⊆ ∂g(x).
Likewise there are similar sufficient conditions. The proof of Lemma 4 is also quite similar to the proofs of the Toland (1979) duality formula (viz. Hiriart-Urruty, 1986)
inf x∈X (g(x)− h(x)) = inf x∗∈X∗
(h∗(x∗)− g∗(x∗)),
which suggests that the principles of Lemma 4 may be more general. A generalisation to a general convex function c satisfying c(0) = 0, would remove the positive homogeneity requirement of Lemma 4, and allow any translation invariant metric in place of c. The assumptions we have made are compatible with metrics which arise from norms, that is, the translation invariant and positively homogeneous metrics. Lemma 4. Assume X is a topological vector space. Let c : X → R̄≥0, and f : X → R. Then for λ > 0 there is
∀x, y ∈ X : f(x)− f(y) ≤ λc(x− y) ⇐⇒ ∂f(X) ⊆ ∂λc(0).
Proof. Suppose λ > 0 is such that f(x) − f(y) ≤ λc(x − y) for all x, y ∈ X . Let x∗ ∈ ∂f(X). Then there is x ∈ X with
∀y ∈ X : 〈y − x, x∗〉 ≤ f(y)− f(x) ≤ λc(y − x) =⇒ ∀y ∈ X : 〈y, x∗〉 ≤ f(y + x)− f(x) ≤ λc(y),
this shows x∗ ∈ ∂λc(0). Next assume λ > 0 satisfies ∂f(X) ⊆ ∂λc(0). Then ∀x ∈ X, ∃x∗ ∈ ∂f(x),∀y ∈ X : f(x)− f(y) ≤ 〈x− y, x∗〉 ≤ λc(x− y),
where the second inequality is because x∗ ∈ ∂λc(0) for all x∗ ∈ ∂f(x).
Theorem 1. Assume X is a separable Fréchet space and fix µ ∈ P(X). Assume c : X → R̄≥0 is sublinear and continuous, and f ∈ L1(X,µ) is upper semicontinuous. Then for all r ≥ 0,
DRR := sup ν∈Bc(µ,r) risk`(f, ν) ≤ r lipc(`f ) + risk`(f, µ).
The tightness of the bound can be quantified as follows. Let ∆(µ) := r lipc(`f ) + risk`(f, µ) − supν∈Bc(µ,r) risk`(f, ν). If lipc(f) <∞ then
∆(µ) ≤ r ( lipc(`f )− [ lipc(co `f )− 1
r
∫ (`f − co `f ) dµ ] + ) ,
where [ · ]+ := max{ · , 0} and 1/0 :=∞, so that when `f is closed convex there is equality in (3).
Proof. Because c is assumed sublinear, it is positively homogeneous and there is c(x, x) = c(x−x) = c(0) = 0 for all x ∈ X . Therefore we can apply Lemma 1 and Lemma 2 to obtain
sup ν∈Bc(µ,r)
∫ `f dν = inf
λ≥0
[ rλ+ ∫ `λcf dµ ] ≤ inf λ≥lipc(`f ) [ rλ+ ∫ `λcf dµ
] = r lipc(`f ) + ∫ `f dµ.
Observing that co `f ≤ `f , applying Lemma 3 and Lemma 4 we find for all x ∈ X
sup λ∈[0,∞)
( `f (x)− `λcf (x)− rλ ) = sup λ∈[0,∞) ( `f (x)− sup y∈X ( `f (y)− λc(x− y) ) − rλ ) = sup λ∈[0,∞) inf y∈X ( `f (x)− `f (y) + λc(x− y)− rλ
) ≤ sup λ∈[0,∞) inf y∈X ( `f (x)− co `f (y) + λc(x− y)− λr
) = sup λ∈[0,∞) ( `f (x)− co `f (x)−∞ Jlipc(co `f ) > λK− λr
) = `f (x)− co `f (x)− r lipc(co `f ). (1)
Similarly, for all x ∈ X there is
sup λ∈[0,∞)
( `f (x)− `λcf (x)− rλ ) ≤ sup λ∈[0,∞) ( `f (x)− `λcf (x) ) + sup λ∈[0,∞) −rλ
= sup λ∈[0,∞)
( `f (x)− `λcf (x) ) = sup λ∈[0,∞) inf y∈X ( `f (x)− `f (y) + λc(x− y)
) ≤ inf y∈X sup λ∈[0,∞) ( `f (x)− `f (y) + λc(x− y)
) = inf y∈X { ∞ c(x− y) > 0 0 c(x− y) = 0
= 0. (2)
Then, using (1) and (2) we find( r lipc(`f ) + ∫ `f dµ ) − inf λ∈[0,∞) ( rλ− ∫ `λcf dµ ) = r lipc(`f ) + sup
λ∈[0,∞)
∫ ( `f − `λcf − λr ) dµ
≤ r lipc(`f ) + ∫
sup λ∈[0,∞)
( `f − `λcf − λr ) dµ
(1),(2) ≤ r lipc(`f ) + min {∫ (`f − co `f ) dµ− r lipc(co `f ), 0 } .
The proof is complete.
Proposition 1. Assume X is a separable Fréchet space with X0 ⊆ X . Assume c : X → R̄≥0 is sublinear and continuous, and `f ∈ ⋂ µ∈P(X0) L1(X,µ) is upper semicontinuous, has lipc(`f ) <∞, and attains its maximum on X0. Then for all r ≥ 0 with 1/0 :=∞,
sup µ∈P(X0)
∆(µ) = r ( lipc(`f )− [ lipc(co `f )− 1
r ρ(`f ) ] + ) .
Proof. Let x0 ∈ X0 be be the point at which `f (x0) = supx∈X0 `f (x). Then
∆(δx0) = r lipc(`f ) + ∫ `f dδx0 − sup
ν∈Bc(δx0 ,r)
∫ `f dδx0
= r lipc(`f ) + ∫ `f dδx0 − ∫ `f dδx0
= r lipc(`f ). (3)
Then there is
r lipc(`f ) (3) ≤ sup
µ∈P(X0) ∆(µ)
(4) ≤ r ( lipc(`f )−max { lipc(co `f )− 1
r ρ(`f ), 0
}) ≤ r lipc(`f ),
which completes the proof.
Remark 2. When f : Rn → R̄ satisfies f 6≡ ∞ and f is minorised by an affine function, there is (cf. Hiriart-Urruty & Lemaréchal, 2010, Prop. 1.5.4)
∀x ∈ Rn : co f(x) = inf  ∑ i∈[n+1] αif(xi) | ∑ i∈[n+1] αi = 1, x = ∑ i∈[n+1] αixi , where the infimum is over all sequences (αi)i∈[n+1] and (xi)i∈[n+1] ⊆ Rn satisfying the conditions above. Consequentially there is the common expression
ρ(f) = sup f( ∑ i∈[n+1] αixi ) − ∑ i∈[n+1] αif(xi) | (αi, xi)i∈[n+1] ⊆ R≥0 ×Rn, ∑ i∈[n+1] αi = 1 . In Lemma 5, by the weak∗ topology on P(Ω) we mean the coarsest topolgoy on P(Ω) that makes the bounded continuous functions on Ω its topological dual space. Likewise ⇀∗ denotes convergence in this topology.
Lemma 5. Assume (Ω, c) is a compact Polish space and µ ∈ P(Ω) is non-atomic. For any ν? ∈ P(Ω) and r > 0 there is a sequence (fi)i∈N ⊆ Aµ(r) := { f ∈ L0(Ω,Ω) | ∫ cd(id, f)#µ ≤ r } with (fi)#µ ⇀∗ ν?.
Proof. Let P (µ, ν) := {f ∈ L0(X,X) | f#µ = ν}. Since µ is non-atomic and c is continuous Pratelli (2007, Thm. B) shows
∀ν ∈ P(Ω) : inf f∈P (µ,ν)
∫ cd(id, f)#µ = costc(µ, ν).
Let r? := costc(µ, ν?), obviously r? ≤ r. Assume r? > 0, otherwise the lemma is trivial. Fix a sequence (εk)k∈N ⊆ (0, r?) with εk → 0. For u ≥ 0 let ν(u) := µ+ u(ν? − µ). Then
costc(µ, ν(0)) = 0 and costc(µ, ν(1)) = r?,
and because costc metrises the weak∗ topology on P(Ω) (Villani, 2008, Thm. 6.9), the mapping u 7→ costc(µ, ν(u)) is continuous. Then by the intermediate value theorem for every k ∈ N there is some uk > 0 with costc(µ, ν(uk)) = r? − εk, forming a sequence (uk)k∈N ⊆ [0, 1]. Then for every k there is a sequence (fjk)j∈N ⊆ P (µ, ν(uk)) so that (fjk)#µ ⇀ ∗ν(k) and
lim j∈N
∫ cd(id, fjk)#µ = inf
f∈P (µ,ν(k))
∫ cd(id, fk)#µ = costc(µ, ν(k)) = r ? − εk.
Therefore for every k ∈ N there exists jk ≥ 0 so that for every j ≥ jk∫ cd(id, fjk)#µ ≤ r?. (3)
Let us pass directly to this subsequence of (fjk)j∈N for every k ∈ N so that (3) holds for all j, k ∈ N. Next by construction we have ν(uk) → ν?. Therefore (fjk)j,k∈N has a subsequence in k so that (fjk)#µ ⇀ ∗ ν?. By ensuring (3) is satisfied, the sequences (fjk)j∈N ⊆ Aµ(r) for every k ∈ N.
Theorem 2. Assume (X, c) is a separable Banach space. Fix µ ∈ P(X) and let Rµ(r) := {g ∈ L0(X,R≥0) | ∫ g dµ ≤ r}. Then for f ∈ L0(Ω, R̄), r > 0 there is
variable-radius risk := sup g∈Rµ(r)
∫ µ(dω) sup
ω′∈Bc(ω,g(ω)) `f (ω ′) ≤ sup ν∈Bc(µ,r) risk`(f, ν) = DRR.
The equality holds in (6) if µ is non-atomically concentrated on a compact subset of X , on which f is continuous with the subspace topology.
Proof. Inequality (6). For g ∈ Rµ(r), let Γg : X ⇒ X denote the set-valued mapping with Γg(x) := Bc(x, g(x)). Let L0(X,Γg) denote the set of Borel a : X → X so that a(x) ∈ Γg(x) for µ-almost all x ∈ X . Let Aµ(r) := ⋃ g∈Rµ(r) L0(X,Γg). Clearly for every a ∈ Aµ(r) there is
r ≥ ∫ c(x, a(x)) dµ = ∫ cd(id, a)#µ,
which shows {a#µ | a ∈ Aµ(r)} ⊆ Bc(µ, r). Then if there is equality in (4), we have
sup g∈Rµ(r)
∫ sup
x′∈Γg(x) f(x) = sup g∈Rµ(r) sup a∈L0(X,Γg)
∫ f da#µ (4)
= sup a∈Aµ(r)
∫ f da#µ
≤ sup ν∈Bc(µ,r)
∫ f da#ν,
which proves the inequality (6).
Equality (4). To complete the proof we will now justify the exchange of integration and supremum in (4). The set L0(X,Γg) is trivially decomposable (Giner, 2009, see the remark at the bottom of p. 323, Def. 2.1). By assumption f is Borel measurable. Since f is measurable, any decomposable subset of L0(X,X) is f -decomposable (Giner, 2009, Prop. 5.3) and f -linked (Giner, 2009, Prop. 3.7 (i)). Giner (2009, Thm. 6.1 (c)) therefore allows us to exchange integration and supremum in (4).
Equality in (6). Under the additional assumptions there exists ν? ∈ P(Ω) with (via Blanchet & Murthy, 2019, Prop. 2) ∫
f dν? = sup ν∈Bc(µ,r)
∫ f dν.
The compact subset where µ is concentrated and non-atomic is a Polish space with the Banach metric. Therefore Using Lemma 5 there is a sequence (fi)i∈N ⊆ Aµ(r) so that
lim i∈N
∫ fi dµ = ∫ f dν? = sup
ν∈Bc(µ,r)
∫ f dν,
proving equality in (6).

heading: PROOFS FOR §4: PROVABLE LIPSCHITZ REGULARISATION FOR KERNEL METHODS
text: Theorem 3. Suppose k0, X0, and ν0 satisfy Assumptions 1 and 2. Let {wsj : s ∈ [n], j ∈ [d]} be sampled i.i.d. from ν0. Then for any f whose coordinate-wise Nyström approximation (11) and (12) satisfy λmax(P̃G) ≤ L2, the Lipschitz condition λmax(G>G) ≤ L2 + ε is met with probability 1− δ, as long as n ≥ Θ̃ ( 1 ε2N 2 εM 2 εQ 2 ε log dNε δ ) , almost independent of d. Here Θ̃ hides all poly-log terms.

heading: PROOFS AND MORE RESULTS FOR §4: KERNEL APPROXIMATION

heading: A.1 RANDOM SAMPLING REQUIRES EXPONENTIAL COST
text: The most natural idea of leveraging the samples is to add the constraints ‖g(ws)‖ ≤ L. For Gaussian kernel, we may sample from N (0, σ2I) while for inverse kernel we may sample uniformly from B. This leads to our training objective:
min f∈H
1
l l∑ i=1 loss(f(xi), yi) + λ 2 ‖f‖2H s.t. ‖g(w s)‖ ≤ L, ∀s ∈ [n].
Unfortunately, this method may require O( 1 εd ) samples to guarantee ∑ j ‖gj‖ 2 H ≤ L
2 + ε w.h.p. This is illustrated in Figure 6, where k is the polynomial kernel with degree 2 whose domain X is the unit ball B, and f(x) = 12 (v
>x)2. We seek to test whether the gradient g(x) = (v>x)v has norm bounded by 1 for all x ∈ B, and we are only allowed to test whether ‖g(ws)‖ ≤ 1 for samples ws that are drawn uniformly at random from B. This is equivalent to testing ‖v‖ ≤ 1, and to achieve it at least one ws must be from the ε ball around v/ ‖v‖ or −v/ ‖v‖, intersected with B. But the probability of hitting such a region decays exponentially with the dimensionality d. The key insight from the above counter-example is that in fact ‖v‖ can be easily computed by∑d s=1(v >w̃s) 2, where {w̃s}ds=1 is the orthonormal basis computed from the Gram–Schmidt process on d random samples {ws}ds=1 (n = d). With probability 1, n samples drawn uniformly from B must span Rd as long as n ≥ d, i.e., rank(W ) = d where W = (w1, . . . , wn). The Gram–Schmidt process can be effectively represented using a pseudo-inverse matrix (allowing n > d) as
‖v‖2 = ∥∥∥(W>W )−1/2W>v∥∥∥
2 ,
where (W>W )−1/2 is the square root of the pseudo-inverse of W>W . This is exactly the intuition underlying the Nyström approximation that we will leveraged.

heading: A.2 SPECTRUM OF KERNELS
text: Let k be a continuous kernel on a compact metric space X , and µ be a finite Borel measure on X with supp[µ] = X . We will re-describe the following spectral properties in a more general way than
in §4. Recall from Chapter 4 of Steinwart & Christmann (2008) that the integral operator for k and µ is defined by
Tk = Ik ◦ Sk : L2(X,µ)→ L2(X,µ) where Sk : L2(X;µ)→ C(X), (Skf)(x) = ∫ k(x, y)f(y)dµ(y), f ∈ L2(X,µ),
Ik : C(X) ↪→ L2(X;µ), inclusion operator.
By the spectral theorem, if Tk is compact, then there is an at most countable orthonormal set (ONS) {ẽj}j∈J of L2(X,µ) and {λj}j∈J with λ1 ≥ λ2 ≥ . . . > 0 such that
Tf = ∑ j∈J λj 〈f, ẽj〉L2(X;µ) ẽj , f ∈ L2(X,µ).
In particular, we have 〈ẽi, ẽj〉L2(X;µ) = δij (i.e., equals 1 if i = j, and 0 otherwise), and T ẽi = λiẽi. Since ẽj is an equivalent class instead of a single function, we assign a set of continuous functions ej = λ −1 j Skẽj ∈ C(X), which clearly satisfies
〈ei, ej〉L2(X;µ) = δij , T ej = λjej .
We will call λj and ej as eigenvalues and eigenfunctions respectively, and {ej}j∈J clearly forms an ONS. By Mercer’s theorem,
k(x, y) = ∑ j∈J λjej(x)ej(y), (5)
and all functions in H can be represented by ∑ j∈J ajej where {aj/ √ λj} ∈ `2(J). The inner
product inH is equivalent to 〈∑ j∈J ajej , ∑ j∈J bjej 〉 H = ∑ j∈J ajbj/λj . Therefore it is easy to see that
ϕj := √ λjej , j ∈ J
is an orthonormal basis of H, with Moreover, for all f ∈ H with f = ∑ j∈J ajej , we have
〈f, ej〉H = aj/λj , 〈f, ϕj〉H = aj/ √ λj , and
f = ∑ j 〈f, ϕj〉H ϕj = ∑ j √ λj 〈f, ej〉H ϕj = ∑ j λj 〈f, ej〉H ej .
Most kernels used in machine learning are infinite dimensional, i.e., J = N. For convenience, we define Φm := (ϕ1, . . . , ϕm) and Λm = diag(λ1, . . . , λm).

heading: A.3 GENERAL SAMPLE COMPLEXITY AND ASSUMPTIONS ON THE PRODUCT KERNEL
text: In this section, we first consider kernels k0 with scalar input, i.e., X0 ⊆ R. Assume there is a measure µ0 on X0. This will serve as the basis for the more general product kernels in the form of k(x, y) = ∏d j=1 k0(xj , yj) defined over X d 0 .
With Assumptions 1 and 2, we now state the formal version of Theorem 3 by first providing the sample complexity for approximating the partial derivatives. In the next subsection, we will examine how three different kernels satisfy/unsatisfy the Assumptions 1 and 2, and what the value of Nε is. For each case, we will specify µ0 on X0, and the measure on Xd0 is trivially µ = µ d 0. Theorem 4. Suppose {ws}ns=1 are drawn iid from µ0 on X0, where µ0 is the uniform distribution on [−v/2, v/2] for periodic kernels or periodized Gaussian kernels. Let Z := (k0(w 1, ·), k0(w2, ·), . . . , k0(wn, ·)), and g1 = 1l ∑l a=1 γag a 1 : X d 0 → R, where ‖γ‖∞ ≤ c1 and
ga1 (y) = ∂ 0,1k(xa, y) = ha1(y1) d∏ j=2 k0(x a j , yj) with h a 1(·) := ∂0,1k0(xa1 , ·).
Given ε ∈ (0, 1], let Φm = (ϕ1, . . . ϕm) where m = Nε. Then with probability 1− δ, the following holds when the sample size n = max(Nε, 53ε2NεQ 2 ε log 2Nε δ ):
‖g1‖2H ≤ 1 l2 γ>K1γ + 3c1
( 1 + 2 √ NεMε ) ε, (6)
where (K1)a,b = (ha1) >Z(Z>Z)−1Z>hb1 d∏ j=2 k0(x a j , x b j).
Then we obtain the formal statement of sample complexity, as stated in the following corollary, by combining all the coordinates from Theorem 4. Corollary 1. Suppose all coordinates share the same set of samples {ws}ns=1. Applying the results in (6) for coordinates from 1 to d and using the union bound, we have that with sample size n = max(Nε, 5 3ε2NεQ 2 ε log 2Nε δ ), the following holds with probability 1− dδ,
λmax(G >G) ≤ λmax(P̃G) + 3c1 ( 1 + 2 √ NεMε ) ε. (7)
Equivalently, if Nε, Mε and Qε are constants or poly-log terms of ε which we treat as constant, then to ensure λmax(G>G) ≤ λmax(P̃G) + ε with probability 1− δ, the sample size needs to be
n = 15
ε2 c21
( 1 + 2 √ NεMε )2 NεQ 2 ε log
2dNε δ .
Remark 1. The first term on the right-hand side of (7) is explicitly upper bounded by L2 in our training objective. In the case of Theorem 1, the values of Qε, Nε, and Mε lead to a Õ( 1ε2 ) sample complexity. If we further zoom into the dependence on the period v, then note that Nε is almost a universal constant while Mε = √ 2π v (Nε − 1). So overall, n depends on v by 1 v2 . This is not surprising because smaller period means higher frequency, hence more samples are needed. Remark 2. Corollary 1 postulates that all coordinates share the same set of samples {ws}ns=1. When coordinates differ in their domains, we can draw different sets of samples for them. The sample complexity hence grows by d times as we only use a weak union bound. More refined analysis could save us a factor of d as these sets of samples are independent of each other.
Proof of Theorem 4. Let ε′ := (1 + 2 √ mMε)ε. Since 〈 ga1 , g b 1 〉 H = 〈 ha1 , h b 1 〉 H0 ∏d j=2 k0(x a j , x b j)
and ∣∣k0(xaj , xbj)∣∣ ≤ 1, it suffices to show that for all a, b ∈ [l],∣∣∣〈ha1 , hb1〉H0 − (ha1)>Z(Z>Z)−1Z>hb1∣∣∣ ≤ 3ε′.
Towards this end, it is sufficient to show that for any h(·) = ϑx∂0,1k0(x, ·) + ϑy∂0,1k0(y, ·) where x, y ∈ X0 and |ϑx|+ |ϑy| ≤ 1, we have∣∣∣h>Z(Z>Z)−1Z>h− ‖h‖2H0∣∣∣ ≤ ε′. (8)
This is because, if so, then∣∣∣〈ha1 , hb1〉H0 − (ha1)>Z(Z>Z)−1Z>hb1∣∣∣ =
∣∣∣∣12(∥∥ha1 + hb1∥∥2H0 − ‖ha1‖2H0 − ∥∥hb1∥∥2H0)− 12 [(ha1 + hb1)>Z(Z>Z)−1Z>(ha1 + hb1) −(ha1)>Z(Z>Z)−1Z>ha1 − (hb1)>Z(Z>Z)−1Z>hb1
]∣∣ ≤1
2 (4ε′ + ε′ + ε′) = 3ε′.
The rest of the proof is devoted to (8). Since n ≥ m, the SVD ofΛ−1/2m Φ>mZ can be written asUΣV >, where UU> = U>U = V >V = Im (m-by-m identity matrix), and Σ = diag(σ1, . . . , σm). Define
α = n−1/2V U>Λ−1/2m Φ > mh.
Consider the optimization problem o(α) := 12 ‖Zα− h‖ 2 H0 . It is easy to see that its minimal objective value is o∗ := 12 ‖h‖ 2 H0 − 1 2h >Z(Z>Z)−1Z>h. So
0 ≤ 2o∗ = ‖h‖2H0 − h >Z(Z>Z)−1Z>h ≤ 2o(α).
Therefore to prove (8), it suffices to bound o(α) = ‖Zα− h‖H0 . Since √ nΦmΛ 1/2UV >α = ΦmΦ > mh, we can decompose ‖Zα− h‖H0 by
‖Zα− h‖H0 ≤ ∥∥(Z − ΦmΦ>mZ)α∥∥H0 + ∥∥∥(ΦmΦ>mZ −√nΦmΛ1/2m UV >)α∥∥∥H0 (9) + ∥∥ΦmΦ>mh− h∥∥H0 .
The last term ∥∥ΦmΦ>mh− h∥∥H0 is clearly below ε because by Assumption 1 and m = Nε∥∥ΦmΦ>mh− h∥∥H0
≤ |ϑx| ∥∥ΦmΦ>m∂0,1k0(x, ·)− ∂0,1k0(x, ·)∥∥H0 + |ϑy|∥∥ΦmΦ>m∂0,1k0(y, ·)− ∂0,1k0(y, ·)∥∥H0
≤(|ϑx|+ |ϑy|)ε ≤ ε. We will next bound the first two terms on the right-hand side of (9).
(i) By Assumption 1, ∥∥k0(ws, ·)− ΦmΦ>mk0(ws, ·)∥∥H0 ≤ ε, hence ∥∥(Z − ΦmΦ>mZ)α∥∥H0 ≤
ε √ n ‖α‖2. To bound ‖α‖2, note all singular values of V U> are 1, and so Assumption 2 implies that for all i ∈ [m],∣∣∣λ−1/2j 〈ϕj , h〉H0∣∣∣ = ∣∣∣〈ej , h〉H0 ∣∣∣ = ∣∣∣〈ej , ϑx∂0,1k0(x, ·) + ϑy∂0,1k0(y, ·)〉H0∣∣∣ (10) ≤ sup x∈X
∣∣∣〈ej , ∂0,1k(x, ·)〉H0∣∣∣ ≤Mε. As a result, ∥∥(Z − ΦmΦ>mZ)αj∥∥H0 ≤ εn1/2 · n−1/2 ∥∥∥Λ−1/2m Φ>mh∥∥∥ ≤ ε√mMε. (ii) We first consider the concentration of the matrix R := 1nΛ −1/2 m Φ>mZZ >ΦmΛ −1/2 m ∈ Rm×m. Clearly,
E {ws} [Rij ] = E {ws}
[ 1
n n∑ s=1 ei(ws)ej(ws)
] = ∫ ei(x)ej(x) dµ(x) = δij .
By matrix Bernstein theorem (Tropp, 2015, Theorem 1.6.2), we have Pr ( ‖R− Im‖sp ≤ ε ) ≥ 1− δ
when n ≥ O(.). This is because ‖(e1(x), . . . , em(x))‖2 ≤ mQ2ε, ∥∥E{ws}[RR>]∥∥sp ≤ mQ2ε/n, and
Pr ( ‖R− Im‖sp ≤ ε ) ≥ 1− 2m exp
( −ε2
mQ2ε n ( 1 + 23ε )) ≥ 1− 2m exp( −ε25mQ2ε 3n ) ≥ 1− δ,
where the last step is by the definition of n. Since R = 1nUΣ 2U>, this means with probability 1− δ,∥∥ 1
nUΣ 2U> − Im ∥∥ sp ≤ ε. So for all i ∈ [m],
∣∣∣∣ 1nσ2i − 1 ∣∣∣∣ ≤ ε which implies ∣∣∣∣ 1√nσi − 1 ∣∣∣∣ < ε ∣∣∣∣ 1√nσi + 1 ∣∣∣∣−1 ≤ ε. (11)
Moreover, λ1 ≤ 1 since k0(x, x) = 1. It then follows that∥∥∥(ΦmΦ>mZ −√nΦmΛ1/2m UV >)α∥∥∥H0 =
∥∥∥∥ΦmΛ1/2m UΣV > 1√nV U>Λ−1/2m Φ>mh−√nΦmΛ1/2m UV > 1√nV U>Λ−1/2m Φ>mh ∥∥∥∥ H0
= ∥∥∥∥Λ1/2m U( 1√nΣ − Im ) U>Λ−1/2m Φ > mh ∥∥∥∥ 2
(because Φ>mΦm = Im)
≤ √ λ1 max
i∈[m] ∣∣∣∣ 1√nσi − 1 ∣∣∣∣ ∥∥∥Λ−1/2m Φ>mh∥∥∥ 2
≤ε √ mMε (by (11), (10), and λ1 ≤ 1).
Combining (i) and (ii), we arrive at the desired bound in (6).
Proof of Corollary 1. Since P̃G approximates G>G only on the diagonal, P̃G −G>G is a diagonal matrix which we denote as diag(δ1, . . . , δd). Let u ∈ Rd be the leading eigenvector of P̃G. Then
λmax(P̃G)− λmax(G>G) ≤ u>P̃Gu− u>G>Gu = u>(P̃G −G>G)u = ∑ j δju 2 j
(by (6)) ≤ 3c1 ( 1 + 2 √ NεMε ) ε.
The proof is completed by applying the union bound and rewriting the results.

heading: A.4 CASE 1: CHECKING ASSUMPTIONS 1 AND 2 ON PERIODIC KERNELS
text: Periodic kernels on X0 := R are translation invariant, and can be written as k0(x, y) = κ(x − y) where κ : R→ R is a) periodic with period v; b) even, with κ(−t) = κ(t); and c) normalized with κ(0) = 1. A general treatment was given by Williamson et al. (2001), and an example was given by David MacKay in MacKay (1998):
k0(x, y) = exp
( − 1 2σ2 sin (π v (x− y) )2) . (12)
We define µ0 to be a uniform distribution on [−v2 , v 2 ], and let ω0 = 2π/v.
Since κ is symmetric, we can simplify the Fourier transform of κ(t)δv(t), where δv(t) = 1 if t ∈ [−v/2, v/2], and 0 otherwise:
F (ω) = 1√ 2π ∫ v/2 −v/2 κ(t) cos(ωt) dt.
It is now easy to observe that thanks to periodicity and symmetry of κ, for all j ∈ Z,
1
v ∫ v/2 −v/2 k0(x, y) cos(jω0y) dy = 1 v ∫ v/2 −v/2 κ(x− y) cos(jω0y) dy
= 1
v ∫ x+v/2 x−v/2 κ(z) cos(jω0(x− z)) dz (note cos(jω0(x− z)) also has period v)
= 1
v ∫ v/2 −v/2 κ(z)[cos(jω0x) cos(jω0z) + sin(jω0x) sin(jω0z)) dz (by periodicity)
= 1
v cos(jω0x) ∫ v/2 −v/2 κ(z) cos(jω0z) dz (by symmetry of κ)
=
√ 2π
v F (jω0) cos(jω0x).
And similarly,
1
v ∫ v/2 −v/2 k0(x, y) sin(jω0y) dy = √ 2π v F (jω0) sin(jω0x).
Therefore the eigenfunctions of the integral operator Tk are
e0(x) = 1, ej(x) := √ 2 cos(jω0x), e−j(x) := √ 2 sin(jω0x) (j ≥ 1)
and the eigenvalues are λj = √ 2π v F (jω0) for all j ∈ Z with λ−j = λj . An important property our proof will rely on is that
e′j(x) = −jω0e−j(x), for all j ∈ Z. Applying Mercer’s theorem in (5) and noting κ(0) = 1, we derive ∑ j∈Z λj = 1.
Checking the Assumptions 1 and 2. The following theorem summarizes the assumptions and conclusions regarding the satisfaction of Assumptions 1 and 2. Again we focus on the case of X ⊆ R. Theorem 1. Suppose the periodic kernel with period v has eigenvalues λj that satisfies
λj(1 + j) 2 max(1, j2)(1 + δ(j ≥ 1)) ≤ c6 · c−j4 , for all j ≥ 0, (13)
where c4 > 1 and c6 > 0 are universal constants. Then Assumption 1 holds with
Nε = 1 + 2 bnεc , where nε := logc4 ( 2.1c6 ε2 max ( 1, v2 4π2 )) . (14)
In addition, Assumption 2 holds with Qε = √ 2 and Mε = 2 √ 2π v bnεc = √ 2π v (Nε − 1).
For example, if we set v = π and σ2 = 1/2 in the kernel in (12), elementary calculation shows that the condition (13) is satisfied with c4 = 2 and c6 = 1.6.
Proof of Theorem 1. First we show that h(x) := ∂0,1k0(x0, x) is in H0 for all x0 ∈ X0. Since k0(x0, x) = ∑ j∈Z λjej(x0)ej(x), we derive
h(x) = ∑ j∈Z λjej(x0)∂ 1ej(x) = ∑ j∈Z λjej(x0)(−jω0e−j(x)) = ω0 ∑ j∈Z λjje−j(x0)ej(x). (15)
h(x) is inH if the sequence λjje−j(x0)/ √ λj is square summable. This can be easily seen by (13):
ω−20 ‖h‖ 2 H0 = ∑ j λjj 2e2−j(x0) = ∑ j∈Z λjj 2e2−j(x0)
= ∑ j∈Z λjj 2e2−j(x0) = λ0 + 2 ∑ j≥1 j2λj ≤ 2c4c5 c4 − 1 .
Finally to derive Nε, we reuse the orthonormal decomposition of h(x) in (15). For a given set of j values A where A ⊆ Z, we denote as ΦA the “matrix” whose columns enumerate the ϕj over j ∈ A. Let us choose
A := { j : λj max(1, j 2)(1 + j2)(1 + δ(j ≥ 1)) ≥ min(1, w−20 ) ε2
2.1
} .
If j ∈ A, then −j ∈ A. LettingN0 = {0, 1, 2, . . .}, we note ∑ j∈N0
1 1+j2 ≤ 2.1. So∥∥h− ΦAΦ>Ah∥∥2H0 = w20 ∑
j∈Z\A
λjj 2e2−j(x0)
= w20 ∑
j∈N0\A
λjj 2 [ (e2j (x) + e 2 −j(x))δ(j ≥ 1) + δ(j = 0) ] = w20
∑ j∈N0\A λjj 2(1 + δ(j ≥ 1))
= w20 ∑
j∈N0\A
{ λjj
2(1 + j2)(1 + δ(j ≥ 1)) 1 1 + j2
}
≤ ε 2
2.1 ∑ j∈N0 1 1 + j2 = ε2 2.1 ∑ j∈N0 1 1 + j2 ≤ ε2.
Similarly, we can bound ∥∥k0(x0, ·)− ΦAΦ>Ak0(x0, ·)∥∥H0 by∥∥k0(x0, ·)− ΦAΦ>Ak0(x0, ·)∥∥2H0
= ∑ j∈Z\A λje 2 j (x0) ≤ ∑ j∈Z\A λj max(1, j 2)e2j (x0)
= ∑
j∈N0\A
λαmax(1, j 2)[ ( e2j (x) + e 2 −j(x) ) δ(j ≥ 1) + δ(j = 0)]
= ∑
j∈N0\A
{ λj max(1, j
2)(1 + j2)(1 + δ(j ≥ 1)) 1 1 + j2 } ≤ 1
2.1 ε2 ∑ j∈N0 1 1 + j2 ≤ ε2.
To upper bound the cardinality of A, we consider the conditions for j /∈ A. Thanks to the conditions in (13), we know that any j satisfying the following relationship cannot be in A:
c6 · c−|j|4 < min(1, w −2 0 )
ε2 2.1 ⇔ c−|j|4 <
1
2.1 · c6 min
( 1, 4π2
v2
) ε2.
So A ⊆ {j : |j| ≤ nε}, which yields the conclusion (14). Finally Qε ≤ √
2, and to bound Mε, we simply reuse (15). For any j with |j| ≤ nε,∣∣〈h, ej〉H∣∣ ≤ ω0 |je−j(x0)| ≤ 2πv √2 bnεc = √ 2π v (Nε − 1).

heading: A.5 CASE 2: CHECKING ASSUMPTIONS 1 AND 2 ON GAUSSIAN KERNELS
text: Gaussian kernels k(x, y) = exp(−‖x− y‖2 /(2σ2)) are obviously product kernels with k0(x1, y1) = κ(x1 − y1) = exp(−(x1 − y1)2/(2σ2)). It is also translation invariant. The spectrum of Gaussian kernel k0 on R is known; see, e.g., Chapter 4.3.1 of Rasmussen & Williams (2006) and Section 4 of Zhu et al. (1998). Let µ be a Gaussian distributionN (0, σ2). Setting ε2 = α2 = (2σ2)−1
in Eq 12 and 13 of E Fasshauer (2011), the eigenvalue and eigenfunctions are (for j ≥ 0):
λj = c −j−1/2 0 , where c0 =
1 2 (3 +
√ 5)
ej(x) = 51/8
2j/2 exp
( − √
5− 1 4 x2 σ2 ) 1√ j! Hj ( 4 √ 1.25 x σ ) ,
where Hj is the Hermite polynomial of order j.
Although the eigenvalues decay exponentially fast, the eigenfunctions are not uniformly bounded in the L∞ sense. Although the latter can be patched if we restrict x to a bounded set, the above closed-form of eigen-pairs will no longer hold, and the analysis will become rather challenging.
To resolve this issue, we resort to the period-ization technique proposed by Williamson et al. (2001). Consider κ(x) = exp(−x2/(2σ2)) when x ∈ [−v/2, v/2], and then extend κ to R as a periodic function with period v. Again let µ be the uniform distribution on [−v/2, v/2]. As can be seen from the discriminant function f = 1l ∑l i=1 γik(x
i, ·), as along as our training and test data both lie in [−v/4, v/4], the modification of κ outside [−v/2, v/2] does not effectively make any difference. Although the term ∂0,1k0(xa1 , w 1 1) in (13) may possibly evaluate κ outside [−v/2, v/2], it is only used for testing the gradient norm bound of κ. With this periodized Gaussian kernel, it is easy to see that Qε = √
2. If we standardize by σ = 1 and set v = 5π as an example, it is not hard to see that (13) holds with c4 = 1.25 and c6 = 50. The expressions of Nε and Mε then follow from Theorem 1 directly.

heading: A.6 CASE 3: CHECKING ASSUMPTIONS 1 AND 2 ON NON-PRODUCT KERNELS
text: The above analysis has been restricted to product kernels. But in practice, there are many useful kernels that are not decomposable. A prominent example is the inverse kernel: k(x, y) = (2−x>y)−1. In general, it is extremely challenging to analyze eigenfunctions, which are commonly not bounded (Lafferty & Lebanon, 2005; Zhou, 2002), i.e., supi→∞ supx |ei(x)| = ∞. The opposite was (incorrectly) claimed in Theorem 4 of Williamson et al. (2001) by citing an incorrect result in König (1986, p. 145), which was later corrected by Zhou (2002) and Steve Smale. Indeed, uniform boundedness is not known even for Gaussian kernels with uniform distribution on [0, 1]d Lin et al. (2017), and (Minh et al., 2006, Theorem 5) showed the unboundedness for Gaussian kernels with uniform distribution on the unit sphere when d ≥ 3. Here we only present the limited results that we have obtained on the eigenvalues of the integral operator of inverse kernels with a uniform distribution on the unit ball. The analysis of eigenfunctions is left for future work. Specifically, in order to drive the eigenvalue λi below ε, i must be at least ddlog2 1 εe+1. This is a quasi-quadratic bound if we view d and 1/ε as two large variables.
It is quite straightforward to give an explicit characterization of the functions in H. The Taylor expansion of z−1 at z = 2 is 12 ∑∞ i=0(− 1 2 ) ixi. Using the standard multi-index notation with
α = (α1, . . . , αd) ∈ (N ∪ {0})d, |α| = ∑d i=1 αi, and x α = xα11 . . . x αd d , we derive
k(x,y) = 1
2− x>y =
1
2 ∞∑ k=0 ( −1 2 )k (−x>y)k = ∞∑ k=0 2−k−1 ∑ α:|α|=k Ckαx αyα
= ∑ α 2−|α|−1C |α|α x αyα,
where Ckα = k!∏d
i=1 αi! . So we can read off the feature mapping for x as
ϕ(x) = {wαxα : α}, where wα = 2− 1 2 (|α|+1)C |α|α ,
and the functions inH are
H = { f =
∑ α ϑαwαx α : ‖ϑ‖`2 <∞
} . (16)
Note this is just an intuitive “derivation” while a rigorous proof for (16) can be constructed in analogy to that of Theorem 1 in Minh (2010).

heading: A.7 BACKGROUND OF EIGENVALUES OF A KERNEL
text: We now use (16) to find the eigenvalues of inverse kernel.
Now specializing to our inverse kernel case, let us endow a uniform distribution over the unit ball B: p(x) = V −1d where Vd = π d/2Γ (d2 + 1) −1 is the volume of B, with Γ being the Gamma function. Then λ is an eigenvalue of the kernel if there exists f = ∑ α ϑαwαx
α such that∫ y∈B k(x,y)p(y)f(y) dy = λf(x). This translates to
V −1d ∫ y∈B ∑ α w2αx αyα ∑ β ϑβwβy β dy = λ ∑ α ϑαwαx α, ∀ x ∈ B.
Since B is an open set, that means wα ∑ β wβqα+βϑβ = λϑα, ∀ α,
where
qα = V −1 d ∫ y∈B yα dy =  2 ∏d i=1 Γ ( 1 2αi+ 1 2 ) Vd·(|α|+d)·Γ ( 1 2 |α|+ d 2 ) if all αi are even 0 otherwise .
In other words, λ is the eigenvalue of the infinite dimensional matrix Q = [wαwβqα+β]α,β,

heading: A.8 BOUNDING THE EIGENVALUES
text: To bound the eigenvalues of Q, we resort to the majorization results in matrix analysis. Since k is a PSD kernel, all its eigenvalues are nonnegative, and suppose they are sorted decreasingly as λ1 ≥ λ2 ≥ . . .. Let the row corresponding to α have `2 norm rα, and let them be sorted as r[1] ≥ r[2] ≥ . . .. Then by Schneider (1953); Shi & Wang (1965), we have
n∏ i=1 λi ≤ n∏ i=1 r[i], ∀ n ≥ 1.
So our strategy is to bound rα first. To start with, we decompose qα+β into qα and qβ via CauchySchwartz:
q2α+β = V −2 d (∫ y∈B yα+β dy )2 ≤ V −2d ∫ y∈B y2α dy · ∫ y∈B y2β dy = q2αq2β.
To simplify notation, we consider without loss of generality that d is an even number, and denote the integer b := d/2. Now Vd = πb/b!. Noting that there are ( k + d− 1
k
) values of β such that
|β| = k, we can proceed by (fix below by changing ( k + d k ) into ( k + d− 1 k ) , or no need
because the former upper bounds the latter)
r2α = w 2 α ∑ β w2βq 2 α+β ≤ w2αq2α ∑ β w2βq2β = w 2 αq2α ∞∑ k=0 2−k−1 ∑ β:|β|=k Ckβq2β
≤ w2αq2α ∞∑ k=0 2−k−1 ( k + d d ) max |β|=k Ckβq2β
= w2αq2α ∞∑ k=0 2−k−1 ( k + d d ) max |β|=k k!∏d i=1 βi! · 2 ∏d i=1 Γ (βi + 1 2 ) Vd · (2k + d) · Γ (k + d2 )
= w2αq2αV −1 d ∞∑ k=0 2−k ( k + d d ) k! (2k + d)Γ (k + d2 ) · max |β|=k d∏ i=1 Γ (βi + 1 2 ) βi!
< w2αq2α · b! πbd! · ∞∑ k=0 2−k−1 (k + d)! (k + b)! (since Γ (βi + 12 ) < Γ (βi + 1) = βi!).
The summation over k can be bounded by
∞∑ k=0 2−k−1 (k + d)! (k + b)! = 1 2 b! ( 2d + ( d b )) ≤ 1 2 ( b!2d + 2b ) ≤ b!2d,
where the first equality used the identity ∑∞ k=1 2 −k ( d+ k b ) = 2d. Letting l := |α|, we can
continue by
r2α < w 2 αq2α ·
b!
πbd! b!2d = 2−l−1 l!∏d i=1 αi!
2 ∏d i=1 Γ ( αi + 1 2 ) Vd · (2l + d) · Γ (l + b) (b!)22d πbd!
≤ 2−l+dπ−2b l!(b!) 3
d!(l + b− 1)!(2l + d) (since Γ (αi + 12 ) < Γ (αi + 1) = αi!)
≤ 2−l+b−1π−2b ( l + b l )−1 (since (b!)2 d! ≤ 2−b).
This bound depends on α, not directly on α. Letting nl = ( l + d− 1
l
) and NL = ∑L l=0 nl =(
d+ L L
) , it follows that
L∑ l=0 lnl = L∑ l=1 l(l + d)! d! · l! = (d+ 1) L∑ l=1
(l + d)!
(d+ 1)!(l − 1)!
=(d+ 1) L∑ l=1 ( l + d d+ 1 ) = (d+ 1) ( L+ d+ 1 d+ 2 ) .
Now we can bound λNL by
λNLNL ≤ NL∏ i=1 λi ≤ L∏ l=0
( 2−l+b−1π−2b ( l + b l )−1)nl
⇒ log λNL ≤ N−1L L∑ l=0 nl ( −(l − b+ 1) log 2− 2b log π − log ( l + b l ))
≤ −N−1L · log 2 · L∑ l=0 lnl (since log 2 < 2 log π as the coefficients of b)
= − ( d+ L+ 1 d+ 1 )−1 · log 2 · (d+ 1) ( d+ L+ 1 d+ 2 ) = −d+ 1
d+ 2 L log 2
≈ −L log 2 ⇒ λNL ≤ 2−L.
This means that the eigenvalue λi ≤ ε provided that i ≥ NL where L = ⌈ log2 1 ε ⌉ . SinceNL ≤ dL+1, that means it suffices to choose i such that
i ≥ ddlog2 1 εe+1.
This is a quasi-polynomial bound. It seems tight because even in Gaussian RBF kernel, the eigenvalues follow the order of λα = O(c−|α|) for some c > 1 (Fasshauer & McCourt, 2012, p.A742).

heading: EXPERIMENTS

heading: A.9 EFFICIENCY OF ENFORCING LIPSCHITZ CONSTANT BY DIFFERENT METHODS
text: The six different ways to train SVMs with Lipschitz regularisation are summarized in Algorithm 1. Figure 7 plots how fast the regularisation on gradient norm becomes effective when more and more points w are added to the constraint set. We call them “samples” although it is not so random in the greedy method, modulo the random initialization of BFGS within the greedy method. The horizontal axis is the loop index i in Algorithm 1, and the vertical axis is L(i) therein, which is the estimation of the Lipschitz constant of the current solution f (i). We used 400 random examples (200 images of digit 1 and 200 images of digit 0) in the MNIST dataset and set L = 3 and RKHS norm ‖f‖H ≤ ∞ for all algorithms. Inverse kernel is used, hence no results are shown for coordinate-wise Nyström.
Clearly the Nyström algorithm is more efficient than the Brute-force algorithm, and the greedy method significantly reduces the number of samples for both algorithms. In fact, Nyström with greedy selection eventually fell below the prespecified L, because of the gap in (9).

heading: A.10 EXTENSION TO `∞-NORM ATTACKS FOR OUR KERNEL BASED METHOD
text: We now extend our kernel based approach to `∞ norm ball attacks. Since most multiclass losses are 1-Lipschitz continuous with respect to `∞ norm on (f1(x), . . . , f10(x)), we will seek
sup x∈X sup u:‖u‖∞≤1
∥∥∥[g1(x), . . . , g10(x)]>u∥∥∥ ∞ ≤ L, where gc(x) := ∇f c(x).
The left-hand side (LHS) can be bounded by
LHS = sup x∈X max 1≤c≤10 ‖gc(x)‖1 ≤ max1≤c≤10 sup‖ϕ‖H≤1 ∥∥G>c ϕ∥∥1 . Given the Nyström approximation G̃c of Gc, we can enforce the convex constraint of
max 1≤c≤10 sup ‖v‖2≤1
∥∥∥G̃>c v∥∥∥ 1 ≤ L.

heading: A.11 MORE RESULTS ON CROSS-ENTROPY ATTACKS
text: A.12 VISUALIZATION OF GRADIENT
A gradient-based attacker tries to decrease the targeted loss by following the negative gradient in Figure 10b, i.e., reduce the pixel value in red area and increase pixel value in blue area.
In order to verify that the robustness of Gauss-Lip is not due to obfuscated gradient, we visualised “large perturbation” adversarial examples, with the `2 norm upper bounded by 6. Figure 11 shows how the PGD attacker uses the gradients to perturb the images step by step. At the end of PGD, there are 46 cases where the original image was successfully attacked, i.e., turned into the target class. This is over 50% of the total of 90 cases, and the resulting images look realistic.
iteration 1 iteration 2 iteration 3 iteration 4 iteration 5
iteration 6 iteration 7 iteration 8 iteration 9 iteration 10
prediction of Gauss-Lip
Figure 11: Gradients and perturbed images at each iteration in a 10-step PGD attack using (targeted) cross-entropy approximation, with the `2 norm upper bounded by 6. Here the classifier is Gauss-Lip (σ = 2). The table in the bottom right presents the final predictions of our trained Gauss-Lip on the perturbed images.
To further look into the attack result, we increased PGD to 100 iterations. As shown in Figure 12, now the number of misclassified cases (i.e., unsuccessful attacks that failed to turn an image into the targeted class) drops from 46 to 22, out of 90 cases. The final images are quite realistic. We will further study these remaining cases in the future.
In the above experiments for Figures 11, 12, and 14, PGD was run on the cross-entropy objective. For example, the row corresponding to class 4 tries to promote the likelihood of the target class 4. Naturally the diagonal is not meaningful, hence left empty.
We further ran PGD for 100 iterations on C&W approximation (an untargeted attack used in Figure 3), and the resulting images after every 10 iterations are shown in Figure 13. Here 9 out of 10 images were eventually turned into a different but untargeted class, and the final images are very realistic.
Another random set of images. To test if the above result is due to the particularly hard images selected, we randomly selected another set of images and its results for 100-step PGD on crossentropy objective and C&W objective are shown in Figures 14 and 15, respectively. Interestingly, C&W attack succeeds on all these images, and cross-entropy attack was only unsuccessful in turning 0 into 1.
Please note that despite the commonality in using the cross-entropy objective, the setting of targeted attack in Figures 11, 12, and 14 is not comparable to that in Figure 8a, where to enable a batch test mode, an untargeted attacker was employed by increasing the cross-entropy loss of the correct class, i.e., decreasing the likelihood of the correct class. This is a common practice.

Review summaries:
Through the lens of Distributional Robust Risk (DRR), this work draws a link between adversarial robustness and Lipschitz constant regularisation. The authors first provide an upper bound of the DRR (with a Wasserstein ball as the ambiguity set) in terms of the true risk and the Lipschitz constant of the loss function under the current model. They show that the standard adversarial risk can be upper bounded by the DRR, emphasizing that the Lipschitz constant regularised loss can be used as a proxy for adversarially robust training.

This paper uses results from distributional robustness to provide bounds of p-norm-constrained adversarial risk which depend on the Lipschitz constant of the underlying classifier. The bulk of the paper focuses on sample-efficient mechanisms to approximate the Lipschitz constants of kernel methods so that a constraint on this Lip constant can be enforced during training. Empirically, the kernel methods are compared to existing deep learning approaches and are shown to be competitive at this scale.