Learning with Cross-Kernels and Ideal PCA

arXiv:1406.2646v1 [cs.LG] 10 Jun 2014

Franz J. Király∗

Martin Kreuzer †

Louis Theran‡

Abstract
We describe how cross-kernel matrices, that is, kernel matrices between the data and a
custom chosen set of ‘feature spanning points’ can be used for learning. The main potential of
cross-kernels lies in the fact that (a) only one side of the matrix scales with the number of data
points, and (b) cross-kernels, as opposed to the usual kernel matrices, can be used to certify
for the data manifold. Our theoretical framework, which is based on a duality involving
the feature space and vanishing ideals, indicates that cross-kernels have the potential to be
used for any kind of kernel learning. We present a novel algorithm, Ideal PCA (IPCA), which
cross-kernelizes PCA. We demonstrate on real and synthetic data that IPCA allows to (a)
obtain PCA-like features faster and (b) to extract novel and empirically validated features
certifying for the data manifold.

1. Introduction
Since their invention by Boser, Guyon and Vapnik [2, 7], kernel methods have had a fundamental
impact on the fields of statistics and machine learning. The major appeal of using kernel methods for learning consists in using the kernel trick, first proposed by Aizerman, Braverman and
Rozonoer [1], which allows to make otherwise costly computations in the feature space implicit
and thus highly efficient for a huge variety of learning tasks – see e.g. [4, 6] for an overview.
Many kernel methods make extensive use of the so-called kernel matrix, a matrix whose
entries are kernel function evaluations k(x i , x j ) at data points x 1 , . . . , x N ∈ Rn . This kernel matrix
KXX = (k(x i , x j ))i j is simultaneously the main source of efficient linearization and the central
computational bottleneck. For instance, the most expensive part of algorithms such as kernel PCA
or kernel ridge regression consists in computing the inverse or a singular value decomposition of
the (N × N ) matrix KXX . We make two observations in this context:
(A) Subsampling KXX , i.e., considering (N × M ) sub-matrices of KXX , is the state-of-the-art in
speeding up the scalability in the number of data points N . While much has been written on the
topic, there seems to be no consensus on how exactly to choose the subsample from the data set.
(B) There seems to be no kernel algorithm which learns the data manifold from which the
x i were obtained, or test whether an unseen data point is on the manifold - in particular none
involving KXX . (Note that kernel PCA outputs projections on predominant data coordinates, but
not the coordinates embedded in data space or the shape of the data manifold.)
∗

Department of Statistical Science, University College, Gower Street, London WC1E 6BT, United Kingdom, f.

kiraly@ucl.ac.uk
†

Faculty of Computer Science and Mathematics, University of Passau, 94030 Passau, Germany, martin.kreuzer@

uni-passau.de
‡

Institut für Mathematik, Freie Universität Berlin, Arnimallee 2, 14195 Berlin, Germany, theran@math.

fu-berlin.de

1

In this paper, we propose an algebraic duality framework which offers both a potential explanation and a potential solution to the above (seemingly unrelated) issues. Our contribution
is centered around cross-kernel matrices KXZ = (k(x i , z j ))i j , where z1 , . . . , z M ∈ Rn , called crosskernel basis points, are not necessarily data points. Assuming that z1 , . . . , z M are chosen suitably
(e.g., randomly) and that M is large enough, we show theoretically and empirically that:
(1) The cross-kernel matrix KXZ can be used to capture both the information contained in
KXX as well as additional information about the data manifold which is provably not contained
in KXX . Practically, we present an algorithm, called IPCA, which is able to learn both coordinates
in the manifold as well as certificates for being contained in the data manifold, addressing issue (B)
above.
(2) The cross-kernel matrix KXZ can be employed in ways completely analogous to a subsampled matrix KXX . As opposed to subsampling, the properties of the points zi can be prescribed by
the experimenter and are thus independent of known or unknown – and potentially detrimental
– properties of the data x i . In particular, KXZ can be used for speeding up kernel learning, while
potentially avoiding data-related subsampling issues, thereby addressing issue (A) above.
In the following we briefly explain why and how considering the cross-kernel KXZ can be
advantageous as compared to the kernel matrix KXX . More detailed technical statements and
explicit algorithms can be found in section 2 and thereafter.
Why the kernel matrix is not enough. The kernel matrix KXX misses information on the data
manifold that can be obtained from KXZ . Let us explain this in the example where the kernel is
the ordinary Euclidean scalar product k(x, y) = 〈x, y〉 in Rn . Suppose the data points x 1 , . . . , x N
all lie in a vector subspace L ⊆ Rn . The entries of the kernel matrix KXX are scalar products of the
type 〈x i , x j 〉 which are invariant under rotation of the coordinate system. Therefore, one obtains
one and the same kernel matrix KXX when considering rotated data points U x 1 , . . . , U x N and the
rotated vector space U · L for any rotation matrix U ∈ Rn×n . Since one and the same kernel
matrix KXX can arise in a non-degenerate way from every vector space L (of that dimension), the
data manifold L can not be obtained back from KXX . Similar considerations hold for arbitrary
kernels: the data manifold can not be learnt from the kernel matrix KXX .
Cross-kernels help! In the above example the vector space L can be identified from a crosskernel matrix KXZ = (k(x i , z j ))i j , if z1 , . . . , z M span Rn . Observe that the rows of KXZ are coordinate representations of the x i . Since z1 , . . . , z M span Rn , the data point x i can be obtained back
from the i-th row of KXZ . Since reconstructing x i from the i-th row of KXZ is a linear reparameterization, an unseen data point x is contained in the data manifold L if and only if the vector
κz (x) = (k(x, z j )) j is contained in the row-span of KXZ . Therefore we can obtain L back from KXZ
and use the row-span of KXZ to efficiently test membership for the data manifold L. Note that the
same reasoning does not work for the matrix KXX , since x 1 , . . . , x N do not span Rn , but only L.
A similar reasoning holds in larger generality: the cross-kernel KXZ contains extensive information
on the data manifold which KXX does not.
From cross-kernels to kernel learning. The “manifold awareness” of the cross-kernel implies further interesting properties which can be used for decomposing, and thus approximating,
>
the original kernel matrix KXX via the potentially much smaller KXZ . Namely, letting KZX = KXZ
−1
and KZZ = (k(zi , z j ))i j , we can show that the equality KXX = KXZ KZZ
KZX holds exactly for polynomial kernels (assuming “feature-spanning” points zi as in the example above, see Theorem 3.4).
We conjecture (and have observed empirically) that the same holds approximately for other kernels. This equality is reminiscent of the central equation in subsampling, with the difference that

