Eigen 高维矩阵运算

本文最后更新于：2024年5月7日下午

Eigen 官方代码仅支持二维矩阵，但其他贡献值提供了高维矩阵处理类 Tensor。

编译速度

准备用 Eigen 做矩阵运算的同学注意这个库在 Release O1 及以上的优化标准下可能需要编译十几分钟！
如果使用禁用高级优化 (\Od) 选项，代码运行速度比 O2 慢2倍多
因此可以在开发时 \Od，在发布时 O2
在 VS 中调整的方法为：属性 -> C/C++ -> 优化 -> 优化下拉选项

Tensor 类

Matrix 和 Array 表示二维矩阵，对于任意维度的矩阵可以使用 Tensor 类（当前最高支持 250 维）
注意：这部分代码是用户提供的，没有获得 Eigen 官方支持，不在官方文档支持的代码包里
官方文档（注明了 unsupported）：https://eigen.tuxfamily.org/dox/unsupported/eigen_tensors.html#title15
仓库链接：https://gitlab.com/libeigen/eigen

引用方法

Tensor 在 eigen\unsupported\Eigen\CXX11 文件夹下，使用时引用：
1
#include <unsupported/Eigen/CXX11/Tensor>
之后可以使用 Tensor 类的相关部分代码。

创建 Tensor 对象

Tensor 也有静态、动态之分，用法和 Matrix、 Array 不同

动态、静态对象

动态 Tensor 语法：

1 2	`Tensor<data_type, rank>(size0, size1, ...) Tensor<data_type, rank>(size_array)`

静态 Tensor 语法：

1	`TensorFixedSize<data_type, Sizes<size0, size1, ...>>`

示例代码：

int a = 2;
Eigen::Tensor<int, 3> t(a, a, 3);
t.setRandom();

cout << t << endl << endl;

Eigen::TensorFixedSize<float, Sizes<2, 2, 2>> f;
f.setConstant(3);
cout << f << endl << endl;

-->
  -84756595  -556960047  -680927012
 2146412165 -1365162692 -1063466897

  -98229169  -283100286 -1721105372
   52398681 -1246925249 -1728135974

3 3
3 3

3 3
3 3

也可以创建 string 类型的 Tensor 对象：

// Create a tensor of strings of rank 2 with sizes 2, 2.
Tensor<string, 2> t_2d({ 2, 2 });
t_2d(0, 0) = "hello";
t_2d(0, 1) = "eigen";
t_2d(1, 0) = "tensor";
t_2d(1, 1) = "world!";
cout << t_2d << endl << endl;

-->
hello  eigen
tensor world!

TensorMap

TensorMap 可以从已经分配内存的数据生成 Tensor 对象

语法：

1	`TensorMap<Tensor<data_type, rank>>(data, size0, size1, ...)`

示例：

int storage[36];
TensorMap<Tensor<int, 4>> t_4d(storage, 2, 3, 2, 3);

-->
1219955604 1258332037  304939008
 343989211  347455365  347499254

1219958005 1219958049 1219955516
         2          0          0

 348393272          0 1219965800
         2          0          0

     32758      22260   67094713
     32761      32761      32761

     32758      32758      32758
     32758          0          0

     32761          0      32758
     32761          0          0

初始化

Tensor 初始化和 Martix 差不多

方法	语法	示例
指定常数初始化	`setConstant(const Scalar& val)`	`a.setConstant(12.3f);` `a.setConstant("yolo");`
0 值初始化	`setZero()`	`a.setZero();`
赋值初始化	`setValues({..initializer_list})`	`Eigen::Tensor<float, 2> a(2, 3); a.setValues({{0.0f, 1.0f, 2.0f}, {3.0f, 4.0f, 5.0f}});`
随机初始化	`setRandom()`	`a.setRandom();`

索引数据

单个数据

语法：
1
<data_type> tensor(index0, index1...)

示例：

int a = 2;
Eigen::Tensor<int, 3> t(a, a, 3);
t.setRandom();

cout << t(1, 2, 0) << endl << endl;

-->
-283100286

切片

当需要引入成块数据时， Tensor 没有 Matrix 类有那么方便的 block 函数，但是支持切片操作
切片需要设置 offset 和 extents 两个变量，offset 表示块左上角坐标，extents 表示块维度

