从拉格朗日乘数到LASSO

拉格朗日乘数法（Lagrange multiplier，以数学家约瑟夫·拉格朗日命名），在数学中的最优化问题中，是一种寻找多元函数在其变量受到一个或多个条件的约束时的局部极值的方法。

1. 最优化问题

最优化，是应用数学的一个分支。主要研究在特定情况下最大化或最小化某一特定函数或变量。对于没有条件约束的一元函数可以求出所有导数为零的点，代入原函数得出目标函数的极值，对于多元函数可以使用梯度下降法求得。而对于有条件约束的目标函数，上面两种方法似乎就不太可行了。对于这类问题主要有两种解体思路。

1.1 转化为无条件极值问题

例：求在条件下的最小值。

解：由易得, 带入到目标函数可得：
.
显然当时, 目标函数有最小值

这类方法显示不适用于约束条件复杂，或者目标函数有很多自变量的问题。

1.2 拉格朗日乘数法

拉格朗日乘数法可以将一个有个变量与个约束条件的最优化问题转换为一个解有个变量的方程组的解的问题。这种方法中引入了一个或一组新的未知数，即拉格朗日乘数，又称拉格朗日乘子，或拉氏乘子，它们是在转换后的方程，即约束方程中作为梯度的线性组合中各个向量的系数。

比如，要求, 在时的局部极值时，我们可以引入新变量拉格朗日乘数 ，这时我们只需要求下列拉格朗日函数的局部极值：

更一般地，对含个变量和个约束的情况，有：

下面是对于和的推导过程：

设是函数在条件下的极小值点.

另外, 由, 得出关于变量的隐函数, 即是的极值点, 可以得到:

