# L1-norm Kernel PCA

We present the first model and algorithm for L1-norm kernel PCA. While L2-norm kernel PCA has been widely studied, there has been no work on L1-norm kernel PCA. For this non-convex and non-smooth problem, we offer geometric understandings through reformulations and present an efficient algorithm where the kernel trick is applicable. To attest the efficiency of the algorithm, we provide a convergence analysis including linear rate of convergence. Moreover, we prove that the output of our algorithm is a local optimal solution to the L1-norm kernel PCA problem. We also numerically show its robustness when extracting principal components in the presence of influential outliers, as well as its runtime comparability to L2-norm kernel PCA. Lastly, we introduce its application to outlier detection and show that the L1-norm kernel PCA based model outperforms especially for high dimensional data.

## Authors

• 2 publications
• 49 publications
• ### 2DR1-PCA and 2DL1-PCA: two variant 2DPCA algorithms based on none L2 norm

In this paper, two novel methods: 2DR1-PCA and 2DL1-PCA are proposed for...
12/23/2019 ∙ by Xing Liu, et al. ∙ 0

• ### Consistent Discretization and Minimization of the L1 Norm on Manifolds

The L1 norm has been tremendously popular in signal and image processing...
09/18/2016 ∙ by Alex Bronstein, et al. ∙ 0

• ### L1-norm Error Function Robustness and Outlier Regularization

In many real-world applications, data come with corruptions, large error...
05/28/2017 ∙ by Chris Ding, et al. ∙ 0

• ### Robust Linear Discriminant Analysis Using Ratio Minimization of L1,2-Norms

As one of the most popular linear subspace learning methods, the Linear ...
06/29/2019 ∙ by Feiping Nie, et al. ∙ 0

• ### Outlier Regularization for Vector Data and L21 Norm Robustness

In many real-world applications, data usually contain outliers. One popu...
06/20/2017 ∙ by Bo Jiang, et al. ∙ 0

• ### The Exact Solution to Rank-1 L1-norm TUCKER2 Decomposition

We study rank-1 L1-norm-based TUCKER2 (L1-TUCKER2) decomposition of 3-wa...
10/31/2017 ∙ by Panos P. Markopoulos, et al. ∙ 0

• ### A Multiphase Image Segmentation Based on Fuzzy Membership Functions and L1-norm Fidelity

In this paper, we propose a variational multiphase image segmentation mo...
04/09/2015 ∙ by Fang Li, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Principal Component Analysis (PCA) is one of the most popular dimensionality reduction techniques [1]. Given a large set of possibly correlated features, it attempts to find a small set of features (principal components

) that retain as much information as possible. To generate such new dimensions, it linearly transforms original features by multiplying

in a way that newly generated features are orthogonal and have the largest variances.

In traditional PCA, variances are measured using the -norm. This has a nice property that although the problem itself is non-convex, the optimal solution can be easily found through matrix factorization. With this property, together with its easy interpretability, PCA has been extensively used in a variety of applications. However, despite of its success, it still has some limitations. First, since it generates new dimensions through a linear combination of features, it is not able to capture non-linear relationships between features. Second, as it uses the -norm for measuring variance, its solutions tend to be substantially affected by influential outliers. To overcome these limitations, the following two approaches have been proposed.

Kernel PCA The idea of kernel PCA is to map original features into a high-dimensional feature space, and perform PCA in that high-dimensional feature space [2]. With non-linear mappings, we can capture non-linear relationships among features, and this computation can be done efficiently using the kernel trick. With the kernel trick, computations of principal components can be done without an explicitly mapping.

L1-norm PCA To alleviate the effects of influential observations, L1-norm PCA uses the L1-norm instead of the L2-norm to measure variances. The L1-norm is more advantageous than the L2-norm when there are outliers having large feature values since it is less influenced by them. By utilizing this property, more robust results can be obtained through the -norm based formulation in the presence of influential outliers.

In this paper, we combine these two approaches for the variance maximization version of -norm PCA (which is not the same as minimizing the reconstruction error with respect to the -norm). In other words, we tackle a kernel version of -norm PCA. Unlike -norm kernel PCA, the kernel version of -norm PCA is a hard problem in that it is not only non-convex but also non-smooth. However, through a reformulation, we make it a geometrically interpretable problem where the goal is to minimize the -norm of a vector subject to a linear constraint involving the -norm terms. For this reformulated problem, we present a ”fixed point” type algorithm that iteratively computes a -1,1 weight for each data point based on the kernel matrix and previous weights. We show that the kernel trick is applicable to this algorithm. Moreover, we prove the efficiency of our algorithm through a convergence analysis. We show that our algorithm converges in a finite number of steps and the objective values decrease with a linear rate. Lastly, we computationally investigate the robustness of our algorithm and illustrate its use for outlier detection.

