论文笔记与复现[156]PARAFAC. tutorial and applications

原文下载：https://www.sciencedirect.com/science/article/abs/pii/S0169743997000324

摘要

本文介绍了PARAFAC的多维分解方法及其在化学计量学中的应用。PARAFAC是PCA向高阶数组的推广，但该方法的一些特性与普通的二维情况截然不同。例如，可以从多维光谱数据(multi-way spectral data)中恢复出纯光谱（pure spectra）。

1 介绍

以交叉方式测量变量，结果的集合为多维数据。

PARAFAC以及二路PCA等方法都是多线性或双线性分解方法，它们将数组分解成分数和负载[16]（loadings）的集合，希望以比原始数据数组更精简的形式描述数据。

主成分分析模型可以被认为是最复杂和最灵活的模型，而PARAFAC是最简单和最受限制的模型。

结构越多，拟合越差，模型越简单。使用多维方法不是为了获得更好的拟合，而是为了获得更充分、更稳健和可解释的模型。

对于组分数为F的I×J×K数组，平行因子模型含有F(I+J+K)个参数。

PARAFAC的一个非常令人讨厌的特性是计算模型所需的时间很长。所使用的算法通常基于交替最小二乘法（ALS），ALS的初始化使用随机值或基于广义特征值问题的直接三线性（trilinear）分解。

在下文中，为了简单起见，讨论将仅限于三维（three-way）数据，但大多数结果对任何（更高）阶的数据和模型都有效。

2 术语

标量：小写斜体

矢量：粗体小写

二维矩阵：粗体大写

三维数组：带下划线的粗体大写字母

xijk：X的第ijk个元素

模式（mode）、way和顺序（order）这三个术语或多或少可以互换使用。

术语因子（factor）和组分（component）之间没有区别。

3 模型

数据被分解为三线性分量（三元组，triads），每个分量由一个分数向量和两个负载向量组成。在三维中，通常不区分分数和负载（因为分数和负载在数学上是同等对待的）。

三维数组的平行因子模型由三个负载矩阵A、B、C组成（其中的元素分别表示为aif、bjf、ckf），建立三线性模型以最小化模型中的残差eijk。三维数组的元素可由负载矩阵的元素与残差计算得到，公式如下：

$$x_{ijk}=\sum _{f=1}^F{a}_{if}{b}_{jf}{c}_{kf}+e_{ijk}(1)$$

图1为公式(1)在二组分情况下的计算示意图。

该模型也可记为：$$\underline{X}=\sum _{f=1}^F{a}_f⨂b_f⨂c_f$$

其中af、bf、cf分别为矩阵A、B、C的第f列。

3.1 唯一性

PARAFAC模型的一个明显优点是解的唯一性。如果数据确实是三线性的，使用了正确数量的分量并且信噪比合适，就能得到真正的潜在光谱。

3.2 多维数组的秩（rank）

秩为1的矩阵可以写成2个向量（分数和负载向量）的外积。这样的组成部分被称为二元组。

三元组是二元组的三线性等价物，即三线性（PARAFAC）分量，是3个向量的积。

4 实现

4.1 交替最小二乘法（Alternating least squares）

PARAFAC模型的解可以通过该方法找到，方法是依次假设两种已知模式下的载荷，然后估计最后一种模式的未知参数集。这也是最初提出的对模型进行估计的方式。

PARAFAC ALS算法的流程：
（0）确定组分数F
（1）初始化B和C
（2）通过最小二乘回归，从X, B, C中估计A
（3）用同样的方法估计B
（4）用同样的方法估计C
（5）从步骤(2)开始往下执行，直到收敛。

ALS算法将在每次迭代中改善模型的拟合。如果算法收敛到全局最小值，则找到模型的最小二乘解。

ALS的优点：确保每次迭代都能优化解；ALS的主要缺点：模型估计时间长，当变量数量很多时，有时需要数百到数千次迭代才能收敛。

6 评估解

6.2 杠杆和残差

杠杆和残差可用于影响和残差分析。

6.3 组分数

提取太多的分量不仅意味着噪声被越来越多地建模，而且真实因素被更多（相关）的分量建模。

确定组分数的主要方法有三种：（1）分半实验，（2）判断残差，（3）与建模数据的外部知识进行比较。

[19]主张使用分半实验。其想法是将数据分为两半，然后在这两半上创建PARAFAC模型。通常情况下，应该以具有足够数量的自变量/样本（independent variables/samples）的模式来分割数据。

9 应用II：稀疏荧光数据的唯一分解

9.1 数据

这个问题是PARAFAC使用非负约束获得唯一分解的一个示例。

样品：含有不同量的酪氨酸、色氨酸和苯丙氨酸的2个样品。

因此，要分解的数组是2×51×201。

在多线性分解中应该避免瑞利散射，有三种方法可以做到这一点：（iii）测量空白，并从样品测量值中减去该测量值。在这个实验中，最初没有采取任何措施来消除瑞利散射。

9.2 结果与讨论

（1）估计的激发光谱如图10（上图）所示。
①三组分PARAFAC溶液的发射负载如图10a所示。从中可以看出，与色氨酸相对应的光谱具有大的负区域。得出的结论是，由于变化性小（两个样品），分解很困难。由于我们知道荧光光谱和浓度应该是正的，所以很自然地将PARAFAC 负载限制在正值 。
②在图10b中，使用非负性约束显示了估计的发射负荷。估计的光谱与分析物的纯光谱非常相似，但对于色氨酸，由于非多重线性瑞利散射，在300mn以下有一个小峰。
③为了避免这种情况，试图将受瑞利散射影响的所有变量设置为缺失值，然后估计相应的PARAFAC模型，结果如图10c所示。
④显然，仅凭这一点不足以确保色氨酸光谱具有良好的曲线分辨率。将缺失元素方法与非负约束相结合，有助于模型关注图中数据的正确方面。在图10d中，估计的发射负载（实线）与纯光谱（虚线，注意区分）一起显示。
（2）估计的激发光谱如图11（下图）所示。

下面使用N-way工具箱的示例数据claus.mat尝试复现图10。

打开Matlab，将数据导入工作区，将EEM数据X重整为5×201×61规格。

绘制5张光谱的填充等高线图：

for i = 1:DimX(1)
    subplot(2,3,i),contourf(squeeze(X(i,:,:)))
end

可见，前3张光谱为纯光谱，后2张光谱为混合光谱。
从X中提取第4-5张光谱：P = X(4:5,:,:)
将工作区保存为claus2.mat

[Factors1] = parafac(P,3);
plotfac(Factors1)

下图还原了图10(a)：

增加非负约束：

[Factors2] = parafac(P,3,0,[2,2,2]);
plotfac(Factors2)

下图还原了图10(b)：

接下来使用drEEM工具箱尝试复现图10d。

ds=assembledataset(X,ExAx,EmAx,'AU')
eemreview(ds)
ds2=smootheem(ds,[18 24],[],[],[],[0 0 0 0],[],3382,0)
eemreview(ds2)
spectralvariance(ds2)
ds3=subdataset(ds2,1:3,[],[])
[Factors] = parafac(ds3.X,3,0,[2 2 2]);
plotfac(Factors)

结果如下图所示。