2

the points z j do not need to be chosen among the data. Suitable choices of the points z j avoid
potential problems in the data while leading to the same speed-up of arbitrary kernel methods.

2. Some dualities for polynomial kernels
The following preview of our approach is a more technically exact and quantitative variant of
the discussion in the introduction. In what follows, we use inhomogenous polynomial kernels.
(A further difficulty, not considered here, lies in the generalization to other kernels and kernel
feature spaces.) The central observation which we to carry over from the introduction is the
certifying property of the row-span of KXZ . In our setup, the data points x i are sampled from a
manifold X ⊂ Rn which is “cut out” by an unknown set of polynomial equations. Our goal is to
relate these equations to the linear span of the images Φ≤d (x i ) under the feature map Φ≤d . This
linear span is then related to the row-span of KXZ by applying several dualities in and on feature
space.
More technically, we first show how to identify R[t]≤d , the set of polynomials of degree
at most d, with linear functionals on the kernel feature space and that every f ∈ R[t]≤d can
be expressed as a kernel decision function. To understand the equations cutting out X (which
are not canonical) it is better to change perspective and consider the set of all polynomials
vanishing on X, called its vanishing ideal (which is canonical). The main results of this section
say, informally, that X can be identified with the intersection of a particular linear space L and
Φ≤d (Rn ), and that the vanishing ideal of X corresponds to L ⊥ , which in turn is isomorphic to the
row-span of KXZ . In other words, we relate a linear duality to the algebraic-geometric concept of
duality between sets of equations and their common vanishing locus.
2.1. Duality of polynomial rings and feature spaces
To begin with, we introduce inhomogeneous polynomial kernel as follows.
Definition 2.1. Let θ ∈ (0, 1) be a fixed real number, and let d ≥ 1. The inhomogeneous polynomial kernel function k≤d : Rn × Rn −→ R is given by k≤d (x, y) = (θ · 〈x, y〉 + 1)d , where 〈., .〉
denotes the standard scalar product.
This definition differs slightly from the usual one which is obtained after dividing by θ d .
Since θ is chosen arbitrarily in (0, 1), no qualitative change is introduced by our convention. It
is however, as we will see, the more natural one.
n+d 
It is well-known that feature space F≤d of k≤d satisfies F≤d ∼
= Rm , where m = d , and
that the feature map Φ≤d : Rn −→ F≤d is given by Φ≤d (x) = (γα x α | |α| ≤ d) (e.g., see [4],
α
α
Sec. 2.1). Here we let x α = c1 1 · · · cn n and |α| = α1 + · · · + αn for x = (c1 , . . . , cn ) and α =
q
d
|α| 
(α1 , . . . , αn ). Moreover, we have γα = θ |α| · |α| · α ,...,α . To define Φ≤d uniquely, we order
1

n

p1 , . . . , pm increasingly with respect to the degree-lexicographic term ordering. The feature map
is characterized by the property that k≤d (x, y) = 〈Φ≤d (x), Φ≤d ( y)〉 for x, y ∈ Rd , i.e., by the fact
that it transforms the kernel function to a standard scalar product on the feature space.
On the algebraic side, the main objects linked via duality are vector spaces of polynomials.
Let t = (t 1 , . . . , t n ) be a tuple of indeterminates, and let R[t] be the polynomial ring in these
indeterminates. For every d ≥ 0, we denote by R[t]≤d the vector space of all polynomials of

3

n+d 
degree at most d. Recall that the dimension of R[t]≤d is d . The inhomogeneous polynomial
kernel function can be extended to a pairing of polynomials as follows.

