Title: CERTIFYING DISTRIBUTIONAL ROBUSTNESS USING LIPSCHITZ REGULARISATION

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.

1. INTRODUCTION:

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) .

2. PRELIMINARIES:

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 L 0 (X, Y ) denote the Borel measurable mappings X → Y , and L 1 (X, µ) denote the Borel functions f ∈ L 0 (X, R) with |f | dµ < ∞ for µ ∈ P(X).

For two measures µ, ν ∈ P(Ω) the set of (µ, ν)-couplings is

Π(µ, ν) := {π ∈ P(Ω × Ω) | µ = π( • , dω), ν = π(dω, • )}. Let c ∈ L 0 (Ω × Ω, R).

The c-transportation cost of µ, ν ∈ P(Ω), and c-transportation cost ball of radius r ≥ 0, centred at µ ∈ P(Ω) are respectively

EQUATION

The least c-Lipschitz constant of f (cf. Cranko et al., 2019) is the infimum over L ≥ 0 satisfying (2), and is denoted by lip c (f ), so that when (X, d) is a metric space lip d (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).

3. CERTIFICATE FOR DISTRIBUTIONAL ROBUSTNESS:

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

) := sup x∈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 = R n 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 ∈ L 1 (X, µ) is upper semicontinuous. Then for all r ≥ 0,

EQUATION

)

The tightness of the bound can be quantified as follows. Let ∆(µ)

EQUATION

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 X 0 ⊆ X, with the collection of probability distributions supported on X 0 . Proposition 1. Assume X is a separable Fréchet space with X 0 ⊆ X. Assume c : X → R≥0 is sublinear and continuous, and

f ∈ µ∈P(X0) L 1 (X, µ) is upper semicontinuous, has lip c ( f ) < ∞,

and attains its maximum on X 0 . Then for all r ≥ 0 with 1/0 := ∞,

sup µ∈P(X0) ∆(µ) = r lip c ( f ) -lip c (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 ) = c X (x -x

) whenever y = y and c(x, y, x , y ) = ∞ whenever y = y , where c X : X → R is symmetric, sublinear, and continuous, and f ∈ L 1 (X × Y, µ) is upper semicontinuous. Then for all r ≥ 0 there is (3).

To see the tightness of the bound, if lip c ( 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).

3.1. DISTRIBUTIONAL ROBUSTNESS AS ADVERSARIAL ROBUSTNESS:

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

EQUATION

where in the equality we extend d to a metric on Ω := X × Y with d((x, y), (x , y )) := d(x, x ) + ∞ y = y .

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 ∈ L 0 (X, R ≥0 ) | g dµ ≤ r}. Then for f ∈ L 0 (Ω, 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 g r ≡ 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 γ i k(x i , •), with x i 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.

4. PROVABLE LIPSCHITZ REGULARISATION FOR KERNEL METHODS:

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 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 ⊆ R d , 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) = 1 l l a=1 γ a k(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 sup x∈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) ,

EQUATION

)

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 {w s } n s=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 g j (x) := ∂ j f (x) be the partial derivative with respect to the j-th coordinate of x, and ∂ i,j k(x, y) be the partial derivative to x i and y j . i or j being 0 means no derivative. Assuming sup x∈X k(x, x) = 1 and g j ∈ H (true for various kernels considered by Assumptions 1 and 2 below), we get a new bound

sup x∈X ∇f (x) 2 2 = sup x∈X d j=1 g j , k(x, •) 2 H ≤ sup ϕ: ϕ H =1 d j=1 g j , ϕ 2 H = λ max (G G), (9)

where λ max evaluates the maximum eigenvalue, and G := (g 1 , . . . , g d ). 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 g i , g j H . Interestingly, λ max (G G) delivers significantly lower (i.e., tighter) value in approximating the Lipschitz constant sup x∈X ∇f (x) 2 , compared with f H max z>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 {w s } 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 g j can be approximated by gj , where gj ∈ R n is the Nyström approximation (Drineas & Mahoney, 2005; Williams & Seeger, 2000) :

a gj := K -1/2 (g j (w 1 ), . . . , g j (w n )) = (Z Z) -1/2 Z g j (noting g j (w 1 ) = g j , k(w 1 , •) H ) (10) where K := [k(w i , w i )] i,i , Z := (k(w 1 , •), k(w 2 , •), . . . , k(w n , •)), G := (g 1 , . . . , gd ).

So to ensure λ max (G G) ≤ L 2 + ε, intuitively we can resort to enforcing λ max ( G G) ≤ L 2 , 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.

4.1. A COORDINATE-WISE NYSTRÖM APPROXIMATION FOR PRODUCT KERNELS:

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 k 0 (x j , y j ) where X = X d 0 and k 0 is a base kernel on X 0 . Let the RKHS of k 0 be H 0 , and let µ 0 be a finite Borel measure with supp

[µ 0 ] = X 0 . Periodic kernels have k 0 (x j , y j ) = exp -sin π v (x j -y j ) 2 /(2σ 2 ) .

