使用Boost有理数的Eigen矩阵仅4维时编译失败问题咨询

问题概述 #

用Boost多精度库的cpp_rational（分子分母为大整数的有理数类型）配合Eigen计算小矩阵的精确行列式时，2、3、5维矩阵均可正常运行，但4维矩阵会触发编译错误。临时改用partialPivLu().determinant()能解决问题，但3维场景的基准测试显示性能会下降4倍，仅适合作为临时 workaround。

最小复现代码 #

#include <boost/multiprecision/cpp_int.hpp>
#include <Eigen/Dense>
#include <iostream>

using Rational = boost::multiprecision::cpp_rational;

template <size_t N>
Rational rational_test() {
    Eigen::Matrix<Rational, N, N> A;
    for (size_t i = 0; i < N; ++i)
        for (size_t j = 0; j < N; ++j)
            A(i, j) = Rational(1, 1 + i + j);
    return A.determinant();
}

int main() {
    std::cout << rational_test<2>() << "\n"
              << rational_test<3>() << "\n"
              // << rational_test<4>() << "\n"  // 取消注释触发编译错误
              << rational_test<5>() << "\n";
    return 0;
}

核心编译错误信息 #

/usr/include/eigen3/Eigen/src/Core/GenericPacketMath.h:956:1: note:   template 
argument deduction/substitution failed:
/usr/include/eigen3/Eigen/src/LU/Determinant.h:96:56: note:   deduced conflicti
g types for parameter ‘const Packet’ (‘boost::multiprecision::detail::expressi
on<boost::multiprecision::detail::negate, boost::multiprecision::number<boost::
multiprecision::backends::rational_adaptor<boost::multiprecision::backends::cpp
_int_backend<> > >, void, void, void>’ and ‘Eigen::internal::determinant_impl<E
igen::Matrix<boost::multiprecision::number<boost::multiprecision::backends::rat
ional_adaptor<boost::multiprecision::backends::cpp_int_backend<> > >, 4, 4, 0, 
4, 4>, 4>::Scalar’ {aka ‘boost::multiprecision::number<boost::multiprecision::b
ackends::rational_adaptor<boost::multiprecision::backends::cpp_int_backend<> > 
>’})
   96 |     return internal::pmadd(m(i0,2), d0, internal::pmadd(-m(i1,2), d1, m
(i2,2)*d2));
      |                                         ~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~
~~~~~~~~~~

原因分析 #

Eigen对4阶矩阵的行列式计算实现依赖internal::pmadd（打包乘加）函数，而Boost cpp_rational的表达式模板类型（比如-m(i1,2)生成的negate表达式类型）与实际的Scalar类型不匹配，导致模板参数推导失败。其他维度（2、3、5）的行列式计算逻辑未使用该打包运算实现，因此可正常编译。

解决方案 #

方案1：禁用Eigen向量化优化（推荐，性能损失小） #

在包含Eigen头文件前定义宏，强制Eigen使用通用行列式计算实现，避免模板推导冲突：

#define EIGEN_DONT_VECTORIZE
#include <Eigen/Dense>
#include <boost/multiprecision/cpp_int.hpp>
#include <iostream>

// 后续代码保持不变
using Rational = boost::multiprecision::cpp_rational;

template <size_t N>
Rational rational_test() {
    Eigen::Matrix<Rational, N, N> A;
    for (size_t i = 0; i < N; ++i)
        for (size_t j = 0; j < N; ++j)
            A(i, j) = Rational(1, 1 + i + j);
    return A.determinant();
}

int main() {
    std::cout << rational_test<2>() << "\n"
              << rational_test<3>() << "\n"
              << rational_test<4>() << "\n"
              << rational_test<5>() << "\n";
    return 0;
}

方案2：补充Eigen自定义标量适配 #

特化Eigen的internal::is_convertible模板，让Eigen识别Boost多精度的表达式类型可转换为cpp_rational：

#include <boost/multiprecision/cpp_int.hpp>
#include <Eigen/Dense>
#include <type_traits>

using Rational = boost::multiprecision::cpp_rational;

// 特化模板，解决表达式类型与Scalar的转换识别问题
namespace Eigen {
namespace internal {
template <typename Expr>
struct is_convertible<boost::multiprecision::detail::expression<Expr, void, void, void>, Rational>
    : std::true_type {};
}
}

// 后续代码保持不变
template <size_t N>
Rational rational_test() {
    Eigen::Matrix<Rational, N, N> A;
    for (size_t i = 0; i < N; ++i)
        for (size_t j = 0; j < N; ++j)
            A(i, j) = Rational(1, 1 + i + j);
    return A.determinant();
}

int main() {
    std::cout << rational_test<2>() << "\n"
              << rational_test<3>() << "\n"
              << rational_test<4>() << "\n"
              << rational_test<5>() << "\n";
    return 0;
}

该方案无需禁用所有向量化优化，性能更接近原生实现。

方案3：使用LU分解（临时 workaround） #

改用A.partialPivLu().determinant()绕过编译错误，但会带来明显性能损耗，仅适合对性能要求不高的场景。

内容的提问来源于stack exchange，提问作者Daniel Shapero