Notation 2.2. Let θ ∈ (0, 1) and d ≥ 1. The inhomogeneous polynomial kernel on R[t] is
the map k≤d : R[t]n × R[t]n −→ R[t] given by k≤d (F, G) = (θ · 〈F, G〉 + 1)d where we define
〈F, G〉 = f1 g1 + · · · + f n g n for F = ( f1 , . . . , f n ) and G = (g1 , . . . , g n ).
The polynomials k≤d (x, t) ∈ R[t] can be used to generate R(t]≤d , as our next proposition
shows.
By the multinomial theorem, they can be expressed as k≤d (x, t) = (θ 〈x, t〉 + 1)d =
P
2 α α
n
|α|≤d γα x t . Recall that points x 1 , . . . , x m ∈ R are called generic for a property P if there
exists a Zariski open subset1 of U ⊆ (Rn )m such that property P holds for all (x 1 , . . . , x m ) ∈ U.
n+d 
Proposition 2.3. Let d ≥ 0, let m0 ≥ d , and let x 1 , . . . , x m0 ∈ Rn be generic. Then we have
R[t]≤d = 〈k≤d (x i , t) | 1 ≤ i ≤ m0 〉.
The next proposition says that R[t]≤d is dual to the feature space F≤d .
Proposition
2.4. (a) For d ≥ 1, the map ϕ : R[t]≤d −→ (F≤d )∨ ∼
= HomR (Rm , R) defined by
P
α
∨
∨
f = |α|≤d cα t 7→ f , where f (eα ) = cα /γα , is an isomorphism of R-vector spaces.
P
P
(b) For d ≥ 1, the map ψ : R[t]≤d −→ F≤d defined by f = |α|≤d cα tα 7→ |α|≤d (cα /γα )eα is an
isomorphism of R-vector spaces.
By dualizing the map ϕ and using the canonical isomorphism between F≤d and its bidual,
we obtain an isomorphism ϕ ∨ : F≤d −→ (R[t]≤d )∨ which maps eα to the R-linear map given by
t β 7→ γαS
δαβ . Another consequence of the preceding proposition is that, if we pass to the union
R[t] = d≥1 R[t]≤d , we see that the polynomial ring contains the duals of all feature spaces
F≤d .
Finally, we interpret the isomorphism ϕ in terms of kernel decision functions. Recall that
the map eval( f ) : Rn −→ R given by x 7→ f (x) is called the polynomial function associated
to f ∈ R[t] and that the polynomial function associated to k≤d (x, t) is called a kernel decision
function. Let us denote the vector space of all kernel decision functions by KDF≤d .
Corollary 2.5. The map Φ∗≤d : (F≤d )∨ −→ KDF≤d given by ` 7→ ` ◦ Φ≤d is an isomorphism.
Now we use the isomorphism ϕ : R[t]≤d −→ (F≤d )∨ to transfer the standard scalar product
on F≤d to R[t]≤d in the natural way. The result is a scalar product 〈., .〉ϕ on R[t]≤d such that
〈tα , tβ 〉ϕ = 〈(1/γα ) eα∨ , (1/γβ ) eβ∨ 〉 = (1/γ2α ) δαβ
for |α| ≤ d. The next proposition provides a basic property of the scalar product 〈., .〉ϕ .
Proposition 2.6. For f ∈ R[t]≤d and x ∈ Rn , we have f (x) = 〈 f , k≤d (x, t)〉ϕ .
The duality expressed by the map ϕ is an algebraic analogue of the theory of reproducing
kernel Hilbert spaces. The associated Hilbert space is the space of polynomial functions f : Rn →
R of degree ≤ d which can be identified with R[t]≤d by replacing the polynomial function f with
the corresponding symbolic polynomial. The equation in the preceding proposition could also be
obtained by combining the Riesz representation with this identification (see [4], Sec. 2.2). In the
next section we go beyond what can be shown using the usual RKHS duality alone.
1

All such sets have full measure under any continuous probability density.

4

2.2. Duality of vanishing ideals and feature spans
Next we show that ideals – a classical concept in algebra – are the proper dual objects of feature
spans of manifolds, in the same way as the polynomial ring is the dual of feature space itself.
Recall that an ideal I in R[t] is a vector subspace such that I · R[t] ⊆ I. Ideals are connected to
subsets of Rn as follows.
Definition 2.7. (a) Given a subset X of Rn , the set of polynomials I(X) = { f ∈ R[t] | f (x) = 0
for all x ∈ X} is an ideal in R[t]. It is called the vanishing ideal of I.
(b) Given an ideal I in R[t], the set of points V(I) = {x ∈ Rn | f (x) = 0 for all f ∈ I} is called
the zero set of I. A subset X of Rn is called an algebraic set if it is the zero set of an ideal in R[t].
Notice that not every ideal in R[t] is a vanishing ideal, since vanishing ideals have the additional property of being radical, i.e., if f ∈ I i for some i ≥ 1 then f ∈ I. Given an ideal I in R[t]
and d ≥ 0, we let I≤d = I ∩ R[t]≤d .
From now on the data manifold X is always assumed to be an algebraic set in Rn . In this
case, the vanishing ideal I(X) is dual to the manifold X in the following sense.
Theorem 2.8. Let X ⊂ Rn be an algebraic set.
(a) We have I(X)≤d = 〈k≤d (x, t) | x ∈ X〉⊥ , where ⊥ is taken with respect to 〈., .〉ϕ .
(b) Under the isomorphism ϕ, the set I(X)≤d corresponds to {` ∈ (F≤d )∨ | `(〈Φ≤d (X)〉) = 0}.
(c) Under the isomorphism ψ, the set I(X)≤d corresponds to the feature span 〈Φ≤d (X)〉.
The space 〈Φ≤d (X)〉 will be called the feature span of X. As discussed above, Theorem 2.8 is
a kernelized version of the usual algebra-geometry duality based on Hilbert’s Nullstellensatz.
Remark 2.9. All the results in this section hold, mutatis mutandis, for homogeneous polynomial
kernels. For other kernels, they can be adapted via kernelizing the polynomial ring, but then
exact statements may need to be replaced with approximate ones. For this reason, and to save
space, we work with inhomogeneous polynomial kernels from now on.

