Dense-Exponential Random Features: Sharp Positive Estimators of the Gaussian Kernel
Authors: Valerii Likhosherstov
Krzysztof Marcin Choromanski
Kumar Avinava Dubey
Frederick Liu
Tamas Sarlos
Adrian Weller
Venue: NeurIPS 2023 poster (2023)
Keywords: Gaussian kernel
softmax kernel
TL;DR: We propose new random features for Gaussian and softmax kernels which have lower variance and efficient computation algorithm.

Abstract: The problem of efficient approximation of a linear operator induced by the Gaussian or softmax kernel is often addressed using random features (RFs) which yield an unbiased approximation of the operator's result. Such operators emerge in important applications ranging from kernel methods to efficient Transformers. We propose parameterized, positive, non-trigonometric RFs which approximate Gaussian and softmax-kernels. In contrast to traditional RF approximations, parameters of these new methods can be optimized to reduce the variance of the approximation, and the optimum can be expressed in closed form. We show that our methods lead to variance reduction in practice (e^{10}-times smaller variance and beyond) and outperform previous methods in a kernel regression task. Using our proposed mechanism, we also present FAVOR#, a method for self-attention approximation in Transformers. We show that FAVOR# outperforms other random feature methods in speech modelling and natural language processing.

Dense-Exponential Random Features: Sharp Positive Estimators of the Gaussian Kernel
Valerii Likhosherstor"'l University of Cambridge
Krzysztof Choromanski"
Avinava Dubey" Google Research
Google DeepMind & Columbia University kchoro@google.com com
v. w.L.hosherstorrgmallcco com
Frederick Liu Google Research
Tamas Sarlos
Adrian Weller University Cambridge & The Alan Turing Institute
Google Research
Abstract
The problem of efficient approximation of linear operator induced by the Gaus- sian or softmax kernel is often addressed using random features features (RFs) which yield unbiased approximation of the operator's result. Such operators emerge in important applications ranging from kernel methods to efficient Transformers. We propose parameterized, positive, non-trigonometric RFs which approximate Gaus- sian and softmax-kerneh In contrast to traditional RF approximations, parameters of these new methods can be optimized reduce the variance of the approximation, and the optimum can be expressed in closed form. We show that our methods lead to variance reduction in practice ttmes smaller variance and beyond) and outperform previous methods in a kernel regression task. Using our proposed mechanism, we also present FAVOR#, a method for self-attention approximation Transformers. We show that FAVOR# outperforms other random feature methods in speech modelling and natural language processing.
Introduction
Random feature decomposition is an important technique for the lincarization of nonlinear kernel functions with theoretical guarantees such as unbiasedness and concentration around the true kernel value. Linearization allows significant reduction in computations from quadratic linear complexity in the size of the operator induced by the kernel. The technique emerged under the name of random 41 and was used in many applications such as kernel SVM
kitchen sinks (RKS) introduced in (39 40
43 29
38 3
dimensionality reduction [18 J neural networks [12] 22 49 24 [4], function-to-
function regression (34), kernel regression 28 2 nonparametric adaptive control 66, differentially-
private ML algorithms [. operator-valued kernels
32 g and semigroup kernels (11. An in-depth
theoretical analysis of random features was performed by 30
An exciting recent application of random features in the area of scalable Transformer networks a low-rank matrix when the sequence is long. However, the RKS family of methods relies on the Fourier transform, resulting in sin and cos types of random features, which were shown be unsuitable for application in Transformers due to negative values in the low-rank matrix. [16] proposed solution in the form of positive-valus random features relying on the exponential function (positive random features, PosRFs), yielding
[16 I3 17 26). where the self-attention matrix is approximated
Equal contribution
Work done while at University of Cambridge. The author is at Waymo now.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

method they called Fast Attention Via Orthogonal positive Random features (FAVOR+) for self- attention approximation. This solution was improved by [3] by means of careful choice of the linear combination parameters under the exponent, and the so-called homogeneity heuristic, which allows choice of one set of parameters for all approximated values. The resulting random features were called generalized exponential random features (GEREs), and the corresponding self-attention approximation method was termed +AVOR++.
Contributions: In this paper, we make a leap forward in the design of positive-valued random features by proposing dense exponential random features (DERFs) which contain both PosREs and GERFs as special cases. Instead of scalar parameters as in GERFs, DERFs rely on matrix parameters and dense quadratic forms inside the exponent. We show how to select parameters of the new random features efficiently without harming the overall subquadratic complexity.
More technically, our contributions
follows:
We show that the homogeneity heuristic of [31] may fact be viewed not as a heuristic, but closed-form optimum of the shifted log-variance objective.
2. We introduce DERF's and three special instantiations: DERFs (SDERFs), and simplified ADERFs (SADERFs). All these instantiations contain GERFs as special case (Figure [ left). For each instantiation we prove that closed-form optimum of the shifted log-variance objective can be found efficiently.
asymmetric DERF' (ADERFs), symmetric
3. We show that new variants result in lower variance than GERFs and other previous methods in practice (e.g. up to e0 times variance improvement as in Figure 1 right). Further, we show that DERFs outperform other random features in kernel regression and Transformer setups (speech modelling and natural language processing). We refer to the DERF-based self-attention approximation method as FAVOR#.
2 Prerequisites
2.1 Scaled softmax kernel and random features
By the scaled softmax kernel K(o), Rd x Rd + (0, +00). where a E R, we denote a mapping defined as K(a) y) = exp(a/|>]]| + y + ally|l²) for all x, Id where || is an Lorno Two important special cases of the scaled softmax kernel are I) the Gaussian kernel (x)-1/2((x))) = exp( -|x - / and 2) the softmax kernel K(0)) yx = exp(x y). For two sets of vectors X = (x(0) E 1L Ji=1 and y = (y()) E J==1
)=) by KC )) E RR we denote a matrix where (( (x,)),, = K(()) (x((() y()) for all Si
In this paper, we will be interested in the problem of computing (((0) (X' C where X' y and
matrix C E Axn are provided as an input. A naive solution requires ( (d + n)) computations for constructing (x,,) (O(T2d)) and computing the matrix multiplication K(aa (x,) )) x C
(O)L'n)). Instead, we will use an efficient Monte Carlo approximation of k(aa) (X,)) ) x C using the following notion of random features (RFs) for the scaled softmax kernel:
cfur30
miist/ifar10
DERF
23-
20
15
SDERF
GERF
SADERF
ADERF
30-
10
10"
Figure 1: (left) Venn diagram of the new types of random features (green) we propose. (right) Logarithm of the relative variance of different random feature maps on pairs of vectors sampled from CIFARIO and MNIST/CIFAR 10. A new random feature map SDERF results in a consistent variance reduction of the previous best method GERF, up to 1 and e times. Figure extends this plot.

---

Definition 2.1. By random features for the scaled softmax kemel k(n,, a E R, we denote
triple T = (v, f f(²) where v is a probability distribution over random objects w E Ω and f(r :P:2x 1 Fd + R. i E {1,2) 2}. are such mappings that, for all x, E IR
(x,) = Ev (x) (w.
(I)
The decomposition of type ① can be used for an efficient unbiased approximation of K(a) C. Let w(l)) w(dd) 92 be i.i.d. samples from v. Define matrices P, E RRLM where for all
1 S i & L.
Pip: M- ) - (w))) y())))
(2)
where P Sy: E Ie are column vectors corresponding the rows of P,S. Then according to 1. R = P is an unbiased Monte Carlo (MC) approximation of (X,,) on M samples. The variance Var Ri, proportional M, hence M is tradeoff parameter between computations and precision. Now, ŽC an unbiased approximation of ((a) J'CC but R PST a rank-M matrix, hence computing EC has (L) complexity. Assuming that sampling each (mm) and computing (-) are (d) operations, which is usually the case, precomputing P and S takes O(LM d) computations, resulting in total (LLM (d + n2) computational complexity. By choosing M < L, we obtain a significant reduction in computations compared the exact variant: O(LM (d + n)) < O (d + n)).
Var, (w,x(() f(2) (w, y()) of this approximation is inversely
- 12/2, emerge in kernel SVM [4). kernel regression 33 47 and in physics in the form of the Gauss transform [50]. Another important application is in the area of efficient Transformers and is discussed in the section [16]
Operations of type K(()) (X' C, especially for the Gaussian kernel
2.2 Random features for efficient Transformers
RFs found a prominent application in the area of efficient formers rely self-attention block for propagating information between elements of the sequence. the sequence length is and input matrices are denoted as Q, K, V E RLxd (queries, keys and values), then self-attention outputs the following matrix:
Transformers [66. Trans-
(3) where 1, IR. is vector of all ones, diag(-) returns diagonal (L x L)-sized matrix with the argument on the diagonal, X (x(*) d 1/4 d R. and y {y d d''/K,, Rd). Hence, substitution of I instead of (0()) (X,3) y))nn in reduces computational complexity from O(L'd) (()) (X, y) is the result of a softmax operation
Y = ( (X,))1L) ( x( (X.
e Rxd
to O(LMd) (n = d + 1). diag(K. (K(() performed on rows of d d-1/2QKT.
2.3 Existing random features for the softmax kernel
Representation is not unique and different RFs be proposed for single K(a), Note that if (v, f(), f()) are RFs for k(0), then (v, f(), f(2)) are RFs for K((a) for a E where flu (w, exp(a||x [2) (() (w, x). Hereafter we focus on the softmax kernel (())) without loss of generality.
15
proposed to use Irigonometric random features (TrigRFs) from 40 for efficient Transformers:
Derig = Rd+1 Virig Uni 2x]) x (O, La)
fe x).
(4)
<[exp[[[x 2/2) cos(W x 0).
(12) ((0 k).y) exp(((lyy1 2 cos(- y y00),
(5)
where w = (0, ()) Unif(.) denotes a uniform distribution on the argument set, (OO,,, La) is multivariate Gaussian distribution with mean 0₂ (vector of d zeros) and covariance matrix I (identity matrix of size d x d).
The next iteration of efficient attention approximators [6 observed problem with TrigRFs 5
The attention matrix ()ag(K(()) (()gg(()))))))))))) (x,)
(X,3) from 3 is right stochastic meaning that its
3

---

entries are nonnegative and each row sums to 1 due to the normalizing term diag (K(()) ( However, since fung can be arbitrary real numbers, P.S and, therefore, R can take negative values. Hence, diag(K1zz)-\ is not right stochastic in general and entries of Æ1, can take very small and/or negative values resulting in unstable behaviour when inverting dlag((K1,))> [16] therefore proposed a new type of positive random features (PosRFs) which have the form:
Opos = Rd Vpos =N(0,
(w, = Spoo [w,
(w, x) = cxp(w x ||x1
It is clear that such fpee only take strictly positive values resulting in the right stochastic diag ( and a stable Transformer training procedure.
31 extended PosREs, proposing generalized exponential random features (GERFs)) for K(0);
AGE RR, VGE AN 1)", fee ()) (2) ( D ex Bw C | /2) (6)
where A, B, C, D are real numbers? satisfying: - 8A > 0, B v1 -4, C -0.5, D
(1 - 444)4// 4/4 31 express B, C, through A and find a closed-form equation for the variance of
Varvar /Gg (w, (c)) x) (2)
(w, (w,y)= - (x00
(GE(A, d'log 1 - 4A1
2(1 - 4A)
+
1 1 8A ||x + y/12 ||x|| llyl1
VI 8A
The minimum variance corresponds to the minimum L(A, x, y) since (x,)(x,,))))) (x, does not depend
on A. Since x.), is defined for single pair of x. y and not for sets X' [31 propose
homogeneity heuristic when they replace ||x yyl ||x |yy in 7 with averages over X; j ( + y Σ, ||x and - E, [| respectively. This heuristic is based on the assumption that {x() and {y() are homogeneous and their statistics are tightly concentrated around the mean. After this, the minimum of (7) with respect to A can be found in closed form.
Σ.
3 Dense-exponential random features (DERFs)
We prove that the homogeneity heuristic (Section 2.3 corresponds a certain minimization problem. Then, we present DERFs which generalize GERFs and provide a tighter solution of that problem.
3.1 The objective minimized by GERFs
Our first contribution is showing that the homogeneity heuristic adopted in GERFs is actually an analytic solution of a certain optimization problem. Define
-(0;X,,,T)= = L-2
Σ
log(Var>[ (w,, (w, + K (x (0 y()))),
(8)
ISHJEL
where T (v, ((). p((2) (2)) are RFs for the kernel (0) and O are their parameters appearing implicitly
in f()), f(2). 8 is a mean log- variveri shifted by K (x(0), y())),, The best possible value of (8)
is log K (x()) y(i)) which corresponds to all variances Var (w) (w,x x()) (() (w,(())) being zero,
meaning that RFs provide exact kernel estimation. Hence, minimization of (8) leads to more precise
estimators on X. We call the loss function L(0; X, T) the shifted log-variance objective.
Since log concave loss L is conceptually similar to relative standard deviation: the smaller ((O) (x(+) (())), the higher Z is for the same amount of variance (l) ((2) (w, y()).
If TGE (VGE. JG:: JGE GE:JGE are taken in d8. then OGE = (A, D} and Z(OGE: X.Y, TGR) Ei, (GE(A;X' , yli). Using
<(1) ((2)
L 2
7
we get: Z(OGE; X.Y, TGE) =
dlog
1 4A
+ 8A
2 - 8A
Σ 1
(x + 20 I Σ 1p
L Σ 1p2
/I 8A
31 define these RF's for K (-1/2) but we adapt them for K (0) using the trick mentioned above.
B1 consider more generalized form when A, B, C, D are complex with an additional parameter = +1,
however only the subfamily (6) with s is proposed for use in the Transformer application.

---

That is, ((A;X')) )) coincides with (7) when ||x + yll,, | ; ||yl| are replaced by their average statistics computed on X, X. Hence, the homogeneity heuristic (Section 2.3 is nothing but minimiza- tion of (8). While in general it's unclear how to find a closed-form optimum of Varλ or Var(KC). the global minimum of (8) is feasible and can be computed in O(1) time. Further, in the case of GERF, the values minimizing (8) lead to very good results in large-scale applications efficient Transformers shown by [31 without knowing about (8) objective. In the next section present extensions of GERFs which aim to minimize (8) in closed form.
Towards DERFs
Dense-exponential random features (DERFs) are an extension of GERFs where scalars A, B, C are replaced with dense matrices. DERFs may be viewed as a generalization that contain the previously introduced classes as special cases. We define DERFs as follows: ΩDE = R VDF = (0, 1)d and for k (1,
A ( Dexp([ Aw + w ( x x x)
where (k)) C(k) Rdxd DE R. A S (a set of x d real symmetrie matrices). Clearly, GERFs with parameters A, B, C, D can be expressed via DERFs with parameters A Ala. B B
BLa C()) C CLa, D is unchanged. Our first theoretical result is giving the conditions when
TDE= = (VDE, JURY ( are valid RFs:
Theorem 4.1. Lel the following conditions hold: 8A La (B()
(I. - La
C(K) = (B ( - D det(La 4A) where [1,2]. Then TDE are RFs
(2) for (()) and for all x, E Rd;; Var PD (w, x yw
detetLLa - 8A) exp 2x (cO) (la
+2yT (C(2)+((())) (L, (L. -IA y+4x ((, -B)))) - (9)
Our ultimate goal is to find optimal parameters A, (k), ((k) and D minimizing the variance of the low-rank approximation of ((() (X')) where sets X, are provided. Our first observation is that we can assume that A Da (a set of d x real diagonal matrices). Indeed, any symmetric A can be expressed as QAQ T where Q E Og (a set of orthogonal matrices {Z E dxd Z = La}) and
Á E Da. Let w ~ NOO,,, Ia). Then, for any x E R k E (1,2),
po (w Dexp(e w QAQ w+ (x))x x TC())
D Dexpl Au w (1)) x x ((0)) = T(k) (w,, x),
where ((k) = Q (() ت Q N'OO,, Ia) since the distribution w ~ N (02z, La) is isometric, ie. colation-invariant and T DERFs with parameters A. (k), ((k), D. We conclude that with any A, ,(k) JDE(W, (w,x) can be expressed as DERFs ((8) with À E D.. Hence, hereafter we only consider
A E ID, without loss of generality.
Since ((), (c)) are dense matrices in general, evaluation of plk) (w. x) takes ((d)) time which is bigger than the (dd) complexity for TrigRFs, PosRFs and GERFs. However, P and S matrices can be still computed in a time subquadratic in L. For that, precompute ( (0)
()) for all A E {1, 2}. 1 < m S M, < i, < L in O(MM + L)d²) time. Then, computing (l (m)) ))), /!2) ( (m)) y()) for all 1 < i,j < L, < m M takes ((LMd) operations. The complexity of constructing (2) then is (LL(M d?) Md?) which is still subquadratic in L.
Our goal is to minimize L(ODE; X, TDE) for θDE = (A, (()) 2 c(1), (( D). However, we find that even for a single pair of x, y it's unclear how to minimize the variance ( in closed form. Hence, below consider special cases where analytic solution is feasible.
5

---

norma!
Schere
hele Togen
05
00
ADC1!
153
-10
-10
1314
107
12
33"
melt
ofbr3
misoocoor
22
10
Figure 2: Log of the relative variance of new and existing RF mechanisms, mean value over multiple samples. 0.1 S a S 1.
4.1 Asymmetric dense-exponential random features
Define RFs TADE = (VADE. SADE. SADE) in the same way TDE with the only difference that A ALa where A E R. We refer to these RFs as asymmeetric dense-exponential RFs (ADERFs) since SADE * (2) in general. The only additional restriction of ADERFs com- pared to DERFs is that all diagonal entries of A ID are the same. The parameters of TADE are ADE (()) (2), C()) C(2), D}. By OADE denote set of all possible O ADE resulting in correct RFs for the kernel k(0), i.e. satisfying Theorem 4.1 with A Ald. The following result gives analytic formula for minimum of L(OADE; X, TADE.. In the theorem, we use notions of SVD and eigendecomposition of symmetrie matrix (5) (all proofs are in the Appendix).
= Theorem 4.2. Let X = fx(() E R y = (yli) E R Let ( = E (x))),
( E (( Suppose that M(), M(2)) E Sa are nonsingular (implying that
2 d). Define y(() = d- (EL, (E)
E R. For E {1,2), let (
Q ((x) A (k) (Q(K)) be eigendecomposition of symmetric ((k) where Q E Os. ( Da has strictly positive diagonal values since (R)) > by definition and (1k) is nonsingular Let UEVT
be SVD decomposition of (A ( Q (2) (A (2)) where U, V 0d, E Da has nonnegative diagonal entries.
One of the solutions OADE A(O N, KKY,TAA TADE) is as follows. Set &b = 2d Σ Σn + 2 (3) and, for E (1,2).
(A,BO) (2) CO), C of minoADRE@AAR
A
[ 26
12
( v1 4AΣ
(1) ()))
16
B
VI - 1AΣ 1/ (A
c(.)_ 1
- 2(1-4A) (+) T (()) D -4)
Further, we have: I(0) X. Y, TADE)
d( log(1 44)
log(1 (8A) 2(1 8A)
d
Eus + (()) tal))) 22x
(10)
2
Theorem 4.2 implies an algorithm for finding 0% OADE efficiently. Namely, compute (() E {1,2) (O(Ld") time) and y (O(Ld) time). Then, perform matrix decompositions to obtain (()) A (1) k E [1, 2). and U, E, in O(d9) time. After that. A, (1)), B22), c(1) c(2) D can be all evaluated
in (dd) time using formulae from Theorem 4.2
The total time complexity of the approximation
scheme is therefore O(L(Md d²) M d
which is subquadratic in L as required.

---

4.2 Symmetric dense-exponential random features
Define TSDE (USDE, JSDE J SDE/ in the same way as TDE with the only difference that B
(2)
( B. From the conditions Theorem 4.1 it follows immediately that also C() C(2) C.
Hence, JADEE = SDD and we refer to these RFs as symmetric dense-exponential RF's (SDERFs).
The parameters of 7SDE are OSDE {A, B, C, D3. By ÔSDE denote set of all possible OSDE'S resulting in correct RFs for the kernel K(O), i.e. satisfying conditions from Theorem 4.1 with
(()) B. C C. k E {1.2]. The following theorem gives an analytic solution for global minimum of Z(OsDEi X TSDE).
and let M ( be defined as in Theorem 4.2 (but without the restriction of nonsingularity) and define μ²() Ei-1 x(9 R. ( A (3) (Q))) be eigendecomposition of a symmetric u - ビ E IRS. Further, let Q
Theorem 4.3. Let X' (x(*) HR y' (y() 1/
positive semidefinite matrix M(1) + (x(4) (1(()) + 4 ( + M where Q e Oa and
(3) E Da with nonnegative diagonal entries. Let the entries on the diagonal of A (3) be sorted in the non-ascending order.
One of the solutions OSDE = [A, D) of minnnnneeeess L(OSDEE X, TSDE) is as follows.
(3) 24 V(2 (9) + +1) 2 + 8A 1.1
A E Da. for all 1 S S d: A,i = 16
B = (Ia
4A) 4A)1/2(0 33) T,, C.C --==--- It D = det (Ia 4A/1/4 Further, we have: Z(OSDE; &, (TDE)) =
L
Σ log(1- -A,1)) 2 log(1 SA..) (1+ A
1 11 I - 11 (11)
Again, Theorem 4.3 implies algorithm for finding OSDE in time subquadratic in L. That is, compute M
1(()) (
14(4, (ن) in (Ld)) total time. Then, perform an eigendecompositiin (3) in O(d3) time. After that, A, B, C, can be computed in Odd) time using
to obtain Q (3) A
formulae from Theorem 4.3) The total time complexity of the approximation scheme is the same as
for ADERFs: (L + d²) + Md + dd) or (L((Md+ + d²) + Md2) if we assume that L 2 d.
4.3 Simplified ADERFs
While having compact and closed-form expression, both ADERFs and SDERFs rely on eigende- composition and SVD decompositions: operations for which implementation has not yet matured in popular deep learning libraries with GPU and TPU support. For this reason, we propose simplified ADERFs (SADERFs) TSADE = (DSADE, SSADE
(()) {(2) JSADe) which extend GERFs but require only basic
unary operations. SADERFS are defined via GERFs as follows: QSADE Rdd USADE = (o, 1)d,
((1) SSADE (w, x) = fGg(W, (w, Tx). SSADE (w = f (w. w "yyy) where T E Da is diagonal ma- trix with nonzero diagonal entries. First of all, TSADE are valid random features for the soft- (1) SSADE x) S)SDA y)] EPCE U (w, Vx) (w, v "y]
((2)
max kernel ((0) since EPSADE
K K(O) [xx y- = ((O)(x, y), where we use (1)0) (x, = exp(x y) by the definition.
We find & by optimizing the objective Z(OGE: X, Y, TGE):
( for TSAdE, the form of which is easily deduced from
Z(OSADE: X,Y, TSADE) dlog
1A
1
8A
1
Σ
+
v1 8A
1
8A
1,3
D Vx "P 20
1
Vx
y
(1)
12
2
(xx()) x0))
(12)
(1A)))
1
1.)
where we move term not depending on L to the lend- side. Since 1 - 8.A > 0, we conclude that minimizing 12 is equivalent to minimizing Ens || [xx(0) + 220011==
Σ>
Vr.x " Σ> (x 3/)) + ())))))
(13)
7

---

82A 1TA0O0N, -176
barknote,
Tight
12
50850
10
16,,00
113,00
1463
12400
Accumcyy 3.
80
30-
19650
100
Chenk 21055
wrerage
62
163,00
103,00
Figure Kernel classification, accuracy (%). The last plot shows average curves over benchmarks. We observe that our proposed method SDERF, shows the best accuracy across most of the (benchmark, M) pairs and also shows the best average performance.
reduces to independent optimization problems with respect L 1 S 1 d. Each problem is convex and the solution is found trivially by setting the derivative to zero: for all
Optimizing 13
1 1 d: Vi (x) )) Σ, (x )))1/4 The solution be computed in O(dL) time, after which the parameters of [GE S be found efficiently described in Section 2.3
easy to see that TSADE are a special case ADERFs (Section 4.1 which explains their name. Furthermore, the case V Ia reduces TSADE to TGE. hence the latter special case of the former. Figure (left) illustrates all the new types of random features a Venn diagram.
Experiments
In this section, we evaluate DERFs experimentally in various machine learning applications. More details about each experiment can be found in Appendix
5.1 Variance comparison
We follow the variance comparison setup from 31: we sample pairs of vectors x,y and compute relative variances of the approximation Vark (x, y) K (x where R(0) denotes the RF approximation and Vark (x, y) is evaluated via (9) We set d 64 as in [31] and take different regimes for sampling x,y: normal where x, are drawn from (O ( sphere where x, y are drawn uniformly on a sphere add-1 heterogen where x, are drawn from different distributions (O,,,, I² and ((al4,0"L)) mnist and cifar10 are where x, y,y are random images from MNIST 19 or CIFARIO (27, resized to x 8, scaled by a > and flattened. Finally, mnist/cifar10 is regime where is drawn in mnist and y is drawn in cifar10.
We do not report SADERFs since they're re special case of ADERF' (Figure 2). SDERFs outperform or are par with other methods in all setups about e times better than GERFs in heterogen, mnist and mmist/cifart and about e10 times better in cifar10. Further, ADERFs outperform GERFs by around e times in mnist/ /ciffrrfl where and y drawn "asymmetricalll".

---

TRANS-8
TRANS-32
0009
0.00
0.04
0200
0.04
240
0325
0279
FANce
RACER
0004
THANSS-22
NS7-8
0.60
0002
0
200
000
0010
Figure 4: Comparison of FAVOR# using SDRF with FAVOR+ Performer for regular Conformer- Transducer training with m random features (TRANS-m) as well as the Noisy Student Training variant with random features (NST-m) on the LibriSpeech corpus. We report commonly used normalized word error rate (NWER) metric.
5.2 Kernel classification
In this experiment, we compare accuracy of different RF methods in kernel classification on benchmarks from UCI (21). following the setup of (31). Kernel regression [3 47] is applied for predicting class probabilities. Training objects denoted u ..... ( E Rdd and their one-hot labels as r(1), (1),,..,,,,,,,, (L) E Rd. During testing, the goal is to predict the class of new object u* as argmax g where r = El= ((-055)ou"" (u(()))r and > is tuned the validation set. With (n)LM) preprocessing, RFs used find unbiased approximation of r' in ((nM) instead of (nL) exact computation. For each benchmark, we range the values of M from 2 2 (Figure 3). We observe that SDERF, which is proposed in this paper, shows the best accuracy across most of the (benchmark, M) pairs and also shows the best average performance.
5.3 DERFs for long-sequence Transformers
In this section, we evaluate DERFs self-attention approximation setups [16]. We refer the DERF-basee self-attention approximation method as FAVOR#
training
5.3.1 Speech modelling
In our first set of experiments, focus speech models. We train Performer-encoderrs and test them on the LibriSpeech corpus [35], commonly used for benchmarking speech models. considered two Transformers architectures/rraining setups: (a) cenffomme--Transduee [33] trained in regular way (TRANS) as well as: (b) the Noisy Student Training (NST) variant introduced in [36). We compare "performized" variants of these architectures, applying FAVOR# with SDERF (since
worked best in the previous setups) as well as FAVOR+ [31].
In the first setting, we see that FAVOR# consistently outperforms FAVOR+ for smaller m (where reduced variance of the softmax-kernd estimation more critical) and both achieve similar scores for larger m. In the NST-experiment, we focused on the smaller m variant, where FAVOR# again beats FAVOR++ All results are presented Fig.
5.3.2 Natural language processing
The General Language Understanding Evaluation (GLUE) benchmark [46] consists different natural language understanding tasks with the sequence length ranging from 32 to 128. We use this to test performance of different low rank attention methods on NL.P tasks. We used the same training parameters mentioned in (20 (see Appendix B.4.2 for details). We warm start all low-rank

---

Table I: GLUE Dev results base sized models. Number of training examples is reported below each task. MCC score is reported for CoLA, FI score is reported for MRPC, Spearman correlation is reported for STS-B, B, and accuracy scores reported for the other tasks. The best result, second best.
System
MNL[(m)
QQP
QNLI
SST-2
CoLA
STS-B
MRPC
RTE
392k
363k
108k
67k
8.5k
5.7k
3.5k
2.5k
Local Attention
37
75.68
86.52
87.31
89.71
61.07
83.65
75.47
59.84
Global + Sliding Wndw
33
81.82
90.35
89.91
92.43
61.81
85.58
87.80
67.15
ELU 26
82.58
90.05
89.81
92.43
58.63
87.91
87.50
67.15
RELU 16
82.49
90.71
89.68
92.32
57.57
88.15
87.25
68,95
FAVOR+ 16
77.69
86.69
89.41
91.80
54.87
83.78
80.73
66.
FAVOR++ 31
82.29
90.43
89.73
92.20
58.85
85.90
88.73
67.63
FAVOR#
82.69
90.68
90.01
92.53
59.33
85.48
87.99
69.68
Transformers with a pre-trained BERT-base model checkpoint 20 thus contrasting how well the low rank methods approximate the softmax kernel.
We compared FAVOR+ [31] FAVOR+ [16], ELU 266 ReLU [6] variants of the Performers [16 and sparse attention 37 against the FAVOR# variant and report the results in Table [ We couldn't use SDERF F in this setup because eigendecompossiton led errors on TPUs due different implementation compared to the speech modelling experiment. For this reason, we used SADERF which doesn't 't require any matrix decompositions. On most tasks we find that FAVOR# is the best or second best performing variant showcasing its effectiveness in modelling the softmax kernel for Transformers.
6 Conclusion
We proposed an extension of generalized exponential random features (GERFs) for the Gaussian and softmax kernels: dense-exponential random features (DERFs). DERFs employ matrix parameters and are more flexible than GERFs. We evaluated DERFs in several applications such as kernel regression and two Transformers training setups, demonstrating downstream performance benefits.
Limitations & broader impact. Optimizing matrix parameters the most general formulation of DERFs (Theorem 4.1l could lead to further variance reductions. It remains open question how to find closed form solution of the optimum which we leave to future work. Our work is primarily theoretical but has broad applications. A prominent application is efficient Transformers which should be used responsibly due to their potential for misuse and significant societal impact, and carbon footprint [48 10).
7 Acknowledgements
L. acknowledges support from the Cambridge Trust and DeepMind. V. was part- employed by Google while a PhD student. A.W. acknowledges support from a Turing Al Fellowship under EPSRC grant EP/V025279/1, The Alan Turing Institute, and the Leverhulme Trust via CFI.
References
[1] Nir Ailon and Edo Liberty. An almost optimal unrestricted fast jownson-lindenstennnne transform. ACM Trans. Algorithms, 93):21:1-21:12, 2013.
[2] Haim Avron, Michael Kapralov, Cameron Musco, Christopher Musco, Ameya Velingker, and Amir Zandieh. Random Fourier features for kernel ridge regression: Approximation bounds and statistical guarantees. In Doina Precup and Yee Whye Teh, editors, Proceedings of the 34th International Conference on Machine Learning, ICML 2017, Sydney, NSW, Australia, 6-11 August 2017, volume 70 of Proceedings of Machine Learning Research, pages 253-262. PMLR, 2017.
[3] Haim Avron, Vikas Sindhwani, Jiyan Yang, and Michael W. Mahoney, Quasi-monte carlo feature maps for shift-invariant kernels. J. Mach. Learn. Res.,
2016.
10

---

[4] Iz Beltagy, Matthew E Peters, and Arman Cohan. Longformer: The long document transformer. arXiv preprint arXiv:2004.05/50, 2020.
[5] Emily M. Bender, Timnit Gebru, Angelina McMillan- McMillan-Major, and Shmargaret Shmitchell, On the dangers of stochastic parrots: Can language models be too big? In Proceedings of the 2021 ACM Conference on Fairness, Accountability, and Transparency, FAccT "21, page 610-623, New York, NY, USA, 2021. Association for Computing Machinery.
[6] Nicholas M. Boffi, Stephen Tu, and Jean-Jacques E. Slotine. Nonparametric adaptive control and prediction: Theory and randomized algorithms. In 60th IEEE Conference on Decision and Control, CDC 2021, Austin, TX, USA, December 14-17, 2021, pages 2935-2942. IEEE, 2021.
[7] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman- Milne,, and Qiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018.
[8] Romain Brault, Markus Heinonen, and Florence Buc. Random fourier features for operator- valued kernels. In Robert Durrant and Kee-Fung Kim, editors, Proceedings The 8th Asian Conference on Machine Learning, volume 63 of Proceedings of Machine Learning Research, pages 110-125, The University of Waikato, Hamilton, New Zealand, 16-18 Nov 2016. PMLR.
[9] R. Brockett. Dynamical systems that sort lists, diagonalize matrices, and solve linear programming problems. Linear Algebra and its Applications, 146:79-91, 1991.
[10] Nicholas Carlini, Florian Tramèr, Eric Wallace, Matthew Jagielski, Ariel Herbert- Voss, Kather- ine Lee, Adam Roberts, Tom B. Brown, Dawn Song, Úlfar Erlingsson, Alina Oprea, and Colin Raffel. Extracting training data from large language models. CoRR, abs/2012.07805, 2020.
[11] Kamalika Chaudhuri, Claire Monteleoni, and Anand Sarwate. Differentially private empirical risk minimization. J. Mach. Learn. Res., 12:1069- 109, 2011.
[12] Youngmin Cho and Lawrence K. Saul. Kernel methods for deep learning. In Yoshua Bengio, Dale Schuurmans, John D. Lafferty, Christopher K. I. Williams, and Aron Culotta, editors, Advances in Neural Information Processing Systems 22: 23rd Anmual Conference on Neural Information Processing Systems 2009. Proceedings of meeting held 7-10 December 2009, Vancouver, British Columbia, Canada, pages 342 350. Curran Associates, Inc., 2009.
[13] Krzysztof Choromanski, Haoxian Chen, Han Lin, Yuanzhe Ma, Arijit Sehanobish, Deepali Jain, Michael S. Ryoo, Jake Varley, Andy Zeng, Valerii Likhosherstov, Dmitry Kalashnikov, Vikas Sindhwani, and Adrian Weller. Hybrid random features. In International Conference on Learning Representations (ICLR), 2022.
[14] Krzysztof Choromanski, Carlton Downey, and Byron Boots. Initialization matters: Orthogonal predictive state recurrent neural networks. In 6th International Conference on Learning Rep- resentations, ICLR 2018, Vancouver, BC, Canada, April 30 May 3, 2018, Conference Track Proceedings. OpenReview.net, 2018.
(15) Krzysztof Choromanski, Valerii Likhosherstov, David Dohan, Xingyou Song, Jared Davis, Tamás Sarlós, David Belanger, Lucy J. Colwell, and Adrian Weller. Masked language modeling for proteins via linearly scalable long-context transformers. CoRR, abs/2006.03555. 2020.
[16] Krzysztof Marcin Choromanski, Valerii Likhosherstov, David Dohan, Xingyou Song, Andreea Gane, Tamas Sarlos, Peter Hawkins, Jared Quincy Davis, Afroz Mohiuddin, Lukasz Kaiser, David Benjamin Belanger, Lucy J Colwell, and Adrian Weller. Rethinking attention with performers. In International Conference on Learning Representattinns 2021.
[17] Sankalan Pal Chowdhury, Adamos Solomou, Avinava Dubey, and Mrinmaya Sachan. On learning the transformer kernel. Transactions of Machine Learning Research, 2022.
[18] Sanjoy Dasgupta and Anupam Gupta. An elementary proof of theorem of johnson and lindenstrauss. Random Struct. Algorithms, 22(1):60-65, 2003.
[19] Li Deng. The mnist database of handwritten digit images for machine learning research. IEEE Signal Processing Magazine, 1916)114 2012.

---

[20] Jacob Devlin, Ming- deep bidirectional transformers for language understanding. arXiv preprint arXiv: 2018.
Ling-Wei Wei Chang, Kenton Lee, and Kristina Toutanova. Bert: Pre-training of
[21] Dheeru Dua and Casey Graff. UCI machine learning repository, 2017.
[22] Lukas Gonon. Random feature neural networks learn Black-Schole type PDEs without curse of dimensionality. CoRR, 0bbs/2106.08000, 2021.
[23] Anmol Gulati, James Qin, Chung-Chene Chiu, Niki Parmar, Yu Zhang, Jiahui Yu, Wei Han, Shibo Wang, Zhengdong Zhang. Yonghui Wu, and Ruoming Pang. Conformer: Convolution- augmented transformer for speech recognition. CoRR, abs/2005.00 18000 2020.
[24] Insu Han, Haim Avron, Neta Shoham, Chaewon Kim, and Jinwoo Shin. Random features for the neural tangent kernel. CoRR. abs/2104.01351, 2021.
[25] Charles R. Harris, Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathamiel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fer- nández, del Río, Mark Wiehe, Pearu Peterson, Pierre Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant. Array programming with NumPy. Nature, 585(7825):357-362, September 2020.
Kevin Sheppard, Tyler
[26] Angelos Katharopoulos, Apoorv Vyas, Nikolaos Pappas, and François Fleuret. Transformers are RNNs: Fast autoregressive transformers with linear attention. In Proceedings of the 37th International Conference on Machine Learning. 2020.
[27] Alex Krizhevsky, Vinod Nair, and Geoffrey Hinton. Cifar-10 (canadian institute for advanced research).
[28] Valero Laparra, Diego Marcos Gonzalez, Devis Tuia, and Gustau Camps-Valts Large-scar scale random features for kernel regression. 2015 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), pages 17-20, 2015.
[29] Quoc V. Le, Tamás Sarlós, and Alexander J. Smola. Fastfood computing hilbert space expansions in loglinear time. In Proceedings of the 30th International Conference on Machine Learning, ICML 2013, Atlanta, GA, USA, 16-21 June 2013, volume 28 of JMLR Workshop and Conference Proceedings, pages 244-252. JMLR.org. 2013.
[30] Zhu Li, Jean-Francois Ton, Dino Oglic, and Dino Sejdinovic. Towards unified analysis of random Fourier features. J. Mach. Learn. Res., 12:108: 1-00:51, 2021.
[31] Valerii Likhosherstov, Krzysztof Choromanski, Avinava Dubey, Frederick Liu, Tamas Sarlos, and Adrian Weller. Chefs' random tables: Non-trigonometrit random features. In Advances in Neural Information Prucessing Systems. Curran Associates, Inc., 2022.
[32] Ha Quang Minh. Operator-valued Bochner theorem, Fourier feature maps for operator-valued kernels, and vector-valued learning. CoRR, abs/1 (668006639, 2016
[33] E. A. Nadaraya. On estimating regression. Theory of Probability & Its Applications, 9(1::141 142, 1964.
[34] Junier B. Oliva, Willie Neiswanger, Barnabás Póczos, Eric P. Xing, Hy Trac, Shirley Ho. and Jeff G. Schneider. Fast function to function regression. In Guy Lebanon and S. V. N. Vishwanathan, editors, Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics, AISTATS 2015, San Diego California, USA, May 9-12, 2015, volume 38 of JMLR Workshop and Conference Proceedings. JMLR.org, 2015.
[35] Vassil Panayotov, Guoguo Chen, Daniel Povey, and Sanjeev Khudanpur. Librispeech: An ASR corpus based on public domain audio books. In 2015 IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP 2015, South Brishane, Queensland, Australia, April 19-24, 2015, pages 5206-5210, IEEE, 2015.
12

---

[36] Daniel S. Park, Yu Zhang, Ye Jia, Wei Han, Chung Cheng Chiu, Bo Li, Yonghui Wu, and Quoc V. Le. Improved noisy student training for automatic speech recognition. In Helen Meng, Bo Xu, and Thomas Fang Zheng, editors, Interspeech 2020, 21 Annual Conference of the International Speech Communication Association, Virtual Event, Shanghai, China, 25-29 October 2020, pages 2817-2821. ISCA, 2020.
[37] Niki J. Parmar, Ashish Vaswani, Jakob Uszkoreit, Lukasz Kaiser, Noam Shazeer, Alexander Ku, and Dustin Tran. Image transformer. In International Conference on Machine Learning (ICML), 2018.
[38] Jeffrey Pennington, Felix X. Yu, and Sanjiv Kumar. Spherical random features for polynomial kernels. In Corinna Cortes, Neil D. Lawrence, Daniel D. Lee, Masashi Sugiyama, and Roman Garnett, editors, Advances in Neural Information Processing Systems 28: Annual Conference on Neural Information Processing Systems 2015, December 7-12, 2015, Montreal, Quebec, Canada, pages 1846-1854, 2015.
[39] A. Rahimi and B. Recht. Uniform approximation of functions with random bases. In 2008 46th Annual Allerton Conference on Communication, Control, and Computing, Los Alamitos, CA, USA, sep 2008. TEEE Computer Society.
[40] Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In J. Platt, Y. Singer, and S. Roweis, editors, Advances in Neural Information Processing Systems, volume 20. Curran Associates, Inc., 2007.
D. Koller,
[41] Ali Rahimi and Benjamin Recht. Weighted sums of random kitchen sinks: Replacing minimiza- tion with randomization in learning. In Daphne Koller, Dale Schuurmans, Yoshua Bengio, and Léon Bottou, editors, Advances in Neural Information Processing Systems 21, Proceedings of the Twenty-Second Annual Conference on Neural Information Processing Systems, Vancouver, British Columbia, Canada, December 8-11, 2008, pages 1313-1320. Curran Associates, Inc., 2008.
[42] Bharath K. Sriperumbudur and Zoltán Szabó. Optimal rates for random Fourier features. In Corinna Cortes, Neil D. Lawrence, Daniel D. Lee, Masashi Sugiyama, and Roman Garnett, editors, Advances in Neural Information Processing Systems 28: Annual Conference on Neural Information Processing Systems 2015, December 7-12, 2015, Montreal, Quebec, Canada, pages 1144-1152, 2015.
[43] Yitong Sun, Anna C. Gilbert, and Ambuj Tewari. But how does it work in theory? Linear SVM with random features. In Samy Bengio, Hanna M. Wallach, Hugo Larochelle, Kristen Grauman, Nicolò Cesa-Bianchi, and Roman Garnett, editors, Advances in Neural Information Processing Systems 31: Annual Conference Neural Information Processing Systems 2018, NeurlPS 2018, December 3-8, 2018, Montréal, Canada, pages 3392- 2018
[44] Danica J. Sutherland and Jeff G. Schneider. On the error of random Fourier features. In Marina Meila and Tom Heskes, editors, Proceedings of the Thirty-First Conference on Uncertainty in Artificial Intelligence, UAI 2015, July 12-16, 2015, Amsterdam, The Netherlands, pages 862 871, AUAI Press, 2015.
[45] Lloyd N. Trefethen and David Bau. Numerical Linear Algebra. SIAM, 1997.
[46] Alex Wang, Amanpreet Singh, Julian Michael, Felix Hill, Omer Levy, and Samuel R Bowman GLUE: A multi-task benchmark and analysis platform for natural language understanding. arXiv preprint
2018.
[47] Geoffrey S. Watson. Smooth regression analysis. Sankhya: The Indian Journal of Statistics, Series A (1961-2002), 26(4):359-372, 1964.
[48] Laura Weidinger, John Mellor, Maribeth Rauh, Conor Griffin, Jonathan Uesato, Po-Sen Huang. Myra Cheng. Mia Glaese, Borja Balle, Atoosa Kasirzadeh, et al. Ethical and social risks of harm from language models. arXiv preprint arXiv:21 2021.
[49] Jiaxuan Xie, Fanghui Liu, Kaijie Wang, and Xiaolin Huang. Deep kernel learning via random Fourier features. CoRR, abs/1910.02660, 2019.
13

---

[50] Yang, Duraiswami, Gumerov, and Davis. Improved fast Gauss transform and efficient kernel density estimation. In Proceedings Ninth IEEE International Conference on Computer Vision, pages 664-671 2003.
[51] Jiyan Yang, Vikas Sindhwani, Quanfu Fan, Haim Avron, and Michael W. Mahoney. Random Laplace feature maps for semigroup kernels on histograms. 2014 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2014, Columbus, OH, USA. June 23-28, 2014, pages 971-978, IEEE Computer Society, 2014.
[52] Tianbao Yang, Yu-Feng Li, Mehrdad Mahdavi, Rong Jin, and Zhi-Hua Zhou. Nyström method vs random Fourier features: theoretical and empirical comparison. In Peter L. Bartlett, Fernando C. N. Pereira, Christopher J. C. Burges, Léon Bottou, and Kilian Q. Weinberger, editors, Advances Neural Information Processing Systems 25: 26th Annual Conference on Neural Information Processing Systems 2012. Proceedings of a meeting held December 3-6, 2012, Lake Tahoe, Nevada, United States, pages 485-493, 2012.
[53] Manzil Zaheer, Guru Prashanth Guruganesh, Avinava Dubey, Joshua Ainslie, Chris Alberti, Santiago Ontanon, Philip Minh Pham, Anirudh Ravula, Qifan Wang, Li Yang, and Amr Mahmoud El Houssieny Ahmed. Big bird: Transformers for longer sequences. 2020.
[54] Yukun Zhu, Ryan Kiros, Rich Zemel, Ruslan Salakhutdinov, Raquel Urtasun, Antonio Torralba, and Sanja Fidler. Aligning books and movies: Towards story-like visual explanations by watching movies and reading books. In IEEE international conference on computer vision, pages 19-27, 2015.
14

---

A Proofs
A.1 Proof of Theorem 4.1
Proof. By the definition of UIDE:/DE, JD::JDE/ we have:
((1) ,(2)
EyDn (DF x) (w) = (2x) D2
D2 / cxp( ( wwll + 2 Aw
+ JuTB(U x + B 2) y) + x xC + T dw
= (2T) D2 exp (x TC() x + TOO0ly
Dxexp
w T (Ia 4A) w ( x dw.
x exp
JRd
Since 8A < Ia, we have 4A < 0.5 < Ia. meaning that I - 4A is positive definite and invertible. The following identity is straightforward to check:
T(a (La 4A) + T (B y) (w 11) - n²) + z/e -
T
2
== (Ia 11a-4A4 -A) μ= E + B)
Therefore, we have:
Evpe (w, p y) (2T) -4/2) exp xT y J.llyy juu 2ut
(-2 (w
(-2(w-1)
- 11) du.
x
exp
Rd
Next, we use the fact that the integral of the multivariate Gaussian distribution with mean μ and variance E is 1:
1 { exp - (w - m1) dw =1
(2m) 4/2 det 1/2
From that we conclude:
1 T N '1))
Eype (D (w. x) So = 122 det (et(())1// exp
x x T
= De 4A) exp x xx(0)x++ xy Tc(2ly
1 [ x x y) (I 4A)
-1 |B x + +B
2
exp x (C()) 1 (B (I. 4A) (B())
D2 det(Ia 4A) 1112
((L
ty ty((00)
1 1 (La T
B 2))
4A) x - (2))
B (1) (Ia 4A) -B))))(2) |B(3)y)
Based on this expression, we conclude that, indeed, EvDE (E (w x) / (w,y) k1t0) (x, y) for all
x,y E Rd if the conditions from theorem's 's statement are satisfied.
Next, we calculate expression for the variance. For any random variable Z, Var Z EZ² - (EZ)².
In particular, if Z = - (w, x) p (w. y), VDE. we get:
(
Var Fp pW (w. x) x2) (w, Evpnn Hf (w, x) (22 (w. y)² (Eunn/liius (w, ())
Eype (1)
(F (w, x)² 22) (w. (w.y)² y)² k-(0) (x,
15

---

We have:
Evps (w, x)² S (w. y)² (2T) ( D2 / exp ||01 + -w Aw 2w (BO) B
2x T C( 2y J c(tly dw (2T) D exp 2x c + +y T c
/ expl z (Ia 8A) + 2w0" (B B(3)y) dw. R
Evaluation of the integral above be done in the same way calculation of (12) Pl (w,, x) ( y), noticing that Ia 8A is positive definite and invertible. The result
Evng is as follows:
EVDR DE
(11) w. x)² (22)
JDE w. y) D det(I. - 8A) 1/2
X exp 2x (ca)) (I)) -A)-1B(1)) x
+2y +
+ (C*)
(La 8A) -A)-
+ 4x xx T T(a - 8A) (
We conclude that the variance expression given in the theorem's 's statement is correct.
A.2 Important lemma
Below, we prove important lemma which is used in the subsequent proofs:
Lemma A.I. Consider function f (-o0, } defined as
f(A) = loogl1 (-0(1-44) 4A) llo(1 8A) 1 8.A
(14)
where o 2 0. Then, the minimum of on (-o0, } is achieved at
A"
1
1- 26
120++1²
(15)
16
Proof. Set 7 (1 8A)- (0, +00). Note that there a one-to-one correspondence between
7 E (0, +oo) and A E (-00, ) Hence, we substitute 7 1 SA and 1 4A4
((1 8A) + 1)/2 = (r-1 1)/2 27
1+a in (14 and equivalently perform minimization with respect
to
1 2
min 7E(0,-00)
1(7)= = log
+ 2
log 7 + 07 log(7 + 1)
log - log 2 + 67.
2-
For h(-)'s derivative, we have:
--1
11(7) Φ + 7+1 27
= +
(16)
27(7 +
20 1) +
207² (26 + 1)7-1
(17)
27(7 +
27(7+1)
Based on (16), we see that 1'(7) + 00 as + and 11(7) > & 2 0 for all > 1. Hence, conclude that fi(-) is bounded from below on (0, +00) and the global minimum 7 on (0, +20)
exists and it is one of the points satisfying H' (>*) 0. Hence, it's one of the positive roots of the polynomial in numerator of 177.
o = 0, there is single root 7 = of the polynomial in the numerator of hence it is a global minimum of ((-). If ó 0, then there are two roots of the polynomial in the numerator of I1
" -(26
(20 + 1)2 86
4
16

---

2 -20 +1 + ((2 + 1)2 + 86
(18)
46
Note that, if Φ > O, then 26 + 1 > 0 and (2Ø + 1)2 + 86 2 (2p + 1) Hence, 22 < 0 and
7 > O. We conclude that >- = 7 is the minimum of h-) on (0, +20). We multiply numerator and
denominator of (8)'s right hand side by (2Ø + 1) + (26 + 1) + 86 0:
+ + 86) (26 (2) (26
2
(19)
46 (2 + 1)+ (2 1)2 + 86) 2Φ + + + (22 1)2 + 86
Note that the right hand side of [99 is equivalent to (8) when > 0 but also holds for the case
when Φ = 0 (i.e. when y" = 1). We conclude that /(.) is minimized at A _ (-*) ) since
7" = (1- -.-)-- It's easy to see that (15 follows from 19 directly.
A.3 Proof of Theorem
Proof. With A = Alg, the conditions from Theorem 4.1 read as
8A 1,
4A ( ( B IIa ((
1
(B()) ((()) D= = (1 - 4A)4/4
2(1 - 4A)
(20)
for k E 21, And the variance expression 99 for all x, E Rd transforms into
Var AADE(W, (WW y) = D²((1 8A)) d/2 exp
2x
((())
1 8A ) x
+2y ((2) + 1
84 (B ( y +
- 8A (B B(3l y) k( (x, 3)
4
We express ((k) through A, (k) and D through A using (0)) in the equation above:
Varpane I (11)
(w, A =
1 1 4A
exp (1 1 8A
1 1 4A x xTBB"")
VI 8A
2
1
y (B((2)) (2)) y
1-84 8A x T (B()) B(O)y - KK(0 (x.))
1 - 8A
4A
Since 1-4A
B (1)) (1))TB(2)
(2) is a full-rank matrix L 20). both (11) and (2) are full-rank. Hence,
can express ((2) = =1 ()) Also, note that
1
8- =(1 (1-88))-(1-
1- 8A
1 4A
(1 -1- 8A)(1 4A)
We rewrite the expression for the variance using the identity above and the formula for ((2);
(
1- 4A 1-44 exp (1 (1-- (A) (11)))))(0-4))=xx((((( T(())
Var (2) ADDE (w. x),
(w) (B)
/1 8A
-11-84) 44)y (BB
(B + 4(1 -8A) '(1) -4)x y (x,
We the expression above rewrite
(8 for (v, ((), {(2)} (VADE, S f(2 follows:
L(OADE; X, TADE) = - log( Vgg(VarD (1) (1 (())) {(0) y (0) + K (x(1), (x00),101))
IStJSL
1
= =dlog(1 - 4A) - log(1 - 8A) + (1 (1-88) 4A)- 1--1 Σ
xx () (B(1)) T )
2
1=1
+(1 (1--8A)-116 (11-)-LL----1
17

---

+4(1 8.)) (1 4A)L -2
Σ (e(')T((i))
(21)
IsUSh
Denote E (B T B Then (1 becomes:
d
A(DADE: N.Y, TADE) dllog11 4A)
log(1 -84)
2
1 ((1-8A) 4.A
Σ (x) T Ex(")
1=1
+(1 -
81) ((1-8A)-14
4A)L Σ
+ 4(1 8A) 4A)L Σ (x(")Ty()
IStJSL
(22)
We next prove the following lemma:
-12)) When A (8A 1) fixed. = E (B(")) B()) minimizes the right hand side of (2 with respect E.
Lemma A.2. Let B(I vI 4A
Proof. We have:
(x(1) (1) T Ex L
Traco((x(1) Exx(t)) (xx())))
---1
Trace Trace(Ex()) (xe())
1=1
Trace
E I. ( (x()))
Σ 1
-
Trace (EM
where we use the cyclic property of trace (-ace(-) and linearity of trace. Analogously, we obtain
1- (ybi)) ( E-1y(0) Trace E (2)), Assuming that A is fixed, optimization of 22
with respect to E reduces to the following minimization problem:
min E=0
(TTaae(EMMM)) F(E) 31 Trace EM ()) + (ETrace(E--MMM ((2))
(23)
where 31 = (1 8A) (1 - 4A) -1 3g = (1- (1 4.A) and the constraint E E Sq. >
follows from the fact that E = (B(1) ( and E is invertible. We have 8.4 >> 0 and
1 -A (1 8. 1/2 > Hence, 81, B2 0. For any E> and any A E Sy there is i R small enough such that E + tB is invertible and the following Neumann series is convergent:
E A -1
[L. + IE 1-
Σ 1=
(-t) (AE-))
We further deduce:
0 ((--))
- M
Trace (E (1)) Trace
(1=00
Σ
(- Trace
E EE-(AE- M2
1=0
and, therefore,
F(E+ta 81 Trace(EE (A)M -2 Σ (-t) Trace E 'AA ((AE-1'M(2)
1=0
= 3 Trace (Mace + 13, Trace (M) + B2 Σ
(-1) Tracee E (AE-A--))' (24)
1
18

---

Further, we have:
①
F(E +1 31 Trace AM! 82
1-1)t Trace (EE
(A
Ot
1m1
82
F E + 14)= = B2
(-1) 2
E (AE
(ōt}²
l=2
8
F(E + +1))
2/3 Trace E AE AE
AE
(25)
(ōt)2
1=
We replace M Q (A (2) 12) 2222
and apply the cyclic property of trace in 25
Ø²
F(B
(A)
-32 Trace (A 1/) E EE'AE" BE-'AEE-N (A (2)))112
(01)2
1=0
= 233 Trace
(TE TTT
where T (A (A Since E is positive definite, E-1 1 is also positive definite and a
TE-TT' at least positive semidefinite. Hence, Trace (TE 'TT Z and also
(#t) F(E
A )l === 2 0. We conclude that F(E) is a convex function on E E Sa E > 0). Since E E
Sa E > 0} is an open set, (every) global minimum E of 23 satisfies two conditions
1) E 0, and 2) VF E
Oaxd
(26)
Set t = 1 and assume that A Sy is small enough by norm so that E + A is invertible and the Neumann series for (Ia AE AE-1)-1 is convergent. Then, (24) holds for
F(E + A) = 31 Trace(EM + 31 Trace AM + 82
Trace
E E
(A'E (AEE)/M())
t=
= F(E) + B1 Trace (M + 82
- Trace E (AE
1=1
F(E) B1 Trace AM - Trace
E P'AEE-IN A
+B²
1) Trace
E ((AE)))(M2))
Clearly, 82 Eiz( (-) E (AE 1) o(|A|]) where || || is an Tg-norm. Also,
using the cyclic property trace, we get:
Trace E EE'AE (E
Trace AE- AE-IM@EE-1)
Therefore, have:
F(E + A) = Trace A ( 3. "ME B-1) + ((I
(27)
Since A, E M() ( Sd, from
27 follows that
VF(E) = ( 3.E (2)
(28)
Let E* = (B(1)) T B())+ 2 0. Note that
B2
(1 -1-88. -1 vI -1-46 44
B
8A) -41-
Since / + 0, Σ. U, A 1/2 ( are full-rank, E" is also full-rank, therefore E and
satisfies condition from (26). Observe that
E (A 1/2/22 V /3// -/) (A TT (A w)
19

---

VA
V "
where we use definitions of E", U, Σ, V and orthogonality of Q Q()) , U, V. Hence, we deduce that
31 E (A
1/1//
(1 "E")
(29)
31
32
Q
(A
11)1//2 (111)
31
-22
(30)
Q(2) A (2) and due orthogonality of V, U. By left- and right- multiplication of 30 by (E*) we deduce that
by the definition of Q() A (1)
(, (1)
or, in other words, VF E"" Odxd and the condition 2 from 26 is also satisfied. We conclude that
the global minimum of 23 is achieved at E
According Lemma A.2 (
(1)+ global minimum of (1 ss right hand side when A is
fixed. Indeed, if there B which leads to smaller value of 21 E (B(") ( would lead
to smaller value of (22 's right hand side. Also, this E is positive definite by definition (note that
(1)) is nonsingular), leading to contradiction with Lemma A.2
Substituting E" instead of E in 22 corresponds to the minimum value of A(OAGE; O, X, TAGE) for a fixed A. Our next step IS to minimize this expression with respect to A. Denote F (A () - (A (1)
(T Then E* = (1 4A)F where F doesn't depend on A.
We substitute E into (2) and get:
d 2 log(1 8.) + (1 8A) -1 4A)
dllo((1 4A)
(1-44)FM
+(1 - 4A)Traco((1 - -AA4-AF-1- (()
+4(1 - 8A) -1 4A)L Σ (x(1) Ty(i)
d
dlog(1 4A)
log(1 8A) + (1 -8A) -TTace (M + - (M
2
+2 (1 + (1 - -8A) (3)
(31)
where we also replace
L- Σ (x - Σ (1)
Σ )فو d)(
ISAJSL
and
(1 - 88) - 4A) = (1- 8A) + = (1+ + (1 -8A)-1)
2(1 8A)
Based on
30) and since F / /B'/BEE', conclude that FM M(2) M
F F- M (2) Using the cyclic property of trace, we get:
Trace(FM
Trace(M Trace(F-'M' (M))
By the definition of F, A (1) Q(1) and using the cyclic property and orthogonality of (1), U, have:
Trace (M) Trace (2"0 (2""0""0""""00000""")))))
20

---

Trace (Q ("""4"))1122330"""""1110"))
Trace ("
d
Trace Σ) =
1=1
Hence, 31 finally becomes:
dogg(1 - -A) log(1 - 8A) + 2(1 -8A)
Σ
En + + (1 - -8A) di
2
1=1
d log(1 - 4A)
1
log(1 8.4) + 2 - 8A)-1
d
1, x(3) + 22((3)
(32)
2
Next, we use Lemma A.1 (Ø = d ビ E + + 2 0) for deriving expression for A which
minimizes 322 This expression coincides with the one Theorem's 's statement. The expressions
for (2) C ((2) follow directly from (20). optimal (1)) =1))* and A. (10) follows from (3)
The proof is concluded.
A.4 Proof of Theorem
Proof. With ( B B and C(1) C() C, the conditions from Theorem 4.1 read as
1
8A < Ia, B (Ia T 4A) - La: C
B (A--4A)-1B
B Id.
D det 1
2
(33)
Denote
(La 4A) 41)-1/1BB 1/2B
Raxd Then, according to
33 QTQ LL, that Q E Og. We
rewrite
using 33 and then substitute B (I. 4A) 1/2 //2
Var SSDE (w, x) X)S ( y) = del 4A) del 8A)
(1)
112 -1/2 exp -||x|12
+2x B T 31( (Ia- --)--A)---- Bx 11y12 +2y B T (L. 8A) By + 4x B (I. 8A) -By-
- -k (x, del(Ia - 4A) (Ia 1/2 exp (-1X/]² 22x -2 Q EQx 11y112
-2y T
EQy T EQy)
(x,)²
(34)
where we denote:
E = - (Ia - 4A) 1/2 (1a -8A) (Ia 4A) - 4A) 8A)
1
2 (II 8A) + L1) (Ia 8A) Is [I. 8A)
(35)
which is in D. since A E D.. Next, we observe:
2x QxQQ_QQ + 2y Q EQy + 4x Q EQy 2(x + y) Q EQ(x y)
We plug this into 34 and use the resulting expression to rewrite 8 for (v. (1), f(2))
(1) (USDE, JSDE: > as follows:
[(OSDE; TTDR) = L-2
Σ
/ (w, x (1)) ) (w, ااشار + K (x(1)x ,
IS
1
= (Tt 1
log del - -8A - L-1 Σ
-- Σ 11
2
1=1
3=1
-2L
Σ
(x(+) + +01 T t0))
EQ(
(36)
1SiJ5L
21

---

Using linearity and cyclic property of trace, we deduce that
1
Σ (x + +y())) Q EQ(x + y())
ISiJSL
1,-
Σ
Trace (x(0 +
+ +y(11)
151J5L
L
Σ Trace (Q
(Q + (x() ((Oiir]
L- Σ (xl)) x(0)) + ((O)f
Trace
Q
QTEQ
Observe that
L- Σ (x(() + xx01)(x(()) + +(
ISUSI
=
Σ
(x)
+ x (5(6)) + (x)) + (x))))))
IStJSL
x (xx (0) T + - 1
L (1
{ yby
+ 1 (
(1T x
T
L
TT
(=1
1=1
j=1
j=1
+L 1 Σ
y (1 (yb)) M( + (1) (5) T + p (5) (x(44)) + (
j=1
Denote N (M + 1 (μ (5)) + (x((4) + ( We conclude that
(10sDE::,)) TSDE) = log (Ta 4A)
1
log det (Ia - 8A) T,
Σ 1=1
10
2
112 -Trace (Q (QNN..).
(37)
With A fixed, we minimize the right hand side of 37 with respect to Q which equivalent to minimizing L(OSDE; X, 3, TSDE) with respect to with fixed A, since there one-to-one correspondence between B and Q This is equivalent maximizing, again using the cyclic property of trace,
Trace (Q T (ETENN Trace (EQNQ")
(38)
with respect to Q. We prove the following lemma first:
Lemma A.3. Suppose that diagonal entries of E are all distinct, and the same holds for A (3), Let
II E {0, 1 be permutation matrix sorting diagonal entries of E (i.e. by applying IEEI in descending order corresponding to a permutation E Ndd, Set Q II (Q(3))T E O.. Then have:
Trace (EQ (Q²))
E Σ [=1
m A (3)
(39)
sup Trace (EQNQT)
(40)
QEta
Proof. First of all, we have:
Trace (EQ) = Trace
EIT (Q (3)) TN 1I)
Trace EEI A 11)
(41)
22

---

Trace IEI (3)
Σ
Em A
(42)
L.e. (39) satisfied.
Optimization for finding supqeos Trace (EQNQ) is a well-studied problem 9 By the definition,
A (3) has eigenvalues of N on the main diagonal and E Da hence contains its eigenvalues on its main diagonal. Then, proven in [9], Q" is indeed global maximum of this problem in the case of distinct eigenvalues for E and N. That is, (40) is proven.
Next, we prove a generalization of Lemma A.3 when diagonal entries and A (3) are not necessarily distinct:
Lemma A.4. Let II E 10.1d be a permutation matrix sorting diagonal entries of (i.e. by applying IEN in any non-ascending order corresponding to a permutation π E N Ser =' (Q (3) E Og. Then we have:
d
Trace (EQ (BQ'N(Q))) = Σ
E Ex A
(3) 1
(43)
1=1
sup Trace (EQNQ)) QeCa
(44)
Proof. In the same way as (4 42 we show that Trace (EQ"N(Q")T) Σi= Em (43) satisfied. Next we prove that for any Od.
A i.c.
(3)
Trace (EQNQT) S
E B A (3)
(45)
which would imply
Our proof is by contradiction. First of all, we can assume that E, A (3) are nonzero matrices since otherwise we have 44 trivially. Since Trace (EQNQT) is continuous function of Q and Og is Trace (EQNQT) is finite. Suppose that there > 0 such that
compact, supqeo,
6 sup Trace (EQNQT) QeOa
(3) Σ E E AL.
(46)
~ (3) Let E, A
Da be matrices with all distinct values on the diagonal such that
E E < min
|
8
| (3)
A (3)
(47)
12
(3)
12
where || denotes Frobenius norm and | ||A (3)
|| 0 since these nonzero matrices.
Further, we assume that diagonal entries of are sorted in a descending order and, addition to A (3) π also sorts entries E in a non-ascending (descending) order. Clearly, such E, X (3)
can
be obtained by small perturbations of E, A (3) Also, denote Ñ = a (3) Ä (3) (Q(3))_ Since Og is a compact closed set and Trace (EQNQ) ) is a continuous function of Q, there exists Q O, such that
Trace [T)) sup Trace (EQNQT).
(48)
QeO,
By the definition of E, N, we have:
Trace EEQ () Trace (TQ')T()))T)
=(Trace ( Trace (EQ")(Q)")T)) + (Traee (TQ))(Q"))T)))
23

---

Trace TEQ"N((C (EQ"NN(Q"")T
Trace EQ" (N N) (Q7)T) + Trace (1-E)QNN(Q~11)
Next, we apply Cauchy-Schwarz inequality to both terms:
Trace EQQ (N - N) (QM)T S II(Q)* E E - (2 |||| Nllr,
Trace (
Q NT))(""))) = (I(QQ*)))))
(BB[)""(_BB)))((((""") = E E [
where we use invariance of the Frobenius norm under multiplications by orthogonal matrices. Using this invariance again, we deduce that
IN
||Q (A
A (Q (3) T | (3) A 1 l|., (3)
[N = 1Q (Q (3) T = | (3)
l|f.
We conclude that
Trace (EQ"N((Q"") < Trace [EQ"N(QQT)) + ||/AA A
1| E
El
E|[r//
(3)
1|--
(49)
A.3 to E = E, (3) (3)
Next, we apply Lemma
and deduce that
d
Trace (EQ"N(Q")T)
<
म. m
A
(3)
=
Σ
Em A + Em. Emg
A
(3)
1=1
1=1
(3)
A.1) )
+ +E m
A
d
E
A
(3)
E E
E
(3)
+
E E
(13)
A1.1 L,5
[=1
1=1
1=1
E A
A
Trace
I[E-E))
Trace
HEI
(A (3)
- A (3))
We apply Cauchy-Schwarc inequality again to the second and the third term:
Trace I EJIT (3)
W E EII
(3)
|E =
E
Elle A ||..
Trace
ПEN (3) A A0)
3
(3)
where we use invariance of the Frobenius norm under column and row permutations, We conclude that
Trace (TQ')(Q")))) < Σ E A +
E
E||||A
+
E |
(3) lI.
We combine this inequality with (49 and obtain:
Trace (BQ"N(QQ"))) S
E m A +2
(3)
|
E/|r||A (3) 1|r. (50)
[=1
Next, we use triangle inequality and deduce that
(3) |F
|
(3)
||
(3) Ir.
Hence, we continue
50
Trace (EQ))N((()
< Σ
E
F....
A
(3)
2 AA
24

---

A
+ 2 -Elle
(|A ||E 1 (3)
|
A 1
E
(3)
+
E 33)
A
|F + + (3) || E
E/F
2
(3)
-3)
1|F
(3) E B 1
2
A
(3)
(3)
|| + 2|| (3)
F |
Elip +
E|||||||||
A 1
|F
where in the last transition we use
E
E I
E
which is according
47
We continue this
chain of inequalities using (4) again:
Trace (EQ"'NQ²))T) S
Σ
2
+
+
2
8
(3) 1,
+
6
12
12
12
E w.m
A
(3) 1,1
+
1=1
This is contradiction with (46) taking into account Q**ss definition (48) Hence,
(44 proven.
π,, (I) II where π, II are defined as in Lemma A.4 statement. That is, π (A)) denotes some permutation which sorts diagonal entries of E in non-asoeee ascendicnndn order, It's a function of since E is function of A defined in 35 In fact, based on (35). we see that π (TA) is some permutation which sorts diagonal entries of in a non-descending order. (I(A) denotes permutation matrix corresponding π (A). That is, diagonal entries of (IIA) are sorted in non-descending order.
Let Q" be defined as in Lemma A.4 statement. Further, we denote T (A)
Let (A) denote the right hand side (37) where we substitute Q = =Q" That is, G(A) is an optimal value of Z(OSDE; X, V, TSDE) with A fixed:
L
L
((A) log det 4A)
log dea SA) --1
(x11
1-1
Σ
111
5=1
d
2
E (T)))
L
L
log det 4A)
log det(Ia -8A) --1
Σ
||x
(1) 1
L-
Σ
11
2
5=1
(1 (A A
(3) -
(51)
where E's definition 35 Let n-- (A) E N' denote permutation inverse T(A). By rearranging terms in the sum, we have:
Σ (L+(I-SAAA)AA))))
(1+ -
(A)x
Therefore, we have:
C(A) log det(I. - -AA) log det 1-1
1
xx 10 1-1
111
1=1
(1+ (1 -8AA) A 1
(52)
25

---

Define new function G(A, A), where A E Dd satisfies 8A x Id, as follows:
g(A, A) = log det (Ia 4A) dea 1 8A) + Σ (1 8A()) AA A (3)
x -A), (A))
By the definition of G(A, A),, we have:
d Σ
((A) G(A. L
Σ
xx
1 --1
Σ
A
(3)
(A)i, I(A),"
Hence, it holds:
G(A) inf G(A,A)-L
Σ
| 12
Σ 10
(3) n (A).. (A)))=((AA)) (53)
SA-L.
Next, we show that there is a closed-form expression for the solution of inf inf AED 8A <t, 9(A, A).
Since A E Dd. we have: log det(Ia 4Â) = Ei Li=l - 4At,), log det(I. 8A)
EL log(1 8Ar,)). We further have:
Σ (lo(1--AA,)) (A.,) log(1 - 8A11) + (1 8Â,,)) -AA (3) - (A)) (54)
G(A A)
From we see that minimization inf AED, sÂ <I, ( A) reduces d independent minimization problems with respect to At, such that 8A1, 1. I'th problem. 1 I S solved using Lemma
A.1 where we set 6 = A - (A)) Let A E Da denote the corresponding solution. Then,
(3)
for all < 1 S d, we have:
A
- 24 (3)
(3) -V (2 (2A +1) + 8A (3) * (A)( (A);
(A)
(55)
From 53
it follows that
L L-1
d Σ
((A) 2 9(A"
A)
L-1
|x(
10
Σ
11
A n 1(A). *(A):
(3)
5=1
1=1
1 log det(LI -AA") - --1
dogdet(L. - 4A**)
xx (1 1-1
Σ
11
2
(1+ (1 (1++1-8A7)) 84 A 33)
(56)
"(A) (())'
1=1
Denote Ex = I 1 (0A-8A")-- Then we have:
Σ
- (1+ A (A).. ((A), 2Trace E A"'I(A))- -AA33 (3) (II(A)-1) 7)
S
2
E A (A))))
(3) L.
(57)
where the second transition follows from Lemma A.4 and the fact that T(A**) sorts diagonal
entries of A ** in non-descending order, hence its sorts diagonal entries of Exx in a non-ascending
order (recall the definition of ( **) and E**). Denote E II(A* >= II(IIA (T. Then
Exs A (A")) (Am)))))) = Ej for all 1 < < d and
A
E
(3) A
(58)
1=1
26

---

Further, we have:
- - - (--)8- - I(A))))
E" = IIAA
2
-
1
[L 8 -(
2
I
[I. 8A")
2
where we denote A = IT A II TT, i.e. A A Ai for all 5/5 d.
Given the definition of A
55
for all 1 < E d we have:
(3) (2A + +1) 8A
A 16
1 - 21
(59)
(3) That is, A is independent of A. Based on 55 we that smaller values of 1 result in (3) bigger values of A Since A 1.1..... A (3)
d.d ordered in non-ascending order, deduce that Ai A....... Ad,a are ordered in a non-descending order. By the definition of (A*)), we then have
A A ((A (A))))))))) for all 1 < d. Therefore,
E = Σ
(1+ -1-881)))))) -A =
8A (1
Combining this with (58). 57 we can continue the chain of inequalities
(56):
L
((A) 2 log det(I,
4A*)
log det (AAI(-8AA") - L-
xx (
L-
Σ
1120
2
i=
j=1
d
+ (1 8A A )( A,
1=1
log (La -44 4A") log det 8A+) - 1-1
1x L
Σ
|
1=1
+ Σ 1 + (1 (1-118A)) 8A (AA))))) 1-1
A = G(A")
(60)
where in the second transition we use the fact that
dct(I. - -AA = det (II(A*) [I. 4A IAA")) det (I2 4A")
definition of g(-) 51. Note that (60) holds for all A E Da such that 8A < Ia and also 8A* < La since SA" La We conclude that, when B. C, D are chosen optimally with given A. the minimum of C(OSDE; X,Y, TSDE) is reached when A = A" As we have already deduced, diagonal entries of A = A are sorted in the non-descending order. Hence, using Lemma notation, diagonal
and, similarly, det(I, 8A") det (I, -A""). In the third transition,
entries of E are already sorted in non-ascending sorting order and π (1,.... d). II La satisfy
requirements of the Lemma. Hence, with A = A", the optimal B has a form (I. 4A 1/2 where
Q' La(Q ) (Q)3) T. Optimal C and D are further determined by
33
follows from
(60 and the fact that, as discussed above, we can replace π (A )) with in
60 S d.
B
Additional experimental details
B.1
Compute resources and implementation
We use NumPy 25 in Google Colaboratory the variance comparison and kernel classification experiment. For the Transformer setups, we use TPU cluster and JAX [7] library. All tested Transformer variants were trained and tested on TPU pods containing TPU v3 chips with JAX and GPUs (V100).
27

---

B.2 Variance comparison
We repeat the setup of [31] closely: we draw pairs of sets {x(") hsis hylss 1024. On each pair, we compute the relative variance for all pairs of points and for all indicated RF methods. Further, the shifted log-variance optimized separately on each pair of sets for GERF, ADERF and SDERE.
We take M 1 since M's value important in this experiment: bigger M would just shift the curves below. The reported curves are means over all pairs of points and over all 5 sets.
Table 2: Hyperparameters for the base models for pre-training for all methods
Parameter
Value
of heads
12
# of hidden layers
12
Hidden layer size
768
# of tokens
512
Batch size
256
256
Pretrain Steps Loss
1M
MLM
Activation layer
gelu
Dropout prob
0.1
Attention dropout prob
0.1
Optimizer
Adam
Learning rate
10 10-4
Compute resources
x TPUv3
Table 3: Dataset used for pre training.
Dataset
# tokens
Avg. doc len.
Books [54]
1.0B
37K
Wikipedia
3.1B
592
B.3 Kernel classification
As in (3I. we obtain training, validation and test splits by shuffling the raw dataset and taking 90%. 5%. 5% objects respectively. The splits are fixed for all RF methods. We tune a logarithmic grid of 10 values on [10-2, 2, 10 For each and each RF type, we try 50 seeds for drawing RFs during validation and testing. Testing is performed for the best a only. Figure reports averages over 50 seeds. We use orthogonal w's for all types of RF' as described in (31), since orthogonal random features work better in practice [16 31].
B.4 DERFs for long-sequence Transformers
B.4.1 Speech modelling
Our Conformer-Transduucee variant was characterized by: 20 conformer layers, model_din = 512, relative position embedding dimensionality rped 512 and h 8 heads. We used batch size bs = 2048 and trained with the adam optimizer on TPUs. For the regular Conformer-Tranne Transducer training, we ablation studies different number of random features: NST setting, we run experiments with m 8. We reported commonly used metric: normalized word error rate (NWER).
8, 32, 128. In the
B.4.2 Natural language processing
We pretrained BERT model on two publicly available datasets (see: Table 3). Following the original BERT training, we mask 15% of tokens in these two datasets, and train to predict the mask. All methods were warm started from exactly the same pre-trained checkpoint after IM iteration of BERT
28

---

pretraining, We used the exact same hyperparameter-setyy for all the baselines (FAVOR++(3) [31], FAVOR+ [o] ELU [2),, ReLU [T]] and FAVOR++ The hyperparameters for pretraining are shown in Table 2 We finetuned on GLUE task, warm-starting with the weights of the pretrained model. The setup is analogous to the one from the original BERT paper.
29


Reviews:
[{"review_id": "hGoK9lTzDH", "reviewer": "Reviewer_r6PB", "summary": "This work focuses on positive linear features for softmax kernels which are relevant for accelerating kernel methods and transformers. The paper observes that the functional form for random Fourier features (GERFs) proposed by prior work in [30] can be seen as optimizing the shift-ed log-variance objective. Building on GERFs, a more expressive and parameterized functional form is proposed, and corresponding shifted-log variance is derived. After that, several simplifications are considered for improved optimization of the log-shifted variance of the random features, which leads to different forms of RFs that are easily computable while keeping subquadratic complexity in the number of sequence lengths. Towards the end, a few empirical results are presented in support of their contributions. ", "questions": "* Is it possible to repeat the proof of Theorem 4.2 for a setting where the log of features comes from a non-linear function, potentially a neural network? The main idea is to extend from the quadratic functions of w and x, use a more general parametric form, and optimize a relaxation of the variance.  Note that to ensure the unbiasedness of the kernel, estimation can be imposed using constraints similar to the proposed work. How would you compare and contrast with the current approach? If not, why is it not possible?\n* Why are time versus accuracy results not reported anywhere, assuming scalability is one of the main motivations? When will these methods accelerate, and when will they not? Theoretical complexity may not account for constants that stem from implementation. \n* What aspects of the implementation may prohibit practitioners interested in the fast self-attention method from considering these proposed methods?\n* Respond to weakness 1.", "limitations": "Please take a look at the weakness and questions.\n", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "* This work studies an important class of problems that has broader applicability. The paper builds on observations/heuristics from the prior work and attempts to justify using theoretical arguments. During this exercise, the proposed methodologies are mathematically elegant and present several exciting results on RFs that are relevant to both kernel methods and transformers.  \n* While I have not examined every mathematical proof, the derivation appears generally accepted, insightful, and praiseworthy.\n", "weaknesses": "* To establish that GERF minimizes the objective, the supporting argument is the solution of 8 matches with the heuristic proposed by the [30]. If I have understood it correctly, more is needed as for an objective to be reasonable, its behavior beyond the maximizer/minimizer needs to be understood. Also, is there any relation between this objective and the actual variance of the GERF estimate? Can Jensen inequality be used? This is one of the main results of this paper, and it needs to be better motivated. \n* Notations could be improved. E.g., use a consistent format for scalars and matrices. L is a scalar. Avoid unnecessary precision; e.g., the upper subscript in lines 157 and 158 could have been avoided and stated as part of the text. Note that this work is built on math and therefore requires extra effort on styling to make it accessible for a larger audience. \n* Pages 5-6-7 are not used carefully and do not present the best aspects of this work. You are considering simple formulations and repeating more or less the same procedure for obtaining close form solution of the objective, which again needs to be justified better in the first place. Also, before going into these details, give an overview of possible choices on the parameters space and justification. Summarizing these results as a corollary instead of Theorem 4.3 and 4.4 might be helpful. \n* In line 217: “operations for which implementation has not yet matured in popular deep learning libraries with GPU and TPU support.” What is the evidence for this claim? Are you suggesting that eigen decomposition can’t be computed efficiently on GPUs using torch and TensorFlow? \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for all the comments and kind words in strengths! We address weaknesses and questions below.\n\n> If I have understood it correctly, more is needed as for an objective to be reasonable, its behavior beyond the maximizer/minimizer needs to be understood.\n\nAn intuitive understanding of this objective is that, if we assume that variances over different pairs of points are tightly concentrated around one point, then the objective is the logarithm of that point. Jensen’s inequality can be used to lower-bound the variance, however not to upper bound it since log is a concave function. Note that our experimental evidence (Figures 1 and 2) suggests that the minimization of (8) indeed leads to variance minimization where we get up to e^10 times variance reduction.\n\n> Notations could be improved.\n> Pages 5-6-7 are not used carefully and do not present the best aspects of this work.\n\nWe will incorporate notation and style suggestions by the Reviewer in the final revision. Thank you!\n\n> Are you suggesting that eigen decomposition can’t be computed efficiently on GPUs using torch and TensorFlow?\n\nWe used the Jax codebase and we experienced errors and unhealthy behavior when using SVD and eigen decomposition in Jax.\n\n> Is it possible to repeat the proof of Theorem 4.2 for a setting where the log of features comes from a non-linear function, potentially a neural network?\n\nThe quadratic nature of the problem is essential for our close form solutions, meaning that we don’t see how our proofs can be extended to arbitrary nonlinear functions. Perhaps approximate solutions are feasible, we leave that to future work since it’s outside of the scope of this paper.\n\n> Why are time versus accuracy results not reported anywhere, assuming scalability is one of the main motivations? When will these methods accelerate, and when will they not?\n\nOur methods’ complexity grows linearly with the sequence length as opposed to quadratic growth for the standard self-attention. Hence, we are guaranteed to get improvements for long sequences (of order 1000 and above). We emphasize that, apart from empirical contributions, our paper is mainly theoretical and provides a significant extension of random features for the Gaussian kernel with nontrivial close form solutions (Theorems 4.2 and 4.3). We ask the Reviewer to take that into account.\n\n> What aspects of the implementation may prohibit practitioners interested in the fast self-attention method from considering these proposed methods?\n\nFor smaller sequence lengths of order d (d is usually 64 and is the dimension of the attention head) our methods won’t give efficiency improvements over standard self-attention since O(LMd) will be close to O(L^2 d).\n"}, {"review_id": "DCgLsfBCgq", "reviewer": "Reviewer_pnU2", "summary": "This paper is the next instalment in a series of works focusing on low-rank transforners, e.g. [15, 30]. They propose Dense Exponential Random Features (DERF) for unbiased Monte Carlo approximation of Gaussian or softmax kernels. This class of features generalizes GERFs from previous work, with the main difference being that the features are parametrized by matrix-valued parameters rather than scalars. Hence, this flexibility allows for capturing a larger class of features for kernel approximation.  The parameters themselves can be analytically optimized in several useful special cases: ADERF, SDERF and SADERF. The parameter optimization is carried out with respect to a so-called shifted log-variance objective averaged over the dataset, and it is shown that optimizing this objective formalizes a \"homogeneity heuristic\" from previous work, which is a nice observation. In particular, this allows to select the same set of parameters across all data points. The experiments are essentially the same as in previous works: variance analysis on synthetic datasets + cifar + mnist, a set of nonparametric kernel regression benchmarks, speech modelling and low-rank uptraining on NLP tasks. Improvements compared to previous work and scalable transformers are achieved.", "questions": "Some questions and remarks:\n- Is $\\sigma$ analogous to the bandwidth parameter of the kernel in Figure 1 (right) and Figure 2? This is not stated explicitly.\n- In line 138-139, it is stated that [30] achieved good results by optimizing the shifted log-variance. Given that this objective is only introduced in this work, maybe I would phrase this differently, since in previous work it was only incorporated as a heuristic for homogenizing the solution, rather than an optimization problem.\n- How significant is the $O(d^3)$ computation cost associated with SDERF and ADERF? I guess the previous layer could easily have an overall dimension in the 1000s, where this could become a bottleneck on GPUs in terms of memory? \n- In line 214, the authors make the remark that assuming $L \\geq d$, the $O(Ld^2)$ cost dominates the $O(d^3)$ cost. Is this not already assumed in Theorem 4.2 (and hence, I assume in Theorem 4.3 from the way it is phrased), since this is necessary for the nonsingularity of $M^{(1)}$ and $M^{(2)}$?\n- If the main improvement of SDERF compared to ADERF comes from allowing a non-isotropic $A$, would it be beneficial to define extensions of GERF or SADERF, which add this flexibility to the parametrization? What might be the computational cost of such an RF construction? This is just a hunch, but if this can be solved for without matrix decompositions, then we might have the best of both worlds? No problem if working this out is difficult, just curious.\n- What is the range of sequence lengths in the speech modelling benchmark? Does this count as a long-range task?\n- In the NLP task, what is the number of random features (number of MC samples) used? Is there any intuition, investigation that the authors could report about how to choose this hyperparameter?\n- One more question; in the Appendix it is stated that replacing the standard normal distribution on $\\omega$ with an orthogonal one works better in practice. is this achieved by QR decomposing a set of Gaussian vectors? Is there any intuition about why this improves the performance?", "limitations": "The paper is mainly theoretical, but it has direct implications regarding Transformers, which carry with themselves a variety of societal and environmental impacts, but this seems to be addressed.  ", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper is overall well written, it does really well at communicating the main issues and ideas on which it is based. There is both a theoretical and experimental component, both of which have some interesting results. The main idea is novel, and the proposed approach for generalizing previous random feature methods is non-trivial; I was impressed that the analytical calculations could also be carried out in the matrix-valued case. I believe it is of interest to both the kernel and efficient transformer communities to have some scalable and expressive random features, and the paper definitely takes another step in this direction.", "weaknesses": "Overall there are not that many weaknesses in my opinion. The impact is slightly less significant than of its predecessor [30], where the main novelty was allowing for low-rank uptraining of pre-trained Transformers, which was already achieved by FAVOR++ (efficient Transformer based on GERFs). Nevertheless, the improved experimental performance should be of interest.\n\nI was missing some more comparisons between the variations: for one, SADERFs were not included in any of them, and only used for the last NLP task. Computationally, this variant seems to be the most favourable since it does not require the eigendecomposition of a $d \\times d$ matrix, hence foregoing an $O(d^3)$ computational cost, which is associated with ADERF and SDERF. It would be interesting to see this variation included in the comparisons to have an idea about the trade-off between the performance and this computational saving. \n\nI was also missing some intuition or hypothesis regarding why SDERF seems to perform best. Is it because it allows for a non-isotropic matrix $A$, hence allowing to adapt the variance of $\\omega$ on a coordinate-by-coordinate basis in the quadratic form containing $\\omega$, if this makes sense? It would be interesting to know which additional parameters result in the biggest improvement compared to GERF, so that focus could be placed on optimizing the scalability-performance trade-off. (See questions for more)\n\nAnother question that seems unaddressed to me is how much information is lost compared to a full Transformer on truly long-range tasks. At the moment, we do not know the \"price\" we pay for the subquadratic scalability in sequence length, and it would also be interesting to include some more long-range tasks, which compares full Transformer vs RF methods. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to sincerely thank the Reviewer for all the comments and a favorable score. We address weaknesses and questions below.\n\n> I was missing some more comparisons between the variations: for one, SADERFs were not included in any of them, and only used for the last NLP task.\n\nNote that we report ADERFs in Figure 2 which can be thought as a stronger version of SADERFs. It doesn’t perform worse than the previous best variant (GERFs), however SDERFs is even better. For that reason, we evaluated SDERFs in the speech modelling experiment, however we weren’t able to use it in the NLP setup because we experienced unhealthy behaviour when using SVD and eigen decomposition in Jax. That’s why we chose SADERFs in this setup and in general we recommend using it when matrix decomposition is infeasible.\n\n> Another question that seems unaddressed to me is how much information is lost compared to a full Transformer on truly long-range tasks.\n\nWe refer to [14] which is the original paper proposing RFs in the context of Transformers and demonstrating superior performance compared to full Transformers. We note that RFs proposed in our paper are a stronger variant of FAVOR+ from [14].\n\n> Is sigma analogous to the bandwidth parameter of the kernel in Figure 1 (right) and Figure 2?\n\nWe don’t parametrize our Gaussian kernel definition (line 58) but sigma in Figures 1, 2 is equivalent to the inverse bandwidth parameter in the standard Gaussian kernel exp(-||x - y||^2 / (2 bandwidth^2)) since the arguments to (our) kernel are sigma x, sigma y.\n\n> Given that this objective is only introduced in this work, maybe I would phrase this differently…\n\nThanks for the suggestion, we will incorporate it in the final version to mitigate a potential confusion. What we meant is that inherently [30] optimize Eq. (8) without knowing it and get good results which suggests that (8) is a good variance proxy.\n\n>  I guess the previous layer could easily have an overall dimension in the 1000s, where this could become a bottleneck on GPUs in terms of memory?\n\nNote that d is not the dimension of the previous layer but the dimension of the attention head which is much smaller and is typically 64.\n\n> Is this not already assumed in Theorem 4.2 (and hence, I assume in Theorem 4.3 from the way it is phrased) … ?\n\nThat’s correct, thank you for your observation! We will add this clarification to the Theorem statements.\n\n> If the main improvement of SDERF compared to ADERF comes from allowing a non-isotropic A, would it be beneficial to define extensions of GERF or SADERF, which add this flexibility to the parametrization?\n\nThank you for suggesting this, this could be a nice extension of our methods and we leave it to the future work.\n\n> What is the range of sequence lengths in the speech modelling benchmark?\n\nFor speech modeling the max sequence length was ~900.\n\n> In the NLP task, what is the number of random features (number of MC samples) used? Is there any intuition, investigation that the authors could report about how to choose this hyperparameter?\n\nThe number of RFs is M=128.  Bigger M can only decrease the variance of the approximation, hence we recommend setting M as higher as possible according to the given computational limitations. In our NLP experiments we tried 64, 128 and 256 and found that 128, 256 perform similar to each other while being a bit better to M=64. Moreover in Figure 3 we have shown plots for how accuracy changes across different dataset as we vary M. \n\n> One more question; in the Appendix it is stated that replacing the standard normal distribution on omega  with an orthogonal one works better in practice. Is this achieved by QR decomposing a set of Gaussian vectors? Is there any intuition about why this improves the performance?\n\nThere is a thread of work about orthogonal random features and why they reduce variance for kernel estimation: see “Orthogonal Random Features”, Yu et al. 2016, [13, 14, 15] and their references/citations.\n"}, {"review_id": "13M3Yy3VPM", "reviewer": "Reviewer_LgLu", "summary": "The authors propose new random features for Gaussian and softmax kernels, and apply the approximation to learning scalable Transformer networks.", "questions": "How feasible are the conditions in Theorem 4.1?\n", "limitations": "N/A", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The proposed dense exponential random features (DERFs) generalize the current positive random features (PosRFs) and  generalized exponential random features (GERFs).\n2. The authors show stronger theoretical results for the proposed DERFs in scalable Transformer networks where the self-attention matrix is approximated as a low-rank matrix when the sequence is long.\n3. The utility of the method is demonstrated on a variety of datasets with competitive results.\n", "weaknesses": "1. While the authors argue that significant variance reduction can be achieved, it is not clear whether the proposed approximation is unbiased.\n2. The method could be computationally expensive---the term D involves computing the determinant.\n3. Theorem 4.1 makes a few assumptions, restrictions of these conditions are not fully discussed.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the review! We address weaknesses and questions below.\n\n> While the authors argue that significant variance reduction can be achieved, it is not clear whether the proposed approximation is unbiased.\n\nThe proposed approximation is precisely unbiased according to Theorem 4.1.\n\n> The method could be computationally expensive---the term D involves computing the determinant.\n\nNote that A is a dxd diagonal matrix in all our DERF variants, hence computing its determinant takes O(d) time.\n\n> Theorem 4.1 makes a few assumptions, restrictions of these conditions are not fully discussed.\n\nAll our close-form optimal solutions (Theorems 4.2, 4.3) satisfy constraints from Theorem 4.1 as mentioned in lines 176-177 for Theorem 4.2 and 198-200 for Theorem 4.1.\n\n> How feasible are the conditions in Theorem 4.1?\n\nAs we mentioned above, these conditions are satisfied automatically for our optimal solutions.\n"}, {"review_id": "x12jfGWHoP", "reviewer": "Reviewer_PE4H", "summary": "In this paper, the GERFs are genearilzied to dense exponential random features (DERFs). The paper shows that GERFs and PosRFs are the specific situation of DERFs, from which if follows that with a suitable parameter estimator, DERFs can have better performance.", "questions": "In different experiments, several random features are compared when the number of random features are the same. How about the real calculation time? Since there could be difference on calculating different random features in the inference. \n", "limitations": "not applicable", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "As far as I know, the DERFs are novel and are a good extension from GERFs.\n\nThere are sufficient discussions and experiments from many aspects, all showing the good performance of the proposed method. \n\nThe paper is well organized and clearly represented.", "weaknesses": "I did not find significant weakness. One possible is that because this paper include many aspects, it may make the reader lose the focus. \n\nIt is better to clearly express the suggestion: in which case DERFs is recommended.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the Reviewer for a high score and kind words in “Strengths”! We address weaknesses and questions below.\n\n> It is better to clearly express the suggestion: in which case DERFs is recommended.\n\nWe recommend using them whenever the sequence length is too prohibitive for using exact kernel matrix/self-attention.\n\n> How about the real calculation time? Since there could be difference on calculating different random features in the inference.\n\nThe computational complexity for all RF methods scales linearly with the sequence length, hence, for a large sequence length, we consider all these methods as efficient and only compare the downstream metric performance.\n"}, {"review_id": "Fz6lLxAPnC", "reviewer": "Reviewer_z18t", "summary": "This paper studies the problem of computing the matrix \"KC\". Here, K is the kernel matrix (L-by-L, with L very large, could be reproducing kernel, but OK otherwise), and C is a known constant matrix. The kernel is restricted to be \"scaled softmax kernel\": exp(x-dot-y)F(x)F(y) (for some function F not difficult to handle). The contribution of this paper, is the design of a class of random features, of which the correlation is exactly K, but to some extent, the variance of replacing the mean (for the correlation) by sample mean, can be minimized with closed form. This research, to the best of my knowledge, is innovative, inspiring, and interesting.", "questions": "1. Lines 76--77. Could you please provide more details on why the computational complexity of hat-K-times-C is O(LMn)?\n2. In Eq. (8), the variable \"theta\" does not show up in the right-hand side, and is not elsewhere defined. Better if this variable can be defined.\n3. It remains very unclear even after a careful study of Section 3, what are \"homogeneity heuristic\" and \"a certain optimization problem\" mentioned in Lines 124 and 125. Better if these can be explicitly pointed out.\n4. Consider the mission of computing the matrix KC, as elaborated in Section 2.1, whether adopting DERFs would reduce computational complexity? To be specific, whether the amount of variance reduced, deserves the computational cost incurred by optimizing the random features via the approaches introduced in Section 4? Can one achieve the same accuracy and computation complexity simply by using a larger M? For example, Theorem 4.2.", "limitations": "foundational research, limitations well managed.", "rating": 3, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "In this paper, the idea, and the derived formulation of optimizing sample variance for improved kernel estimation is innovative. Besides the notation problem raised below, the paper has a clear and friendly structure that is easy to follow. The main results are presented clearly. These suggest excellent writing quality and clarity. We do have concerns on significance, but we would like to raise such concerns below, instead of here.", "weaknesses": "1. Concerns about the writing. I feel that the notations are unnecessarily too heavy. For example, I can not see any reason that blocks one from rewriting:\nf^(k) --> f^k (or just f and g), B^(k) --> B_k, C^(k) --> C_k or C^k, x^(i) --> x_i, y^(i) --> y_i, M^(k), Q^(k), Lambda^(k), .... By the rewriting, the number of notations is reduced by 50%! Some notations only appear a few times and can just be replaced by English, e.g., O_d (the set of orthogonal matrices).\n2. The significance of the paper relies on the significance of the exponential kernel exp(x-dot-y). If computing the transformer attention is the only application, the vanilla computational cost should be compared. Otherwise, I feel not certain if this is a widely used kernel family, and (e.g., for Gaussian kernel) in the application scenarios (e.g., regularized least squares, or support vector machines), whether computing KC is really the bottleneck.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to sincerely thank the Reviewer for all the comments. We address weaknesses and questions below.\n\n> Concerns about the writing.\n\nWe will address this concern by simplifying and de-densifying notation in the revision of the paper.\n\n> Otherwise, I feel not certain if this is a widely used kernel family, and … whether computing KC is really the bottleneck.\n\nNote that our derivations, as mentioned in the title and in e.g. lines 57-58, apply for both softmax and Gaussian kernels. Gaussian kernel is a very important object of research in the kernel methods literature and a lot of efforts have been made to speed it up in situations where the kernel matrix is big (e.g. in SVM or kernel regression where the number of data points is above a few thousands) – see [38,39,40] and thousands of follow ups of these papers.\n\n> Lines 76--77. Could you please provide more details on why the computational complexity of hat-K-times-C is O(LMn)?\n\nAs mentioned in the text, hat{K} = P S^T. Hence, hat{K} C = P S^T C = P (S^T C). The complexity of computing U = S^T C is O(LMn) since S^T, C are of shapes MxL and Lxn. The complexity of computing P U is also O(LMn) since P is of shape LxM and U is of shape Mxn.\n\n> In Eq. (8), the variable \"theta\" does not show up in the right-hand side, and is not elsewhere defined. Better if this variable can be defined.\n\nTheta is defined in line 126. We will expand its definition to make it clearer in the final revision.\n\n> It remains very unclear even after a careful study of Section 3, what are \"homogeneity heuristic\" and \"a certain optimization problem\" mentioned in Lines 124 and 125. \n\n“homogeneity heuristic” is defined in line 116 and the optimization problem is the minimization of (8) as mentioned in line 129. We will clarify these notations in the final version.\n\n> Can one achieve the same accuracy and computation complexity simply by using a larger M?\n\nAs we demonstrate in our experimental results (Figure 1), our RF variants result in up to e^10 variance improvements over the previous best variant. To achieve that by using a larger M means taking e^20 * M random features (since variance is proportional to M^{-1/2}) which would be prohibitively expensive.\n"}]

Author rebuttal:

Meta-review: The submission is about scaling up kernel methods using the random feature (RF) approach. Particularly, the authors concentrate on the problem of efficient multiplication/approximation with the Gram matrix associated to the softmax kernel (or more generally to the scaled softmax kernel defined in line 55), which is a commonly-appearing operation in various kernel algorithms. They propose the dense-exponential random feature class (DERF, line 145-146) extending the A,B,C scalar parameters of the family GERF (generalized exponential RF) family to matrices. After showing that these features are RFs for the softmax kernel (meant in the sense of Def. 2.1) and providing analytical expression for their variance (Theorem 4.1), they specialize the construction to ADERF (asymmetric dense exponential RFs), SDERF (symmetric dense-exponential RFs) and SADERF (simplified ADERF) where one can get the variance-minimizing parameters in closed form (Theorem 4.2, 4.3). The approach is illustrated in kernel regression and self-attention approximation.

Providing principled techniques to scale up kernel methods is an important problem in contemporary data science; the considered task is of clear interest to the NeurIPS community. The presented approach is novel, with relevance from both empirical and theoretical point of view, as it was elaborated by the reviewers.

Note: line 22: random features for operator-valued kernels: The authors might wish to cite "Romain Brault, Markus Heinonen, Florence d'Alché-Buc. Random Fourier Features For Operator-Valued Kernels, Asian Conference on Machine Learning (ACML), pages 110-125, 2016" as the first work in the area.
Decision: Accept (poster)