Our work has the following contributions.

1. We are the first to present a model and an algorithm for -norm kernel PCA. While -norm kernel PCA has been widely used, a kernel version of -norm PCA has never been studied before. In this work, we show that the kernel trick which made -norm kernel PCA successful is also applicable for -norm kernel PCA.

2. We provide a rate of convergence analysis for our -norm kernel PCA algorithm. Although many algorithms have been proposed for -norm PCA, none of them provided a rate of convergence analysis. The work shows that our algorithm achieves a linear rate of convergence by exploiting the structure of the problem.

3. We introduce a methodology based on -norm kernel PCA for outlier detection.

In what follows, we always refer to the variance maximization version of

-norm PCA and we assume that every variable in the input data is standardized with a mean of 0 and standard deviation of 1. This paper is organized as follows. Section

2 reviews recent works on -norm PCA, and points out how our work is distinguishable from them. Section 3 covers various formulations of the kernel version of -norm PCA. Through reformulations, we offer an understanding of the problem and based on these understandings we present our algorithm in Section 4. Section 5 gives a convergence analysis for our algorithm. Lastly, we show the robustness of our algorithm and its application to outlier detection in Section 6.

## 2 Related Work

-norm PCA has been shown to be NP-hard in [3] and [4]. Nevertheless, an algorithm finding a global optimal solution is proposed in [3]. Utilizing the auxiliary-unit-vector technique [5], it computes a global optimal solution with complexity where is the number of observations, is the rank of the data matrix, and is the desired number of principal components. Assuming and are fixed, the runtime of this algorithm is polynomial in . However, if are large its computation time can be prohibitive. Rather than finding a global optimal solution which is intractable for general problems, our work focuses on developing an efficient algorithm finding a local optimal solution.

Recognizing the hardness of -norm PCA, an approximation algorithm is presented in [4] based on the known Nesterov’s theorem [6]. In this work, -norm PCA is relaxed to a semi-definite programming (SDP) problem and alternatively the SDP relaxation is considered. After solving the relaxed problem, it generates a random vector and uses a randomized rounding to produce a feasible solution. This randomized algorithm is a

-approximate algorithm in expectation. To achieve this approximation guarantee with high probability, it performs multiple times of randomized roundings and takes the one having the best objective value. Instead of providing an approximation guarantee by solving a relaxed problem, our work directly considers the

-norm kernel PCA problem, and develops an efficient algorithm finding a local optimal solution.

Another approach using a known mathematical programming model is introduced in [7]

. Specifically, it proposes an iterative algorithm that solves a mixed integer programming (MIP) problem in each iteration. Given an orthonormal matrix of loading vectors, it perturbs the matrix slightly in a way that the resulting matrix yields the largest objective value. After perturbation, it uses singular value decomposition to recover orthogonality. The algorithm is completely different from the one proposed herein, the objective values of the iterates do not necessarily improve over iterations. As opposed to it, our work shows monotone convergence of the objective values as well as a linear rate of convergence to a local optimal solution.

On the other hand, a simple numerical algorithm finding a local optimal solution is proposed in [8], and its extended version that finds multiple loading vectors at the same time is presented in [9]. In the former work, the optimal solution is assumed to have a certain form, and parameters involved in that form are updated at each iteration improving the objective values. For a linear kernel, our algorithm has the same form as this algorithm. However, while the algorithm in [8] is derived without any justification, we provide an understanding behind the algorithm. Moreover, we prove a rate of convergence analysis, and introduce a kernel version of it.

## 3 Kernel-based L1-norm PCA Formulations

We consider -norm PCA in a high-dimensional feature space . Suppose we map data vectors , into a feature space by a possibly non-linear mapping . Assuming for every , the kernel version of -norm PCA can be formulated as follows.

 maximizex∈F f(x)=n∑i=1|Φ(ai)Tx| (1) subject to ∥x∥2=1

As shown in (1), we only consider extracting the first loading vector. This assumption is justifiable since subsequent loading vectors can be found by iteratively running the same algorithm. Specifically, each time a new loading vector x is obtained, we update the kernel matrix defined by by orthogonally projecting onto the space orthogonal to the most recently obtained loading vector, and then run the same algorithm on the updated kernel matrix .

