Simplicial principal component analysis for density functions in Bayes spaces

目录

• 问题
• $\mathcal{B}^2(I)$
• $\mathcal{B}^2(I)$上的PCA

Hron K, Menafoglio A, Templ M, et al. Simplicial principal component analysis for density functions in Bayes spaces[J]. Computational Statistics & Data Analysis, 2016: 330-350.

问题

我们知道一般的PCA，其数据是$x \in \mathbb{R}^n$的，事实上，已经有很多关于函数类数据的PCA了.

一般的函数型PCA是定义在$L^2$空间上的. 假设$x_1, x_2, \ldots, x_N \in L^2(I)$, 并假设是中心化的. 我们希望找到一个$\xi$最大化:

$$\frac{1}{N} \sum_{i=1}^N \langle x_i, \xi \rangle_2^2, \mathrm{s.t.} \: \|\xi\|_2=1.$$

其中$\langle x, y \rangle=\int_I xy \: \mathrm{d}t$.
假设:

$$\xi = \sum_{i=1}^N v_i x_i.$$

并记:

$$M \in \mathbb{R}^{N \times N}, M_{i,j}=\langle x_i, y_j \rangle_2$$

则最初的式子可以表示为:

$$\frac{1}{N} v^TM^TMv, \quad \mathrm{s.t.} \: \|Xv\|_2=1.$$

可以证明，KKT条件为:

$$M^2v=\lambda Mv$$

显然，$v$是$M$的首特征向量(当然$\|v\|=1$不一定成立).
类似的，其它的载荷向量也是如此求得. 上面有一点存疑的地方是:

$$\xi = \sum_{i=1}^N v_i x_i.$$

在$\mathbb{R}^n$中是绝对没问题的是，问题是在$L^2$，是否可以分解一个元素呢? 可以的，绝对是可以的.

作者是将一般的函数的PCA，限定在密度函数的PCA，我们知道，密度函数$f$满足:

$$f \ge 0, \\ \int_If\mathrm{d}t=1.$$

显然$\xi = \sum_{i=1}^N v_i x_i$并不一定能够满足上面的性质，为此，作者引入了一个新的贝叶斯空间$\mathcal{B}^2(I)$.

$\mathcal{B}^2(I)$

假设$I=[a,b]$，我们的工作是构造一个空间，使得上面的元素其线性运算能够保持密度函数的性质.
首先说明，$\mathcal{B}^2(I)$里的元素为$\{f|\int_I f(t) \mathrm{d}t=1, f\ge 0, t\in I\}$.

记$\eta=b-a$，后续我们会发现，$1/\eta$是这个空间的零元素.

首先定义加法和数乘法，使其称为一个向量空间.

$$(f \oplus g) (t)=\frac{f(t)g(t)}{\int_If(s)g(s) \mathrm{d}s}, \quad t \in I,$$

可以发现$\oplus$是保持密度函数的性质的(只要$f,g$在$I$上满足).

$$(\alpha \odot f)(t)=\frac{f(t)^{\alpha}}{\int_I f(s)^{\alpha} \mathrm{d}s}, \quad t \in I,$$

显然也是保持的.
并且，容易证明(利用类似核方法的思想):

$$f \oplus g = g \oplus f, \\ f \oplus g \oplus h=f \oplus (g \oplus h), \\ \alpha \odot (f \oplus g) = (\alpha \odot f) \oplus (\alpha \odot g), \\ (\alpha \cdot \beta) \odot f= \alpha \odot (\beta \odot f), \\ (\alpha + \beta) \odot f= (\alpha \odot f) \oplus (\beta \odot f).$$

注意到:
令$g(t)=1/\eta, t\in I$

$$f \oplus g=f, \quad 0 \odot f = \frac{1}{\eta}$$

所以$1/\eta$是零元素，那么可以如此定义差:

$$f \ominus g= f \oplus [(-1) \odot g],$$

易得:

$$f \ominus f= 1 /\eta.$$

再定义内积，使其成为一个内积空间:

$$\langle f, g \rangle_{\mathcal{B}} = \frac{1}{2\eta} \int_I \int_I \ln \frac{f(t)}{f(s)} \ln \frac{g(t)}{g(s)} \mathrm{d}t \mathrm{d}s, \quad, f, g \in \mathcal{B}^2(I).$$

则，我们可以定义其上的范数为:

$$\|f\|_{\mathcal{B}} = [\frac{1}{2\eta} \int_I \int_I \ln^2 \frac{f(t)}{f(s)} \mathrm{d}{t} \mathrm{d}s]^{1/2}.$$

下证其为一范数:
非负性是显然的, 首先证明其是正定的，即，零元素的大小为0:

$$\|1/\eta\|_{\mathcal{B}} = [\frac{1}{2\eta} \int_I \int_I \ln^2 1 \mathrm{d}{t} \mathrm{d}s]^{1/2}=0.$$

其次，证明其是其次的，即$\|\alpha \odot f\|_{\mathcal{B}}=|\alpha|\|f\|_{\mathcal{B}}$:

