Tangent Hyperplane Kernel Principal Component Analysis for Denoising

13 Pages·2012·1.53 MB·English

644 IEEETRANSACTIONSONNEURALNETWORKSANDLEARNINGSYSTEMS,VOL.23,NO.4,APRIL2012
Tangent Hyperplane Kernel Principal
Component Analysis for Denoising
Joon-Ku Im, Daniel W. Apley, and George C. Runger
Abstract—Kernel principal component analysis (KPCA) is subspace is defined as the span of the dominant eigenvectors
a method widely used for denoising multivariate data. Using oftheM×M covariancematrixof{φ(x ),φ(x ),...,φ(x )}.
1 2 N
geometric arguments, we investigate why a projection operation The projected feature vector, which we denote by Pφ(x),
inherent to all existing KPCA denoising algorithms can some-
is a point in the feature space. In order to denoise x, we
times cause very poor denoising. Based on this, we propose
a modification to the projection operation that remedies this must find a point xˆ in the input space that best corresponds
problem and can be incorporated into any of the existing to Pφ(x). The exact preimage of Pφ(x) (under the map φ)
KPCA algorithms. Using toy examples and real datasets, we generally does not exist. That is, there generally exists no xˆ
show that the proposed algorithm can substantially improve suchthatφ(xˆ)= Pφ(x)exactly.Thisistheso-calledpreimage
denoising performance and is more robust to misspecification
problemin KPCA denoising[3].By preimageapproximation,
of an important tuning parameter.
we mean that one must find a value xˆ that minimizes some
Index Terms—Denoising, kernel, kernel principal component appropriately defined measure of dissimilarity between xˆ [or
analysis (KPCA), preimage problem.
φ(xˆ)] and Pφ(x).
AlloftheaforementionedKPCAdenoisingmethodsemploy
I. INTRODUCTION the same orthogonal projection operation. Their distinctions,
KERNEL principal component analysis (KPCA) is a which we discuss in more detail in Section II, lie in how
nonlinear generalization of linear principal component they perform the preimage approximation. For example, the
analysis (PCA) and is widely used for denoising multivariate original method of Schölkopf et al. [3] finds the point xˆ that
datatoidentifyunderlyingpatternsandforotherpurposes[1]. minimizes a dissimilarity measure defined as the Euclidean
Let x be a multivariate observation vector having n compo- distance between φ(xˆ) and Pφ(x). That is, they take, xˆ =
nents,andsupposethedataconsistof N suchobservations{x1, argminz||φ(z)−Pφ(x)||2 as the denoised x, where the norm
x2,...,xN}.InKPCA,onedefinesasetof M features{φ1(x), is the standard Euclidean 2-norm.Fig. 2 illustrates the results
φ2(x), ..., φM(x)}, each of which is some appropriately of this method for a toy example with N = 100 observations
chosen nonlinear function of x, and forms a feature vector of n = 2 variables, using a polynomial kernel of degree 2.
φ(x) = [φ1(x), φ2(x), ..., φM(x)](cid:2). One then conducts PCA Fig. 2(a) shows the original data, and Fig. 2(b) shows the
on the set of feature vectors {φ(x1), φ(x2), ..., φ(xN)}. If denoised data using the method of Schölkopf et al. [3]. The
we considerRn the “inputspace” in which x lies and RM the salientfeatureofFig.2(b)isthatthepointsaroundthemiddle
“feature space” in which φ(x) lies, we can view the feature ofthecurve(e.g.,x )arenotdenoisednearlyaswellasthose
A
map as φ: Rn →RM. In practice, the feature map is usually near the ends of the curve (e.g., x ).
B
defined implicitly by some kernel function for computational Inthispaper,weinvestigatethisphenomenonandshowthat
reasons [2], [3]. it is a consequence of the projection operation, as opposed to
Several different KPCA denoising methods have been pro- the preimage approximation. Consequently, this undesirable
posed[2]–[9],which we reviewandcontrastmorethoroughly characteristic (much poorer denoising in some regions of
inSectionII.Theircommonalityisthattheyallshareageneral the input space than in other regions) is inherent to all the
framework that involves two primary operations, which we aforementioned KPCA algorithms, and the solution to this
refertoasprojectionandpreimageapproximation.By projec- problem that we propose can be used to enhance any of the
tion, we mean that, to denoise an observation x via KPCA, these algorithms.
oneprojectsitsfeaturevectorφ(x)ontotheprincipalsubspace An intuitive explanation for the poor denoising in some
in the feature space, as illustrated in Fig. 1. The principal regions is as follows. Consider the image of Rn under the
map φ, which is the subset of the feature space for which
ManuscriptreceivedDecember2,2010;acceptedDecember28,2011.Date
each point has an exact preimage. This subset having exact
of publication February 21, 2012; date of current version March 6, 2012.
This work was supported in part by the National Science Foundation under preimages forms a manifold in RM, which we refer to as the
GrantCMMI-0826081 andGrantCMMI-0825331. full manifold. As operations in RM, we can view projection
J.-K.ImandD.W.ApleyarewiththeDepartmentofIndustrialEngineering ontotheprincipalsubspaceastakingthepointφ(x)awayfrom
andManagementSciences,NorthwesternUniversity,Evanston,IL60208USA
(e-mail: [email protected]; [email protected]). thefullmanifoldtoobtainPφ(x),andpreimageapproximation
G.C.RungeriswiththeSchoolofComputing, Informatics, andDecision as taking the projection Pφ(x) back to some point on the full
Systems Engineering, Arizona State University, Tempe, AZ 85287 USA manifold that is as similar to Pφ(x) as possible, to obtain
(e-mail: [email protected]).
Digital ObjectIdentifier 10.1109/TNNLS.2012.2185950 φ(xˆ), or equivalently xˆ. As we will show, how well a point
2162–237X/$31.00 ©2012IEEE
IMetal.: TANGENTHYPERPLANEKERNELPCAFORDENOISING 645
and the covariance matrix of the feature vectors as
(cid:3)N
S(cid:3) = N−1 φ˜(xj)φ˜(cid:2)(xj)= N−1(cid:3)˜(cid:2)(cid:3)˜.
j=1
KPCA therefore involves the eigenvalues and eigenvectors of
the M ×M matrix S(cid:3).
Given that direct calculation of the eigenvectors of S(cid:3) is
usuallycomputationallyinfeasible(M isoftenlarge),thenovel
enabling aspect of KPCA is the kernel trick, by which we
can solve this eigenproblem without actually computing the
Ffeiagt.u1re. vDecetnoorisφin(xg)aonntoobstehrevaptriionncixpailnvsoulbvsepsactweoanodpepraretiiomnasg:epraopjpecrotixoinmoaftioitns high-dimensionalfeaturesor their covariancematrix S(cid:3). This
tofindanxˆ thatminimizessomemeasureofdissimilaritybetweenxˆ(orφ(xˆ)) is possible because in KPCA the feature map φ is defined
and Pφ (x). implicitly via its inner products (cid:4)φ(x),φ(y)(cid:5) = φ(cid:2)(x)φ(y) =
K(x,y) for some appropriately chosen positive definite kernel
function K(•,•) (see [2], [3]; for comprehensivetreatments of
x is denoised is strongly influenced by how far from the full kernel functions, see [13]), and the inner products are all that
manifold its projection Pφ(x) falls. For points such as x in are required in the computations. Common choices for the
A
Fig. 2, the projection step takes them quite far away from the kernel are the polynomial kernel K(x,y) = (1 + x(cid:2)y)d of
fullmanifold.Hence, Pφ(x )mustbemovedalargedistance degreed, andthe Gaussianradialbasisfunction(RBF) kernel
A
togetbacktothefullmanifoldinthepreimageapproximation K(x,y) = exp(–||x − y||2/ρ) where the parameters d and ρ
step, which results in poor denoising. On the other hand, for are chosenbythe users. A polynomialkernelof degree2 was
pointssuch as x , the projectionstep keepsthem close to the used in the Fig. 2 example.
B
full manifold, which results in more effective denoising. More specifically, define the centered kernel function
We propose a modification to the projection step which K˜(x,y) = (cid:4)φ˜(x),φ˜(y)(cid:5)=(cid:4)φ(x)−φ¯,φ(y)−φ¯(cid:5)
is designed to keep the “projected” feature vector closer to
= (cid:4)φ(x),φ(y)(cid:5)−(cid:4)φ(x),φ¯(cid:5)−(cid:4)φ(y),φ¯(cid:5)+(cid:4)φ¯,φ¯(cid:5)
the full manifold, thereby improving denoising performance,
especially at points such as x in Fig. 2, which are more (cid:3)N (cid:3)N
difficult to denoise. We refer toA the approach as the tangent = K(x,y)−N−1 K(x,xi)−N−1 K(y,xi)
hyperplane KPCA algorithm for reasons that will become i=1 i=1
(cid:3)N (cid:3)N
apparent later. The modified projection can be used in con-
junction with the preimage approximation algorithms of any +N−2 K(xi,xj)
of the aforementioned KPCA approaches. i=1 j=1
The remainder of this paper is organized as follows. = K(x,y)−N−11(cid:2)Kx−N−11(cid:2)Ky+N−21(cid:2)K1
SectionIIbrieflyreviewsthegeneralKPCAapproachandvar-
where 1 denotes a column vector of ones, and K =
x
ious preimage approximation algorithms. Section III provides [K(x,x ), K(x,x ),...,K(x,x )](cid:2). Also define the centered
a geometric explanation of why certain points such as xA in kernel 1matrix K˜2= (cid:3)˜(cid:3)˜(cid:2), whoNse (i, j)th element is K˜ =
ij
Fig. 2 are denoisedpoorlyusing the regularprojectionstep in K˜(x ,x ). Note that K˜ is an N × N matrix that can be
i j
KPCA. In Section IV, we introduce our tangent hyperplane
calculatedvia evaluationsof the kernelfunction.Thisis com-
KPCA approach as a solution to this problem. Section V
putationally significant, because the eigenvector/eigenvalue
presents denoising results for the toy example of Fig. 2 and pairs {(vk, λk): k = 1, 2, ..., L0} of S(cid:3) with nonzero
forrealimagedatasets,togetherwithadiscussionofthechar-
eigenvalue (we refer to such an eigenvector as a nonzero
acteristics of the approach. Section VI concludes this paper.
eigenvector, and assume there exist L of them) and the
0
nonzero eigenvector/eigenvalue pairs {(αk = [αk,1, αk,2, …,
αk,N](cid:2), ηk): k = 1,2,...,L0} o(cid:2)f K˜ are closely related (see
II. REVIEWOF PRIOR KPCA WORK AND THE [1] for details): vk = (cid:3)˜(cid:2)αk = iN=1(cid:3)˜ (xi)αk,i, and ηk =
PREIMAGE PROBLEM Nλ , where the convention for scaling the eigenvectors is to
k
take ||v || = 1, which implies ||α || = (Nλ )−1/2 = η−1/2.
Regular(linear)PCAonthen-dimensionaldata{x1,x2,…, Using thkis, one can solve the eigeknproblemkfor the M ×k M
x }findsthesetoforthonormalbasisvectorsalongwhichthe
N matrixS(cid:3) bysolvingthemuchlesscomputationallyexpensive
componentsof xhavelargestvariance.Thesebasisvectorsare eigenproblem for the N ×N matrix K˜ (N << M typically).
theeigenvectorsofthen×n (sample)covariancematrixofthe
Asdiscussedintheintroduction,theexistingKPCAdenois-
data [10]–[12]. KPCA is PCA on the M-dimensional feature
ing methods share a common projection operation. That is,
vectors, {φ(x ), φ(x ), …, φ(x )}, where M >> n usually.
φD¯e=finNe−th1e(cid:2)ceNjn=1t1eφre(dxjf2)e,attuhreefveeactutNorer amsaφt˜r(ixx)a=s φ(x)−φ¯ where pgrivoejenctainonobservation x, the(cid:3)yL find (implicitly) the orthogonal
P(cid:3)˜(x)= (cid:4)φ˜(x),v (cid:5)v
(cid:3)˜ =[φ˜(x1) φ˜(x2) ··· φ˜(xN)](cid:2) k=1 k k
646 IEEETRANSACTIONSONNEURALNETWORKSANDLEARNINGSYSTEMS,VOL.23,NO.4,APRIL2012
Fig.2. (a)Original data(whitemarkers)distributed approximately alongaquadratic curve(solidcurve). (b)Denoiseddata(whitemarkers)obtained using
themethodof[3].PointssuchasxA nearthemiddleofthecurvearedenoised lesseffectively thanpoints suchasxB neartheendsofthecurve.
ofφ˜(x)ontotheprincipalsubspacespannedbythe L principal where K˜x =[K˜(x,x1), K˜(x,x2), ··· , K˜(x,xN)](cid:2). Thus, the
eigenvectors {v1, v2, ..., vL} of S(cid:3). This is followed by the quantity to be minimized (with respect to z, for a given x
preimage approximationfor finding the xˆ that is most similar to be denoised) in the preimage approximation of [3] can be
to Pφ˜(x), which is where the KPCA denoising methods written as
mentioned in the introduction differ. Specifically, they differ ||φ˜(z)− Pφ˜(x)||2
in how they define the dissimilarity measure and/or how they
= (cid:4)φ˜(z),φ˜(z)(cid:5)−2(cid:4)φ˜(z),Pφ˜(x)(cid:5)+(cid:4)Pφ˜(x),Pφ˜(x)(cid:5)
implementtheoptimizationalgorithmtominimizeit,including
the variables over which they optimize. = K˜(z,z)−2β(cid:2)β +β(cid:2)β
z x x x
In [2] and [3], the authors suggested minimization of the
where
Euclidean distance (cid:6) (cid:7) (cid:6) (cid:7)
||φ˜(z)− Pφ˜(x)||2 (1) βx = v1 v2 ··· vL (cid:2)φ˜ (x)= α1 α2 ··· αL (cid:2)K˜x
is the vector of principal components scores of φ˜(x).
overz.Inotherwords,thedissimilaritymeasurebetweenxˆ (or
φ˜(xˆ)) and Pφ˜(x) is the Euclidean distance between φ˜(xˆ) and In [2] and [3], the authors suggested using a fixed-point
Pφ˜(x) in the feature space. The computational kernel trick iterative algorithm, related to a gradient descent algorithm,
to minimize (1) with respect to z over the entire input space
that is central to all KPCA methods is based on the fact that
Rn.TakahashiandKurita[4]modifiedthefixed-pointiterative
for any x ∈Rn, z ∈Rn, and nonzero eigenvector v , we can
k
algorithm by updating certain weighting coefficients at each
evaluate the following inner products via the kernel:
(cid:4) (cid:5) iteration, based on updatingthe feature vector to be projected
(cid:3)N ontotheprincipalsubspace.Afixedpointoftheiralgorithmis
(cid:4)φ˜(x),vk(cid:5) = φ˜(x), φ˜(xi)αk,i more consistent with minimizing ||φ˜(z)− Pφ˜(z)||2. Teixeira
i=1 et al. [5] proposed an alternative initial guess for the fixed-
(cid:3)N
= (cid:4)φ˜(x),φ˜(xi)(cid:5)αk,i =K˜(cid:2)xαk poWinthiitleerathtievepraelcgeodriinthgmm.ethodssharea commondissimilarity
i=1 measure(1),othersusedifferentdissimilaritymeasures.Kwok
and (cid:4) (cid:5) and Tsang [6] proposed a method akin to multidimensional
(cid:3)L scaling. To denoise x, they found xˆ so that the distances
(cid:4)φ˜(z),Pφ˜(x)(cid:5) = φ˜(z), (cid:4)φ˜(x),vk(cid:5)vk ||xi − xˆ|| (i = 1,2,...,n) are as close as possible to
k=1 the distances ||φ˜(x ) − Pφ˜(x)||2 (i = 1,2,...,n). Their
i
(cid:3)L dissimilarity measure is the distance between the vector of
= (cid:4)φ˜(z),vk(cid:5)(cid:4)φ˜(x),vk(cid:5) ||xi − xˆ||2’s and the vector of ||φ˜(xi) − Pφ˜(x)||2’s for
k=1 values of i such that φ˜(x ) is a neighbor of Pφ˜(x), and
i
= (cid:3)L (K˜(cid:2)α )(K˜(cid:2)α ) to reduce computational expense they restricted xˆ to being
z k x k a linear combination of the nearest neighbors. Bakır et al.
k=1 [7] suggested a kernel regression method to estimate the
= K˜(cid:2)[α ··· α ][α ··· α ](cid:2)K˜ map (cid:8)(Pφ˜(x)) = [(cid:8) (Pφ˜(x)),(cid:8) (Pφ˜(x)),...,(cid:8) (Pφ˜(x))](cid:2),
z 1 L 1 L x 1 2 n
IMetal.: TANGENTHYPERPLANEKERNELPCAFORDENOISING 647
(a) (b)
Fig.3. (a)Datawithnoiseintheinputspace.Thesolidcurveisthetruepatternthatthedatawouldfollowwithnonoise.(b)Centeredfeaturevectorsand
truepatternmappedtothefeature space.
(a) (b)
Fig. 4. (a) Result of [3] in the input space. Approximate preimages (black markers). (b) Result of [3] in the feature space. Projections (gray marker) of
the features of the noisy data (white marker) using the simple 3-D feature and the method of [3]. The approximate preimages in (a) are supposed to best
correspondtotheprojections in(b).Wecanobservethat xB isdenoisedmuchmorethanxA.
(cid:2)
where one views the components of Pφ˜(x) as the predictor Hence, they used a dissimilarity measure ||φ˜( N w x )−
i=1 i i
variables and the components of xˆ as the response variables. Pφ˜(x)||2,andchosetheweightstominimizethisundercertain
T(cid:2)hey assume the kernel regression model (cid:8)j(Pφ˜(x)) = constraints.
N βjκ(Pφ˜(x),Pφ˜(x )), where κ is a different kernel Nguyen and De la Torre [9] approached the problem in a
opei=r(cid:2)a1tori, and minimize tihe sum of squared regression residu- different way. They proposed to solve xˆ = argminz||φ˜ (z)−
als N ||x −(cid:8)(Pφ˜(x ))||2 with respect to the βj(cid:2)s. Their Pφ˜ (z)||2+C ·d(z,x) where d(z,x) is any distance measure
i=1 i i i
approach applies mainly to the situation in which noiseless between z and x, and C is a weight that is chosen to balance
observations are available for training. Zheng et al. [8] noted betweenthetwoterms.Althoughthismethodslightlydeviates
that a denoised observation of x obtained by most of the from the general framework we have portrayed, it is similar
preceding metho(cid:2)ds is a weighted sum of {x1, x2, ..., xN}. in the sense thatitstill involvesanorthogonalprojectiononto
That is, xˆ = N w x for some weights w , ..., w . the principal subspace.
i=1 i i 1 N
648 IEEETRANSACTIONSONNEURALNETWORKSANDLEARNINGSYSTEMS,VOL.23,NO.4,APRIL2012
(a)
Fig. 5. Feature vectors for the data (open dots), true pattern (solid curve),
and full manifold (surface) from the toy example, in the feature space. All
feature vectors mustresideexactly inthefullmanifold.
III. GEOMETRICINTERPRETATIONOF
KPCA DENOISING PERFORMANCE
In thissection,we providea geometricinterpretationof the
inconsistent denoisingperformanceseen in Fig. 2 (i.e., points
suchasx aredenoisedlesseffectivelythanpointssuchasx )
A B
and relate it to the orthogonal projection procedure common
to all of the KPCA methods. We argue that the orthogonal Fig.6. (a)Methodof[3]inthefeaturespace.Depictioninthefeaturespace
projection has undesirable properties and that modifying this oftheprincipalprojection Pφ˜(x)andpreimageprojectionφ˜(xˆ)forpointsxA
andxB,togetherwiththefullmanifoldandprincipalsubspace.(b)Expanded
procedure (as described in the next section) can substantially callout for xA.(c)Expandedcallout for xB.
improve denoising performance.
To illustrate, consider again the toy example introduced in
Section I, for which the data are distributed along a quadratic The reason behind the poor denoising of x and the good
A
curve in R2 with Gaussian noise added. In order to better denoisingof x becomesmoreapparentwhenwevisualizethe
B
visualizethefeaturespace(tofacilitatetheillustration),instead full manifold and the projection and preimage approximation
of using a feature map induced by a polynomial kernel, we steps in the featurespace. Recall thatthe fullmanifoldφ˜(Rn)
use the simple feature map φ(x) = (x ,x ,x2)(cid:2) which is is the image of Rn under the feature map, which is an n-
1 2 1
sufficiently rich to capture the quadratic pattern. Fig. 3(a) dimensional manifold in M-dimensional space. Fig. 5 plots
shows the noisy data and the true quadratic pattern in the the full manifold for the toy example. Notice that the feature
inputspace,andFig.3(b)showsthesamedataandtruepattern vectorforanydatapointmustlieinthisfullmanifold,regard-
mapped to the 3-D (centered) feature space. The data points less of the level of noise. More generally, the feature vector
are thefeaturevectors,andthe solidcurveisthe imageof the for any pointin Rn must lie on the full manifold,and the full
true pattern under the feature map. Ideally, we would like the manifold is the set of all points in the feature space that have
denoised observation to lie close to the solid curve. exact preimages. Consequently, after finding the projection
Consider the orthogonal projection Pφ˜(x) of a feature Pφ˜(x)ontotheprincipalsubspace,wecanviewφ˜(xˆ)fromthe
vector onto the principal subspace. For the toy example with preimageapproximationxˆ =argmin ||φ˜(z)−Pφ˜(x)||2 of[3]
z
L = 2, the projections of the feature vectors are plotted in as an orthogonal projection of Pφ˜(x) onto the full manifold
Fig.4(b)andtheirapproximatepreimagesusingthemethodof in the feature space. To distinguish between these two feature
[3]areplottedinFig.4(a).AswithlinearPCA,wewouldhope spaceprojections,wereferto Pφ˜(x)astheprincipalprojection
that the projections in the feature space are close to the true and φ˜(xˆ) as the preimage projection.
pattern.However,fromFig.4(b),theprojectionsofpointsnear Fig. 6 shows the principal projections and preimage pro-
the middleofthe curveare muchfartherfromthetrue pattern jections for pointsx and x , together with the full manifold
A B
than are the projectionsof pointsnear both endsof the curve. andtheprincipalsubspace.Thereasonwhyx isdenoisedless
A
Thisisconsistentwiththevariabledenoisingperformancethat effectively than x has much more to do with the principal
B
we see in Fig. 4(a) and discussed in the context of Fig. 2. projection than the preimage projection. Because of the local
IMetal.: TANGENTHYPERPLANEKERNELPCAFORDENOISING 649
geometryaroundφ˜(x ) of the full manifold and the principal Then we take the principal projection to be
A
subspace [see Fig. 6(b)], the principal projection Pφ˜(x ) for
A (cid:3)L
thepointxA fallsquitefarawayfromthefullmanifold.Thatis P φ˜(x)= Pφ˜(x)+ b∗v (2)
to say, it is quite far from having a preimage. Hence, in order T k k
k=1
to get back to the full manifold when finding the preimage
projection, we must choose a point φ˜(xˆA) that is quite far where a∗ = [a1∗, a2∗, …, an∗](cid:2) and b∗ = [b1∗, b2∗, …, b∗L](cid:2)
from the principalsubspace. In turn, this means that φ˜(xˆ ) is denote the coefficient vectors that minimize the distance. The
A
quitefar fromtheimage ofthe truepattern,whichmeansthat subscript T on the projection operator indicates it is for the
the denoised point xˆ is quite far from the true pattern (i.e., tangent hyperplane method.
A
poor denoising). ThesolutionisobtainedbydefiningU=[u1,u2,…,un]and
The situation for the point xB is very different. Because V = [v1, v2, …, vL], the sets of basis vectors for the tangent
of the local geometry around φ˜(x ) [see Fig. 6(c)], the hyperplane and for the principal subspace, respectively. The
B
principal projection Pφ˜(x ) remains relatively close to the optimal coefficient vector minimizes
B
fculollsemtaonitfhoeldp,risnocitphaalt sthuebsppraeciem.aTgheeprreosjueclttiiosnthφ˜a(txˆxˆB) riesmcaloinses (cid:10)(cid:10)(cid:10)(cid:11)φ˜(x)+Ua(cid:12)−(cid:11)Pφ˜(x)+Vb(cid:12)(cid:10)(cid:10)(cid:10)2
B
to the true pattern (i.e., good denoising). (cid:10)(cid:10)(cid:11) (cid:12) (cid:6) (cid:7)(cid:13) (cid:14)(cid:10)(cid:10)2
= (cid:10)(cid:10) φ˜(x)− Pφ˜(x) + U, −V a (cid:10)(cid:10)
b
(cid:10)(cid:11) (cid:12) (cid:10)
IV. TANGENTHYPERPLANE METHODFOR = (cid:10)(cid:10) φ˜(x)− Pφ˜(x) +Ac(cid:10)(cid:10)2
KPCA DENOISING
with respect to c = [a(cid:2), b(cid:2)](cid:2), and A = [U, −V], the solution
The arguments in the previous section indicate that the to which is the standard linear least squares solution.
denoisingperformanceisstronglyinfluencedbytheprojection Because the solution c∗ may not be unique (if the columns
step. The geometric interpretation in Fig. 6 motivates our of A are linearlydependent)ormay be poorlyconditioned(if
proposed modification to the projection step, which can be the columns of A are nearly linearly dependent), we add a
incorporated into the preimage approximation criteria for any regularization penalty and choose c∗ so that
of the KPCA denoising algorithms discussed previously. For
eachpointxtobedenoised,insteadofprojectingorthogonally c∗ =argmin||(φ˜(x)− Pφ˜(x))+Ac||2+δ||c||2 (3)
c
onto the principal subspace, we incorporate information on
the local geometry of the full manifold in the neighborhood for some small user-selected δ as in [14]. In Section V, we
of φ˜(x), so that the principal projection remains closer to the compare the results with different levels of δ and discuss
full manifold. This is accomplished as follows. its appropriate choice. The regularization penalty has other
First note that any po(cid:2)int in the principal subspace can be intuitive appeal. It pena(cid:2)lizes solutions for which the tangent
φrb˜eL(px]r(cid:2))es∈aesnRttehLde. anFs-odrPimaφ˜ne(yxn)sxi+o,nwaelkLh=dye1pfibenkrevpkltahnfeoert(asinnogmtehenetbMh=y-pd[eibmr1pe,lanbns2ie,o…nfoar,l hathnyedpr/eoarrwpltaahyneefprproiomnicniPtpφφa˜˜l((xsxu))+b.sHpeancnjce=e1p,aoinijnutthjPeisφe˜xf(utxrr)et+hmee(cid:2)rcakLaw=sae1ybofkfrvoδkm→isφ˜f∞(uxr)-,
feature space) that is tangent to the full manifold at φ˜(x). our principal projection coincides with the standard principal
More precisely, th(cid:2)e tangent hyperplane for φ˜(x) is the set of projection Pφ˜(x). Moreover, the regularization penalty keeps
all points φ˜(x)+ nj=1ajuj for a = [a1, a2, …, an](cid:2) ∈Rn, the solution well defined in the case where A is singular. For
where example,inthespecialcaseoflinearPCA[whichcorresponds
uj = ∂∂φ˜x(xj) =(cid:8) ∂φ∂˜1x(jx) ∂φ∂˜2x(jx) ··· ∂φ˜∂Mx(jx)(cid:9)(cid:2). tRpornitn,hcesiopfaeltahtsauutrbetshpmeaacpdeiφsatn(axdn)cte=heixst]a,znetghreoenttbaehntygwpeeenertpnhlayanpneey.rpIplnaocnilneutdisionnaglltthhoeef
regularization penalty ensures that the solution is unique and
In the following, we refer to the quantity that will replace coincides with the standard linear PCA solution in this case.
Pφ˜(x) in the projection step as the principal projection, It is straightforward to show that the solution to (3) is
althoughtechnicallyitisnotaprojection.Rather,theprincipal
projection in our algorithm is defined as the point on the c∗ =(A(cid:2)A+δI)−1A(cid:2)(Pφ˜(x)−φ˜(x)) (4)
principal subspace that is closest to the tangent hyperplane
for φ˜(x). This is equivalent to finding the minimum distance from which we obtain our principal projection as defined in
(2). One can then use the preimage approximation criterion
between the two subspaces (the tangent hyperplane and the
from any of the existing KPCA denoising methods to find
principalsubspace).Asin[14],whichcomputestheminimum the approximate preimage of P φ˜(x), which serves as the
distance between two hyperplanes in a different context, the T
denoised point xˆ.
solution is obtained by finding the coefficientvectorsa and b
Itisimportanttonotethatthecomputationalkerneltrickcan
to minimize the distance between
be used in all steps of the algorithm. Hence, no calculations
(cid:3)L (cid:3)n needtobeperformedinthe M-dimensionalfeaturespace.The
Pφ˜(x)+ b v and φ˜(x)+ a u . penalized least squares solution for the principal projection
k k j j
k=1 j=1 PTφ˜(x) requires only inner products in the feature space,
650 IEEETRANSACTIONSONNEURALNETWORKSANDLEARNINGSYSTEMS,VOL.23,NO.4,APRIL2012
and6),x isdenoisedsubstantiallybetter,and x isdenoised
A B
slightly better. In fact, x is now denoised as well as x .
A B
The reason for the dramatic improvement for the denoising
of points such as x is illustrated geometrically in Fig. 7(a),
A
which also shows the tangent hyperplanes for the two points
φ˜(x )andφ˜(x ).Theapproachforcestheprincipalprojection
A B
to be close to the tangent hyperplane [compare P φ˜(x ) in
T A
Fig. 7(b)with Pφ˜(x ) in Fig. 6(b)],whichalso keepsitclose
A
to the full manifold. If the principal projection is close to the
full manifold, then it is close to having an exact preimage,
so that its preimage projection φ˜(xˆ) is more likely to remain
closeto theprincipalsubspace.Theresultis thatthe denoised
point is more likely to be close to the true pattern.
V. RESULTSAND DISCUSSION
In this section, we apply the proposed tangent hyperplane
method (using the preimage approximation criterion of [3])
to three different datasets: the toy data that we have used
to illustrate the approach and two different image datasets
in which each multivariate observation represents a different
image. For the toy dataset, we use a feature map consisting
of all polynomial terms of degree ≤ 2, which corresponds
to the kernel K(x, y) = (1 + x(cid:2)y)2. For the image datasets,
we used the Gaussian RBF kernel K(x,y) = exp(–||x −
y||2/ρ) with ρ = nc, where c equals 0.5 times the average
of component variances, which we chose because it is more
appropriatethan the polynomialkernelson both methods.We
compare the tangent hyperplane approach with that of [3].
As we have discussed, the tangent hyperplane approach is
Fig. 7. Proposed method in the feature space. Geometric depiction of the a modification of the principal projection step that can be
tangent hyperplane method for denoising points xA and xB. Forcing the used with any of the preimage approximation criteria. For
principalprojection PTφ˜(x)tobeclosetothetangenthyperplanealsokeeps thepreimageapproximation,weusedthequasi-NewtonBFGS
ithtactlobseecatuosethtehefupllremimaangifeolpdroajnedcticolnosseφ˜t(oxˆAh)avainndgφ˜an(xˆeBx)acatrepsreoimclaogsee.toNothtiecier algorithmtominimize||φ˜(z)−Pφ˜(x)||2or||φ˜(z)−PTφ˜(x)||2.
respective principal projections PTφ˜(xˆA) and PTφ˜(xˆB), it is difficult to
visuallydistinguish theminthefigure.(b),(c)Expandedcallouts ofthetwo
A. Results for the Toy Example
regions in(a).
Thedataforthetoyexampleweregeneratedbystandardiz-
and the Appendix shows how all these can be calculated via ingasimulatedtrainingsetofsize N =100bivariateobserva-
the kernel trick. Because the principal projection is a linear tionsoftherandomvectorx=[x ,x ](cid:2)=[t,t2](cid:2)+[w ,w ](cid:2),
1 2 1 2
combinationoftheeigenvectors,thekerneltrickcanbeapplied where t follows a uniform distribution over the interval [−1,
as usual in any of the preimage approximation approaches. 1],andw andw areindependentGaussianrandomvariables
1 2
Regarding the computational complexity, because the tan- with zero mean and standard deviation 0.1. We can view the
genthyperplanemethodsimplymodifiestheprojectionstepof quadratic curve {[t, t2](cid:2): t ∈ [−1, 1]} as the true pattern, and
existingKPCAmethods,itscomputationalexpenseisthesame each observation x as a random point along the true pattern
asthatoftheexistingmethodsbutwiththeadditionalexpense plusGaussiannoise.Fig.8(a)showsthedenoisingresultusing
of the modified projection step. The most computationally themethodof[3],whichdiffersslightlyfromFig.4(a)because
expensivepartof the modified projectionstep is the inversion ofslightlydifferentchoiceoffeaturemap:φ(x)=(x ,x ,x2)(cid:2)
1 2 1
of(n+L)×(n+L) matrix A(cid:2)A + δI. Thisis generallymuch in F√ig. 4(√a) versus the√full quadratic feature map φ(x) =
less expensive than the remainder of the KPCA algorithm (1, 2x , 2x ,x2,x2, 2x x )(cid:2)inFig.8(a).Fig.8(b)shows
1 2 1 2 1 2
(especially if one uses a matrix inversion algorithm that takes the denoising result using the proposed tangent hyperplane
advantage of the orthogonality of the columns of A that are method.Forbothmethods, L =4 principaleigenvectorswere
associatedwith V).Hence,themodifiedprojectionstep ofthe retained. L =4 providedthe best denoising for both methods
tangenthyperplanemethoddoesnotsubstantiallyincrease the (results for other L are omitted for brevity).
computational expense of the entire KPCA algorithm. ComparingFig.8(a)and(b),itisclearthatforthisexample
In Fig. 7, we illustrate the tangent hyperplane method for the modified principal projection of the tangent hyperplane
the same toy example with 3-D feature space. The preimage method provides much better denoising than the conventional
approximation is still that of [3]. Relative to the denoising orthogonal principal projection. We also conducted simula-
results using the conventional orthogonal projection (Figs. 4 tionsforamodifiedtoyexamplewithhigherdimensionaldata.
IMetal.: TANGENTHYPERPLANEKERNELPCAFORDENOISING 651
(a)
(a) (b) 2 Method of [3]
δ =10−4 Our Proposed Method
0 δ =1
1.5 0
Distance 1
MS
R
0.5
00 50 100 150 200 250 300
L
(b)
3.5 Method of [3]
3 Our Proposed Method
2.5
F(ai)g.th8e. mFeothrotdheofto[y3]exaanmdp(leb,)athceomtapnagreinsotnhyopfetrhpeladneenomiseitnhgodr,eswuhltischushinags Distance 2
substantially better performance. MS 1.5
R
1
δ =10−4 δ =1
(a) 0.5 0 0
0
0 50 100 150 200 250 300
L
Fig. 10. RMS distance for the method of [3] (dashed) and the tangent
hyperplane method (solid), as a function of L for (a) the images with low
noiseand(b)theimageswithhighnoiselevels.Forthelownoise,themethod
of [3] shows much better performance than the tangent hyperplane method.
For high noise, the tangent hyperplane method with δ0 = 1 has the best
optimal performance for this example, and the tangent hyperplane method
(b) withbothsettings ismuchmorerobusttomisspecification of L.
Gaussian noise at two different noise levels. We standardized
the data before applying the denoising methods. The training
data consisted of N = 500 noisy images. Examples of the
original clean images and the noisy images for the two
different noise levels are shown in Fig. 9.
(c) As the denoising performance depends on choice of L
(the number of retained principal components) and δ [the
regularization parameter in (3)], we investigated the effects
ofthese parameters.Fig. 10plots, as a functionof L, the root
mean square (RMS) distance
(cid:15) (cid:16)
1
(cid:3)N 2
N−1 ||xi,true−xˆi||2
i=1
Fig. 9. Examples of (a) clean images, (b) images with low noise, and (c)
images with high noise. The noisy ones are the clean ones with Gaussian between the denoised images (xˆ ) and the original clean
i
noiseaddedwithdifferent standarddeviations (b)0.04and(c)0.12. images (xi,true) over the set of N training images. After
generating the noisy images, the original clean images were
The true pattern was a 1-D quadraticcurve in a 2-D subspace
usedonlyforevaluationpurposesbutnotwithinthealgorithms
ofR20,rotatedsothatithadcomponentsalongalleachnatural
in any way.
basis vector. The denoising results for the higher dimensional
Fig. 10(a) shows the denoising results for the images with
examplewerequitesimilartothoseshowninFig.8,especially
low noise. The dotted horizontal line is the RMS distance
the relative performance of the two methods. For brevity, we
between the noisy images and the original clean images,
omit the results for this example and instead focus on the
representing the performance with no denoising. The dashed
higher dimensional example in the following section.
curve is the RMS distance for the method of [3], and the two
solid curves are the RMS distance for the proposed tangent
B. Results for the Image Examples hyperplane method with two different values for δ, which
One of the two image datasets was obtained from the we take to be δ times the maximum eigenvalue of A(cid:2)A for
0
website http://cs.nyu.edu/∼roweis/data.html.Each image is of δ =1 and 10−4. Fig. 10(b) is analogous to Fig. 10(a) but
0
thesameperson’sfacebutwithdifferentfacialexpressionsand for the images with high noise. Notice that the method of [3]
angles.Theelementsofeachoriginalmultivariateobservation corresponds to the tangent hyperplane method with δ =∞.
0
x arethegrayscalelevels(overtherange[0,1])forn = 28× A number of notable characteristics can be observed from
20 = 560 pixels. The data to be denoised were contaminated Fig. 10. First, when the noise level is low, the method of [3]
versions of the images for which we added independent performs better, in terms of RMS distance, than the tangent
652 IEEETRANSACTIONSONNEURALNETWORKSANDLEARNINGSYSTEMS,VOL.23,NO.4,APRIL2012
(a) (a)
(b) (b)
(c) (c)
Fig. 11. Sample ofthedenoised image result forthelow noisecase using Fig.12. Sampleofthedenoisedimageresultforthehighnoisecaseusing
(a) method of [3] with L = 200, and (b) and (c) our method with L = (a) method of [3] with L = 30, and (b) and (c) our method with L = 75,
300,δ0=1and L =150,δ0=10−4,respectively. AlthoughRMSdistance δ0=1andL=100,δ0=10−4,respectively.Ourmethodunderbothsettings
corresponding to (c) is high, the denoised images are quite clean and retain showsbetterdenoisingperformanceoverall,whichisconsistentwithwhatwe
keycharacteristics/patterns intheoriginal data. observefromFig.10(b).
hyperplanemethodwithfiniteδ0.However,theRMSdistance for optimal choices of L per each method for the same 16
measure does not necessarily capture all relevant aspects of samples shown in Fig. 11. However, quite significantly, our
performance. For example, consider Fig. 11, which visually method is much more robust to the choice of L. For any
compares the denoised images for the two approaches for a L between roughly 50 and 200 for δ = 10−4, our method
0
sample of 16 images. For each method, L was chosen to be performs nearly as well as for its optimal L. In contrast, the
optimal,i.e., L =200wasusedforthemethodof[3],L =300 “sweet spot” for the method of [3] is much smaller, with
for the tangenthyperplanemethod with δ0 = 1, and L =150 performancerapidlydegradingfor L outsidetherange15–75,
for the tangenthyperplanemethod with δ0 = 10−4. Although roughly.This robustness to choice of L is important, because
thetangenthyperplanemethodresultsinworseRMSdistance, in practice it is difficult to determine the optimal L. When
it clearly provides more denoising and produces images that usingthe polynomialkernel, the same trendwas observedbut
appear cleaner than the images produced by [3]. The upside the overall performances of all methods were not as good.
of this is that the facial expressions are generally cleaner in Notice that plots such as Fig. 10 cannot be used in practice,
Fig. 11(b) and (c), than in Fig. 11(a). But the downside is because the plotted RMS distances were calculated using the
that, in addition to removing more of the noise, the tangent true noiseless images, which are not available in practice.
hyperplane method has removed some of the real patterns The second image dataset, which are also of faces but of
in the images. For example, although some of the facial a different nature than in the previous example, was obtained
expressionpatternsareretainedinFig.11(b)and(c),othersare from the website http://isomap.stanford.edu. As described in
notretained,andthepatternresultingfromthesubjectlooking [15], the dataset consists of a sequence of 4096-dimensional
to his right has been “denoised” out of the data [as apparent vectors, representing the grayscale levels (over the range [0,
from the second image in the top rows of Fig. 11(a)–(c)]. 1])forn =64×64=4096pixels.Eachimageisananimated
This is perhaps because the subject’s rotation to the far right facewithdifferentposesandlightingdirections.Thedatatobe
is something of an outlier expression. denoised were again contaminated versions of the images for
The results are quite different when the noise level is which we added independent Gaussian noise at two different
higher. Fig. 10(b) shows that, in this case, the tangent hyper- noise levels. As the data dimension (4096) is quite high, for
plane method with δ0 = 1 has the best optimal performance dimensionality reduction purpose we preprocessed the data
according the RMS distance measure, although the optimal using PCA before applyingthe denoising methods.The train-
performances of the other two methods are not far behind. ingdataconsistedof N = 698noisyimages.Examplesofthe
ThiscanalsobeseeninFig.12,whichagaincomparesresults originalcleanimagesandthenoisyonesareshowninFig.13.
IMetal.: TANGENTHYPERPLANEKERNELPCAFORDENOISING 653
(a) (a)
(b) (b)
(c) (c)
Fig. 13. Examples of (a) clean, (b) and (c) noisy image data. The noisy Fig.15. Sample ofthe denoised imageresult forthe low noisecase using
ones are the clean ones with Gaussian noise added with standard deviations (a)methodof[3]with L =250,and(b)and(c)ourmethodwith L =300,
(b)0.04and(c)0.12. δ0 =1and L =75,δ0 =10−4,respectively.
(a) (a)
12
10 δ 0 =10−4 δ 0 =1
Distance 68
MS
R 4
(b)
2 Method of [3]
Our Proposed Method
00 50 100 150 200 250 300
L
(b)
12
10 δ =1
0
S Distance 68 δ =10−4 (c)
RM 4 0
2 Method of [3]
Our Proposed Method
0
0 50 100 150 200 250 300
L
Fig. 14. RMS distance for the method of [3] (dashed) and the tangent
hyperplane method (solid), as a function of L for both (a) low noise and Fig.16. Sampleofthedenoisedimageresultforthehighnoisecaseusing
(b)highnoisecases.Asintheprevious example,themethodof[3]isbetter (a)method of[3] with L = 75, and (b) and(c) our method with L = 100,
inthelow noisecasewhile ourmethodbetter inhighnoisecase androbust δ0 =1and L =200,δ0 =10−4,respectively.
tothechoice of L.
somewhere between the low and high noise cases and the
As in the previousexample,we applied boththe methodof δ = 1 and δ =10−4 cases shown previously,we omit them
[3] and the tangent hyperplanemethod with δ = 1 and δ = 0 0
0 0 for brevity.
10−4. Figs. 14–16,which are analogousto Figs. 10–12,show
the results. The key observationsare quite consistentwith the
C. Choice of the Parameter δ and Kernels
previous example. In terms of RMS distance, our method is
outperformedbythemethodof[3]whenthenoiselevelislow. From the previous experimental results, it can be observed
The tangent hyperplane method appears to result in greater that δ substantially influences the denoising results. Although
“denoising,”butoflegitimatepatternsaswellasnoise.Athigh we envision no clear strategy for optimally choosing δ, we
noise level, on the other hand, the tangenthyperplanemethod do recommend that it be chosen to at least accomplish its
outperforms the method of [3] in terms of RMS distances as intendedpurpose:namely,to preventerratic results caused by
well as robustness to L. poorconditioningof A whenits columnsare(nearly)linearly
We also performedsimulationsformediumnoisecasesand dependent. As a rule of thumb, we suggest choosing δ =
for δ = 10−2. Because the results were predictable, falling 10−4 times the maximum eigenvalue of A(cid:2)A. One could also
0