arXiv:1608.07179v1 [cs.LG] 25 Aug 2016

Minimizing Quadratic Functions in Constant Time

Kohei Hayashi
National Institute of Advanced Industrial Science and Technology
hayashi.kohei@gmail.com
Yuichi Yoshida
National Institute of Informatics and Preferred Infrastructure, Inc.
yyoshida@nii.ac.jp

Abstract
A sampling-based optimization method for quadratic functions is proposed.
Our method approximately solves the following n-dimensional quadratic minimization problem in constant time, which is independent of n: z ∗ =
minv∈Rn hv, Avi + nhv, diag(d)vi + nhb, vi, where A ∈ Rn×n is a matrix and
d, b ∈ Rn are vectors. Our theoretical analysis specifies the number of samples
k(δ, ǫ) such that the approximated solution z satisfies |z − z ∗ | = O(ǫn2 ) with
probability 1 − δ. The empirical performance (accuracy and runtime) is positively
confirmed by numerical experiments.

1 Introduction
A quadratic function is one of the most important function classes in machine learning, statistics,
and data mining. Many fundamental problems such as linear regression, k-means clustering, principal component analysis, support vector machines, and kernel methods [14] can be formulated as a
minimization problem of a quadratic function.
In some applications, it is sufficient to compute the minimum value of a quadratic function rather
than its solution. For example, Yamada et al. [21] proposed an efficient method for estimating the
Pearson divergence, which provides useful information about data, such as the density ratio [18].
They formulated the estimation problem as the minimization of a squared loss and showed that the
Pearson divergence can be estimated from the minimum value. The least-squares mutual information [19] is another example that can be computed in a similar manner.
Despite its importance, the minimization of a quadratic function has the issue of scalability. Let
n ∈ N be the number of variables (the “dimension” of the problem). In general, such a minimization
problem can be solved by quadratic programming (QP), which requires poly(n) time. If the problem
is convex and there are no constraints, then the problem is reduced to solving a system of linear
equations, which requires O(n3 ) time. Both methods easily become infeasible, even for mediumscale problems, say, n > 10000.
Although several techniques have been proposed to accelerate quadratic function minimization, they
require at least linear time in n. This is problematic when handling problems with an ultrahigh
dimension, for which even linear time is slow or prohibitive. For example, stochastic gradient
descent (SGD) is an optimization method that is widely used for large-scale problems. A nice
property of this method is that, if the objective function is strongly convex, it outputs a point that
is sufficiently close to an optimal solution after a constant number of iterations [5]. Nevertheless,
in each iteration, we need at least O(n) time to access the variables. Another technique is lowrank approximation such as Nyström’s method [20]. The underlying idea is the approximation
of the problem by using a low-rank matrix, and by doing so, we can drastically reduce the time
29th Conference on Neural Information Processing Systems (NIPS 2016), Barcelona, Spain.

complexity. However, we still need to compute the matrix–vector product of size n, which requires
O(n) time. Clarkson et al. [7] proposed sublinear-time algorithms for special cases of quadratic
function minimization. However, it is “sublinear” with respect to the number of pairwise interactions
of the variables, which is O(n2 ), and their algorithms require O(n logc n) time for some c ≥ 1.
Our contributions: Let A ∈ Rn×n be a matrix and d, b ∈ Rn be vectors. Then, we consider the
following quadratic problem:
minimize
pn,A,d,b (v), where pn,A,d,b (v) = hv, Avi + nhv, diag(d)vi + nhb, vi.
(1)
n
v∈R

