手写数据集非线性最小二乘二分类实现_最小二乘数据集

---

一、非线性最小二乘

对于非线性等式
$$f_i(x)=g_i(x)-y_i=0,i=1,...,m$$
写作向量形式$f(x)=(f_1(x),...,f_m(x))$
$$f(x)=0$$
对任意x,我们称f_i(x)为第i个残差，我们不能直接得到x的数值解，但是可以通过使残差平方和最小，来求一个近似解。
$$f_1(x{)}^2+...+f_m(x{)}^2=||f(x)|{|}^2$$
这意味着找到对任意x具有$||f(x)|{|}^2\ge ||f(\tilde{x})|{|}^2$的$\tilde{x}$,我们把这一点$\tilde{x}$称为式（1）的最小二乘近似解。

将博客手写数据集的线性最小二乘二分类问题转化为非线性问题(这里将原文中的参数换为$\beta$,未知数为x不然不好看)：
$$g={\hat{x}}^T\hat{\beta }$$
$$f_i(x)=g_i(x)-y_i$$
将$f_i(x)$做归一化处理
$$\phi (u)=\frac{e^u-e^{-u}}{e^u+e^{-u}}$$
我们通过求解非线性最小二乘问题来选择x
$$\sum _{i=1}^m\phi (g_i(x))-y_i$$

二、Levenberg-Marquardt

1.牛顿法

对于一维$f(x)=0$求根，牛顿法是一种利用泰勒展开的线性化方法

对于优化问题
$$minf(x)$$
问题转化为
$$f^{{}^{\prime }}(x)=0$$
牛顿法将问题转化为求
$$f^{{}^{\prime }}(x_0)+(x-x_0)f^{{}^{\prime \prime }}(x_0)=0$$
同样用迭代法求解x
$$x_{n+1}=x_n-\frac{f^{{}^{\prime }}(x_n)}{f^{{}^{\prime \prime }}(x_n)}$$
对于多维的$f(x)=0$问题，这时我们就可以利用雅克比，海赛矩阵之类的来表示高维求导函数。
$$f(x)=0，x=(x_0,...,x_n)$$
雅克比矩阵(表示一阶偏导)：
$$J_f=(\frac{\partial f}{\partial x_0},\frac{\partial f}{\partial x_1}.........\frac{\partial f}{\partial x_n})$$
黑赛矩阵表示二阶偏导
H f = [ ∂ 2 f ∂ x 0 2 ∂ 2 f ∂ x 0 x 1 . . . ∂ 2 f ∂ x 0 x n ∂ 2 f ∂ x 1 x 0 ∂ 2 f ∂ x 1 2 . . . ∂ 2 f ∂ x 1 x n . . . . ∂ 2 f ∂ x n x 0 ∂ 2 f ∂ x n x 1 . . . ∂ 2 f ∂ x n 2 ] H_f=

$$\begin{bmatrix} \frac{\partial^2 f}{\partial x_0^2} & \frac{\partial^2 f}{\partial x_0x_1}&...&\frac{\partial^2 f}{\partial x_0x_n}\\ \frac{\partial^2 f}{\partial x_1x_0} & \frac{\partial^2 f}{\partial x_1^2}&...&\frac{\partial^2 f}{\partial x_1x_n}\\ .&.&.&.\\\frac{\partial^2 f}{\partial x_nx_0} & \frac{\partial^2 f}{\partial x_nx_1}&...& \frac{\partial^2 f}{\partial x_n^2}\\ \end{bmatrix}$$ Hf​=⎣⎢⎢⎢⎢⎡​∂x02​∂2f​∂x1​x0​∂2f​.∂xn​x0​∂2f​​∂x0​x1​∂2f​∂x12​∂2f​.∂xn​x1​∂2f​​..........​∂x0​xn​∂2f​∂x1​xn​∂2f​.∂xn2​∂2f​​⎦⎥⎥⎥⎥⎤​
高维牛顿法解最优化问题又可写成
$$x_{n+1}=x_n-H_f(x_n{)}^{-1}{J}_f(x_n)$$

2.高斯牛顿法

对于目标函数
$$min\sum _{i=1}^m{f}_i^2(x)$$
雅克比矩阵(表示一阶偏导)：
J f = [ ∂ f 0 ∂ x 0 ∂ f 0 ∂ x 1 . . . ∂ f 0 ∂ x n ∂ f 1 ∂ x 0 ∂ f 1 ∂ x 1 . . . ∂ f 1 ∂ x n . . . . ∂ f m ∂ x 0 ∂ f m ∂ x 1 . . . ∂ f m ∂ x n ] J_f=

