本文共 8370 字，大约阅读时间需要 27 分钟。

1、Softmax模型
    Softmax回归和逻辑回归的回归模型很类似，为如下条件概率：

$$P(Y=k|x_i;w_k)=\frac{\exp(w_k\cdot x_i)}{\sum_{k=1}^K\exp(w_k\cdot x_i)}\quad$$
其中$x\in R^{n+1}$，$w_k\in R^{n+1}$，$k=1,2,3...,K$，表示有K个类别。
    对于给定的测试输入$x_i$，我们想用条件概率针对每一个类别$k$估算出概率值 $P(Y=k | x_i)$。也就是说，我们想估计$x_i$ 的每一种分类结果出现的概率。因此，我们的条件概率函数将要输出一个 $K$维的向量（向量元素的和为1）来表示这$K$ 个估计的概率值。 具体地说，我们的条件概率函数$P_w(x_i)$如下：
$$P_w(x_i)=\begin{bmatrix}P(Y=1|x_i;w_1)\\ P(Y=2|x_i;w_2)\\ ...\\ P(Y=K|x_i;w_K)\end{bmatrix}= \frac{1}{\sum_{k=1}^K\exp(w_k\cdot x_i)}\quad \begin{bmatrix}\exp(w_1\cdot x_i)\\ \exp(w_2\cdot x_i)\\ ...\\ \exp(w_K\cdot x_i)\end{bmatrix}$$
其中$w_1,w_2,...,w_K\in R^{n+1}$。
    为了方便起见，我们同样使用符号 $w$来表示全部的模型参数。在实现Softmax回归时，将$w$用一个$k*(n+1)$ 的矩阵来表示会很方便，该矩阵是将$w_1,w_2,...,w_K$ 按行罗列起来得到的，如下所示：
$$w=\begin{bmatrix}w_1\\ w_2\\ ...\\ w_K\end{bmatrix}$$

2、模型参数估计
在下面的公式中，$1\{ \}$ 是示性函数，其取值规则为：$1\{ 值为假的表达式 \}=0,1\{ 值为真的表达式 \}=1$。举例来说，表达式$1\{2+2=4\}$ 的值为1 ，$1\{1+1=5\}$的值为 0。我们对数似然函数为：

$$L(w)=\sum_{i=1}^N\sum_{k=1}^K1\{y_i=k\}\log\frac{\exp(w_k\cdot x_i)}{\sum_{l=1}^K\exp(w_l\cdot x_i)}\quad$$
值得注意的是，上述公式是logistic回归代价函数的推广。logistic回归代价函数可以改为：
$$L(w)=\sum_{i=1}^N[y_i\log\pi(x_i)+(1-y_i)\log(1-\pi(x_i))] \\=\sum_{i=1}^N\sum_{k=1}^K1\{y_i=k\}\log P(y_i=k|x_i;w_k)$$
    对于$L(w)$的最小化问题，目前还没有闭式解法。因此，我们使用迭代的优化算法（例如梯度下降法）。经过求导，我们得到梯度公式如下：
$$\frac{\partial L(w)}{\partial w_k}\quad=\sum_{i=1}^N[x_i(1\{y_i=k\}-P(y_i=k|x_i;w_k))]$$
其中$\frac{\partial L(w)}{\partial w_k}\quad$本身是一个向量,它的第$l$个元素 $\frac{\partial L(w)}{\partial w_{kl}}\quad$是$L(w)$对$w_k$的第$l$个分量的偏导数。
    有了上面的偏导数公式以后，我们就可以将它代入到梯度下降法等算法中，来最小化 $L(w)$。 例如，在梯度下降法的标准实现中，每一次迭代需要进行如下更新:
