|
|
Fortran或者matlab老手在用C++以后常会怀恋矩阵的各类操作：包括但不限于实复矩阵转换，切片向量化。遗憾是C++并没有内建矩阵数据类型，因此常借助Eigen/gsl等三方库完成这些功能。但不知为何，这些库语法也不像Fortran/Python/Matlab看齐，用起来十分不爽。
既然这样，不如自己设计一个矩阵库，让它尽量接近Matlab/Python/Fortran的数值计算语法。
基础设计
设计矩阵库的第一步是矩阵类本身。这里我不打算遵循Java/C++的惯例，让线代算法成为矩阵类的成员方法，而是参考C/python的设计，让部分算法成为全局函数。矩阵本身数据成员只有数据本身，它看起来就像这样：
template <std::size_t M, std::size_t N>
class Matrix
{
template <typename T, typename RT = void>
using enable_arith_type_t = typename std::enable_if<std::is_arithmetic<T>::value, RT>::type;
public:
// Constructors
Matrix() : m_data(M * N, 0.0){};
~Matrix() = default;
Matrix(const Matrix &) = default;
Matrix(Matrix &&) = default;
Matrix &operator=(const Matrix &) = default;
Matrix &operator=(Matrix &&) = default;
template <typename T, size_t L, enable_arith_type_t<T> * = nullptr>
Matrix(const std::array<T, L> &list);
template <typename T, enable_arith_type_t<T> * = nullptr>
Matrix(const std::initializer_list<T> &list);
template <typename T, enable_arith_type_t<T> * = nullptr>
Matrix(const std::vector<T> &list);
// Member functions
double *data();
const double *data() const;
std::vector<double> &flat();
const std::vector<double> &flat() const;
std::array<double, N> row(size_t idx) const;
std::array<double, M> col(size_t idx) const;
double &operator()(size_t row, size_t col);
const double &operator()(size_t row, size_t col) const;
void fill(double val);
Matrix<N, M> T() const;
// Overloaded Operators
Matrix operator+(const Matrix &other) const;
Matrix operator-(const Matrix &other) const;
template <size_t L>
Matrix<M, L> operator*(const Matrix<N, L> &other) const;
std::array<double, M> operator*(const std::array<double, N> &x) const;
...
template <typename T>
enable_arith_type_t<T, Matrix<M, N>> operator*(T ele);
template <typename T>
enable_arith_type_t<T, Matrix<M, N>> operator/(T ele);
...
// Generate function.
friend std::ostream &operator<<(std::ostream &os, const Matrix &self);
static Matrix eye();
static size_t row_counts();
static size_t col_counts();
private:
// data member
std::vector<double> m_data;
};
我们使用非类型模板参数在构造时候约束矩阵的大小，选择vector作为矩阵的数据容器。这样做的理由是：实际使用中矩阵的维数往往是预先知道的信息，并且编译期约束了大小，可以检查一些矩阵运算的错误。而vector相比与double*和double**的好处在于省去了数据成员的生命周期管理，并且一维数组的访问效率有保障。不使用array在于其在栈上分配内存，矩阵维度稍大容易在嵌入式设备中爆栈。
这是一个基础的矩阵类接口设计。优点是矩阵对象实质上是扁平化数组，操作方便。缺点在于：无数组按行索引或者按列索引，总有相应的行或列返回的是不连续内存，这使得行列的访问接口不一致。这也是为什么C++矩阵库按行或者按列赋值总要用奇怪的函数，没法像PY或者Fortran一样直接拿到子矩阵索引的原因。
操作代理
上述问题并非无法解决。其它语言储存矩阵的底层也是数组。但是对它的切片进行操作时应该是通过某种代理转化到了数组上的操作。在C++中，我们可以模仿该操作，设计访问代理，实现切片访问功能。
要模仿的语法是这样的：
A(:,1)=C(:,2);
A(1:2,2:3)=[3.3,4.4;5.5,6.6];首先需要重载括号操作符，用pair{a,b}代替a:b作为参数。接受到该类型参数时视为输入了一个切片范围，然后并不返回元素引用，而是返回我们的代理类。
代理类设计如下，我把它放到了Matrix的私有部分：
struct SubPart
{
public:
SubPart(Matrix<M, N> &self, IndexRange r, IndexRange c)
: row_idx(0), row_end(0), col_idx(0), col_end(0), data(self)
{
row_idx = r.first >= 0 ? r.first : 0;
col_idx = c.first >= 0 ? c.first : 0;
row_end = r.second >= 0 ? r.second : M - 1;
col_end = c.second >= 0 ? c.second : N - 1;
}
SubPart(const SubPart &) = delete;
SubPart(SubPart &&) = delete;
SubPart &operator=(const SubPart &) = delete;
SubPart &operator=(SubPart &&) = delete;
template <size_t A, size_t B>
void operator=(const Matrix<A, B> &other);
template <typename T, size_t A, enable_arith_type_t<T> * = nullptr>
void operator=(const std::array<T, A> &list);
template <typename T, enable_arith_type_t<T> * = nullptr>
void operator=(const std::initializer_list<T> &list);
template <typename T, enable_arith_type_t<T> * = nullptr>
void operator=(const std::vector<T> &list);
std::array<double, M> operator*(const std::array<double, N> &x) const;
template <size_t L>
Matrix<M, L> operator*(const Matrix<N, L> &other) const;
private:
size_t row_idx;
size_t col_idx;
size_t row_end;
size_t col_end;
Matrix<M, N> &data;
template <typename T>
void generator_by_list(const T &list);
};
在矩阵中调用括号操作符并且输入是范围时，将返回SubPart对象，继续调用对SubPart的函数就完成了对矩阵本身操作（我只实现了赋值操作符，其它同理），这样完成了切片的引用语义。
切片的值语义通过全局函数slice实现：
template <size_t A, size_t B, size_t M, size_t N>
Matrix<A, B> slice(const Matrix<M, N> &m, size_t row_start, size_t col_start)有了这两个功能以后，矩阵类就可以完成其它语法一般的切片操作了：
Matrix<4, 4> a = {1, 2, 3, 4, 5, 6, 7, 8};
PRINT_SINGLE_ELEMENTS(a);
a({2, 3}, {2, 3}) = Matrix<2, 2>{1, 1, 1, 1};
a( 0, {0, -1}) = {89.0, 23.0, 44.0, 9.8};
PRINT_SINGLE_ELEMENTS(a, &#34;a = &#34;);
==================================================
Matrix<4,4>:
1 5 0 0
2 6 0 0
3 7 0 0
4 8 0 0
a = Matrix<4,4>:
89 23 44 9.8
2 6 0 0
3 7 1 1
4 8 1 1乘法、行列式与逆
乘法
稠密矩阵乘法叫做GEMM，是专门的领域。简单地说：想要发挥乘法的最大性能，需要硬件逻辑配合，针对专门的平台提供相应的汇编指令或者开发接口，然后再数据上做分块、对齐、喂给并行指令。这样做绑定了硬件，如X86的AVX，CUDA的cuBLAs等等，没有移植性。感兴趣的可以通过下面链接入门：
抛开汇编指令层，矩阵乘法实际上是个三重积： c_{ij}=a_{ik}b_{kj} 。最朴素实现：
for (size_t i = 0; i < M; i++)
{
for (size_t j = 0; j < L; j++)
{
for (size_t k = 0; k < N; k++)
{
result(i, j) += (*this)(i, k) * other(k, j);
}
}
}所有优化都是参考上述进行对比的。我这里简单的取个巧，将优化抛给STL底层库：
Matrix<N, M> AT = this->T();
auto iter_A = AT.data();
auto iter_B = other.data();
for (size_t i = 0; i < M; i++)
{
iter_B = other.data();
for (size_t j = 0; j < L; j++)
{
result(i, j) = std::inner_product(iter_A, iter_A + N, iter_B, 0.0);
iter_B += N;
}
iter_A += N;
}这么做浪费了一倍的空间。好处在于，标准库的内积可以在两块连续的内存上滑动。具体优化就看各家编译器实现。
在MSVC和GCC上benchmark，测500x500稠密矩阵乘法（引入了Eigen作为参考）：
优化以前:
2022-11-13T13:47:05+08:00
Running C:\Users\Desktop\matrix\build\Release\matrix_bench.exe
Run on (12 X 1616.6 MHz CPU s)
CPU Caches:
L1 Data 32 KiB (x6)
L1 Instruction 32 KiB (x6)
L2 Unified 256 KiB (x6)
L3 Unified 12288 KiB (x1)
------------------------------------------------------------
Benchmark Time CPU Iterations
------------------------------------------------------------
BM_MatrixMul 7126967 ns 7118056 ns 90
BM_MatrixMulEigen 1058715 ns 1074219 ns 640
优化以后:
2022-11-13T13:45:23+08:00
Running C:\Users\Desktop\matrix\build\Release\matrix_bench.exe
Run on (12 X 1608.38 MHz CPU s)
CPU Caches:
L1 Data 32 KiB (x6)
L1 Instruction 32 KiB (x6)
L2 Unified 256 KiB (x6)
L3 Unified 12288 KiB (x1)
------------------------------------------------------------
Benchmark Time CPU Iterations
------------------------------------------------------------
BM_MatrixMul 2188033 ns 2148438 ns 320
BM_MatrixMulEigen 1071652 ns 1035392 ns 498标准的三重循环性能大概只有Eigen乘法的1/7，简单的用空间换时间以后大约是Eigen的1/2。
更多性能提升需要对齐数据做掩模调用平台汇编手动分配寄存器，非常耗费人工。
行列式
行列式的定义是从代数余子式来的，因此一种求法是递归的求解代数余子式：
代码大概是这样，先求余子式：
template <size_t M, size_t N>
Matrix<M - 1u, N - 1u> cofactor(const Matrix<M, N> &mat, size_t p, size_t q)
{
Matrix<M - 1u, N - 1u> result{};
size_t i = 0, j = 0;
for (size_t row = 0; row < M; row++)
{
for (size_t col = 0; col < N; col++)
{
if (row != p && col != q)
{
result(i, j++) = mat(row, col);
if (j == N - 1)
{
j = 0;
i++;
}
}
}
}
return result;
}然后递归求伴随矩阵与行列式：
template <size_t M>
double determinant(const Matrix<M, M> &mat)
{
double D = 0;
int sign = 1;
for (size_t f = 0; f < M; f++)
{
D += sign * mat(0, f) * determinant(cofactor(mat, 0, f));
sign = -sign;
}
return D;
}
template <>
double determinant(const Matrix<1, 1> &mat)
{
return mat.flat()[0];
}
template <size_t M>
Matrix<M, M> adjoint(const Matrix<M, M> &mat)
{
Matrix<M, M> result{};
for (int i = 0; i < M; i++)
{
for (int j = 0; j < M; j++)
{
auto sign = (i + j) % 2 == 0 ? 1 : -1;
result(j, i) = sign * (determinant(cofactor(mat, i, j)));
}
}
return result;
}
template <>
Matrix<1, 1> adjoint(const Matrix<1, 1, nullptr> &mat)
{
return {1};
}这样做有蛮大的缺陷。通过代数余子式求行列式的过程是模板递归的，对维数N的矩阵调用一次行列式，将会实例化N...1维的所有矩阵模板，这根本无法接受。因此我不使用这种方法实现行列式计算。
另一个方法是对矩阵做LU分解： A=PLU 。由于LU是三角矩阵，行列式等于对角线乘积 det(A)=det(P)det(diag(U_{11}...U_{nn})) ，得到LU矩阵就求出了行列式。另一个好处在于LU分解求方阵逆差不多是计算量最小方法。实现PLU能一次性求出行列式与逆。
Pivot LU
LU分解，前后迭代（forward/backward substitute）这些常规操作我在专栏稀疏求解里面已经说过了，这里不赘述。直接给出公式：
\beta_{ij}=a_{ij}-\sum_{k=1}^{i-1}\alpha_{ik}\beta_{kj}\tag{1} \\ \alpha_{ij}=\frac{1}{\beta_{jj}}(\alpha_{ij}-\sum_{k=1}^{j-1}\alpha_{ik}\beta_{kj})\tag{2} \\ 这个是朴素LU，我们需要选主元保证稳定性。选主元就是交换行的过程，交换单位变换矩阵存在列数组中：
template <size_t N>
Matrix<N, N> ludcmp(Matrix<N, N> A, std::array<int, N> &indx, bool &even)
{
even = true;
for (int i = 0; i < N; i++)
{
indx = i;
}
for (int k = 0; k < N - 1; k++)
{
auto valmax = fabs(A(k, k));
auto ip = k;
for (int row = k + 1; row < N; row++)
{
double tmp = fabs(A(row, k));
if (valmax < tmp)
{
valmax = tmp;
ip = row;
}
}
if (valmax < gl_rep_eps)
{
return {};
}
if (ip != k)
{
for (int col = k; col < N; col++)
{
std::swap(A(ip, col), A(k, col));
}
std::swap(indx[ip], indx[k]);
for (int col = 0; col < k; col++)
{
std::swap(A(k, col), A(ip, col));
}
even = !even;
}
for (int row = k + 1; row < N; row++)
{
double weight = A(row, k) / A(k, k);
A(row, k) = weight;
for (int col = k + 1; col < N; col++)
{
A(row, col) -= weight * A(k, col);
}
}
}
return A;
}求行列式直接取U的对角线乘积，然后根据选主元结果交换次数乘正负号。
求逆就是拿LU去解按列解单位阵，不过顺序按P交换：
template <size_t N>
void ludbksb(const Matrix<N, N> &A, const std::array<int, N> &indx, double *b)
{
std::array<double, N> y{};
y[0] = b[indx[0]];
for (int row = 1; row < N; row++)
{
double sum = 0.0;
for (int col = 0; col < row; col++)
{
sum += A(row, col) * y[col];
}
y[row] = b[indx[row]] - sum;
}
int n = N;
b[n - 1] = y[n - 1] / A(n - 1, n - 1);
for (int row = n - 2; row >= 0; row--)
{
auto id = row + 1;
double sum = 0.0;
for (int col = id; col < n; col++)
{
sum += A(row, col) * b[col];
}
b[row] = (y[row] - sum) / A(row, row);
}
}这里LU的接口参考LAPACK设计。这种实现应该是非并行情况下最优的。
benchmark测一下500x500稠密矩阵逆（这里引入了Eigen作为参考）：
2022-11-13T14:23:39+08:00
Running C:\Users\Desktop\matrix\build\Release\matrix_bench.exe
Run on (12 X 1617.98 MHz CPU s)
CPU Caches:
L1 Data 32 KiB (x6)
L1 Instruction 32 KiB (x6)
L2 Unified 256 KiB (x6)
L3 Unified 12288 KiB (x1)
------------------------------------------------------------
Benchmark Time CPU Iterations
------------------------------------------------------------
BM_MatrixInv 2030 ns 2002 ns 320000
BM_MatrixEigenInv 2846 ns 2846 ns 280000稠密矩阵逆稍微比Eigen快一些。
经过这样一番简单操作以后，我们得到了一个能像Fortran/Matlab一样加减乘、逐元、切片的矩阵类。并且这个矩阵类性能相对还不错（稠密乘法约Eigen一半，逆速度与LAPACK相当）
这个矩阵类大约500行，都放在头文件中：
之后我会基于这个做一些李群代数的代码工作，都放在这个库里。 |
|
发表于 2022-12-1 19:06:04