一步一步学习Deep Learning (3) — Softmax Regression (C++)

前言

关于SR理论公式以及求导过程可以参考本博客另外一篇博文 一步一步学习Deep Learning (2) — Softmax Regression (Matlab)。本博文主要是利用C++来实现SR,并在MNIST数据库上做测试(关于MNIST数据库的loading问题,可以参看博文MNIST loading)。


实现

首先利用SR的Matlab版在MNIST数据库上做测。这里有一个实现上的trick可以参考,其中一个是梯度检测,在实现优化的过程中,人为的去debug这些梯度部分会比较麻烦,梯度检测旨在利用有限梯度\(\frac{f(x+delta(x)) + f(x-delta(x))}{2*delta(x)}\)去检查自己的实现有没有问题,这里给出别人实现的一个代码:

function average_error = grad_check(fun, theta0, num_checks, varargin)

delta = 1e-3; 
sum_error = 0;

fprintf(' Iter       i             err');
fprintf('           g_est               g               f\n')

for i=1:num_checks
    T = theta0;
    j = randsample(numel(T), 1);
    T0=T; T0(j) = T0(j) - delta;
    T1=T; T1(j) = T1(j) + delta;

    [f, g] = fun(T, varargin{:});
    f0 = fun(T0, varargin{:});
    f1 = fun(T1, varargin{:});

    g_est = (f1 - f0) / (2 * delta);
    error = abs(g(j) - g_est);

    fprintf('% 5d  % 6d % 15g % 15f % 15f % 15f\n', ...
            i, j, error, g(j), g_est, f);

    sum_error = sum_error + error;
end

average=sum_error/num_checks;

当然,在运行的时候记得将梯度检测以及任何后续可能会涉及的检测全部关掉,这些检测可是非常耗时间的。在检测无误后,我们可以进行MNIST数据的测试了,这里选择的优化算法为LBFGS,优化工具箱为minFunc结果如下图所示:

softmax

根据UFLDL,结果只要不低于91%,且不高于99%的话,应该都算正常。我们来看看C++的实现是否正确。

C++实现

本次实现主要是为了方便理解SR的一些基本原理,包括练习一下C++语言,所以并没有考虑太多的整体架构,可能使用起来不是特别方便,日后改进。

首先,我们来定义一个Softmax Regression的类,如下:

class SMaxReg
{
public:
    SMaxReg() {}
    SMaxReg(const int epochs, const int numClasses, const double learningRate,
            const double moment, const double regularizer, const double epsilon);
    ~SMaxReg();

    // setter and getter
    void setInitWeight(const Mat &weight);

    Mat getWeight() const;

    // train
    void train(const Mat &data, const Mat &label, int batchSize);

    // prediction
    void predict(const Mat &data, const Mat &label, Mat &predictLabel, 
                 Mat &prob, double &accuracy);

protected:
    void realLabel2Matrix(const Mat &label, Mat &biLabel);
    void randperm(vector<int> &index);
    void getBatchData(const Mat &data, const Mat &label, const vector<int> &index,
                      const int batchSize, const int batchIdx, 
                      Mat &batchData, Mat &batchLabel);

private:
    Mat weight;      
    int epochs;      
    int numClasses;  
    double epsilon;  
    double learningRate;
    double moment;
    double regularizer;

private:
    SMaxReg(const SMaxReg &rhs); // do not allow copy constructor
    const SMaxReg &operator = (const SMaxReg &); // nor assignement operator
};

C++的实现中,我们主要采用Stochastic Gradient Descent (SGD)算法来进行优化,主要过程为:

1. 设定迭代次数或者周期数,并在每一次的周期过程中随机化数据的顺序
2. 每次读取一定量的数据batch,可以是100个,当然最朴素的随机情况是只读取一个
3. 利用这一个batch进行学习更新现有的权值,并计算出目标函数值
4. 重复2,3直到一个周期结束
5. 重复1,2,3直到达到停止条件

关于SGD部分的深入理解,会在另外一篇博客中单独讨论。这里仅仅将其作为一个工具来使用。基于上述过程的指导,先简单解释private某些变量部分,其中weight就是SR分类器的权值,是需要学习得到的东西,epochs是最大的周期数,一周期数等于一次扫描过程,也就是N次迭代,这里的N为数据的数量。epsilon是其中的一个停止条件,它的含义是如果两次迭代的目标函数值低于这个设定值的话,则认为优化收敛,即可以停止优化。learningRate是梯度下降中的一个学习率。关于moment的技术,可以参考UFLDLMomentum的解释部分。regularizer主要是控制权值weight在学习的过程中,避免出现特别大或者特别小的值,以减少过拟合现象。接下来我们一一来解释几个重要的函数。

void SMaxReg::realLabel2Matrix(const Mat &label, Mat &labelMatrix)
{
    assert(label.type() == CV_32FC1);

    labelMatrix = Mat::zeros(label.cols, numClasses, label.type());
    float *lmPtr = (float *)labelMatrix.data;
    float *lPtr = (float *)label.data;
    for (int i = 0; i < label.cols; i++) {
        lmPtr[i*labelMatrix.cols + (uchar)(*lPtr++)] = 1;
    }

    lmPtr = NULL;
    lPtr = NULL;
}

