arXiv:1711.01596v1 [] 5 Nov 2017

Is Input Sparsity Time Possible for
Kernel Low-Rank Approximation?
Cameron Musco
MIT
cnmusco@mit.edu

David P. Woodruff
Carnegie Mellon University
dwoodruf@cs.cmu.edu

November 7, 2017
Abstract
Low-rank approximation is a common tool used to accelerate kernel methods: the n × n
kernel matrix K is approximated via a rank-k matrix K̃ which can be stored in much less space
and processed more quickly. In this work we study the limits of computationally efficient lowrank kernel approximation. We show that for a broad class of kernels, including the popular
Gaussian and polynomial kernels, computing a relative error k-rank approximation to K is
at least as difficult as multiplying the input data matrix A ∈ Rn×d by an arbitrary matrix
C ∈ Rd×k . Barring a breakthrough in fast matrix multiplication, when k is not too large, this
requires Ω(nnz(A)k) time where nnz(A) is the number of non-zeros in A. This lower bound
matches, in many parameter regimes, recent work on subquadratic time algorithms for lowrank approximation of general kernels [MM16, MW17], demonstrating that these algorithms are
unlikely to be significantly improved, in particular to O(nnz(A)) input sparsity runtimes. At
the same time there is hope: we show for the first time that O(nnz(A)) time approximation
is possible for general radial basis function kernels (e.g., the Gaussian kernel) for the closely
related problem of low-rank approximation of the kernelized dataset.

1

Introduction

The kernel method is a popular technique used to apply linear learning and classification algorithms
to datasets with nonlinear structure. Given training input points a1 , ..., an ∈ Rd , the idea is to
replace the standard Euclidean dot product hai , aj i = aTi aj with the kernel dot product ψ(ai , aj ),
where ψ : Rd × Rd → R+ is some positive semidefinite function. Popular kernel functions include
2
e.g., the Gaussian kernel with ψ(ai , aj ) = e−kai −aj k /σ for some bandwidth parameter σ and the
polynomial kernel of degree q with ψ(ai , aj ) = (c + aTi aj )q for some parameter c.
Throughout this work, we focus on kernels where ψ(ai , aj ) is a function of the dot products
aTi ai = kai k2 , aTj aj = kaj k2 , and aTi aj . Such functions encompass many kernels used in practice,
including the Gaussian kernel, the Laplace kernel, the polynomial kernel, and the Matern kernels.
Letting F be the reproducing kernel Hilbert space associated with ψ(·, ·), we can write ψ(ai , aj ) =
hφ(ai ), φ(aj )i where φ : Rd → F is a typically non-linear feature map. We let Φ = [φ(a1 ), ..., φ(an )]T
denote the kernelized dataset, whose ith row is the kernelized datapoint φ(ai ).
There is no requirement that Φ can be efficiently computed or stored – for example, in the case
of the Gaussian kernel, F is an infinite dimensional space. Thus, kernel methods typically work with
the kernel matrix K ∈ Rn×n with Ki,j = ψ(ai , aj ). We will also sometimes denote K = {ψ(ai , aj )}
1

to make it clear which kernel function it is generated by. We can equivalently write K = ΦΦT . As
long as all operations of an algorithm only access Φ via the dot products between its rows, they
can thus be implemented using just K without explicitly computing the feature map.
Unfortunately computing K is expensive, and a bottleneck for scaling kernel methods to large
datasets. For the kernels we consider, where ψ depends on dot products between the input points,
we must at least compute the Gram matrix AAT , requiring Θ(n2 d) time in general. Even if A
is sparse, this takes Θ(nnz(A)n) time. Storing K then takes Θ(n2 ) space, and processing it for
downstream applications like kernel ridge regression and kernel SVM can be even more expensive.

1.1

Low-rank kernel approximation

For this reason, a vast body of work studies how to efficiently approximate K via a low-rank
surrogate K̃ [SS00, AMS01, WS01, FS02, RR07, ANW14, LSS13, BJ02, DM05, ZTK08, BW09,
CKS11, WZ13, GM13]. If K̃ is rank-k, it can be stored in factored form in O(nk) space and
operated on quickly – e.g., it can be inverted in just O(nk 2 ) time to solve kernel ridge regression.
One possibility is to set K̃ = Kk where Kk is K’s best k-rank approximation – the projection
onto its top k eigenvectors.
over all rank-k K̃, the error kK − K̃kF , where kM kF is
P K2k minimizes,
1/2
the Frobenius norm: ( i,j Mi,j ) . It in fact minimizes error under any unitarily invariant norm,
e.g., the popular spectral norm. Unfortunately, Kk is prohibitively expensive to compute, requiring
Θ(n3 ) time in practice, or nω in theory using fast matrix multiplication, where ω ≈ 2.373 [LG14].
The idea of much prior work on low-rank kernel approximation is to find K̃ which is nearly
as good as Kk , but can be computed much more quickly. Specifically, it is natural to ask for K̃
fulfilling the following relative error guarantee for some parameter ǫ > 0:
kK − K̃kF ≤ (1 + ǫ)kK − Kk kF .

(1)

Other goals, such as nearly matching the spectral norm error kK − Kk k or approximating K entrywise have also been considered [RR07, GM13]. Of particular interest to our results is the closely
related goal of outputting an orthonormal basis Z ∈ Rn×k satisfying for any Φ with ΦΦT = K:
kΦ − ZZ T ΦkF ≤ (1 + ǫ)kΦ − Φk kF .

(2)

(2) can be viewed as a Kernel PCA guarantee – its asks us to find a low-rank subspace Z such
that the projection of our kernelized dataset Φ onto Z nearly optimally approximates this dataset.
Given Z, we can approximate K using K̃ = ZZ T ΦΦT ZZ T = ZZ T KZZ T . Alternatively, letting P
be the projection onto the row span of ZZ T Φ, we can write K̃ = ΦP ΦT , which can be computed
efficiently, for example, when P is a projection onto a subset of the kernelized datapoints [MM16].

1.2

Fast algorithms for relative-error kernel approximation

Until recently, all algorithms achieving the guarantees of (1) and (2) were at least as expensive as
computing the full matrix K, which was needed to compute the low-rank approximation [GM13].
However, recent work has shown that this is not required. Avron, Nguyen, and Woodruff
[ANW14] demonstrate that for the polynomial kernel, Z satisfying (2) can be computed in O(nnz(A)q)+
n poly(3q k/ǫ) time for a polynomial kernel with degree q.
Musco and Musco [MM16] give a fast algorithm for any kernel, using recursive Nyström sampling, which computes K̃ (in factored form) satisfying kK − K̃k ≤ λ, for input parameter λ. With
2