Here, h·, ·i denotes the inner product and diag(d) denotes the matrix whose diagonal entries are
specified by d. Note that a constant term can be included in (1); however, it is irrelevant when
optimizing (1), and hence we ignore it.
Let z ∗ ∈ R be the optimal value of (1) and let ǫ, δ ∈ (0, 1) be parameters. Then, the main goal of
this paper is the computation of z with |z − z ∗ | = O(ǫn2 ) with probability at least 1 − δ in constant
time, that is, independent of n. Here, we assume the real RAM model [6], in which we can perform
basic algebraic operations on real numbers in one step. Moreover, we assume that we have query
accesses to A, b, and d, with which we can obtain an entry of them by specifying an index. We note
that z ∗ is typically Θ(n2 ) because hv, Avi consists of Θ(n2 ) terms, and hv, diag(d)vi and hb, vi
consist of Θ(n) terms. Hence, we can regard the error of Θ(ǫn2 ) as an error of Θ(ǫ) for each term,
which is reasonably small in typical situations.
Let ·|S be an operator that extracts a submatrix (or subvector) specified by an index set S ⊂ N; then,
our algorithm is defined as follows, where the parameter k := k(ǫ, δ) will be determined later.
Algorithm 1
Input: An integer n ∈ N, query accesses to the matrix A ∈ Rn×n and to the vectors d, b ∈ Rn ,
and ǫ, δ > 0
1: S ← a sequence of k = k(ǫ, δ) indices independently and uniformly sampled from
{1, 2, . . . , n}.
2
n
2: return n
k2 minv∈R pk,A|S ,d|S ,b|S (v).
In other words, we sample a constant number of indices from the set {1, 2, . . . , n}, and then solve
the problem (1) restricted to these indices. Note that the number of queries and the time complexity
are O(k 2 ) and poly(k), respectively. In order to analyze the difference between the optimal values
of pn,A,d,b and pk,A|S ,d|S ,b|S , we want to measure the “distances” between A and A|S , d and d|S ,
and b and b|S , and want to show them small. To this end, we exploit graph limit theory, initiated by
Lovász and Szegedy [11] (refer to [10] for a book), in which we measure the distance between two
graphs on different number of vertices by considering continuous versions. Although the primary
interest of graph limit theory is graphs, we can extend the argument to analyze matrices and vectors.
Using synthetic and real settings, we demonstrate that our method is orders of magnitude faster than
standard polynomial-time algorithms and that the accuracy of our method is sufficiently high.
Related work: Several constant-time approximation algorithms are known for combinatorial optimization problems such as the max cut problem on dense graphs [8, 13], constraint satisfaction
problems [1, 22], and the vertex cover problem [15, 16, 25]. However, as far as we know, no such
algorithm is known for continuous optimization problems.
A related notion is property testing [9, 17], which aims to design constant-time algorithms that
distinguish inputs satisfying some predetermined property from inputs that are “far” from satisfying
it. Characterizations of constant-time testable properties are known for the properties of a dense
graph [2, 3] and the affine-invariant properties of a function on a finite field [23, 24].
Organization In Section 2, we introduce the basic notions from graph limit theory. In Section 3,
we show that we can obtain a good approximation to (a continuous version of) a matrix by sampling a
constant-size submatrix in the sense that the optimizations over the original matrix and the submatrix
are essentially equivalent. Using this fact, we prove the correctness of Algorithm 1 in Section 4. We
show our experimental results in Section 5.
2

