softmax理论及代码解读——UFLDL

前言

看了各种softmax以后迷迷糊糊的，还是研究一下UFLDL的教程稳点。当然还是得参考挺多教程的：UFLDL-softmax 、Softmax的理解与应用 、Logistic 分类器与 softmax分类器 、详解softmax函数以及相关求导过程 、Exercise:Softmax Regression 。
【UFLDL】网址有所变动：戳这里

其实，在最初比较疑惑的问题是：softmax的回归和分类到底是不是一个解法，或者它俩是不是一个知识？然后在一些博客中发现，softmax有时候也称为softmax分类回归器，而且在机器学习的定义中，回归问题通常用来预测一个值，如预测房价或者未来的天气状况等；分类用于将事物打上一个标签，通常结果为离散值，但是分类通常建立在回归上的，最常见的分类方法是逻辑回归或者叫逻辑分类。详细请戳回归(regression)与分类(classification)的区别

还是将logistic回归和softmax回归放一起看看。

Logistic回归

需要注意的是，logistic回归的中文叫做对数回归 ，我们经常也称他为逻辑回归等

广义线性模型

上面说了分类通常建立在回归中，那么这个“建立”方法是什么？周志华老师《机器学习》给出了一个答案说的非常好：“只需要找一个单调可微函数将分类任务的真实标记与线性回归模型的预测值联系起来，上公式
$$y=g^{-1}({\mathbf{w}}^T\mathbf{x}+b)$$
这样得到的一个模型$y$称为"广义线性模型(generalized linear model)",其中$g(\cdot )$是单调可微函数，称为联系函数(link function)。对数广义函数$g(\cdot )=ln(\cdot )$是广义线性函数的一个特例

logistic回归(对数回归)

