函数型数据分析
Published:
函数型数据分析针对于数据为函数的情况，而函数可以看作无穷维向量，因此函数型数据分析是有限维模型向无限维模型的推广。本文可以作为函数型数据分析的入门级导论，用于展现函数型数据分析的某些特质。
Functional PCA
PCA是一种常见的降维手段，在向量型数据中，PCA的任务是找到主要的几个方向向量，令数据在该方向上投影的误差最小，或者说方差最大，也即信息最大。对于PCA不熟悉的读者可以参考 PCA , 其中介绍了PCA，核PCA，概率PCA等PCA相关的降维方法。
Maximize Variance
在函数型数据中，协方差矩阵是一个无穷维的矩阵，更准确的说法其是一个核函数，
\[\sigma(s,t) = \sum_{i=1}^N X_i(s) X_i(t)\]从最大化方差的角度来理解PCA，将PCA中主成分的寻找看作一个类似贪心的过程，也即从最大的主成分开始，寻找第二个主成分，依次进行寻找，第一个主成分的方程可以写为（其中$\Sigma$为样本协方差矩阵）：
\[\begin{align} \max w^T \Sigma w ,\text{ s.t.} w^T w = 1 \end{align}\]利用Lagrange乘子法可以给出上述问题的解，
\[\Sigma w = \lambda w\]也即$u$是协方差矩阵$\Sigma$的特征向量，使得原问题最大则需要选取最大的特征值所对应的特征向量，
在给定第一个主成分的基础上，第二个主成分需要在于第一个主成分正交的前提下最大化方差，
\[\max v^T \Sigma v, \text{ s.t. } v^Tv=1, v^Tw =0\]同样可以使用Lagrange乘子法，
\[\begin{align} L &= \frac{1}{2} v^T \Sigma v- \frac{1}{2} \lambda(v^Tv-1) - \mu(v^Tw) \\ \frac{dL}{dv} &= \Sigma v - \lambda v - \mu w = 0 \\ \end{align}\]有趣的是Lagrange乘子$\mu$实际上为0，
\[\begin{align} \Sigma v &= \lambda v + \mu w \\ w^T \Sigma v &= \lambda w^Tv + \mu w^Tw, \\ \rho w^T v &= \lambda w^T v + \mu w^Tw,\text{With }\Sigma w = \rho w \\ \mu &= 0, \text{With } w^T v = 0 \end{align}\]因此可以得到第二个主成分也是一个特征向量，
\[\Sigma v = \lambda v\]在正交性的前提下，第二个主成分就是第二大特征值所对应的特征向量，逐步递推可以求解所有的$k$个主成分。
下面我们将上述的过程推广到函数型PCA上，对矩阵的梯度替换为对函数的梯度即可，
此时的主成分为一个函数，对于第一个主成分需要最大化其方差，
\[\max \iint w(t) \sigma(s,t) w(s) dsdt ,\text{ s.t.} \int w(s)^2 ds = 1\]利用Lagrange乘子法并且对函数求梯度可以得到和普通PCA一致的结果，
\[\int \sigma(s,t) w(t) dt = \lambda w(s)\]对于其他的主成分，加入正交约束可以类似地求解，由于和普通PCA类似，此处暂不赘述。
Minimize Reconstruction Error
更有意思的角度是从最小化重构误差的角度理解PCA，此时的PCA类似于一个自编码器（Auto-Encoder）。
对于普通的PCA，根据 Eckart-Young-Mirsky 定理，下式的结果由PCA给出，
\[\begin{align} \min \sum_i \Vert X_i- (Z Q)_i \Vert_2^2 , \text{s.t.} Q^T Q = I_k \end{align}\]而在 Eckart-Young-Mirsky 定理及其泛函版本简证 中，上式结果可以推广为，
\[\begin{align} \min \sum_i \int (X_i(t) - \sum_{j=1}^K S_{ij} \phi_j(t))^2 dt ,\text{s.t.} \int \phi_j(t) \phi_k(t) dt = \delta_{jk} \end{align}\]也即在寻找$K$个正交基的前提下，最小化重构误差的结果由PCA给出
该角度相较于最大化方差的角度是更有意思的，因为最大化方差的角度只能逐次求解主成分，相当于一种贪心的方法。而最小化重构误差的角度可以同时求出所有的主成分，并且其相当于证明了最大化方差中所给出的贪心算法实际上可以找到全局最优解。
Computing Functional PCA
对于函数型数据，关键环节在于如何求解一个无限维的问题，其中一个关键的思路是使用基函数展开，将一个函数在有限个基函数下进行展开，这个基函数可以是傅里叶基函数，也可以是B-样条基函数等。由于在有限个基函数下展开产生的误差，相较于过程中需要计算的数值积分、矩阵运算等的误差相比可能是微不足道的，因此选取有限组基展开并不会影响问题的求解。
将函数型PCA中涉及到的函数都在一组基函数$\phi(t) = [\phi_1(t),…\phi_K(t)]^T$下展开，
\(\begin{align} w(s) &= \sum_{j=1}^K\phi_j(s) b_j = \phi(s)^T b = b^T \phi(s)\\ \sigma(s,t) &= \sum_{i=1}^N X_i(s)X_i(t) = X^T(s)X(t) = \phi^T(s) C^TC \phi(s) \end{align}\) 上式中的$C,b$代表基前面的系数，利用上述关系式对其求解的等式进行化简，
\[\begin{align} \int \sigma(s,t) w(t) dt &= \lambda w(s) ,\text{ s.t.} \int w(s)^2 ds =1 \\ \int \phi^T(s) C^T C \phi(t) \phi(t)^T b dt &= \lambda \phi(s)^T b ,\text{ s.t. } \int b^T \phi(s) \phi(s)^T bds=1 \\ \phi^T(s) \int C^T C \phi(t) \phi(t)^T b dt &= \lambda \phi(s)^T b ,\text{ s.t. } \int b^T \phi(s) \phi(s)^T bds=1 \\ \int C^T C \phi(t) \phi(t)^T b dt &= \lambda b ,\text{ s.t. } \int b^T \phi(s) \phi(s)^T bds=1 \\ C^T C \int \phi(t) \phi(t)^T dt b &= \lambda b ,\text{ s.t. } b^T \int\phi(s) \phi(s)^T dsb=1 \\ C^TC W b &=\lambda b, \text{ s.t. } b^TWb = 1,\text{ Let } W = \int\phi(s) \phi(s)^Tds \end{align}\]定义上述的矩阵$W$之后，可以将式子化简为简单的矩阵表达式，而由于选择的基函数$\phi(s)$是有限个，因此其内积构成的矩阵$W$也是有限维的矩阵，上述问题可以转化为一个特征值问题求解。由于$W$为内积构成的矩阵，其为正定矩阵，因此可以将其分解为$W = W^{\frac{1}{2}}W^{\frac{1}{2}}$ ,基于上述分解并且利用换元，可以进行转化，
\[\begin{align} C^TC W b &=\lambda b, \text{ s.t. } b^TWb = 1 \\ W^{\frac{1}{2}}C^TC W^{\frac{1}{2}} W^{\frac{1}{2}} b &=\lambda W^{\frac{1}{2}}b, \text{ s.t. } b^TW^{\frac{1}{2}}W^{\frac{1}{2}}b = 1 \\ \tilde W u &= \lambda u,\text{ s.t.} u^T u=1 \\ \text{Let } \tilde W &= W^{\frac{1}{2}}C^TC W^{\frac{1}{2}}, u = W^{\frac{1}{2}}b \end{align}\]到此为止，我们将函数型PCA问题利用基函数展开的方式，转化为了一个常见的矩阵的特征值求解问题。
Functional Linear Regression
函数型线性回归模型是普通线性回归模型的推广，根据自变量和因变量是否为函数或向量，有不同的几个模型，
| 自变量 | 因变量 | 模型 |
|---|---|---|
| 向量 | 向量 | X - Y |
| 函数 | 向量 | X(t) - Y |
| 向量 | 函数 | X - Y(t) |
| 函数 | 函数 | X(t) - Y(t) |
可以看到第一种模型是最简单的模型，也即普通的线性回归模型，而其他三种模型是函数型线性回归模型所需要解决的问题。本节只关注于函数-向量线性回归模型，也即上表中的第二类模型。
不失一般性，我们假设因变量为一个标量，采用最小二乘法建立线性回归模型，
\[\min \mathcal{L} = \min \sum_i (y_i - \int \beta(t) x_i(t) )^2\]由于此时的系数$\beta(t)$为函数，上式的零解通常很容易取到，因此通常需要对$\beta(t)$做某些限制，例如限制其光滑性，这也是实际中所需要满足的性质，其他的限制还包括周期性质等。如果采用二阶梯度的$L_2$范数对光滑性进行约束，损失函数可以写为，
\[\min \mathcal{L} = \min \sum_i (y_i - \int \beta(t) x_i(t) )^2 + \lambda \int (\nabla^2\beta(t))^2 dt\]同样基于基函数展开的思想，将$\beta(t)$在一组基下进行展开，
\[\begin{align} \beta(t) & = \sum_{i=1}^K c_i \phi_i(t) = c^T \phi(t)\\ \nabla^2 \beta(t) &= \sum_{i=1}^K c_i \nabla^2\phi_i(t) = c^T \nabla^2 \phi(t) \end{align}\]此时的损失函数可以进行化简为矩阵形式，
\[\begin{align} &\min \sum_i (y_i - \int \beta(t) x_i(t) )^2 + \lambda \int (\nabla^2\beta(t))^2 dt \\ =& \sum_i (y_i - c^T \int \phi(t) x_i(t) )^2 + \lambda c^T\int \nabla^2\phi(t)\nabla^2\phi(t)^T dt c \\ =& \sum_i(y_i - c^T z_i)^2+ \lambda c^T R c \\ =& (y- Z c)^T(y-Zc) +\lambda c^T Rc \\ \text{With} & z_i = \int \phi(t) x_i(t) , R = \int \nabla^2\phi(t)\nabla^2\phi(t)^T dt \end{align}\]最后的等式和岭回归（Ridge Regression）是一致的， \(\begin{align} L &= (y- Z c)^T(y-Zc) +\lambda c^T Rc \\ \frac{dL}{dc} &= Z^T(Zc - y) +\lambda R c = 0 \end{align}\)
上式的解同样与岭回归的解类似，
\[\begin{align} Z^T(Zc-y) + \lambda Rc &= 0 \\ (Z^T Z +\lambda R)c &= Z^T y \\ (Z^T Z +\lambda R)^{-1} Z^T y &=c \end{align}\]限于篇幅，本节仅用该例子展现将函数型线性回归模型转化为普通线性回归模型的方法，对于其他类型的函数型线性回归模型的建立和求解此处暂略。