The key benefit of this decomposition is that the derivative ∂ 0,1 k(x, y) can be written as ∂ 0,1 k 0 (x 1 , y 1 ) d j=2 k 0 (x j , y j ). Since k 0 (x j , y j ) can be easily dealt with, approximation will be needed only for ∂ 0,1 k 0 (x 1 , y 1 ). Applying this idea to

EQUATION

g 1 , g 2 H = l -2 l a,b=1 γ a γ b ∂ 0,1 k 0 (x a 1 , x b 1 )∂ 0,1 k 0 (x b 2 , x a 2 ) d j=3 k 0 (x a j , x b j ).

So the off-diagonal entries of G G can be computed exactly. To approximate the diagonal, we sample

{w 1 1 , . . . , w n 1 } from µ 0 , set Z 1 = (k 0 (w 1 1 , •), . . . , k 0 (w n 1 , •))

, and apply Nyström approximation:

EQUATION

and analogously for Z 1 ∂ 0,1 k 0 (x b 1 , •). We will denote this approximation of G G as PG . Clearly,

λ max ( PG ) ≤ L 2 is a convex constraint on γ, based on i.i.d. samples {w s j : s ∈ [n], j ∈ [d]} from µ 0 .

It is now important to analyse how many samples w s j are needed, such that

λ max ( PG ) ≤ L 2 =⇒ λ max (G G) ≤ L 2 + ε with high probability.

4.2. GENERAL SAMPLE COMPLEXITY AND ASSUMPTIONS ON THE PRODUCT KERNEL:

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 k 0 with a scalar input, i.e., X 0 ⊆ R. Recall from Steinwart & Christmann (2008, Chapter 4 ) that the integral operator for k 0 and µ 0 is defined by

T k0 = I • S : L 2 (X 0 , µ 0 ) → L 2 (X 0 , µ 0 )

where S :

EQUATION

I: C(X 0 ) → L 2 (X 0 ; µ 0 ) is the inclusion operator. By the spectral theorem, if T k0 is compact, then there is an at most countable orthonormal set {ẽ j } j∈J of L 2 (X 0 , µ 0 ) and {λ j } j∈J with λ 1 ≥ λ 2 ≥ . . . > 0 such that T k0 f = j∈J λ j f, ẽj L2(X0;µ0) ẽj for all f ∈ L 2 (X 0 , µ 0 ).

