梯度下降C++实现

梯度下降是机器学习的基础内容，而机器学习又是人工智能这个领域非常基础的内容。作为一个自认为已经入门机器学习的小菜鸡，就来讲述一下已经可以熟练应用的梯度下降该怎么用C++来实现。

首先需要描述一下梯度下降一般会应用在什么情景中。

比如说，我们现在想要求一个w = f(x, y, z)的这个f函数。经过理论推导，我们已经知道了f(x, y) = ax + by + cz，但是参数a和b是由实际材料或者其他因素决定的参数。这个时候，如果我们可以通过实际测量知道每一个(x, y, z)对应的w，我们就可以采用梯度下降的办法来拟合出a,b和c的近似值。当然，对于上述这种简单的线性函数表达式，我们远远不必使用梯度下降这么高大上的方法。但是，当函数的表达式非常非常复杂，或者说，我们根本无法写出这个函数的表达式的时候，我们就需要采用梯度下降的办法了。只要我们可以通过给定预测的参数值算出这个预测参数值对应的函数结果，就可以使用梯度下降。很多复杂的函数往往可以用程序通过传入参数而给出返回值，却无法直接写成一个表达式，这就是之前说的写不出表达式的情况。现在，为了介绍原理，我们就使用这个简单的ax + by + cz函数。想要拟合出a,b和c，我们就需要有很多很多组(x, y, z, (ax + by +cz))这样的数据。例程里，就让x和y和z从(0, 1)这个范围里均匀取值了。

假设a,b和c的真实值是0.8,0.2和0.3。我们最开始预测的a,b和c的值都为0.5，然后一步一步向真实值靠近。

