1. 偏最小二乘回归

线性回归算法通常是根据最小均方误差MSE准则计算权重$w$，然后利用权重对输入数据进行线性加权，得到预测结果。

通常认为学习到的权重$w$表示特征各维度的重要性。$w$中正的权重表示该维度对预测结果起正向作用；负的权重表示该维度对预测结果起反向作用。

然而上述认知并不总是正确的，产生错误的原因可能是训练数据量太少，或特征维度之间的相关性过大。

比如对于数据$X=(x_1,x_2) = { (1,2), (2,4), (3,6), (4,8) }$，标签$y= { 1, 2, 3, 4 }$。通过线性回归得到的结果为$y=-x_1+x_2\propto -x_1$；然而直接观察可知存在关系$y \propto x_1$。产生上述矛盾的原因是数据的两个特征维度是相关的$x_2=2x_1$。

为了降低数据维度之间的相关性，同时保留数据和标签之间的相关性，引入偏最小二乘回归(Partial Least Squares, PLS)。

偏最小二乘回归的实现过程为寻找一组权重$w_1,…,w_d$，把每个数据$x^{n}$的特征维度线性组合为一个新的特征$t^{n}$:

\[t^{n} = w_1x^{n}_1+ w_2x^{n}_2+ \cdots w_dx^{n}_d\]

若共寻找$K$组权重$W \in \Bbb{R}^{d\times K}$，则把数据$X\in \Bbb{R}^{N\times d}$映射为特征维度为$K$的新数据$T=XW\in \Bbb{R}^{N\times K}$，再对数据$T$和标签$Y\in \Bbb{R}^{N\times 1}$进行线性回归：

\[Y = TC, C \in \Bbb{R}^{K\times 1}\]

寻找权重$W$时，要求构造的新数据$T$的特征维度之间尽量不相关，并且$T$与标签$Y$之间存在相关性。

2. 非线性迭代偏最小二乘回归

下面介绍求解偏最小二乘回归的非线性迭代(nonlinear iterative)算法NIPALS。

为了保证构造的新数据$T$的特征维度之间尽量不相关，循环地构造$T$的每一个特征维度（k=1:K）。额外引入一个投影矩阵$P\in \Bbb{R}^{d\times K}$，用于根据现有的新数据$T$重构原始数据$X=TP^T$。在构造完成每一个维度的新数据$T_k$后，从原始数据$(X,Y)$中减去已构造新数据的影响，从而使得新数据的下一个维度尽可能与已有维度不相关。

为了保证$T$与标签$Y$之间存在相关性，构造新数据的权重$W$通过原始数据$X$和标签$Y$的内积生成：

\[W_k = \frac{(X^{(k)})^TY^{(k)}}{||(X^{(k)})^TY^{(k)}||_2} \in \Bbb{R}^{d\times 1}\]

根据权重$W_k$构造一组新数据$T_k=X^{(k)}W_k\in \Bbb{R}^{N\times 1}$，通过线性回归求解$Y^{(k)}=T_kC_k$:

\[C_k = (T_k^TT_k)^{-1}T_k^TY^{(k)} = \frac{T_k^TY^{(k)}}{||T_k||_2^2} \in \Bbb{R}^{1}\]

如果标量$C_k$数值比较小，说明新构造的特征维度对回归的影响已经比较小，则可以退出算法。

否则把新数据$T_k$投影到原始数据空间$\hat{X}^{(k)}=T_kP_k^T$，仍然通过线性回归求解：

\[P_k = ((T_k^TT_k)^{-1}T_k^TX^{(k)})^T = \frac{(X^{(k)})^TT_k}{||T_k||_2^2} \in \Bbb{R}^{d\times 1}\]

把新构造的数据分别投影到原始数据空间和标签空间后从中减去，以消除现有特征对新特征的相关性：

\[\begin{aligned} X^{(k+1)} &= X^{(k)} - T_kP_k^T \\ Y^{(k+1)} &= Y^{(k)} - T_kC_k \end{aligned}\]

训练完成后，对于新的数据$Z\in \Bbb{R}^{M\times d}$，可以迭代地预测结果：

\[\begin{aligned} k=&1:K \\ &T_k = Z^{(k)}W_k \\ &Z^{(k+1)} = Z^{(k)} - T_kP_k^T \\ Y =& TC \end{aligned}\]

NIPALS算法的程序实现如下：

class PLS():
    def __init__(self, k=2):
        # 最终的隐变量的成分数
        self.components = k

    def fit(self, X, Y):
        # 获取特征维度
        N, D = np.shape(X)
        K = self.components
        # 对X进行降维时，K个基的系数 
        W = np.empty([D, K])
        # 利用隐变量对X进行回归的系数
        P = np.empty([D, K])
        # 存储变换后的隐变量
        T = np.empty([N, K])
        # 隐变量T对Y的回归系数
        c = np.empty([K, 1])

        X_k = X # [N, d]
        Y_k = Y # [N, 1]
        for k  in range(K):
            # 计算X每个维度上的特征与Y的相关性
            # 并利用这个相关性作为初始权重
            w_k = X_k.T @ Y_k # [d, 1]
            w_k /= np.linalg.norm(w_k, 2)
            # 对X进行加权求和得到t
            t_k = X_k @ w_k # [N, 1]
            # 利用t对Y进行回归得到系数c    
            c_k = (t_k.T @ Y_k) / (t_k.T @ t_k) # [1, ]

            if c_k < 1e-6:
                self.components = k
                break
        
            # 利用t对X进行回归得到回归系数P
            p_k = (X_k.T @ t_k) / (t_k.T @ t_k) # [d, 1]
            # 利用t,P计算X的残差
            X_k = X_k - t_k @ p_k.T # [n, d]
            # 利用t,c计算Y的残差
            Y_k = Y_k - t_k * c_k # [n, 1]

            # 中间结果存储
            W[:, k] = w_k[:, 0]
            P[:, k] = p_k[:, 0]
            T[:, k] = t_k[:, 0]
            c[k, :] = c_k
            # 利用X,Y残差进行下一轮的迭代
        
        self.W = W[:, 0:self.components]
        self.P = P[:, 0:self.components]
        self.T = T[:, 0:self.components]
        self.c = c[0:self.components]

    def predict(self, Z):
        N, _ = np.shape(Z)
        t = np.empty((N, self.components))
        Z_k = Z

        for k in range(self.components):
            w_k = np.expand_dims(self.W[:, k], axis=-1) # [d, 1]
            p_k = np.expand_dims(self.P[:, k], axis=-1) # [d, 1]
            t_k = Z_k @ w_k
            Z_k = Z_k - t_k @ p_k.T
            t[:, k]=t_k[:, 0]

        result = t @ self.c
        return result


pls = PLS(2)
X = np.array(
    [[1,2],[2,4],[3,6],[4,8]],
    dtype=np.float64)
Y = np.array(
    [[1],[2],[3],[4]],
    dtype=np.float64)
pls.fit(X, Y)
print(pls.predict(X))