$$\begin{bmatrix} \frac{\partial f_0}{\partial x_0} & \frac{\partial f_0}{\partial x_1}&...&\frac{\partial f_0}{\partial x_n}\\ \frac{\partial f_1}{\partial x_0} & \frac{\partial f_1}{\partial x_1}&...&\frac{\partial f_1}{\partial x_n}\\ .&.&.&.\\ \frac{\partial f_m}{\partial x_0} & \frac{\partial f_m}{\partial x_1}&...&\frac{\partial f_m}{\partial x_n}\\ \end{bmatrix}$$ Jf​=⎣⎢⎢⎢⎡​∂x0​∂f0​​∂x0​∂f1​​.∂x0​∂fm​​​∂x1​∂f0​​∂x1​∂f1​​.∂x1​∂fm​​​..........​∂xn​∂f0​​∂xn​∂f1​​.∂xn​∂fm​​​⎦⎥⎥⎥⎤​
黑赛矩阵表示二阶偏导(将二次偏导省略)
$$H=2J_f^T{J}_f$$
代入牛顿法高维迭代方程的基本形式，得到高斯牛顿法迭代方程
$$x_{n+1}=x_n-（J_f^T{J}_f{）}^{-1}{J}_f(x_n)f(x_n)$$

3.Levenberg-Marquardt算法

高斯牛顿法的H矩阵并不一定可逆，Levenberg-Marquardt对其进行改进
$$x_{n+1}=x_n-（J_f^T{J}_f+\lambda I{）}^{-1}{J}_f(x_n)f(x_n)$$
然后Levenberg-Marquardt方法的好处就是在于可以调节:
如果下降太快，使用较小的λ，使之更接近高斯牛顿法
如果下降太慢，使用较大的λ，使之更接近梯度下降法
细节：Levenberg-Marquardt算法浅谈
因此在解决了高斯-牛顿法基础之上，LM算法的重点就是如何确定λ值。λ值更新
算法流程

代码实现

其实有很多地方没有详尽的推导，代码也可能不正确，后面重新捋一遍（给自己说）
参考博客LM算法C++实现
手写数据集非线性最小二乘二分类实现

#include <vector> 
#include <string> 
#include <fstream> 
#include <iostream>
#include <Eigen/Dense> 
#include <ctime>
#include <windows.h>
using namespace std;
using namespace Eigen;

const double DERIV_STEP = 1e-5;
const int MAX_ITER = 100;