还需要介绍的是loss这个概念。因为我们已经知道了每个(x, y, z)对应的w，所以我们可以计算当前预测参数值和真实参数值之间的差距大小。我们将预测参数值带入f函数，对每个(x, y, z)求出预测参数对应的w，这里记作w'。我们先定义一种比较简单的loss function，就是Σ(w-w')*(w-w')。可想而知，当预测参数跟实际参数之间的差距越大时，这个w和w的差距一般情况下就会越大。所以，我们的目标就是寻找loss的最小值。

真正应用梯度下降的时候，我们是需要事先获得大量数据的，这里就直接用程序生成这些数据了

vector<float> func(vector<float>& para) {
	vector<float> ret;
	for (float x = 0.f; x < 1.f; x += .1f) {
		for (float y = 0.f; y < 1.f; y += .1f) {
			for (float z = 0.f; z < 1.f; z += .1f) {
				ret.push_back(para[0] * x + para[1] * y + para[2] * z);
			}
		}
	}
	return ret;
}

vector<float> p;
p.push_back(.8f);
p.push_back(.2f);
p.push_back(.3f);
vector<float> s = func(p);

这里传入的para参数是实际真实值，等到进入梯度下降的过程的时候，就需要传入预测值了，因为在解决问题的时候我们是不知道最终结果的。

这样我们就模拟好了一个完整的问题情境。重新描述一下现在这个问题：已知f(x, y, z) = ax + by + cz，并且知道很多组(x, y, z, f(x, y, z))的数据，要求出a, b, c的近似值。思路就如上文所说，求出loss的最小值。

我们这里是需要找loss(x, y, z)的最小值，要阐述为什么梯度下降可以找到答案，先用一维函数y = f(x)举例。假定x等于0.4的时候，y取到最小值。那么现在预测x为0.5，在0.5出求y对x的导数，很大概率这个点的导数是正数，也就是说dy/dx>0，因为只有这样，当x从0.5减小到0.4的时候，y会减小到最小值。正是应用这个原理，我们求loss对所有变量的导数，然后沿着导数为负的方向移动一段距离。在大部分情况下，每次移动之后，loss是会比前一次的loss值要小。形象一点地解释，就比如一个大坑，我们在坑的边上，要到坑的最底部。我们可以仅考虑脚下的坡度，往坡下走一小步。每次都往坡下走一小步，总能走到最低点。

原理解释完毕，于是可以阅读一下最主要的伪代码，就是梯度下降部分。

while (1) {
	tmpLoss = loss(para)
	for each para {
		dp[i] = d(loss)/d(para[i])
	}
	normalize(dp)
	for each para {
		para[i] -= dp[i]
	}
}

伪代码可以很明确地表示，我们沿着使loss减小的方向移动para这个向量(我们可以把数组看作多维向量)。不过这个while(1)什么时候退出呢？大部分情况是只要我们的loss减小到可接受的程度，我们就可以break了。

基本原理介绍完毕，有关梯度下降的优化以后再继续补充。

下附完整代码：

gradient.h

#pragma once
#include <vector>
#include <iostream>

#define RATIO .8f

using std::cout;
using std::endl;
using std::vector;

float variant(vector<float> &std, vector<float> &cal);
float clamp(float floor, float ceil, float d);

class Grad {
public:
	Grad(int paraNum, float initStep, float funcDiff, vector<float>(*predFunc)(vector<float> ¶),
		bool debug, bool randGen = true, int maxRound = INT_MAX,
		float(*lossFunc)(vector<float> &, vector<float> &) = variant,
		float(*miniStep)(float var, float step) = NULL, bool(*breakCond)(float var) = NULL);
	
	void prepare();
	void standard(vector<float> &std);
	long desc();
	void print();

private:
	bool randGen;
	bool debug;

	int paraNum;
	float funcDiff;
	float descStep;
	int tmpRound, maxRound;
	float tmpLoss, minLoss = INFINITY;

	vector<float> tmpPara, resPara;
	vector<float> tmpDiff, oldDiff;
	vector<float> tmpVel, oldVel;
	vector<int> mostShrink;

	vector<float> stdOpt, calOpt;
	vector<float> (*predFunc)(vector<float> &);
	float(*lossFunc)(vector<float> &, vector<float> &);
	float(*miniStep)(float var, float step);
	bool(*breakCond)(float var);
};

gradient.cpp

#include "gradient.h"
#include "time.h"

Grad::Grad(int paraNum, float initStep, float funcDiff, vector<float>(*predFunc)(vector<float> ¶), bool debug, bool randGen, int maxRound,
	float(*lossFunc)(vector<float> &, vector<float> &), float(*miniStep)(float var, float step), bool(*breakCond)(float var)) {
	this->paraNum = paraNum;
	this->descStep = initStep;
	this->funcDiff = funcDiff;
	this->predFunc = predFunc;
	this->randGen = randGen;
	this->debug = debug;
	this->maxRound = maxRound;
	this->lossFunc = lossFunc;
	this->miniStep = miniStep;
	this->breakCond = breakCond;

	tmpPara = vector<float>(paraNum);
	resPara = vector<float>(paraNum);
	tmpDiff = vector<float>(paraNum);
	oldDiff = vector<float>(paraNum);
	tmpVel = vector<float>(paraNum);
	oldVel = vector<float>(paraNum);
	mostShrink = vector<int>(paraNum);
}

void Grad::prepare() {
	if (randGen)
		for (float &i : tmpPara)i = float(rand() % 11) / 10;
	else
		for (float &i : tmpPara)i = .5f;

	for (int i = 0; i < paraNum; i++) {
		tmpDiff[i] = oldDiff[i] = tmpVel[i] = oldVel[i] = resPara[i] = .0f;
		mostShrink[i] = 0;
	}

	tmpRound = 0;
}
void Grad::standard(vector<float> &std) {
	stdOpt = std;
}
long Grad::desc() {
	clock_t t = clock();

	while (1) {
		if(debug)print();
		calOpt = predFunc(tmpPara);
		tmpLoss = lossFunc(stdOpt, calOpt);
		if (tmpLoss < minLoss) {
			minLoss = tmpLoss;
			resPara = tmpPara;
			tmpRound = 0;
		}
		else {
			tmpRound++;
			if (tmpRound > maxRound)break;
		}

		for (int i = 0; i < paraNum; i++) {
			tmpPara[i] += funcDiff;
			calOpt = predFunc(tmpPara);
			tmpDiff[i] = (lossFunc(stdOpt, calOpt) - tmpLoss) / funcDiff;
			tmpPara[i] -= funcDiff;
			if (tmpPara[i] == 0.f&&tmpDiff[i] > 0)tmpDiff[i] = 0.f;
			if (tmpPara[i] == 1.f&&tmpDiff[i] < 0)tmpDiff[i] = 0.f;
			tmpDiff[i] /= 1 << (mostShrink[i] * 2);
		}

		float len = 0;
		for (float &d : tmpDiff)len += d*d;
		if (len == 0)break;
		len = sqrt(len);
		for (float &d : tmpDiff)d /= len;

		for (int i = 0; i < paraNum; i++) {
			if (oldDiff[i] > 0.5&&tmpDiff[i] < -0.5&&abs(oldDiff[i] + tmpDiff[i])<.01f)
				mostShrink[i]++;
			if (oldDiff[i] < 0.5&&tmpDiff[i] > -0.5&&abs(oldDiff[i] + tmpDiff[i])<.01f)
				mostShrink[i]++;

			oldDiff[i] = tmpDiff[i];
			oldVel[i] = tmpVel[i];
		}

		for (int i = 0; i < paraNum; i++) {
			tmpVel[i] = RATIO * tmpVel[i] - descStep * tmpDiff[i];
			tmpPara[i] += tmpVel[i] + RATIO * (tmpVel[i] - oldVel[i]);
			tmpPara[i] = clamp(0.f, 1.f, tmpPara[i]);
		}

		if(miniStep != NULL)descStep = miniStep(tmpLoss, descStep);
		if (breakCond != NULL && breakCond(tmpLoss))break;
	}
	return clock() - t;
}
void Grad::print() {
	for (float p : tmpPara)cout << p << " ";
	cout << endl;
}

float variant(vector<float> &std, vector<float> &cal) {
	float res = 0.f;

#pragma omp parallel for
	for (unsigned int i = 0; i < std.size(); i++) {
		res += (std[i] - cal[i])*(std[i] - cal[i]);
	}
	return res;
}
float clamp(float floor, float ceil, float d) {
	if (d < floor)return floor;
	if (d > ceil)return ceil;
	return d;
}

main.cpp

#include "gradient.h"

vector<float> func(vector<float>& para) {
	vector<float> ret;
	for (float x = 0.f; x < 1.f; x += .1f) {
		for (float y = 0.f; y < 1.f; y += .1f) {
			for (float z = 0.f; z < 1.f; z += .1f) {
				ret.push_back(para[0] * x + para[1] * y + para[2] * z);
			}
		}
	}
	return ret;
}
float minimize(float var, float step) {
	if (var<.1f&&step>.001f)step /= 10;
	if (var<.001f&&step>.0001f)step /= 10;
	if (var<.00001f&&step>.00001f)step /= 10;
	return step;
}
bool quit(float var) {
	if (var < .00000001f)return true;
	else return false;
}

int main() {
	vector<float> p;
	p.push_back(.8f);
	p.push_back(.2f);
	p.push_back(.3f);
	vector<float> s = func(p);

	Grad grad(3, .01f, .0001f, func, true, true, 100, variant, minimize);
	grad.standard(s);
	grad.prepare();
	cout<<grad.desc();

	system("pause");
}