《机器学习（周志华）》学习笔记（三）

Q：机器学习中最简单的学习算法是什么？

最简单的机器学习算法莫过于线性回归算法了。线性回归算法的基本形式如下：

$f(x)=w_1 x_1+w_2 x_2+···+w_n x_n + b$

如果说把 $b$ 换成 $w_0$ ，我们能得到一个更简洁的公式

$f(x)=w_1 x_1 +w_2 x_2+···+w_n x_n + w_0$

其中x是自变量，w是参数。我们训练的目的就是找到一组 w，使得我们输入一组x的值后得到的f(x)，也就是我们想要预测的量（房子价格、西瓜类别等）。一个常见的例子是房价预测例子。通过住房面积、房间个数两个特征来预测房价，则f(x)就是代表房价，x0=1,x1=住房面积，x2=房间个数。

Q：线性回归模型的参数w和b怎样确定？

一个预测模型的效果越好，意味着其预测误差越小。那我们的任务就是找出使得预测值f(x)与真实值y的差最小的参数w。那么我们可以建立一个函数：

$J(w) = \sum_i [f(x_i|w)-y_i ]^2$

这个其实是在取定一组参数w和b的情况下，把每一次的预测值与真实值的差的平方加起来，得到了总的均方误差。令这个函数取最小值，那么就得到了能让f(x)与y的差最小，亦即我们的预测值与真实值的总误差最小的参数w和b，我们的回归模型效果也就最好了。

现在问题转化为最小化上面的损失函数 $J(w)$ 。那么这个函数的最小化又要怎么做呢？

高中的数学告诉我们，求一个函数的最小值，可以先求导，令导数为0，就找到了一堆极值点。然后把极值点代入原函数，算出函数值最小的点就是对应的最小值点了。

求导找最值

w的最优值

b的最值

因为损失函数 $J(w)$ 是一个凸函数，所以求出来的解直接就是最优解了。如此我们就获得了一个可以使用的线性回归模型。

我们可以把上面的情况写成更一般，更简洁的形式。考察下面线性回归模型的一般形式

$f(x)=w_1 x_1 +w_2 x_2+···+w_n x_n + w_0$

把w写成向量，x写成矩阵，则有

$f(\textbf{X}) = \textbf{X}\vec{w}$

对于训练样本，我们可以得出下面等式，其中X时训练数据各种特征值，y时训练样本对应的真实值，w是待求的参数：

$\textbf{X}\vec{w} = \vec{y}$

根据我们大一的线性代数知识，可以知道当X为方阵时，有

$\vec{w} = \textbf{X}^{-1}\vec{y}$

更一般的，即便X不为方阵，也有：

$\vec{w} = (\textbf{X}^T\textbf{X})^{-1}\textbf{X}^T\vec{y}$

如此我们便得出了直接求出了w的最优值。

这样直接使用微积分求导或者解线性方程组一次性得到结果的方法被称为analytical方法。与之构成对比的是numerical方法，这种方法通常需要通过迭代不断逼近，而且得到的结果往往只是近似值，比如将会在下面提到的梯度下降法。

Q：线性回归模型是用来做预测的，但做分类任务时应该用哪一个模型？

上面提到的线性模型是最基本的线性模型，把这个基本线性模型加以改造，就可以得到用于分类任务的回归模型。

假设我们输入一组住房面积和房间个数的数据，预测到了一组住房的房价，比如A房80万，B房135万，C房150万等。如果我们认为画一条线，认为140万以下为低档房，140万以上为高档房，那我们就可以根据住房面积和房间个数来给住房分类了。也就是说在预测的基础上加上一个判断条件，就构成了分类算法。

对数几率回归模型，或者称逻辑回归（Logit Regression）就是这样一种基于回归的分类算法。

$h_w(x) = g(w^Tx) = \frac{1}{1+e^{-(w_0+w_1x_1...+w_nx_n)}}$

一般来说中文对Logistic Regression的翻译多是“逻辑斯蒂回归”，也就是音译，不过周志华老师这里的“对数几率回归”则是意译（又有对数，又有概率解释，还是回归方法）。然而我自己也没找到Logistic Regression这个词的内涵，维基百科也是：Verhulst did not explain the choice of the term "logistic", but it is presumably in contrast to the logarithmic curve, and by analogy with arithmetic and geometric.

