原理基于 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;
}
测试结果：
当然代码也还有可以优化的地方，以及这些都是争对方阵而编写的一段矩阵、向量运算相关的简陋代码。
