GZ Cloudhome Logo

偏最小二乘：基于协方差的有监督隐变量模型

发布于：2023 年 6 月 9 日

长久以来，我对偏最小二乘（Partial Least Squares，PLS） 整个模型都是处于一个“知道一点概念，但是细说说不上来”的状态。因此，最近花了几天的时间研究了一下 PLS 模型，这篇也算是对最近的一些研究的总结。

引言

PLS 尝试对两个矩阵 X∈RN×p\mathbf{X} \in \mathbb{R}^{N \times p}，Y∈RN×q\mathbf{Y} \in \mathbb{R}^{N \times q} 之间的关系进行建模。其中 NN 是样本数，pp 和 qq 分别是 X\mathbf{X} 中的变量个数和 Y\mathbf{Y} 中的变量个数。PLS 会尝试寻找 X\mathbf{X} 各列的线性组合 ξ\boldsymbol{\xi}，和 Y\mathbf{Y} 的各列的线性组合 ω\boldsymbol{\omega}，使得 ξ\boldsymbol{\xi} 能够最大程度上解释 ω\boldsymbol{\omega} 中的方差（Variance）。换言之，我们期望 ξ\boldsymbol{\xi} 和 ω\boldsymbol{\omega} 之间的协方差（Covariance）最大。用一个最优化问题来表述，即为

max⁡u1,v1ξ1Tω1=u1TXTYv1s.t.∥u1∥=∥v1∥=1\begin{align*} \underset{\boldsymbol{u}_1, \boldsymbol{v}_1}{\max} && & \boldsymbol{\xi}_1^T \boldsymbol{\omega}_1 = \boldsymbol{u}_1^T \mathbf{X}^T \mathbf{Y} \boldsymbol{v}_1\\ \mathrm{s.t.} && & \| \boldsymbol{u}_1 \| = \|\boldsymbol{v}_1\| = 1 \end{align*}

这里 u1∈Rp\boldsymbol{u}_1 \in \mathbb{R}^p 和 v1∈Rq\boldsymbol{v}_1 \in \mathbb{R}^q 即为 X\mathbf{X} 和 Y\mathbf{Y} 中各列的线性组合系数。因为 u1TXTYv1\boldsymbol{u}_1^T \mathbf{X}^T \mathbf{Y} \boldsymbol{v}_1 的结果与 u1\boldsymbol{u}_1 和 v1\boldsymbol{v}_1 的模有关，我们不希望系数向量的模影响到协方差的计算，所以有限定条件 ∥u1∥=∥v1∥=1\| \boldsymbol{u}_1 \| = \|\boldsymbol{v}_1\| = 1，即限定它们的模长均为 1。

基于拉格朗日乘子法，得到增广拉格朗日函数：

L=u1TXTYv1−12λuu1Tu1−12λvv1Tv1.L = \boldsymbol{u}_1^T \mathbf{X}^T \mathbf{Y} \boldsymbol{v}_1 - \frac{1}{2} \lambda_u \boldsymbol{u}_1^T \boldsymbol{u}_1 - \frac{1}{2} \lambda_v \boldsymbol{v}_1^T \boldsymbol{v}_1.

分别基于 LL 对 u1\boldsymbol{u}_1 和 v1\boldsymbol{v}_1 求偏导，令之为 0 可得

∂L∂u1=XTYv1−λuu1=0∂L∂v1=YTXu1−λvv1=0.\begin{align*} \frac{\partial L}{\partial \boldsymbol{u}_1} &= \mathbf{X}^T \mathbf{Y} \boldsymbol{v}_1 - \lambda_u \boldsymbol{u}_1 = 0 \\ \frac{\partial L}{\partial \boldsymbol{v}_1} &= \mathbf{Y}^T \mathbf{X} \boldsymbol{u}_1 - \lambda_v \boldsymbol{v}_1 = 0. \end{align*}

上式第一行左乘 u1T\boldsymbol{u}_1^T，第二行左乘 v1T\boldsymbol{v}_1^T，可以得到

u1TXTYv1=λu=λv.\boldsymbol{u}_1^T \mathbf{X}^T \mathbf{Y} \boldsymbol{v}_1 = \lambda_u = \lambda_v.

因此，记 λ=λu=λv=u1TXTYv1\lambda = \lambda_u = \lambda_v = \boldsymbol{u}_1^T \mathbf{X}^T \mathbf{Y} \boldsymbol{v}_1 （注意 λ\lambda 即为 ξ1\boldsymbol{\xi}_1 和 ω1\boldsymbol{\omega}_1 的协方差），上式可以化简为

XTYv1=λu1(XTY)Tu1=λv1\begin{align*} \mathbf{X}^T \mathbf{Y} \boldsymbol{v}_1 &= \lambda \boldsymbol{u}_1\\ (\mathbf{X}^T \mathbf{Y})^T \boldsymbol{u}_1 &= \lambda \boldsymbol{v }_1 \end{align*}

