Sep 082014
 

拉普拉斯方法(Laplace’s Method)

拉普拉斯方法又称为拉普拉斯近似(Laplace Approximation)。它可以用来计算一元或多元积分[1]。

举例来说,假设 \(f(x)\) 是一维函数,我们要计算
\( \int_{-\infty}^{\infty} f(x) \mathrm{d} x \).

如果\(f(x)\)形式很复杂,我们往往找不到定积分的公式(Close form)。
如果用数值方法来计算积分,计算量又很大。
所以要想个办法,得到一个比较精确的积分结果。

拉普拉斯方法可以适用于这种情况。我们先用泰勒展开(Taylor Expansion):

$$ f(x) \approx f(x_0) + f'(x_0)(x-x_0) + \frac{1}{2}f”(x_0) (x-x_0)^2 $$
如果选取\(x_0\)使得\(f'(x_0) = 0\),则可以进一步简化为
$$ f(x) \approx f(x_0) + \frac{1}{2}f”(x_0) (x-x_0)^2 $$

再引入一个假设,\(\int_{-\infty}^{\infty} f(x) \mathrm{d} x \)
可以被写成 \(\int_{-\infty}^{\infty} e^{f(x)} \mathrm{d} x \) 的形式,那么就有:

$$ \begin{align}
\int_{-\infty}^{\infty} e^{f(x)} \mathrm{d} x
&\approx \int_{-\infty}^{\infty} e^{f(x_0) + \frac{1}{2}f”(x_0) (x-x_0)^2} \mathrm{d} x \\
&= e^{f(x_0)} \int_{-\infty}^{\infty} e^{- \frac{1}{2}|f”(x_0)| (x-x_0)^2} \mathrm{d} x \\
&= e^{f(x_0)} \sqrt{\frac{2\pi}{|f”(x_0)|}}
\end{align}
$$

注意上式中\(f”(x_0)\)要取绝对值。因为\(f'(x_0) = 0\),\(x_0\)是\(f(x)\)的极值点。
对于概率密度函数而言,一般它也是最大值点(mode), 因此\(f”(x_0) < 0[/latex]. 从几何上讲,拉普拉斯方法是要用一个[latex]e^{-x^2}[/latex]形式的函数近似[latex]f(x)[/latex]. 或者说,要用一个长的像正态函数的函数取近似,这样的函数服从[latex]N(x_0, 1/f''(x_0))[/latex]. 下面举个例子来计算 $$\int_{-\infty}^{\infty} x^2 e^{-x^2/2} \mathrm{d} x $$. 这个积分的结果是[latex]\sqrt{2\pi}[/latex],这可以这样计算出来。 假设[latex]Z[/latex]是一个标准正态分布,那么: $$\int_{-\infty}^{\infty} x^2 e^{-x^2/2} \mathrm{d} x = \sqrt{2\pi} E(Z^2) = \sqrt{2\pi} Var(Z) = \sqrt{2\pi} = 2.507$$. 用拉普拉斯方法: $$\int_{-\infty}^{\infty} x^2 e^{-x^2/2} \mathrm{d} x = \int_{-\infty}^{\infty} e^{-x^2/2 + 2 \log{x}} \mathrm{d} x $$. 对[latex]f(x)[/latex]求导: $$ \begin{align} f(x) &= -\frac{x^2}{2} + 2 \log{x} \\ f'(x) &= -x + \frac{2}{x} \\ f''(x) &= -1 - \frac{2}{x^2} \end{align} $$ 令[latex]f'(x) = 0[/latex],可以解出[latex]x_0 = \sqrt{2}[/latex]. 我们有[latex]f(x_0) = -1 + \log{2}[/latex] 和 [latex]f''(x_0) = -2[/latex]. 积分结果为: $$ e^{f(x_0)} \sqrt{\frac{2\pi}{|f''(x_0)|}} = e^{-1 + \log{2}} * \sqrt{\frac{2\pi}{|-2|}} = \frac{2}{e} \sqrt{\pi} = 1.304 $$ 比较2.507 和1.304,我们发现理论结果和拉普拉斯方法的结果几乎差了一倍。 为啥差别那么大呢? 看下面的图: [caption id="attachment_564" align="alignnone" width="300"]Demo Laplace Method Demo Laplace Method[/caption]

我们发现要求的积分值是黑色曲线下的面积。拉普拉斯方法想用灰色曲线([latex]e^{-\frac{x^2}{2}}\))通过拉伸得到红色曲线,然后用红色曲线下的面积来近似积分。但是由于黑色曲线有两个峰值,这个近似显然不算成功。
所以拉普拉斯方法是有局限性的:被积分的函数有一个峰值,并且和正态曲线长的像。这种情况下的近似才能比较精确。

对于概率密度函数,(无论一维还是多维),大部分都是一个峰值。或者因为中心极限定理,统计量均值的分布和正态函数很像,也是一个峰值。因此拉普拉斯方法的用处还是很多的。此外,拉普拉斯方法计算很快。比如在线性混合效果模型中,拉普拉斯方法是用的最广的方法,也可能是唯一能在实际中使用的方法[2-4]。

[1]维基百科 http://en.wikipedia.org/wiki/Laplace’s_method
[2]Approximate Inference in Generalized Linear Mixed Models. N. E. Breslow and D. G. Clayton. Journal of the American Statistical Association Vol. 88, No. 421 (Mar., 1993) , pp. 9-25
[3]Variance component testing in generalised linear models with random effects. XIHONG LIN. Biometrika (1997) 84 (2): 309-326.
[4]lme4 package Douglas Bates et al. http://cran.r-project.org/web/packages/lme4/lme4.pdf (see nAGQ parameter)

Jul 062014
 

卡方检验
Chi-square test

卡方检验的全称是Pearson’s Chi-square test,这是统计中最重要的检验。
从用途来讲, 卡方检验可以做两件事:
(1)适配度检定(Goodness of fit),检查样本是否符合某种随机分布;
(2)独立性检定(Independence test),检查多个变量之间是否独立。

这篇不讲怎么用卡方检验,而是讲讲卡方检验的公式,推导,以及和似然检验的联系,目的是温故知新。

从公式来讲卡方检验:

$$ \chi^2 = \sum_i \frac{(O_i – E_i)^2}{E_i} $$
这里\(O_i\)是观察到的次数,\(E_i\)是期望的次数,
统计量\(\chi^2\) 服从自由度\(df\)的卡方分布(\(\chi^2\) distribution)。

假设有一个有可能是多项分布的随机变量\(X \sim Multinomial(N, p_1, p_2, \ldots, p_k)\)。
观察到\(k\)类的次数分别是\(x_1, x_2, \ldots,x_k\).
为了统计检验\(X\)是否服从这个分布,列出零假设\(H_0\) 为\(X\)服从多项分布,对立假设为\(X\)不服从这个多项分布,可以得到:

$$ \chi^2 = \sum_i \frac{(O_i – E_i)^2}{E_i} = \sum_i \frac{(O_i – N \cdot p_i)^2}{N \cdot p_i}$$

注意这个每一项的分母是\( N \cdot p_i \), 但根据中心极限定理(Central Limit Theorem),
$$\frac{O_i – N \cdot p_i}{\sqrt{N \cdot p_i \cdot (1-p_i)}} \rightarrow Normal(0, 1)$$

同卡方检验的公式相比较,这里面有两个问题:
1. 为什么卡方检验不是\(N\)个\(\frac{O_i – N \cdot p_i}{\sqrt{N \cdot p_i \cdot (1-p_i)}}\)的平方和?
2. 为什么\( \chi^2 \)的自由度是\( n-1 \),而不是 \( n \)?

回答这两个问题的一个方法是推导一下 \( \chi^2 \) 的分布。
根据[1],推导的思路如下:

1. 令\( Z_i = \frac{O_i – N * p_i}{\sqrt{N*p_i}} \),因此\( \chi^2 = \sum_i Z_i^2 \).

2.
根据中心极限定理,
\( Z_i \rightarrow \sqrt{1-p_i} N(0, 1) = N(0, 1-p_i) \)
\( Cov(Z_i, Z_j) = \frac{Cov(O_i, O_j)}{n*\sqrt{p_i*p_j}} = -\sqrt{p_i * p_j} \)

3. 为了得到\( Z_1^2 + Z_2+ \ldots + Z_k^2 \) 的极限分布,
我们先假设\( \vec{g} = (g_1, g_2, \ldots, g_r) \) 是i.i.d.的正态分布,
再假设单位向量 \( \vec{p} = (\sqrt{p_1}, \sqrt{p_2}, \ldots, \sqrt{p_k}) \).
当\( \vec{g} \) 向量 减去 它在\( \vec{p} \)的投影时,
\( \vec{g^p} = \vec{g} – ( \vec{g} \cdot \vec{p} ) \cdot \vec{p} \) 的分布和 \( Z_1, Z_2, \ldots, Z_k \) 是 一样的。
为了证明这一点, 可以证明 \( E(Z_i^2) \rightarrow (g^p_i)^2 \) 和 \( E(Z_i Z_j) \rightarrow g^p_i \cdot g^p_j \) 。
有了\( \vec{g^p} \), 就可以把\( \chi^2 \) 表示为 \( \chi^2 = \sum_i (\vec{g^p}_i)^2 \) 即 \( \vec{g^p} \) 长度的平方。

4. 下面计算\( \vec{g^p} \) 的分布
\( \vec{g} \) 的每个分量是i.i.d的正态分布。
根据正态分布的性质,对\( \vec{g} \) 的坐标进行正交变换后,在新坐标下,每个坐标仍然是i.i.d的正态分布 (见[2])
因此我们构造出一个特殊的正交变换使得\( \vec{g} \) 变换之后有一个坐标系和\( \vec{p} \)同一个方向,
因此可以把\( \vec{g} \) 表示为\( (g’_1, g’_2, \ldots, g’_k) \),并且
\( \vec{g^p} = (g’_1, g’_2, \ldots, g’_{k-1}, 0) \)。
因此\( \chi^2 = \sum_{i=1}^{k-1} (g’_i)^2 = \chi^2_{k-1} \)。

以上四个步骤可以从理论上证明卡方检验。另外一种思路是把似然检验(G-test)和卡方检验联系起来:
$$ \begin{align} G
&= 2 \sum_i O_i \log(\frac{O_i}{E_i}) \\
&= 2 \sum_i O_i \log(1 + \frac{O_i – E_i}{E_i}) \\
&\sim 2 \sum_i O_i \frac{O_i – E_i}{E_i} \\
&= 2 \sum_i \frac{O_i^2 – E_i * O_i}{E_i}\\
&= 2 \sum_i \frac{O_i^2 – 2 * E_i * O_i + E_i^2}{E_i}\\
&= 2 \sum_i \frac{(O_i – E_i)^2} {E_i}
\end{align}
$$
上面的计算中用到了\( \log(1+x) \sim x \) 和 \( \sum_i O_i = \sum_i E_i \)。
可见两种检验基本上是等价的(上面推导来自[3],但更简洁)。

然而,似然检验和卡方检验相比,似然检验需要计算对数,在皮尔森的年代,这是很不方便的。
卡方检验则便于手算,这或许是得到广泛应用的原因之一吧。

[1]Pearson’s Theorem http://ocw.mit.edu/courses/mathematics/18-443-statistics-for-applications-fall-2003/lecture-notes/lec23.pdf
[2]Orthogonal Transformation of Standard Normal Sample http://ocw.mit.edu/courses/mathematics/18-443-statistics-for-applications-fall-2003/lecture-notes/lec15.pdf
[3]The Two-Way Likelihood Ratio (G) Test http://arxiv.org/pdf/1206.4881v2.pdf

Sep 262013
 

Ridge regression 相关的资料

以前学Ridge regression,大概记住了如何求回归系数 \(\beta^{ridge}\),但有些细节不了解。我找到下面两个slide感觉又学到了新东西。

 

第一个讲义 【Link or RidgeRegression – Ryan Tibshirani

http://www.stat.cmu.edu/~ryantibs/datamining/lectures/16-modr1.pdf

Slide 6:Shrink之后的回归系数看起来什么样?

Slide 7:为什么要center和normalize predictors?

Slide 13:当真正的回归系数很大的时候,Ridge regression还有好的表现么?或者说当true coefficient很大的时候,应该向什么方向调整\(\lambda\)可以减少prediction error?

 

第二个讲义 【Link or RidgeRegression-Rudyregularization

http://www-stat.stanford.edu/~owen/courses/305/Rudyregularization.pdf

Slide 17:Ridge regression的回归系数的bias具体是多少?

Slide 23-24:Ridge regression和Principal component analysis (PCA)的关系是什么?

Slide 37:选择Ridge regression的tuning parameter时一般用cross validation。如果用GCV(Generalized cross validation)的话,可以有close form.

May 092011
 

本文介绍如何使用R软件来分析一维随机变量。分析的内容包括如何查找一维数据的分布类型,如何估计分布参数以及如何用假设检验来测试一维数据的分布类型。
How to find, fit, test the distribution of univariate variable in R?
我们经常见到一维随机变量,比如线性模型的响应,我们通常需要检验它是否是正态分布来决定模型中直接用Y还是用log(Y),或者其他的transformation。
本文主要参考【1】,我会介绍一些基本的方法,但建议读者参考原文获得更多的信息。

1. 画密度图,CDF图

直方图:history(x)
密度图:plot(density(x))
CDF图:plot(ecdf(x))

检查是否是正态分布:

z= (x-mean(x))/sd(x)
qqnorm(z)
abline(0,1)

类似的可以检查其他分布(先构造一个理论分布,再qqnorm)

x.wei <- rweibull(200, shape=2.1, scale=1.1)
x.teo <- rweibull(200, shape=2.1, scale=1.0)
qqplot(x.teo, x.wei)
abline(0,1)

http://www.statsoft.com/textbook/distribution-fitting/

2. 利用矩估计猜测分布类型
主要是standardize之后计算一二三四阶矩(moment),然后对比下面网页列举的常见分布,猜出到底是哪一种分布:
NIST 1.3.5.11. Measures of Skewness and Kurtosis

3. 估计分布参数
当我们知道分布类型后,可以估计分布参数,常见的有矩估计和最大似然估计。
矩估计相对简单,可以用mean,var函数计算,但可能不具有无偏的性质。
最大似然估计有
1) mle() 在 stats4 包里
2) fitdistr() 在 MASS 包里
1)的方法显然更基本,但能适用于各种分布,2)的方法使用简单,对Gamma, Weibull, Normal等分布只需要一个命令,例如:

fitdistr(x.norm,"normal") ## fitting gaussian pdf parameters 
mean	sd
9.9355373 2.0101691 
(0.1421404) (0.1005085)

4. 检查分布是否合适?
在做Goodness of fit tests之前,可以先画出直方图和理论密度分布图。
之后,可以利用卡方检验来做Goodness of fit tests。具体来讲:
i) 对于Poisson, binomial, negative binomail, 我们可以使用vcd包中的goodfit函数。
ii) 对于一般的分布,可以把变量归类,然后利用卡方检验公示计算观察到变量数量和理论值之间的差异,然后计算pvalue
iii) 对于一般的分布,也可以使用Kolmogorov-Smirnov test来做统计检验

对第三种情况举例如下:

> x.wei <- rweibull(n=200, shape=2.1, scale = 1.1)
> ks.test(x.wei, "pweibull", shape=2, scale= 1)

	One-sample Kolmogorov-Smirnov test

data:  x.wei 
D = 0.1042, p-value = 0.02591
alternative hypothesis: two-sided 

特别的,我们需要检查数据是否是正态分布。
最常用的是Shapiro-Wilk test:shapiro.test()
此外,R里面有一个package nortest,提供了另外5种检查正态分布的函数:
i) Shapiro-Francia test: sf.test()
ii) Anderson-Darling test: ad.test()
iii) Cramer-Von Mises test: cvm.test()
iv) Lilliefors test: lillie.test() 适用于小样本,参数未知的正态分布
v) pearson.test: pearson.test()
这5种test各有细致的差异,使用的时候需自己区分。