$$\begin{align*} & \frac{\textrm{dz}}{\textrm{d}x}|_{x=x_0} = f_x(x_0,y_0) + f_y(x_0, y_0) \cdot \frac{\textrm{d}y}{\textrm{d}x}|_{x=x_0} = 0 \\ & \because \frac{\textrm{d}y}{\textrm{d}x} = -\frac{\varphi_x}{\varphi_y} \\ & \therefore \frac{\textrm{dz}}{\textrm{d}x}|_{x=x_0} = f_x(x_0,y_0) - f_y(x_0, y_0) \cdot \frac{\varphi_x(x_0, y_0)}{\varphi_y(x_0, y_0)}|_{x=x_0} = 0 \\ & \frac{\textrm{dz}}{\textrm{d}x}|_{x=x_0} = f_x(x_0,y_0) - \varphi_x(x_0, y_0) \cdot \frac{f_y(x_0, y_0)}{\varphi_y(x_0, y_0)}|_{x=x_0} = 0 \\ & \textrm{make.}\, \varphi_x(x_0, y_0) \cdot \frac{f_y(x_0, y_0)}{\varphi_y(x_0, y_0)}|_{x=x_0} = - \lambda \\ & \left\{ \begin{array}{**lr**} f_x + \lambda \varphi_x = 0 \\ f_y+\lambda \varphi_y=0 \\ \varphi = 0 \end{array} \right. \end{align*}$$

上述方程组中的被称为”拉格朗日乘子“，满足此方程组的点即可能为目标函数在约束条件下的一个极值点。为此可构建函数:

$$\mathcal {L}(x,y,\lambda)=f(x,y)+\lambda \varphi(x,y)$$

满足的点为潜在的极值点。

e.g.

$$\begin{align*} & \textrm{make}. \ \mathcal{L}(x,y,\lambda) = xy - \lambda(x^2y^2 - x^4 - y^2) = 0 \\ & \left\{ \begin{array}{**lr**} \mathcal{L}_x=y-\lambda(2xy^2 - 4x^3)=0 &\ \ \ (1) \\ \mathcal{L}_y = x-\lambda(2yx^2-2y)=0 &\ \ \ (2) \\ \mathcal{L}_{\lambda} = x^2y^2 - x^4 - y^2 =0 &\ \ \ (3) \\ \end{array} \right. \\ & (1) \times (2yx^2-2y) - (2) \times (2xy^2 - 4x^3) = y(2yx^2-2y) - x(2xy^2 - 4x^3) = 0 \\ & \Rightarrow 2x^4= y^2 \\ & (3) = 2x^6 - x^4 - 2x^4 = 0 \end{align*}$$

之后求上式的导数，后面过程省略。

2. 线性回归模型

2.1 多元线性回归

线性模型的基本结构为：

$$f(x)=w_1x_1 + w_2x_2 +...+w_dx_d + b \tag{2.1}$$

写成向量形式为：

$$f(x)=\boldsymbol{w x} + b \tag{2.2}$$

其中。定义损失函数为：

$$\sum\limits_{i=1}^m(y_i-\hat{y_i})^2 = \sum\limits_{i=1}^m(y_i - \boldsymbol{x}_i\boldsymbol{w})^2 \tag{2.3}$$

其实就是范数的平方结果，损失函数的求解目标可以转化为：

$$\min\limits_w{\Vert \boldsymbol{y} - \bf{X} \boldsymbol{w} \Vert}_2^2 \tag{2.4}$$

对求偏导：

$$\begin{align*} \frac {\partial{\Vert \boldsymbol{y} - \bf{X} \boldsymbol{w} \Vert}_2^2} {\partial \boldsymbol{w}} & = \frac {\partial{(\boldsymbol{y} - \bf{X}\boldsymbol{w})}^\mathrm{T}(\boldsymbol{y} - \bf{X}\boldsymbol{w})} {\partial \boldsymbol{w}} \\ & = \frac {\partial(\boldsymbol{y}^\mathrm{T}\boldsymbol{y} - \boldsymbol{w}^\mathrm{T}\bf{X}^\mathrm{T}\boldsymbol{y} - \boldsymbol{y}^\mathrm{T}\bf{X}\boldsymbol{w}+\boldsymbol{w}^\mathrm{T}\bf{X}^\mathrm{T}\bf{X}\boldsymbol{w})} {\partial \boldsymbol{w}} \\ & = 0 - \bf{X}^\mathrm{T}\boldsymbol{y} - \bf{X}^\mathrm{T}\boldsymbol{y} + 2\bf{X}^\mathrm{T}\bf{X}\boldsymbol{w} \\ & = 2 \bf{X}^\mathrm{T}\bf{X}\boldsymbol{w} - 2\bf{X}^\mathrm{T}\boldsymbol{y} \end{align*} \tag{2.5}$$

令导数结果为0，得到：

$$\bf{X}^\mathrm{T}\bf{X}\boldsymbol{w} = \bf{X}^\mathrm{T}\boldsymbol{y} \tag{2.6}$$

左乘一个：

$$\boldsymbol{w} = (\bf{X}^\mathrm{T}\bf{X})^{-1}\bf{X}^\mathrm{T}\boldsymbol{y} \tag{2.7}$$

令，将代入，得：

$$f(\hat{x}_i) = \hat{x}_i (\bf{X}^\mathrm{T}\bf{X})^{-1}\bf{X}^\mathrm{T}y \tag{2.8}$$

import numpy as np
from numpy.typing import ArrayLike


class MLR:

    def __init__(self):
        self._w = None

    @property
    def w(self):
        return self._w

    def fit(self, X: ArrayLike, y: ArrayLike):
        size = X.shape[0]
        ones = np.ones((size, 1))
        X_ones = np.concatenate((ones, X), axis=1)

        X_mat = np.mat(X_ones)
        y_mat = np.mat(y.reshape(-1, 1))

        temp = X_mat.T * X_mat

        self._w = temp.I * X_mat.T * y_mat

    def predict(self, X_test: ArrayLike):
        assert self._w is not None

        _X_test = np.array(X_test)
        m, n = X_test.shape

        assert n == (self.w.shape[0] - 1)

        ones = np.ones((m, 1))
        X_mat = np.mat(np.concatenate((ones, X_test), axis=1))

        yi = X_mat * self.w
        
        return np.array(yi).flatten()

import matplotlib.pyplot as plt
from sklearn.datasets import make_regression

_X, _y = make_regression(n_features=1, noise=5, random_state=777)
xMin = _X.min()
xMax = _X.max()

_noise_X = np.linspace(1, 2, 5).reshape(-1, 1)
_noise_y = np.array([350, 380, 410, 430, 480])
_Xn = np.r_[_X, _noise_X]
_yn = np.r_[_y, _noise_y]

f, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 9))