容易看出，u1\boldsymbol{u}_1 和 v1\boldsymbol{v}_1 分别为 XTY\mathbf{X}^T \mathbf{Y} 对应奇异值 λ\lambda 的左奇异向量和右奇异向量。而我们期望协方差 λ\lambda 最大化，所以令 u1\boldsymbol{u}_1 和 v1\boldsymbol{v}_1 对应 XTY\mathbf{X}^T \mathbf{Y} 的最大奇异值，即可最大化协方差。

我们往往不期望就此止步。前面我们只找到了协方差最大的一组组合 ξ1\boldsymbol{\xi}_1 和 ω1\boldsymbol{\omega}_1；那么是否还有其它组合呢？答案是肯定的，因为显然 X\mathbf{X} 和 Y\mathbf{Y} 中的协方差并没有挖掘完毕。如何进一步寻找其它的组合？这就涉及了 PLS 的几种不同的变种：

三种不同类型的 PLS

PLS-W2A

PLS-W2A 的思路是，ξ1,ξ2,…\boldsymbol{\xi}_1, \boldsymbol{\xi}_2, \dots 之间，以及 ω1,ω2,…\boldsymbol{\omega}_1, \boldsymbol{\omega}_2, \dots 之间，都需要线性无关。为了达到这一点，我们需要从 X\mathbf{X} 的列空间和 Y\mathbf{Y} 的列空间中分别除去和 ξ1\boldsymbol{\xi}_1 和 ω1\boldsymbol{\omega}_1 相关的部分，然后再对剩下的部分做奇异值分解。这可以用基于最小二乘的线性回归来解决。

基于上面的思路，当我们提取第 rr 组隐变量 ξr\boldsymbol{\xi}_r 和 ωr\boldsymbol{\omega}_r 时，

X(r)=X(r−1)−ξr−1γr−1TY(r)=Y(r−1)−ωr−1δr−1T,\begin{align*} \mathbf{X}^{(r)} &= \mathbf{X}^{(r-1)} - \boldsymbol{\xi}_{r-1} \boldsymbol{\gamma}_{r-1}^T \\ \mathbf{Y}^{(r)} &= \mathbf{Y}^{(r-1)} - \boldsymbol{\omega}_{r-1} \boldsymbol{\delta}_{r-1}^T, \end{align*}

注意， X(1)=X\mathbf{X}^{(1)} = \mathbf{X}，Y(1)=Y\mathbf{Y}^{(1)} = \mathbf{Y}，且

γr=ξrTXξrTξδr=ωrTYωrTωr,\begin{align*} \boldsymbol{\gamma}_r &= \frac{\boldsymbol{\xi}_r^T \mathbf{X}}{\boldsymbol{\xi}_r^T \boldsymbol{\xi}} \\ \boldsymbol{\delta}_r &= \frac{\boldsymbol{\omega}_r^T \mathbf{Y}}{\boldsymbol{\omega}_r^T \boldsymbol{\omega}_r}, \end{align*}

这里面 γr\boldsymbol{\gamma}_r 和 δr\boldsymbol{\delta}_r 分别是用 ξr\boldsymbol{\xi}_r 回归 X\mathbf{X} 的线性回归系数和用 ωr\boldsymbol{\omega}_r 回归 Y\mathbf{Y} 的线性回归系数。这样的操作又叫做一阶 deflation。

这样，Xr\mathbf{X}^{r} 的各列将垂直于 ξ1,…,ξr−1\boldsymbol{\xi}_1, \dots, \boldsymbol{\xi}_{r-1}，Y(r)\mathbf{Y}^{(r)} 也同理。基于 X(r)\mathbf{X}^{(r)} 和 Y(r)\mathbf{Y}^{(r)} 采取同样的奇异值分解策略，就可以计算得到 ξr\boldsymbol{\xi}_r 和 ωr\boldsymbol{\omega}_r。

PLS-W2A 相当于在引言中的优化目标的基础上，加入了额外的约束。新的优化目标为：

max⁡ur,vrξrTωr=urT(X(r))TY(r)vrs.t.∥ur∥=∥vr∥=1ξr⊥ξlωr⊥ωl,      1≤l<r.\begin{align*} \underset{\boldsymbol{u}_r, \boldsymbol{v}_r}{\max} && & \boldsymbol{\xi}_r^T \boldsymbol{\omega}_r = \boldsymbol{u}_r^T (\mathbf{X}^{(r)})^T \mathbf{Y}^{(r)} \boldsymbol{v}_r \\ \mathrm{s.t.} && & \| \boldsymbol{u}_r \| = \|\boldsymbol{v}_r\| = 1 \\ && & \boldsymbol{\xi}_r \perp \boldsymbol{\xi}_l \\ && & \boldsymbol{\omega}_r \perp \boldsymbol{\omega}_l, ~~~~~~ 1 \le l < r. \end{align*}

或者等价地，