语法：

1	`slice(const StartIndices& offsets, const Sizes& extents)`

示例：

int a = 2;
Eigen::Tensor<int, 3> t(a, a, 3);
t.setConstant(2);
cout << t << endl << endl;

const Eigen::Tensor<int, 1>::Dimensions d(2);
Eigen::array<Eigen::Index, 3> offsets = { 0, 0, 0};
Eigen::array<Eigen::Index, 3> extents = { 2, 2, 2};

t.slice(offsets, extents).setZero();

cout << t;

-->
2 2 2
2 2 2

2 2 2
2 2 2

0 0 2
0 0 2

0 0 2
0 0 2

可以看到，原始数据以 $[0,0,0]$ 起始的 $2 \times2\times2$ 的区域内都被切片设置成了0，说明切片起了作用，而且切片的数据是引用。

切片也可以相互赋值：

Eigen::array<Eigen::Index, 3> offsets1 = {0, 0, 0};
Eigen::array<Eigen::Index, 3> offsets2 = { 1, 1, 0 };
Eigen::array<Eigen::Index, 3> extents = {200, 200, 1};

float_tensor.slice(offsets1, extents) = float_tensor.slice(offsets2, extents);

另外一种特殊的切片叫做 chip

语法：

1	`slice(const StartIndices& offsets, const Sizes& extents)`

示例：

Eigen::Tensor<int, 2> a(4, 3);
a.setValues({{0, 100, 200}, {300, 400, 500},
             {600, 700, 800}, {900, 1000, 1100}});
Eigen::array<Eigen::Index, 2> offsets = {1, 0};
Eigen::array<Eigen::Index, 2> extents = {2, 2};
Eigen::Tensor<int, 2> slice = a.slice(offsets, extents);
cout << "a" << endl << a << endl;

-->
a
   0   100   200
 300   400   500
 600   700   800
 900  1000  1100
cout << "slice" << endl << slice << endl;

-->
slice
 300   400
 600   700

引用

如果只需要从表达式的值中访问几个元素，那么可以通过使用 TensorRef 避免在完整张量中具体化该值。
TensorRef 是任何特征操作的小包装类。它为()操作符提供重载，允许您访问表达式中的各个值。TensorRef 很方便，因为 Operation 本身不提供访问单个元素的方法。

只有在需要表达式值的子集时才使用 TensorRef。TensorRef 只计算您访问的值。但是请注意，如果你要访问所有的值，Tensor 计算将会更快一些。

// Create a TensorRef for the expression.  The expression is not
// evaluated yet.
TensorRef<Tensor<float, 3> > ref = ((t1 + t2) * 0.2f).exp();

// Use "ref" to access individual elements.  The expression is evaluated
// on the fly.
float at_0 = ref(0, 0, 0);
cout << ref(0, 1, 0);

数组表示

Tensor 可以将数据部分以数组指针形式传出来，在数组中修改数据

Eigen::Tensor<float, 2> a(3, 4);
float* a_data = a.data();
a_data[0] = 123.45f;
cout << "a(0, 0): " << a(0, 0);
=> a(0, 0): 123.45

表达式计算规则

Tensor 的表达式可以仅以操作符而不是计算的结果的形式存在
计算表达式最常用的方法是将其赋给张量，而使用 auto 则可以保留运算过程而不计算结果。

auto 保留计算

在下面的示例中，auto 声明使中间值为 Operations，而不是 Tensors，并且不会导致计算表达式。对张量结果的赋值将导致对所有操作的计算。

auto t3 = t1 + t2;             // t3 is an Operation.
auto t4 = t3 * 0.2f;           // t4 is an Operation.
auto t5 = t4.exp();            // t5 is an Operation.
Tensor<float, 3> result = t5;  // The operations are evaluated.

如果知道 Operation 值的等级和大小，可以将 Operation 赋值为 TensorFixedSize 而不是张量，这样计算效率能高一点。
1
2
// We know that the result is a 4x4x2 tensor! TensorFixedSize<float, Sizes<4, 4, 2>> result = t5;

eval 强制计算