参考文献:
【1】 http://cran.r-project.org/doc/contrib/Ricci-distributions-en.pdf

Feb 212011
 

We will briefly explain why we would be interested in implementing exact logistic regression, then provides C++ and R codes.

1. Why exact test?

Since we want to have a clear mind of how likely/unlikely the realization we observed. In the classic example of 2×2 table without covariates, especially the 2×2 table has very few (<5) occurrence, Fisher exact tests are often applied, and large sample theory cannot give an accurate estimation.

2. Why exact logistic regression?

Fisher’s exact test cannot applied to logistic model. For example, when we have covariates in the model, we want BOTH estimate the effect size and get its exact p-value. In this case, only exact logistic regression provides solution.The theoretical background is provided in reference [1].

My implementation:

Download: exactLogisticRegression.tar

1. I verified the results with SAS.
2. The speed is comparable to, or faster than SAS.

Cons:
1. Only 1 interested parameter conditioning on all other parameter is supported for now.
2. I have not implemented the confidence limit parts, as it’s a bit more tedious.

R binding : mypackage_1.0.tar

See mypackage/R/rcpp_hello_world.R, I wrote a R function to wrap the C++ function.

R binding is helped by RCpp package. It greatly reduced the workload of exchanging date (in the form of matrix, list, vector) between C++ and R. A quick tutorial can be found from RCpp homepage(http://dirk.eddelbuettel.com/code/rcpp.html). For experienced Rcpp user, the quick-ref documentation (http://dirk.eddelbuettel.com/code/rcpp/Rcpp-quickref.pdf) is helpful.

【1】 Exact Logistic Regression, Robert E. Derr, SAS Institute Inc., Cary, NC http://support.sas.com/rnd/app/da/new/daexactlogistic.html

【2】Rcpp: Rcpp: Seamless R and C++ Integration dirk.eddelbuettel.com/code/rcpp.html

Feb 082011
 

利用矩阵的分解和分析图是个很有意思的话题。当我们能用这个技术来改进PCA的时候,或者降维的时候,我们有可能相信有意思的结果会蹦出来。这里主要参考了文献【1】和Pluskid的blog【2】。其中【1】给出了推导过程:目标函数是二次型的矩阵,约束同样是二次型的;还有详细的Algorithm:里面最关键的一步是Generalized eigenvector problem(wiki有非常简短的介绍),理论上可以用Golub Matrix Computation Chapter 8 的方法(我没读,差不多忘了Numerial Method课的知识了),但我并没有使用。另外【2】里的文字流畅,言简意赅,是入门的好文章。

简单来讲,当利用k-neighbor 或\( \epsilon \)方法构造临街矩阵,利用Simpled minded(0 or 1)或者Heat Kernel来构建Weight矩阵后,我们的问题是求解:

$$ L f = \lambda D f \quad st. \quad f^T D f = 1 \quad \mathrm{and} \quad f^T D \mathbf{1} = 0 $$

我们已经知道D是一个对角阵,所以求解$$ D^{-1} L f = \lambda f $$即可。

考虑到约束条件,我们只需对每一个\(f[\latex]做如下变换:

$$ f^T D f = 1 \Rightarrow f = f / \sqrt{\sum_i f_i^2 d_i} $$

$$ f^T D \mathbf{1} = 0 \Rightarrow f = f – \frac{ \sum_i f_i d_i } { \sum_i d_i} $$

R code

LaplacianEigenmap.R Link

LaplacianEigenmapTest.R Link

结果:

从 [latex] 0 = \lambda_0 \le \lambda_1 \le \lambda_2 \ldots \lambda_n \) 中取最小的两个非零的特征根\( \lambda_1\, \lambda_2 \)对应的特征向量,仿照paper,得到结果如下:

注意,这里的图案和paper不符,但我的验算,检查是否是特征值、约束条件,表明我的计算过程应该是正确的。

另外Pluskid 的Blog上的图案中,Laplacian Eigenmap的结果是一个彩带, 我认为有可能是使用了\( \lambda_0, \lambda_1 \)对应的特征向量。

参考:

【1】Mikhail Belkin and Partha Niyogi, “Laplacian Eigenmaps for Dimensionality Reduction and Data Representation,” Neural Computation 15, no. 6 (February 6, 2011): 1373-1396.

【2】漫谈 Clustering (番外篇): Dimensionality Reduction http://blog.pluskid.org/?p=290

Feb 052011
 

主成份(Principal Component Analysis)分析是降维(Dimension Reduction)的重要手段。每一个主成分都是数据在某一个方向上的投影,在不同的方向上这些数据方差Variance的大小由其特征值(eigenvalue)决定。一般我们会选取最大的几个特征值所在的特征向量(eigenvector),这些方向上的信息丰富,一般认为包含了更多我们所感兴趣的信息。当然,这里面有较强的假设:(1)特征根的大小决定了我们感兴趣信息的多少。即小特征根往往代表了噪声,但实际上,向小一点的特征根方向投影也有可能包括我们感兴趣的数据; (2)特征向量的方向是互相正交(orthogonal)的,这种正交性使得PCA容易受到Outlier的影响,例如在【1】中提到的例子(3)难于解释结果。例如在建立线性回归模型(Linear Regression Model)分析因变量(response)和第一个主成份的关系时,我们得到的回归系数(Coefficiency)不是某一个自变量(covariate)的贡献,而是对所有自变量的某个线性组合(Linear Combination)的贡献。

在Kernel PCA分析之中,我们同样需要这些假设,但不同的地方是我们认为原有数据有更高的维数,我们可以在更高维的空间(Hilbert Space)中做PCA分析(即在更高维空间里,把原始数据向不同的方向投影)。这样做的优点有:对于在通常线性空间难于线性分类的数据点,我们有可能再更高维度上找到合适的高维线性分类平面。我们第二部分的例子就说明了这一点。

本文写作的动机是因为作者没有找到一篇好的文章(看了wikipedia和若干google结果后)深层次介绍PCA和Kernel PCA之间的联系,以及如何以公式形式来解释如何利用Kernel PCA来做投影,特别有些图片的例子只是展示了结果和一些公式,这里面具体的过程并没有涉及。希望这篇文章能做出较好的解答。

1. Kernel Principal Component Analysis 的矩阵基础

我们从解决这几个问题入手:传统的PCA如何做?在高维空间里的PCA应该如何做?如何用Kernel Trick在高维空间做PCA?如何在主成分方向上投影?如何Centering 高维空间的数据?

1.1 传统的PCA如何做?

让我先定义如下变量: \( X= [x_1, x_2, \ldots, x_N] \) 是一个\( d \times N \)矩阵,代表输入的数据有\( N\) 个,每个sample的维数是\( d \)。我们做降维,就是想用\( k \)维的数据来表示原始的\(d\)维数据(\( k \le d \))。
当我们使用centered的数据(即\(\sum_i x_i = 0\))时,可定义协方差矩阵\(C\)为:

$$C=\frac{1}{N} x_i x_i^T = \frac{1}{N} X X^T$$
做特征值分解,我们可以得到:
$$CU = U \Lambda \Rightarrow C = U \Lambda U^T = \sum_a \lambda_a u_a u_a^T$$
注意这里的\(C, U , \Lambda\)的维数都是\(d \times d\), 且\( U=[u_1, u_2, \ldots, u_d] \), \(\Lambda = diag(\lambda_1, \lambda_2, \ldots, \lambda_d) \)。
当我们做降维时,可以利用前\(k\)个特征向量\(U_k = [u_1, u_2, \ldots, u_k] \)。则将一个\(d\)维的\(x_i\)向\(k\)维的主成分的方向投影后的\(y_i = U_k^T x_i\) (这里的每一个\(u_i\)都是\(d\)维的,代表是一个投影方向,且\(u_i^T u_i = 1\),表示这是一个旋转变量)

1.2 在高维空间里的PCA应该如何做?

高维空间中,我们定义一个映射\(\Phi : X^d \rightarrow \mathcal{F}\),这里\(\mathcal{F}\)表示Hilbert泛函空间。
现在我们的输入数据是\(\Phi(x_i), i = 1, 2, …n\), 他们的维数可以说是无穷维的(泛函空间)。
在这个新的空间中,假设协方差矩阵同样是centered,我们的协方差矩阵为:
$$\bar{C}=\frac{1}{N} \Phi(x_i) \Phi(x_i)^T = \frac{1}{N} \Phi(X) \Phi(X)^T$$
这里有一个陷阱,我跳进去过:
在对Kernel trick一知半解的时候,我们常常从形式上认为\(\bar{C} \)可以用\(K_{i,j} = K(x_i, x_j)\)来代替,
因此对\(K=(K_{ij})\)做特征值分解,然后得到\(K=U \Lambda U^T\),并且对原有数据降维的时候,定义\(Y_i=U_k^T X_i\)。
但这个错误的方法有两个问题:一是我们不知道矩阵\(\bar{C}\)的维数;二是\(U_k^T X_i\)从形式上看不出是从高维空间的\(\Phi(X_i)\)投影,并且当有新的数据时,我们无法从理论上理解\( U_k^T X_{\mathrm{new}} \)是从高维空间的投影。
如果应用这种错误的方法,我们有可能得到看起来差不多正确的结果,但本质上这是错误的。
正确的方法是通过Kernel trick将PCA投影的过程通过内积的形式表达出来,详细见1.3

1.3 如何用Kernel Trick在高维空间做PCA?

在1.1节中,通过PCA,我们得到了\(U\)矩阵。这里将介绍如何仅利用内积的概念来计算传统的PCA。
首先我们证明\(U\)可以由\(x_1, x_2, \ldots, x_N\)展开(span):
$$ C u_a = \lambda_a u_a $$
$$ \begin{align}
u_a &= \frac{1}{\lambda_a} C u \\
&= \frac{1}{\lambda_a} (\sum_i x_i x_i^T )u \\
&= \frac{1}{\lambda_a} \sum_i x_i (x_i^T u) \\
&= \frac{1}{\lambda_a} \sum_i (x_i^T u) x_i \\
&= \sum_i \frac{x_i^T u}{\lambda_a} x_i \\
&= \sum_i \alpha_i^a x_i
\end{align}$$
这里定义\(\alpha_i^a=\frac{x_i^T u}{\lambda_a} \)。
因为\(x_i^T u\) 是一个标量(scala),所以\(\alpha_i^a\)也是一个标量,因此\(u_i\) 是可以由\(x_i\)张成。

进而我们显示PCA投影可以用内积运算表示,例如我们把\(x_i\)向任意一个主成分分量\(u_a\)进行投影,得到的是\(u_a^T x_i\),也就是\(x_i^T u_a\) 。作者猜测写成这种形式是为了能抽出\(x_i^T x_j = \mathrm{<} x_i, x_j\mathrm{>} \)的内积形式。

$$\begin{align}
x_i^T C u_a &= \lambda_a x_i^T u_a \\
x_i^T \frac{1}{N} \sum_j x_j x_j^T \sum_k \alpha_k^a x_k &= \lambda_a x_i^T \sum_k \alpha_k^a x_k \\
\sum_j \alpha_k^a \sum_k (x_i^T x_j) (x_j^T x_k) &= N \lambda_a \sum_k \alpha_k^a (x_i^T x_k) \\
\end{align}$$
当我们定义\( K_{ij} = x_i^T x_j \)时,上式可以写为\( K^2 \alpha = N \lambda_a K \alpha^a\)
(这里\(\alpha^a\)定义为\([\alpha_1^a, \alpha_2^a, \ldots, \alpha_N^a]^T\).)
进一步,我们得到解为:
$$K \alpha = \tilde{\lambda}_a \alpha \quad \mathrm{with} \quad \tilde{\lambda}_a = N \lambda_a$$
\(K\)矩阵包含特征值\(\tilde{\lambda}\)和\(\alpha^a\),我们可以通过\(\alpha\)可以计算得到\(u_a\),
注意特征值分解时Eigendecomposition,\(\alpha^a\)只代表一个方向,它的长度一般为1,但在此处不为1。
这里计算出\(\alpha_a\)的长度(下面将要用到):
因为\(u_a\)的长度是1,我们有:
$$\begin{align}
1 &= u_a^T u_a \\
&= ( \sum_i \alpha_i^a x_i) ^T (\sum_j \alpha_j^a x_j ) \\
&= \sum_i \sum_j \alpha_i^a \alpha_j^a x_i^T x_j^T \\
&=(\alpha^a)^T K \alpha_a \\
&=(\alpha^a)^T ( N \lambda_a \alpha_a) \\
&=N \lambda_a ({\alpha^a}^T \alpha_a )\\
&\Rightarrow \quad \lVert \alpha^a \rVert = 1/\sqrt{N \lambda_a} = 1/\sqrt{\tilde{\lambda}_a}
\end{align}$$

在上面的分析过程中,我们只使用了内积。因此当我们把\(K_{ij} = x_i^T x_j\)推广为\(K_{ij} = <\Phi(x_i), \Phi(x_j> = \Phi(x_i)^T \Phi(x_j)\)时,上面的分析结果并不会改变。

1.4 如何在主成分方向上投影?

投影时,只需要使用\(U\)矩阵,假设我们得到的新数据为\(t\),那么\(t\)在\(u_a\)方向的投影是:
$$u_a^T t = \sum_i \alpha_i^a x_i^T t = \sum_i \alpha_i^a ( x_i^T t)$$
对于高维空间的数据\(\Phi(x_i),\Phi(t)\),我们可以用Kernel trick,用\(K(x_i, t)\)来带入上式:
$$u_a^T t = \sum_i \alpha_i^a K(x_i, t)$$

1.5 如何Centering 高维空间的数据?

在我们的分析中,协方差矩阵的定义需要centered data。在高维空间中,显式的将\(\Phi(x_i)\)居中并不简单,
因为我们并不知道\(\Phi\)的显示表达。但从上面两节可以看出,所有的计算只和\(K\)矩阵有关。具体计算如下:
令\(\Phi_i = \Phi(x_i)\),居中\(\Phi_i^C = \Phi_i – \frac{1}{N}\sum_k \Phi_k\)
$$\begin{align}
K_{ij}^C
&= <\Phi_i^C \Phi_j^C> \\
&= (\Phi_i – \frac{1}{N}\sum_k \Phi_k)^T (\Phi_j – \frac{1}{N}\sum_l \Phi_l ) \\
&=\Phi_i^T\Phi_j – \frac{1}{N}\sum_l \Phi_i^T \Phi_l – \frac{1}{N}\sum_k \Phi_k^T \Phi_j + \frac{1}{N^2}\sum_k \sum_l \Phi_k^T \Phi_l \\
&=K_{ij} – \frac{1}{N}\sum_l K_{il} – \frac{1}{N}\sum_k K_{kj} + \frac{1}{N^2}\sum_k \sum_l K_{kl}
\end{align}$$
不难看出,
$$K^C = K – 1_N K – K 1_N + 1_N K 1_N$$
其中\(1_N\) 为\(N \times N \)的矩阵,其中每一个元素都是\(1/N\)
对于新的数据,我们同样可以
$$\begin{align}
K(x_i, t)^C
&= <\Phi_i^C \Phi_t^C> \\
&= (\Phi_i – \frac{1}{N}\sum_k \Phi_k)^T (\Phi_t – \frac{1}{N}\sum_l \Phi_l ) \\
&=\Phi_i^T\Phi_t – \frac{1}{N}\sum_l \Phi_i^T \Phi_l – \frac{1}{N}\sum_k \Phi_k^T \Phi_t + \frac{1}{N^2}\sum_k \sum_l \Phi_k^T \Phi_l \\
&=K(x_i,t) – \frac{1}{N}\sum_l K_{il} – \frac{1}{N}\sum_k K(x_k,t) + \frac{1}{N^2}\sum_k \sum_l K_{kl}
\end{align}$$

2. 演示 (R code)


首先我们应该注意输入数据的格式,一般在统计中,我们要求\(X\)矩阵是\(N \times d\)的,但在我们的推导中,\(X\)矩阵是\(d \times N\)。
这与统计上的概念并不矛盾:在前面的定义下协方差矩阵为\(X^T X\),而在后面的定义中是\(X X^T\)。另外这里的协方差矩阵是样本(Sample)的协方差矩阵,我们的认为大写的X代表矩阵,而不是代表一个随机变量。
另外,在下面的结果中,Gaussian 核函数(kernel function)的标准差(sd)为2。在其他取值条件下,所得到的图像是不同的。

KPCA图片:

R 源代码(Source Code):链接到完整的代码 KernelPCA

Kernel PCA部分代码:

# Kernel PCA
# Polynomial Kernel
# k(x,y) = t(x) %*% y + 1
k1 = function (x,y) { (x[1] * y[1] + x[2] * y[2] + 1)^2 }
K = matrix(0, ncol = N_total, nrow = N_total)
for (i in 1:N_total) {
  for (j in 1:N_total) {
    K[i,j] = k1(X[i,], X[j,])
}}
ones = 1/N_total* matrix(1, N_total, N_total)
K_norm = K - ones %*% K - K %*% ones + ones %*% K %*% ones
res = eigen(K_norm)

V = res$vectors
D = diag(res$values)

rank = 0
for (i in 1:N_total) {
	if (D[i,i] < 1e-6) { break }
      V[,i] = V[,i] / sqrt (D[i,i])
	rank = rank + 1
}
Y = K_norm %*%  V[,1:rank]
plot(Y[,1], Y[,2], col = rainbow(3)[label], main = "Kernel PCA (Poly)"
, xlab="First component", ylab="Second component")

3. 主要参考资料


【1】A Tutorial on Principal Component Analysis ,Jonathon Shlens, Shlens03

【2】Wikipedia: http://en.wikipedia.org/wiki/Kernel_principal_component_analysis

【3】 Original KPCA Paper:Kernel principal component analysis,Bernhard Schölkopf, Alexander Smola and Klaus-Robert Müller http://www.springerlink.com/content/w0t1756772h41872/fulltext.pdf

【4】Max Wellings’s classes notes for machine learning Kernel Principal Component Analaysis http://www.ics.uci.edu/~welling/classnotes/papers_class/Kernel-PCA.pdf