我要提问
ARTICLE DETAIL

资讯详情

前沿编程新知与开发实战干货的深度解读。

Eigen矩阵逐元素操作与按规则生成:unaryExpr和NullaryExpr实战解析

Eigen矩阵逐元素操作与按规则生成:unaryExpr和NullaryExpr实战解析 写过一段时间的 Eigen应该能感觉到一个“手感差异”加减乘除、转置、求逆、特征分解这些常规操作Eigen 都给你封装好了表达式链随便写性能也很能打。可一旦你想对矩阵的每个元素做点自定义逻辑——比如把角度批量转弧度、做一次分段函数裁剪、或者根据行号列号按公式生成一个初值矩阵——默认的接口就不太顺手了。这篇教程讲的就是专门处理这类场景的两个家伙unaryExpr和NullaryExpr。前者负责“把同一个函数作用到矩阵的每个元素上”后者负责“按任意规则现场生成矩阵”。适合两类人看一是用 Eigen 写数值求解器、做数据预处理的人二是被 Eigen 表达式模板绕晕、想弄明白“为什么它能这么写”的初学者。1. 为什么需要这两个特殊的表达式类1.1 表达式模板机制的一点背景先聊个背景。Eigen 里a b返回的不是一个全新的Matrix对象而是一个“加法表达式”——一个轻量的模板对象里面只保存了对a和b的引用。你真正把它赋值给某个矩阵或者塞进下一个运算表达式时Eigen 才会把求值动作真正执行一遍。这是整个库的性能核心叫“表达式模板”。unaryExpr和NullaryExpr也是这一类机制。unaryExpr构造一个“一元运算表达式”等元素真正被访问时才调用你的函数NullaryExpr构造一个“生成器表达式”它不拥有任何数据每个元素都是现场算出来的。理解这一点很重要因为它决定了两件事第一这些接口不产生额外的中间矩阵内存开销比“临时构造一个矩阵再逐元素改”要低第二它们可以和Matrix、Array、甚至其他表达式直接做运算结果仍然是一条表达式链直到最终赋值才落地。1.2 unaryExpr 和 NullaryExpr 的分工一句话区分unaryExpr是“映射”。给你一个已有矩阵对每个元素应用f(x)得到一个新矩阵。NullaryExpr是“生成”。指定行列数给出一个规则f(i, j)每次都实时算出该位置的元素。它们解决的是同一个问题让“逐元素自定义操作”和“按规则初始化”不依赖手写for循环而是成为一流公民可以参与表达式组合、利用编译期优化。2. unaryExpr对矩阵每个元素施加你的自定义函数2.1 从最简单的 lambda 开始假设你要把矩阵里的每个元素都过一遍 sigmoid。数学上就是s(x) 1 / (1 exp(-x))。如果手写循环你得这样Eigen::MatrixXd m Eigen::MatrixXd::Random(3, 3); Eigen::MatrixXd r(3, 3); for (int i 0; i m.rows(); i) for (int j 0; j m.cols(); j) r(i, j) 1.0 / (1.0 std::exp(-m(i, j)));代码本身不复杂但如果你后续还要对这个结果做转置、乘一个系数、再加另一个矩阵循环就没办法自然地接进表达式链里只能一步一步写临时变量。用unaryExpr就是一行auto sig [](double x) { return 1.0 / (1.0 std::exp(-x)); }; Eigen::MatrixXd r m.unaryExpr(sig);直接传给unaryExpr的 lambda 会被 Eigen 存储为表达式的一部分等到r被赋值时才逐元素调用。编译器一般能把 lambda 内联掉生成的循环和手写版几乎没有差别。2.2 函数指针、函数对象与带状态变换lambda 不是唯一选择。你可以传一个普通函数指针double twice(double x) { return 2.0 * x; } Eigen::MatrixXd r m.unaryExpr(twice);也可以传一个结构体仿函数特别是当变换需要额外参数时。比如你想做一个仿射变换y scale * x offsetstruct AffineTransform { double scale; double offset; double operator()(double x) const { return scale * x offset; } }; Eigen::MatrixXd r m.unaryExpr(AffineTransform{2.0, 1.0});为什么要用仿函数而不是 lambda两种都能写。但仿函数有一个优势它把“状态”和“操作”打包成一个可拷贝的自定义类型Eigen 在多次拷贝表达式时行为更可控。lambda 在 C11 里如果捕获了变量拷贝语义也是明确的但operator()默认是const如果你想让每次调用改变内部状态比如统计调用次数就必须用mutable捕获。而仿函数可以手动把计数器声明为mutable意图更清晰。实际场景里这种带状态的变换不算少做归一化时把均值和方差打包成仿函数传给unaryExpr做非线性回归时把带固定参数的模型表达式传进去代码都比一长串for循环好维护得多。2.3 求值时机、类型推导和一些“手感”细节unaryExpr返回的是一个表达式不是一个已经算好的矩阵。什么时候真正算赋值的时候或者参与其他运算的时候。auto expr m.unaryExpr(sig); // 这里不计算 Eigen::MatrixXd r expr; // 这里才真正遍历执行 sig Eigen::MatrixXd r2 expr m; // 这里也一样逐元素先 sig再相加类型方面有一个容易踩的坑lambda 的参数类型最好和矩阵的标量类型严格一致。如果你的矩阵是floatlambda 却写double func(double)编译报错会发生在 Eigen 模板深处信息很长新手容易懵。稳妥的做法是让 lambda 参数匹配矩阵标量类型Eigen::MatrixXd m ...; // double auto f [](double x) { return x * x; }; auto r m.unaryExpr(f); // OK如果确实不想写死类型可以用泛型 lambdaauto f [](const auto x) { return x * x; };不过矩阵标量类型是double这种简单类型时直接写具体类型可读性更好。注意unaryExpr的仿函数operator()最好声明为const。Eigen 内部有时会在const表达式上调用它如果operator()不是const会报一个很长的编译错误提示“discards qualifiers”。这是我最常碰到的报错后面会有专门一节讲排坑。使用std::function也不是不行但我不推荐。std::function会做类型擦除每次调用都有一层虚函数间接跳转性能比函数指针和 lambda 差不少。更重要的是std::function存在堆分配的可能在 Eigen 这种注重性能的代码里属于“非必要不引入”。3. NullaryExpr不占数据空间按规则实时生成矩阵3.1 从“生成一个矩阵”说起NullaryExpr的官方定义是“一个没有输入参数的表达式”但实际上它有两套形态两参版本NullaryExpr(rows, cols, func)func接受行索引和列索引返回该位置的元素。无参版本NullaryExpr(rows, cols, func)func不接受参数每次访问返回一个值。无参版本能做的事比较有限除了用来生成常量矩阵或者接入带状态的随机数生成器一般用得不多。真正高频使用的是两参版本因为大多数“规则生成”本质上都是“根据坐标算元素”。比如生成一个i j分布的矩阵Eigen::MatrixXd A Eigen::MatrixXd::NullaryExpr(4, 4, [](Eigen::Index i, Eigen::Index j) { return i * 10.0 j; });你得到的就是0 1 2 3 10 11 12 13 20 21 22 23 30 31 32 33NullaryExpr不会先分配一整块内存等你读取A(i,j)的时候才调用 lambda 算出那个位置的值。如果你只是把这个表达式继续往下传比如再加一个矩阵Eigen 会在最终赋值时逐元素计算两个操作数中间不会生成临时矩阵。3.2 三角矩阵、三对角矩阵这种“结构矩阵”怎么造数值计算里经常需要测试矩阵。比如我想生成一个五阶对角占优三对角矩阵主对角元2(i1)两条副对角元1其余位置0。手写三重循环非常啰嗦用NullaryExpr就清晰得多int n 5; auto tridiag [n](Eigen::Index i, Eigen::Index j) - double { if (i j) return 2.0 * (i 1); if (std::abs(i - j) 1) return 1.0; return 0.0; }; Eigen::MatrixXd M Eigen::MatrixXd::NullaryExpr(n, n, tridiag);这个 lambda 里用了分支判断。NullaryExpr对分支没有特殊优化但好处在于你不再需要两三层循环和一堆条件判断来填充矩阵代码结构跟数学定义一一对应后面别人接手时读起来非常轻松。同样地单位矩阵也可以不调用Identity()而用i j的规则生成——虽然没必要但有助于理解机制Eigen::MatrixXd I Eigen::MatrixXd::NullaryExpr(4, 4, [](Eigen::Index i, Eigen::Index j) - double { return i j ? 1.0 : 0.0; });3.3 与 Random、LinSpaced 的区别和组合能力Eigen 自己也有几个“生成类”接口比如MatrixXd::Random(3, 4)、ArrayXd::LinSpaced(100, 0.0, 1.0)。它们和NullaryExpr的区别是Random和LinSpaced是内置好的特化规则效率当然没问题NullaryExpr则是“规则生成器”的通用入口你想生成什么规律都行自由度完全在你的 lambda 里。更关键的是组合能力。NullaryExpr返回的是一个表达式所以可以直接参与算术Eigen::MatrixXd B Eigen::MatrixXd::Random(3, 3); Eigen::MatrixXd C Eigen::MatrixXd::NullaryExpr(3, 3, [](Eigen::Index i, Eigen::Index j) { return i - j; }) 2.0 * B;这条语句里NullaryExpr生成一个“坐标差矩阵”表达式然后和一个随机矩阵的 2 倍做逐元素相加最终结果赋给C。整个过程没有多余的中间矩阵一次遍历完成全部计算。如果某个规则生成的结果你会反复使用那建议先eval()出来存成真正的矩阵。否则每次访问都重算一遍 lambda在迭代算法里会造成大量无用功Eigen::MatrixXd M Eigen::MatrixXd::NullaryExpr(...).eval(); // 落成真矩阵4. 我踩过的坑和实测结论4.1 常见问题速查表unaryExpr和NullaryExpr用起来其实不难难的是编译报错和表达式拷贝时的状态问题。我把实战中遇到的高频问题整理成了一张表问题现象根本原因解决方案编译错误passing ... as this argument discards qualifiers仿函数或 lambda 的operator()不是constEigen 在const表达式上调用它给operator()加const需要改动的成员加mutable编译错误no matching function for call to unaryExprlambda 参数类型和矩阵标量类型不匹配让参数类型严格等于矩阵标量类型或用const autostd::function版本运行明显慢类型擦除、函数指针跳转、堆分配用 lambda、函数指针或仿函数替换表达式被拷贝多次后mutablelambda 的状态“分裂”Eigen 表达式对象本身会被拷贝内部保存的仿函数也被拷贝计数器各算各的尽量用仿函数并把表达式立即eval()不要在auto上长时间持有带状态的表达式用NullaryExpr生成大矩阵后反复访问程序变慢每次访问都重新执行 lambda先用.eval()转成真矩阵再使用第 4 条我吃过亏值得多说两句。Eigen 的表达式对象是可拷贝的auto expr m.unaryExpr(counterFunctor)之后如果把expr传进函数、放进容器甚至只是拷贝了一份内部的仿函数状态就和原来的分家了。最后赋值求值时到底哪一份状态生效你根本想不到。解决思路就两条要么仿函数里不要写可变状态要么从构造到赋值都直来直去别把表达式存下来到处传。4.2 unaryExpr 到底会不会拖慢性能我对这三种写法做过简单实测手写双重循环、unaryExpr传 lambda、Array直接表达式在-O2优化下耗时几乎完全一致。因为 lambda 是内联的Eigen 的逐元素访问在 256 位 SIMD 下也能自动向量化手写循环能走的路unaryExpr一样能走。真正有差距的是std::function版本因为每次元素访问都要经过一次间接调用。虽然现代 CPU 有分支预测虚调用在某些热点循环里可能被优化掉一部分但在我本机测试中普遍比 lambda 版本慢 3 到 5 倍具体看平台。结论很简单能用 lambda 或仿函数就别用std::function。4.3 关于四元数的一个澄清热词榜里有“eigen库四元数”这里插一句。Quaterniond和unaryExpr其实是两个维度的东西unaryExpr是针对矩阵中每一个标量元素做映射不是针对每个“对象元素”。如果你想给一批四元数做归一化、共轭或按某个公式生成四元数矩阵unaryExpr并不能直接优雅地替你做。更合适的做法是把四元数放进std::vector用std::transform写显式循环或者用Eigen::Map把四元数数组按分量映射成一个Matrixdouble, 4, N再把分量操作写清楚。这个和本主题无关但看到关键词就提前说清楚免得大家把两个概念搅在一起。5. 新手补丁在 VS Code 里把 Eigen 环境配好讲完进阶用法最后回归新手问题怎么在 VS Code 里加载 Eigen。Eigen 是 header-only 的库没有.dll、.so、.lib需要链接。你要做的只有一件事让编译器能找到Eigen/Dense这套头文件所在目录。如果用的 CMake 构建最标准的写法是cmake_minimum_required(VERSION 3.16) project(eigen_demo) find_package(Eigen3 REQUIRED) add_executable(demo main.cpp) target_link_libraries(demo Eigen3::Eigen)Eigen3::Eigen这个 target 本质上是头文件搜索路径不涉及链接。在 Ubuntu 系系统里用apt install libeigen3-dev安装后头文件在/usr/include/eigen3CMake 的find_package通常也能直接找到。如果不用 CMake而是直接在 VS Code 里编译单文件可以打开c_cpp_properties.json在includePath里加上 Eigen 的路径{ configurations: [ { name: Linux, includePath: [ ${workspaceFolder}/**, /usr/include/eigen3 ], defines: [], compilerPath: /usr/bin/g, cStandard: c17, cppStandard: c17, intelliSenseMode: linux-gcc-x64 } ], version: 4 }这样 IntelliSense 不报红编译命令里也得带上同样的-I参数。最省事的是写个简单的tasks.json里面用-I/usr/include/eigen3编译或者干脆用 CMake Tools 插件一条CMake: Configure搞定所有路径问题。如果出现fatal error: Eigen/Dense: No such file or directory先确认三点头文件确实在那个目录、includePath或-I参数写的目录正确、目录层级对不对——Eigen 的头文件是Eigen/Dense所以路径必须指到包含Eigen文件夹的上一级。最后说点实在的。我平时写数值代码的体会是unaryExpr和NullaryExpr属于那种“知道的人觉得顺手不知道的人一辈子用不上”的 API。它们不会改变算法复杂度但能让代码离数学定义更近、少写几层循环排查问题的时候一眼能看穿。如果你刚开始接触 Eigen建议先把Matrix和Array的基础操作玩熟等你开始觉得“这一步如果不写for循环就好了”的时候就是这两个表达式类派上用场的时机了。
返回列表