$$\|\alpha \odot f\|_{\mathcal{B}} = [\frac{1}{2\eta} \int_I \int_I \ln^2 \frac{f^{\alpha}(t)}{f^{\alpha}(s)} \mathrm{d}{t} \mathrm{d}s]^{1/2} = |\alpha|[\frac{1}{2\eta} \int_I \int_I \ln^2 \frac{f(t)}{f(s)} \mathrm{d}{t} \mathrm{d}s]^{1/2} = |\alpha|\|f\|_{\mathcal{B}}.$$

最后证其满足三角不等式:

$$\begin{array}{ll} \|f \oplus g\|_{\mathcal{B}}&=[\frac{1}{2 \eta}\int_I \int_I \ln^2 \frac{f(t)g(t)}{f(s)g(s)}\mathrm{d}t \mathrm{d}s]^{1/2} = [\frac{1}{2 \eta}\int_I \int_I \ln^2 \frac{f(t)g(t)}{f(s)g(s)}\mathrm{d}t \mathrm{d}s]^{1/2}\\ &= [\frac{1}{2 \eta}\int_I \int_I \ln^2 \frac{f(t)}{f(s)}\mathrm{d}t \mathrm{d}s + \frac{1}{2 \eta}\int_I \int_I \ln^2 \frac{g(t)}{g(s)}\mathrm{d}t \mathrm{d}s]^{1/2} \\ & \le [\frac{1}{2 \eta}\int_I \int_I \ln^2 \frac{f(t)}{f(s)}\mathrm{d}t \mathrm{d}s]^{1/2} + [\frac{1}{2 \eta}\int_I \int_I \ln^2 \frac{g(t)}{g(s)}\mathrm{d}t \mathrm{d}s]^{1/2} \\ &= \|f\|_{\mathcal{B}}+\|g\|_{\mathcal{B}}. \end{array}$$

证毕.

定义一个$\mathcal{B}^2(I) \rightarrow L^2(I)$上的函数:

$$\mathrm{clr} (f)(t) = f_c(t) = \ln f(t) - \frac{1}{\eta} \int_I \ln f(s) \mathrm{d}s.$$

为什么要定义一个这样的函数等等再讲，先来看看它的性质——不仅仅是等距映射.

$$\mathrm{clr} (f \oplus g)(t)=f_c(t)+g_c(t), \quad \mathrm{clr} (\alpha \odot f)(t) =\alpha \cdot f_c(t), \quad \langle f, g \rangle_{\mathcal{B}}=\langle f_c, g_c \rangle_2=\int_I f_c(t) g_c(t) \mathrm{d}t.$$

这些性质的证明是容易的.

还需要注意的一个性质，不应该称之为限制条件才对:

$$\int_I f_c \mathrm{d}t=\int_I \ln f(t) \mathrm{d}t - \int_I \ln f(s) \mathrm{d}s=0.$$

这就意味着，只有$L^2(I)$中满足积分为0的函数才能在$\mathcal{B}^2(I)$中有原像.

接下来解释为什么要弄这样一个映射. 因为一般情况下，我们首先面对的都是一些离散的数据，然后利用某些方法进行拟合，比如论文中提到的$B-$样条，但是拟合出来的函数往往并不是密度函数，所以便有了$\mathrm{clr}$变化，这个变化可以帮助我们有效利用已有的函数，利用已有函数的积分等性质来应对$\mathcal{B}^2(I)$中的一些计算.
当然这也给函数逼近增加了难度，就是在区间$I$上积分和需要为1，这个问题在另一篇文章中进行了详细的讨论.

$\mathcal{B}^2(I)$上的PCA

假设$x_i, i=1,2,\ldots, N\in \mathcal{B}^2(I)$, 那么令:

$$\xi = \sum_{i=1}^N v_i \odot x_i = (v_1 \odot x_1) \oplus (v_2 \odot x_2) \oplus \cdots \oplus (v_N \odot x_N).$$

令矩阵$M$其元素$M_{ij}=\langle x_i, x_j \rangle_{\mathcal{B}}= \langle \mathrm{clr}(x_i), \mathrm{clr}(x_J) \rangle_2$. 则有类似的公式:

$$M^2v = \lambda Mv, \|Xv\|_{\mathcal{B}}=1.$$

转化为$L^2(I)$上的PCA是类似的:

$$\mathrm{clr}(\xi) = \sum_i^N v_i\mathrm{clr}(x_i),$$

$$M^2v = \lambda Mv, \|\mathrm{clr}(\xi)\|_2=1.$$

在实际情况中$\mathrm{clr}(x_i)$是通过函数逼近得到的，假设为:

$$\mathrm{clr}(x_i)=\Phi c_i, \Phi=[\phi_1, \ldots, \phi_K].$$

则

$$\mathrm{clr}(X)=\Phi C,$$

假设$M'_{ij} = \langle \phi_i, \phi_j \rangle_2$, 则:

$$M = C^TM'C$$

又

$$\mathrm{clr}(\xi) = \Phi Cv$$

令$b = Cv$, 可得:

$$Mv = C^T\Phi^T \Phi Cv = C^T\Phi^T \Phi b = \lambda v,$$

两边同乘以$C$可得:

$$CC^T \Phi^T \Phi b = \lambda b$$

解得$b$, 可知:

$$\mathrm{clr}(\xi) = \Phi b \Rightarrow \xi = \mathrm{clr}^{-1}(\Phi b).$$

注意: $\int_I \mathrm{clr}(x_i) \mathrm{d}t=0.$