在计算大型复合表达式时，有时需要告诉 Eigen 表达式树中的一个中间值值得提前计算。这是通过插入对表达式 Operation 的 eval()方法的调用来完成的。

// The previous example could have been written:
Tensor<float, 3> result = ((t1 + t2) * 0.2f).exp();

// If you want to compute (t1 + t2) once ahead of time you can write:
Tensor<float, 3> result = ((t1 + t2).eval() * 0.2f).exp();

可以避免重复计算：

Broadcasting ()表达式导致对 X.max ()表达式进行多次计算:
1
2
3
Tensor<...> X ...; Tensor<...> Y = ((X - X.maximum(depth_dim).reshape(dims2d).broadcast(bcast)) * beta).exp();
在 maximum() 和 reshape() 调用之间插入 eval () 调用可以保证maximum() 只计算一次，从而大大加快执行速度。

控制计算设备

张量库提供了诸如收缩和卷积等各种运算的几种实现。这些实现针对不同的环境进行了优化: CPU 上的单线程，CPU 上的多线程，或者使用 Cuda 的 GPU。

可以在指定设备上运行计算功能：

1 2	`Eigen::Tensor<float, 2> c(30, 40); c.device(...) = a + b;`

多线程计算

使用线程池进行计算

// Create the Eigen ThreadPool
Eigen::ThreadPool pool(8 /* number of threads in pool */)

// Create the Eigen ThreadPoolDevice.
Eigen::ThreadPoolDevice my_device(&pool, 4 /* number of threads to use */);

// Now just use the device when evaluating expressions.
Eigen::Tensor<float, 2> c(30, 50);
c.device(my_device) = a.contract(b, dot_product_dims);

属性获取

示例数据：

1 2	`Eigen::Tensor<int, 3> t(3, 3, 3); t.setConstant(2);`

属性	语法	示例
所有维度尺寸	`.dimensions()`	`t.dimensions() --> [3, 3, 3]`
维度数	`.NumDimensions`	`t.NumDimensions --> 3`
指定维度指定	`.dimension(Index n)`	`t.dimension(1) --> 3`
数据数量	`.size()`	`t.size() --> 27`

正经的 Tensor 对象是可以获取上述属性的，但是 Operation 就不一定了
比较好的办法是用 TensorRef 指向Tensor 对象，以在没有计算时获取其属性。

常用操作

矩阵运算

操作	语法	示例
生成和当前矩阵一样大的常数矩阵	`constant(const Scalar& val)`	`a.constant(2.0f);`
生成和当前矩阵一样大的随机数矩阵	`random()`	`a.random();`
逐元素加、减、乘、除、负号	`+, -, *, /, -`	`a + b; -a`
逐元素开方	`sqrt()`	`a.sqrt()`
逐元素开方取倒数	`rsqrt()`	`a.rsqrt()`
逐元素平方	`square()`	`a.square()`
逐元素取倒数	`inverse()`	`a.inverse()`
逐元素自然指数	`exp()`	`a.exp()`
逐元素自然对数	`log()`	`a.log()`
逐元素幂运算	`pow(p)`	`a.pow(2)`
两个矩阵中逐元素取最大值	`cwiseMax(q)`	`a.cwiseMax(q)`
两个矩阵中逐元素取最小值	`cwiseMin(q)`	`a.cwiseMin(q)`
根据真假选择矩阵	`select(const ThenDerived& thenTensor, const ElseDerived& elseTensor)`	`f.select(then, else);`
所有元素求和	`sum()`	`a.sum()`
按指定维度求和	`sum(const Dimensions& new_dims)`	`a.sum(Eigen::array<int, 2>({0, 1}))`
所有元素求均值	`mean()`	`a.mean()`
按指定维度求均值	`mean(const Dimensions& new_dims)`	`a.mean(Eigen::array<int, 2>({0, 1}))`
所有元素最大值	`maximum()`	`a.maximum()`
按指定维度求最大值	`maximum(const Dimensions& new_dims)`	`a.maximum(Eigen::array<int, 2>({0, 1}))`
所有元素最小值	`minimum()`	`a.minimum()`
按指定维度求最小值	`minimum(const Dimensions& new_dims)`	`a.minimum(Eigen::array<int, 2>({0, 1}))`
迹	`trace()`	`a.trace()`
指定维度迹	`trace(const Dimensions& new_dims)`	`a.trace(Eigen::array<int, 2>({0, 1}))`
指定维度积分图	`cumsum(const Index& axis)`	`a.cumsum(1)`
指定维度乘积图	`cumprod(const Index& axis)`	`a.cumprod(1)`
矩阵反向	`reverse()`	`a.reverse()`
矩阵复制	`broadcast(const Broadcast& broadcast)`	`a.broadcast(bcast);`
pad 0	`pad()`	`a.pad(paddings)`
转换数据类型	`cast<type>()`	`a.cast<float>()`
按维度取 argmax	`argmax(int dim)`	`a.argmax(2)`
按维度取 argmin	`argmin(int dim)`	`a.argmin(2)`

