牛顿迭代法及其应用

牛顿迭代法及其应用

考虑无约束优化问题
$$\underset{\mathbf{x}}{\min }f(\mathbf{x})$$
假设$f(\mathbf{x})$二阶可微，在第$k$次的迭代值为${\mathbf{x}}^k$，可将$f(\mathbf{x})$在${\mathbf{x}}^k$处泰勒展开
$$f(\mathbf{x})=f({\mathbf{x}}^k)+(\mathbf{x}-{\mathbf{x}}^k{)}^T\nabla f({\mathbf{x}}^k)+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}^k{)}^{T}H({\mathbf{x}}^k)(\mathbf{x}-{\mathbf{x}}^k)$$
其中$H({\mathbf{x}}^k)={\nabla }^{2}f({\mathbf{x}}^k)$，为使下一次迭代使函数值变小，对泰勒展开式求导得
$$\nabla f(\mathbf{x})=\nabla f({\mathbf{x}}^k)+H({\mathbf{x}}^k)(\mathbf{x}-{\mathbf{x}}^k)$$
令$\nabla f(\mathbf{x})=0$得
$$\mathbf{x}={\mathbf{x}}^k-H({\mathbf{x}}^k{)}^{-1}\nabla f({\mathbf{x}}^k)$$
因此第$k+1$次的迭代公式为
$${\mathbf{x}}^{k+1}={\mathbf{x}}^k-H({\mathbf{x}}^k{)}^{-1}\nabla f({\mathbf{x}}^k)$$
接下来考虑使用牛顿迭代法求解对数几率回归，对数几率回归可以写成如下优化问题
$$\underset{\mathbf{\beta }}{\min }\sum _{i=1}^m(-y_i{\mathbf{\beta }}^T{\hat{\mathbf{x}}}_i+ln(1+e^{{\mathbf{\beta }}^T{\hat{\mathbf{x}}}_i}))$$
令$l(\mathbf{\beta })={\sum }_{i=1}^m(-y_i{\mathbf{\beta }}^T\hat{\mathbf{x}}+ln(1+e^{{\mathbf{\beta }}^T\hat{\mathbf{x}}}))$，首先求梯度的表达式,
$$\frac{\partial l(\mathbf{\beta })}{\partial \mathbf{\beta }}=\sum _{i=1}^m(-y_i{\hat{\mathbf{x}}}_i+\frac{{\hat{\mathbf{x}}}_i{e}^{{\mathbf{\beta }}^T{\hat{\mathbf{x}}}_i}}{1+e^{{\mathbf{\beta }}^T{\hat{\mathbf{x}}}_i}})$$
其次计算$Hessian$矩阵的表达式
$$\frac{{\partial }^{2}l(\mathbf{\beta })}{\partial \mathbf{\beta }\partial {\mathbf{\beta }}^T}=\sum _{i=1}^m\frac{e^{{\mathbf{\beta }}^T{\hat{\mathbf{x}}}_i}}{(1+e^{{\mathbf{\beta }}^T{\hat{\mathbf{x}}}_i}{)}^2}{\hat{\mathbf{x}}}_i{\hat{\mathbf{x}}}_i^T$$
最后给出使用牛顿迭代法求解对数几率回归的python实现

import numpy as np
from sklearn import datasets
import matplotlib.pyplot as plt