The problem (1) has a convex non-smooth objective function to maximize. Moreover, the feasibility set is non-convex. To better understand this problem and derive an efficient algorithm, we reformulate (1) in the following way.

 minimizex∈F g(x)=∥x∥2 (2) subject to n∑i=1|Φ(ai)Tx|=1

To prove the equivalence of two formulations, we show that an optimal solution of one formulation can be derived from an optimal solution of the other formulation.

###### Proposition 1.

The following holds.

1. [nolistsep]

2. If is optimal to (1), then is an optimal solution to (2).

3. If is optimal to (2), then is an optimal solution to (1).

###### Proof.

It is easy to check that is feasible to (2). If is not optimal to (2), there exists z such that . From its feasibility, it is obvious that holds. Then, for , we have

 f(w)=n∑i=1|Φ(ai)Tw|=(1/∥z∥2)n∑i=1|Φ(ai)Tz|=1/∥z∥2.

In the same way, we have

 f(x∗)=1/∥^x∗∥2

due to . From ,

 f(x∗)=1/∥^x∗∥2

must hold, which contradicts the assumption that is an optimal solution of (1).

Again, it is easy to check is feasible to (1). If is not optimal to (1), then there exists w such that

 n∑i=1|Φ(ai)T^y∗|

Since , for , we have

 g(z)=(1/n∑i=1|Φ(ai)Tw|)∥w∥2=1/n∑i=1|Φ(ai)Tw|.

On the other hand,

follows from resulting in

 g(y∗)=(1/n∑i=1|Φ(ai)T^y∗|)>g(z)=1/n∑i=1|Φ(ai)Tw|.

This contradicts the assumption that is optimal to (2). ∎

To understand formulation (2), we first examine the constraint set . Geometrically, this constraint set is symmetric with respect to the origin and represents a boundary of a polytope . It is easy to check that is a polytope since it can be represented by a finite set of linear constraints as

 P={x|(n∑i=1Φ(ai)ci)Tx≤1 where ci∈{−1,1},i=1,…,n}.

Therefore, formulation (2) is to find the closest point to the origin from the boundary of the polytope . The following proposition shows an optimal solution must be perpendicular to one of the faces.

###### Proposition 2.

An optimal solution is perpendicular to the face which it lies on.

###### Proof.

Suppose that an optimal solution of (2) is . Assuming , consider the face

 F={x|(n∑i=1Φ(ai)c∗i)Tx=1}∩P.

If is not perpendicular to face , then

 w=n∑i=1Φ(ai)c∗i/∥n∑i=1Φ(ai)c∗i∥22

is the closest point to the origin from

 {x|(n∑i=1Φ(ai)c∗i)Tx% =1}

having . Now, let us define its scalar multiple

 z=(1/n∑i=1|Φ(ai)Tw|)% w.

By construction, z is a feasible solution to (2) and has the objective value of

 g(z)=(1/n∑i=1|Φ(ai)Tw)∥w∥2.

From

 n∑i=1|Φ(ai)T(n∑j=1Φ(aj)c∗j)|−∥n∑i=1Φ(ai)c∗i∥22 =n∑i=1|Φ(ai)T(n∑j=1Φ(aj)c∗j)|−n∑i=1Φ(ai)Tc∗i(n∑j=1Φ(aj)c∗j) =n∑i=1(|Φ(ai)T(n∑j=1Φ(aj)c∗j)|−Φ(ai)Tc∗i(n∑j=1Φ(aj)c∗j))≥0,

we have

 n∑i=1|Φ(ai)Tw|=n∑i=1|Φ(ai)T(n∑j=1Φ(aj)c∗j)|∥n∑i=1Φ(ai)c∗i∥22≥1.

As a result,

 g(z) =(1/n∑i=1|Φ(ai)Tw|)∥w∥2 ≤g(w)=∥w∥2

follows. This contradicts the assumption that is an optimal solution to (2). Therefore, the optimal solution must be perpendicular to face . ∎

From Proposition 2, we can easily derive the following Corollary 1.

###### Corollary 1.

An optimal solution of (2) must have the form of where and .

The form that is a scalar multiple of is assumed in [8] for the linear kernel case without any justification but by Corollary 1, it follows that it is essentially the right form for the optimal solution. Moreover, from

 ∥x∗∥2 =∥y∗∥2/n∑i=1|Φ(ai)T% y∗|=∥y∗∥2/n∑i=1c∗iΦ(ai)Ty∗ =∥y∗∥2/n∑i=1n∑j=1c∗ic∗jΦ(ai)TΦ(aj)=∥y∗∥2/∥y∗∥22 =1/∥n∑i=1Φ(ai)c∗i∥2 (3)

