函数型数据主成分分析

函数型数据主成分分析

基本思想

函数型主成分分析（FPCA，Functional Principal Components Analysis）是传统的PCA的一种推广。

考虑我们已经从数据中得到拟合曲线$x_i(s),s\in \mathcal{T},i=1,\cdots ,n$，所谓的第一主成分，就是我们希望能找到一个模为1的函数$\beta (s)$，使得 { x i } \{x_i\} {xi​}在$\beta$上的投影（$L_2$內积） { ξ i } \{\xi _i\} {ξi​}的方差达到最大，方差最大其实也就体现 { x i } \{x_i\} {xi​}整体到$\beta$的距离达到最小。$\beta$一般就叫做权重函数（可以理解为“坐标轴”单位长度量）。

我们管各个函数到$\beta$上的投影叫做观测曲线的主成分得分：
$${\xi }_i={\int }_{\mathcal{T}}\beta (s)x_i(s)ds,i=1,\cdots ,n$$

故而，求解第一个主成分就变成了求解一个优化问题：
$$\begin{array}{rl}\max \frac{1}{n}\sum _{i=1}^n{\xi }_i^2& =\max \frac{1}{n}\sum _{i=1}^n{\left({\int }_{\mathcal{T}}\beta (s)x_i(s)ds\right)}^2\\ \text{ s.t. }& \Vert \beta {\Vert }^2={\int }_T\beta (s)\beta (s)ds=1\end{array}$$
求解这个优化问题，我们就得到了第一主成分${\beta }^1(s)$。

第$k$主成分无非就是在满足和前面$k-1$个主成分权重函数垂直的基础上，求解上述优化问题而已，即求解

$$\begin{array}{l}\max \frac{1}{n}{\sum }_{i=1}^n{\xi }_i^2=\max \frac{1}{n}{\sum }_{i=1}^n{\left({\int }_{\mathcal{T}}\beta (s)x_i(s)ds\right)}^2\\ \text{ s.t. }\Vert \beta {\Vert }^2={\int }_T\beta (s)\beta (s)ds=1\\ {\int }_T\beta (s){\beta }^l(s)ds=0,l=1,\cdots ,k-1\end{array}$$

这个优化问题的解可以表述如下。记协方差函数：

$$v(s,t)=\frac{1}{n-1}\sum _{i=1}^n\left(x_i(s)-\overline{x}(s)\right)\left(x_i(t)-\overline{x}(t)\right)$$

那么权重函数满足特征方程：

$${\int }_{\mathcal{T}}v(s,t)\beta (t)dt=\lambda \beta (s)$$

定义积分变换：

$$V\beta (s)={\int }_{\mathcal{T}}v(s,t)\beta (t)dt$$

这里的$V$称为协方差算子，它将函数$\beta$变成一个函数。那么，我们有：

$$V\beta (s)=\lambda \beta (s)$$

我们也类比PCA，使用特征值的累积贡献率来衡量主成分所占比例：

$$\mathrm{F}\mathrm{V}\mathrm{E}=\sum _{i=1}^K{\lambda }_i/\sum _{i=1}^{n-1}{\lambda }_i$$

这里之所以对$\lambda$只累计到$n$是因为协方差算子$V$的秩为样本数量减一个，则非零特征根的个数最多为$n-1$个。

特征问题求解

由上述已知，我们求解主成分最后归结为求解一个特征值问题。

求解这个问题，目前比较流行的有三种方法：

• 对函数进行SVD离散化
• 对函数进行基函数展开
• 运用一般性的数值积分方法

我们最后需要的是特征函数，为了避免插值而带来更大的误差，我选用对基函数进行展开的方法。下面简单介绍一个对函数进行基函数展开的基本思路。

我们的样本基函数$x_i$可以通过基函数展开，如下：
$$X_i(s)=\sum _{k=1}^K{c}_{ik}{\Phi }_k(s),i=1,2,\dots ,N$$

我们记$$X={\left(x_1,x_2,\dots ,x_N\right)}^{\prime },\Phi ={\left({\Phi }_1,\dots ,{\Phi }_k\right)}^{\prime },C={\left(c_{ik}\right)}_{N\times K}$$

那么样本函数就可以写为等价的矩阵形式$X=C\Phi$。那么协方差函数就可以写为（假设已经标准化）：

$$v(s,t)=\frac{1}{n-1}{\Phi }^{\prime }(s)C^{\prime }C\Phi (t)$$

定义K阶对称矩阵

$$W=\int \Phi {\Phi }^{\prime }$$

当选择正交基的时候，比如说正交傅里叶基，这就是一个单位矩阵。关于这个基如何选取，我们后面还会详谈。

同样地，将特征函数进行展开：