2 Preliminaries
For an integer n, let [n] denote the set {1, 2, . . . , n}. The notation a = b ± c means that b − c ≤ a ≤
b + c. In this paper, we only consider functions and sets that are measurable.
Let S = (x1 , . . . , xk ) be a sequence of k indices in [n]. For a vector v ∈ Rn , we denote the
restriction of v to S by v|S ∈ Rk ; that is, (v|S )i = vxi for every i ∈ [k]. For the matrix A ∈ Rn×n ,
we denote the restriction of A to S by A|S ∈ Rk×k ; that is, (A|S )ij = Axi xj for every i, j ∈ [k].
2.1 Dikernels
Following [12], we call a (measurable) function f : [0, 1]2 → R a dikernel. A dikernel is a generalization of a graphon [11], which is symmetric and whose range is bounded in [0, 1]. We can regard a
dikernel as a matrix whose index is specified by a real value in [0, 1]. We stress that the term dikernel
has nothing to do with kernel methods.
R1
For two functions f, g : [0, 1] → R, we define their inner product as hf, gi = 0 f (x)g(x)dx. For a
dikernel W : [0, 1]2 → R and a function f : [0, 1] → R, we define a function W f : [0, 1] → R as
(W f )(x) = hW (x, ·), f i.
Let W : [0, 1]2 → R be a dikernel. The Lp norm kW kp for p ≥ 1 and the cut norm kW k of W are
R R
1/p
R R
1 1
defined as kW kp = 0 0 |W (x, y)|p dxdy
and kW k = supS,T ⊆[0,1] S T W (x, y)dxdy ,
respectively, where the supremum is over all pairs of subsets. We note that these norms satisfy the
triangle inequalities and kW k ≤ kW k1 .
Let λ be a Lebesgue measure. A map π : [0, 1] → [0, 1] is said to be measure-preserving, if
the pre-image π −1 (X) is measurable for every measurable set X, and λ(π −1 (X)) = λ(X). A
measure-preserving bijection is a measure-preserving map whose inverse map exists and is also
measurable (and then also measure-preserving). For a measure preserving bijection π : [0, 1] →
[0, 1] and a dikernel W : [0, 1]2 → R, we define the dikernel π(W ) : [0, 1]2 → R as π(W )(x, y) =
W (π(x), π(y)).
2.2 Matrices and Dikernels
Let W : [0, 1]2 → R be a dikernel and S = (x1 , . . . , xk ) be a sequence of elements in [0, 1]. Then,
we define the matrix W |S ∈ Rk×k so that (W |S )ij = W (xi , xj ).
b : [0, 1]2 → R from the matrix A ∈ Rn×n as follows. Let I1 =
We can construct the dikernel A
1
1 2
[0, n ], I2 = ( n , n ], . . . , In = ( n−1
n , . . . , 1]. For x ∈ [0, 1], we define in (x) ∈ [n] as a unique
b y) = Ai (x)i (y) . The main motivation for creating a
integer such that x ∈ Ii . Then, we define A(x,
n
n
dikernel from a matrix is that, by doing so, we can define the distance between two matrices A and
b − Bk
b .
B of different sizes via the cut norm, that is, kA
We note that the distribution of A|S , where S is a sequence of k indices that are uniformly and
b S , where S is a sequence of
independently sampled from [n] exactly matches the distribution of A|
k elements that are uniformly and independently sampled from [0, 1].

3 Sampling Theorem and the Properties of the Cut Norm
In this section, we prove the following theorem, which states that, given a sequence of dikernels
W 1 , . . . , W T : [0, 1]2 → [−L, L], we can obtain a good approximation to them by sampling a
sequence of a small number of elements in [0, 1]. Formally, we prove the following:
Theorem 3.1. Let W 1 , . . . , W T : [0, 1]2 → [−L, L] be dikernels. Let S be a sequence of k
elements uniformly and independently sampled from [0, 1]. Then, with a probability of at least
1 − exp(−Ω(kT / log2 k)), there exists a measure-preserving bijection π : [0, 1] → [0, 1] such that,
for any functions f, g : [0, 1] → [−K, K] and t ∈ [T ], we have


p
[
t | )gi| = O LK 2 T / log k .
|hf, W t gi − hf, π(W
S
2
3

We start with the following lemma, which states that, if a dikernel W : [0, 1]2 → R has a small cut
norm, then hf, W f i is negligible no matter what f is. Hence, we can focus on the cut norm when
proving Theorem 3.1.
Lemma 3.2. Let ǫ ≥ 0 and W : [0, 1]2 → R be a dikernel with kW k ≤ ǫ. Then, for any functions
f, g : [0, 1] → [−K, K], we have |hf, W gi| ≤ ǫK 2 .
Proof. For τ ∈ R and the function h : [0, 1] → R, let Lτ (h) := {x ∈ [0, 1] | h(x) = τ } be the level
set of h at τ . For f ′ = f /K and g ′ = g/K, we have
Z 1Z 1
Z
Z
2
′
′
2
|hf, W gi| = K |hf , W g i| = K
τ1 τ2
W (x, y)dxdydτ1 dτ2
−1

≤ K2

Z

1

−1

≤ ǫK 2

Z

Z

1

−1

1

|τ1 ||τ2 |

−1

Z

Z

−1

Lτ1 (f ′ )

1

Z

Lτ1 (f ′ )

Lτ2 (g′ )

W (x, y)dxdy dτ1 dτ2

Lτ2 (g′ )

|τ1 ||τ2 |dτ1 dτ2 = ǫK 2 .

−1

To introduce the next technical tool, we need several definitions. We say that the partition Q is a
refinement of the partition P = (V1 , . . . , Vp ) if Q is obtained by splitting each set Vi into one or more
parts. The partition P = (V1 , . . . , Vp ) of the interval [0, 1] is called an equipartition if λ(Vi ) = 1/p
for every i ∈ [p]. For the dikernel W : [0, 1]2 → R and the equipartition P = (V1 , . . . , Vp ) of [0, 1],
we define WP : [0, 1]2 → R as the function obtained by averaging each Vi × Vj for i, j ∈ [p]. More
formally, we define
Z
Z
1
WP (x, y) =
W (x′ , y ′ )dx′ dy ′ ,
W (x′ , y ′ )dx′ dy ′ = p2
λ(Vi )λ(Vj ) Vi ×Vj
Vi ×Vj
where i and j are unique indices such that x ∈ Vi and y ∈ Vj , respectively.
The following lemma states that any function W : [0, 1]2 → R can be well approximated by WP
for the equipartition P into a small number of parts.
Lemma 3.3 (Weak regularity lemma for functions on [0, 1]2 [8]). Let P be an equipartition of [0, 1]
into k sets. Then, for any dikernel W : [0, 1]2 → R and ǫ > 0, there exists a refinement Q of P with
2
|Q| ≤ k2C/ǫ for some constant C > 0 such that
kW − WQ k ≤ ǫkW k2.
Corollary 3.4. Let W 1 , . . . , W T : [0, 1]2 → R be dikernels. Then, for any ǫ > 0, there exists an
2
equipartition P into |P| ≤ 2CT /ǫ parts for some constant C > 0 such that, for every t ∈ [T ],
kW t − WPt k ≤ ǫkW t k2 .
Proof. Let P 0 be a trivial partition, that is, a partition consisting of a single part [n]. Then, for each
t ∈ [T ], we iteratively apply Lemma 3.3 with P t−1 , W t , and ǫ, and we obtain the partition P t into
2
at most |P t−1 |2C/ǫ parts such that kW t − WPt t k ≤ ǫkW t k2 . Since P t is a refinement of P t−1 ,
we have kW i − WPi t k ≤ kW i − WPi t−1 k for every i ∈ [t − 1]. Then, P T satisfies the desired
2
2
property with |P T | ≤ (2C/ǫ )T = 2CT /ǫ .
d
As long as S is sufficiently large, W and W
|S are close in the cut norm:

Lemma 3.5 ((4.15) of [4]). Let W : [0, 1]2 → [−L, L] be a dikernel and S be a sequence of k
elements uniformly and independently sampled from [0, 1]. Then, we have
−

8L
2L
d
≤ ES k W
|S k − kW k < 1/4 .
k
k

Finally, we need the following concentration inequality.
4

Lemma 3.6 (Azuma’s inequality). Let (Ω, A, P ) be a probability space, k be a positive integer, and
C > 0. Let z = (z1 , . . . , zk ), where z1 , . . . , zk are independent random variables, and zi takes
values in some measure space (Ωi , Ai ). Let f : Ω1 × · · · × Ωk → R be a function. Suppose that
|f (x) − f (y)| ≤ C whenever x and y only differ in one coordinate. Then
h
i
2
Pr |f (z) − Ez [f (z)]| > λC < 2e−λ /2k .
Now we prove the counterpart of Theorem 3.1 for the cut norm.
Lemma 3.7. Let W 1 , . . . , W T : [0, 1]2 → [−L, L] be dikernels. Let S be a sequence of k
elements uniformly and independently sampled from [0, 1]. Then, with a probability of at least
1 − exp(−Ω(kT / log2 k)), there exists a measure-preserving bijection π : [0, 1] → [0, 1] such that,
for every t ∈ [T ], we have
 p

[
t | )k = O L T / log k .
kW t − π(W
S 
2

Proof. First, we bound the expectations and then prove their concentrations. We apply Corollary 3.4
2
to W 1 , . . . , W T and ǫ, and let P = (V1 , . . . , Vp ) be the obtained partition with p ≤ 2CT /ǫ parts
such that
kW t − WPt k ≤ ǫL.
for every t ∈ [T ]. By Lemma 3.5, for every t ∈ [T ], we have
8L
t
t
t
\
[
t
ES k W
\
P − W )|S k ≤ ǫL + 1/4 .
P |S − W |S k = ES k(W
k
Then, for any measure-preserving bijection π : [0, 1] → [0, 1] and t ∈ [T ], we have
t
t
\
\
t | )k ≤ kW t − W t k + E kW t − π(W
t
[
[
ES kW t − π(W
S 
S
P 
P
P |S )k + ES kπ(WP |S ) − π(W |S )k
8L
t
\
≤ 2ǫL + 1/4 + ES kWPt − π(W
(2)
P |S )k .
k

Thus, we are left with the problem of sampling from P. Let S = {x1 , . . . , xk } be a sequence of
independent random variables that are uniformly distributed in [0, 1], and let Zi be the number of
points xj that fall into the set Vi . It is easy to compute that
1
k
k
1
E[Zi ] =
and Var[Zi ] =
− 2 k< .
p
p p
p
The partition P ′ of [0, 1] is constructed into the sets V1′ , . . . , Vp′ such that λ(Vi′ ) = Zi /k and λ(Vi ∩
t

Vi′ ) = min(1/p, Zi /k). For each t ∈ [T ], we construct the dikernel W : [0, 1] → R such that the
t
t
value of W on Vi′ × Vj′ is the same as the value of WPt on Vi × Vj . Then, W agrees with WPt on
S
t
t
\
the set Q = i,j∈[p] (Vi ∩Vi′ )× (Vj ∩Vj′ ). Then, there exists a bijection π such that π(W
P |S ) = W
for each t ∈ [T ]. Then, for every t ∈ [T ], we have
t
t
t
t
t
\
kWPt − π(W
P |S )k = kWP − W k ≤ kWP − W k1 ≤ 2L(1 − λ(Q))
 1 Z 

 1 Z 2 

X
X
i
i
min ,
≤ 4L 1 −
min ,
= 2L 1 −
p k
p k
i∈[p]

i∈[p]

 X  1 Z 2 1/2
X 1 Zi
i
= 2L
,
−
−
≤ 2L p
p
k
p
k
i∈[p]

i∈[p]

which we rewrite as
2
2
t
\
kWPt − π(W
P |S )k ≤ 4L p

P

X1

i∈[p]

p

−

Zi 2
.
k

The expectation of the right hand side is (4L2 p/k 2 ) i∈[p] Var(Zi ) < 4L2 p/k. By the Cauchyp
t | )k ≤ 2L p/k.
\
Schwartz inequality, EkW t − π(W
S 
P

P

5

Inserted this into (2), we obtain

r
2
8L
8L
2L
p
[
t
≤ 2ǫL + 1/4 + 1/2 2CT /ǫ .
EkW − π(W |S )k ≤ 2ǫL + 1/4 + 2L
k
k
k
k
p
p
Choosing ǫ = CT /(log2 k 1/4 ) = 4CT /(log2 k), we obtain the upper bound
s
s

8L
4CT
T 
2L
[
t | )k ≤ 2L
+ 1/4 + 1/4 = O L
.
EkW t − π(W
S 
log2 k k
log2 k
k
t

t | )k changes by at most O(L/k) if one element in S changes, we
[
Observing that kW t − π(W
S  p
apply Azuma’s inequality with λ = k T / log2 k and the union bound to complete the proof.

The proof of Theorem 3.1 is immediately follows from Lemmas 3.2 and 3.7.

4 Analysis of Algorithm 1
In this section, we analyze Algorithm 1. Because we want to use dikernels for the analysis, we
introduce a continuous version of pn,A,d,b (recall (1)). The real-valued function Pn,A,d,b on the
functions f : [0, 1] → R is defined as
d⊤ 1i + hf, b1
d
⊤ 1i,
b i + hf 2 , d1
Pn,A,d,b (f ) = hf, Af

where f 2 : [0, 1] → R is a function such that f 2 (x) = f (x)2 for every x ∈ [0, 1] and 1 : [0, 1] → R
is the constant function that has a value of 1 everywhere. The following lemma states that the
minimizations of pn,A,d,b and Pn,A,d,b are equivalent:
Lemma 4.1. Let A ∈ Rn×n be a matrix and d, b ∈ Rn×n be vectors. Then, we have
min

v∈[−K,K]n

pn,A,d,b (v) = n2 ·

inf

f :[0,1]→[−K,K]

Pn,A,d,b (f ).

for any K > 0.
Proof. First, we show that n2 · inf f :[0,1]→[−K,K] Pn,A,d,b (f ) ≤ minv∈[−K,K]n pn,A,d,b (v). Given
a vector v ∈ [−K, K]n , we define f : [0, 1] → [−K, K] as f (x) = vin (x) . Then,
X Z Z
1 X
1
b i=
hf, Af
Aij f (x)f (y)dxdy = 2
Aij vi vj = 2 hv, Avi,
n
n
i,j∈[n] Ii Ij
i,j∈[n]
X Z Z
XZ
1
1 X
d⊤ 1i =
di vi2 = hv, diag(d)vi,
hf 2 , d1
di f (x)2 dx =
di f (x)2 dxdy =
n
n
i∈[n]
i,j∈[n] Ii Ij
i∈[n] Ii
X Z Z
XZ
X
1
1
d
⊤ 1i =
hf, b1
bi f (x)dx =
bi f (x)dxdy =
bi vi = hv, bi.
n
n
Ii Ij
Ii
i,j∈[n]

i∈[n]

i∈[n]

2

Then, we have n Pn,A,d,b (f ) ≤ pn,A,d,b (v).

Next, we show that minv∈[−K,K]n pn,A,d,b (v) ≤ n2 · inf f :[0,1]→[−K,K] Pn,A,d,b (f ). Let f :
[0, 1] → [−K, K] be a measurable function. Then, for x ∈ [0, 1], we have
XZ
XZ
∂Pn,A,d,b (f (x))
=
Aiin (x) f (y)dy +
Ain (x)j f (y)dy + 2din (x) f (x) + bin (x) .
∂f (x)
Ii
Ij
i∈[n]

j∈[n]

Note that the form of this partial derivative only depends on in (x); hence, in the optimal solution
f ∗ : [0, 1] → [−K, K], we can assume f ∗ (x) = f ∗ (y) if in (x) = in (y). In other words, f ∗
is constant on each of the intervals I1 , . . . , In . For such f ∗ , we define the vector v ∈ Rn as
vi = f ∗ (x), where x ∈ [0, 1] is any element in Ii . Then, we have
X Z Z
X
b ∗ i,
Aij vi vj = n2
Aij f ∗ (x)f ∗ (y)dxdy = n2 hf ∗ , Af
hv, Avi =
i,j∈[n]

Ii

i,j∈[n]

6

Ij

hv, diag(d)vi =

X

di vi2

=n

i∈[n]

i∈[n]

hv, bi =

X

XZ

bi vi = n

i∈[n]

XZ

i∈[n]

Ii

Ii

dT 1i,
di f ∗ (x)2 dx = nh(f ∗ )2 , d1

dT 1i.
bi f ∗ (x)dx = nhf ∗ , b1

Finally, we have pn,A,d,b (v) ≤ n2 Pn,A,d,b (f ∗ ).

Now we show that Algorithm 1 well-approximates the optimal value of (1) in the following sense:
Theorem 4.2. Let v ∗ and z ∗ be an optimal solution and the optimal value, respectively, of prob2
lem (1). By choosing k(ǫ, δ) = 2Θ(1/ǫ ) + Θ(log 1δ log log δ1 ), with a probability of at least
1 − δ, a sequence S of k indices independently and uniformly sampled from [n] satisfies the following: Let ṽ ∗ and z̃ ∗ be an optimal solution and the optimal value, respectively, of the problem
minv∈Rk pk,A|S ,d|S ,b|S (v). Then, we have
n2 ∗
z̃ − z ∗ ≤ ǫLK 2 n2 ,
k2
where K = max{maxi∈[n] |vi∗ |, maxi∈[n] |ṽi∗ |} and L = max{maxi,j |Aij |, maxi |di |, maxi |bi |}.
2
b
Proof. We instantiate Theorem 3.1 with k = 2Θ(1/ǫ ) + Θ(log δ1 log log δ1 ) and the dikernels A,
d
d
⊤
⊤
d1 , and b1 . Then, with a probability of at least 1 − δ, there exists a measure preserving bijection
π : [0, 1] → [0, 1] such that
n
o ǫLK 2
2
d⊤
d
d
\
\
⊤
⊤
⊤
b − π(A|
max |hf, (A
S ))f i|, |hf , (d1 − π(d1 |S ))1i|, |hf, (b1 − π(b1 |S ))1i| ≤
3

for any function f : [0, 1] → [−K, K]. Then, we have

z̃ ∗ = min pk,A|S ,d|S ,b|S (v) =
v∈Rk

= k2 ·

min

v∈[−K,K]k

pk,A|S ,d|S ,b|S (v)

(By Lemma 4.1)
Pk,A|S ,d|S ,b|S (f )

2
2 d
d⊤
d
\
⊤
⊤
b
b
hf, (π(A|
= k2 ·
inf
S ) − A)f i + hf, Af i + hf , (π(d1 |S ) − d1 )1i + hf , d1 1i+
f :[0,1]→[−K,K]

d
d
\
⊤ | ) − b1
⊤ )1i + hf, b1
⊤ 1i
hf, (π(b1
S


d⊤ 1i + hf, b1
d
⊤ 1i ± ǫLK 2
b i + hf 2 , d1
hf, Af
≤ k2 ·
inf
inf

f :[0,1]→[−K,K]

f :[0,1]→[−K,K]

2

=

k
· min pn,A,d,b (v) ± ǫLK 2 k 2 .
n2 v∈[−K,K]n

=

k2
k2
· minn pn,A,d,b (v) ± ǫLK 2 k 2 = 2 z ∗ ± ǫLK 2 k 2 .
2
n v∈R
n

(By Lemma 4.1)

Rearranging the inequality, we obtain the desired result.
We can show that K is bounded when A is symmetric and full rank. To see this, we first note
that we can assume A + ndiag(d) is positive-definite, as otherwise pn,A,d,b is not bounded and
the problem is uninteresting. Then, for any set S ⊆ [n] of k indices, (A + ndiag(d))|S is again
positive-definite because it is a principal submatrix. Hence, we have v ∗ = (A + ndiag(d))−1 nb/2
and ṽ ∗ = (A|S + ndiag(d|S ))−1 nb|S /2, which means that K is bounded.

5 Experiments
In this section, we demonstrate the effectiveness of our method by experiment. All experiments
were conducted on an Amazon EC2 c3.8xlarge instance. Error bars indicate the standard deviations
over ten trials with different random seeds.
7

n=200

0.10
0.05

Table 1: Pearson divergence: runtime (second).

0.05
0.00
0.10

1000

|z − z∗| n2

500

0.10

0.05
0.00

2000

0.10
0.05
0.00
10

20

40

80

160

k

Nyström Proposed

0.00

k
20
40
80
160
20
40
80
160

n = 500
0.002
0.003
0.007
0.030
0.005
0.010
0.022
0.076

1000
0.002
0.003
0.007
0.030
0.012
0.022
0.049
0.116

2000
0.002
0.003
0.008
0.033
0.046
0.087
0.188
0.432

5000
0.002
0.003
0.008
0.035
0.274
0.513
0.942
1.972

Figure 1: Numerical simulation: absolute approximation error scaled by n2 .

Nyström Proposed

Table 2: Pearson divergence: absolute approximation error.
k
20
40
80
160
20
40
80
160

n = 500
0.0027 ± 0.0028
0.0018 ± 0.0023
0.0007 ± 0.0008
0.0003 ± 0.0003
0.3685 ± 0.9142
0.3549 ± 0.6191
0.0184 ± 0.0192
0.0143 ± 0.0209

1000
0.0012 ± 0.0012
0.0006 ± 0.0007
0.0004 ± 0.0003
0.0002 ± 0.0001
1.3006 ± 2.4504
0.4207 ± 0.7018
0.0398 ± 0.0472
0.0348 ± 0.0541

2000
0.0021 ± 0.0019
0.0012 ± 0.0011
0.0008 ± 0.0008
0.0003 ± 0.0003
3.1119 ± 6.1464
0.9838 ± 1.5422
0.2056 ± 0.2725
0.0585 ± 0.1112

5000
0.0016 ± 0.0022
0.0011 ± 0.0020
0.0007 ± 0.0017
0.0002 ± 0.0003
0.6989 ± 0.9644
0.3744 ± 0.6655
0.5705 ± 0.7918
0.0254 ± 0.0285

Numerical simulation We investigated the actual relationships between n, k, and ǫ. To this end,
we prepared synthetic data as follows. We randomly generated inputs as Aij ∼ U[−1,1] , di ∼ U[0,1] ,
and bi ∼ U[−1,1] for i, j ∈ [n], where U[a,b] denotes the uniform distribution with the support [a, b].
After that, we solved (1) by using Algorithm 1 and compared it with the exact solution obtained by
QP.1 The result (Figure 1) show the approximation errors were evenly controlled regardless of n,
which meets the error analysis (Theorem 4.2).
Application to kernel methods Next, we considered the kernel approximation of the Pearson
divergence [21]. The problem is defined as follows. Suppose we have the two different data sets
′
x = (x1 , . . . , xn ) ∈ Rn and x′ = (x′1 , . . . , x′n′ ) ∈ Rn where n, n′ ∈ N. Let H ∈ Rn×n
Pn′
Pn
1−α
′
′
be a gram matrix such that Hl,m = α
j=1 φ(xj , xl )φ(xj , xm ),
i=1 φ(xi , xl )φ(xi , xm ) + n′
n
n
where φ(·, ·) is aP
kernel function and α ∈ (0, 1) is a parameter. Also, let h ∈ R be a vector
n
such that hl = n1 i=1 φ(xi , xl ). Then, an estimator of the α-relative Pearson divergence between
the distributions of x and x′ is obtained by − 21 − minv∈Rn 12 hv, Hvi − hh, vi + λ2 hv, vi. Here,
λ > 0 is a regularization parameter. In this experiment, we used the Gaussian kernel φ(x, y) =
exp((x − y)2 /2σ 2 ) and set n′ = 200 and α = 0.5; σ 2 and λ were chosen by 5-fold cross-validation
as suggested in [21]. We randomly generated the data sets as xi ∼ N (1, 0.5) for i ∈ [n] and
x′j ∼ N (1.5, 0.5) for j ∈ [n′ ] where N (µ, σ 2 ) denotes the Gaussian distribution with mean µ and
variance σ 2 .
λ
We encoded this problem into (1) by setting A = 21 H, b = −h, and d = 2n
1n , where 1n denotes
the n-dimensional vector whose elements are all one. After that, given k, we computed the second
step of Algorithm 1 with the pseudoinverse of A|S +kdiag(d|S ). Absolute approximation errors and
runtimes were compared with Nyström’s method whose approximated rank was set to k. In terms of
accuracy, our method clearly outperformed Nyström’s method (Table 2). In addition, the runtimes
of our method were nearly constant, whereas the runtimes of Nyström’s method grew linearly in k
(Table 1).

1

We used GLPK (https://www.gnu.org/software/glpk/) for the QP solver.

8

6 Acknowledgments
We would like to thank Makoto Yamada for suggesting a motivating problem of our method. K. H. is
supported by MEXT KAKENHI 15K16055. Y. Y. is supported by MEXT Grant-in-Aid for Scientific
Research on Innovative Areas (No. 24106001), JST, CREST, Foundations of Innovative Algorithms
for Big Data, and JST, ERATO, Kawarabayashi Large Graph Project.

References
[1] N. Alon, W. F. de la Vega, R. Kannan, and M. Karpinski. Random sampling and approximation of MAXCSP problems. In STOC, pages 232–239, 2002.
[2] N. Alon, E. Fischer, I. Newman, and A. Shapira. A combinatorial characterization of the testable graph
properties: It’s all about regularity. SIAM Journal on Computing, 39(1):143–167, 2009.
[3] C. Borgs, J. Chayes, L. Lovász, V. T. Sós, B. Szegedy, and K. Vesztergombi. Graph limits and parameter
testing. In STOC, pages 261–270, 2006.
[4] C. Borgs, J. T. Chayes, L. Lovász, V. T. Sós, and K. Vesztergombi. Convergent sequences of dense graphs
i: Subgraph frequencies, metric properties and testing. Advances in Mathematics, 219(6):1801–1851,
2008.
[5] L. Bottou. Stochastic learning. In Advanced Lectures on Machine Learning, pages 146–168. 2004.
[6] V. Brattka and P. Hertling. Feasible real random access machines. Journal of Complexity, 14(4):490–526,
1998.
[7] K. L. Clarkson, E. Hazan, and D. P. Woodruff. Sublinear optimization for machine learning. Journal of
the ACM, 59(5):23:1–23:49, 2012.
[8] A. Frieze and R. Kannan. The regularity lemma and approximation schemes for dense problems. In
FOCS, pages 12–20, 1996.
[9] O. Goldreich, S. Goldwasser, and D. Ron. Property testing and its connection to learning and approximation. Journal of the ACM, 45(4):653–750, 1998.
[10] L. Lovász. Large Networks and Graph Limits. American Mathematical Society, 2012.
[11] L. Lovász and B. Szegedy. Limits of dense graph sequences. Journal of Combinatorial Theory, Series B,
96(6):933–957, 2006.
[12] L. Lovász and K. Vesztergombi. Non-deterministic graph property testing. Combinatorics, Probability
and Computing, 22(05):749–762, 2013.
[13] C. Mathieu and W. Schudy. Yet another algorithm for dense max cut: go greedy. In SODA, pages 176–182,
2008.
[14] K. P. Murphy. Machine learning: a probabilistic perspective. The MIT Press, 2012.
[15] H. N. Nguyen and K. Onak. Constant-time approximation algorithms via local improvements. In FOCS,
pages 327–336, 2008.
[16] K. Onak, D. Ron, M. Rosen, and R. Rubinfeld. A near-optimal sublinear-time algorithm for approximating the minimum vertex cover size. In SODA, pages 1123–1131, 2012.
[17] R. Rubinfeld and M. Sudan. Robust characterizations of polynomials with applications to program testing.
SIAM Journal on Computing, 25(2):252–271, 1996.
[18] M. Sugiyama, T. Suzuki, and T. Kanamori. Density Ratio Estimation in Machine Learning. Cambridge
University Press, 2012.
[19] T. Suzuki and M. Sugiyama. Least-Squares Independent Component Analysis. Neural Computation,
23(1):284–301, 2011.
[20] C. K. I. Williams and M. Seeger. Using the nyström method to speed up kernel machines. In NIPS, 2001.
[21] M. Yamada, T. Suzuki, T. Kanamori, H. Hachiya, and M. Sugiyama. Relative density-ratio estimation for
robust distribution comparison. In NIPS, 2011.
[22] Y. Yoshida. Optimal constant-time approximation algorithms and (unconditional) inapproximability results for every bounded-degree CSP. In STOC, pages 665–674, 2011.
[23] Y. Yoshida. A characterization of locally testable affine-invariant properties via decomposition theorems.
In STOC, pages 154–163, 2014.
[24] Y. Yoshida. Gowers norm, function limits, and parameter estimation. In SODA, pages 1391–1406, 2016.
[25] Y. Yoshida, M. Yamamoto, and H. Ito. Improved constant-time approximation algorithms for maximum
matchings and other optimization problems. SIAM Journal on Computing, 41(4):1074–1093, 2012.

9

