Eigen自定义标量类型：4x4矩阵求逆可编译，5x5矩阵报错

为啥会这样？ #

Eigen处理矩阵求逆分两种逻辑：

• 4x4及更小矩阵：用硬编码的解析公式直接计算逆，只需要自定义标量支持加减乘除基础运算，所以你之前的代码能正常运行。
• 5x5及更大矩阵：默认走LU分解的数值求逆流程，这个过程会调用abs()（判断主元大小）、</>（选择主元）这些操作，但你的int_mod<3>类型没重载这些函数/运算符，直接导致编译失败。

怎么解决？ #

给int_mod<3>补上缺失的运算符和函数重载即可，分以下几步：

1. 重载比较运算符 #

有限域元素可以直接基于底层存储的整数值实现比较：

template<int MOD>
bool operator<(const int_mod<MOD>& a, const int_mod<MOD>& b) {
    return a.val < b.val;
}

template<int MOD>
bool operator>(const int_mod<MOD>& a, const int_mod<MOD>& b) {
    return a.val > b.val;
}

// 顺带补上==、!=、<=、>=，避免后续其他场景出问题

2. 重载abs()函数 #

有限域元素的绝对值直接返回其非负的模值即可：

template<int MOD>
int_mod<MOD> abs(const int_mod<MOD>& a) {
    return int_mod<MOD>(std::abs(a.val));
}

3. 可选：强制使用解析求逆（不推荐大矩阵） #

如果不想补全重载，也可以强制Eigen调用解析求逆接口，但大矩阵这么做效率极低，仅适合测试场景：

// 强制调用解析求逆方法
Matrix<int_mod<3>, 5, 5> inv_mat = mat.inverseUsingAdjoint();

补全后的完整示例 #

假设你的int_mod<3>已有基础运算重载，补完上述内容后的可运行代码：

#include <Eigen/Dense>
#include <cmath>

template<int MOD>
struct int_mod {
    int val;

    int_mod(int v = 0) : val((v % MOD + MOD) % MOD) {}

    // 基础加减乘除运算
    int_mod operator+(const int_mod& other) const {
        return int_mod(val + other.val);
    }

    int_mod operator-(const int_mod& other) const {
        return int_mod(val - other.val);
    }

    int_mod operator*(const int_mod& other) const {
        return int_mod(val * other.val);
    }

    // 有限域除法（求乘法逆元）
    int_mod operator/(const int_mod& other) const {
        // 用扩展欧几里得算法计算模MOD的逆元
        int a = other.val, m = MOD;
        int m0 = m, y = 0, x = 1;
        if (m == 1) return int_mod(0);
        while (a > 1) {
            int q = a / m;
            int t = m;
            m = a % m;
            a = t;
            t = y;
            y = x - q * y;
            x = t;
        }
        if (x < 0) x += m0;
        return int_mod(val * x);
    }
};

// 补全比较运算符
template<int MOD>
bool operator<(const int_mod<MOD>& a, const int_mod<MOD>& b) {
    return a.val < b.val;
}

template<int MOD>
bool operator>(const int_mod<MOD>& a, const int_mod<MOD>& b) {
    return a.val > b.val;
}

template<int MOD>
bool operator==(const int_mod<MOD>& a, const int_mod<MOD>& b) {
    return a.val == b.val;
}

// 补全abs函数
template<int MOD>
int_mod<MOD> abs(const int_mod<MOD>& a) {
    return int_mod<MOD>(std::abs(a.val));
}

// 给Eigen做类型特化，确保正确适配自定义标量
namespace Eigen {
    template<int MOD>
    struct NumTraits<int_mod<MOD>> : NumTraits<int> {
        typedef int_mod<MOD> Real;
        typedef int_mod<MOD> NonInteger;
        typedef int_mod<MOD> Nested;
        enum {
            IsComplex = 0,
            IsInteger = 0, // 有限域元素支持除法逆元，不属于整数类型
            IsSigned = 1,
            RequireInitialization = 1,
            ReadCost = 1,
            AddCost = 1,
            MulCost = 1
        };
    };
}

int main() {
    Eigen::Matrix<int_mod<3>, 5, 5> mat;
    // 手动赋值一个可逆矩阵，避免随机生成不可逆矩阵
    mat << 1,2,0,0,0,
           0,1,2,0,0,
           0,0,1,2,0,
           0,0,0,1,2,
           2,0,0,0,1;
    Eigen::Matrix<int_mod<3>, 5, 5> inv_mat = mat.inverse();
    return 0;
}

额外注意点 #

• 有限域中矩阵可逆的条件是行列式≠0（模3），测试时尽量手动构造可逆矩阵，别依赖setRandom，大概率会生成不可逆矩阵。
• 补充NumTraits特化是为了让Eigen更好地识别自定义标量，避免其他潜在的适配问题。

内容的提问来源于stack exchange，提问作者Harry Reed