3. Vanishing ideals and cross-kernel matrices
In this section we discuss structural and algebraic properties of cross-kernel matrices between
data points and further sampled points. As in the previous section we use inhomogeneous polynomial kernels. Similar statements can be proven for homogeneous polynomial kernels, and
even non-polynomial kernels. In the latter case our exact statements have to be transformed into
spectral approximation results.
The following setting is used throughout the section. There is a manifold X ⊆ Rn from which
data points are sampled, and a manifold Z ⊆ Rn from which further points are sampled in such
a way that they feature-span Z. We assume that both X and Z are algebraic sets and that there is
a degree d ≥ 1 such that X and Z are cut out by polynomials of degree at most d. Furthermore,
we assume X ⊆ Z. The data manifold X is considered to be fixed and unknown, while Z can
be chosen by the experimenter and is typically given by Z = Rn . The sampled data points are
denoted by x 1 , . . . , x n ∈ X, and the points samples from Z are denoted by z1 , . . . , z M ∈ Z. They
are the rows of the matrices X ∈ RN ×n and Z ∈ R M ×n , respectively. We begin by introducing
kernel and cross-kernel matrices.

5

Definition 3.1. (a) Given d ≥ 1 and the inhomogeneous polynomial kernel k≤d : Rn × Rn −→ R,
we denote the matrix of size N × M whose entry in position (i, j) is k≤d (x i , z j ) by K≤d (X , Z). This
matrix is called the cross-kernel matrix between X and Z.
(b) The matrix K≤d (X , X ) is simply called the kernel matrix of X .
In the following we study properties of the matrices K≤d (X , X ), K≤d (X , Z) and K≤d (Z, Z), and
we prove duality statements between cross-kernel matrices and vanishing ideals. An important
condition is genericity of the points sampled from Z. In our case, the appropriate genericity
condition is that the points feature-span Z. It is defined as follows.
Definition 3.2. Let d ≥ 1, let Φ≤d : Rn −→ F≤d be the feature map, and let z1 , . . . , z M ∈ Z.
(a) For a set S ⊆ Rn , the R-linear span fspan(S) = 〈Φ≤d (s) | s ∈ S〉 is called the feature span of S,
and the number frk S = dimR (fspan(S)) is called the feature rank of S.
(b) We say that z1 , . . . , z M ∈ Z feature-span the algebraic set Z if fspan({z1 , . . . , z M }) = fspan(Z).
By a slight abuse of notation, we also say that the matrix Z feature-spans Z in this case, and
we write frk(Z) for frk({z1 , . . . , z M }). Moreover, we say that z1 , . . . , z M are feature independent if
Φ≤d (z1 ), . . . , Φ≤d (z M ) are linearly independent in F≤d .
The matrix Z is also called the feature generating matrix. Notice that, by elementary linear
algebra, a set of points {z1 , . . . , z M } feature-spans Z if and only if frk(Z) = frk(Z). A traditional
method to analyse the kernel matrix K≤d (X , X ) is to subsample X , i.e., to choose the points zi in
{x 1 , . . . , x N }. This entails the common problem that one cannot assume that Z feature-spans X.
Our next theorem characterizes the matrices Z which avoid this problem.
Theorem 3.3. In the above setting, assume that fspan(X ) ⊆ fspan(Z).
(a) (Cross-Kernel Rank) We have rank K≤d (X , Z) = frk(X ). Hence, if X feature-spans X then
rank K≤d (X , Z) = frk(X).
(b) (Cross-Kernel Data Manifold Certificate) Let c ∈ Rn such that Φ≤d (c) ∈ fspan(Z). Then we
have Φ≤d (c) ∈ fspan(X ) if and only if (k≤d (c, z1 ), . . . , k≤d (c, z M )) is contained in the row span
of K≤d (X , Z). If X feature-spans X then this implies c ∈ X.
(c) (Cross-Kernel Nullspace) Let f ∈ I(Z)⊥
and suppose that X feature-spans X and Z feature≤d
PM
spans Z. Then we have f ∈ I(X)≤d if and only if f is of the form f = i=1 ci k≤d (zi , t) with a
coefficient tuple c = (c1 , . . . , c M ) ∈ R M such that K≤d (X , Z) · c > = 0.
(d) (Cross-Kernel Range) Let f ∈ I(Z)⊥
. Suppose that X feature-spans X, Z feature-spans Z, and
≤d
z1 , . . . , z M are feature-independent. Then f ∈ I(X)⊥
if and only if there is a vector (c1 , . . . , c M ) in
≤d
PM
−1/2
the left range of K≤d (X , Z) · K≤d (Z, Z)
such that f = j=1 c j k≤d (z j , t).
The statements of this theorem can be turned into statements about K≤d (X , X ) and K≤d (Z, Z)
by taking X = Z and X = Z. Let us denote the Moore-Penrose pseudoinverse of a matrix A by A+ .
The cross-kernel matrix allows us to reconstruct the kernel matrix K≤d (X , X ) as follows.
Theorem 3.4. In the above setting we have K≤d (X , X ) = K≤d (X , Z) · K≤d (Z, Z)+ · K≤d (Z, X ) if and
only if fspan(X ) ⊆ fspan(Z). In particular, equality holds if Z feature-spans Z.
The equality in Theorem 3.4 is an exact matrix decomposition of the kernel matrix of X which
does not require Z to be subsampled from X . This differs markedly from Nyström type methods.