这个函数的目的是将简单的数值label,例如\((1, 2, 2, 1, 2, 4, 3, 3, ...)\)转换成对应的二值化的label情况,例如,本MNIST中,类别数为\(10\),则\(1 = 1000000000\)\(2 = 0100000000\) 等等。这样做的目的是为了利用OpenCV的矩阵操作,并且同时计算多类的梯度,以及利用多类梯度进行权值的同时更新过程。当然,如果是针对二分类的情况,也可以直接将某一类看成正类,其余的看成负类,依次交替,训练多次。

void SMaxReg::randperm(vector<int> &index)
{
    srand((uchar)time(NULL));
    std::random_shuffle(index.begin(), index.end());
}

void SMaxReg::getBatchData(const Mat &data, const Mat &label, const vector<int> &index,
                           const int batchSize, const int batchIdx, 
                           Mat &batchData, Mat &batchLabel)
{
    int batchStart = batchIdx * batchSize;
    int batchEnd = min((batchIdx + 1)*batchSize, data.rows);
    for (int i = batchStart; i < batchEnd; ++i) {
        data.row(index[i]).copyTo(batchData.row(i - batchStart));
        label.row(index[i]).copyTo(batchLabel.row(i - batchStart));
    }
}

randperm 是将数据的顺序打乱,即随机梯度过程的第一个过程。而getBatchData则是得到数据的其中一个小的batch。

void SMaxReg::train(const Mat &data, const Mat &label, int batchSize)
{
    assert(data.type() == CV_32FC1 && label.type() == CV_32FC1);
    assert(batchSize >= 1);

    srand((uchar)time(NULL));

    Mat labelMatrix;
    realLabel2Matrix(label, labelMatrix);

    // initialize weight and bias
    int numData = data.rows;
    int dimData = data.cols;
    if (weight.empty() == true) {
        weight = Mat::zeros(dimData, numClasses, data.type());
        randu(weight, 0, 1); weight *= 0.001;
    }
    else {
        if (weight.rows != dimData || weight.cols != numClasses) {
            printf("Initial weight dimension wrong: weight.rows == dimData && weight.cols == numClasses!\n");
            abort();
        }
    }
    Mat velocity = Mat::zeros(weight.size(), weight.type());

    vector<int> index(numData, 0);
    for (int i = 0; i < numData; ++i) {
        index[i] = i;
    }

    Mat batchData, batchLabel; // batch data and label
    Mat rsp, maxRsp; // feed-forward response
    Mat prob, sumProb, logProb; // softmax prediction probability
    Mat gradient, ww; // update

    bool isConverge = false; 
    double prevCost = 0;
    double currCost = 0; // cost
    double mom = 0.5; // moment
    int t = 0; // iteration
    int numBatches = floor(numData / batchSize);
    for (int ei = 0; ei < epochs; ++ei) {
        randperm(index);
        batchData = Mat::zeros(batchSize, dimData, data.type());
        batchLabel = Mat::zeros(batchSize, numClasses, labelMatrix.type());
        for (int batchIdx = 0; batchIdx < numBatches; ++batchIdx) {
            // batch SGD
            t++;
            if (t == 30)
                mom = moment;

            getBatchData(data, labelMatrix, index, batchSize, batchIdx, batchData, batchLabel);


            // compute softmax response
            rsp = batchData * weight;
            reduce(rsp, maxRsp, 1, REDUCE_MAX, rsp.type());
            rsp -= repeat(maxRsp, 1, numClasses);
            exp(rsp, prob);
            reduce(prob, sumProb, 1, REDUCE_SUM, prob.type());
            prob = prob / repeat(sumProb, 1, numClasses);

            // compute gradient
            gradient = batchLabel - prob;
            gradient = batchData.t() * gradient;
            gradient = -gradient / batchData.rows;
            gradient += regularizer * weight;

            // update weight and bias
            velocity = mom * velocity + learningRate * gradient;
            weight -= velocity;

            // compute objective cost
            log(prob, logProb);
            logProb = batchLabel.mul(logProb);
            currCost = -(sum(logProb)[0] / batchData.rows);
            pow(weight, 2, ww);
            currCost += sum(ww)[0] * 0.5 * regularizer;
            printf("epoch %d: processing batch %d / %d cost %f\n", ei + 1, batchIdx + 1, numBatches, currCost);

            if (abs(currCost - prevCost) < epsilon && ei != 0) {
                printf("objective cost variation less than pre-defined %f\n", epsilon);
                isConverge = true;
                break;
            }
            prevCost = currCost;
        }

        batchData.release();
        batchLabel.release();
        if (isConverge == true) break;
    }

    if (isConverge == false) {
        printf("stopped by reaching maximum number of iterations\n");
    }

    // free and destroy space
    index.clear();
    labelMatrix.release();
    batchData.release(); 
    batchLabel.release(); 
    rsp.release();
    maxRsp.release(); 
    prob.release(); 
    sumProb.release();
    logProb.release();
    gradient.release();
    ww.release();
    velocity.release();
}