class LogitRegression:
    def __init__(self):
        self._beta = None
        self.coef_ = None
        self.interception_ = None

    def _sigmoid(self, x):
        return 1. / (1. + np.exp(-x))

    def fit(self, X, y, eta=0.1, n_iters=1e3):
        assert X.shape[0] == y.shape[0]

        def dL(beta, X_b, y):
            try:
                vec = beta.dot(X_b.T)
                coef = np.exp(vec) / (1 + np.exp(vec))
                return (-y.dot(X_b) + X_b.T.dot(coef)) / len(X_b)
            except Exception:
                return float('inf')

        def Hessian(beta, X_b):
            H = np.zeros(shape=(X_b.shape[1], X_b.shape[1]))
            for x in X_b:
                e = np.exp(beta.dot(x))
                H += e / (e ** 2) * x[:, np.newaxis].dot(x[np.newaxis])
            return H

        def newton_method(X_b, y, n_iters, initial_beta, eta, eps):
            beta = initial_beta
            i_iter = 0
            while i_iter < n_iters:
                gradient = dL(beta, X_b, y)
                if np.sum(np.abs(gradient)) < eps:
                    break
                H = Hessian(beta, X_b)
                beta = beta - eta * np.linalg.inv(H).dot(gradient)
                i_iter += 1
            return beta

        X_b = np.hstack([np.ones((len(X), 1)), X])
        initial_beta = np.zeros(X_b.shape[1])
        eps = X.shape[1] / 10
        self._beta = newton_method(X_b, y, n_iters, initial_beta, eta, eps)
        self.coef_ = self._beta[1:]
        self.interception_ = self._beta[0]
        return self

    def predict_probability(self, X_predict):
        X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
        return self._sigmoid(X_b.dot(self._beta))

    def predict(self, X_predict):
        assert self.coef_ is not None and self.interception_ is not None
        assert X_predict.shape[1] == len(self.coef_)
        probability = self.predict_probability(X_predict)
        return np.array(probability >= 0.5, dtype='int')

    def score(self, X_test, y_test):
        y_predict = self.predict(X_test)
        return np.sum(y_predict == y_test) / len(y_test)

    def __repr__(self):
        return "LogitRegression()"


if __name__ == "__main__":
    iris = datasets.load_iris()
    X = iris.data[:, :2]
    y = iris.target
    X = X[y < 2]
    y = y[y < 2]
    logit = LogitRegression().fit(X, y)
    y_predict = logit.predict(X)
    print(np.sum(y_predict == y))
    w, b = logit.coef_, logit.interception_

    def x2(w, b, x1):
        return -(w[0] * x1 + b) / w[1]

    plt.figure()
    plt.scatter(X[y == 0, 0], X[y == 0, 1], c='r')
    plt.scatter(X[y == 1, 0], X[y == 1, 1], c='b')
    xmin, xmax = np.min(X[:, 0]), np.max(X[:, 0])
    plt.plot([xmin, xmax], [x2(w, b, xmin), x2(w, b, xmax)], label='boundary')
    plt.title('Boundary of LogitRegression')
    plt.legend()
    plt.show()

扩展：逆牛顿法

牛顿迭代法每一次都要计算$Hessian$矩阵的逆，计算量较大，并且$Hessian$矩阵可能不可逆，逆牛顿法的基本思想是使用一个矩阵${\mathbf{G}}_k=G({\mathbf{x}}^{\mathbf{k}})$来近似代替$H({\mathbf{x}}^k{)}^{-1}$。

牛顿迭代法中的$H({\mathbf{x}}^k)$满足以下条件
$$\nabla f({\mathbf{x}}^{k+1})-\nabla f({\mathbf{x}}^k)=H({\mathbf{x}}^k)({\mathbf{x}}^{k+1}-{\mathbf{x}}^k)$$
等价于
$$H({\mathbf{x}}^k{)}^{-1}(\nabla f({\mathbf{x}}^{k+1})-\nabla f({\mathbf{x}}^k))=({\mathbf{x}}^{k+1}-{\mathbf{x}}^k)$$
这称为拟牛顿条件，逆牛顿法将${\mathbf{G}}_k$作为$H({\mathbf{x}}^k{)}^{-1}$的近似，要求矩阵${\mathbf{G}}_k$正定，并且满足
$${\mathbf{G}}_{k+1}(\nabla f({\mathbf{x}}^{k+1})-\nabla f({\mathbf{x}}^k))=({\mathbf{x}}^{k+1}-{\mathbf{x}}^k)$$
不同的逆牛顿法，区别就在于如何确定$\mathbf{G}$。