the proper setting of λ, it can output Z satisfying (2) (see Section C.3 of [MM16]). Computing Z
requires evaluating Õ(k/ǫ) columns of the kernel matrix along with Õ(n(k/ǫ)ω−1 ) additional time
for other computations. Assuming the kernel is a function of the dot products between the input
points, the kernel evaluations require Õ(nnz(A)k/ǫ) time. The results of [MM16] can also be used
√
to compute K̃ satisfying (1) with ǫ = n in Õ(nnz(A)k+nkω−1 ) time (see Appendix A of [MW17]).
Woodruff and Musco [MW17] show that for any kernel, and for any ǫ > 0, it√is possible top
achieve
2
(1) in Õ(nnz(A)k/ǫ)+n poly(k/ǫ) time plus the time needed to compute an Õ( nk/ǫ )× Õ( nk/ǫ)
submatrix of K. If A has uniform row sparsity – i.e., nnz(ai ) ≤ c nnz(A)/n for
√ some constant c and
all i, this step can be done in Õ(nnz(A)k/ǫ2.5 ) time. Alternatively, if d ≤ ( nk/ǫ2 )α for α < .314
this can be done in Õ(nk/ǫ4 ) = Õ(nnz(A)k/ǫ4 ) time using fast rectangular matrix multiplication
[LG12, GU17] (assuming that there are no all zero data points so n ≤ nnz(A).)

1.3

Our results

The algorithms of [MM16, MW17] make significant progress in efficiently solving (1) and (2) for general kernel matrices. They demonstrate that, surprisingly, a relative-error low-rank approximation
can be computed significantly faster than the time required to write down all of K.
A natural question is if these results can be improved. Even ignoring ǫ dependencies and
typically lower order terms, both algorithms use Ω(nnz(A)k) time. One might hope to improve this
to input sparsity, or near input sparsity time, Õ(nnz(A)), which is known for computing a low-rank
approximation of A itself [CW13]. The work of Avron et al. affirms that this is possible for the
kernel PCA guarantee of (2) for degree-q polynomial kernels, for constant q. Can this result be
extended to other popular kernels, or even more general classes?
1.3.1

Lower bounds

We show that achieving the guarantee of (1) significantly more efficiently than the work of [MM16,
MW17] is likely very difficult. Specifically, we prove that for a wide class of kernels, the kernel
low-rank approximation problem is as hard as multiplying the input A ∈ Rn×d by an arbitrary
C ∈ Rd×k . We have the following result for some common kernels to which our techniques apply:
Theorem 1 (Hardness for low-rank kernel approximation). Consider any polynomial kernel ψ(mi , mj ) =
2
(c + mTi mj )q , Gaussian kernel ψ(mi , mj ) = e−kmi −mj k /σ , or the linear kernel ψ(mi , mj ) =
mTi mj . Assume there is an algorithm which given M ∈ Rn×d with associated kernel matrix
K = {ψ(mi , mj )}, returns N ∈ Rn×k in o(nnz(M )k) time satisfying:
kK − N N T k2F ≤ ∆kK − Kk k2F
for some approximation factor ∆. Then there is an o(nnz(A)k) + O(nk2 ) time algorithm for multiplying arbitrary integer matrices A ∈ Rn×d , C ∈ Rd×k .
The above applies for any approximation factor ∆. While we work in the real RAM model,
ignoring bit complexity, as long as ∆ = poly(n) and A, C have polynomially bounded entries, our
reduction from multiplication to low-rank approximation is achieved using matrices that can be
represented with just O(log(n + d)) bits per entry.
√
Theorem 1 shows that the runtime of Õ(nnz(A)k + nk ω−1 ) for ∆ = n achieved by [MM16] for
general kernels cannot be significantly improved without advancing the state-of-the-art in matrix
3