6

4. The ideal-kernel duality and kernel learning
We demonstrate how kernel-ideal duality and cross-kernel matrices can be employed for common
learning tasks. They enable us to improve computational cost, stability, and to obtain novel
information about the data manifold. In the setting of the preceding section, we collect some
basic observations on the matrix K := K≤d (X , Z) · K≤d (Z, Z)−1/2 under the assumption that Z
feature-spans Z (or Rn ).
Remark 4.1. (a) Assuming that arithmetic and kernel function evaluation can be done at cost
O(1), it takes O(M 3 + M N ) steps to compute K and O(M 2 N + M 3 ) steps to obtain a singular
value decomposition of K. Thus both tasks require asymptotically linear time in N .
(b) By Theorem 3.4, the equality K≤d (X , X ) = K · K > holds exactly. Thus a SVD for K≤d (X , X ) can
be obtained from a SVD K = USV > of K via K≤d (X , X ) = US 2 U > . By (a), the cost of this method
is only O(M 2 N + M 3 ), which is linear in N , while the direct algorithm costs O(N 3 ).
(c) If Z is chosen to be a subsample of X , one recovers common subsampling strategies. However, the feature generating matrix Z can be chosen completely independently from X . Hence
the matrix K≤d (Z, Z)−1/2 needs to be precomputed only once for a given kernel degree and data
dimension. In particular, it is not necessary to subsample X . For instance, we can choose Z randomly. Or, if desired, the matrix Z can be chosen to contain both a subsample of X and random
points.
(d) As mentioned in Theorem 3.3.b, the right singular vectors of K characterize the data manifold X from which X is sampled. These right singular vectors cannot be obtained from K≤d (X , X )
alone.
The Ideal PCA (IPCA) algorithm
The kernel-ideal duality can be used to both speed up kernel PCA and additionally yield feature
functions cutting out the data manifold. By using the matrix K≤d (X , Z) instead of K≤d (X , X ),
we obtain the same feature projections and principal components as in ordinary kernel PCA
(both in terms of the X - and the Z-basis), the corresponding singular values, and novel feature
projections which cut out the data manifold (which we call “right principal components”). All of
this is achieved by Algorithm 1 in linear time and without the necessity of subsampling.
For dimension reduction tasks, the important quantities are evaluations of the principal components. The left principal components agree with the ones from PCA, while the orthogonal of
the space generated by the right principal components is used in Algorithm 2 to project onto the
data manifold.
The vector u contains exactly the principal feature evaluations of x one would obtain from
kernel PCA. The vector v is the IPCA-type right analogue, while v ⊥ can be seen as a measure for
how closely x lies on the data manifold. Neither v nor v ⊥ can be obtained from any algorithm
involving only the matrix K≤d (X , X ).

5. Experimental Validation
We validate the main claims of our paper experimentally. We show, in this section: (a) the
left principal IPCA features obtained from the cross-kernel are equal to the principal features
obtained from kernel PCA [5] and the kernel matrix KXX (5) but faster to obtain, (b) the vanishing
7

Algorithm 1 IPCA Computes left and right principal vectors of the cross-kernel.
Input: a degree d ≥ 1, a data matrix X ∈ RN ×n , a matrix Z ∈ R M ×n . A threshold " > 0 or a
cut-off m.
Output: m left principal components as columns of a matrix U ∈ RN ×m and m right principal
components as columns of a matrix V ∈ R M ×m , as well as the corresponding singular values on
the diagonal of a matrix S ∈ Rm×m .
Compute the matrices KXZ = K≤d (X , Z) ∈ RN ×M and KZZ = K≤d (Z, Z) ∈ R M ×M .
−1/2
2: Compute K = KXZ · KZZ .
1
3: If centering is desired, replace K by K − 1K, where 1 ∈ RN ×N denotes the all-ones matrix.
N
4: Compute the truncated singular value decomposition USV > of K, with U ∈ RN ×m , V ∈ R M ×m ,
and S = diag(σ1 , . . . , σm ). (If " is given, then σm is the smallest singular value with σm ≥ ".)

1:

5:

Output U, V, S.

Algorithm 2 eval-IPCA Computes IPCA features of the data.
Input: a degree d ≥ 1, a data matrix X ∈ RN ×n , a matrix Z ∈ R M ×n , a threshold " > 0, and a test
data point x t ∈ Rn .
Output: left principal features u ∈ Rm (PCA-like), right principal features v ∈ Rm (PCA-like) and
certifying features v ⊥ ∈ Rm (IPCA-like). The entries of u and v are local, whitened coordinates on
the data manifold X (u in X - and v in Z-basis), and the entries of v ⊥ are orthogonal coordinates
on X (i.e., they should be close to zero for x t ∈ X).
1:
2:
3:
4:
5:

Apply the IPCA algorithm to obtain matrices U, S, V . The columns of U ∈ RN ×m contain m
left principal vectors and the columns of V ∈ R M ×m contain m right principal vectors.
Compute the evaluation vectors κX = K(x t , X ) and κ Z = K(x t , Z).
−1/2
Compute (or pass over from IPCA) the matrix KZZ = K≤d (Z, Z)−1/2 .
Compute P = I − V V > , where I is the identity matrix of size M .
−1/2
Output: u = κX · U · S −1 v = κ Z · V · S −1 v ⊥ = κ Z · KZZ · P