we can further show that the optimal solution of formulation (2) can be derived from the optimal solution of the following binary problem.

 maximizec∈{−1,1}n ∥n∑i=1Φ(ai)ci∥22. (4)
###### Proposition 3.

Let an optimal solution of binary formulation (4) be . Then, satisfies , for . Therefore, it follows that is an optimal solution for formulation (2).

###### Proof.

Since is an optimal solution of (4), flipping the sign of any must not improve the objective value, .

To deduce a contradiction, let us assume that there exists some such that for .

Then, for any , flipping the sign of gives

 ∥y∗−2Φ(aj)c∗j∥22 =∥y∥22−4yT(Φ(aj)c∗j)+4∥Φ(aj)∥22 =∥y∥22+4|yT(Φ(aj))|+4∥Φ(aj)∥22 >∥y∥22=∥n∑i=1Φ(ai)c∗i∥22

which contradicts the assumption that is an optimal solution to (4). Therefore, must satisfy

 c∗i=sgn(Φ(ai)Ty∗) for i=1,…,n.

Since maximizes ,

 x∗=(1/n∑i=1|Φ(ai)Ty∗|)y∗

is a minimizer for (2) due to Corollary 1 and (3). ∎

The following result has been shown in [3] for the linear kernel case but here we generalize it.

###### Corollary 2.

Formulation (2) is equivalent to formulation (4).

###### Proof.

Due to Corollary 1 and (3), formulation (2) comes down to

 maximizec∈{−1,1}n ∥n∑i=1Φ(ai)ci∥22 subject to y=n∑i=1Φ(ai)ci ci=sgn(Φ(ai)Ty),i=1,...,n.

Since an optimal solution to (4) satisfies the constraints by Proposition 3, the two formulations are essentially the same. ∎

Interestingly, the binary formulation (4) is NP-hard. Since expanding gives

 n∑i=1n∑j=1Φ(ai)TΦ(aj)+∑i

assuming be the weight of the complete graph having nodes , we can show the equivalence of the quadratic binary program (4) and the max-cut problem.

## 4 Algorithm

In this section, we develop an efficient algorithm that finds a local optimal solution to problem (2) based on the findings in Section 3. Before giving the details of the algorithm, we first provide an idea behind the algorithm.

The main idea of the algorithm is to move along the boundary of so that the -norm of the iterate decreases. Figure 1 graphically shows a step of Algorithm 1. Starting with iterate

, we first identify hyperplane

which current iterate lies on. After identifying the equation of , we find the closest point to the origin from , which we denote by . After that, we obtain by projecting to the constraint set by multiplying it by an appropriate scalar. We repeat this process until iterate converges.

Next, we develop an algorithm based on the above idea. Let . From Corollary 1, we know that the optimal solution has the form of

 x∗ =(1/n∑i=1|Φ(ai)Ty∗|)y∗=(1/n∑i=1ciΦ(ai)Ty∗)y∗ =(1/(c∗)TKc∗)n∑i=1Φ(ai)c∗i.

Utilizing the fact that the optimal solution is characterized by the sign vector , we characterize the initial iterate with the sign vector as

 x0=(1/(c0)TKc0)n∑i=1Φ(ai)c0i.

With , the equation of the hyperplane is represented by

 (yk)T(x−xk)=0.

The closest point to the origin among the points in the hyperplane has the form of . By plugging into , we have

 s=(yk)T(xk)(yk)T(yk), zk=(yk)T(xk)(yk)T(yk)yk.

We multiply by to make it feasible and thereby get

 xk+1=(1/n∑i=1|Φ(ai)Tzk|)zk. (5)

Utilizing

 (yk)T(xk)=n∑i=1Φ(ai)Txkcki=n∑i=1|Φ(ai)Txk|=1, (6)

we get the followings.

 yk=n∑i=1Φ(ai)cki, (7) zk=yk∥yk∥22, (8) xk+1=(1/n∑i=1|Φ(ai)Tyk|)yk. (9)

By plugging (7) into (9), can be represented as

 xk+1 =(1/n∑i=1|Φ(ai)Tyk|)yk =(1/(ck)TKck)n∑i=1Φ(ai)cki.

This implies that we only need to update at each iteration.
From

 ck+1i =sgn((Φ(ai))Txk+1)=sgn((Φ(ai))Tyk) =sgn(n∑j=1ckjKi,j),

we only require to compute

 ck+1i=sgn(n∑j=1ckjKi,j)

at each iteration.
From

 ∥xk+1−xk∥22=0⟺(ck−ck+1)TK(ck−ck+1)=0,

