矩阵求逆（LUP分解算法）

　　原理基于 gaussian jordan elimination 方法，考虑求解如下线性方程组：

$$ egin{aligned} mathbf{A}mathbf{x} &= mathbf{b} end{aligned} $$

　　我们设 $ L $、$ U $、$ P $ 满足：

$$ egin{aligned} mathbf{PA} &= mathbf{LU} end{aligned} $$

　　所以只要找到这样的 $ L $、$ U $、$ P $， 就可以通过：

$$ egin{aligned} mathbf{Ax} &= mathbf{b} \ mathbf{PAx} &= mathbf{Pb} \ mathbf{LUx} &= mathbf{Pb} \ mathbf{Ux} &= mathbf{y} \ mathbf{Ly} &= mathbf{Pb} end{aligned} $$

　　求出 $ x $ 了，注意这里 $ x $ 是一个列向量。所以该方法可以用来求矩阵的逆。

　　根据矩阵与其逆矩阵的关系 $ mathbf{AA^{-1} = I} $， $ mathbf{I} $ 是一个单位矩阵，所以只需要把每行（或列）作为 $ mathbf{A}mathbf{x} = mathbf{b} $ 中的 $ mathbf{b} $ 就可以求出逆矩阵$ mathbf{{A^{-1}}^{T}} $ 的对应行的向量元素，然后转置即可（实际上不用转置，这个是由于没有实现一个自己用的矩阵库的原因。），也可以考虑使用 opencv 自带的矩阵库的 inv 方法或者 boost 的矩阵库实现。

　　这里 $ L $、$ U $、$ P $ 都是我们在做 gaussian jordan elimination 方法时交换行所产生的矩阵（代码中 ```lup_decomposition``` 函数中实现了该过程。）

　　代码实现：

#include <vector>
#include <cassert>
#include <ostream>
#include <iomanip>
#include <iostream>

using std::vector;
using std::ostream;
using std::cout;
using std::swap;

bool lup_decomposition(vector<vector<double>> a, vector<vector<double>>& l, vector<vector<double>>& u, vector<double>& p)
{
	for (int i = 0; i < a.size(); i++)
		p[i] = (double)i;
	for (int i = 0; i < a.size(); i++)
	{
		l[i][i] = 1;
		for (int j = i + 1; j < a.size(); j++)
			l[i][j] = 0;
	}
	for (int i = 1; i < a.size(); i++)
		for (int j = 0; j < i; j++)
			u[i][j] = 0;
	for (int i = 0; i < a.size(); i++)
	{
		double n = 0; int i_;
		for (int j = i; j < a.size(); j++)
			if (abs(a[j][i]) > n)
				n = abs(a[j][i]), i_ = j;
		if (n == 0) return false;
		swap(p[i], p[i_]);
		for (int j = 0; j < a.size(); j++)
			swap(a[i][j], a[i_][j]);
		u[i][i] = a[i][i];
		for (int j = i + 1; j < a.size(); j++)
		{
			l[j][i] = a[j][i] / u[i][i];
			u[i][j] = a[i][j];
		}
		for (int j = i + 1; j < a.size(); j++)
			for (int k = i + 1; k < a.size(); k++)
				a[j][k] -= l[j][i] * u[i][k];
	}
	return true;
}

void lup_solve(vector<vector<double>> a, vector<double> b, vector<vector<double>> l, vector<vector<double>> u, vector<double> p, vector<double>& x)
{
	vector<double> y(a.size(), 0);
	for (int i = 0; i < a.size(); i++)
	{
		double sum = 0;
		for (int j = 0; j < i; j++)
			sum += l[i][j] * y[j];
		y[i] = b[p[i]] - sum;
	}
	for (int i = a.size() - 1; i >= 0; i--)
	{
		double sum = 0;
		for (int j = i + 1; j < a.size(); j++)
			sum += u[i][j] * x[j];
		x[i] = (y[i] - sum) / u[i][i];
	}
}

vector<vector<double>> inverse(vector<vector<double>> a)
{
	vector<vector<double>>
		unit_matrix(a.size(), vector<double>(a.size(), 0)),
		x(a.size(), vector<double>(a.size(), 0)),
		l(a.size(), vector<double>(a.size(), 0)),
		u(a.size(), vector<double>(a.size(), 0));
	vector<double> p(a.size());
	for (int i = 0; i < a.size(); i++)
		unit_matrix[i][i] = 1;
	lup_decomposition(a, l, u, p);
	for (int i = 0; i < a.size(); i++)
		lup_solve(a, unit_matrix[i], l, u, p, x[i]);
	return x;
}