$$\beta (s)=\sum _{k=1}^K{b}_k{\Phi }_k(s)={\Phi }^{\prime }(s)b$$

将其代入

$${\int }_{\mathcal{T}}v(s,t)\beta (t)dt=\lambda \beta (s)$$

就可以得到（$N=n-1$）：

$$\frac{1}{N}{\Phi }^{\prime }(s)C^{\prime }CWb=\lambda {\Phi }^{\prime }(s)b$$

进一步能得到$\frac{1}{N}{C}^{\prime }CWb=\lambda b$，由特征向量正交和单位长度的约束要求，有$$b_k^{\prime }W{b}_k=1,b_k^{\prime }W{b}_m=0,k̸=m$$

对$W$做cholesky分解，可得$W=L{L}^{\prime }$。

定义$u=L^{\prime }b$，那么上述问题就变成了对称矩阵的代数特征值问题：

$$\frac{1}{N}{L}^{\prime }{C}^{\prime }CLu=\lambda u$$

据此可以求得$u$，进而求得$b$，最后求得特征函数$\beta$。

关于基函数的选取

如上所述，可以看出，如果所选的基函数是正交的，本质上和PCA的以拟合系数为坐标点的函数空间PCA推广是实际上是一样的。若基函数不是正交的，无非就是在此基础上对要求特征值的矩阵得多乘一个$W=\int \Phi {\Phi }^{\prime }$，再求特征向量，以及进行$W$意义下对特征向量进行单位化而已（不单位话也没事，只不过权重函数$\beta$不再是模长为1的而已）。

常用的基函数有傅里叶基函数和B样条基函数，傅里叶基函数适用于周期性函数数据，B样条基函数适用于非周期函数数据，当然，也可以用多项式基函数。

B样条基函数的递归定义为：

$$\begin{array}{c}{B}_{j,0}(x)=\{\begin{array}{l}1,t_j\le x<t_{j+1}\\ 0,\text{else}\end{array}\\ B_{i,k}(x)=\frac{x-t_i}{t_{i+k}-t_i}{B}_{i,k-1}(x)+\frac{t_{i+k+1}-x}{t_{i+k+1}-t_{i+1}}{B}_{i+1,k-1}(x),k>0\end{array}$$

程序代码

下面附一段简单的以多项式为基的matlab代码。

clc
clear
close all
format short
digits(5)
%% 设定维数
p = 8;
d = 2;
%% 读入数据，保存
D = dir('../data/*.csv');
D = struct2cell(D);%结构体不好调用，转成元胞数组来使用
csvs_name = D(1,:);
N = 0;%这个表示数据的个数
%csvs_name = csvs_name(5:6);
for csv_name = csvs_name %cpicName = 'a001.bmp'
    N = N+1;
    csvName = csv_name{1};
    Data{N} = load(['../data/' csvName]);
end
%clear csv_name csvN csvName csvs_name D
%% 定义原空间基，生成模型函数，由模型函数通过拟合得到点的坐标表示
%p = 15;%原空间维数
%%根据p自动生成模型函数
eval_poly = 'a(1)';
for i = 2:p
    eval_poly = [eval_poly '+a(' num2str(i) ').*x.^' num2str(i-1)];
end
eval_poly = ['modelfun = @(a,x) x.^(0.5).*(1-x).^(0.5).*(' eval_poly ');'];
eval(eval_poly);

%modelfun = @(a,x) x.^(0.5).*(1-x).^(0.5).*(a(1)+a(2).*sin(x)+a(3).*cos(x)+a(4).*1./x+a(5).*x+a(6).*x.^2);
%   modelfun = @(a,x)(x.^(0.5).*(1-x).^(0.5).*(a(1)+a(2).*x+a(3).*x.^2+...
%       a(4).*x.^3))%+a(5).*x.^4+a(6)*x.^5))+...
%       a(7)*x.^6+a(8).*x.^7+a(9).*x.^8+a(10).*x.^9));
%         modelfun = @(a,x)((a(1)+a(2).*x.^2+a(3).*x.^2+...
%         a(4).*x.^3+a(5).*x.^4+a(6)*x.^5+...
%         a(7)*x.^6+a(8).*x.^7+a(9).*x.^8+a(10).*x.^9));

%%做一波拟合，求出系数，进行拟合，拟合结果为A
%%A中存的是数据点的坐标表示
for k = 1:N
    data = Data{k};
    data(:,2) = abs(data(:,2));%取绝对值
    x = data(:,1);
    y = data(:,2);
    %y = y./(x.*(1-x));%%%%%%%%如果将x(1-x)项挪到左边，拟合会出问题
    a0 = ones(1,p);
    a = nlinfit(x,y,modelfun,a0);
    A(:,k) = a';
end
%clearvars -except A modelfun x y