constant

可以相当加、减、乘、除常数操作

Eigen::Tensor<float, 2> a(2, 3);
a.setConstant(1.0f);
Eigen::Tensor<float, 2> b = a + a.constant(2.0f);
Eigen::Tensor<float, 2> c = b * b.constant(0.2f);
cout << "a" << endl << a << endl << endl;
cout << "b" << endl << b << endl << endl;
cout << "c" << endl << c << endl << endl;

-->
a
1 1 1
1 1 1

b
3 3 3
3 3 3

c
0.6 0.6 0.6
0.6 0.6 0.6

random

Eigen::Tensor<float, 2> a(2, 3);
a.setConstant(1.0f);
Eigen::Tensor<float, 2> b = a + a.random();
cout << "a" << endl << a << endl << endl;
cout << "b" << endl << b << endl << endl;
=>
a
1 1 1
1 1 1

b
1.68038   1.5662  1.82329
0.788766  1.59688 0.395103

sum

指定维度求和，两种写法

Eigen::Tensor<float, 3> t(3, 3, 3);
t.setConstant(2);

const Eigen::Tensor<int, 2>::Dimensions d(0, 1); // 将矩阵前两维坍缩成和
cout << t.sum(d) << endl;

cout << t.sum(Eigen::array<int, 2>({ 0, 1 }));

-->
18 18 18
18 18 18

cumsum

列方向求和

// Create a tensor of 2 dimensions
Eigen::Tensor<int, 2> a(2, 3);
a.setValues({{1, 2, 3}, {4, 5, 6}});
// Scan it along the second dimension (1) using summation
Eigen::Tensor<int, 2> b = a.cumsum(1);
// The result is a tensor with the same size as the input
cout << "a" << endl << a << endl << endl;
cout << "b" << endl << b << endl << endl;
=>
a
1 2 3
4 5 6

b
1  3  6
4  9 15

可以原地分别求各个维度的和，最终得到积分图

// Create a tensor of 2 dimensions
Eigen::Tensor<int, 2> a(2, 3);
a.setValues({ {1, 2, 3}, {4, 5, 6} });
// Scan it along the second dimension (1) using summation
cout << "a" << endl << a << endl << endl;
a = a.cumsum(1);
// The result is a tensor with the same size as the input
cout << "a" << endl << a << endl << endl;
a = a.cumsum(0);
// The result is a tensor with the same size as the input
cout << "a" << endl << a << endl << endl;

-->
a
1 2 3
4 5 6

a
 1  3  6
 4  9 15

a
 1  3  6
 5 12 21

broadcast

Eigen::Tensor<int, 2> a(2, 3);
a.setValues({{0, 100, 200}, {300, 400, 500}});
Eigen::array<int, 2> bcast({3, 2});
Eigen::Tensor<int, 2> b = a.broadcast(bcast);
cout << "a" << endl << a << endl << "b" << endl << b << endl;
=>
a
   0   100   200
 300   400   500
b
   0   100   200    0   100   200
 300   400   500  300   400   500
   0   100   200    0   100   200
 300   400   500  300   400   500
   0   100   200    0   100   200
 300   400   500  300   400   500

pad

Eigen::Tensor<int, 2> a(2, 3);
a.setValues({{0, 100, 200}, {300, 400, 500}});
Eigen::array<pair<int, int>, 2> paddings;
paddings[0] = make_pair(0, 1); // 左右
paddings[1] = make_pair(2, 3); // 上下
Eigen::Tensor<int, 2> b = a.pad(paddings);
cout << "a" << endl << a << endl << "b" << endl << b << endl;
=>
a
   0   100   200
 300   400   500