features obtained from IPCA can be used for manifold learning with kernels, and (c) the IPCA
vanishing features outperform kernel PCA features in the USPS classification scenario in terms of
compactness, computational cost, and accuracy.
Data sets. Our synthetic data below are: (i) points uniformly from a circle of radius 10, then
perturbed by independent 2-dimensional standard, Gaussian random vectors, and (ii) 2 circles
on a sphere of radius 5. For each circle, 200 points are sampled uniformly, then 3-dimensional
Gaussian noise with variance 0.1 is added. The real world data set used below is (iii) the USPS
handwritten digits data set [8]. All experiments below are done with the inhomogenous polynomial kernel of degree 2, with θ = 1.
Kernel PCA vs IPCA. We compared kernel PCA features and IPCA left features on the (ii)
2 circles data. We fixed M = 12 as the number of feature spanning points, and computed
principal components for N = 10, 20, . . . , 1000 data points. The six non-degenerate principal
vectors and their principal values were equal to the IPCA left feature vectors and the squares
of the corresponding singular values, up to machine precision. Figure 1 shows log-runtimes in
comparison. The runtime for IPCA is lower by orders of magnitude.

8

2

0
IPCA
kPCA

log(time/sec)

−2

−4

−6

−8

−10

0

100

200

300

400
500
600
number of data points (N)

700

800

900

1000

Figure 1: Runtime of kernel PCA and IPCA. Thick line is mean, dotted lines are 0.1 and 0.9
quantiles over 20 repetitions.
Manifold learning with cross-kernel and IPCA. We used IPCA to learn the data manifold
from 200 samples in the (i) circle and the (ii) 2 circle example. The manifold is the approximated
by the points with low evaluation of the normalized IPCA certifying features, the result can be
seen in figure 2. The singular value threshold was fixed to m = 5 for (i) and m = 9 for (ii).
The certifying features were evaluated on a fine grid, grid points with evaluation norm near the
training data were considered in-manifold, with 0.1-quantile for (i) and 0.01-quantile for (ii).
6

6

4

4

2

2

0

0

−2

−2

−4

−4

−6

−6

10

10
5

10
5

0

5

10
5

0

0
−5

−5

−5
−10

(a)

0
−5
−10

−10

−10

(b)

Figure 2: Manifold estimation with IPCA. The crosses are the data points. The estimated manifolds are the light green ring in 2(a) and the two blue rings in 2(b).
Classification with IPCA certifying vectors. We used the IPCA features for classification on
the USPS data set. In all experiments, 700 training examples are used, and the rest of the USPS
data set is used for evaluation.
t
Denote the different training classes X 1 , . . . , X t and the pooled data by X train := ∪i=1
X i . The
natural form of a one-vs-all classifier using IPCA certifying features is as follows: (1) to train,
compute IPCA(X i , Z) for each of the X i with the same fixed Z; to classify a set of points X test ,
compute the right-evaluation vectors using Algorithm 2 for each class and assign the label corresponding to the minimum norm evaluation. The correctness of this construction follows directly
from Theorem 3.3. All the stability and convergence guarantees of IPCA go through without
modification.
We compare three choices for the feature generating points Z: (1) choose the featuregenerating matrix Z with standard normal Gaussian entries; (2) choose it as a uniform random sub-sample of X train ; (3) choose it as a degenerate sub-sample of X train ; we simulate a high
9

number of repetitions by sub-sampling from a sub-sample of size N /4.
Figure 3(a) summarizes the results. Non-degenerate sub-sampling (2) is best, since it takes
advantage of the implicit dimension-reduction of expressing the Xi in a coordinate system derived X ⊂ Rn , which can be of much lower dimension than n. However, when the sampling is
degenerate as in (3), e.g., with many repetitions, a random feature-generating matrix of the same
size as in (1) gives better performance. Eventually, with growing M all choices of Z yield comparable results. As it can not be distinguished a-priori whether the training sample is degenerate
or not, random Z can be preferrable.
We also compared the IPCA certifying features to the IPCA left and right principal features
(PCA-like), by extracting m = 32 features, then performing linear classification via LIBLINEAR
[3]. The feature-generating matrix Z is chosen randomly with standard normal entries. Figure
3(b)–(c) show that m = 32 certifying features yields more accurate and faster results than the
same number of principal features. The results reported in [5] for classical kernel PCA are
comparable to the PCA-like features, which all are outperformed by the IPCA certifying features.
1

12

1

0.9

IPCA
PCA−left
PCA−right

0.95

subsample
random
repetition

0.8

10
0.9

IPCA
PCA−left
PCA−right

0.85

0.7

8

0.5

training time

accuracy

accuracy

0.8

0.6

0.75

6

0.7

0.4

4

0.65

0.3

0.6

0.2

0.55

2

0.1
50

100

150

200
250
feature geneating points (M)

(a)

300

350

400

0.5
50

100

150

200
250
feature generating points (M)

(b)

300

350

400

0
50

100

150

200
250
feature generating points (M)

300

350

400

(c)