$$w_k = w_k - \alpha\frac{\partial L(w)}{\partial w_k}\quad( k=1,\ldots,K）$$

3、Softmax回归 vs 二元回归
    当类别数 $k = 2$ 时，softmax 回归退化为 logistic 回归。这表明 softmax 回归是 logistic 回归的一般形式。如果你在开发一个音乐分类的应用，需要对k种类型的音乐进行识别，那么是选择使用 softmax 分类器呢，还是使用 logistic 回归算法建立 k 个独立的二元分类器呢？
    这一选择取决于你的类别之间是否互斥，例如，如果你有四个类别的音乐，分别为：古典音乐、乡村音乐、摇滚乐和爵士乐，那么你可以假设每个训练样本只会被打上一个标签（即：一首歌只能属于这四种音乐类型的其中一种），此时你应该使用类别数 k = 4 的softmax回归。（如果在你的数据集中，有的歌曲不属于以上四类的其中任何一类，那么你可以添加一个“其他类”，并将类别数 k 设为5。）
    如果你的四个类别如下：人声音乐、舞曲、影视原声、流行歌曲，那么这些类别之间并不是互斥的。例如：一首歌曲可以来源于影视原声，同时也包含人声 。这种情况下，使用4个二分类的 logistic 回归分类器更为合适。这样，对于每个新的音乐作品 ，我们的算法可以分别判断它是否属于各个类别。
    现在我们来看一个计算视觉领域的例子，你的任务是将图像分到三个不同类别中。(i) 假设这三个类别分别是：室内场景、户外城区场景、户外荒野场景。你会使用sofmax回归还是 3个logistic 回归分类器呢？ (ii) 现在假设这三个类别分别是室内场景、黑白图片、包含人物的图片，你又会选择 softmax 回归还是多个 logistic 回归分类器呢？
    在第一个例子中，三个类别是互斥的，因此更适于选择softmax回归分类器 。而在第二个例子中，建立三个独立的 logistic回归分类器更加合适。

4、Softmax回归的实现：
    这里用的数据集仍然为上一篇逻辑回归中的数据集：

#coding=utf-8from numpy import *import matplotlib.pyplot as pltclass SoftmaxRegression:    def __init__(self):        self.dataMat = []   # 数据集        self.labelMat = []  # 类标签        self.weights = []   # 权重        self.M = 0          # 样本点个数        self.N = 0          # 样本特征个数+1        self.K = 0          # 保存累别的总数        self.alpha = 0.001  # 学习率    #加载数据集    def loadDataSet(self):        for line in open('softmax.txt','r'):            items = line.strip().split('    ')            self.dataMat.append([1.0, float(items[0]), float(items[1])])            self.labelMat.append(int(items[2]))        self.K = len(set(self.labelMat))                self.dataMat = mat(self.dataMat)        self.labelMat = mat(self.labelMat).transpose()        self.M,self.N = shape(self.dataMat)        self.weights = mat(ones((self.N,self.K)))  # 权重初始化为全1的 N*K的矩阵,每一列为一个权重向量        # self.weights = [[-1.19792777,6.05913226,-4.44164147,3.58043698],        #                 [ 1.78758743,0.47379819,0.63335518,1.1052592 ],        #                 [ 1.48741185,-0.18748907,1.79339685,0.90668037]]    def likelihoodfunc(self):        likelihood = 0.0        for i in range(self.M):                    # 对所有样本点迭代            t = exp(self.dataMat[i]*self.weights)            likelihood += log(t[0,self.labelMat[i,0]]/sum(t))        print likelihood    # 梯度下降法，每次更新权重时选择所有样本点    def gradientAscent(self):        for l in range(10):                        # 进行迭代10次            error = exp(self.dataMat*self.weights) # 得到一个矩阵，每个元素代表一个样本点和一个权值的内积            rowsum = -error.sum(axis=1)            # 矩阵每一行的所有分量相加，得到一列。每个元素对应条件概率公式中的分母            rowsum = rowsum.repeat(self.K, axis=1) # 对rowsum进行 列方向上的扩展，结果为 N行K列的矩阵          #求出每个样本点属于每个类别的条件概率，即条件概率模型的公式，矩阵第几行代表第几个样本点，第几列代表第几类            error = error/rowsum                              for m in range(self.M):        # 对所有样本点进行迭代                error[m,self.labelMat[m,0]] += 1   # 如果某个样本点属于某个类，对应元素+1，这对应公式中1{ }示性函数            self.weights = self.weights + self.alpha * self.dataMat.transpose()* error  # 跟新权重向量            self.likelihoodfunc()        print self.weights   # 每一列为一个权重向量    # 随机梯度上升法version0，每次按顺序选一个样本点更新权重    def stochasticGradientAscent_V0(self):        for l in range(500):               # 进行迭代500次            for i in range(self.M):        # 对每个样本进行遍历                error = exp(self.dataMat[i]*self.weights)                rowsum = -error.sum(axis=1)                rowsum = rowsum.repeat(self.K, axis=1)                error = error/rowsum                error[0,self.labelMat[i,0]] += 1                self.weights = self.weights + self.alpha * self.dataMat[i].transpose()* error                  # self.likelihoodfunc()        print self.weights    # 随机梯度下降法version1，每次随机选一个样本点更新权重    def stochasticGradientAscent_V1(self):        for l in range(500):            idxs = range(self.M)            for i in range(self.M):                alpha = 4.0/(1.0+l+i)+0.01                rdmidx = int(random.uniform(0,len(idxs)))                error = exp(self.dataMat[rdmidx]*self.weights)                rowsum = -error.sum(axis=1)                rowsum = rowsum.repeat(self.K, axis=1)                error = error/rowsum                error[0,self.labelMat[rdmidx,0]] += 1                self.weights = self.weights + alpha * self.dataMat[rdmidx].transpose()* error                del(idxs[rdmidx])                # self.likelihoodfunc()        print self.weights    def classify(self,X):        p = X * self.weights        return p.argmax(1)[0,0]    def test(self):        xcord0 = []        ycord0 = []        xcord1 = []        ycord1 = []        xcord2 = []        ycord2 = []        xcord3 = []        ycord3 = []        for i in range(50):            for i in arange(80):                x = random.uniform(-3.0, 3.0)                y = random.uniform(0.0, 15.0)                c = self.classify(mat([[1.0,x,y]]))                if c==0:                    xcord0.append(x)                    ycord0.append(y)                if c==1:                    xcord1.append(x)                    ycord1.append(y)                if c==2:                    xcord2.append(x)                    ycord2.append(y)                if c==3:                    xcord3.append(x)                    ycord3.append(y)        # for i in range(self.M):            # if self.labelMat[i]==0:                # xcord0.append(self.dataMat[i,1]); ycord0.append(self.dataMat[i,2])            # elif self.labelMat[i]==1:                # xcord1.append(self.dataMat[i,1]); ycord1.append(self.dataMat[i,2])            # elif self.labelMat[i]==2:                # xcord2.append(self.dataMat[i,1]); ycord2.append(self.dataMat[i,2])            # else:                # xcord3.append(self.dataMat[i,1]); ycord3.append(self.dataMat[i,2])        fig = plt.figure()        ax = fig.add_subplot(111)        ax.scatter(xcord0, ycord0, s=20, c='yellow', marker='s')        ax.scatter(xcord1, ycord1, s=20, c='blue')        ax.scatter(xcord2, ycord2, s=20, c='red')        ax.scatter(xcord3, ycord3, s=20, c='black')        plt.title('inference')        # plt.title('train data')        plt.xlabel('X1')        plt.ylabel('X2');        plt.show()if __name__=='__main__':        myclassification = SoftmaxRegression()    myclassification.loadDataSet()       # myclassification.gradientAscent()    myclassification.stochasticGradientAscent_V0()       # myclassification.stochasticGradientAscent_V1()    myclassification.test()