b
   0     0     0    0
   0     0     0    0
   0   100   200    0
 300   400   500    0
   0     0     0    0
   0     0     0    0
   0     0     0    0

argmax

argmax 和 argmin 是当前（2023.1）的 Tensor 隐藏功能

为啥说他隐藏因为文档上没写
参考：https://stackoverflow.com/questions/62280277/argmax-method-in-c-eigen-library
Tensor 和 Matrix 、Vector 等正统模板一样拥有 arg 功能，可以对某一维度数据求取最值下标

#include <Eigen/Dense>
#include <unsupported/Eigen/CXX11/Tensor>
#include <iostream>

#define STR_(x)  #x
#define STR(x)  STR_(x)
#define PRINT(x)  std::cout << STR(x) << ":\n" << (x) << std::endl

int main()
{
    using namespace Eigen;
    using T = int;
    using S = Sizes<2, 3>;

    S const sizes{};

    T constexpr data[S::total_size]{
        8, 4,
        1, 6,
        9, 2,
    };

    Map<MatrixX<T> const> const matrix(data, sizes[0], sizes[1]);
    PRINT(matrix);

    RowVector2<Index> argmax{};
    matrix.maxCoeff(&argmax.x(), &argmax.y());
    PRINT(argmax);

    VectorX<Index> argmax0{matrix.cols()};
    for (Index col = 0; col < matrix.cols(); ++col)
        matrix.col(col).maxCoeff(&argmax0[col]);
    PRINT(argmax0);

    VectorX<Index> argmax1{matrix.rows()};
    for (Index row = 0; row < matrix.rows(); ++row)
        matrix.row(row).maxCoeff(&argmax1[row]);
    PRINT(argmax1);

    TensorMap<Tensor<T const, S::count>> const tensor(data, sizes);
    PRINT(tensor);
    PRINT(tensor.argmax());
    PRINT(tensor.argmax(0));
    PRINT(tensor.argmax(1));

    // Note that tensor.argmax() is the index for a 1D view of the data:
    Index const matrix_index = sizes.IndexOfColMajor(std::array{argmax.x(), argmax.y()});
    Index const tensor_index = Tensor<Index, 0>{tensor.argmax()}();
    PRINT(matrix_index == tensor_index);
}

->
matrix:
8 1 9
4 6 2
argmax:
0 2
argmax0:
0
1
0
argmax1:
2
1
tensor:
8 1 9
4 6 2
tensor.argmax():
4
tensor.argmax(0):
0
1
0
tensor.argmax(1):
2
1
matrix_index == tensor_index:
1

经过 arg 操作的返回结果为 Index 数据，可以这样拿到：

Tensor<float, 3> a(100, 100, 8);
a.setRandom();
Tensor<Index, 2> b;
b = a.argmax(2);

Reduction 操作

sum 或 maximum 和 minimum 函数的结果可以 cout 出来，但事实上他们并不是真正的 double 等普通类型
他们的输出事实上是某种 Tensor 的操作，此时虽然可以看见但是无法使用该数据
如果我们需要使用此类数据，需要赋值给一个 Tensor
参考：https://stackoverflow.com/questions/45940313/eigentensor-double-contraction-to-scalar-value
在 Eigen 文档中也有类似示例使用方法：

// Create a tensor of 2 dimensions
Eigen::Tensor<int, 2> a(2, 3);
a.setValues({{1, 2, 3}, {6, 5, 4}});
// Reduce it along the second dimension (1)...
Eigen::array<int, 1> dims({1 /* dimension to reduce */});
// ...using the "maximum" operator.
// The result is a tensor with one dimension.  The size of
// that dimension is the same as the first (non-reduced) dimension of a.
Eigen::Tensor<int, 1> b = a.maximum(dims);
cout << "a" << endl << a << endl << endl;
cout << "b" << endl << b << endl << endl;
=>
a
1 2 3
6 5 4

b
3
6

对于非矩阵的数值结果，也需要一个 Tensor 对象承接，可以使用空尺寸的 Eigen::TensorFixedSize<double, Eigen::Sizes<>>：