vector<double> sum(vector<vector<double>> a, int axis = 1)
{
	if (axis == 1)
	{
		vector<double> e(a.size());
		for (int i = 0; i < a.size(); i++)
			for (int j = 0; j < a[0].size(); j++)
				e[i] += a[i][j];
		return e;
	}
	else if (axis == -1)
	{
		vector<double> e(a[0].size());
		for (int i = 0; i < a[0].size(); i++)
			for (int j = 0; j < a.size(); j++)
				e[i] += a[i][j];
		return e;
	}
	else {
		vector<double> e(1);
		for (int i = 0; i < a.size(); i++)
			for (int j = 0; j < a[0].size(); j++)
				e[0] += a[i][j];
		return e;
	}
}

double sum(vector<double> a)
{
	double e = 0;
	for (auto i : a)
		e += i;
	return e;
}

vector<double> add(vector<double> a, vector<double> b)
{
	assert(a.size() == b.size());
	vector<double> c(a.size());
	for (int i = 0; i < a.size(); i++)
		c[i] = a[i] + b[i];
	return c;
}

vector<vector<double>> add(vector<vector<double>> a, vector<vector<double>> b)
{
	assert(a.size() == b.size());
	assert(a[0].size() == b[0].size());
	vector<vector<double>> c(a.size(), vector<double>(a[0].size(), 0));
	for (int i = 0; i < a.size(); i++)
		for (int j = 0; j < a[0].size(); j++)
			c[i][j] = a[i][j] + b[i][j];
	return c;
}

vector<double> sub(vector<double> a, vector<double> b)
{
	assert(a.size() == b.size());
	vector<double> c(a.size());
	for (int i = 0; i < a.size(); i++)
		c[i] = a[i] - b[i];
	return c;
}

vector<vector<double>> sub(vector<vector<double>> a, vector<vector<double>> b)
{
	assert(a.size() == b.size());
	assert(a[0].size() == b[0].size());
	vector<vector<double>> c(a.size(), vector<double>(a[0].size(), 0));
	for (int i = 0; i < a.size(); i++)
		for (int j = 0; j < a[0].size(); j++)
			c[i][j] = a[i][j] - b[i][j];
	return c;
}

vector<double> mul(vector<vector<double>> a, vector<double> b)
{
	assert(a[0].size() == b.size());
	vector<double> c(b.size(), 0);
	for (int i = 0; i < a.size(); i++)
		for (int j = 0; j < a[0].size(); j++)
			c[i] += a[i][j] * b[j];
	return c;
}

vector<vector<double>> mul(vector<vector<double>> a, vector<vector<double>> b)
{
	assert(a[0].size() == b.size());
	vector<vector<double>> r(a.size(), vector<double>(b[0].size(), 0));
	for (int i = 0; i < r.size(); i++)
		for (int j = 0; j < r[0].size(); j++)
			for (int k = 0; k < b.size(); k++)
				r[i][j] += a[i][k] * b[k][j];
	return r;
}

double dot(vector<double> a, vector<double> b)
{
	assert(a.size() == b.size());
	double c = 0.;
	for (int i = 0; i < a.size(); i++)
		c += a[i] * b[i];
	return c;
}

vector<vector<double>> transpose(vector<vector<double>> a)
{
	for (int i = 0; i < a.size(); i++)
		for (int j = i + 1; j < a[0].size(); j++)
			swap(a[i][j], a[j][i]);
	return a;
}

ostream& operator<<(ostream& os, vector<vector<double>> a)
{
	os << "[
";
	for (auto rows : a)
	{
		for (auto col : rows)
			os << std::right << std::setw(12) << col << ' ';
		os << '
';
	}
	os << "]
";
	return os;
}

ostream& operator<<(ostream& os, vector<double> a)
{
	for (auto i : a) os << i << '
';
	return os;
}

int main()
{
	vector<vector<double>> m{ {1, 2}, {2, 3} };

	cout << m << '
';
	cout << transpose(inverse(m)) << '
';
	
	return 0;
}

　　测试结果：

　　当然代码也还有可以优化的地方，以及这些都是争对方阵而编写的一段矩阵、向量运算相关的简陋代码。