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
angle_2^2, mathrm{s.t.} : |xi|_2=1.
]
其中(langle x, y
angle=int_I xy : mathrm{d}t).
假设:
[xi = sum_{i=1}^N v_i x_i.
]
并记:
[M in mathbb{R}^{N imes N}, M_{i,j}=langle x_i, y_j
angle_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_Ifmathrm{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, fge 0, tin 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 eta) odot f= alpha odot (eta odot f), \
(alpha + eta) odot f= (alpha odot f) oplus (eta odot f).
]
注意到:
令(g(t)=1/eta, tin 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
angle_{mathcal{B}} = frac{1}{2eta} 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}{2eta} 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}{2eta} 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}{2eta} int_I int_I ln^2 frac{f^{alpha}(t)}{f^{alpha}(s)} mathrm{d}{t} mathrm{d}s]^{1/2} = |alpha|[frac{1}{2eta} int_I int_I ln^2 frac{f(t)}{f(s)} mathrm{d}{t} mathrm{d}s]^{1/2} = |alpha||f|_{mathcal{B}}.
]
最后证其满足三角不等式:
[egin{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)
ightarrow 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
angle_{mathcal{B}}=langle f_c, g_c
angle_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, Nin 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
angle_{mathcal{B}}= langle mathrm{clr}(x_i), mathrm{clr}(x_J)
angle_2). 则有类似的公式:
[M^2v = lambda Mv, |Xv|_{mathcal{B}}=1.
]
转化为(L^2(I))上的PCA是类似的:
[mathrm{clr}(xi) = sum_i^N v_imathrm{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
angle_2), 则:
[M = C^TM'C
]
又
[mathrm{clr}(xi) = Phi Cv
]
令(b = Cv), 可得:
[Mv = C^TPhi^T Phi Cv = C^TPhi^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.)