对数几率回归函数输出的值是0到1范围的数。对于一个二分类任务来说，假如只分正例和反例（好瓜和坏瓜），那么我们可以认为函数输出的值越大，则这个样本是正例的几率越大，值越小，样本是反例的可能越大。这就是“对数几率回归”名字的由来。我们可以定一个临界值，一般是0.5，当输出值大于0.5时认为该样本是正例，否则为反例。

Q：对数几率回归的参数w和b又是怎样确定的？

我们知道，对数几率回归函数的输出值是一个概率。那么一个对率回归模型越好，意味着当输入模型的样本为正例时，模型的输出值越大（越接近1），当输入模型的样本为反例时，模型的输出值越小（越接近0）。

所以我们可以构建这个函数，具体的推论参考这篇文章
$l(w)=\sum_{i=1}^{m}[y_iln(h_w(x))+(1-y_i)ln(1-h_w(x))]$

上式中 $h_w(x)$ 是对率模型的输出值，亦即x为正例的概率值，yi是样本的真实标记，当样本为正例时y取1，为反例时y取0。这整一个函数的含义就是把每个样本属于其真实标记的概率加起来。若我们使 l(w,b) 最大，那就可以找到一组θ使得每个样本属于其真实标记的概率最大，亦即使得正例的函数值尽可能接近1，反例的函数值尽可能接近0。现在问题又转化成了最优化问题，只是这一次的函数没有上面线性回归函数那么好，很难用求导的方式确定最优解。

其实最优化问题是一类经典的工程数学问题，现在已经有了许多种成熟的更泛用解决方案，比如随机游走，比如遗传算法，但在机器学习领域最常用的还是梯度下降方法。可能是因为梯度下降比较容易理解，效果也比较确定吧？

“梯度下降”中的梯度又是一个什么东西？梯度是一个向量，其方向表示函数变化最快的方向，其大小（模）表示函数在此方向上的变化率。计算函数在某一点上的梯度，就是要算出函数在该点上的所有维度上的偏导数。简单理解就是，梯度是一元函数下的导数在多元函数下的推广（虽然并不准确）。

weikipedia image

梯度下降法是一个迭代进行的方法，每一轮都朝着梯度方向走动一小步，逐渐逼近最优值。

// psudocode for gradient descent
procedure grad_desc(X, y, w, learning_rate, precision)
    prev_w = w
    do
        w = w + learning_rate * gradient(w, X, y)
    while f(prev_w) - f(w) < precision
    return w
end

如此我们就得到了一个优秀的对数几率回归模型。

Q：对数几率回归模型是处理二分类任务的，那处理多分类任务时怎么办?

多分类任务可以凭借二分类任务为基础得到解决。解决思路一般有一对一，一对多，和多对多三种。

一对一策略：举个例子，现在有A、B、C、D四个种类，要确定某个样本属于四类中的哪一类。那么我们可以事先训练好六个二分类的分类器——A/B、A/C、A/D、B/C、B/D、C/D。然后把要确定类别的样本分别放入这六个分类器。假设分类结果分别是A、A、A、B、D、C。可以知道六个分类器中有三个认为这个样本属于A，认为是B、C、D的各有一个。所以我们可以认为这个样本就是属于A类的。

一对多策略：举个例子，现在有A、B、C、D四个种类，要确定某个样本属于四类中的哪一类。那么我们可以事先训练好四个二分类的分类器——A/BCD、B/ACD、C/ABD、D/ABC,分类器输出的是一个函数值。然后把要确定类别的样本分别放入这四个分类器。假设四个分类器的结果分别是“属于A的概率是0.9”，“属于B的概率是0.8”、“属于C的概率是0.7”、“属于B的概率是0.6”。那我们可以认为这个样本是属于A。

多对多策略：每次将若干个类作为正类，若干个其他类作为反类。其中需要这种策略需要特殊的设计，不能随意取。常用的技术：纠错输出码。工作步骤分为：