base_model = MLR()
base_model.fit(_X, _y)
base_pred = base_model.predict(_X)

base_model_ = MLR()
base_model_.fit(_Xn, _yn)
base_pred_ = base_model_.predict(_Xn)

ax1.scatter(_X, _y.reshape(-1, 1), color='y', label="data without outlier")
ax2.scatter(_Xn, _yn.reshape(-1, 1), color='y', label="data with outlier")

ax1.plot(_X, base_pred, lw=2, label='MLR without outlier')
ax2.plot(_Xn, base_pred_, lw=2, label='MLR with outlier')

ax1.legend()
ax2.legend()
plt.show()
    

在上面的例子中，除了随机生成了100个线性相关的点，还增加了5个离群值，并分别训练了2个模型。可以看到，使用离群值训练的模型受到了较大的影响。

2.2 岭回归

然而，现实任务中往往不是满秩矩阵。例如在许多任务中我们会遇到大量的变量，其数目甚至超过样本例数，导致的列数多于行数，显然不满秩。此时可解出多个，它们都能使均方差误差最小化。选择哪一个解作为输出，将由学习算法的归纳偏好决定，常见的做法是引入正则化（regularization）项。若使用，则有：

$$\min\limits_w{ \Vert \boldsymbol{y} - \boldsymbol{X} \boldsymbol{w} \Vert_2^2} + \lambda \Vert \boldsymbol{w} \Vert_2^2 \tag{2.9}$$

式被称为“岭回归”（ridge regression），其中正则化参数。相较于式，多了正则化项 。式同样对求偏导：

$$\begin{align*} \frac {\partial ({\Vert \boldsymbol{y} - \bf{X} \boldsymbol{w} \Vert}_2^2 + \lambda \Vert \boldsymbol{w} \Vert_2^2 )} {\partial \boldsymbol{w}} & = \frac {\partial{(\boldsymbol{y} - \bf{X}\boldsymbol{w})}^\mathrm{T}(\boldsymbol{y} - \bf{X}\boldsymbol{w})} {\partial \boldsymbol{w}} + \frac {\partial \lambda \Vert \boldsymbol{w} \Vert_2^2} {\partial \boldsymbol{w}} \\ & = 2 \bf{X}^\mathrm{T}\bf{X}\boldsymbol{w} - 2 \bf{X}^\mathrm{T}\boldsymbol{y} + 2 \lambda \boldsymbol{w} \\ & = 2(\bf{X}^\mathrm{T}\bf{X} + \lambda \bf{E}) \boldsymbol{w} - 2 \bf{X}^\mathrm{T}\boldsymbol{y} \end{align*} \tag{2.10}$$

令式等于0：

$$(\bf{X}^\mathrm{T}\bf{X} + \lambda \bf{E}) \boldsymbol{w} = \bf{X}^\mathrm{T}\boldsymbol{y} \tag{2.11}$$

式左乘:

$$\boldsymbol{w} = (\bf{X}^\mathrm{T}\bf{X} + \lambda \bf{E}) ^{-1} \bf{X}^\mathrm{T}\boldsymbol{y} \tag{2.12}$$

同样，令，将代入，得：

