在Python中使用MNIST较便捷的方法是安装MNIST工具包，在这里可以下载安装。

安装后便可如下直接导入训练数据和测试数据：

<code class="language-python">import matplotlib.pyplot as plt
import numpy as np
from mnist import MNIST
%matplotlib inline
mndata = MNIST('Path to your data/MNIST')
images, labels=mndata.load_training()
images_test, labels_test=mndata.load_testing()
</code>

计算Loss function，同样得解决指数函数数据溢出的问题：

<code class="language-python">def calculate_loss(w,x,y):
    theta=np.dot(np.transpose(w),x)
    ### first part
    l1=np.multiply(theta,y)
    ### second part
    thetaMax=np.max(theta,axis=0)
    theta-=thetaMax
    l2=np.log(np.sum(np.exp(theta),axis=0))
    ### sum together
    L=np.sum(-l1+l2+thetaMax)
    return L
</code>

计算导数：

<code class="language-python">def calculate_gradient(w,x_batch,y_batch):
    theta=np.dot(np.transpose(w),x_batch)
    theta-=np.max(theta,axis=0)
    mu=np.exp(theta)
    muSum=np.sum(mu,axis=0)
    mu=mu/muSum
    mu-=y_batch
    dL=np.dot(x_batch,np.transpose(mu))
    return dL
</code>

用SGD训练：

<code class="language-python">def train(x,y,batch_sz,lr,max_iter,loss_thresh):
    ### bias trick ###
    batch_sz=100
    d_len=x.shape[1]
    d_dim=x.shape[0]+1
    d_class=np.max(y)-np.min(y)+1
    x_b=np.concatenate((x,np.ones((1,d_len))),axis=0)
    yc=np.zeros((d_class,d_len))
    yc[y,np.arange(d_len)]=1;
    w=np.zeros((d_dim,d_class),dtype=np.float)
    
    Loss_old=0
    Loss=[]
    stepCnt=0
    ### Run SGD ###
    for iter in range(max_iter):
        ### sample a mini batch ###
        batch=np.arange(d_len)
        np.random.shuffle(batch)
        x_batch=x_b[:,batch[:batch_sz]]
        y_batch=yc[:,batch[:batch_sz]]
        ### update weight ###
        dL=calculate_gradient(w,x_batch,y_batch)
        w-=lr*dL
        ### record loss changes ###
        Loss.append(calculate_loss(w,x_b,yc))
        ### learning rate annealing ###
        stepCnt+=1
        if stepCnt==10:
            stepCnt=0
            lr*=0.8
        ### Check if converge ###
        if abs(Loss[-1]-Loss_old)<loss_thresh: break loss_old="Loss[-1]" return w,loss < code></loss_thresh:></code>

主函数：

<code class="language-python">x_train=np.transpose(np.matrix(images))
y_train=np.array(labels)
lr=0.5
batch_sz=50
Max_iter=200
loss_thresh=1e-3
w,Loss = train(x_train,y_train,batch_sz,lr,Max_iter,loss_thresh)
print w
plt.plot(Loss)
</code>

Loss function的收敛情况如下：

得到的权重W是一个 785x10 的矩阵，其中10表示有10个不同的类别。

单独拿出其中某一行矩阵，拿掉最后一列bias，是一个784维的向量\(w_c\)，即是对应于该行数字的预测权重。

这个权重在机器视觉上有何种意义呢？

下面变一个小魔术，如果你还记得我们的MNIST数据中的单个训练数据\(x_i\)，是将一张28x28的数字图片，转换成一个784维的向量得到的。于是我们这里反过来把784维的\(w_c\)转换成28x28的图片，得到了以下的10张图片：

分别对应0～9这10个数字，从中依稀可以辨别出来。所以，这里分类字符的原理，其实可以看作是模板匹配（template matching），训练过程是为了从训练集中综合出最合适的数字模板。

简要对测试集做了一下测试，以下补上prediction部分的代码：

<code class="language-python">def predict(w,x):
    score=np.dot(np.transpose(w),x)
    pred=np.argmax(score)
    return pred
</code>

随机抽取40个测试数据做了预测，这里仅定性说明softmax分类器训练结果的有效性，更加具体的测试数据以后更新。令人excited的一点是，虽然这个实验中我们只是用最简单的思路做了一些微小的工作，不加处理的对像素数据进行了regression，准确率已经比较可观了，这很大程度上可以解释为mnist数字，从分层理解的模型中来看，层次是很浅的，大约就是 像素->边缘->基本几何形状 这三层。

引用 novakon:
这个例子所用的算法，很接近传统OCR。但由于MNIST都是手写字体，字形变化很大，所以直接用784维向量而不加处理（比如卷积！）的话，loss func 的收敛情况其实不佳（从模板中8与9的相像程度来……