C++元编程——线性回归

腾昵猫

已于 2023-05-12 19:48:57 修改

阅读量278

点赞数

文章标签：线性回归算法回归

于 2023-05-11 18:24:51 首次发布

本文链接：https://blog.csdn.net/Dr_Jack/article/details/130626488

版权

如下就是线性回归的公式：

$\vec{y}=\vec{x}\cdot\vec{w}$

$\overrightarrow{w}=(w_{0},w_{1}, ...,w_{n},b)^{'}$

$\overrightarrow{x}=(x_{0},x_{1},...,x_{n},1)$

如果将误差函数设置为：

$E=\left \| y-\vec{x}\cdot\vec{w} \right \|_{2}^{2}$

则其最小值点即为

$\vec{y}=\vec{x}\cdot\vec{w}$

$\vec{x}^{'}\cdot\vec{y}=\vec{x}^{'}\cdot\vec{x}\cdot\vec{w}$

也即

$\vec{w}=(\vec{x}^{'}\cdot\vec{x})^{-1} \vec{x}^{'}\cdot y$

所以根据这个公式可以写代码如下：

矩阵定义

#ifndef _MAT_HPP_
#define _MAT_HPP_
#include <memory>
#include <climits>
#include <float.h>

#include <map>
#include <iostream>
#include <iomanip>
#include "ht_memory.h"
#ifdef WITH_BOOST
#include <boost/pool/pool.hpp>
#endif


template<typename val_t>
val_t max_and_swap(const val_t& v1, const val_t& v2)
{
	return v1 < v2 ? v2 : v1;
}

template<typename type>
struct destoryer 
{
	static void do_destory(type& v)
	{
		//printf("%s\r\n", typeid(type).name());
	}
};

template<int i_size, typename val_t>
struct mat_m
{
#ifdef WITH_BOOST
	static boost::pool<> s_pool;
#else
	val_t sz_ele[i_size];
#endif
	val_t* p;

	mat_m() :p(nullptr)
	{
#ifndef WITH_BOOST
		//p = (val_t*)malloc(i_size * sizeof(val_t));
		p = sz_ele;
#else
		p = (val_t*)(s_pool.malloc());
#endif
		for (int i = 0; i < i_size; ++i)
		{
			new(p + i) val_t(0);
			//p[i] = 0;
		}
	}
	~mat_m()
	{
		if (p)
		{
#ifndef WITH_BOOST
			/*
			for (int i = 0; i < i_size; ++i)
				destoryer<val_t>::do_destory(p[i]);
			free(p);
			*/

#else
			for (int i = 0; i < i_size; ++i)
				destoryer<val_t>::do_destory(p[i]);
				//p[i].~val_t();
			s_pool.free(p);
#endif
			p = 0;
		}
	}
	val_t& get(const int& len_1d, const int& i_1d_idx, const int& i_2d_idx)
	{
		val_t& ret = p[i_2d_idx + len_1d * i_1d_idx];
		return ret;
	}

	val_t max_abs() const
	{
		double d = -1*DBL_MAX;
		for (int i = 0; i < i_size; ++i) 
		{
			d = d < abs(p[i]) ? abs(p[i]) : d;
		}
		return d;
	}

	val_t max() const
	{
		val_t d = -1 * DBL_MAX;
		for (int i = 0; i < i_size; ++i)
		{
			//d = d < (p[i]) ? (p[i]) : d;
			d = max_and_swap(p[i], d);
		}
		return d;
	}

	val_t sum() const 
	{
		val_t d_sum = 0.;
		for (int i = 0; i < i_size; ++i)
		{
			d_sum = d_sum + p[i];
		}
		return d_sum;
	}

	template<int len_1d, int i_1d_idx, int i_2d_idx>
	inline val_t& get_val()
	{
		static_assert((i_2d_idx + len_1d * i_1d_idx) < i_size, "ERROR:mat_m over flow！！！");
		return p[i_2d_idx + len_1d * i_1d_idx];
	}

	template<int len_1d, int i_1d_idx, int i_2d_idx>
	inline val_t get_val() const
	{
		return p[i_2d_idx + len_1d * i_1d_idx];
	}
};

#ifdef WITH_BOOST
template<int i_size, typename val_t>
boost::pool<> mat_m<i_size, val_t>::s_pool(i_size * sizeof(val_t));
#endif

template<int row_num, int col_num, typename val_t = double>
struct mat
{
	using t_type = mat<col_num, row_num, val_t>;
	using type = val_t;
	typedef val_t vt;
	static constexpr int r = row_num;
	static constexpr int c = col_num;
	using mat_m_t = mat_m<row_num * col_num, val_t>;
	std::shared_ptr<mat_m_t> pval;
	bool b_t;

	mat():b_t(false)
	{
		pval = std::make_shared<mat_m_t>();
	}
	mat(const mat<row_num, col_num, val_t>& other) :pval(other.pval), b_t(other.b_t)
	{
	}
	~mat() 
	{
		pval.reset();
	}

	mat(const val_t&& v):b_t(false)
	{
		pval = std::make_shared<mat_m_t>();
		for (int i = 0; i < row_num; ++i)
		{
			for (int j = 0; j < col_num; ++j)
			{
				pval->get(col_num, i, j) = v;
			}
		}
	}
	mat(const val_t& v) :b_t(false)
	{
		pval = std::make_shared<mat_m_t>();
		for (int i = 0; i < row_num; ++i)
		{
			for (int j = 0; j < col_num; ++j)
			{
				pval->get(col_num, i, j) = v;
			}
		}
	}
#if 0
	template<typename val_other_t>
	mat(const val_other_t& v) :b_t(false)
	{
		pval = std::make_shared<mat_m_t>();
		for (int i = 0; i < row_num; ++i)
		{
			for (int j = 0; j < col_num; ++j)
			{
				pval->get(col_num, i, j) = static_cast<val_t>(v);
			}
		}
	}
#endif
	mat(const std::initializer_list<val_t>& lst):b_t(false)
	{
		pval = std::make_shared<mat_m_t>();
		auto itr = lst.begin();
		for (int i = 0; i < row_num; ++i)
		{
			for (int j = 0; j < col_num; ++j)
			{
				if (itr == lst.end())return;
				pval->get(col_num, i, j) = *itr;
				itr++;
			}
		}
	}

	val_t& get(const int& i_row, const int& i_col)
	{
		if (!b_t)
			return pval->get(col_num, i_row, i_col);
		else
			return pval->get(row_num, i_col, i_row);
	}

	val_t get(const int& i_row, const int& i_col) const
	{
		if (!b_t)
			return pval->get(col_num, i_row, i_col);
		else
			return pval->get(row_num, i_col, i_row);
	}

	template<int i_1d_idx, int i_2d_idx>
	inline val_t& get_val()
	{
		if (!b_t)
			return pval->template get_val<col_num, i_1d_idx, i_2d_idx>();
		else
			return pval->template get_val<row_num, i_2d_idx, i_1d_idx>();
	}

	template<int i_1d_idx, int i_2d_idx>
	inline val_t get_val() const
	{
		static_assert(i_1d_idx < row_num && i_2d_idx < col_num, "ERROR: mat::get_val overflow!!!!!");
		if (!b_t)
			return pval->template get_val<col_num, i_1d_idx, i_2d_idx>();
		else
			return pval->template get_val<row_num, i_2d_idx, i_1d_idx>();
	}

	mat<col_num, row_num, val_t> t() const
	{
		mat<col_num, row_num, val_t> ret;
		ret.pval = pval;
		ret.b_t = !b_t;
		return ret;
	}

	val_t max_abs() const
	{
		return pval->max_abs();
	}

	val_t max() const 
	{
		return pval->max();
	}

	val_t sum() const 
	{
		return pval->sum();
	}

	void print() const
	{
		std::cout << "[" << std::endl;
		for (int i = 0; i < row_num; ++i)
		{
			std::cout << std::setw(3) << "[";
			for (int j = 0; j < col_num; ++j)
			{
				std::cout << (j != 0 ? "," : "") << std::setw(10) << get(i, j);
			}
			std::cout << std::setw(3) << "]" << std::endl;
		}
		std::cout << "]" << std::endl;
	}

	template<int other_col_num>
	mat<row_num, other_col_num, val_t> dot(const mat<col_num, other_col_num, val_t>& mt) const
	{
		return ::template dot<row_num, col_num, col_num, other_col_num, val_t>(*this, mt);
	}

	mat<row_num, col_num, val_t> rot180() const
	{
		mat<row_num, col_num, val_t> ret;
		for (int r = 0; r < row_num; ++r) 
		{
			for (int c = 0; c < col_num; ++c)
			{
				ret.get(r, c) = get(row_num-1-r, col_num-1-c);
			}
		}
		return ret;
	}

	template<int row_base, int col_base, int row_num_other, int col_num_other>
	void assign(const mat<row_num_other, col_num_other, val_t>& mt_other) 
	{
		/* 这里不麻烦了，直接写成运行时 */
		for (int r = 0; r < row_num_other; ++r) 
		{
			for (int c = 0; c < col_num_other; ++c) 
			{
				if (r + row_base < 0 || c + col_base < 0)
				{
					continue;
				}
				if (r + row_base >= row_num || c + col_base >= col_num) 
				{
					break;
				}
				get(r + row_base, c + col_base) = mt_other.get(r, c);
			}
		}
	}

	template<int top_pad, int left_pad, int right_pad, int bottom_pad>
	mat<row_num + top_pad + bottom_pad, col_num + left_pad + right_pad, val_t>
		pad() const
	{
		using mat_ret_t = mat<row_num + top_pad + bottom_pad, col_num + left_pad + right_pad, val_t>;
		mat_ret_t mt_ret;
		mt_ret.template assign<top_pad, left_pad>(*this);
		return mt_ret;
	}

	template<int row_span, int col_span>
	mat<row_num + row_span*(row_num - 1), col_num + col_span*(col_num-1)>
		span() const
	{
		using mat_ret_t = mat<row_num + row_span * (row_num - 1), col_num + col_span * (col_num - 1)>;
		mat_ret_t mt_ret;
		for (int r = 0; r < row_num; ++r) 
		{
			for (int c = 0; c < col_num; ++c) 
			{
				mt_ret.get(r*(row_span + 1), c*(col_span + 1)) = get(r, c);
			}
		}
		return mt_ret;
	}

	template<int row_base, int col_base, int row_len, int col_len>
	val_t region_max(int& i_row, int& i_col) const 
	{
		static_assert(row_base < row_num && col_base < col_num, "region_max overflow!!!");
		val_t d_max = -1. * DBL_MAX;
		for (int r = row_base; r < row_base + row_len && r < row_num; ++r) 
		{
			for (int c = col_base; c < col_base + col_len && c < col_num; ++c) 
			{
				if (d_max < get(r, c)) 
				{
					i_row = r, i_col = c;
					d_max = get(r, c);
				}
			}
		}
		return d_max;
	}

	mat<row_num*col_num, 1, val_t> one_col() const 
	{
		mat<row_num*col_num, 1, val_t> ret;
		ret.pval = pval;
		return ret;
	}

	template<typename t>
	static void print_sub_type(t) 
	{
		std::cout << typeid(t).name();
	}

	template<int r, int c, typename t>
	static void print_sub_type(mat<r, c, t>)
	{
		mat<r, c, t>::print_type();
	}

	static void print_type() 
	{
		printf("<matrix %d * %d type: ", row_num, col_num);
		print_sub_type<val_t>(val_t());
		printf(">\r\n");
	}

	const val_t& operator[](const int& idx) const
	{
		return pval->p[idx];
	}

	val_t& operator[](const int& idx)
	{
		return pval->p[idx];
	}

	friend std::ostream& operator<<(std::ostream& ofs, const mat<row_num, col_num, val_t>& mt)
	{
		std::cout << "[" ;
		for (int i = 0; i < row_num; ++i)
		{
			std::cout << std::setw(3) << "[";
			for (int j = 0; j < col_num; ++j)
			{
				std::cout << (j != 0 ? "," : "") << mt.get(i, j);
			}
			std::cout << std::setw(3) << "]";
		}
		std::cout << "]";
		return ofs;
	}

	int size() const
	{
		return row_num * col_num;
	}

	mat<row_num - 1, col_num - 1, val_t> algebraic_complement_val(const int& r, const int& c) const
	{
		mat<row_num - 1, col_num - 1, val_t> mt_ret;
		int retr = 0, retc = 0;
		for (int mr = 0; mr < row_num; mr++)
		{
			if (mr == r) continue;
			for (int mc = 0; mc < col_num; mc++)
			{
				if (mc == c) continue;
				mt_ret.get(retr, retc++) = get(mr, mc);
			}
			retc = 0;
			retr++;
		}
		return mt_ret;
	}

	mat<row_num, 1, val_t> col(const int& c) const 
	{
		mat<row_num, 1, val_t> ret;
		for (int i = 0; i < row_num; ++i) 
		{
			ret.get(i, 0) = this->get(i, c);
		}
		return ret;
	}

	mat<1, col_num, val_t> row(const int& r) const 
	{
		mat<1, col_num, val_t> ret;
		for (int i = 0; i < col_num; ++i) 
		{
			ret.get(0, i) = this->get(r, i);
		}
		return ret;
	}
};

template<typename type>
struct mat_size
{
	static constexpr int num = 1;
};

template<int target_col, int target_row, typename target_val_t>
struct mat_size<mat<target_col, target_row, target_val_t> >
{
	static constexpr int num = target_col * target_row * mat_size<target_val_t>::num;
};

template<typename unite_type>
struct mat_unite_type 
{
	using type = unite_type;
};

template<int row_num, int col_num, typename val_t>
struct mat_unite_type<mat<row_num, col_num, val_t> >
{
	using type = typename mat_unite_type<val_t>::type;
};

template<typename mat_t>
struct unite_mat
{
	using type = mat<mat_size<mat_t>::num, 1, typename mat_unite_type<mat_t>::type>;
};

template<typename mat_t>
using unite_mat_t = typename unite_mat<mat_t>::type;

template<typename target_type>
struct stretch_to_unite
{
	template<typename type1>
	static target_type cal(type1& v)
	{
		auto k = v.one_col();
		return stretch_to_unite<target_type>::cal(k);
	}

	static target_type cal(target_type& v)
	{
		return v;
	}
};

template<int i1, int i2, typename val_t>
mat<i1, i2, val_t> max_and_swap(const mat<i1, i2, val_t>& mt_max, const mat<i1, i2, val_t>& mt_optional)
{
	using ret_type = mat<i1, i2, val_t>;
	ret_type ret;
	for (int i = 0; i < i1; ++i) 
	{
		for (int j = 0; j < i2; ++j) 
		{
			ret.get(i, j) = max_and_swap(mt_max.get(i,j), mt_optional.get(i, j));
		}
	}
	return ret;
}

template<typename val_t>
val_t max_and_choose(const val_t& v1, const val_t& v2, const val_t& v3, const val_t& v4)
{
	return v1 < v2 ? v3 : v4;
}

template<int i1, int i2, typename val_t>
mat<i1, i2, val_t> max_and_choose(const mat<i1, i2, val_t>& mt_judge1, const mat<i1, i2, val_t>& mt_judge2, const mat<i1, i2, val_t>& mt_optional1, const mat<i1, i2, val_t>& mt_optional2)
{
	using ret_type = mat<i1, i2, val_t>;
	ret_type ret;
	for (int i = 0; i < i1; ++i)
	{
		for (int j = 0; j < i2; ++j)
		{
			ret.get(i, j) = max_and_choose(mt_judge1.get(i, j), mt_judge2.get(i, j), mt_optional1.get(i, j), mt_optional2.get(i, j));
		}
	}
	return ret;
}

template<int row_num, int col_num, typename val_t>
struct destoryer<mat<row_num, col_num, val_t> >
{
	using type = mat<row_num, col_num, val_t>;
	static void do_destory(type& v)
	{
		v.~type();
	}
};

#endif

矩阵基础运算

#ifndef _BASE_FUNCTION_HPP_
#define _BASE_FUNCTION_HPP_

#include <algorithm>
#include <math.h>
#include "base_logic.hpp"

// * 计算梯度
template<typename func_t, typename v_t>
auto derivative(func_t&& f, const v_t& v, const v_t& small_step = 0.001, const v_t& div_step = 0.002, const v_t& mul_step = 0.001)
{
	return (f(v + small_step) - f(v - small_step)) / div_step * mul_step;
}

/* 点乘运算 */
template<int r1, int c1, int r2, int c2, typename imatt1, typename imatt2, typename vt = double>
inline typename imatt1::type n_dot(const imatt1& mt1, const imatt2& mt2)
{
	static_assert(c1 == r2, "[matrix dot error]\tleft matrix column number do not match right matrix's row number.");
	if constexpr (c1 != 0 || r2 != 0)
	{
		return mt1.template get_val<r1, c1>() * mt2.template get_val<r2, c2>() + n_dot<r1, c1 - 1, r2 - 1, c2>(mt1, mt2);
	}
	if constexpr (c1 == 0 && r2 == 0)
	{
		return mt1.template get_val<r1, c1>() * mt2.template get_val<r2, c2>();
	}
}

template<int r, int c>
class v_dot
{
public:
	template<typename imatt1, typename imatt2, typename vt = double>
	static typename imatt1::type cal(const imatt1& mt1, const imatt2& mt2)
	{
		return n_dot<r, mt1.c - 1, mt2.r - 1, c>(mt1, mt2);
	}
};

template<int row_num1, int col_num1, int row_num2, int col_num2, typename val_t = double>
mat<row_num1, col_num2, val_t> dot(const mat<row_num1, col_num1, val_t>& mt1, const mat<row_num2, col_num2, val_t>& mt2)
{

	using omatt = mat<row_num1, col_num2, val_t>;
	using imatt1 = mat<row_num1, col_num1, val_t>;
	using imatt2 = mat<row_num2, col_num2, val_t>;
	omatt mt_ret;
	col_loop<col_num2 - 1, v_dot>(mt_ret, mt1, mt2);
	return mt_ret;
}

/* 加法运算 */
template<int r, int c>
class v_add
{
public:
	template<typename imatt1, typename imatt2>
	static typename imatt1::type cal(const imatt1& mt1, const imatt2& mt2)
	{
		return mt1.template get_val<r, c>() + mt2.template get_val<r, c>();
	}
};

template<int r, int c>
class n_add
{
public:
	template<typename imatt2, typename vt = double>
	static vt cal(const vt& mt1, const imatt2& mt2)
	{
		return mt1 + mt2.template get_val<r, c>();
	}
};

template<int r, int c>
struct c_add
{
	template<int row_num, int col_num, typename vt>
	static vt cal(const mat<row_num, col_num, vt>& mt, const mat<row_num, 1, vt>& v) 
	{
		return mt.template get_val<r, c>() + v.template get_val<r, 0>();
	}
};

template<int r, int c>
struct r_add
{
	template<int row_num, int col_num, typename vt>
	static vt cal(const mat<row_num, col_num, vt>& mt, const mat<1, col_num, vt>& v)
	{
		return mt.template get_val<r, c>() + v.template get_val<0, c>();
	}
};

template<int row_num, int col_num, typename val_t = double>
mat<row_num, col_num, val_t> operator+(const mat<row_num, col_num, val_t>& mt1, const mat<row_num, col_num, val_t>& mt2)
{
	using omatt = mat<row_num, col_num, val_t>;
	omatt mt_ret;
	col_loop<col_num - 1, v_add>(mt_ret, mt1, mt2);
	return mt_ret;
}

template<int row_num, int col_num, typename val_t = double>
mat<row_num, col_num, val_t> operator+(const val_t& v, const mat<row_num, col_num, val_t>& mt)
{
	using omatt = mat<row_num, col_num, val_t>;
	omatt mt_ret;
	col_loop<col_num - 1, n_add>(mt_ret, v, mt);
	return mt_ret;
}

template<int row_num, int col_num, typename val_t = double>
mat<row_num, col_num, val_t> operator+(const mat<row_num, col_num, val_t>& mt, const val_t& v)
{
	using omatt = mat<row_num, col_num, val_t>;
	omatt mt_ret;
	col_loop<col_num - 1, n_add>(mt_ret, v, mt);
	return mt_ret;
}

#if 0
template<int row_num, int col_num, typename val_t_other, typename val_t = double>
mat<row_num, col_num, val_t> operator+(const mat<row_num, col_num, val_t>& mt, const val_t_other& v)
{
	return mt + static_cast<val_t>(v);
}

template<int row_num, int col_num, typename val_t_other, typename val_t = double>
mat<row_num, col_num, val_t> operator+(const val_t_other& v, const mat<row_num, col_num, val_t>& mt)
{
	return mt + static_cast<val_t>(v);
}
#endif

template<int row_num, int col_num, typename val_t = double>
void spread_add(mat<row_num, col_num, val_t>& mt_ret, const mat<row_num, col_num, val_t>& mt, const mat<row_num, 1, val_t>& v)
{
	col_loop<col_num - 1, c_add>(mt_ret, mt, v);
}

template<int row_num, int col_num, typename val_t = double>
void spread_add(mat<row_num, col_num, val_t>& mt_ret, const mat<row_num, col_num, val_t>& mt, const mat<1, col_num, val_t>& v)
{
	col_loop<col_num - 1, r_add>(mt_ret, mt, v);
}

template<int row_num, int col_num, typename val_t = double>
void spread_add(mat<row_num, col_num, val_t>& mt_ret, const mat<row_num, col_num, val_t>& mt, const mat<1, 1, val_t>& v)
{
	col_loop<col_num - 1, n_add>(mt_ret, v.template get_val<0, 0>(), mt);
}

template<int row_num, int col_num, typename val_t = double>
void spread_add(mat<row_num, col_num, val_t>& mt_ret, const mat<row_num, 1, val_t>& v, const mat<row_num, col_num, val_t>& mt)
{
	col_loop<col_num - 1, c_add>(mt_ret, mt, v);
}

template<int row_num, int col_num, typename val_t = double>
void spread_add(mat<row_num, col_num, val_t>& mt_ret, const mat<1, col_num, val_t>& v, const mat<row_num, col_num, val_t>& mt)
{
	col_loop<col_num - 1, r_add>(mt_ret, mt, v);
}

template<int row_num, int col_num, typename val_t = double>
void spread_add(mat<row_num, col_num, val_t>& mt_ret, const mat<1, 1, val_t>& v, const mat<row_num, col_num, val_t>& mt)
{
	col_loop<col_num - 1, n_add>(mt_ret, v.template get_val<0, 0>(), mt);
}

template<typename val_t = double>
void spread_add(mat<1, 1, val_t>& mt_ret, const mat<1, 1, val_t>& v, const mat<1, 1, val_t>& mt)
{
	col_loop<0, n_add>(mt_ret, v.template get_val<0, 0>(), mt.template get_val<0, 0>);
}

/* 减法运算 */
template<int r, int c>
class v_minus
{
public:
	template<typename imatt1, typename imatt2, typename vt = double>
	static typename imatt1::type cal(const imatt1& mt1, const imatt2& mt2)
	{
		return mt1.template get_val<r, c>() - mt2.template get_val<r, c>();
	}
};

template<int r, int c>
class n_minus
{
public:
	template<int mat_r, int mat_c, typename vt = double>
	static vt cal(const vt& v, const mat<mat_r, mat_c, vt>& mt2)
	{
		return v - mt2.template get_val<r, c>();
	}

	template<int mat_r, int mat_c, typename vt = double>
	static vt cal(const mat<mat_r, mat_c, vt>& mt2, const vt& v)
	{
		return mt2.template get_val<r, c>() - v;
	}
};

template<int row_num, int col_num, typename val_t = double>
mat<row_num, col_num, val_t> operator-(const mat<row_num, col_num, val_t>& mt1, const mat<row_num, col_num, val_t>& mt2)
{
	using omatt = mat<row_num, col_num, val_t>;
	omatt mt_ret;
	col_loop<col_num - 1, v_minus>(mt_ret, mt1, mt2);
	return mt_ret;
}

template<int row_num, int col_num, typename val_t = double>
mat<row_num, col_num, val_t> operator-(const val_t& v, const mat<row_num, col_num, val_t>& mt)
{
	using omatt = mat<row_num, col_num, val_t>;
	omatt mt_ret;
	col_loop<col_num - 1, n_minus>(mt_ret, v, mt);
	return mt_ret;
}

template<int row_num, int col_num, typename val_t = double>
mat<row_num, col_num, val_t> operator-(const mat<row_num, col_num, val_t>& mt, const val_t& v)
{
	using omatt = mat<row_num, col_num, val_t>;
	omatt mt_ret;
	col_loop<col_num - 1, n_minus>(mt_ret, mt, v);
	return mt_ret;
}
/*
template<int row_num, int col_num, typename val_t_other, typename val_t = double>
mat<row_num, col_num, val_t> operator-(const mat<row_num, col_num, val_t>& mt, const val_t_other& v)
{
	return mt - val_t(v);
}

template<int row_num, int col_num, typename val_t_other, typename val_t = double>
mat<row_num, col_num, val_t> operator-(const val_t_other& v, const mat<row_num, col_num, val_t>& mt)
{
	return (val_t(v) - mt);
}
*/
/* 乘法运算 */
template<int r, int c>
class n_mul
{
public:
	template<typename imatt2>
	static typename imatt2::type cal(const typename imatt2::type& mt1, const imatt2& mt2)
	{
		return mt1 * mt2.template get_val<r, c>();
	}
};

template<int r, int c>
class v_mul
{
public:
	template<typename imatt1, typename imatt2>
	static typename imatt1::type cal(const imatt1& mt1, const imatt2& mt2)
	{
		return mt1.template get_val<r, c>() * mt2.template get_val<r, c>();
	}
};

template<int row_num, int col_num, typename val_t = double>
mat<row_num, col_num, val_t> operator*(const val_t& v, const mat<row_num, col_num, val_t>& mt)
{
	using omatt = mat<row_num, col_num, val_t>;
	omatt mt_ret;
	col_loop<col_num - 1, n_mul>(mt_ret, v, mt);
	return mt_ret;
}

template<int row_num, int col_num, typename val_t = double>
mat<row_num, col_num, val_t> operator*(const mat<row_num, col_num, val_t>& mt, const val_t& v)
{
	using omatt = mat<row_num, col_num, val_t>;
	omatt mt_ret;
	col_loop<col_num - 1, n_mul>(mt_ret, v, mt);
	return mt_ret;
}

template<int row_num, int col_num, typename val_t = double>
mat<row_num, col_num, val_t> operator*(const mat<row_num, col_num, val_t>& mt1, const mat<row_num, col_num, val_t>& mt2)
{
	using omatt = mat<row_num, col_num, val_t>;
	omatt mt_ret;
	col_loop<col_num - 1, v_mul>(mt_ret, mt1, mt2);
	return mt_ret;
}

/* 除法 */
template<int r, int c>
class n_div
{
public:
	template<int row_num, int col_num, typename vt = double>
	static vt cal(const mat<row_num, col_num, vt>& mt, const vt& v)
	{
		return mt.template get_val<r, c>() / v;
	}

	template<int row_num, int col_num, typename vt = double>
	static vt cal(const vt& v, const mat<row_num, col_num, vt>& mt)
	{
		return v / mt.template get_val<r, c>();
	}
};

template<int r, int c>
class v_div
{
public:
	template<int row_num, int col_num, typename vt = double>
	static vt cal(const mat<row_num, col_num, vt>& mt1, const mat<row_num, col_num, vt>& mt2)
	{
		return mt1.template get_val<r, c>() / mt2.template get_val<r, c>();
	}
};

template<int row_num, int col_num, typename val_t>
mat<row_num, col_num, val_t> operator/(const mat<row_num, col_num, val_t>& mt, const val_t& v)
{
	using omatt = mat<row_num, col_num, val_t>;
	omatt mt_ret;
	col_loop<col_num - 1, n_div>(mt_ret, mt, v);
	return mt_ret;
}

template<int row_num, int col_num, typename val_t>
mat<row_num, col_num, val_t> operator/(const val_t& v, const mat<row_num, col_num, val_t>& mt)
{
	using omatt = mat<row_num, col_num, val_t>;
	omatt mt_ret;
	col_loop<col_num - 1, n_div>(mt_ret, v, mt);
	return mt_ret;
}

template<int row_num, int col_num, typename vt>
mat<row_num, col_num, vt> operator/(const mat<row_num, col_num, vt>& mt1, const mat<row_num, col_num, vt>& mt2)
{
	using omatt = mat<row_num, col_num, vt>;
	omatt mt_ret;
	col_loop<col_num - 1, v_div>(mt_ret, mt1, mt2);
	return mt_ret;
}

double sqrtl(double d)
{
    return sqrt(d);
}

template<int i1, int i2, typename val_t>
mat<i1, i2, val_t> sqrtl(const mat<i1, i2, val_t>& mt)
{
	using type = mat<i1, i2, val_t>;
	type ret;
	for (int i = 0; i < i1; ++i) 
	{
		for (int j = 0; j < i2; ++j) 
		{
			ret.get(i, j) = sqrtl(mt.get(i,j));
		}
	}
	return ret;
}

template<int r, int c>
class n_sqrt
{
public:
	template<typename imatt>
	static typename imatt::type cal(const imatt& mt)
	{
		return sqrtl(mt.template get_val<r, c>());
	}
};

template<int row_num, int col_num, typename val_t = double>
mat<row_num, col_num, val_t> sqrtm(const mat<row_num, col_num, val_t>& mt)
{
	using omatt = mat<row_num, col_num, val_t>;
	omatt mt_ret;
	col_loop<col_num - 1, n_sqrt>(mt_ret, mt);
	return mt_ret;
}

/* exp运算 */
template<int r, int c>
struct n_exp 
{
	template<typename imatt>
	static typename imatt::type cal(const imatt& mt ) 
	{
		return exp(mt.template get_val<r, c>());
	}
};

template<int row_num, int col_num, typename val_t = double>
mat<row_num, col_num, val_t> exp(const mat<row_num, col_num, val_t>& mt)
{
	using omatt = mat<row_num, col_num, val_t>;
	omatt mt_ret;
	col_loop<col_num - 1, n_exp>(mt_ret, mt);
	return mt_ret;
}

/* 卷积运算 */
template<int row_base, int col_base, int row_delta, int col_delta, typename imat_origin, typename imat_tpl>
inline auto col_loop_mul(const imat_origin& mt_origin, const imat_tpl& mt_tpl)
{
	if constexpr (col_delta != 0)
	{
		return mt_origin.template get_val<row_base + row_delta, col_base + col_delta>() * mt_tpl.template get_val<row_delta, col_delta>()
			+ col_loop_mul<row_base, col_base, row_delta, col_delta - 1, imat_origin, imat_tpl>(mt_origin, mt_tpl);
	}
	if constexpr (col_delta == 0)
	{
		return mt_origin.template get_val<row_base + row_delta, col_base + col_delta>() * mt_tpl.template get_val<row_delta, col_delta>();
	}
}

template<int row_base, int col_base, int row_delta, int col_delta, typename imat_origin, typename imat_tpl>
inline auto row_loop_add(const imat_origin& mt_origin, const imat_tpl& mt_tpl)
{
	if constexpr (row_delta != 0)
	{
		return col_loop_mul<row_base, col_base, row_delta, col_delta>(mt_origin, mt_tpl)
			+ col_loop_mul<row_base, col_base, row_delta - 1, col_delta>(mt_origin, mt_tpl);
	}
	if constexpr (row_delta == 0)
	{
		return col_loop_mul<row_base, col_base, row_delta, col_delta>(mt_origin, mt_tpl);
	}
}

template<int r, int c>
struct v_inner_conv
{
	template<typename imat_origin_t, typename imat_tpl_t>
	inline static auto cal(const imat_origin_t& mt_origin, const imat_tpl_t& mt_tpl)
	{
		return row_loop_add<r, c, imat_tpl_t::r - 1, imat_tpl_t::c - 1>(mt_origin, mt_tpl);
	}
};

constexpr int get_step_inner_size(int i_origin, int i_tpl, int i_step)
{
	return (i_origin - i_tpl) / i_step + 1;
}
class def_config;
class same_pad;

template<typename op_t = def_config>
constexpr int get_pad_size(int i_origin, int i_tpl, int i_step) 
{
	return (((i_origin - i_tpl) / i_step) + (((i_origin - i_tpl) % i_step) == 0 ? 0 : 1))*i_step - (i_origin - i_tpl);
}

template<>
constexpr int get_pad_size<same_pad>(int i_origin, int i_tpl, int i_step)
{
	//static_assert(i_step == 1, "ERROR:get_pad_size same pad step must equal 1.");
	return i_tpl - 1;			// a + x - b + 1 = a; x = b-1;
}

constexpr int get_ceil_div(int i_origin, int i_tpl)
{
	return (i_origin / i_tpl + ((i_origin%i_tpl)==0?0:1));
}

template<int input_row, int intput_col, int tpl_row, int tpl_col, int row_step, int col_step, typename op_t = def_config>
struct pad_size_t
{
	static constexpr int top = get_pad_size<op_t>(input_row, tpl_row, row_step) / 2;				// 上方的扩充数量
	static constexpr int left = get_pad_size<op_t>(intput_col, tpl_col, col_step) / 2;			// 左方的扩充数量
	static constexpr int right = get_pad_size<op_t>(intput_col, tpl_col, col_step) - left;		// 右方的扩充数量
	static constexpr int bottom = get_pad_size<op_t>(input_row, tpl_row, row_step) - top;			// 下方的扩充数量
};

// * 内部卷积
template<int row_step, int col_step, int row_num, int col_num, int tpl_row, int tpl_col, typename val_t>
inline mat<get_step_inner_size(row_num, tpl_row, row_step), get_step_inner_size(col_num, tpl_col, col_step), val_t>
inner_conv(const mat<row_num, col_num, val_t>& mt_origin, const mat<tpl_row, tpl_col, val_t>& mt_tpl)
{
	using ret_type = mat<get_step_inner_size(row_num, tpl_row, row_step), get_step_inner_size(col_num, tpl_col, col_step), val_t>;
	ret_type mt_ret;
	col_loop<ret_type::c - 1, v_inner_conv>(mt_ret, mt_origin, mt_tpl);
	return mt_ret;
}

template<typename mat_t, typename ...mat_ts>
struct st_one_col 
{
	static constexpr int all_size = (mat_t::r*mat_t::c) + st_one_col<mat_ts...>::all_size;
};

template<typename mat_t>
struct st_one_col<mat_t>
{
	static constexpr int all_size = (mat_t::r*mat_t::c);
};

template<typename mat_t, typename ...mat_ts>
void concat_mat(typename mat_t::type* p, const mat_t& mt, const mat_ts... mts) 
{
	constexpr int cpy_size = mat_t::r*mat_t::c;
	//memcpy(p, mt.pval->p, cpy_size*sizeof(mat_t::type));
	for (int i = 0; i < cpy_size; ++i) 
	{
		p[i] = mt.pval->p[i];
	}
	if constexpr(0!=sizeof...(mat_ts))
		concat_mat(p + cpy_size, mts...);
}

template<typename mat_t, typename ...mat_ts>
mat<st_one_col<mat_t, mat_ts...>::all_size, 1> stretch_one_col(const mat_t& mt, const mat_ts&...mts)
{
	using ret_type = mat<st_one_col<mat_t, mat_ts...>::all_size, 1>;
	ret_type ret;
	concat_mat(ret.pval->p, mt, mts...);
	return ret;
}


template<typename mat_t, typename ...mat_ts>
void split_mat(typename mat_t::type* p, const mat_t& mt, const mat_ts... mts)
{
	constexpr int cpy_size = mat_t::r*mat_t::c;
	//memcpy(mt.pval->p, p, cpy_size * sizeof(mat_t::type));
	for (int i = 0; i < cpy_size; ++i)
	{
		mt.pval->p[i] = p[i];
	}
	if constexpr (0 != sizeof...(mat_ts))
		split_mat(p + cpy_size, mts...);
}

template<typename mat_t, typename ...mat_ts>
void split_one_mat(const mat_t& mt, const mat_ts&...mts)
{
	split_mat(mt.pval->p, mts...);
}

template<int i1, int i2, typename val_t>
mat<i1, i2, val_t> abs(const mat<i1, i2, val_t>& mt)
{
	mat<i1, i2, val_t> ret;
	for (int i = 0; i < i1; ++i) 
	{
		for (int j = 0; j < i2; ++j) 
		{
			ret.get(i, j) = abs(mt.get(i, j));
		}
	}
	return ret;
}

#include <vector>
/* 输出改成0均值1均方差的 */
template<int row_num, int col_num, typename val_t>
inline std::vector<mat<row_num, col_num, val_t> > normalize(const std::vector<mat<row_num, col_num, val_t> >& vec_input, mat<row_num, col_num, val_t>& mt_mean, mat<row_num, col_num, val_t>& mt_div)
{
	if (1 == vec_input.size())return vec_input;
	using type = mat<row_num, col_num, val_t>;
	std::vector<type> vec_ret;
	//type mt_mean, mt_div;
	for (int i = 0; i < vec_input.size(); ++i)
	{
		mt_mean = mt_mean + vec_input[i] / static_cast<val_t>(vec_input.size());
	}

	for (int i = 0; i < vec_input.size(); ++i)
	{
		auto delta = vec_input[i] - mt_mean;
		mt_div = mt_div + delta * delta / static_cast<val_t>(vec_input.size());
	}

	auto mt_s = sqrtl(mt_div);

	for (int i = 0; i < vec_input.size(); ++i)
	{
		vec_ret.push_back((vec_input[i] - mt_mean) / mt_s);
	}
	return vec_ret;
}

#include "ht_memory.h"

template<typename val_t>
void write_file(const val_t& vt, ht_memory& mry) 
{
	mry << vt;
}

template<int row_num, int col_num, typename val_t>
void write_file(const mat<row_num, col_num, val_t>& mt, ht_memory& mry)
{
	for (int r = 0; r < row_num; ++r)
	{
		for (int c = 0; c < col_num; ++c)
		{
			write_file(mt.get(r, c), mry);
		}
	}
}

template<typename val_t>
void read_file(ht_memory& mry, val_t& vt) 
{
	mry >> vt;
}

template<int row_num, int col_num, typename val_t>
void read_file(ht_memory& mry, mat<row_num, col_num, val_t>& mt)
{
	for (int r = 0; r < row_num; ++r)
	{
		for (int c = 0; c < col_num; ++c)
		{
			read_file(mry, mt.get(r, c));
		}
	}
}

template<int N>
double sub_mul_cal(const mat<N, N, double>& mt, const std::vector<int>& vec_cur)
{
	double dmul = 1.;
	for (size_t siz_add_itr = 0; siz_add_itr < vec_cur.size(); ++siz_add_itr)
	{
		dmul *= mt.get(siz_add_itr, vec_cur[siz_add_itr]);
		//printf("%4.2lf ", mt.get(siz_add_itr, vec_cur[siz_add_itr]));
	}
	//printf("=%4.2lf", dmul);
	return dmul;
}

template<int N>
void det_cal(double& dret, const mat<N, N, double>& mt, const std::vector<int>& vec_idx, const size_t& siz_cur, const double& dflag = 1.)
{
	//printf("\n%4.2lf ", dflag);
	dret += (sub_mul_cal(mt, vec_idx)*dflag);			// * 把自己算进去
	for (size_t siz_itr = siz_cur; siz_itr < vec_idx.size() - 1 ; ++siz_itr)
	{
		std::vector<int> vec_cur(vec_idx);
		double dcurflag = dflag;
		for (size_t siz_in_itr = siz_itr + 1; siz_in_itr < vec_cur.size(); ++siz_in_itr)
		{
			std::swap(vec_cur[siz_itr], vec_cur[siz_in_itr]);
			dcurflag *= -1.;
			//dret += (sub_mul_cal(mt, vec_cur)*dflag);
			det_cal(dret, mt, vec_cur, siz_itr + 1, dcurflag);
		}
	}
}

template<int N>
double det(const mat<N, N, double>& mt)
{
	std::vector<int> vec(N, 0);
	int idx = 0;
	std::generate(vec.begin(), vec.end(), [&idx]()
	{
		return idx++;
	});
	double dret = 0.;
	det_cal(dret, mt, vec, 0, 1.);
	return dret;
}

template<int N, typename val_t>
mat<N, N, val_t> algebraic_complement(const mat<N, N, val_t>& mt)
{
	mat<N, N, val_t> mtret;
	double drflag = 1.;
	for (int i = 0; i < N; ++i)
	{
		double dcflag = 1.;
		for (int j = 0; j < N; ++j)
		{
			mtret.get(i, j) = (drflag * dcflag * det(mt.algebraic_complement_val(i, j)));
			dcflag *= -1.;
		}
		drflag *= -1.;
	}
	return mtret.t();
}

template<int N, typename val_t>
mat<N, N, val_t> inverse(const mat<N, N, val_t>& mt)
{
	return algebraic_complement(mt)/det(mt);
}

#endif

试验程序

#include "mat.hpp"
#include "base_function.hpp"


int main(int argc, char** argv) 
{
	mat<4, 3, double> x = {1., 2., 1., 2., 1., 1., 3., 1., 1., 3., 1., 1.};
	mat<4, 1, double> y = {12., 9.5, 11.4, 11.6};
	auto w = inverse(x.t().dot(x)).dot(x.t()).dot(y);		// W=inv(X.t·X)·X.t·Y
	w.print();
	x.dot(w).print();
	return 0;
}

运行结果如下：

由于并非平面，所以展现出来的就是最后两个值线性拟合在中间位置。

用梯度下降法解决

#include "mat.hpp"
#include "base_function.hpp"
#include <functional>

template<typename mt_t>
double regression_func(const mt_t& x, const decltype(x.t())& w) 
{
	return exp(x.dot(w)[0]);
}

template<typename mt_t>
double error_func(double y, const mt_t& x, const decltype(x.t())& w)
{
	double y_ = regression_func(x, w);
	return (y - y_)*(y - y_);
}


int main(int argc, char** argv) 
{
	mat<3, 3, double> x = { 1., 2., 1., 2., 1., 1.};
	mat<3, 1, double> y = { 12., 9.5 };
	mat<3, 1, double> w = {0, 0, 0};
	using namespace std;
	for (int i = 0; i < 200; ++i) 
	{
		for (int j = 0; j < 2; ++j)
		{
			function<double(const mat<3, 1, double>&)> wgrad_f = std::bind(error_func<mat<1, 3, double>>, y[j], x.row(j), std::placeholders::_1);
			
			mat<3, 1, double> div_step = 0.0002;
			for (int k = 0; k < div_step.size(); ++k)
			{
				mat<3, 1, double> step, mul_step;
				step[k] = 0.0001;
				mul_step[k] = 1.;
				mat<3, 1, double> delta = derivative(wgrad_f, w, step, div_step, mul_step);
				w = w - 0.001*delta;
			}
		}		
	}
	w.print();
	printf("%lf, %lf\r\n", regression_func(x.row(0), w), regression_func(x.row(1), w));
	return 0;
}

运行结果如下：