we get the following termination criteria:

 (ck−ck+1)TK(ck−ck+1)=0,

resulting in Algorithm 1.

After getting the final from Algorithm 1, we can compute principal scores without explicit mapping . For example, the principal component of observation can be computed by

 Φ(ai)Tx∗∥x∗∥2 =Φ(ai)Ty∗∥y∗∥2 =n∑j=1Φ(ai)TΦ(aj)c∗j√n∑i=1n∑j=1Φ(ai)TΦ(aj)c∗ic∗j =n∑j=1Ki,jc∗j√n∑i=1n∑j=1Ki,jc∗ic∗j.

We can also proceed to find more principal components without explicit mapping . As computing a loading vector and principal components only require the kernel matrix , we only need to update the kernel matrix each time a new loading vector is found. We can update the kernel matrix without explicit mapping by

 ˜Ki,j =(Φ(ai)−Φ(ai)Tx∗∥% x∗∥22x∗)T(Φ(aj)−Φ(aj)Tx∗∥x∗∥22x∗) =Φ(ai)TΦ(aj)−Φ(ai)Tx∗Φ(aj)Tx∗∥x∗∥22 =Ki,j−(n∑k=1Ki,kc∗k)(n∑k=1Kj,kc∗k)(n∑i=1n∑j=1Ki,jc∗ic∗j).

## 5 Convergence Analysis

In this section, we provide a convergence analysis of Algorithm 1. We first prove finite convergence, and then provide a rate of convergence analysis.

Before proving the convergence of the algorithm, we first show that the sequence is non-increasing.

We have .

###### Proof.

Inequality follows from

 ∥xk∥22−∥zk∥22 =∥xk∥22−1∥yk∥22 =∥xk∥22−((yk)T(xk))2∥yk∥22 =∥xk∥22∥yk∥22−((y% k)T(xk))2∥yk∥22≥0. (10)

Here, the second equality is from (6) and the last inequality holds by Cauchy-Schwarz inequality where the equality holds if is a scalar multiple of .
Next, we have

 n∑i=1|Φ(ai)Tzk| =n∑i=1|Φ(ai)T(yk)|(yk)T(yk) =1(yk)T(yk)n∑i=1|Φ(ai)Tyk| =n∑i=1|Φ(ai)T(n∑j=1Φ(aj)ckj)|n∑i=1n∑j=1Φ(ai)TΦ(aj)ckickj =n∑i=1|(Φ(ai)cki)T(n∑j=1Φ(aj)ckj)|n∑i=1n∑j=1Φ(ai)TΦ(aj)ckickj =n∑i=1|n∑j=1Φ(ai)TΦ(aj)ckickj|n∑i=1n∑j=1Φ(ai)TΦ(aj)ckickj≥1. (11)

Finally, follows from

 ∥xk+1∥22=(1/n∑i=1|Φ(ai)Tzk|)2∥zk∥22≤∥zk∥22.

###### Lemma 2.

If , we have , and .

###### Proof.

From Lemma 1, we have . Then, from (10), is a scalar multiple of . Assuming , follows from (6) resulting in

 xk=yk∥yk∥22.

We can show in the same manner. As a result, holds by (8). From , we have

 xk+1=(1/n∑i=1|Φ(ai)Tzk|)zk=(1/n∑i=1|Φ(ai)Txk|)x% k=xk.

###### Theorem 1.

The sequence converges in a finite number of steps.

###### Proof.

Suppose the sequence does not converge. As vector is solely determined by , the number of possible is finite. Therefore, if the sequence does not converge, then some vectors appear more than once in the sequence .

Without loss of generality, let . By Lemma 1, we have

 ∥xl+m∥2=∥xl∥2≥∥xl+1∥2≥...≥∥xl+m∥2

forcing us to have . Now, by Lemma 2, must hold, which contradicts the assumption that the sequence does not converge. In other words, the algorithm stops at iteration . Therefore, the sequence generated by Algorithm 1 converges in a finite number of steps. ∎

Next, we prove that the sequence generated by Algorithm 1 converges with a linear rate.

###### Theorem 2.

Let Algorithm 1 start from and terminate with at iteration . Then we have where for all .

###### Proof.

From (5), we have

 ∥xk∥2=(1/n∑i=1|Φ(ai)Tzk−1|)∥zk−1∥2.

Since by Lemma 1, we have

 ∥xk∥2≤(1/n∑i=1|Φ(ai)Tzk−1|)∥xk−1∥2. (12)

Now, we show

 ∥xk∥2−∥x∗∥2≤(1/