max⁡pr,qrξrTωr=prTXTYqrs.t.∥pr∥=∥qr∥=1prTXTXpl=0qrTYTYql=0,      1≤l<r.\begin{align*} \underset{\boldsymbol{p}_r, \boldsymbol{q}_r}{\max} && & \boldsymbol{\xi}_r^T \boldsymbol{\omega}_r = \boldsymbol{p}_r^T \mathbf{X}^T \mathbf{Y} \boldsymbol{q}_r \\ \mathrm{s.t.} && & \| \boldsymbol{p}_r \| = \|\boldsymbol{q}_r\| = 1 \\ && & \boldsymbol{p}_r^T \mathbf{X}^T \mathbf{X} \boldsymbol{p}_l = 0 \\ && & \boldsymbol{q}_r^T \mathbf{Y}^T \mathbf{Y} \boldsymbol{q}_l = 0, ~~~~~~ 1 \le l < r. \end{align*}

注意：pr,qr\boldsymbol{p}_r, \boldsymbol{q}_r 和 ur,vr\boldsymbol{u}_r, \boldsymbol{v}_r 并不相同，但是具有一定的转化关系：

P=U(ΓTU)−1,\mathbf{P} = \mathbf{U}\big( \mathbf{\Gamma}^T \mathbf{U} \big)^{-1},

其中 P\mathbf{P} 各列正比于 pr\boldsymbol{p}_r，Γ\mathbf{\Gamma} 各列为 γr\boldsymbol{\gamma}_r，同样地也可以得出

Q=V(ΔTV)−1.\mathbf{Q} = \mathbf{V} \big( \mathbf{\Delta}^T \mathbf{V} \big)^{-1}.

PLS-SVD

PLS-SVD 的思路是逐步拆解 X\mathbf{X} 和 Y\mathbf{Y} 的协方差矩阵 XTY\mathbf{X}^T \mathbf{Y}。设 XTY\mathbf{X}^T \mathbf{Y} 的奇异值分解为（不失一般性，假设 q<pq < p）

XTY=ADBT=[α1…αp][d1⋱dq][β1T⋮βqT]=∑r=1qαrdrβrT,\begin{align*} \mathbf{X}^T \mathbf{Y} &= \mathbf{A} \mathbf{D} \mathbf{B}^T \\ &= \left[ \begin{matrix} \boldsymbol{\alpha}_1 & \dots & \boldsymbol{\alpha}_p \end{matrix} \right] \left[ \begin{matrix} d_1 & & \\ & \ddots & \\ & & d_q \\ && \end{matrix} \right] \left[ \begin{matrix} \boldsymbol{\beta}_1^T \\ \vdots \\ \boldsymbol{\beta}_q^T \end{matrix} \right] \\ &= \sum_{r=1}^q \boldsymbol{\alpha}_r d_r \boldsymbol{\beta}_r^T, \end{align*}

那么我们一共可以提取出 qq 组 ξr\boldsymbol{\xi}_r 和 ωr\boldsymbol{\omega}_r（1≤r≤q1 \le r \le q）：

ξr=Xαrωr=Yβr,\begin{align*} \boldsymbol{\xi}_r &= \mathbf{X} \boldsymbol{\alpha}_r \\ \boldsymbol{\omega}_r &= \mathbf{Y} \boldsymbol{\beta}_r, \end{align*}

注意它和 PLS-W2A 的不同之处：PLS-SVD 得到的 ξ1,ξ2,…\boldsymbol{\xi}_1, \boldsymbol{\xi}_2, \dots 之间，以及 ω1,ω2,…\boldsymbol{\omega}_1, \boldsymbol{\omega}_2, \dots 之间，不一定是相互垂直的。

PLS2 和 PLS1

和 PLS-W2A 相比，PLS2 中，在得到 ξr\boldsymbol{\xi}_r 和 ur\boldsymbol{u}_r 之后，仍对 X(r)\mathbf{X}^{(r)} 做相同的一阶 deflation，但是 Y(r)\mathbf{Y}^{(r)} 则根据 ξ1\boldsymbol{\xi}_1 来做 deflation。这是因为 PLS2 一般看作一个回归模型，我们需要逐步基于 X\mathbf{X} 中的隐变量来预测 Y\mathbf{Y}，因此每步只需要计算 ξr\boldsymbol{\xi}_r 即可；甚至，我们不需要对 Y(r)\mathbf{Y}^{(r)} 做 deflation，因为不管做不做 deflation，最终得到的结果都一样。

PLS1 则是 PLS2 的一种特殊情况——当 Y\mathbf{Y} 只有一列时，PLS2 退化为 PLS1。

NIPALS 算法

前面介绍的 PLS-W2A、PLS2 和 PLS1 都是基于奇异值分解来计算第 rr 对的隐变量的；事实上，它们也可以通过迭代的方法求得，对应的迭代法即为幂法。事实上，这也是 NIPALS 算法中用到的。