It is easy to see that ϕ j := λ j e j is an orthonormal basis of H 0 (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 k 0 (x, •) and ∂ 0,1 k 0 (x, •) "concentrates" on the leading eigenfunctions. Assumption 1. Suppose k 0 (x, x) = 1 and ∂ 0,1 k 0 (x, •) ∈ H 0 for all x ∈ X 0 . For all ε > 0, there exists N ε ∈ N such that the tail energy of ∂ 0,1 k 0 (x, •) beyond the N ε -th eigenpair is less than ε, uniformly for all x ∈ X 0 . That is, denoting Φ m := (ϕ 1 , . . . , ϕ m ),

N ε := inf m ∂ 0,1 k 0 (x, •) -Φ m Φ m ∂ 0,1 k 0 (x, •) H0 < ε for all x ∈ X 0 and k 0 (x, •) -Φ m Φ m k 0 (x, •) H0 < ε for all x ∈ X 0 < ∞.

The second assumption asserts the smoothness and range of eigenfunctions in a uniform sense.

Assumption 2. Under Assumption 1, {e j (x) : j ∈ N ε } is uniformed bounded over x ∈ X 0 , and the RKHS inner product of ∂ 0,1 k 0 (x, •) with {e j : j ∈ N ε } is also uniformly bounded over x ∈ X 0 :

EQUATION

Q ε := sup x∈X0 max j∈[Nε] |e j (x)| < ∞.

Theorem 3. Suppose k 0 , X 0 , and µ 0 satisfy Assumptions 1 and 2. Let {w s j : 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 ( PG ) ≤ L 2 , the Lipschitz condition λ max (G G) ≤ L 2 + ε is met with probability 1 -δ, as long as n ≥ Θ 1 ε 2 N 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 Õ(1) values of N ε , M ε , and Q ε . 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.

5. EXPERIMENTAL RESULTS:

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 {w s }, 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 ( PG ) 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 {w s } (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 a k(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 → (f 1 (x), . . . , f 10 (x)) by 

EQUATION

)

where

G(x) := [∇f 1 (x), • • • , ∇f 10 (x)] = [G 1 k(x, •), . . . , G 10 k(x, •)] with G c := (g c 1 , . . . , g c d )

. 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. Gauss-Lip 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.

Return if L (i) := max x∈X ∇f (i) (x) falls below L.

Finding the exact argmax x∈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) = 1 l l a=1 γ a k(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(x a , •), denoted as φ(x a ) ∈ R 1000 . Then f (x) can be written as 1 l l a=1 γ a φ(x a ) φ(x). Defining w = 1 l 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 g j 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 γ 1 i k(x i , •); . . . ; n i=1 α 10 i k(x i , •) s.t. sup ϕ H ≤1 λ max 10 c=1 G c ϕϕ G c ≈ sup v 2 ≤1 λ max 10 c=1 G c vv Gc ≤ L 2 ,

where (F (x), y) is the Crammer & Singer loss, and the constraint is derived from ( 14) by using its Nyström approximation Gc = [g c 1 , . . . , gc 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 Gc = sup v 2 ≤1, u 2 ≤1 u 10 c=1 G c vv Gc 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 g 1 (x), . . . , g 10 (x) ∞ = sup x∈X max 1≤c≤10 g c (x) 1 = max 1≤c≤10 sup x∈X g c (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

EQUATION

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 .

A PROOFS OF RESULTS:

The following appendix contains the complete set of proofs and auxiliary results.

PROOFS FOR §3: CERTIFICATE FOR DISTRIBUTIONAL ROBUSTNESS:

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 ∈ L 1 (Ω, µ) is upper semicontinuous. Then for all r ≥ 0 there is

EQUATION

)

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) ⊆ ∂λc(0).

Proof. Because f is closed convex, it is equal to its biconjugate (Zȃlinescu, 2002, Thm. 2.3 .3), because c is sublinear and lower semicontinuous λc(x) = sup x * ∈ ∂λc(0) x, x * for all x ∈ X (Zȃ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 (Zȃlinescu, 2002, Thm. 2.4.9) , and so we can apply a minimax theorem (Zȃ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 * = x * = sup x * ∈ ∂f (X) y, x * -f * (x * ) x * ∈ ∂λc(0) ∞ x * / ∈ ∂λc(0) = f (y) ∂f (X) ⊆ ∂λc(0) ∞ ∂f (X) ⊆ ∂λ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 ∈ L 1 (X, µ) is upper semicontinuous. Then for all r ≥ 0,

DRR := sup ν∈Bc(µ,r) risk (f, ν) ≤ r lip c ( f ) + risk (f, µ).

The tightness of the bound can be quantified as follows. Let ∆(µ)

:= r lip c ( f ) + risk (f, µ) - sup ν∈Bc(µ,r) risk (f, ν). If lip c (f ) < ∞ then ∆(µ) ≤ r lip c ( f ) -lip c (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λ + λc f dµ ≤ inf λ≥lip c ( f ) rλ + λc f dµ = r lip c ( f ) + f dµ.

Observing that co f ≤ f , applying Lemma 3 and Lemma 4 we find for all x ∈ X

EQUATION

Similarly, for all x ∈ X there is

EQUATION

Then, using (1) and ( 2) we find

r lip c ( f ) + f dµ -inf λ∈[0,∞) rλ - λc f dµ = r lip c ( f ) + sup λ∈[0,∞) f -λc f -λr dµ ≤ r lip c ( f ) + sup λ∈[0,∞) f -λc f -λr dµ (1),(2) ≤ r lip c ( f ) + min ( f -co f ) dµ -r lip c (co f ), 0 .

The proof is complete.

Proposition 1. Assume X is a separable Fréchet space with X 0 ⊆ X. Assume c : X → R≥0 is sublinear and continuous, and f ∈ µ∈P(X0) L 1 (X, µ) is upper semicontinuous, has lip c ( f ) < ∞, and attains its maximum on X 0 . Then for all r ≥ 0 with 1/0 := ∞,

sup µ∈P(X0) ∆(µ) = r lip c ( f ) -lip c (co f ) - 1 r ρ( f ) + .

Proof. Let x 0 ∈ X 0 be be the point at which f (x 0 ) = sup x∈X0 f (x). Then

EQUATION

)

Then there is

r lip c ( f ) (3) ≤ sup µ∈P(X0) ∆(µ) (4) ≤ r lip c ( f ) -max lip c (co f ) - 1 r ρ( f ), 0 ≤ r lip c ( f ),

which completes the proof.

Remark 2. When f : R n → R satisfies f ≡ ∞ and f is minorised by an affine function, there is (cf. Hiriart-Urruty & Lemaréchal, 2010, Prop. 1.5.4)

∀x ∈ R n : co f (x) = inf    i∈[n+1] α i f (x i ) | i∈[n+1] α i = 1, x = i∈[n+1] α i x i    ,

where the infimum is over all sequences (α i ) i∈[n+1] and (x i ) i∈[n+1] ⊆ R n satisfying the conditions above. Consequentially there is the common expression

ρ(f ) = sup    f i∈[n+1] α i x i - i∈[n+1] α i f (x i ) | (α i , x i ) i∈[n+1] ⊆ R ≥0 × R n , 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

(f i ) i∈N ⊆ A µ (r) := f ∈ L 0 (Ω, Ω) | c d(id, f ) # µ ≤ r with (f i ) # µ * ν . Proof. Let P (µ, ν) := {f ∈ L 0 (X, X) | f # µ = ν}.

Since µ is non-atomic and c is continuous Pratelli (2007, Thm. B) shows

∀ν ∈ P(Ω) : inf f ∈P (µ,ν) c d(id, f ) # µ = cost c (µ, ν). Let r := cost c (µ, ν ), 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 cost c (µ, ν(0)) = 0 and cost c (µ, ν(1)) = r ,

and because cost c metrises the weak * topology on P(Ω) (Villani, 2008, Thm. 6 .9), the mapping u → cost c (µ, ν(u)) is continuous. Then by the intermediate value theorem for every k ∈ N there is some

u k > 0 with cost c (µ, ν(u k )) = r -ε k , forming a sequence (u k ) k∈N ⊆ [0, 1]. Then for every k there is a sequence (f jk ) j∈N ⊆ P (µ, ν(u k )) so that (f jk ) # µ * ν(k) and lim j∈N c d(id, f jk ) # µ = inf f ∈P (µ,ν(k)) c d(id, f k ) # µ = cost c (µ, ν(k)) = r -ε k .

Therefore for every k ∈ N there exists j k ≥ 0 so that for every j ≥ j k c d(id, f jk ) # µ ≤ r .

(3)

Let us pass directly to this subsequence of (f jk ) j∈N for every k ∈ N so that (3) holds for all j, k ∈ N.

Next by construction we have ν(u k ) → ν . Therefore (f jk ) j,k∈N has a subsequence in k so that (f jk ) # µ * ν . By ensuring (3) is satisfied, the sequences (f jk ) 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 ∈ L 0 (X, R ≥0 ) | g dµ ≤ r}. Then for f ∈ L 0 (Ω, 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) := B c (x, g(x)). Let L 0 (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) L 0 (X, Γ g ). Clearly for every a ∈ A µ (r) there is

r ≥ c(x, a(x)) dµ = c d(id, a) # µ, which shows {a # µ | a ∈ A µ (r)} ⊆ B c (µ, 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 L 0 (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 L 0 (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

(f i ) i∈N ⊆ A µ (r) so that lim i∈N f i dµ = f dν = sup ν∈Bc(µ,r) f dν,

proving equality in (6). 5) the perturbation size at each point is at most r (blue), however with the variable-radius risk the expected perturbation size is at most r.

PROOFS FOR §4: PROVABLE LIPSCHITZ REGULARISATION FOR KERNEL METHODS:

Theorem 3. Suppose k 0 , X 0 , and ν 0 satisfy Assumptions 1 and 2. Let {w s j : 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 ( PG ) ≤ L 2 , the Lipschitz condition λ max (G G) ≤ L 2 + ε is met with probability 1 -δ, as long as n ≥ Θ 1 ε 2 N 2 ε M 2 ε Q 2 ε log dNε δ , almost independent of d.

Here Θ hides all poly-log terms.

PROOFS AND MORE RESULTS FOR §4: KERNEL APPROXIMATION:



A.1 RANDOM SAMPLING REQUIRES EXPONENTIAL COST:

The most natural idea of leveraging the samples is to add the constraints g(w s ) ≤ L. For Gaussian kernel, we may sample from N (0, σ 2 I) 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 (x i ), y i ) + λ 2 f 2 H s.t. g(w s ) ≤ L, ∀s ∈ [n].

Unfortunately, this method may require O( 1 ε d ) samples to guarantee j g j 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) = 1 2 (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(w s ) ≤ 1 for samples w s that are drawn uniformly at random from B. This is equivalent to testing v ≤ 1, and to achieve it at least one w s 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 

v 2 = (W W ) -1/2 W 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.

A.2 SPECTRUM OF KERNELS:

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

𝜖𝜖 𝑣𝑣 𝑤𝑤 2 𝑤𝑤 1 = � 𝑤𝑤 1 � 𝑤𝑤 2

𝜖𝜖:

Figure 6 : Suppose we use a polynomial kernel with degree 2, and

f (x) = 1 2 (v x) 2 for x ∈ B. Then g(x) = (v x)v.

If we want to test whether sup x∈B g(x) 2 ≤ 1 by evaluating g(w) 2 on w that is randomly sampled from B such as w 1 and w 2 , we must sample within the ε balls around the intersection of B and the ray along v (both directions). See the blue shaded area. The problem, however, becomes trivial if we use the orthonormal basis { w1 , w2 }. in §4. Recall from Chapter 4 of Steinwart & Christmann (2008) that the integral operator for k and µ is defined by

T k = I k • S k : L 2 (X, µ) → L 2 (X, µ) where S k : L 2 (X; µ) → C(X), (S k f )(x) = k(x, y)f (y)dµ(y), f ∈ L 2 (X, µ), I k : C(X) → L 2 (X; µ), inclusion operator.

By the spectral theorem, if T k is compact, then there is an at most countable orthonormal set (ONS) {ẽ j } j∈J of L 2 (X, µ) and {λ j } j∈J with λ 1 ≥ λ 2 ≥ . . . > 0 such that

T f = j∈J λ j f, ẽj L2(X;µ) ẽj , f ∈ L 2 (X, µ).

In particular, we have ẽi , ẽj L2(X;µ) = δ ij (i.e., equals 1 if i = j, and 0 otherwise), and T ẽi = λ i ẽi .

Since ẽj is an equivalent class instead of a single function, we assign a set of continuous functions e j = λ -1 j S k ẽj ∈ C(X), which clearly satisfies e i , e j L2(X;µ) = δ ij , T e j = λ j e j .

We will call λ j and e j as eigenvalues and eigenfunctions respectively, and {e j } j∈J clearly forms an ONS. By Mercer's theorem,

EQUATION

and all functions in H can be represented by j∈J a j e j where {a j / λ j } ∈ 2 (J). The inner product in H is equivalent to j∈J a j e j , j∈J b j e j H = j∈J a j b j /λ j . Therefore it is easy to see that ϕ j := λ j e j , j ∈ J is an orthonormal basis of H, with Moreover, for all f ∈ H with f = j∈J a j e j , we have f, e j H = a j /λ j , f, ϕ j H = a j / λ j , and

f = j f, ϕ j H ϕ j = j λ j f, e j H ϕ j = j λ j f, e j H e j .

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 ).

A.3 GENERAL SAMPLE COMPLEXITY AND ASSUMPTIONS ON THE PRODUCT KERNEL:

In this section, we first consider kernels k 0 with scalar input, i.e., X 0 ⊆ R. Assume there is a measure µ 0 on X 0 . This will serve as the basis for the more general product kernels in the form of k(x, y) = d j=1 k 0 (x j , y j ) 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 X 0 , and the measure on X d 0 is trivially µ = µ d 0 . Theorem 4. Suppose {w s } n s=1 are drawn iid from µ 0 on X 0 , where µ 0 is the uniform distribution on [-v/2, v/2] for periodic kernels or periodized Gaussian kernels. Let Z := (k 0 (w 1 , •), k 0 (w 2 , •), . . . , k 0 (w n , •)), and

g 1 = 1 l l a=1 γ a g a 1 : X d 0 → R, where γ ∞ ≤ c 1 and g a 1 (y) = ∂ 0,1 k(x a , y) = h a 1 (y 1 ) d j=2 k 0 (x a j , y j ) with h a 1 (•) := ∂ 0,1 k 0 (x a 1 , •). Given ε ∈ (0, 1], let Φ m = (ϕ 1 , . . . ϕ m )

where m = N ε . Then with probability 1 -δ, the following holds when the sample size

EQUATION

where

(K 1 ) a,b = (h a 1 ) Z(Z Z) -1 Z h b 1 d j=2 k 0 (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 {w s } n s=1 . Applying the results in (6) for coordinates from 1 to d and using the union bound, we have that with sample size

EQUATION

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 ( PG ) + ε with probability 1 -δ, the sample size needs to be

n = 15 ε 2 c 2 1 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 L 2 in our training objective. In the case of Theorem 1, the values of Q ε , N ε , and M ε lead to a Õ( 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 v 2 . 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 {w s } n s=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 g a 1 , g b 1 H = h a 1 , h b 1 H0 d j=2 k 0 (x a j , x b j ) and k 0 (x a j , x b j ) ≤ 1, it suffices to show that for all a, b ∈ [l], h a 1 , h b 1 H0 -(h a 1 ) Z(Z Z) -1 Z h b 1 ≤ 3ε .

Towards this end, it is sufficient to show that for any h(

EQUATION

This is because, if so, then

h a 1 , h b 1 H0 -(h a 1 ) Z(Z Z) -1 Z h b 1 = 1 2 h a 1 + h b 1 2 H0 -h a 1 2 H0 -h b 1 2 H0 - 1 2 (h a 1 + h b 1 ) Z(Z Z) -1 Z (h a 1 + h b 1 ) -(h a 1 ) Z(Z Z) -1 Z h a 1 -(h b 1 ) Z(Z Z) -1 Z h b 1 ≤ 1 2 (4ε + ε + ε ) = 3ε .

The rest of the proof is devoted to (8). Since n ≥ m, the SVD of Λ -1/2 m Φ m Z can be written as U ΣV , where U U = U U = V V = I m (m-by-m identity matrix), and Σ = diag(σ 1 , . . . , σ m ). Define

EQUATION

+ Φ m Φ m h -h H0 . The last term Φ m Φ m h -h H0 is clearly below ε because by Assumption 1 and m = N ε Φ m Φ m h -h H0 ≤ |ϑ x | Φ m Φ m ∂ 0,1 k 0 (x, •) -∂ 0,1 k 0 (x, •) H0 + |ϑ y | Φ m Φ m ∂ 0,1 k 0 (y, •) -∂ 0,1 k 0 (y, •) H0 ≤(|ϑ x | + |ϑ y |)ε ≤ ε.

We will next bound the first two terms on the right-hand side of (9).

(i) By Assumption 1, k 0 (w s , •) -Φ m Φ m k 0 (w s , •) H0 ≤ ε, hence (Z -Φ m Φ m Z)α 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/2 j ϕ j , h H0 = e j , h H0 = e j , ϑ x ∂ 0,1 k 0 (x, •) + ϑ y ∂ 0,1 k 0 (y, •) H0 (10) ≤ sup x∈X e j , ∂ 0,1 k(x, •) H0 ≤ M ε .

As a result,

(Z -Φ m Φ m Z)α j H0 ≤ εn 1/2 • n -1/2 Λ -1/2 m Φ m h ≤ ε √ mM ε .

(ii) We first consider the concentration of the matrix R :=

1 n Λ -1/2 m Φ m ZZ Φ m Λ -1/2 m ∈ R m×m . Clearly, E {ws} [R ij ] = E {ws} 1 n n s=1

e i (w s )e j (w s ) = e i (x)e j (x) dµ(x) = δ ij .

By matrix Bernstein theorem (Tropp, 2015, Theorem 1.6 .2), we have Pr R -I m sp ≤ ε ≥ 1 -δ when n ≥ O(.). This is because (e 1 (x), . . . , e m (x))

2 ≤ mQ 2 ε , E {ws} [RR ] sp ≤ mQ 2 ε /n, and Pr R -I m sp ≤ ε ≥ 1 -2m exp -ε 2 mQ 2 ε n 1 + 2 3 ε ≥ 1 -2m exp -ε 2 5mQ 2 ε 3n ≥ 1 -δ,

where the last step is by the definition of n. Since R = 1 n U Σ 2 U , this means with probability 1 -δ,

EQUATION

Moreover, λ 1 ≤ 1 since k 0 (x, x) = 1. It then follows that 11), ( 10), and λ 1 ≤ 1).

(Φ m Φ m Z - √ nΦ m Λ 1/2 m U V )α H0 = Φ m Λ 1/2 m U ΣV 1 √ n V U Λ -1/2 m Φ m h - √ nΦ m Λ 1/2 m U V 1 √ n V U Λ -1/2 m Φ m h H0 = Λ 1/2 m U 1 √ n Σ -I m U Λ -1/2 m Φ m h 2 (because Φ m Φ m = I m ) ≤ λ 1 max i∈[m] 1 √ n σ i -1 Λ -1/2 m Φ m h 2 ≤ε √ mM ε (by (

Combining (i) and (ii), we arrive at the desired bound in (6).

Proof of Corollary 1. Since PG approximates G G only on the diagonal, PG -G G is a diagonal matrix which we denote as diag(δ 1 , . . . , δ d ). Let u ∈ R d be the leading eigenvector of PG . Then

λ max ( PG ) -λ max (G G) ≤ u PG u -u G Gu = u ( PG -G G)u = j δ j u 2 j (by (6)) ≤ 3c 1 1 + 2 N ε M ε ε.

The proof is completed by applying the union bound and rewriting the results.

A.4 CASE 1: CHECKING ASSUMPTIONS 1 AND 2 ON PERIODIC KERNELS Periodic kernels on X 0 := R are translation invariant, and can be written as k 0 (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):

EQUATION

)

We define µ 0 to be a uniform distribution on [-v 2 , 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 k 0 (x, y) cos(jω 0 y) dy = 1 v v/2 -v/2 κ(x -y) cos(jω 0 y) 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ω 0 x) cos(jω 0 z) + sin(jω 0 x) sin(jω 0 z)) dz (by periodicity)

= 1 v cos(jω 0 x) v/2 -v/2 κ(z) cos(jω 0 z) dz (by symmetry of κ) = √ 2π v F (jω 0 ) cos(jω 0 x).

And similarly,

1 v v/2 -v/2 k 0 (x, y) sin(jω 0 y) dy = √ 2π v F (jω 0 ) sin(jω 0 x).

Therefore the eigenfunctions of the integral operator T k are e 0 (x) = 1, e j (x) := √ 2 cos(jω 0 x), e -j (x) := √ 2 sin(jω 0 x) (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ω 0 e -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

EQUATION

where c 4 > 1 and c 6 > 0 are universal constants. Then Assumption 1 holds with

EQUATION

)

In addition, Assumption 2 holds with

Q ε = √ 2 and M ε = 2 √ 2π v n ε = √ 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 c 4 = 2 and c 6 = 1.6.

Proof of Theorem 1. First we show that h(x) := ∂ 0,1 k 0 (x 0 , x) is in H 0 for all x 0 ∈ X 0 . Since k 0 (x 0 , x) = j∈Z λ j e j (x 0 )e j (x), we derive h(x) = j∈Z λ j e j (x 0 )∂ 1 e j (x) = j∈Z λ j e j (x 0 )(-jω 0 e -j (x)) = ω 0 j∈Z λ j je -j (x 0 )e j (x). (15) h(x) is in H if the sequence λ j je -j (x 0 )/ λ j is square summable. This can be easily seen by ( 13):

ω -2 0 h 2 H0 = j λ j j 2 e 2 -j (x 0 ) = j∈Z λ j j 2 e 2 -j (x 0 ) = j∈Z λ j j 2 e 2 -j (x 0 ) = λ 0 + 2 j≥1 j 2 λ j ≤ 2c 4 c 5 c 4 -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 + j 2 )(1 + δ(j ≥ 1)) ≥ min(1, w -2 0 ) ε 2 2.1 . If j ∈ A, then -j ∈ A. Letting N 0 = {0, 1, 2, . . .}, we note j∈N0 1 1+j 2 ≤ 2.1. So h -Φ A Φ A h 2 H0 = w 2 0 j∈Z\A λ j j 2 e 2 -j (x 0 ) = w 2 0 j∈N0\A λ j j 2 (e 2 j (x) + e 2 -j (x))δ(j ≥ 1) + δ(j = 0) = w 2 0 j∈N0\A λ j j 2 (1 + δ(j ≥ 1)) = w 2 0 j∈N0\A λ j j 2 (1 + j 2 )(1 + δ(j ≥ 1)) 1 1 + j 2 ≤ ε 2 2.1 j∈N0 1 1 + j 2 = ε 2 2.1 j∈N0 1 1 + j 2 ≤ ε 2 . Similarly, we can bound k 0 (x 0 , •) -Φ A Φ A k 0 (x 0 , •) H0 by k 0 (x 0 , •) -Φ A Φ A k 0 (x 0 , •) 2 H0 = j∈Z\A λ j e 2 j (x 0 ) ≤ j∈Z\A λ j max(1, j 2 )e 2 j (x 0 ) = j∈N0\A λ α max(1, j 2 )[ e 2 j (x) + e 2 -j (x) δ(j ≥ 1) + δ(j = 0)] = j∈N0\A λ j max(1, j 2 )(1 + j 2 )(1 + δ(j ≥ 1)) 1 1 + j 2 ≤ 1 2.1 ε 2 j∈N0 1 1 + j 2 ≤ ε 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:

c 6 • c -|j| 4 < min(1, w -2 0 ) ε 2 2.1 ⇔ c -|j| 4 < 1 2.1 • c 6 min 1, 4π 2 v 2 ε 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, e j H ≤ ω 0 |je -j (x 0 )| ≤ 2π v √ 2 n ε = √ 2π v (N ε -1).

A.5 CASE 2: CHECKING ASSUMPTIONS 1 AND 2 ON GAUSSIAN KERNELS Gaussian kernels k(x, y) = exp(-x -y 2 /(2σ 2 )) are obviously product kernels with k 0 (x 1 , y 1 ) = κ(x 1 -y 1 ) = exp(-(x 1 -y 1 ) 2 /(2σ 2 )). It is also translation invariant. The spectrum of Gaussian kernel k 0 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 distribution N (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 c 0 = 1 2 (3 + √ 5) e j (x) = 5 1/8 2 j/2 exp - √ 5 -1 4 x 2 σ 2 1 √ j! H j 4 √ 1.25 x σ ,

where H j 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(-x 2 /(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 = 1 l l i=1 γ i k(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,1 k 0 (x a 1 , 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 c 4 = 1.25 and c 6 = 50. The expressions of N ε and M ε then follow from Theorem 1 directly.

A.6 CASE 3: CHECKING ASSUMPTIONS 1 AND 2 ON NON-PRODUCT KERNELS:

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., sup i→∞ sup x |e i (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 d log 2 1 ε +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 1 2 ∞ i=0 (-1 2 ) i x i . Using the standard multi-index notation with α = (α 1 , . . . , α d ) ∈ (N ∪ {0}) d , |α| = d i=1 α i , and x α = x α1 1 . . . 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 C k α x α y α = α 2 -|α|-1 C |α| α x α y α , where C k α = 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 in H are

EQUATION

)

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).

A.7 BACKGROUND OF EIGENVALUES OF A KERNEL:

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 -1 d where V d = π d/2 Γ ( d 2 + 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 -1 d y∈B α w 2 α 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 V d •(|α|+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 α+β ] α,β ,

A.8 BOUNDING THE EIGENVALUES:

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 Schneider (1953) ; Shi & Wang (1965) , we have

[1] ≥ r [2] ≥ . . .. Then by

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 Cauchy-Schwartz:

q 2 α+β = V -2 d y∈B y α+β dy 2 ≤ V -2 d y∈B y 2α dy • y∈B y 2β dy = q 2α q 2β .

To simplify notation, we consider without loss of generality that d is an even number, and denote 

r 2 α = w 2 α β w 2 β q 2 α+β ≤ w 2 α q 2α β w 2 β q 2β = w 2 α q 2α ∞ k=0 2 -k-1 β:|β|=k C k β q 2β ≤ w 2 α q 2α ∞ k=0 2 -k-1 k + d d max |β|=k C k β q 2β = w 2 α q 2α ∞ k=0 2 -k-1 k + d d max |β|=k k! d i=1 β i ! • 2 d i=1 Γ (β i + 1 2 ) V d • (2k + d) • Γ (k + d 2 ) = w 2 α q 2α V -1 d ∞ k=0 2 -k k + d d k! (2k + d)Γ (k + d 2 ) • max |β|=k d i=1 Γ (β i + 1 2 ) β i ! < w 2 α q 2α • b! π b d! • ∞ k=0 2 -k-1 (k + d)! (k + b)! (since Γ (β i + 1 2 ) < Γ (β i + 1) = β i !).

The summation over k can be bounded by

∞ k=0 2 -k-1 (k + d)! (k + b)! = 1 2 b! 2 d + d b ≤ 1 2 b!2 d + 2 b ≤ b!2 d ,

where the first equality used the identity 

r 2 α < w 2 α q 2α • b! π b d! b!2 d = 2 -l-1 l! d i=1 α i ! 2 d i=1 Γ α i + 1 2 V d • (2l + d) • Γ (l + b) (b!) 2 2 d π b d! ≤ 2 -l+d π -2b l!(b!) 3 d!(l + b -1)!(2l + d) (since Γ (α i + 1 2 ) < Γ (α 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 α. 

= - 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 ⇒ λ N L ≤ 2 -L .

This means that the eigenvalue λ i ≤ ε provided that i ≥ N L where L = log 2 1 ε . Since N L ≤ d L+1 , that means it suffices to choose i such that

i ≥ d log 2 1 ε +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) . 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).

A.10 EXTENSION TO ∞ -NORM ATTACKS FOR OUR KERNEL BASED METHOD:

We now extend our kernel based approach to ∞ norm ball attacks. Since most multiclass losses are 1-Lipschitz continuous with respect to ∞ norm on (f 1 (x), . . . , f 10 (x)), we will seek 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.

:



Figures and Tables:

Figure1: Empirical evaluation of the sum of the gaps from Theorems 1 and 2. The Lipschitz constants sup x∈X ∇f (x) q (left: p = 2, right: p = ∞, 1/p + 1/q = 1) were estimated by BFGS.

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.

Figure 3: Test accuracy under PGD attacks on the C&W approximation with 2 norm bound

Figure5: In the classical adversarial risk (5) the perturbation size at each point is at most r (blue), however with the variable-radius risk the expected perturbation size is at most r.

d s=1 (v ws ) 2 , where { ws } d s=1 is the orthonormal basis computed from the Gram-Schmidt process on d random samples {w s } d s=1 (n = d). With probability 1, n samples drawn uniformly from B must span R d as long as n ≥ d, i.e., rank(W ) = d where W = (w 1 , . . . , w n ). The Gram-Schmidt process can be effectively represented using a pseudo-inverse matrix (allowing n > d) as

the integer b := d/2. Now V d = π 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)

-k d + k b = 2 d . Letting l := |α|, we can continue by

Letting n l = l + d -1 l and N L = log λ N L ≤ N -1 L L l=0 n l -(l -b + 1) log 2 -2b log π -log l 2 < 2 log π as the coefficients of b)

Figure 7: Comparison of efficiency in enforcing Lipschitz constant by various methods

∞ ≤1 g 1 (x), . . . , g 10 (x) u ∞ ≤ L, where g c (x) := ∇f c (x).The left-hand side (LHS) can be bounded byLHS = sup x∈X max 1≤c≤10 g c (x) 1 ≤ max 1≤c≤10 sup ϕ H ≤1 G c ϕ 1 .Given the Nyström approximation Gc of G c , we can enforce the convex constraint of max

Figure11: 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.

Figure 12: Left: perturbed images at the end of 100-step PGD attack using (targeted) cross-entropy approximation. Right: classification on the perturbed image given by the trained Gauss-Lip. They are quite consistent with human's perception on the left images.

Figure13: Perturbed images at the end of 100-step PGD attack using (untargeted) C&W approximation. The 11 rows show the images after 0, 10, 20, ..., 100 steps of PGD.

Figure 14: (Another random trial) Left: perturbed images at the end of 100-step PGD attack using (targeted) cross-entropy approximation. Right: classification on the perturbed image given by the trained Gauss-Lip. They are quite consistent with human's perception on the left images.

Figure 15: (Another random trial) Perturbed images at the end of 100-step PGD attack using (untargeted) C&W approximation. The 11 rows show the images after 0, 10, 20, ..., 100 steps of PGD.

Peer review (iclr2020):
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. The authors then apply this idea to kernel methods and aim to minimise the true risk under a Lipschitz constant constraint. By using a bound between the Lipschitz constant and the largest eigenvalue of the Gram matrix, and by approximating the Gram matrix using Nystrom approximation, they express the Lipschitz constant constraint as a convex constraint. In the case of multiplicative kernels, this approach is shown to be efficient when only using a polynomial number of sample points for the Nystrom approximation. This method is tested on 3 standard datasets and is shown to outperform state-of-the-art deep learning based approaches over a wide range of perturbation in both l_2 and l_infinity norm. I have concerns about the novelty of the two main contributions of this work, which are Theorems 1 and 2: * Theorem 1 is a direct implication of Kantorovich duality, well known in optimal transport. * Theorem 2 is a rewriting of Proposition 3.1 in [1], which is not cited. My second concern is that the value of n (number of sample points used in Nystrom approximation) has not been specified for the presented experiments (unless I missed it). This is an important parameter to specify since the algorithm requires to invert an n by n matrix when applying Nystrom approximation, and scalability is usually the main criticism to kernel methods. When comparing the methods, would it be possible to estimate the Lipschitz constants of the models trained by the various algorithms ? This would enable to observe the true effect of the methods on the Lipschitz constant. Due to these two concerns above, I tend to reject this submission. A few comments about clarity: * It would help to provide the full kernel method algorithm (with more detail than Algorithm 1 in Appendix), e.g., explicit the convex constraint on gamma when training the SVM. * epsilon notation is used both for the accuracy in the Lipschitz constraint (as in Theorem 3) and the adversarial perturbation (as in equation 3 or Figures 4 and 5). In page 7, paragraph *Attacks*, delta is used once for denoting the perturbation, but epsilon is reused just after. More consistency in the notation would be greatly appreciated. * It is not clear in the main text that L is a bound on the gradient norms at sample points. It is made clearer in the appendix, so one sentence from the appendix should be added to the main text. Typos: * p.6 top: *and analogously for ...(x_1^a,.)* instead of *(x_1^b,.)* [1] Distributionally Robust Deep Learning as a Generalization of Adversarial Training, M. Staib, S. Jegelka, NIPS, 2017