这是完整的SR的训练函数,重点看三个部分,第一个是计算SR的响应部分:

// compute softmax response
rsp = batchData * weight;
reduce(rsp, maxRsp, 1, REDUCE_MAX, rsp.type());
rsp -= repeat(maxRsp, 1, numClasses);
exp(rsp, prob);
reduce(prob, sumProb, 1, REDUCE_SUM, prob.type());
prob = prob / repeat(sumProb, 1, numClasses);

这里利用OpenCVreduce以及repeat函数来进行矩阵操作,免得自己写各种for循环。repeat的目的比较简单,即将maxRsp按照[rows, cols]来复制。例如maxRsp\(1000\) x \(1\)的矩阵,numClasses = 10,则repeat(maxRsp, 1, numClasses)之后,输出为\(1000\) x \(10\)的矩阵。reduce的目的是按照不同的准则规约,例如REDUCE_MAX则是求出每一行或者列的最大值。具体函数的应用可以参考reduce。其实不停的运用repeat以及reduce的目的是为了利用高效的OpenCV矩阵计算,比自己写出来的可能要快不少。接下来我们来看看梯度的计算部分:

// compute gradient
gradient = batchLabel - prob;
gradient = batchData.t() * gradient;
gradient = -gradient / batchData.rows;
gradient += regularizer * weight;

因为OpenCV对于类Mat的操作做了各种优化,其中包括各种运算符的重载,使得这一部分的代码与Matlab的实现没有什么区别。接下来是梯度的更新部分:

// update weight and bias

// original
thisLr = learningRate / (1 + regularizer * learningRate * t);
weight -= thisLr * gradient;

// momentum
velocity = mom * velocity + learningRate * gradient;
weight -= velocity;

// after one epochs
learningRate = 0.5 * learningRate;

这里给出了两种梯度更新的准则,其中original是最基本的更新准则,直接根据当前的梯度进行更新。而关于momentum的更新模式是在最近的deep learning的框架中常用的一种方法,据说针对大数据能够比较有效快速的收敛。在MNIST数据库上,两种更新模式都测试过,感觉无论是收敛速度还是准确度都差不多。原因可能是60000的训练数据比较小吧。整个训练函数最核心的三部分解释完了,接下来看看预测部分:

void SMaxReg::predict(const Mat &data, const Mat &label, Mat &predictLabel, 
                      Mat &prob, double &accuracy)
{
    prob = data * weight;
    predictLabel = Mat::zeros(1, data.rows, CV_32FC1);
    int *maxLoc = (int *)calloc(2, sizeof(int));
    float *pPtr = (float *)predictLabel.data;
    float *lPtr = (float *)label.data;
    accuracy = 0;
    for (int i = 0; i < data.rows; ++i) {
        minMaxIdx(prob.row(i), NULL, NULL, NULL, maxLoc);
        *pPtr++ = maxLoc[1];
        if ((*lPtr++) == maxLoc[1])
            accuracy++;
    }
    accuracy /= data.rows;

    if (maxLoc != NULL) free(maxLoc); 
    maxLoc = NULL;
    lPtr = NULL;
}

预测就比较简单了,可以直接参考原始的Matlab的实现,给出来的目的是自己做一个记录。


MNIST测试

有了上面的实现,接下来进行MNIST数据的测试,完整的代码包括测试的参数请参看SoftMax Regression。这里给出一个结果截图:

softmax

结果在可接受的范围,由于是优化过程中有随机存在,所以进行多次运行测试,发现结果都在可接受的范围类。这里有几个经验需要记住:

1. 随机梯度需要有足够多的迭代次数才能够接近最优解,所以对于本例60000个训练数据,个人设定的batchSize为100~1000都可以

2. 学习率也是影响结果的一个重要因素,当然简单一点的想法是,设定足够多的训练迭代次数,学习率尽可能的低一点

3. 随机梯度进行学习的过程会出现不断的摇摆,但可以通过观察目标函数的整体趋势判断是否在正确的道路上

4. 随机梯度达到合理预测值的训练速度非常的快,但是会经常性的在最优值附近摆动,因此可能与最优的结果有一定程度上的偏差

5. 随机梯度下降的学习率衰减非常的重要,经过一些测试,如果设定学习率衰减的比较快的话,意味这再后面的迭代更新过程中,当前梯度带来的影响变小,类似曲线曲率在最优解附近降低的过程。具体来说,如果将 `learningRate = factor * learningRate`中,factor设定为0.1的话,得到的结果基本上在92%附近。

后记

利用C++实现机器学习算法最重要的是能够合理的解出优化过程,然后可以利用各种现有的工具包如OpenCV以及armadillo来进行快速的矩阵操作,其训练速度不比Matlab矩阵化后的慢。针对大数据下的机器学习模型,可能C++更合适一点。例如深度学习工具包cxxnet, caffe 以及 MatConv里面核心的算法都是c++来实现。所以,日后可能会更多的自己实现一些机器学习的算法吧!

Comments