$$f(\hat{x}_i) = \hat{x}_i (\bf{X}^\mathrm{T}\bf{X} + \lambda \bf{E})^{-1}\bf{X}^\mathrm{T}y \tag{2.13}$$

多元线性回归中的未必是满秩的，由于在主对角线上加上了，在多数情况下此矩阵会是满秩的，在很大程度上可以减轻（消除）数据中多重共线性对偏最小二乘的影响，但会以牺牲模型性能为代价。

import numpy as np
from numpy.typing import ArrayLike


class RidgeRegr:

    def __init__(self, alpha: float = 0):  # lambda is keyword
        self._w = None
        self._alpha = alpha

    @property
    def w(self):
        return self._w

    @property
    def alpha(self):
        return self._alpha

    def fit(self, X: ArrayLike, y: ArrayLike):
        X = np.array(X)
        size = X.shape[0]

        ones = np.ones((size, 1))
        X_hat = np.concatenate((ones, X), axis=1)

        X_mat = np.mat(X_hat)

        temp = X_mat.T * X_mat
        eye = np.mat(np.eye(temp.shape[1]))

        y_mat = np.mat(y.reshape(-1, 1))

        self._w = (temp + self._alpha * eye).I * X_mat.T * y_mat

    def predict(self, X_test: ArrayLike):
        assert self._w is not None

        _X_test = np.array(X_test)
        m, n = X_test.shape

        assert n == (self.w.shape[0] - 1)

        ones = np.ones((m, 1))
        X_mat = np.mat(np.concatenate((ones, X_test), axis=1))

        yi = X_mat * self.w

        return np.array(yi).flatten()

import matplotlib.pyplot as plt
from sklearn.datasets import make_regression

_X, _y = make_regression(n_features=1, noise=5, random_state=777)
xMin = _X.min()
xMax = _X.max()

_noise_X = np.linspace(1, 2, 5).reshape(-1, 1)
_noise_y = np.array([350, 380, 410, 430, 480])
_Xn = np.r_[_X, _noise_X]
_yn = np.r_[_y, _noise_y]

f, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 9))

ax1.scatter(_X, _y.reshape(-1, 1), color='y', label="data without outlier")
ax2.scatter(_Xn, _yn.reshape(-1, 1), color='c', label="data with outlier")

for _alpha in [0, 10, 100, 200, 1000]:
    ridge1 = RidgeRegr(alpha=_alpha)
    ridge1.fit(_X, _y)
    y_pred1 = ridge1.predict(_X)
    ax1.plot(_X, y_pred1, lw=2, label=fr'$Ridge \ with \ outlier (\lambda = {_alpha})$')

    ridge2 = RidgeRegr(alpha=_alpha)
    ridge2.fit(_Xn, _yn)
    y_pred2 = ridge2.predict(_Xn)
    ax2.plot(_Xn, y_pred2, lw=2, label=fr'$Ridge \ with \ outlier (\lambda = {_alpha})$')

    ax1.legend()
    ax2.legend()
    plt.show()
    

可以看到，随着正则化参数的增加，虽然在没有添加离群值的数据集中模型的欠拟合现象愈发严重。但在真实的数据集中，离群值是很难避免的，岭回归由于添加了正则化项，模型泛化能力得到提升。

2.3 LASSO

和岭回归一样，LASSO（Least Absolute Shrinkage and Selection Operator）也是被 创作用于多重共线性问题的算法，不过LASSO使用的是正则化，其损失函数为：

$$\min\limits_w{ \Vert \boldsymbol{y} - \boldsymbol{X} \boldsymbol{w} \Vert_2^2} + \lambda \Vert \boldsymbol{w} \Vert_1 \tag{2.14}$$

对式中的求偏导：

