OC实现Softmax识别手写数字

简介

Softmax回归模型是logistic回归模型在多分类问题上的推广，在多分类问题中，类标签 y 可以取两个以上的值。Softmax模型运用广泛，很多复杂精细的训练模型最后一步都会用softmax来分配概率。

公式

Softmax回归最主要的代价函数如下：

代价函数

其中，1｛-｝是示性函数，其取值规则为：1｛值为真的表达式｝＝ 1。
根据代价函数，分别对k个分类进行求导，就可以得到：

k分类偏导函数

其中，Xi为类别j的概率：

Xi为类别j的概率

有了上面的偏导数公式以后，我们就可以用梯度下降算法更新theta值，每一次迭代需要进行如下更新:

j = 1,...,k

以上公式并没有涉及到bias优化，加上bias的公式形如：

softmax-regression-vectorequation.png

其实bias的迭代更新和theta一样，只需要找到bias的偏导数就行，我的代码中包含了对bias求导优化。

数据

训练所用数据集是MNIST，MNIST数据集的官网是Yann LeCun's website。这个数据集包含60000行的训练数据集和10000行的测试数据集。

文件	内容
train-images-idx3-ubyte.gz	训练集图片 - 60000 张训练图片
train-labels-idx1-ubyte.gz	训练集图片对应的数字标签
t10k-images-idx3-ubyte.gz	测试集图片 - 10000 张图片
t10k-labels-idx1-ubyte.gz	测试集图片对应的数字标签

其中每一张图片都是28*28的，矩阵表示如下：

MNIST-Matrix.png

在数据集中，每张图片都是一个展开的向量，长度是 28x28 = 784。因此，在MNIST训练数据集中，train-images-idx3-ubyte解析后是一个形状为 [60000, 784]的张量。

mnist-train-xs.png

而train-labels-idx1-ubyte则是一个长度60000的向量，每一个值对应train-images-idx3-ubyte中代表的数字范围0～9。
测试数据与训练数据结构一样，只是数量只有10000张。

下载过后读取与解析到数据代码：

//
//  MLLoadMNIST.m
//  MNIST
//
//  Created by Jiao Liu on 9/23/16.
//  Copyright © 2016 ChangHong. All rights reserved.
//

#import "MLLoadMNIST.h"

@implementation MLLoadMNIST

int reverseInt(int input)
{
    unsigned char ch1, ch2, ch3, ch4;
    ch1=input&255;
    ch2=(input>>8)&255;
    ch3=(input>>16)&255;
    ch4=(input>>24)&255;
    return((int)ch1<<24)+((int)ch2<<16)+((int)ch3<<8)+ch4;
}

double **readImageData(const char *filePath)
{
    FILE *file = fopen(filePath, "rb");
    double **output = NULL;
    if (file) {
        int magic_number=0;
        int number_of_images=0;
        int n_rows=0;
        int n_cols=0;
        fread((char*)&magic_number, sizeof(magic_number), 1, file);
        magic_number= reverseInt(magic_number);
        fread((char*)&number_of_images, sizeof(number_of_images), 1, file);
        number_of_images= reverseInt(number_of_images);
        fread((char*)&n_rows, sizeof(n_rows), 1, file);
        n_rows= reverseInt(n_rows);
        fread((char*)&n_cols, sizeof(n_cols), 1, file);
        n_cols= reverseInt(n_cols);
        output = (double **)malloc(sizeof(double) * number_of_images);
        for(int i=0;i<number_of_images;++i)
        {
            output[i] = (double *)malloc(sizeof(double) * n_rows * n_cols);
            for(int r=0;r<n_rows;++r)
            {
                for(int c=0;c<n_cols;++c)
                {
                    unsigned char temp=0;
                    fread((char*)&temp, sizeof(temp), 1, file);
                    output[i][(n_rows*r)+c]= (double)temp;
                }
            }
        }
    }
    fclose(file);
    return output;
}

int *readLabelData(const char *filePath)
{
    FILE *file = fopen(filePath, "rb");
    int *output = NULL;
    if (file) {
        int magic_number=0;
        int number_of_items=0;
        fread((char*)&magic_number, sizeof(magic_number), 1, file);
        magic_number= reverseInt(magic_number);
        fread((char*)&number_of_items, sizeof(number_of_items), 1, file);
        number_of_items= reverseInt(number_of_items);
        output = (int *)malloc(sizeof(int) * number_of_items);
        for(int i=0;i<number_of_items;++i)
        {
            unsigned char temp=0;
            fread((char*)&temp, sizeof(temp), 1, file);
            output[i]= (int)temp;
        }
    }
    fclose(file);
    return output;
}

Softmax实现