编码：对N个类别做M次划分，每次划分将一部分类别作为正类，一部分划分为反类，从而形成一个二分类训练集；这样一共产生M个训练集，可以训练出M个分类器。
解码：M个分类器分别对测试样本进行预测，这些预测标记组成一个编码，将这个预测编码与每个类别各自的编码进行比较，返回其中距离最小的类别作为最终预测结果。

多对多：编码

以原书的例子做一个详细的演示：

假如现在有一个训练数据集，可以分为四个类——C1, C2, C3, C4

再具体一点可以想象为——西瓜可以分为一等瓜、二等瓜、三等瓜、四等瓜，要训练一个分类系统，使之能判断一个西瓜的等级。

我们对训练数据集做五次划分

第一次，标记C2为正例，其他为反例，训练出一个二分类的分类器f1
第二次，标记C1、C3为正例，其他为反例，训练出一个二分类的分类器f2
第三次，标记C3、C4为正例，其他为反例，训练出一个二分类的分类器f3
第四次，标记C1、C2、C4为正例，其他为反例，训练出一个二分类的分类器f4
第五次，标记C1、C4为正例，其他为反例，训练出一个二分类的分类器f5

根据这五次划分的过程，每一个类都获得了一个编码（向量）：

C1：（-1，1，-1，1，1）
C2：（1，-1，-1，1，-1）
C3：（-1，1，1，-1，1）
C4：（-1，-1，1，1，-1）

若现在有一个测试样本，五个分类器对应的累计结果为

f1：反， f2：反， f3：正， f4：反， f5：正

即该测试样本对应的编码/向量为（-1，-1，1，-1，1）

那么分别计算这个测试样本的编码与四个类别的编码的向量距离，可以使用欧氏距离，算的与C3类的距离是最小的。因此判定该测试样本属于C3类。

Talk is cheap, show me the code!

先来一段Linear Regression的代码


"""
Simple linear regression

:file: supervised.py
:author: Richy Zhu
:email: rickyzhu@foxmail.com
"""
class MyLinearRegression:

    def __init__(self):
        self.w = None

    def _loss_function(self, X, y, w):
        '''mean square as loss function'''
        y_hat = X.dot(w)
        return np.sum((y - y_hat)**2)

    def _gradient_descent(self, X, y, w, eta):
        '''compute new weights by batch gradient descent'''
        y_hat = X.dot(w)
        diff = y - y_hat
        grad = eta * (y - y_hat).reshape([-1,1]) * X * 1.0/len(X)
        w += grad.sum(axis=0)
        return w

    def fit_analitical(self, X, y):
        '''
        Train the logistic regression model, by analytical method

        Parameters
        ----------
        X: ndarray of shape (m, n)
            sample data where row represent sample and column represent feature
        y: ndarray of shape (m,)
            labels of sample data

        Returns
        -------
        self
            trained model
        '''
        X, y = np.array(X), np.array(y)
        X = np.hstack([X, np.ones([len(X),1])]) # add bias
        w = np.linalg.inv(X.T.dot(X)).dot(X.T.dot(y))
        self.b = w[-1]
        self.w = w[:-1]
        return self


    def fit(self, X, y, eta=1e-8, epslon=1e-6, max_iter=100000):
        '''
        Train the simple linear regression model, by gradient descent

        Parameters
        ----------
        X: ndarray of shape (m, n)
            sample data where row represent sample and column represent feature
        y: ndarray of shape (m,)
            labels of sample data
        eta: float
            learning rate
        epslon: float
            terminate when prev_loss - this_loss < epslon
        max_iter: int
            terminate when number of iterations excceed max_iter

        Returns
        -------
        self
            trained model
        '''
        X = np.hstack([X, np.ones([len(X),1])]) # add bias
        # w = np.random.random(X.shape[1]) 
        w = np.ones(X.shape[1])  # standard practice is to randomly initialize
        prev_loss = self._loss_function(X, y, w)
        for i in range(max_iter):
            w = self._gradient_descent(X, y, w, eta)
            this_loss =  self._loss_function(X, y, w)
            if prev_loss - this_loss < epslon:
                break
            prev_loss = this_loss
        self.b = w[-1]
        self.w = w[:-1]
        return self
    
    def predict(self, X):
        '''
        Make prediction by the trained model.

        Parameters
        ----------
        X: ndarray of shape (m, n)
            data to be predicted, the same shape as trainning data

        Returns
        -------
        C: ndarray of shape (m,)
            Predicted value per sample.
        '''
        return X.dot(self.w) + self.b