$$\begin{align*} \frac {\partial ({\Vert \boldsymbol{y} - \bf{X} \boldsymbol{w} \Vert}_2^2 + \lambda \Vert \boldsymbol{w} \Vert_1 )} {\partial \boldsymbol{w}} & = \frac {\partial{(\boldsymbol{y} - \bf{X}\boldsymbol{w})}^\mathrm{T}(\boldsymbol{y} - \bf{X}\boldsymbol{w})} {\partial \boldsymbol{w}} + \frac {\partial \lambda \Vert \boldsymbol{w} \Vert_1} {\partial \boldsymbol{w}} \\ & = 2 \bf{X}^\mathrm{T}\bf{X}\boldsymbol{w} - 2 \bf{X}^\mathrm{T}\boldsymbol{y} + \lambda \\ \end{align*} \tag{2.15}$$

令等于零，左乘，化简：

$$\boldsymbol{w} = (\bf{X}^\mathrm{T}\bf{X})^{-1} (\bf{X}^\mathrm{T}\boldsymbol{y} - \frac{\lambda\bf{E}}{2}) \tag{2.16}$$

通过增大，可以为计算增加一个负向，从而限制参数估计中的大小，而防止多重共线性问题引起的参数被估计过大而导致模型失准的问题。另外，范数和范数正则化都有助于降低过拟合风险，但前者还会带来一个额外的好处：它比后者更易于获得”稀疏“（sparse）解，即它求得的会有更少的非零分量，因此经常被用于特征选择。在scikit-learn中求解LASSO使用的方法为坐标下降法：

$$h_{\theta}(x)=\sum\limits_{i=1}^n \theta_ix_i=\bf{\theta}^\mathrm{T}\boldsymbol{x} \tag{2.17}$$ $$J(\theta) = \frac{1}{2m} \sum\limits_{i=1}^{m} \left( h_{\theta} (x^{(i)}) -y^{(i)} \right)^2 +\lambda \vert\theta\vert_1 \tag{2.18}$$

将损失函数调整为：

$$J(\theta) = \sum\limits_{i=1}^{m} \left( y^{(i)} - \sum\limits_{j=0}^n \theta_jx^{(i)}_j \right)^2 + \lambda \sum\limits_{j=0}^n\vert \theta \vert \tag{2.19}$$

其中，为第个样本的第个特征，由于式存在绝对值项，采取残差求导的方式，将前半部分的函数记为，其对参数的偏导数为：

$$\begin{align*} \frac {\partial} {\partial \theta_k} RSS(\theta) & = 2\sum\limits_{i=1}^{m} ( y^{(i)} - \sum\limits_{j=0}^n\theta_jx^{(i)}_j) \cdot \frac{\partial}{\partial\theta_k}\sum\limits_{k=0}^n \theta_kx^{(i)}_k \\ & = 2 \sum\limits_{i=1}^{m} ( y^{(i)} - \sum\limits_{j=0}^n\theta_jx^{(i)}_j) \cdot -x^i_k \\ & = -2 \sum\limits_{i=1}^{m} \left(x^{(i)}_k y^{(i)} - x^{(i)}_k \sum\limits_{j=0}^n\theta_jx^{(i)}_j\right) \\ & = -2 \sum\limits_{i=1}^{m} \left(x^{(i)}_k y^{(i)} - x^{(i)}_k (\sum\limits_{j=0,j \neq k}^n\theta_jx^{(i)}_j + \theta_kx^{(i)}_k)\right) \\ & = -2 \sum\limits_{i=1}^{m} \left(x^{(i)}_k y^{(i)} - x^{(i)}_k\sum\limits_{j=0,j \neq k}^n\theta_jx^{(i)}_j - \theta_k{x^{(i)}_k}^2\right) \\ & = -2 \sum\limits_{i=1}^{m} \left(x^{(i)}_k (y^{(i)} - \sum\limits_{j=0,j \neq k}^n\theta_jx^{(i)}_j)\right) + 2\theta_k \sum\limits_{i=1}^m {x^{(i)}_k}^2 \end{align*} \tag{2.20}$$

令：