首先我这里选用的梯度下降算法去逼近最值，所以需要一个迭代次数，代码中默认的是500次。输入训练图片是60000，如果每次迭代都全部用上，训练会花去很多时间，所以每次迭代默认随机取100张图片进行训练。
下降梯度的速率默认是0.01。
代码中实现两种梯度下降逼近，一种是每次迭代依次使用每张图片去更新所有分类变量，另一种是每次迭代顺序更新每个分类变量，更新每个分类时候使用所有随机图片数据。两种方法其实效率一样，但是测试时候发现第二种方法的正确率略高。

代码如下：

//
//  MLSoftMax.m
//  MNIST
//
//  Created by Jiao Liu on 9/26/16.
//  Copyright © 2016 ChangHong. All rights reserved.
//

#import "MLSoftMax.h"

@implementation MLSoftMax

- (id)initWithLoopNum:(int)loopNum dim:(int)dim type:(int)type size:(int)size descentRate:(double)rate
{
    self = [super init];
    if (self) {
        _iterNum = loopNum == 0 ? 500 : loopNum;
        _dim = dim;
        _kType = type;
        _randSize = size == 0 ? 100 : size;
        _bias = malloc(sizeof(double) * type);
        _theta = malloc(sizeof(double) * type * dim);
        for (int i = 0; i < type; i++) {
            _bias[i] = 0;
            for (int j = 0; j < dim; j++) {
                _theta[i * dim +j] = 0.0f;
            }
        }
        
        _descentRate = rate == 0 ? 0.01 : rate;
    }
    return  self;
}

- (void)dealloc
{
    if (_bias != NULL) {
        free(_bias);
        _bias = NULL;
    }
    
    if (_theta != NULL) {
        free(_theta);
        _theta = NULL;
    }
    
    if (_randomX != NULL) {
        free(_randomX);
        _randomX = NULL;
    }
    
    if (_randomY != NULL) {
        free(_randomY);
        _randomY = NULL;
    }
}

#pragma mark - SoftMax Main

- (void)randomPick:(int)maxSize
{
    long rNum = random();
    for (int i = 0; i < _randSize; i++) {
        _randomX[i] = _image[(rNum+i) % maxSize];
        _randomY[i] = _label[(rNum+i) % maxSize];
    }
}
/*
- (double *)MaxPro:(double *)index
{
    long double maxNum = index[0];
    for (int i = 1; i < _kType; i++) {
        maxNum = MAX(maxNum, index[i]);
    }
    
    long double sum = 0;
    for (int i = 0; i < _kType; i++) {
        index[i] -= maxNum;
        index[i] = expl(index[i]);
        sum += index[i];
    }
    
    for (int i = 0; i < _kType; i++) {
        index[i] /= sum;
    }
    return index;
}

- (void)updateModel:(double *)index currentPos:(int)pos
{
    double *delta = malloc(sizeof(double) * _kType);
    for (int i = 0; i < _kType; i++) {
        if (i != _randomY[pos]) {
            delta[i] = 0.0 - index[i];
        }
        else
        {
            delta[i] = 1.0 - index[i];
        }
        
        _bias[i] -= _descentRate * delta[i];
        
        for (int j = 0; j < _dim; j++) {
            _theta[i * _dim +j] += _descentRate * delta[i] * _randomX[pos][j] / _randSize;
        }
    }
    
    if (delta != NULL) {
        free(delta);
        delta = NULL;
    }
}

- (void)train
{
    _randomX = malloc(sizeof(double) * _randSize);
    _randomY = malloc(sizeof(int) * _randSize);
    double *index = malloc(sizeof(double) * _kType);
    
    for (int i = 0; i < _iterNum; i++) {
        [self randomPick:_trainNum];
        for (int j = 0; j < _randSize; j++) {
            // calculate wx+b
            vDSP_mmulD(_theta, 1, _randomX[j], 1, index, 1, _kType, 1, _dim);
            vDSP_vaddD(index, 1, _bias, 1, index, 1, _kType);
            
            index = [self MaxPro:index];
            [self updateModel:index currentPos:j];
        }
    }
    if (index != NULL) {
        free(index);
        index = NULL;
    }
}
*/

- (int)indicator:(int)label var:(int)x
{
    if (label == x) {
        return 1;
    }
    return 0;
}

- (double)sigmod:(int)type index:(int) index
{
    double up = 0;
    for (int i = 0; i < _dim; i++) {
        up += _theta[type * _dim + i] * _randomX[index][i];
    }
    up += _bias[type];
    
    double *down = malloc(sizeof(double) * _kType);
    double maxNum = -0xfffffff;
    vDSP_mmulD(_theta, 1, _randomX[index], 1, down, 1, _kType, 1, _dim);
    vDSP_vaddD(down, 1, _bias, 1, down, 1, _kType);
    
    for (int i = 0; i < _kType; i++) {
        maxNum = MAX(maxNum, down[i]);
    }
    
    double sum = 0;
    for (int i = 0; i < _kType; i++) {
        down[i] -= maxNum;
        sum += exp(down[i]);
    }
    
    if (down != NULL) {
        free(down);
        down = NULL;
    }
    
    return exp(up - maxNum) / sum;
}

