c++实现线性回归（高斯消元）（附python实现）

最新推荐文章于 2022-05-30 12:43:19 发布

luke2834

最新推荐文章于 2022-05-30 12:43:19 发布

阅读量564

点赞数

分类专栏： ACM_线性代数

本文链接：https://blog.csdn.net/luke2834/article/details/82745035

版权

ACM_线性代数专栏收录该内容

1 篇文章 0 订阅

订阅专栏

前言

写这次blog的契机是上次笔试的时候，遇到了这个问题
当时以为numpy库是可以用的，就先写了个python版，结果并不能用。。
最后愤然写了个c++版
不过最后一个小问题导致我差了两分钟没交上去代码，所以这一版源码只是通过了案例但没有提交ac。。如果大家发现有什么bug也欢迎大家帮我指出

题意：

就是给定一个X矩阵，一个y向量，分别表示数据特征矩阵（n*m），ground truth向量（n*1）
输入就是一个大矩阵，最后一列是y，前面是X，同时原题里n，m是先不给定的，需要根据矩阵的行数和列数判断
要求返回系数向量，包含偏置b，所以返回的是一个n+1维的vector
我现在找不到原题了，找了一个类似的简化的题，记在这里：链接

思路：

实现的话，其实就是通过高斯消元，解一个线性方程组：
$X^TXw = X^Ty$
这里把推导过程简单说一下
首先我们知道线性回归的loss function：
$min_ w L = ||y - Xw||$
然后我们求关于w的梯度:
$L = (y - Xw)^T(y - Xw) \\ = (y^Ty - w^TX^Ty - y^TXw + w^TX^TXw)$
$\nabla_w L = -X^Ty - X^Ty + 2X^TXw$
令梯度为0，容易推出最后结论：
$X^TXw = X^Ty$
最后附上矩阵求导参考文献：链接

c++实现

#include <bits/stdc++.h>
using namespace std;

const double eps = 1e-8;
typedef vector<double> vec;
typedef vector<vec> mat;

vec gauss_jordan(const mat& A, const vec& b){
	int n = A.size();
	mat B(n, vec(n+1));
	for (int i = 0; i < n; i++){
		for (int j = 0; j < n; j++){
			B[i][j] = A[i][j];
		}
	}
	for (int i = 0; i < n; i++){
		B[i][n] = b[i];
	}
	for (int i = 0; i< n; i++){
		int pivot = i;
		for (int j = i; j < n; j++){
			if (abs(B[j][i]) > abs(B[pivot][i]))
				pivot = j;
		}
		swap(B[i], B[pivot]);
	
		if (abs(B[i][i]) < eps)
			return vec();
	
		for (int j = i + 1; j <= n; j++){
			B[i][j] /= B[i][i];
		}
		for (int j = 0; j < n; j++){
			if (i != j){
				for (int k = i + 1; k <= n; k++)
					B[j][k] -= B[j][i] * B[i][k];
			}
		}
	}
	vec x(n);
	for (int i = 0; i < n; i++)
		x[i] = B[i][n];
	return x;
}

mat mul(mat& A, mat& B){
	mat C(A.size(), vec(B[0].size()));
	for (int i = 0; i < A.size(); i++){
		for (int k = 0; k < B.size(); k++){
			for (int j = 0; j < B[0].size(); j++){
				C[i][j] += A[i][k] * B[k][j];
			}
		}
	}
	return C;
}

mat trans(mat& A){
	mat B(A[0].size(), vec(A.size()));
	for (int i = 0; i < A.size(); i++){
		for (int j = 0; j < A[i].size(); j++){
			B[j][i] = A[i][j];
		}
	}
	return B;
}
vec tovec(mat& A){
	vec B(A.size());
	for (int i = 0; i < A.size(); i++){
		B[i] = A[i][0];
	}
	return B;
}

void SplitString(const string& s, vector<double>& v, const string& c)
{
  string::size_type pos1, pos2;
  pos2 = s.find(c);
  pos1 = 0;
  while(string::npos != pos2)
  {
    v.push_back(stod(s.substr(pos1, pos2-pos1)));
 
    pos1 = pos2 + c.size();
    pos2 = s.find(c, pos1);
  }
  if(pos1 != s.length())
    v.push_back(stod(s.substr(pos1)));
}
int main(){
	string s;
	ios::sync_with_stdio(false);
	mat X, Y;
	for (int i = 0; i < 4; i++){
		cin>>s;
		vector<double> v;
		SplitString(s, v, ",");
		vector<double> v1(v.begin(), v.end() - 1);
		v1.insert(v1.begin(), 1.0);
		vector<double> v2(1, *v.rbegin());
		X.push_back(v1);
		Y.push_back(v2);
	}
	for (int i = 0; i < X.size(); i++){
		for (int j = 0; j < X[i].size(); j++){
			cout << X[i][j] << ", ";
		}
		cout << endl;
	}
	for (int i = 0; i < Y.size(); i++){
		for (int j = 0; j < Y[i].size(); j++){
			cout << Y[i][j] << ", ";
		}
		cout << endl;
	}
	auto Xt = trans(X);
	auto A = mul(Xt, X);
	auto B = mul(Xt, Y);
	auto C = tovec(B);
	auto ans = gauss_jordan(A, C);
	for (int i = 0; i < ans.size() - 1; i++){
		cout << fixed<<setprecision(2)<< ans[i] << ",";
	}
	cout << fixed<<setprecision(2) << *ans.rbegin() << endl;
	return 0;
}

python实现

import sys, math
import numpy as np


def solve(X,Y):
    A = np.dot(X.T, X)
    B = np.dot(X.T, Y.T)
    W = np.squeeze(np.linalg.solve(A, B).T).tolist()
    return W


if __name__=='__main__':
    X = []
    Y = []
    while True:
        line = sys.stdin.readline().strip('\r\n')

        if line == '':
            break

        nums = line.split(',')
        nums = [float(n) for n in nums]
        X.append([1.0] + nums[:-1])
        Y.append(nums[-1])

    X = np.array(X, dtype = np.float32)
    Y = np.array(Y, dtype = np.float32, ndmin = 2)
    coefs = solve(X,Y)

    formatted_coefs = []
    for coef in coefs:
        coef = math.floor(coef*100)/100
        formatted_coefs.append('%.2f' %coef)
    print(','.join(formatted_coefs))