测试一下


import numpy as np
from sklearn.datasets import load_boston, load_iris
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, accuracy_score
from supervised import MyLinearRegression

print('Linear Regression')
print('---------------------------------------------------------------------')
boston = load_boston()
X_train, X_test, y_train, y_test = train_test_split(boston.data, boston.target)

mylr = MyLinearRegression()
mylr.fit(X_train, y_train)
y_pred = mylr.predict(X_test)
print('My MSE:', mean_squared_error(y_test, y_pred))

mylr.fit_analitical(X_train, y_train)
y_pred = mylr.predict(X_test)
print('My MSE:', mean_squared_error(y_test, y_pred), '(analytical method)')

from sklearn.linear_model import LinearRegression as SKLinearRegression
lr = SKLinearRegression()
lr.fit(X_train, y_train)
y_pred = lr.predict(X_test)
print('Sk MSE:', mean_squared_error(y_test, y_pred))

测试的输出结果如下


$ python supervised_example.py

Linear Regression
---------------------------------------------------------------------
My MSE: 491.6360044772676
My MSE: 17.8423520615301 (analytical method)
Sk MSE: 17.84235206153754

从输出结果可见，我们手写的使用analytical方法的代码居然高度逼近scikit-learn包的结果，可以猜测sklearn的线性回归模型用的也是analytical的方法。在通常情况下，实现线性回归模型也应当优先考虑使用analytical方法。

下面是Logistic Regression的代码


"""
Simple Logistic Regression

:file: supervised.py
:author: Richy Zhu
:email: rickyzhu@foxmail.com
"""

class MyLogisticRegression:

    def __init__(self):
        self.w = None
        self.b = None
    
    def _loss_function(self, X, y, w): 
        '''cross entropy as loss function'''
        y_hat = sigmoid(X.dot(w))
        ll = np.sum(y * np.log(y_hat) - (1-y) * np.log(1-y_hat))
        return ll

    def _gradient_descent(self, X, y, w, eta):
        '''compute new weights by batch gradient descent'''
        y_hat = sigmoid(X.dot(w))
        diff = y - y_hat
        grad = eta * diff.reshape([-1,1]) * X * 1.0/len(X)
        w += grad.sum(axis=0)
        return w

    def fit(self, X, y, eta=1e-5, epslon=1e-6, max_iter=100000):
        '''
        Train the logistic regression classifier model

        Parameters
        ----------
        X: ndarray of shape (m, n)
            sample data where row represent sample and column represent feature
        y: ndarray of shape (m,)
            labels of sample data
        eta: float
            learning rate
        epslon: float
            terminate when prev_loss - this_loss < epslon
        max_iter: int
            terminate when number of iterations excceed max_iter

        Returns
        -------
        self
            trained model
        '''
        X = np.hstack([X, np.ones([len(X),1])]) # add bias
        # standard practice is to randomly initialize
        w = np.ones(X.shape[1])  
        prev_loss =  self._loss_function(X, y, w)
        for i in range(max_iter):
            self.w = self._gradient_descent(X, y, w, eta)
            this_loss = self._loss_function(X, y, w)
            if prev_loss - this_loss < epslon:
                break
            prev_loss = this_loss
        self.b = w[-1]
        self.w = w[:-1]
        return self
    
    def predict(self, X):
        '''
        Make prediction by the trained model.

        Parameters
        ----------
        X: ndarray of shape (m, n)
            data to be predicted, the same shape as trainning data

        Returns
        -------
        C: ndarray of shape (m,)
            Predicted class label per sample.
        '''
        if self.w is None:
            raise Exception("Model haven't been trained!")
        # X = np.hstack([X, np.ones([len(X), 1])])
        return binarize(sigmoid(X.dot(self.w)+self.b))