最小二乘回归R绘图 最小二乘回归的假设

文章目录

• 前言
• 一、最小二乘法的概念
• 二、数字集MNIST二分类

• 1.题目介绍
• 2.分析
• 2.最小二乘法思路

• 三、代码实现
• 总结

前言

残差：残差在数理统计中是指实际观察值与估计值（拟合值）之间的差。“残差”蕴含了有关模型基本假设的重要信息。如果回归模型正确的话， 我们可以将残差看作误差的观测值。

一、最小二乘法的概念

普通最小二乘（Ordinary least squares ）估计通常用于分析实验和观测数据，OLS方法最小化残差的平方和，并导致未知参数向量β的估计值的封闭形式表达式。
$\hat{\beta}=(X^TX)^{-1}X^T\hat{y}\tag{1}$
线性最小二乘(Linear least squares)是线性函数对数据的最小二乘近似。它是一套解决线性回归中涉及的统计问题的公式，包括普通(未加权)、加权和广义(相关)残差的变量。线性最小二乘的数值求解方法包括正交分解法和正交方程组的矩阵求逆法。

二、数字集MNIST二分类

1.题目介绍

对于手写数字集MNIST进行处理，实现0和非0元素的分类，手写数字集MNIST中包含60000个训练集和10000个测试集，共包括70000 张 28 ×28的手写数字灰度图像，图像数据已经被转换为28 × 28 = 784维的向量形式存储，标签对应的为10维向量存储，如数字3对应的标签为[0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0.0.0]

2.分析

对于手写数字的图片来说，数字大多数在图像的中央，而周围有很多冗余的区域，对于图像识别来说意义不大，每张图片28×28= 784个存放像素的位置，故在60000张图片中，如果有600张图片在该位置像素值不等于0，我们就提取这个位置作为特征，经过这种方法提取后，我们可以算出来有493个特征位置。

$用作特征的像素的位置$

2.最小二乘法思路

线性回归，最小二乘，矩阵条件数，QR分解的理解联系在这篇文章中讲得很清楚https://zhuanlan.zhihu.com/p/60743204

我们定义$y=f(x_1,x_2,...,x_n)$，$y$为实数，标签为0图片$y=1$，标签非0图片$y=-1$，$\hat{x}=x_1,x_2,...,x_{493}$为特征位置上的像素值。
考虑这样一个问题：给定60000个样本点$(\hat{x_1},y_1)，(\hat{x_2},y_2)...(\hat{x}_{60000},y_{60000})$,利用样本点建立如下关系（线性回归模型），让我们根据输入图片，来预测$y$的值
$\overline{y}=\beta_0+\sum_{i=1}^{493}\beta_ix_1$
为方便表示，记权重向量 $\hat{\beta}=[\beta_0,\beta_1,\beta_2...\beta_{493}]^T$，输入向量$\hat{x}=[1,x_1,x_2...x_{493}]^T$，则线性回归的，目标即找到如下关系：
$\overline{y}=\hat{\beta}^T\hat{x}$
解这个模型我们希望$\hat{\beta}$使预测误差$y-\overline{y}$尽可能的小，即求其最小二乘解，形成如下无约束的优化问题(线性最小二乘问题)：
$min\sum_{i=1}^{493}(\hat{\beta}^T\hat{x_i}-y_i)^2$
解最小二乘问题常用方法（推导这张开头的链接中有重点提到）：
1.正规方程 (normal equations) & 伪逆（式（1）的来源）
2.最小二乘问题的QR分解(full size & economy size)
当我还没有查资料时，第一次看到伪逆这个公式的时候充满不相信,所以手算了一下

最小二乘问题：$I（\beta_0,\beta_1,\beta_2,...,\beta_{493}）=min\sum_{j=1}^{60000}(\hat{\beta}^T\hat{x_j}-y_j)^2$
$I$对$\beta_k$求导
$\frac{\partial I}{\partial \beta}=2\sum_{j=1}^{60000}(\hat{\beta}^T\hat{x_j}-y_j)x_{jk}=0$
$2\sum_{j=1}^{60000}\hat{\beta}^T\hat{x_j}x_{jk}=2\sum_{j=1}^{60000}\sum_{i=0}^{493}\beta_ix_{ji}x_{jk}==2\sum_{i=0}^{493}\beta_i\sum_{j=1}^{60000}x_{ji}x_{jk}=2\sum_{j=1}^{60000}y_jx_{jk}$
可以写成方阵的形式
$\begin{bmatrix} \sum_{j=1}^{60000}x_{j,0}x_{j,0} & \sum_{j=1}^{60000}x_{j,1}x_{j0} & ...& \sum_{j=1}^{60000}x_{j,493}x_{j,0} \\ \\ \ \sum_{j=1}^{60000}x_{j,0}x_{j,1} & \sum_{j=1}^{60000}x_{j,1}x_{j,1} & ...& \sum_{j=1}^{60000}x_{j,493}x_{j,1} \\&......... ...& \\ \sum_{j=1}^{60000}x_{j,0}x_{j,493} & \sum_{j=1}^{60000}x_{j,1}x_{j,493} & ...& \sum_{j=1}^{60000}x_{j,493}x_{j,493} \\ \end{bmatrix} \begin{bmatrix} \beta_0 \\\\ \beta_1 \\...\\\beta_{493} \\ \end{bmatrix}=\begin{bmatrix} \sum_{j=1}^{60000}y_jx_{j,0} \\\\ \sum_{j=1}^{60000}y_jx_{j,2} \\...\\\sum_{j=1}^{60000}y_jx_{j,493}\\ \end{bmatrix}$
令
$X=\begin{bmatrix}x_{1,1} & x_{1,1}& x_{1,2}&...& x_{1,493} \\ x_{2,1}& x_{2,1}& x_{2,2}&...& x_{2,493} \\&&... \\x_{60000,1}& x_{60000,1}& x_{60000,2}&...& x_{60000,493} \\ \end{bmatrix}$
矩阵形式可写为
$\hat{\beta}=(X^TX)^{-1}X^T\hat{y}$

三、代码实现

#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, ¶ms1, x, i);
	float p2 = func_i(labels1, ¶ms2, 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="<< 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;
				cout << i << endl;
			}
			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;
		}
		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;
}