1 2	`Eigen::TensorFixedSize<double, Eigen::Sizes<>> b = a.maximum(); double c = b(0);`

逻辑、比较运算

操作	语法	示例
逐元素与 (bool 型 Tensor 对象)	`&&`	`a && b`
逐元素或 (bool 型 Tensor 对象)	`
逐元素大于	`>`	`a > b`
逐元素不小于	`>=`	`a >= b`
逐元素小于	`<`	`a < b`
逐元素不大于	`<=`	`a <= b`
逐元素等于	`==`	`a == b`
逐元素不等于	`!=`	`a != b`
所有元素为 True	`all()`	`a.all()`
指定维度所有元素为 True	`all(const Dimensions& new_dims)`	`a.all(Eigen::array<int, 2>({0, 1}))`
存在元素为 True	`any()`	`a.any()`
指定维度存在元素为 True	`any(const Dimensions& new_dims)`	`a.any(Eigen::array<int, 2>({0, 1}))`

其他操作

卷积

// Compute convolution along the second and third dimension.
Tensor<float, 4, DataLayout> input(3, 3, 7, 11);
Tensor<float, 2, DataLayout> kernel(2, 2);
Tensor<float, 4, DataLayout> output(3, 2, 6, 11);
input.setRandom();
kernel.setRandom();

Eigen::array<ptrdiff_t, 2> dims({1, 2});  // Specify second and third dimension for convolution.
output = input.convolve(kernel, dims);

for (int i = 0; i < 3; ++i) {
  for (int j = 0; j < 2; ++j) {
    for (int k = 0; k < 6; ++k) {
      for (int l = 0; l < 11; ++l) {
        const float result = output(i,j,k,l);
        const float expected = input(i,j+0,k+0,l) * kernel(0,0) +
                               input(i,j+1,k+0,l) * kernel(1,0) +
                               input(i,j+0,k+1,l) * kernel(0,1) +
                               input(i,j+1,k+1,l) * kernel(1,1);
        VERIFY_IS_APPROX(result, expected);
      }
    }
  }
}

reshape

语法：
1
reshape(const Dimensions& new_dims)

示例：

// Increase the rank of the input tensor by introducing a new dimension
// of size 1.
Tensor<float, 2> input(7, 11);
array<int, 3> three_dims{{7, 11, 1}};
Tensor<float, 3> result = input.reshape(three_dims);

// Decrease the rank of the input tensor by merging 2 dimensions;
array<int, 1> one_dim{{7 * 11}};
Tensor<float, 1> result = input.reshape(one_dim);

shuffle

我运行时报错，怀疑文档和代码不匹配

语法：
1
shuffle(const Shuffle& shuffle)

示例：

// Shuffle all dimensions to the left by 1.
Tensor<float, 3> input(20, 30, 50);
// ... set some values in input.
Tensor<float, 3> output = input.shuffle({1, 2, 0})

eigen_assert(output.dimension(0) == 30);
eigen_assert(output.dimension(1) == 50);
eigen_assert(output.dimension(2) == 20);

stride

语法：
1
stride(const Strides& strides)

示例：

Eigen::Tensor<int, 2> a(4, 3);
a.setValues({{0, 100, 200}, {300, 400, 500}, {600, 700, 800}, {900, 1000, 1100}});
Eigen::array<Eigen::DenseIndex, 2> strides({3, 2});
Eigen::Tensor<int, 2> b = a.stride(strides);
cout << "b" << endl << b << endl;
=>
b
   0   200
 900  1100

printing

Eigen::Tensor<float, 3> tensor3d = {4, 3, 2};
tensor3d.setValues( {{{1, 2}, {3, 4}, {5, 6}}, {{7, 8}, {9, 10}, {11, 12}}, {{13, 14}, {15, 16}, {17, 18}}, {{19, 20}, {21, 22}, {23, 24}}} );
std::cout << tensor3d.format(Eigen::TensorIOFormat::Plain()) << std::endl;
==>
    1  2 
    3  4 
    5  6 

    7  8 
    9 10
    11 12

    13 14
    15 16
    17 18

    19 20
    21 22
    23 24