multiplication. Currently no general algorithm is known for multiplying integer A ∈ Rn×d , C ∈
Rd×k in o(nnz(A)k) time, except when k ≥ nα for α < .314 and A is dense. In this case, AC can
be computed in O(nd) time using fast rectangular matrix multiplication
[LG12, GU17].
√
As discussed, when A has uniform row sparsity or when d ≤ ( nk/ǫ2 )α , the runtime of [MW17]
for ∆ = (1 + ǫ), ignoring ǫ dependencies and typically lower order terms, is Õ(nnz(A)k), which is
also nearly tight.
In recent work, Backurs et al. [BIS17] give lower bounds for a number of kernel learning
problems, including kernel PCA for the Gaussian kernel. However, their strong bound, of Ω(n2 )
time, requires very small error ∆ = exp(−ω(log2 n), whereas ours applies for any relative error ∆.
1.3.2

Improved algorithm for radial basis function kernels

In contrast to the above negative result, we demonstrate that achieving the alternative Kernel PCA
guarantee of (2) is possible in input sparsity time for any shift and rotationally invariant kernel
– e.g., any radial basis function kernel where ψ(xi , xj ) = f (kxi − xj k). This result significantly
extends the progress of Avron et al. [ANW14] on the polynomial kernel.
Our algorithm is based off a fast implementation of the random Fourier features method [RR07],
which uses the fact that that the Fourier transform of any shift invariant kernel is a probability
distribution after appropriate scaling (a consequence of Bochner’s theorem). Sampling frequencies
from this distribution gives an approximation to ψ(·, ·) and consequentially the matrix K.

We employ a new analysis of this method [AKM+ 17], which shows that sampling Õ ǫ2nλ random
Fourier features suffices to give K̃ = Φ̃Φ̃T satisfying the spectral approximation guarantee:
(1 − ǫ)(K̃ + λI)  K + λI  (1 + ǫ)(K̃ + λI).
If we set λ ≤ σk+1 (K)/k, we can show that Φ̃ also gives a projection-cost preserving sketch
[CEM+ 15] for the kernelized dataset Φ. This ensures that any Z satisfying kΦ̃ − ZZ T Φ̃k2F ≤
(1 + ǫ)kΦ̃ − Φ̃k k2F also satisfies kΦ − ZZ T Φk2F ≤ (1 + O(ǫ))kΦ
− Φk k2F and thus achieves (2).


Our algorithm samples s = Õ ǫ2nλ = Õ ǫ2 σ nk (K) random Fourier features, which naively rek+1
quires O(nnz(A)s) time. We show that this can be accelerated to O(nnz(A))+poly(n, s) time, using
a recent result of Kapralov et al. on fast multiplication by random Gaussian matrices [KPW16].
Our technique is analogous to the ‘Fastfood’ approach to accelerating random Fourier features
using fast Hadamard transforms [LSS13]. However, our runtime scales with nnz(A), which can
be significantly smaller than the Õ(nd) runtime given by Fastfood when A is sparse. Our main
algorithmic result is:
Theorem 2 (Input sparsity time kernel PCA). There is an algorithm that given A ∈ Rn×d along
with shift and rotation-invariant kernel function ψ : Rd × Rd → R+ with ψ(x, x) = 1, outputs, with
probability 99/100, Z ∈ Rn×k satisfying:
kΦ − ZZ T Φk2F ≤ (1 + ǫ)kΦ − Φk k2F
for any Φ with ΦΦT = K = {ψ(ai , aj )} and any ǫ > 0. Letting σk+1 denote the (k + 1)th largest
eigenvalue of Kand ω < 2.373 be the exponent
of fast matrix multiplication, the algorithm runs in
ω−1.5 

k
time.
O(nnz(A)) + Õ nω+1.5 · σk+1
ǫ2
4

We note that the runtime of our algorithm is O(nnz(A)) whenever n, k, 1/σk+1 , and 1/ǫ are
not too large. Due to the relatively poor dependence on n, the algorithm is relevant for very
high dimensional datasets with d ≫ n. Such datasets are found often, e.g., in genetics applications
[HDC+ 01, JDMP11]. While we have dependence on 1/σk+1 , in the natural setting, we only compute
a low-rank approximation up to an error threshold, ignoring very small eigenvalues of K, and so
σk+1 will not be too small. We do note that if we apply Theorem 2 to the low-rank approximation
instances given by our lower bound construction, σk+1 can be very small, ≤ 1/ poly(n, d) for
matrices with poly(n) bounded entries. Thus, removing this dependence is an important open
question in understanding the complexity of low-rank kernel approximation.
We leave open the possibility of improving our algorithm, achieving O(nnz(A)) + n · poly(k, ǫ)
runtime, which would match the state-of-the-art for low-rank approximation of non-kernelized
matrices [CW13]. Alternatively, it is possible that a lower bound can be shown, proving the that
high n dependence, or the 1/σk+1 term are required even for the Kernel PCA guarantee of (2).

2

Lower bounds

Our lower bound proof argues that for a broad class of kernels, given input M , a low-rank approximation of the associated kernel matrix K achieving (1) can be used to obtain a close approximation
to the Gram matrix M M T . We write ψ(mi , mj ) as a function of mTi mj (or kmi − mj k2 for distance kernels) and expand this function as a power series. We show that the if input points are
appropriately rescaled, the contribution of degree-1 term mTi mj dominates, and hence our kernel
matrix approximates M M T , up to some easy to compute low-rank components.
We then show that such an approximation can be used to give a fast algorithm for multiplying
any two integer matrices A ∈ Rn×d and C ∈ Rd×k . The key idea is to set M = [A, wC] where w is
a large weight. We then have:


AAT
wAC
T
MM =
.
wC T AT w2 C T C
Since w is very large, the AAT block is relatively very small, and so M M T is nearly rank-2k – it
has a ‘heavy’ strip of elements in its last k rows and columns. Thus, computing a relative-error
rank-2k approximation to M M T recovers all entries except those in the AAT block very accurately,
and importantly, recovers the wAC block and so the product AC.

2.1

Lower bound for low-rank approximation of MM T .

We first illustrate our lower bound technique by showing hardness of direct approximation of M M T .
Theorem 3 (Hardness of low-rank approximation for M M T ). Assume there is an algorithm A
which given any M ∈ Rn×d returns N ∈ Rn×k such that kM M T − N N T k2F ≤ ∆1 kM M T −
(M M T )k k2F in T (M, k) time for some approximation factor ∆1 .
For any √
A ∈ Rn×d and C ∈ Rd×k each with integer entries in [−∆2 , ∆2 ], let B = [AT , wC]T
where w = 3 ∆1 ∆22 nd. It is possible to compute the product AC in time T (B, 2k) + O(nk ω−1 ).
Proof. We can write the (n + k) × (n + k) matrix BB T as:

AAT
T
T
T
BB = [A , wC] [A, wC] =
wC T AT
5


wAC
.
w2 C T C

Let Q ∈ Rn×2k be an orthogonal span for the columns of the n × 2k matrix:


0
wAC
V w2 C T C
where V ∈ Rk×k spans the columns of wC T AT ∈ Rk×n . The projection QQT BB T gives the best
Frobenius norm approximation to BB T in the span of Q. We can see that:
kBB T − (BB T )2k k2F ≤ kBB T − QQT BB T k2F ≤



AAT
0

0
0



2
F

≤ ∆42 n2 d2

(3)

since each entry of A is bounded in magnitude by ∆2 and so each entry of AAT is bounded by d∆22 .
Let N be the matrix returned by running A on B with rank 2k. In order to achieve the
approximation bound of kBB T − N N T k2F ≤ ∆1 kBB T − (BB T )2k k2F we must have, for all i, j:
(BB T − N N T )2i,j ≤ kBB T − N N T k2F ≤ ∆1 ∆42 n2 d2

√
where the last inequality is from (3). This gives |BB T − N N T |i,j ≤ ∆1 ∆22 nd. Since A and
√ C have
T is an integer multiple of w = 3 ∆ ∆2 nd.
integer entries, each entry in the submatrix wAC of
BB
1 2
√
Since (N N T )i,j approximates this entry to error ∆1 ∆22 nd, by simply rounding (N N T )i,j to the
nearest multiple of w, we obtain the entry exactly. Thus, given N , we can exactly recover AC in
O(nkω−1 ) time by computing the n × k submatrix corresponding to AC in BB T .
Theorem 3 gives our main bound Theorem 1 for the case of the linear kernel ψ(mi , mj ) = mTi mj .
Proof of Theorem 1 – Linear Kernel. We apply Theorem 3 after noting that for B = [AT , wC]T ,
nnz(B) ≤ nnz(A) + nk and so T (B, 2k) = o(nnz(A)k) + O(nk 2 ).
We show in Appendix A that there is an algorithm which nearly matches the lower bound of
Theorem 1 for any ∆ = (1 + ǫ) for any ǫ > 0. Further, in Appendix B we show that even just
outputting an orthogonal matrix Z ∈ Rn×k such that K̃ = ZZ T M M T is a relative-error low-rank
approximation of M M T , but not computing a factorization of K̃ itself, is enough to give fast
multiplication of integer matrices A and C.

2.2

Lower bound for dot product kernels

We now extend Theorem 3 to general dot product kernels – where ψ(ai , aj ) = f (aTi aj ) for some
function f . This includes, for example, the polynomial kernel.
Theorem 4 (Hardness of low-rank approximation for dot product kernels). Consider any kernel
ψ : Rd P
× Rd → R+ with ψ(ai , aj ) = f (aTi aj ) for some function f which can be expanded as
q
q−1 and for all q ≥ 2 and some G ≥ 1.
f (x) = ∞
q=0 cq x with c1 6= 0 and |cq /c1 | ≤ G
Assume there is an algorithm A which given M ∈ Rn×d with kernel matrix K = {ψ(mi , mj )},
returns N ∈ Rn×k satisfying kK − N N T k2F ≤ ∆1 kK − Kk k in T (M, k) time.
For any A ∈ Rn×d , C ∈ Rd×k with integer entries in [−∆2 , ∆2 ], let B = [w1 AT , w2 C]T with
1
w1 = 12√∆w2∆2 nd , w2 = 4√Gd∆
. Then it is possible to compute AC in time T (B, 2k+1)+O(nk ω−1 ).
1

2

2

6

Proof. Using our decomposition of ψ(·, ·), we can write the kernel matrix for B and ψ as:




1 1
w1 w2 AC
w12 AAT
K = c0
+ c1
+ c2 K (2) + c3 K (3) + ...
1 1
w1 w2 C T AT w22 C T C

(4)

(q)

where Ki,j = (bTi bj )q and 1 denotes the all ones matrix of appropriate size. The key idea is to
show that the contribution of the K (q) terms is small, and so any relative-error rank-(2k + 1)
approximation to K must recover an approximation to BB T , and thus the product AC as in
Theorem 3.
1
, the fact that w1 < w2 , and our bound on the entries of A and
By our setting of w2 = 4√Gd∆
2

C, we have for all i, j, |bTi bj | ≤ w22 d∆22 <
∞
X
q=2

(q)
cq Ki,j

≤

c1 |bTi bj |


·

∞
X

1
16G .

q−1

G

q=2

Thus, for any i, j, using that |cq /c1 | ≤ Gq−1 :

|bTi bj |q−1

≤

c1 |bTi bj |

∞
X
q=2

Gq−1
1
≤ c1 |bTi bj |.
q−1
(16G)
12

(5)




1 1
, with its top right n × n block set to 0. K̄ just has its
1 1
last k columns and rows non-zero, so has rank ≤ 2k. Let Q ∈ Rn×2k+1 be an orthogonal span for
the columns K̄ along with the all ones vector of length n. Let N be the result of running A on B
with rank 2k + 1. Then we have:
Let K̄ be the matrix

K − c0

kK − N N T k2F ≤ ∆1 kK − K2k+1 k2F ≤ ∆1 kK − QQT Kk2F


(c1 w12 AAT + c2 K̂ (2) + ...) 0
≤ ∆1
0
0

2

(6)
F

where K̂ (q) denotes the top left n × n submatrix of K (q) . By our bound on the entries of A and (5):



13
c1 w12 AAT i,j ≤ 2c1 w12 d∆22 .
c1 w12 AAT + c2 K̂ (2) + c3 K̂ (3) + ...
≤
12
i,j
Plugging back into (6) and using w1 =

√ w2
,
12 ∆1 ∆22 nd

this gives for any i, j:

(K − N N T )i,j ≤ kK − N N T kF ≤

p
√

∆1 n2 · 2c1 w12 d∆22

∆1 n · 2c1 d∆22
√
· w1 w2
12 ∆1 ∆22 nd
w1 w2 c1
.
≤
6

≤

(7)

Since A and C have integer entries, each entry of c1 w1 w2 AC is an integer multiple of c1 w1 w2 . By
the decomposition of (4) and the bound of (5), if we subtract c0 from the corresponding entry of
K and round it to the nearest multiple of c1 w1 w2 , we will recover the entry of AC. By the bound
of (7), we can likewise round the corresponding entry of N N T . Computing all nk of these entries
given N takes time O(nk ω−1 ), giving the theorem.
Theorem 4 lets us lower bound the time to compute a low-rank kernel approximation for any
kernel function expressible as a reasonable power expansion of aTi aj . As a straightforward example,
it gives the lower bound for the polynomial kernel of any degree stated in Theorem 1.
7

Proof of Theorem 1 – Polynomial Kernel. We applyPTheorem 4, noting that ψ(mi , mj ) = (c +
mTi mj )q can be written as f (mTi mj ) where f (x) = qj=0 cj xj with cj = cq−j qj . Thus c1 6= 0 and
|cj /c1 | ≤ Gj−1 for G = (q/c). Finally note that nnz(B) ≤ nnz(A) + nk giving the result.

2.3

Lower bound for distance kernels

We finally extend Theorem 4 to handle kernels like the Gaussian kernel whose value depends on
the squared distance kai − aj k2 rather than just the dot product aTi aj . We prove:
Theorem 5 (Hardness of low-rank approximation for distance kernels). Consider any kernel function ψ : RdP
× Rd → R+ with ψ(ai , aj ) = f (kai − aj k2 ) for some function f which can be expanded
q
q−1 and for all q ≥ 2 and some G ≥ 1.
as f (x) = ∞
q=0 cq x with c1 6= 0 and |cq /c1 | ≤ G
Assume there is an algorithm A which given input M ∈ Rn×d with kernel matrix K = {ψ(mi , mj )},
returns N ∈ Rn×k satisfying kK − N N T k2F ≤ ∆1 kK − Kk k in T (M, k) time.
For any A ∈ Rn×d , C ∈ Rd×k with integer entries in [−∆2 , ∆2 ], let B = [w1 AT , w2 C]T with
1√
w1 = 36√∆w2∆2 nd , w2 = (16Gd2 ∆4 )(36
. It is possible to compute AC in T (B, 2k + 3) +
∆ ∆2 nd)
1

2

O(nkω−1 ) time.

2

1

2

The proof of Theorem 5 is similar to that of Theorem 4, and relegated to Appendix C. The
key idea is to write K as a polynomial in the distance matrix D with Di,j = kbi − bj k22 . Since
kbi − bj k22 = kbi k22 + kbj k22 − 2bTi bj , D can be written as −2BB T plus a rank-2 component. By
setting w1 , w2 sufficiently small, as in the proof of Theorem 4, we ensure that the higher powers of
D are negligible, and thus that our low-rank approximation must accurately recover the submatrix
of BB T corresponding to AC. Theorem 5 gives Theorem 1 for the popular Gaussian kernel:
Proof of Theorem 1 – Gaussian Kernel. ψ(mi , mj ) can be written as f (kmi − mj k2 ) where f (x) =
P
(−1/σ)q q
x . Thus c1 6= 0 and |cq /c1 | ≤ Gq−1 for G = 1/σ. Applying Theorem 5 and
e−x/σ = ∞
q=0
q!
bounding nnz(B) ≤ nnz(A) + nk, gives the result.

3

Input sparsity time kernel PCA for radial basis kernels

Theorem 1 gives little hope for achieving o(nnz(A)k) time for low-rank kernel approximation.
However, the guarantee of (1) is not the only way of measuring the quality of K̃. Here we show
that for shift/rotationally invariant kernels, including e.g., radial basis kernels, input sparsity time
can be achieved for the kernel PCA goal of (2).

3.1

Basic algorithm

Our technique is based on the random Fourier features technique [RR07]. Given any shift-invariant
kernel, ψ(x, y) = ψ(x − y) with ψ(0) = 1 (we will assume this w.l.o.g. as the function can always
be scaled), there is a probability density function p(η) over vectors in Rd such that:
Z
T
e−2πiη (x−y) p(η)dη.
(8)
ψ(x − y) =
Rd

p(η) is just the (inverse) Fourier transform of ψ(·), and is a density function by Bochner’s theorem.
Informally, given A ∈ Rn×d if we let Z denote the matrix with columns z(η) indexed by η ∈ Rd .
8

T

z(η)j = e−2πiη aj . Then (8) gives ZP Z ∗ = K where P is diagonal with Pη,η = p(η), and Z ∗ denotes
the Hermitian transpose.
The idea of random Fourier features is to select s frequencies η1 , ..., ηs according to the density
p(η) and set Z̃ = √1s [z(η1 ), ...z(ηs )]. K̃ = Z̃ Z̃ T is then used to approximate K.
In recent work, Avron et al. [AKM+ 17] give a new analysis of random Fourier features. Extending prior work on ridge leverage scores in the discrete setting [AM15, CMM17], they define the
ridge leverage function for parameter λ > 0:
τλ (η) = p(η)z(η)∗ (K + λI)−1 z(η)

(9)

As part of their results, which seek K̃ that spectrally approximates K, they prove the following:
Lemma 6. For all η, τλ (η) ≤ n/λ.
While simple, this bound is key to our algorithm. It was shown in [CMM17] that if the columns
of a matrix are sampled by over-approximations to their ridge leverage scores (with appropriately
set λ), the sample is a projection-cost preserving sketch for the original matrix. That is, it can be
used as a surrogate in computing a low-rank approximation. The results of [CMM17] carry over to
the continuous setting giving, in conjunction with Lemma 6:
Lemma 7 (Projection-cost preserving sketch via random Fourier features). Consider any A ∈ Rn×d
and shift-invariant kernel ψ(·) with ψ(0) = 1, with associated kernel matrix K = {ψ(ai − aj )} and
P
for sufficiently
kernel Fourier transform p(η). For any 0 < λ ≤ k1 ni=k+1 σi (K), let s = cn log(n/δλ)
ǫ2 λ
1
large c and let Z̃ = √s [z(η1 ), ..., z(ηs )] where η1 , ..., ηs are sampled independently according to p(η).
Then with probability ≥ 1 − δ, for any orthonormal Q ∈ Rn×k and any Φ with ΦΦT = K:
(1 − ǫ)kQQT Z̃ − Z̃k2F ≤ kQQT Φ − Φk2F ≤ (1 + ǫ)kQQT Z̃ − Z̃k2F .

(10)

By (10) if we compute Q satisfying kQQT Z̃ − Z̃k2F ≤ (1 + ǫ)kZ̃ − Z̃k k2F then we have:
(1 + ǫ)2
kUk UkT Φ − Φk2F
1−ǫ
= (1 + O(ǫ))kΦ − Φk k2F

kQQT Φ − Φk2F ≤ (1 + ǫ)2 kZ̃ − Z̃k k2F ≤

where Uk ∈ Rn×k contains the top k column singular vectors of Φ. By adjusting constants on ǫ by
making c large enough, we thus have the relative error low-rank approximation guarantee of (2).
It remains to show that this approach can be implemented efficiently.

3.2

Input sparsity time implementation

Given Z̃ sampled as in Lemma 7, we can find a near optimal subspace Q using any input sparsity
time low-rank approximation algorithm (e.g., [CW13, NN13]). We have the following Corollary:
Corollary 8. Given Z̃ sampled as in Lemma 7 with s = Θ̃( ǫ2 σ nk (K) ), there is an algorithm running
k+1

in time

2
Õ( ǫ2 σ n k(K) )
k+1

that computes Q satisfying with high probability, for any Φ with ΦΦT = K:
kQQT Φ − Φk2F ≤ (1 + ǫ)kΦ − Φk k2F .

With Corollary 8 in place the main bottleneck to our approach becomes computing Z̃.
9

3.2.1

Sampling Frequencies

To compute Z̃, we first sample η1 , ..., ηs according to p(η). Here we use the rotational invariance of
ψ(·). In this case, p(η) is also rotationally invariant [LSS13] and so, letting p̂(·) be the distribution
over norms of vectors sampled from p(η) we can sample η1 , ..., ηn by first selecting s random Gaussian vectors and then rescaling them to have norms distributed according to p̂(·). That is, we can
write [η1 , ..., ηn ] = GD where G ∈ Rd×s is a random Gaussian matrix and D is a diagonal rescaling
m
with m ∼ p̂. We will assume that p̂ can be sampled from in O(1) time.
matrix with Dii = kG
ik
This is true for many natural kernels – e.g., for the Gaussian kernel, p̂ is just a Gaussian density.
3.2.2

Computing Z̃

Due to our large sample size, s > n, even writing down G above requires Ω(nd) time. However,
to form Z̃ we do not need G itself: it suffices to compute for m = 1, ..., s the column z(ηm ) with
T
z(ηm )j = e−2πiηm aj . This requires computing AGD, which contains the appropriate dot products
aTj ηm for all m, j. We use a recent result [KPW16] which shows that this can be performed
approximately in input sparsity time:
4

3 ω−1.5

Lemma 9 (From Theorem 1 of [KPW16]). There is an algorithm running in O(nnz(A)+ log dnδ s
)
time which outputs random B whose distribution has total variation distance at most δ from the
distribution of AG where G ∈ Rd×s is a random Gaussian matrix. Here, ω < 2.373 is the exponent
of fast matrix multiplication.
Proof. Theorem 1 of [KPW16] shows that for B to have total variation distance δ from the distribution of AG it suffices to set B = ACG′ where C is a d × O(log4 dn2 s1/2 /δ) CountSketch matrix
and G′ is an O(log4 dn2 s1/2 /δ) × s random Gaussian matrix. Computing AC requires O(nnz(A))
4
3 s1.5
) time if fast matrix multiplication is
time. Multiplying the result by G′ then requires O( log dn
δ
not employed. Using fast matrix multiplication, this can be improved to O( log

4

dn3 sω−1.5
).
δ

Applying Lemma 9 with δ = 1/200 lets us compute random BD with total variation distance
1/200 from AGD. Thus, the distribution of Z̃ generated from this matrix has total variation
distance ≤ 1/200 from the Z̃ generated from the true random Fourier features distribution. So,
by Corollary 8, we can use Z̃ to compute Q satisfying kQQT Φ − Φk2F ≤ (1 + ǫ)kΦ − Φk k2F with
probability 1/100 accounting for the the total variation difference and the failure probability of
Corollary 8. This yields our main algorithmic result, Theorem 2.

3.3

An alternative approach

We conclude by noting that near input sparsity time Kernel PCA can also be achieved for a broad
class of kernels using a very different approach. We can approximate ψ(·, ·) via an expansion into
polynomial kernel matrices as is done in [CKS11] and then apply the sketching algorithms for the
polynomial kernel developed in [ANW14]. As long as the expansion achieves high accuracy with
low degree, and as long as 1/σk+1 is not too small – since this will control the necessary approximation factor, this technique can yield runtimes of the form Õ(nnz(A)) + poly(n, k, 1/σk+1 , 1/ǫ),
giving improved dependence on n for some kernels over our random Fourier features method. Improving the poly(n, k, 1/σk+1 , 1/ǫ) term in both these methods, and especially removing the 1/σk+1
dependence and achieving linear dependence on n is an interesting open question for future work.
10

4

Conclusion

In this work we have shown that for a broad class of kernels, including the Gaussian, polynomial, and
linear kernels, given data matrix A, computing a relative error low-rank approximation to A’s kernel
matrix K (i.e., satisfying (1)) requires at least Ω(nnz(A)k) time, barring a major breakthrough in
the runtime of matrix multiplication. In the constant error regime, this lower bound essentially
matches the runtimes given by recent work on subquadratic time kernel and PSD matrix low-rank
approximation [MM16, MW17].
We show that for the alternative kernel PCA guarantee of (2), a potentially faster runtime of
O(nnz(A)) + poly(n, k, 1/σk+1 , 1/ǫ) can be achieved for general shift and rotation-invariant kernels.
Practically, improving the second term in our runtime, especially the poor dependence on n, is
an important open question. Generally, computing the kernel matrix K explicitly requires O(n2 d)
time, and so our algorithm only gives runtime gains when d is large compared to n – at least
Ω(nω−.5 ), even ignoring k, σk+1 , and ǫ dependencies. Theoretically, removing the dependence on
σk+1 would be of interest, as it would give input sparsity runtime without any assumptions on
the matrix A (i.e., that σk+1 is not too small). Resolving this question has strong connections to
finding efficient kernel subspace embeddings, which approximate the full spectrum of K.

References
[AKM+ 17] 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 Proceedings of the 34th International
Conference on Machine Learning (ICML), 2017.
[AM15]

Ahmed Alaoui and Michael W Mahoney. Fast randomized kernel ridge regression
with statistical guarantees. In Advances in Neural Information Processing Systems
28 (NIPS), pages 775–783, 2015.

[AMS01]

Dimitris Achlioptas, Frank Mcsherry, and Bernhard Schölkopf. Sampling techniques
for kernel methods. In Advances in Neural Information Processing Systems 14 (NIPS),
2001.

[ANW14]

Haim Avron, Huy Nguyen, and David Woodruff. Subspace embeddings for the polynomial kernel. In Advances in Neural Information Processing Systems 27 (NIPS), pages
2258–2266, 2014.

[BIS17]

Arturs Backurs, Piotr Indyk, and Ludwig Schmidt. On the fine-grained complexity
of empirical risk minimization: Kernel methods and neural networks. In Advances in
Neural Information Processing Systems 30 (NIPS), 2017.

[BJ02]

Francis Bach and Michael I. Jordan. Kernel independent component analysis. Journal
of Machine Learning Research, 3(Jul):1–48, 2002.

[BW09]

Mohamed-Ali Belabbas and Patrick J. Wolfe. Spectral methods in machine learning:
New strategies for very large datasets. Proceedings of the National Academy of Sciences
of the USA, 106:369–374, 2009.

11

[CEM+ 15] Michael B. Cohen, Sam Elder, Cameron Musco, Christopher Musco, and Madalina
Persu. Dimensionality reduction for k-means clustering and low rank approximation.
In Proceedings of the 47th Annual ACM Symposium on Theory of Computing (STOC),
pages 163–172, 2015.
[CKS11]

Andrew Cotter, Joseph Keshet, and Nathan Srebro. Explicit approximations of the
Gaussian kernel. arXiv:1109.4603, 2011.

[CMM17]

Michael B. Cohen, Cameron Musco, and Christopher Musco. Input sparsity time lowrank approximation via ridge leverage score sampling. In Proceedings of the 28th Annual
ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 1758–1777, 2017.

[CW13]

Kenneth L Clarkson and David P Woodruff. Low rank approximation and regression
in input sparsity time. In Proceedings of the 45th Annual ACM Symposium on Theory
of Computing (STOC), pages 81–90, 2013.

[CW17]

Kenneth L. Clarkson and David P. Woodruff. Low-rank PSD approximation in inputsparsity time. In Proceedings of the 28th Annual ACM-SIAM Symposium on Discrete
Algorithms (SODA), pages 2061–2072, 2017.

[DM05]

Petros Drineas and Michael W Mahoney. On the Nyström method for approximating a Gram matrix for improved kernel-based learning. Journal of Machine Learning
Research, 6:2153–2175, 2005.

[FS02]

Shai Fine and Katya Scheinberg. Efficient SVM training using low-rank kernel representations. Journal of Machine Learning Research, 2:243–264, 2002.

[FT07]

Shmuel Friedland and Anatoli Torokhti. Generalized rank-constrained matrix approximations. SIAM Journal on Matrix Analysis and Applications, 29(2):656–659, 2007.

[GM13]

Alex Gittens and Michael Mahoney. Revisiting the Nyström method for improved
large-scale machine learning. In Proceedings of the 30th International Conference on
Machine Learning (ICML), pages 567–575, 2013. Full version at arXiv:1303.1849.

[GU17]

François Le Gall and Florent Urrutia. Improved rectangular matrix multiplication using
powers of the Coppersmith-Winograd tensor. arXiv:1708.05622, 2017.

[HDC+ 01] Ingrid Hedenfalk, David Duggan, Yidong Chen, Michael Radmacher, Michael Bittner,
Richard Simon, Paul Meltzer, Barry Gusterson, Manel Esteller, Mark Raffeld, et al.
Gene-expression profiles in hereditary breast cancer. New England Journal of Medicine,
344(8):539–548, 2001.
[JDMP11] Asif Javed, Petros Drineas, Michael W Mahoney, and Peristera Paschou. Efficient
genomewide selection of PCA-correlated tSNPs for genotype imputation. Annals of
Human Genetics, 75(6):707–722, 2011.
[KPW16]

Michael Kapralov, Vamsi Potluru, and David Woodruff. How to fake multiply by a
Gaussian matrix. In Proceedings of the 33rd International Conference on Machine
Learning (ICML), pages 2101–2110, 2016.

12

[LG12]

François Le Gall. Faster algorithms for rectangular matrix multiplication. In Proceedings of the 53rd Annual IEEE Symposium on Foundations of Computer Science
(FOCS), pages 514–523, 2012.

[LG14]

François Le Gall. Powers of tensors and fast matrix multiplication. In Proceedings of the
39th international symposium on symbolic and algebraic computation, pages 296–303.
ACM, 2014.

[LSS13]

Quoc Le, Tamás Sarlós, and Alexander Smola. Fastfood - Computing Hilbert space
expansions in loglinear time. In Proceedings of the 30th International Conference on
Machine Learning (ICML), pages 244–252, 2013.

[MM16]

Cameron Musco and Christopher Musco. Recursive sampling for the Nyström method.
In Advances in Neural Information Processing Systems 30 (NIPS), 2016.

[MW17]

Cameron Musco and David P Woodruff. Sublinear time low-rank approximation of
positive semidefinite matrices. In Proceedings of the 58th Annual IEEE Symposium on
Foundations of Computer Science (FOCS), 2017.

[NN13]

Jelani Nelson and Huy L Nguyên. OSNAP: Faster numerical linear algebra algorithms
via sparser subspace embeddings. In Proceedings of the 54th Annual IEEE Symposium
on Foundations of Computer Science (FOCS), pages 117–126, 2013.

[RR07]

Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In
Advances in Neural Information Processing Systems 20 (NIPS), pages 1177–1184, 2007.

[SS00]

Alex J Smola and Bernhard Schökopf. Sparse greedy matrix approximation for machine
learning. In Proceedings of the 17th International Conference on Machine Learning
(ICML), pages 911–918, 2000.

[WS01]

Christopher Williams and Matthias Seeger. Using the Nyström method to speed up
kernel machines. In Advances in Neural Information Processing Systems 14 (NIPS),
pages 682–688, 2001.

[WZ13]

Shusen Wang and Zhihua Zhang. Improving CUR matrix decomposition and the
Nyström approximation via adaptive sampling. Journal of Machine Learning Research,
14:2729–2769, 2013.

[ZTK08]

Kai Zhang, Ivor W. Tsang, and James T. Kwok. Improved Nyström low-rank approximation and error analysis. In Proceedings of the 25th International Conference on
Machine Learning (ICML), pages 1232–1239, 2008.

13

A

Fast low-rank approximation of AAT

We give an algorithm which matches the lower bound of Theorem 3.
Theorem 10. There is an algorithm, which given A ∈ Rn×d computes N ∈ Rn×k in O(nnz(A)k) +
n · poly(k/ǫ) time such that probability 99/100:
kAAT − N N T k2F ≤ (1 + ǫ)kAAT − (AAT )k k2F .
Proof. It is known (see Lemma 11 of [CW17]) that there exists a distribution over random matrices
R, S ∈ Rn×O(k/ǫ) which can be applied to A in O(nnz(A)) + n · poly(k/ǫ) time such that with
probability 199/200, setting
Y∗ =

arg min
Y ∈O(k/ǫ)×O(k/ǫ) with rank k

kAAT RY S T AAT − AAT k2F

we have:
kAAT RY ∗ S T AAT − AAT k2F ≤ (1 + ǫ)kAAT − (AAT )k k2F .
We can solve for an approximately optimal Ỹ by further sketching our problem on the left and right
(similar to the technique used in Lemma 15 of [CW17]). Specifically, if we let TL , TR ∈ Rn×poly(k/ǫ)
be drawn from the Count Sketch distribution, we can solve:
Ỹ =

arg min
Y ∈O(k/ǫ)×O(k/ǫ) with rank k

kTLT AAT RY S T AAT TR − TLT AAT TR k2F

and are guaranteed that with probability 99/100,
kAAT RỸ S T AAT − AAT k2F ≤ (1 + 2ǫ)kAAT − (AAT )k k2F .

(11)

Computing Ỹ requires forming TLT A, AT R, S T A, and AT TR and then multiplying the appropriate matrices together. This takes O(nnz(A)) + n poly(k/ǫ) time. Once TLT AAT R, S T AAT TR
and TLT AAT TR have been formed we can solve for Ỹ in poly(k/ǫ) time using the formula of [FT07].
Finally, since Ỹ is rank-k we can factor Ỹ = V V T for V ∈ RO(k/ǫ)×k using the SVD. We can
then compute N1 = AAT RV ∈ Rn×k and N3 = AAT SV ∈ Rn×k which satisfy kAAT − N1 N2T k2F ≤
(1 + 2ǫ)kAAT − (AAT )k k2F with probability 99/100 by (11).
N1 and N2 both require O(nnz(A)k) + n · poly(k/ǫ) time to compute. The theorem follows
from adjusting constants on ǫ and noting that we can symmetrize N1 N2T to form N N T if desired
in n · poly(k/ǫ) time.

B

Hardness of outputting a low-rank subspace

Theorem 3 shows a lower bound on outputting a relative-error low-rank approximation to M M T .
Here we show that this hardness extends to the possibly easier problem of just outputting a lowrank span that contains a relative-error low-rank approximation. This result extends analogously
to the other kernel lower bounds discussed in Section 2.
14

Theorem 11 (Hardness of low-rank span for M M T ). Assume there is an algorithm A which given
any M ∈ Rn×d returns orthonormal Z ∈ Rn×k such that kM M T − ZZ T M M T k2F ≤ ∆1 kM M T −
(M M T )k k2F in T (M, k) time for some approximation factor ∆1 .
For any √
A ∈ Rn×d and C ∈ Rd×k each with integer entries in [−∆2 , ∆2 ], let B = [AT , wC]T
where w = 3 ∆1 ∆22 nd. It is possible to compute the product AC in time T (B, 2k)+ Õ((n+d)kω−1 ).
Proof. ZZ T M M T is the projection of M onto the column span of Z. This projection can be
performed approximately using standard leverage score sampling techniques (see e.g., [CW13]).
Let S ∈ Rs×n be a sampling matrix sampling rows of Z by its row norms (its leverage scores since
it is orthonormal) where s = c(k log k) or some sufficiently large constant c. Let R be the n × k
matrix which selects the last k columns of M M T .
Letting X ∗ = arg minX∈k×k kZX T − M M T Rk2F and X = arg minX∈k×k kSZX T − SM M T Rk2F
we have by a well known leverage score approximate regression result with high probability in k:
kZX T − M M T Rk2F = O(1) · kZ(X ∗ )T − M M T Rk2F

= O(1) · kZZ T M M T R − M M T Rk2F

= O(∆1 )kM M T − (M M T )k k2F .

Further, computing X requires Õ(dk ω−1 ) time to compute the O(k log k) × k submatrix SM M T R
as well as Õ(kω ) = Õ(nkω−1 ) to perform the regression. This gives the result via Theorem
3 since computing Z with rank-2k ZX T gives a low-rank approximation of M M T with error
O(∆1 )kM M T − (M M T )2k k2F measured on the last k columns of M . Small error on these columns
is all that is needed to recover AC accurately (see proof of Theorem 3).

C

Additional lower bound proofs

We now prove our hardness result for kernels depending on the squared distance kai − aj k22 .
× Rd → R+ with ψ(ai , aj ) = f (kai − aj k2 ) for
Theorem 5. Consider any kernel function ψ : Rd P
q
q−1 and
some function f which can be expanded as f (x) = ∞
q=0 cq x with c1 6= 0 and |cq /c1 | ≤ G
for all q ≥ 2 and some G ≥ 1.
Assume there is an algorithm A which given input M ∈ Rn×d with kernel matrix K = {ψ(mi , mj )},
returns N ∈ Rn×k satisfying kK − N N T k2F ≤ ∆1 kK − Kk k in T (M, k) time.
For any A ∈ Rn×d , C ∈ Rd×k , with integer entries in [−∆2 , ∆2 ], let B = [w1 AT , w2 C]T
1 √
where w1 = 36√∆w2∆2 nd and w2 = (16Gd2 ∆4 )(36
. It is possible to compute AC in time
∆ ∆2 nd)
1

2

T (B, 2k + 3) + O(nk ω−1 ).

1

2

2

Proof. Define the distance matrix D ∈ Rn+k×n+k with Di,j = kbi − bj k2 . Using the fact that
kbi − bj k2 = kbi k2 + kbi k2 − 2bTi bj we have D = E + E T − 2BB T where E is a rank-1 matrix with
all rows equal to [kb1 k22 , ..., kbn+k k22 ]. We can write the kernel matrix for B and k as:




1 1
w1 w2 AC
w12 AAT
T
K = c0
+ c1 (E + E ) − 2c1
+ c2 D (2) + c3 D (3) + ...
(12)
1 1
w1 w2 C T AT w22 C T C
(q)

where Di,j = kbi − bj k2q . Let K̄ be K − c0 · 1 − c1 (E + E T ) , with its top n × n block set to 0. K̄
has rank at most 2k and if we set Q ∈ Rn×2k+3 to be a matrix with columns spanning the columns
15

of K̄, the all ones vector, E and E T , then letting N be the result of running A on B with rank
2k + 3:

 2
−2c1 w12 AAT + c2 D̂ (2) + ... 0
T 2
T
2
(13)
kK − N N kF ≤ ∆1 kK − QQ KkF ≤ ∆1
0
0 F
where D̂ (q) denotes the top left n × n submatrix of D (q) .
By our bounds on the entries of A, for i, j ≤ n, kbi − bj k2 ≤ 4d∆22 w12 and by our setting of
w1 , w2 , plugging into (13) we have for all i, j:
|(K − N N T )i,j | ≤ kK − N N T kF


∞
X
p
cq (4d∆22 w12 )q 
≤ ∆1 n 2c1 d∆22 w12 +

(14)

q=2

≤

p




∞
X
(4Gd∆22 w12 )q−1 
∆1 nc1 d∆22 w12 2 +

(Since |cq /c1 | ≤ Gq−1 )

q=2

p
w1 w2 c1
≤ 3 ∆1 nc1 d∆22 w12 ≤
12

(15)

where the second to last bound follows from the fact that w1 < w2 and w2 is set small enough so
(4Gd∆22 ) · w22 ≪ 1/2 so the series converges to a sum < 1. Additionally, for i ≤ n and j ≤ k (i.e.,
considering the entries of K corresponding to AC) we have:
T

Ki,n+j = c0 + c1 (E + E )i,n+j − 2c1 w1 w2 (AC)i,j +

∞
X

(q)

cq Di,n+j .

q=2

This last sum can be bounded by:
∞
X
q=2

(q)

cq Di,n+j ≤ c1

∞
X

Gq−1 (4∆22 dw22 )q

(By assumption |cq /c1 | ≤ Gq−1 )

q=2

≤ c1 w1 w2
≤ c1 w1 w2

∞
X

q=2
∞
X
q=2

2(q−1) w2

Gq−1 w2

w1

q
4∆22 d

q
Gq−1 w22q−3 4∆22 d

1
≤ c1 w1 w2 .
3

(Using w2 ≤

(Using
1/4
16G∆42 d2

w2
w1

≤

1
w2 .)

so the series converges.)

T
If we set v = N Ni,n+j
− c0 − c1 (E + E T )i,n+j we thus have combining with (14) for i ≤ n, j ≤ k

|v + 2c1 w1 w2 (AC)i,j | ≤

5c1 w1 w2
12

and so we can compute (AC)i,j exactly by rounding v to the nearest integer multiple of c1 w1 w2 .
This gives the theorem since we can compute the required entries of N N T and E in O(nk ω−1 )
time.

16