- (double *)fderivative:(int)type
{
    double *outP = malloc(sizeof(double) * _dim);
    for (int i = 0; i < _dim; i++) {
        outP[i] = 0;
    }
    
    double *inner = malloc(sizeof(double) * _dim);
    for (int i = 0; i < _randSize; i++) {
        long double sig = [self sigmod:type index:i];
        int ind = [self indicator:_randomY[i] var:type];
        double loss = -_descentRate * (ind - sig) / _randSize;
        _bias[type] += loss * _randSize;
        vDSP_vsmulD(_randomX[i], 1, &loss, inner, 1, _dim);
        vDSP_vaddD(outP, 1, inner, 1, outP, 1, _dim);
    }
    if (inner != NULL) {
        free(inner);
        inner = NULL;
    }
    
    return outP;
}

- (void)train
{
    _randomX = malloc(sizeof(double) * _randSize);
    _randomY = malloc(sizeof(int) * _randSize);
    for (int i = 0; i < _iterNum; i++) {
        [self randomPick:_trainNum];
        for (int j = 0; j < _kType; j++) {
            double *newTheta = [self fderivative:j];
            for (int m = 0; m < _dim; m++) {
                _theta[j * _dim + m] = _theta[j * _dim + m] - _descentRate * newTheta[m];
            }
            if (newTheta != NULL) {
                free(newTheta);
                newTheta = NULL;
            }
        }
    }
}

- (void)saveTrainDataToDisk
{
    NSFileManager *fileManager = [NSFileManager defaultManager];
    NSString *thetaPath = [[NSSearchPathForDirectoriesInDomains(NSCachesDirectory, NSUserDomainMask, YES) objectAtIndex:0] stringByAppendingString:@"/Theta.txt"];
//    NSLog(@"%@",thetaPath);
    NSData *data = [NSData dataWithBytes:_theta length:sizeof(double) *  _dim * _kType];
    [fileManager createFileAtPath:thetaPath contents:data attributes:nil];
    
    NSString *biasPath = [[NSSearchPathForDirectoriesInDomains(NSCachesDirectory, NSUserDomainMask, YES) objectAtIndex:0] stringByAppendingString:@"/bias.txt"];
    data = [NSData dataWithBytes:_bias length:sizeof(double) * _kType];
    [fileManager createFileAtPath:biasPath contents:data attributes:nil];
}

- (int)predict:(double *)image
{
    double maxNum = -0xffffff;
    int label = -1;
    double *index = malloc(sizeof(double) * _kType);
    vDSP_mmulD(_theta, 1, image, 1, index, 1, _kType, 1, _dim);
    vDSP_vaddD(index, 1, _bias, 1, index, 1, _kType);
    for (int i = 0; i < _kType; i++) {
        if (index[i] > maxNum) {
            maxNum = index[i];
            label = i;
        }
    }
    return label;
}

- (int)predict:(double *)image withOldTheta:(double *)theta andBias:(double *)bias
{
    double maxNum = -0xffffff;
    int label = -1;
    double *index = malloc(sizeof(double) * _kType);
    vDSP_mmulD(theta, 1, image, 1, index, 1, _kType, 1, _dim);
    vDSP_vaddD(index, 1, bias, 1, index, 1, _kType);
    for (int i = 0; i < _kType; i++) {
        if (index[i] > maxNum) {
            maxNum = index[i];
            label = i;
        }
    }
    return label;
}

@end

最后训练结果：

Simulator Screen Shot

不断改变循环次数，下降速率等参数，都会带来正确率的变化。测试结果发现默认参数能带来接近最好的识别的正确率，90%左右👏。

结语

90％的正确率并非达到最优，因为这仅仅是个简单的模型，可以加上卷积神经网络来改善效果。
关于卷积神经网络实现，我也实现了一版，但由于效率与正确率不高，后面优化后再分享😊。

最后编辑于：2017.12.04 04:24:56

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 194,390评论 5赞 459
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 81,821评论 2赞 371
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 141,632评论 0赞 319
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 52,170评论 1赞 263
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 61,033评论 4赞 355
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 46,098评论 1赞 272
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 36,511评论 3赞 381
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 35,204评论 0赞 253
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 39,479评论 1赞 290
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 34,572评论 2赞 309
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 36,341评论 1赞 326
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 32,213评论 3赞 312
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 37,576评论 3赞 298
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 28,893评论 0赞 17
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 30,171评论 1赞 250
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 41,486评论 2赞 341
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 40,676评论 2赞 335

OC实现Softmax识别手写数字

简介

公式

数据

Softmax实现

结语

推荐阅读更多精彩内容