Figure 3: (a) Classification with IPCA certifying features, comparing random Z and subsampled
Z from generic and degenerate sampling; (b) performance comparison of the three IPCA features
- left and right principal features are closely related to kernel PCA features; (c) speed comparison
of IPCA feature based classifiers.
Above comparing different modes of instanciating IPCA, the experiments also show that IPCA
is capable of extracting competitive features for considerably less and better scaling computational cost than kernel PCA, and achieves comparable performance on a much smaller number
of relevant features.
Acknowledgments LT is supported by the European Research Council under the European
Union’s Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement no 247029SDModels. This research was conducted while MK and LT were guests of FK at Mathematisches
Forschungsinstitut Oberwolfach, supported by FK’s Oberwolfach Leibniz Fellowship.
References
[1] Mark A. Aizerman, Emmanuel M. Braverman, and Lev I. Rozonoer. Theoretical foundations of the potential function method in pattern recognition learning. In Automation and Remote Control,, number 25
in Automation and Remote Control„ pages 821–837, 1964.

10

[2] Bernhard E. Boser, Isabelle M. Guyon, and Vladimir N. Vapnik. A training algorithm for optimal
margin classifiers. In Proceedings of the 5th Annual ACM Workshop on Computational Learning Theory,
pages 144–152. ACM Press, 1992.
[3] R.-E. Fan, K.-W. Chang, C.-J. Hsieh, X.-R. Wang, and C.-J. Lin. LIBLINEAR: A library for large linear
classification. Journal of Machine Learning Research 9:1871–1874, 2008.
[4] Bernhard Schölkopf and Alexander J Smola. Learning with kernels. MIT Press, 2002.
[5] Bernhard Schölkopf, Alexander Smola, and Klaus-Robert Müller. Nonlinear component analysis as a
kernel eigenvalue problem. Neural computation, 10(5):1299–1319, 1998.
[6] John Shawe-Taylor and Nello Cristianini. Kernel Methods for Pattern Analysis. Cambridge University
Press, New York, 2004.
[7] Vladimir N. Vapnik. The Nature of Statistical Learning Theory. Springer Verlag, New York, 1995.
[8] USPS handwritten digits data set A Database for Handwritten Text Recognition Research, J. J. Hull,
IEEE PAMI 16(5) 550-554, 1994.

11

Appendix
Proof: (of Proposition 2.3)
P
n+d 
Let m = d . Under the isomorphism R[t]≤d → Rm given by α cα tα 7→ (cα )|α|≤d , the polynomial k≤d (x, t) maps to (γ2α x α )|α|≤d for every x ∈ Rn . Up to scaling the coordinates by γ2α , the
image of this isomorphism is a Veronese variety. Now the claim follows from the fact that this
Veronese variety is irreducible and non-degenerate.
Proof: (of Proposition 2.4)
First we show (a). Both vector spaces have the same dimension. The map ϕ is R-linear and sends
the basis vector tα to the basis vector (1/γα ) · eα∗ , i.e., to the map defined by eβ 7→ (1/γα ) δαβ .
From this the claim follows.
To prove (b), we use the standard scalar product on F≤d to identify it with (F≤d )∨ . This (noncanonical) isomorphism identifies the standard basis vector eα with eα∗ . Thus the claim follows
from (a).
Proof: (of Corollary 2.5)
It is easy to check that Φ∗≤d = eval ◦ϕ −1 . Since both maps eval and ϕ −1 are isomorphisms, the
claim follows.
Proof: (of Proposition 2.6) P
This follows by writing f = α cα tα with cα ∈ R and computing
P
P
P
P
〈 f , k≤d (x, t)〉ϕ = 〈 α cα tα , β γ2β x β tβ 〉ϕ = α cα γ2α x α 〈tα , tα 〉ϕ = α cα x α = f (x)