$$\begin{align*} & p_k = \sum\limits_{i=1}^{m} (x^{(i)}_k (y^{(i)} - \sum\limits_{j=0,j \neq k}^n\theta_jx^{(i)}_j)) \\ & z_k = \sum\limits_{i=1}^m {x^{(i)}_k}^2 \end{align*} \tag{2.21}$$

因此有：

$$\frac{\partial RSS(\theta)}{\partial \theta_k} = -2p_k + 2 \theta_k z_k \tag{2.22}$$

对式中的正则项求偏导有：

$$\begin{align*} \frac{\partial R (\theta)}{\partial \theta_k} = \left\{ \begin{array}{**lr**} -\lambda &\ \ \ \theta_k <0 \\ [-\lambda, \lambda] &\ \ {\theta_k = 0} \\ \lambda &\ \ \ \theta_k > 0 \end{array} \right. \\ \end{align*} \tag{2.23}$$

将式, 代入到中：

$$\begin{align*} \frac{\partial J(\theta)}{\partial \theta_k} = \left\{ \begin{array}{**lr**} -2p_k + 2 \theta_k z_k -\lambda &\ \ \ \theta_k <0 \\ [-2p_k -\lambda, -2p_k +\lambda] &\ \ {\theta_k = 0} \\ -2p_k + 2 \theta_k z_k + \lambda &\ \ \ \theta_k > 0 \end{array} \right. \\ \end{align*} \tag{2.24}$$

令式为0：

$$\begin{align*} \theta_k = \left\{ \begin{array}{**lr**} (p_k + \lambda/2)/z_k &\ \ \ p_k < -\lambda/2 \\ 0 &\ \ -\lambda/2 \leq \ p_k \leq \lambda/2 \\ (p_k - \lambda/2)/z_k &\ \ \ p_k > \lambda/2 \end{array} \right. \\ \end{align*} \tag{2.25}$$

import numpy as np
from numpy.typing import ArrayLike


class LASSORegr:

    def __init__(self, alpha: float = 0):
        self._alpha = alpha
        self._w = None

    def __str__(self):
        return f"LassoRegr(alpha={self._alpha})"

    @property
    def alpha(self):
        return self._alpha

    @property
    def w(self):
        return self._w

    def fit(self, X_train: ArrayLike, y_train: ArrayLike, epochs: int = 100):

        m, n = X_train.shape  # m is the size of sample, n is the number of feat
        ones = np.ones((m, 1))
        X_hat = np.concatenate((ones, X_train), axis=1)
        X_mat = np.mat(X_hat)

        y_mat = np.mat(y_train.reshape(-1, 1))

        w = np.matrix(np.zeros((n + 1, 1)))

        for _ in range(epochs):

            for k in range(n + 1):

                _w = w.copy()
                _w[k] = 0

                pk = (X_mat[:, k].T * (y_mat - X_mat * _w))[0, 0]
                zk = (X_mat[:, k].T * X_mat[:, k])[0, 0]

                if pk < -self.alpha / 2:
                    wk = (pk + self.alpha / 2) / zk
                elif pk > self.alpha / 2:
                    wk = (pk - self.alpha / 2) / zk
                else:
                    wk = 0

                w[k, 0] = wk

        self._w = w

        print(model)

    def predict(self, X_test: ArrayLike):

        _X_test = np.array(X_test)
        m, n = X_test.shape

        assert n == (self.w.shape[0] - 1)

        ones = np.ones((m, 1))
        X_mat = np.mat(np.concatenate((ones, X_test), axis=1))

        yi = X_mat * self.w
        return np.array(yi).flatten()
    

3. 小结

岭回归和LASSO并不属于严格意义上的“算法”，而是对多元线性回归中潜在的多重共线性问题的一种优化方法。在原有损失函数的基础上，通过加入正则化项，限制模型参数的范围。

4. 参考资料

1. 机器学习 - 周志华
2. 拉格朗日乘数法的推导
3. Ridge回归求解
4. LASSO回归求解