出现原因是，在二分类问题中，线性预测$z={\mathbf{w}}^T\mathbf{x}+b$的输出值$z$为实值，转换为二分类最理想的方法是使用阶跃函数
$$y=\{\begin{array}{l}0,z<0\\ 0.5,z=0\\ 1,z>0\end{array}$$
但是阶跃函数明显是不可导的，需要采用近似的函数代替，这样对数几率函数就出现了(logistic function)
$$y=\frac{1}{1+e^{-z}}$$
它与我们经常说的sigmoid函数一样。带入到广义线性模型中就得到了对数回归(logistic regression)，说是回归，其实是分类，而且对数回归的函数也经常被当做样本$\mathbf{x}$类别为1的概率：
$$P(y=1|\mathbf{x};\mathbf{w},b)=h_{\theta }(\mathbf{x})=\frac{1}{1+e^{-({\mathbf{w}}^T\mathbf{x}+b)}}$$
而且这个函数上下同时乘以$e^{(w^{T}x+b)}$就能够变形成
$$P(y=1|\mathbf{x};\mathbf{w},b)=\frac{e^{({\mathbf{w}}^T\mathbf{x}+b)}}{1+e^{({\mathbf{w}}^T\mathbf{x}+b)}}$$
据此，在UFLDL中定义了对数回归的损失函数
$$J(\theta )=-\frac{1}{m}\left[\sum _{i=1}^m{y}^{(i)}\log h_{\theta }(x^{(i)})+(1-y^{(i)})\log (1-h_{\theta }(x^{(i)}))\right]$$
更逼格一点的写法是
J ( θ ) = − 1 m [ ∑ i = 1 m ∑ j = 0 1 1 { y ( i ) = j } log ⁡ P ( y ( i ) = j ∣ x ( i ) ; θ ) ] J(\theta)=-\frac{1}{m}\left[\sum_{i=1}^m\sum_{j=0}^11\{y^{(i)}=j\}\log P(y^{(i)}=j|x^{(i)};\theta) \right] J(θ)=−m1​[i=1∑m​j=0∑1​1{y(i)=j}logP(y(i)=j∣x(i);θ)]
式中的 1 { ⋅ } 1\{\cdot\} 1{⋅}称为示性函数，大括号里面的值为真的时候表达式值取1，反之取0；其中$i$表示第几个样本，样本总数为$m$ ，因为是有监督训练训练更新嘛，那么$y^{(i)}$就是指第$i$个样本的标签了，由于是二分类问题，所以 y ( i ) ∈ { 0 , 1 } y^{(i)}\in\{0,1\} y(i)∈{0,1} ,更新办法就是直接用梯度更新就OK了。

Softmax 回归

梯度求解

由于softmax就是将简单的二分类扩充为多分类， 所以只需要对logistic回归代价函数中的两类问题改为多分类即可，改法就是观察logistic回归逼格写法中的第二个累加和的范围是从0到1，那么多分类就是从0到类别总数了，这样就可以得到softmax的代价函数
J ( θ ) = − 1 m [ ∑ i = 1 m ∑ j = 1 k 1 { y ( i ) = j } log ⁡ P ( y ( i ) = j ∣ x ( i ) ; θ ) ] J(\theta)=-\frac{1}{m}\left[\sum_{i=1}^m\sum_{j=1}^k1\{y^{(i)}=j\}\log P(y^{(i)}=j|x^{(i)};\theta) \right] J(θ)=−m1​[i=1∑m​j=1∑k​1{y(i)=j}logP(y(i)=j∣x(i);θ)]
由于softmax 输出的是对每一个类别估算的概率值，为了保证概率和为1，需要进行归一化，所以条件概率函数,也就是知道样本输入和模型参数的时候，这个样本属于各个标签的概率为：
$$P(y^{(i)}=l|x^{(i)};\theta )=\frac{e^{{\theta }_l^T{x}^{(i)}}}{{\sum }_{j=1}^k{e}^{{\theta }_j^T{x}^{(i)}}}$$
需要注意的是，这里的$\theta$代表全部模型参数，大小是$k\ast (n+1)$
$$\theta ={\left[{\theta }_1^T,{\theta }_2^T,\cdots ,{\theta }_k^T\right]}^T$$
其实这样一看，就感觉应该是每一个类别都有一组模型参数，而每一组模型参数大小就是输入维度或者权重维度$n$再加上偏置的维度1,跟简单的神经网络连接方法一样，右边的每一个标签输出与左边的输入都有权重连接。

然后便可以得到softmax更具体的代价函数
J ( θ ) = − 1 m [ ∑ i = 1 m ∑ j = 1 k 1 { y ( i ) = j } log ⁡ e θ l T x ( i ) ∑ j = 1 k e θ j T x ( i ) ] J(\theta)=-\frac{1}{m}\left[\sum_{i=1}^m\sum_{j=1}^k1\{y^{(i)}=j\}\log \frac{e^{\theta_l^Tx^{(i)}}}{\sum_{j=1}^ke^{\theta_j^Tx^{(i)}}}\right] J(θ)=−m1​[i=1∑m​j=1∑k​1{y(i)=j}log∑j=1k​eθjT​x(i)eθlT​x(i)​]
插曲：上式是多类分类情况，如果利用sigmoid激活到(0,1)范围，获得最后一层值，如果改成多标签分类的表达式就是:
$$loss=\sum _{j=1}^m\sum _{i=1}^n-y_{ji}log(\widehat{y_{ji}})-(1-y_{ji})log(1-\widehat{y_{ji}})$$
其中$y^{(i)}$是真实标签，$P(y^{(i)})$是预测输出的概率，通常为经过sigmoid激活函数归一化到(0,1)后的值。
关于多标签分类与单标签分类的区别，戳这里

直接令其对权重或者偏置的偏导数为零并不能求出闭式解(解析解)，解析解就是有严格的公式可以直接根据输入得到输出，就像一元二次方程一样，有确定的求解公式，这一点可以参考百度百科。那么也就只能找近似解了，称为数值解，采用逼近的方法求解，一般就是梯度下降法或者变种方法了。

说到这个梯度，可以参考Softmax公式推导 ，或者看看本渣的推导，可能比他那个简单，但正确与否谢谢大家指正。

求导的时候主要注意一点，是针对某一个类别输出求导，也就是说对应k个$\theta$ 中的一个，所以在对某个权重求解偏导$\frac{\partial J(\theta )}{{\theta }_l}$的时候可以直接对大括号后面一部分求导，即
$$\begin{array}{rl}& \frac{\partial \log \frac{e^{{\theta }_l^T{x}^{(i)}}}{{\sum }_{j=1}^k{e}^{{\theta }_j\cdot x^{(i)}}}}{\partial {\theta }_l}\\ =& \frac{{\sum }_{j=1}^k{e}^{{\theta }_j\cdot x^{(i)}}}{e^{{\theta }_l^T{x}^{(i)}}}\cdot \frac{\partial \frac{e^{{\theta }_l^T{x}^{(i)}}}{{\sum }_{j=1}^k{e}^{{\theta }_j\cdot x^{(i)}}}}{\partial {\theta }_l}\\ =& \frac{{\sum }_{j=1}^k{e}^{{\theta }_j\cdot x^{(i)}}}{e^{{\theta }_l^T{x}^{(i)}}}\cdot \frac{e^{{\theta }_l^T{x}^{(i)}}\cdot x^{(i)}\cdot {\sum }_{j=1}^k{e}^{{\theta }_j\cdot x^{(i)}}-e^{{\theta }_l^T{x}^{(i)}}\cdot e^{{\theta }_l^T{x}^{(i)}}\cdot x^{(i)}}{({\sum }_{j=1}^k{e}^{{\theta }_j\cdot x^{(i)}}{)}^2}\\ =& \frac{x^{(i)}\cdot ({\sum }_{j=1}^k{e}^{{\theta }_j\cdot x^{(i)}}-e^{{\theta }_l^T{x}^{(i)}})}{{\sum }_{j=1}^k{e}^{{\theta }_j\cdot x^{(i)}}}\\ =& x^{(i)}\cdot (1-P(y^{(l)}=1|x^{(i)},k,\theta ))\end{array}$$
然后带入到$\frac{\partial J(\theta )}{{\theta }_l}$ 中，就可以得到UFLDL中的那个代价函数
∂ J ( θ ) ∂ θ j = − 1 m [ ∑ i = 1 m x ( i ) ⋅ ( 1 { y ( i ) = j } − P ( y ( i ) = j ∣ x ( i ) ； θ ) ) ] \frac{\partial J(\theta)}{\partial\theta_j}=-\frac{1}{m}\left[\sum_{i=1}^mx^{(i)}\cdot (1\{y^{(i)}=j\}-P(y^{(i)}=j|x^{(i)}；\theta))\right] ∂θj​∂J(θ)​=−m1​[i=1∑m​x(i)⋅(1{y(i)=j}−P(y(i)=j∣x(i)；θ))]
梯度下降法的时候就直接使用
$${\theta }_j={\theta }_j-\alpha \frac{\partial J(\theta )}{\partial {\theta }_j}$$

权重衰减

先看UFLDL中的一个式子
$$\begin{array}{rl}P(y^{(i)}=j|x^{(i)};\theta )& =\frac{e^{({\theta }_j-\psi {)}^T\cdot x^{(i)}}}{{\sum }_{l=1}^k{e}^{({\theta }_j-\psi {)}^T\cdot x^{(i)}}}\\ & =\frac{e^{({\theta }_j{)}^T\cdot x^{(i)}}}{{\sum }_{l=1}^k{e}^{({\theta }_j{)}^T\cdot x^{(i)}}}\end{array}$$
可以看出权重减去一个向量与不减去这个向量得到的结果一样，说明这个求得的权重参数是冗余的，也就是说最小化代价函数$J(\theta )$的解不唯一,如果$\theta$是解，那么$\theta -\psi$也是一个最优解。这就引出了“权重衰减”，用于解决softmax回归的参数冗余带来的数值问题。

新的代价函数为
J ( θ ) = − 1 m [ ∑ i = 1 m ∑ j = 1 k 1 { y ( i ) = j } log ⁡ P ( y ( i ) = j ∣ x ( i ) ; θ ) ] + 1 2 λ ∑ i = 1 k ∑ j = 0 n θ i j 2 J(\theta)=-\frac{1}{m}\left[\sum_{i=1}^m\sum_{j=1}^k1\{y^{(i)}=j\}\log P(y^{(i)}=j|x^{(i)};\theta) \right]+\frac{1}{2}\lambda\sum_{i=1}^k\sum_{j=0}^n\theta_{ij}^2 J(θ)=−m1​[i=1∑m​j=1∑k​1{y(i)=j}logP(y(i)=j∣x(i);θ)]+21​λi=1∑k​j=0∑n​θij2​
这样加入了权重衰减以后，能够保证代价函数$J(\theta )$是严格的凸函数，并且有唯一解。
∂ J ( θ ) ∂ θ j = − 1 m [ ∑ i = 1 m x ( i ) ⋅ ( 1 { y ( i ) = j } − P ( y ( i ) = j ∣ x ( i ) ； θ ) ) ] + λ θ j \frac{\partial J(\theta)}{\partial\theta_j}=-\frac{1}{m}\left[\sum_{i=1}^mx^{(i)}\cdot (1\{y^{(i)}=j\}-P(y^{(i)}=j|x^{(i)}；\theta))\right]+\lambda\theta_j ∂θj​∂J(θ)​=−m1​[i=1∑m​x(i)⋅(1{y(i)=j}−P(y(i)=j∣x(i)；θ))]+λθj​
何时使用softmax？

当所有的类别互斥的时候就选择softmax分类器，比如：古典音乐、乡村音乐、摇滚乐和爵士乐；但是如果是人声音乐、舞曲、影视原声、流行歌曲，那就最好还是用四个二元logistic分类器了，因为一首歌曲可以来源于影视原声，同时也包含人声。

代码解读

参考UFLDL上的顺序编程，文末会给出全部程序。

(1) 四准备

• 手写数字数据文件
• 读取数据的MATLAB程序
• 需要补全的softmax分类器程序
• computeNumericalGradient源码

(2) 初始化参数和读数据

需要初始化的参数是输入大小、总类别数目、权重衰减、初始模型参数$\theta$

inputSize = 28 * 28; % Size of input vector (MNIST images are 28x28)
numClasses = 10;     % Number of classes (MNIST images fall into 10 classes)
lambda = 1e-4; % Weight decay parameter
% Randomly initialise theta
theta = 0.005 * randn(numClasses * inputSize, 1);
theta = reshape(theta, numClasses, inputSize);

读数据有提供实例，主要注意的是将标签0改为10，读取数据的时候已经自动归一化到

addpath ../data/mnist
images = loadMNISTImages('train-images-idx3-ubyte');
labels = loadMNISTLabels('train-labels-idx1-ubyte');
labels(labels==0) = 10; % Remap 0 to 10

(3) 损失函数的书写

在官方代码中，损失函数并未补全，但是提供了实现技巧:

• 计算真实值矩阵(ground truth matrix )

其实就是每个样本的示性函数，matlab中有自带的函数可以做到这一点，利用sparse可以将对应样本的对应标签标注出来，比如labels是$60000\times 1$的标签向量，标记了60000个样本分属标签，numCases是总样本数目60000，那么sparse(labels,1:numCases,1)表示的就是一个60000行的稀疏矩阵，最后一行是(8,60000) 1意思是第60000个样本对应第八个标签；为了将稀疏矩阵转换为$10\times 60000$的单热度(one hot)编码，那就需要对稀疏矩阵调用一个full函数，其实计算真实值矩阵就是一句话

groundTruth = full(sparse(labels, 1:numCases, 1));

• 计算梯度函数

∂ J ( θ ) ∂ θ j = − 1 m [ ∑ i = 1 m x ( i ) ⋅ ( 1 { y ( i ) = j } − P ( y ( i ) = j ∣ x ( i ) ； θ ) ) ] + λ θ j = − 1 m [ ∑ i = 1 m x ( i ) ⋅ ( 1 { y ( i ) = j } − e ( θ j ) T ⋅ x ( i ) ∑ l = 1 k e ( θ j ) T ⋅ x ( i ) ] + λ θ j \begin{aligned} \frac{\partial J(\theta)}{\partial\theta_j}&=-\frac{1}{m}\left[\sum_{i=1}^mx^{(i)}\cdot (1\{y^{(i)}=j\}-P(y^{(i)}=j|x^{(i)}；\theta))\right]+\lambda\theta_j\\ &=-\frac{1}{m}\left[\sum_{i=1}^mx^{(i)}\cdot (1\{y^{(i)}=j\}-\frac{e^{(\theta_j)^T\cdot x^{(i)}}}{\sum_{l=1}^k e^{(\theta_j)^T\cdot x^{(i)}}}\right]+\lambda\theta_j \end{aligned} ∂θj​∂J(θ)​​=−m1​[i=1∑m​x(i)⋅(1{y(i)=j}−P(y(i)=j∣x(i)；θ))]+λθj​=−m1​[i=1∑m​x(i)⋅(1{y(i)=j}−∑l=1k​e(θj​)T⋅x(i)e(θj​)T⋅x(i)​]+λθj​​

式子是分别对连接到每一个输出标签的权重求梯度，这里直接用矩阵把所有的权重梯度一次性求出来，先看看参数$\theta$是$10\times 784$的矩阵，data是$784\times 60000$的矩阵，第一维是每张图片的大小$28\times 28$被拉长为一个向量，groundTruth是$10\times 60000$的矩阵，梯度求解三个步骤：求分母，求分数项，求这个损失函数梯度(可以不用人工计算梯度，方法看下一个小步骤)

M = exp(theta*data);
M = bsxfun(@rdivide, M, sum(M));
thetagrad = (-1/size(data,2)).*(((groundTruth-(M))*data') + (lambda*theta));

顺便也可以计算一下损失函数

cost = (-1/size(data,2)).*sum(sum(groundTruth.*log(M)));
cost = cost + (lambda/2).*sum(sum(theta.^2));

为了后续方便，为这一步骤定义一个函数称为

function [cost, grad] = softmaxCost(theta, numClasses, inputSize, lambda, data, labels)

(4) matlab自动梯度优化

随后可能与理论部分有点区别了，正常求解是直接用当前权重减去步长和梯度的乘积就是新的模型参数也就是${\theta }_j={\theta }_i-\alpha \nabla {\theta }_i$，但是呢UFLDL中提供了另一种快速解法，直接使用MATLAB的优化函数minFunc，关于此函数的解释可以戳这里 ,所以我们可以利用此函数对上面定义的损失函数的最小值，而且minFunc提供了几个很好地求梯度方法，这里使用大规模优化算法(LBFGS)代替上面人工计算的损失函数的梯度。

if ~exist('options', 'var')
    options = struct;
end
if ~isfield(options, 'maxIter')
    options.maxIter = 400;
end
% initialize parameters
theta = 0.005 * randn(numClasses * inputSize, 1);
% Use minFunc to minimize the function
% addpath minFunc/
options.Method = 'lbfgs'; % Here, we use L-BFGS to optimize our cost
                          % function. Generally, for minFunc to work, you
                          % need a function pointer with two outputs: the
                          % function value and the gradient. In our problem,
                          % softmaxCost.m satisfies this.
options.display = 'on';
% options.displayImage = false;
[softmaxOptTheta, cost] = minFunc( @(p) softmaxCost(p, ...
                                   numClasses, inputSize, lambda, ...
                                   inputData, labels), ...                     
                                   theta, options);
% Fold softmaxOptTheta into a nicer format
softmaxModel.optTheta = reshape(softmaxOptTheta, numClasses, inputSize);
softmaxModel.inputSize = inputSize;
softmaxModel.numClasses = numClasses;

这里优化的时候会自动优化(3)中输出项的第一个cost，而非第二项grad

这样就得到了模型参数了，可以进入测试阶段了

(5)测试阶段

先读数据

images = loadMNISTImages('t10k-images-idx3-ubyte');
labels = loadMNISTLabels('t10k-labels-idx1-ubyte');
labels(labels==0) = 10; % Remap 0 to 10
inputData = images;

预测标签

M = exp(theta*inputdata);
M = bsxfun(@rdivide, M, sum(M));
[temp pred] = max(M);

计算准确率

acc = mean(labels(:) == pred(:));
fprintf('Accuracy: %0.3f%%\n', acc * 100);

logistic回归和softmax关系

其实logistic回归就是具有两种类别的softmax, 直接由softmax推导看看(注意softmax具有参数冗余特性, 靠的是权重衰解决的):
$$\begin{array}{rl}& \frac{1}{e^{{\theta }_{1}x}+e^{{\theta }_{2}x}}\left[\begin{array}{c}{e}^{{\theta }_1\ast x}\\ e^{{\theta }_2\ast x}\end{array}\right]\\ =& \frac{1}{e^{0\ast x}+e^{({\theta }_2-{\theta }_1)x}}\left[\begin{array}{c}{e}^{0\ast x}\\ e^{({\theta }_2-{\theta }_1)\ast x}\end{array}\right]\\ =& \frac{1}{1+e^{{\theta }^{\prime }\ast x}}\left[\begin{array}{c}1\\ e^{{\theta }^{\prime }\ast x}\end{array}\right]\end{array}$$

结语

以上便是整个流程了，可运行代码下载戳此处

给各位大佬们提个建议哈，大部分代码不是能够一次性运行成功的，如果是计算机专业，强烈建议自行分析整个程序，嘿嘿，加油。若有错误，多多指正

本文已经同步到微信公众号中，公众号与本博客将持续同步更新运动捕捉、机器学习、深度学习、计算机视觉算法，敬请关注