Proof: (of Theorem 2.8)
Part (a) follows from the definition of I(X) and Prop. 2.6.
To prove (b), we first use (a) in order
X)≤d isomorphically into (F≤d )∨ . The result
P to map I(
is the orthogonal complement of U = 〈 |α|≤d γα x α eα∨ | x ∈ X〉. Since we are using the standard
scalar products on F≤d and its dual, this orthogonal complement equals
P the annihilator of the
preimage of U in F≤d . That preimage is generated by the elements |α|≤d γα x α eα = Φ≤d (x)
with x ∈ X. In other words, the vector space U ⊥ is the annihilator of the feature span of X, as
claimed.
Finally, we note that (c) follows from (b) by identifying F≤d with its dual via the standard
scalar product.
Proof: (of Theorem 3.3)
First we prove claim (a). Letting vi = Φ≤d (x i ) ∈ Rm and w j = Φ≤d (z j ) ∈ Rm for i = 1, . . . , N
and j = 1, . . . , M , we have to show that the rank of the matrix G = (〈vi , w j 〉)i, j equals the
dimension of 〈v1 , . . . , vN 〉. By removing columns and rows of G, we may assume that {v1 , . . . , vN }
and {w1 , . . . , w M } are linearly independent. Using the hypothesis 〈v1 , . . . , vN 〉 ⊆ 〈w1 , . . . , w M 〉, we
see that M ≥ N and may assume vi = w i for i = 1, . . . , N . Then the first N columns of G are the
Gram matrix of V = (v1 , . . . , vN ) and the claim follows from the fact that the rank of the Gram
matrix of V equals the rank of V .
Next we show (b). The implication “⇐” is trivially true. To prove the converse,
we note
PN
that the hypothesis implies that there are a1 , . . . , aN ∈ R such that Φ≤d (c) − i=1 ai Φ≤d (x i )
is orthogonal to fspan(Z). Since this vector is contained in fspan(Z), it is zero. To prove the
additional claim, we note that Φ≤d (c) ∈ fspan(X ) implies that all linear forms in the annihilator
12

of fspan(X ) vanish at c. Via ϕ −1 , these linear forms correspond to the polynomials in I(X)≤d .
Since X is cut out by these polynomials, we get c ∈ X, as claimed.
For the proof of part (c) we use the isomorphism ϕ to write the hypothesis f ∈ I(Z)⊥
as
≤d
PM
f = j=1 c j k≤d (z j , t) for some c j ∈ R. We have f ∈ I(X) if and only if ϕ( f ) is orthogonal to
〈ϕ(k≤d (x i , t))〉. Therefore we require that we have
PM P
P
0 = 〈ϕ( f ), ϕ(k≤d (x i , t))〉 = 〈 j=1 c j α γα z αj eα∨ , α γα x iα eα∨ 〉
P

P
M
α
c
z
x iα
= α γ2α
j=1 j j
for i = 1, . . . , N . On the other hand, the entry in position (i, j) of K≤d (X , Z) is k≤d (x i , z j ) =
P
P
P
〈Φ≤d (x i ), Φ≤d (z j )〉 = 〈 α γα x iα eα , α γα z αj eα 〉 = α γ2α x iα z αj . Hence (c1 , . . . , c M )> is in the
PM P
nullspace of K≤d (X , Z) if and only if j=1 α γ2α x iα z αj c j = 0 for i = 1, . . . , N . As we have seen,
this is equivalent to f ∈ I(X).
Finally we show (d). Using Theorem 2.8, we see that the hypothesis yields I(X)⊥ = 〈k≤d (x i , t) |
i = 1, . . . , N 〉 and I(Z)⊥ = 〈k≤d (z j , t) | j = 1, . . . , M 〉. Therefore there exist numbers c j ∈ R such
PM
that we have f = j=1 c j k≤d (z j , t), and the question is whether there exists a representation
PN
f = i=1 ai k≤d (x i , t) with ai ∈ R. By the hypothesis, the matrix K≤d (Z, Z) is invertible, since
it is the Gram matrix of the full rank matrix (Φ≤d (z j )). So, knowing (c1 , . . . , c M ) is equivalent to
knowing
 
P 2 P M
α
(c1 , . . . , c M ) · K≤d (Z, Z) =
γ
c
z
z αj j
α α
`=1 ` `
To prove the implication “⇒” we can write the latter tuple as
 
P 2 PN
α
γ
z αj j = (a1 , . . . , an ) · K≤d (X , Z)
a
x
α α
i=1 i i
and conclude that c = (a1 , . . . , an ) · K≤d (X , Z) · K≤d (Z, Z)−1 . For the implication “⇐”, we are
given the equalities


P 2 P M
P 2 PN
α
α
α
γ
c
z
z
=
γ
a
x
z αj
α α
α α
`=1 ` `
j
i=1 i i
PM
PN
for j = 1, . . . , M . This means that, for j = 1, . . . , M , the vectors `=1 γα c` z`α eα − i=1 γα ai x iα eα
PM
and j=1 γα z αj eα are orthogonal with respect to the standard scalar product in F≤d . Thus the first
vector is
orthogonal
Pboth
PN to fspan(Z) and contained in fspan(Z), i.e., it is zero. Consequently,
M
we get `=1 c` z`α = i=1 ai x iα , and as we have seen, this yields f ∈ I(X)⊥
, as claimed.
≤d
Proof: (of Theorem 3.4)
Let Y ∈ R(N +M )×n be the matrix obtained by row concatenating X and Z. Observe that, by
construction, there are projection matrices P and Q such that K(X , X ) = P · K(Y, Y ) · P > and
K(X , Z) = P · K(Y, Y ) · Q> , as well as K(Z, Z) = Q · K(Y, Y ) · Q> .
To prove the implication “⇒”, we suppose that fspan(X ) * fspan(Z). By the rank-nullity theorem, this implies frk(Y ) > frk(Z). Thus Theorem 3.3.a yields rank K≤d (Z, Z) < rank K≤d (Y, Y ).
The claim then follows by applying Lemma 5.1 to A = K≤d (Y, Y ) and B = Q and C = Q> .
For the reverse implication, we first note that fspan(X ) ⊆ fspan(Z) implies fspan(Y ) =
fspan(Z). Thus we have frk(Y ) = frk(Z) and Theorem 3.3.a shows rank K≤d (Z, Z) = rank K≤d (Y, Y ).
Now the claim follows by applying Lemma 5.1 in the same way as before.
13

The following Lemma provides the key ingredient for the preceding proof.
Lemma 5.1. Let A ∈ Rm×n , let B ∈ Rn×k , and let C ∈ R`×m . Then the matrix equality A =
AB(CAB)+ CA holds if and only if rank A = rank CAB.
Proof: First we show the implication “⇐”. Since rank CAB = rank A, we have rank CA = rank AB =
rank A, and therefore colspan(CA) = colspan(A) as well as rowspan(AB) = rowspan(A). This implies A = C + CA = ABB + . The definition of the pseudoinverse yields CAB = CAB(CAB)+ CAB.
Multiplying this equality by C + and B + , we get A = C + CABB + = C + CAB(CAB)+ CABB + =
AB(CAB)+ CA, as claimed.
To prove the reverse implication, suppose that rank A 6= rank CAB. Since the rank of a matrix
does not increase by matrix multiplication, we must have rank A > rank CAB. Since the pseudoinverse of CAB satisfies rank(CAB)+ = rank CAB, we find rank A > rank CAB = rank(CAB)+ ≥

rank AB(CAB)+ CA , in contradiction to the hypothesis.

14