int ReversInt(int i)//保证数据在不同主机之间传输时能够被正确解释
{
	unsigned char ch1, ch2, ch3, ch4;//无符号版本和有符号版本的区别就是有符号类型需要使用一个bit来表示数字的正负
	ch1 = i & 255;
	ch2 = (i >> 8) & 255;
	ch3 = (i >> 16) & 255;
	ch4 = (i >> 24) & 255;
	return((int)ch1 << 24) + ((int)ch2 << 16) + ((int)ch3 << 8) + ch4;
}
void read_Mnist_label(string filename, vector<float>& labels)
{
	ifstream file(filename, ios::binary);
	if (file.is_open())
	{
		int magic_number = 0;
		int number_of_images = 0;
		file.read((char*)&magic_number, sizeof(magic_number));
		file.read((char*)&number_of_images, sizeof(number_of_images));
		magic_number = ReversInt(magic_number);
		number_of_images = ReversInt(number_of_images);
		cout << endl;
		cout << "\t" << "magic number = " << magic_number << endl;
		cout << "\t" << "number of images = " << number_of_images << endl;
		for (int i = 0; i < number_of_images; i++)
		{
			unsigned char label = 0;//一个字节存储8位无符号数，储存的数值范围为0-255
			file.read((char*)&label, sizeof(label));
			labels.push_back((float)label);
		}
	}
}
void read_Mnist_Images(string filename, vector<vector<float>>& images)
{
	ifstream file(filename, ios::binary);//直接调用了其默认的打开方式
	if (file.is_open())
	{
		int magic_number = 0;
		int number_of_images = 0;
		int n_rows = 0;
		int n_cols = 0;
		// unsigned char label;

		file.read((char*)&magic_number, sizeof(magic_number));
		file.read((char*)&number_of_images, sizeof(number_of_images));
		file.read((char*)&n_rows, sizeof(n_rows));
		file.read((char*)&n_cols, sizeof(n_cols));
		magic_number = ReversInt(magic_number);
		number_of_images = ReversInt(number_of_images);
		n_rows = ReversInt(n_rows);
		n_cols = ReversInt(n_cols);
		cout << endl;
		cout << "\t" << "magic number = " << magic_number << endl;
		cout << "\t" << "number of images = " << number_of_images << endl;//60000
		cout << "\t" << "rows = " << n_rows << endl;//28
		cout << "\t" << "cols = " << n_cols << endl;//28
		for (int i = 0; i < number_of_images; i++)
		{
			vector<float>tp;
			for (int r = 0; r < n_rows; r++)
			{
				for (int c = 0; c < n_rows; c++)
				{
					unsigned char image = 0;//一个字节存储8位无符号数，储存的数值范围为0-255
					file.read((char*)&image, sizeof(image));
					tp.push_back(image);
				}
			}
			images.push_back(tp);//60000×784
		}
	}
}
float f_i(VectorXd* beta, MatrixXf* x, int i)//f_i
{
	float result = 0.0f;
	float sum = 0.0f;
	for (int j = 0; j < 494; j++)
		sum += ((*x)(i, j) * (*beta)(j));
	result += (exp(sum) - exp(-sum) / (exp(sum) + exp(-sum)));
	return result;
}
float func_i(VectorXf* labels1, VectorXd* beta, MatrixXf* x, int i)//f_i
{
	float result = 0.0f;
	float sum = 0.0f;
	for (int j = 0; j < 494; j++)
		sum += ((*x)(i, j) * (*beta)(j));
	result += (exp(sum) - exp(-sum) / (exp(sum) + exp(-sum)) - (*labels1)[i]);
	return result;
}
float FUNC(VectorXf* labels1, VectorXd* beta, MatrixXf* x)
{
	float sum = 0;
	for (int i = 0; i < (*labels1).size(); i++)
	{
		sum += func_i(labels1, beta, x, i);
	}
	return sum;
}
float Deriv(VectorXf* labels1, VectorXd* beta, MatrixXf* x, int n, int i)//beta_n,f_i
{
	// Assumes input is a single row matrix
	//clone方法返回，一个新的相同的对象被创建
	// Returns the derivative of the nth parameter
	VectorXd params1(*beta);
	VectorXd params2(*beta);

	// Use central difference  to get derivative
	params1[n] -= DERIV_STEP;
	params2[n] += DERIV_STEP;

	// c += a 等效于 c = c + a 
	float p1 = func_i(labels1, &params1, x, i);
	float p2 = func_i(labels1, &params2, x, i);

	float d = (p2 - p1) / (2 * DERIV_STEP);

	return d;
}
int main(int argc, char* argv[])
{
	Eigen::initParallel();
	//训练
//60000×784，组成一个矩阵
	vector<vector<float>>images;
	//60000×494，组成一个矩阵
	vector<int> effective_col;
	//一个60000维标签向量
	vector<float>labels;
	//一个494维b向量
	cout << endl << "训练集图片信息：" << endl;
	read_Mnist_Images("train-images.idx3-ubyte", images);
	cout << endl << "训练集标签信息：" << endl;
	read_Mnist_label("train-labels.idx1-ubyte", labels);
	//60000中的600张
	for (int i = 0; i < images[0].size(); i++)
	{
		int s = 0;
		for (int j = 0; j < images.size(); j++)
		{
			if (images[j][i] > 0)
			{
				s++;
			}
			if (s == 600)
			{
				effective_col.push_back(i);//得到有效列
				break;
			}
		}
	}
	auto m = images.size();//训练矩阵的行数60000
	auto n = images[0].size();//矩阵列数
	auto b = labels.size();//矩阵行数
	auto num_key = effective_col.size() + 1;//有效列数


	//***************************************建立训练集60000*494矩阵***************************************
	MatrixXf Characteristic_matrix(m, num_key);//其每个元素都是float类型。
	for (int i = 0; i < m; i++)
	{
		for (int j = 0; j < num_key - 1; j++)
		{
			Characteristic_matrix(i, j) = images[i][effective_col[j]];
		}
		Characteristic_matrix(i, num_key - 1) = 1;//最后一个特征是1
	}
	images.~vector<vector<float>>();//析构函数，释放内存
	int tr = Characteristic_matrix.rows();
	int tc = Characteristic_matrix.cols();
	VectorXf actual_Y(b);//初始化训练集的标签向量
	for (int i = 0; i < b; i++)
	{
		labels[i] == 0 ? actual_Y(i) = 1 : actual_Y(i) - 1;
	}
	labels.~vector<float>();//析构函数，释放内存

	//初始化beta矩阵
	VectorXd beta(num_key);
	for (int i = 0; i < 494; i++)
	{
		beta(i) = 0.5f;
	}
	VectorXd tempbeta(num_key);
	MatrixXd J(m, num_key);//雅各比矩阵
	VectorXd G(num_key);//g矩阵
	float u = 1.0f;//阻尼因子
	float v = 2.0f;
	MatrixXd I(num_key, num_key);
	I = I.Identity(num_key, num_key);//单位矩阵
	float q = 1.0f;
	MatrixXd H(num_key, num_key);//黑塞矩阵
	float mse;
	float mse_temp;
	VectorXd f(m);
	for (int k = 0; k < 10; k++)
	{
		cout << k << endl;
		if (q > 0)
		{
			for (int i = 0; i < m; i++)
			{
				f(i) = func_i(&actual_Y, &beta, &Characteristic_matrix, i);
			}
			for (int i = 0; i < m; i++)
			{
				for (int j = 0; j < 494; j++)
					J(i, j) = Deriv(&actual_Y, &beta, &Characteristic_matrix, j, i);
				G = J.transpose() * f;
			}
			H = J.transpose() * J;
		}
		tempbeta = beta - ((H + u * I).inverse()) * G;
		mse = FUNC(&actual_Y, &beta, &Characteristic_matrix);
		mse_temp = FUNC(&actual_Y, &tempbeta, &Characteristic_matrix);
		q = 2 * (mse - mse_temp) / (beta.transpose() * (u * beta - J.transpose() * f));
		if (q > 0)
		{
			mse = mse_temp;
			beta = tempbeta;
			float a = 1 / 3;
			u = u * max(a, 1 - (2 * q - 1) * (2 * q - 1) * (2 * q - 1));
			v = 2;
		}
		else
		{
			u = u * v;
			v = v * 2;
		}
		cout << beta << endl;
		if (fabs(mse - mse_temp) < 1e-8)
			break;
	}

	vector<vector<float>>images1;
	vector<float>labels1;
	std::cout << endl << "测试集图片信息：" << endl;
	read_Mnist_Images("t10k-images.idx3-ubyte", images1);
	std::cout << endl << "测试集标签信息：" << endl;
	read_Mnist_label("t10k-labels.idx1-ubyte", labels1);
	auto m1 = images1.size();
	std::cout << m1 << endl;

	MatrixXf Characteristic_matrix1(m1, num_key);//其每个元素都是float类型。
	for (int i = 0; i < m1; i++)
	{
		for (int j = 0; j < num_key - 1; j++)
		{
			Characteristic_matrix1(i, j) = images1[i][effective_col[j]];
		}
		Characteristic_matrix1(i, num_key - 1) = 1;
	}
	images1.~vector<vector<float>>();//析构函数，释放内存
	VectorXf actual_y1(m1);//初始化测试集的标签向量
	for (int i = 0; i < m1; i++)
	{
		labels1[i] == 0 ? actual_y1(i) = 1 : actual_y1(i) = -1;
	}
	labels1.~vector<float>();//析构函数，释放内存


	//***************************************建立测试集10000*494矩阵***************************************


	//写判断函数
	int tp = 0, fp = 0, tn = 0, fn = 0;
	for (int i = 0; i < m1; i++)
	{
		int c = (int)(f_i(&tempbeta, &Characteristic_matrix1, i) + 0.5);
		std::cout << c << endl;
		if (c == actual_y1(i))
		{
			tp++;
		}
		else
		{
			fp++;
		}
	}
	std::cout << endl;
	std::cout << "测试集预测结果：" << endl;
	std::cout << "\t" << "tp = " << tp << endl
		<< "\t" << "fp = " << fp << endl;
	std::cout << "\t" << "错误率 = " << 1.0 * (fp) / 10000 << endl;
	return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298