函数型数据主成分分析

函数型数据主成分分析

基本思想

函数型主成分分析(FPCA,Functional Principal Components Analysis)是传统的PCA的一种推广。

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

我们管各个函数到 β \beta β上的投影叫做观测曲线的主成分得分:
ξ i = ∫ T β ( s ) x i ( s ) d s , i = 1 , ⋯   , n \xi_{i}=\int_{\mathcal{T}} \beta(s) x_{i}(s) d s, \quad i=1, \cdots, n ξi=Tβ(s)xi(s)ds,i=1,,n

故而,求解第一个主成分就变成了求解一个优化问题:
max ⁡ 1 n ∑ i = 1 n ξ i 2 = max ⁡ 1 n ∑ i = 1 n ( ∫ T β ( s ) x i ( s ) d s ) 2  s.t.  ∥ β ∥ 2 = ∫ T β ( s ) β ( s ) d s = 1 \begin{aligned} \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) d s\right)^{2} \\ \text { s.t. } &\|\beta\|^{2}=\int_{T} \beta(s) \beta(s) d s=1 \end{aligned} maxn1i=1nξi2 s.t. =maxn1i=1n(Tβ(s)xi(s)ds)2β2=Tβ(s)β(s)ds=1
求解这个优化问题,我们就得到了第一主成分 β 1 ( s ) \beta^1(s) β1(s)

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

max ⁡ 1 n ∑ i = 1 n ξ i 2 = max ⁡ 1 n ∑ i = 1 n ( ∫ T β ( s ) x i ( s ) d s ) 2  s.t.  ∥ β ∥ 2 = ∫ T β ( s ) β ( s ) d s = 1 ∫ T β ( s ) β l ( s ) d s = 0 , l = 1 , ⋯   , 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) d s\right)^{2}} \\ {\text { s.t. }\|\beta\|^{2}=\int_{T} \beta(s) \beta(s) d s=1} \\ {\int_{T} \beta(s) \beta^{l}(s) d s=0, l=1, \cdots, k-1}\end{array} maxn1i=1nξi2=maxn1i=1n(Tβ(s)xi(s)ds)2 s.t. β2=Tβ(s)β(s)ds=1Tβ(s)βl(s)ds=0,l=1,,k1

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

v ( s , t ) = 1 n − 1 ∑ i = 1 n ( x i ( s ) − x ‾ ( s ) ) ( x i ( t ) − x ‾ ( t ) ) 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) v(s,t)=n11i=1n(xi(s)x(s))(xi(t)x(t))

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

∫ T v ( s , t ) β ( t ) d t = λ β ( s ) \int_{\mathcal{T}} v(s, t) \beta(t) d t=\lambda \beta(s) Tv(s,t)β(t)dt=λβ(s)

定义积分变换:

V β ( s ) = ∫ T v ( s , t ) β ( t ) d t V \beta(s)=\int_{\mathcal{T}} v(s, t) \beta(t) d t Vβ(s)=Tv(s,t)β(t)dt

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

V β ( s ) = λ β ( s ) V \beta(s)=\lambda \beta(s) Vβ(s)=λβ(s)

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

F V E = ∑ i = 1 K λ i / ∑ i = 1 n − 1 λ i \mathrm{FVE}=\sum_{i=1}^{K} \lambda_{i} / \sum_{i=1}^{n-1} \lambda_{i} FVE=i=1Kλi/i=1n1λi

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

特征问题求解

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

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

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

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

我们的样本基函数 x i x_i xi可以通过基函数展开,如下:
X i ( s ) = ∑ k = 1 K c i k Φ k ( s ) , i = 1 , 2 , … , N X_{i}(s)=\sum_{k=1}^{K} c_{i k} \Phi_{k}(s), i=1,2, \ldots, N Xi(s)=k=1KcikΦk(s),i=1,2,,N

我们记 X = ( x 1 , x 2 , … , x N ) ′ , Φ = ( Φ 1 , … , Φ k ) ′ , C = ( c i k ) N × K X=\left(x_{1}, x_{2}, \ldots, x_{N}\right)^{\prime}, \Phi=\left(\Phi_{1}, \ldots, \Phi_{k}\right)^{\prime}, C=\left(c_{i k}\right)_{N \times K} X=(x1,x2,,xN),Φ=(Φ1,,Φk),C=(cik)N×K

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

v ( s , t ) = 1 n − 1 Φ ′ ( s ) C ′ C Φ ( t ) v(s, t)=\frac{1}{n-1} \Phi^{\prime}(s) C^{\prime} C \Phi(t) v(s,t)=n11Φ(s)CCΦ(t)

定义K阶对称矩阵

W = ∫ Φ Φ ′ W=\int \Phi \Phi^{\prime} W=ΦΦ

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

同样地,将特征函数进行展开:

β ( s ) = ∑ k = 1 K b k Φ k ( s ) = Φ ′ ( s ) b \beta(s)=\sum_{k=1}^{K} b_{k} \Phi_{k}(s)=\Phi^{\prime}(s) b β(s)=k=1KbkΦk(s)=Φ(s)b

将其代入

∫ T v ( s , t ) β ( t ) d t = λ β ( s ) \int_{\mathcal{T}} v(s, t) \beta(t) d t =\lambda \beta(s) Tv(s,t)β(t)dt=λβ(s)

就可以得到( N = n − 1 N=n-1 N=n1):

1 N Φ ′ ( s ) C ′ C W b = λ Φ ′ ( s ) b \frac{1}{N} \Phi^{\prime}(s) C^{\prime} C W b=\lambda \Phi^{\prime}(s) b N1Φ(s)CCWb=λΦ(s)b

进一步能得到 1 N C ′ C W b = λ b \frac{1}{N} C^{\prime} C W b=\lambda b N1CCWb=λb,由特征向量正交和单位长度的约束要求,有 b k ′ W b k = 1 , b k ′ W b m = 0 , k ≠ m b_{k}^{\prime} W b_{k}=1, b_{k}^{\prime} W b_{m}=0,k \neq m bkWbk=1,bkWbm=0,k̸=m

W W W做cholesky分解,可得 W = L L ′ W=LL' W=LL

定义 u = L ′ b u=L'b u=Lb,那么上述问题就变成了对称矩阵的代数特征值问题:

1 N L ′ C ′ C L u = λ u \frac{1}{N} L' C^{\prime} C L u=\lambda u N1LCCLu=λu

据此可以求得 u u u,进而求得 b b b,最后求得特征函数 β \beta β

关于基函数的选取

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

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

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

B j , 0 ( x ) = { 1 , t j ≤ x &lt; t j + 1 0 , else B i , k ( x ) = x − t i t i + k − t i B i , k − 1 ( x ) + t i + k + 1 − x t i + k + 1 − t i + 1 B i + 1 , k − 1 ( x ) , k &gt; 0 \begin{array}{c}{B_{j, 0}(x)=\left\{\begin{array}{l}{1, t_{j} \leq x&lt;t_{j+1}} \\ {0, \text {else}}\end{array}\right.} \\ {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&gt;0}\end{array} Bj,0(x)={1,tjx<tj+10,elseBi,k(x)=ti+ktixtiBi,k1(x)+ti+k+1ti+1ti+k+1xBi+1,k1(x),k>0

程序代码

下面附一段简单的以多项式为基的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)
©️2020 CSDN 皮肤主题: 代码科技 设计师:Amelia_0503 返回首页