%% 开始奇异值分解降维，得到新空间的原点和坐标轴，坐标表示
%%由奇异值分解做主成分分析，A的每一列为原空间下的系数，U的每一列为主成分
[p,N] = size(A);
%d = 5;%新空间维数
x0 = mean(A,2);%计算样本均值
A_shift = A - repmat(x0,1,N);%中心平移每个训练样本，A_shift为原数据点的平移表示

%%%%对于A_shift的补充处理
syms x;
base_vec_T = x.^(0.5).*(1-x).^(0.5).*x.^[0:p-1];
base_vec = conj((base_vec_T))';
W = int(base_vec*base_vec_T,[0,1]);
L_prime = chol(W);
L_prime_A_shift = L_prime*A_shift;

%%%%
[U,S,V] = svd(L_prime_A_shift);
B = L_prime\U;
if size(S,2) == 1
    S_diag = S(1);
else
    S_diag = diag(S);
end
R = cumsum(S_diag)./sum(S_diag);
Coord_new = B(:,1:d);


%% 新空间的坐标轴和原点的函数表示
%%生成新的空间中的基函数，找到d个新空间基函数及原点坐标
%%Coord_new x0中存的是新坐标系，fi存的是新的函数空间的原点和坐标轴
for k = 1:d;%取出第k个基函数
    syms x;
    c = B(:,k);
    exp = modelfun(c,x);
    exp = simplify(exp);%第k个基函数的表达式
    f{k} = exp;
    digits(2);
    latex(vpa(exp,2));%写成latex表达式黏贴到解释器中
end

f0 = modelfun(x0,x);
f0 = simplify(f0);
latex(vpa(f0,2));
% for i = 1:d
%     f_fun{i} = matlabFunction(f{i});
% end
% f0_fun = matlabFunction(f0);

%% 求原坐标系中的投影后在原空间中的位置



%A_new_coord = Coord_new'*A_shift;%A_new_coord表示原数据点在新空间中的表示

syms x;
for k = 1:N
    a = A_shift(:,k);
    pk_f = modelfun(a,x);
    for j = 1:d
        A_new_coord(j,k) = int(pk_f.*f{j},0,1);
    end
end

%%函数运算
for k = 1:N
    A_shadow_f = f0;
    for i = 1:d
        A_shadow_f = A_shadow_f+f{i}*A_new_coord(i,k);
    end
    A_shadow_f = simplify(A_shadow_f);
    A_shadow_fs{k} = A_shadow_f;
    A_shadow_fs_fun{k} = matlabFunction(A_shadow_f);
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%坐标运算，可以注释掉
% A_shadow_coord = Coord_new*A_new_coord+repmat(x0,1,N);%A_shadow_coord表示投影坐标点在原坐标系下的表示
% syms x;
% for k = 1:N
%     a = A_shadow_coord(:,k);
%     A_shadow_coord_f = modelfun(a,x);
%     A_shadow_coord_fs{k} = simplify(A_shadow_coord_f);
% end
%
% %%不考虑舍入误差的情况下，可以证明坐标运算和函数运算的得到的函数表达式是一样的，验证代码如下
% k=5;
% diff = A_shadow_coord_fs{k}-A_shadow_fs{k}
% diff = simplify(diff)
% digits(2);
% latex(vpa(diff,2))
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 画图和求误差
Errors = [];
NN = 1000;
h = 1/NN;
xx = linspace(0,1,NN+1);
close all;
for k = 1:N
    %%原数据点
    data = Data{k};
    data(:,2) = abs(data(:,2));%取绝对值
    x = data(:,1);
    y = data(:,2);
    %%拟合曲线
    a = A(:,k);
    data_fit_f = @(x) modelfun(a,x);
    %%投影函数
    A_shadow_f_fun = A_shadow_fs_fun{k};
    %%原点函数
    f0_fun = matlabFunction(f0);
    
    %%开始画图
    figure(k);
    scatter(x,y,'.','MarkerEdgeColor','b',...
              'MarkerFaceColor','c',...
              'LineWidth',0.5);
    hold on;
    
    y_fit = data_fit_f(xx);
    y_shadow = A_shadow_f_fun(xx);
    y_f0 = f0_fun(xx);
    plot(xx,y_fit,'g');
    plot(xx,y_shadow,'black');
    plot(xx,y_f0,'r');
    legend('数据','拟合','投影','原点');
    
    %%求拟合曲线和投影函数之间的L2误差
    %Errors(end+1) = norm((y_fit-y_shadow),2);
    Errors(end+1) = sum(abs(y_fit-y_shadow)*h);
end

Errors_mean = mean(Errors)



%% test
% close all;
% plot(xx,xx);
% hold on;
% plot(xx,xx.^9)