Eigen无矩阵BiCGSTAB报错：THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY

问题背景 #

需要求解规模约100k×100k的稠密浮点矩阵方程Ax=y（已知y求x），计划采用Eigen库的BiCGSTAB迭代算法。因内存限制无法存储完整矩阵A，参考Eigen官方示例实现了模拟单位矩阵的自定义无矩阵包装类。代码在不调用BiCGSTAB时可正常运行，但调用该算法时触发编译错误：

static assertion failed: THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY

错误原因 #

Eigen的迭代求解器（包括BiCGSTAB）要求传入的矩阵类型具备左值特性，你的自定义矩阵类存在两处关键问题：

• 在Eigen::internal::traits<MyMatrixReplacement>中，Flags被设为0，缺少Eigen::LvalueBit标记，导致Eigen无法识别该类为可作为左值的类型。
• 迭代求解器内部逻辑会尝试对矩阵进行左值绑定操作，缺少该标记会触发静态断言阻止编译。

解决办法 #

1. 修正traits中的Flags配置 #

在traits结构体的枚举中添加必要标记位，至少包含LvalueBit和布局标记：

enum { 
  Flags = Eigen::LvalueBit | Eigen::RowMajorBit, // 可根据矩阵布局选择ColMajorBit
  RowsAtCompileTime = Eigen::Dynamic,
  ColsAtCompileTime = Eigen::Dynamic,
  MaxRowsAtCompileTime = Eigen::Dynamic,
  MaxColsAtCompileTime = Eigen::Dynamic,
  IsVectorAtCompileTime = false,
  Alignment = Eigen::Aligned16, // 使用标准对齐而非自定义值1
  outer_stride_at_compile_time = Eigen::Dynamic,
  inner_stride_at_compile_time = Eigen::Dynamic
};

2. 完善自定义矩阵类的接口 #

• 移除不必要的非constdata()方法，或改为返回const空指针：

  const double* data() const { return nullptr; }
  

• 调整innerStride()和outerStride()的返回值，符合Eigen的 stride 规则：

  Index innerStride() const { return Index(1); }
  Index outerStride() const { return Index(rows()); }
  

3. 修正evaluator的Flags #

在evaluator<MyMatrixReplacement>的枚举中添加访问标记：

enum { 
  CoeffReadCost = 1,
  Flags = Eigen::PacketAccessBit | Eigen::LinearAccessBit
};

替代方案建议 #

针对100k×100k规模的矩阵求解，除BiCGSTAB外可考虑以下方案：

• 利用矩阵结构优化：若矩阵具备稀疏、对称正定、带状等特性，选择对应迭代器（如对称正定矩阵用ConjugateGradient，稀疏矩阵用SparseLU）。
• 分块求解：将大矩阵拆分为小分块，分批次计算矩阵-向量乘积，降低内存占用。
• 预处理技术：为BiCGSTAB添加高效预处理器（如IncompleteLUT、DiagonalPreconditioner），加速收敛。
• 第三方库：若Eigen无法满足需求，可尝试PETSc、Trilinos等专门针对大规模线性系统的数值计算库。

修正后的完整测试代码 #

#include <iostream>
#include <eigen3/Eigen/Core>
#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/IterativeLinearSolvers>
#include <eigen3/unsupported/Eigen/IterativeSolvers>

class MyMatrixReplacement;

namespace Eigen {
  namespace internal {
    template<>
    struct traits<MyMatrixReplacement>
    {
      typedef Eigen::Dense StorageKind;
      typedef Eigen::MatrixXpr XprKind;
      typedef double Scalar;
      typedef double RealScalar;
      typedef int StorageIndex;
      typedef Eigen::Index Index;
      enum { 
        Flags = Eigen::LvalueBit | Eigen::RowMajorBit,
        RowsAtCompileTime = Eigen::Dynamic,
        ColsAtCompileTime = Eigen::Dynamic,
        MaxRowsAtCompileTime = Eigen::Dynamic,
        MaxColsAtCompileTime = Eigen::Dynamic,
        IsVectorAtCompileTime = false,
        Alignment = Eigen::Aligned16,
        outer_stride_at_compile_time = Eigen::Dynamic,
        inner_stride_at_compile_time = Eigen::Dynamic
      };
    };
  }
}


class MyMatrixReplacement : public Eigen::MatrixBase<MyMatrixReplacement>
{
public:
  MyMatrixReplacement(){}
  MyMatrixReplacement(int x, int y){}

  typedef typename Eigen::internal::ref_selector<MyMatrixReplacement>::type Nested; 

  Index rows() const {return Index(4);};
  Index cols() const {return Index(4);};
 
  template<typename Rhs>
  Eigen::Product<MyMatrixReplacement,Rhs,Eigen::AliasFreeProduct> operator*(const Eigen::MatrixBase<Rhs>& x) const {
    return Eigen::Product<MyMatrixReplacement,Rhs,Eigen::AliasFreeProduct>(*this, x.derived());
  }

  const double* data() const {
        return nullptr;
      }

  Index innerStride() const {return Index(1);}
  Index outerStride() const {return Index(rows());}

  MyMatrixReplacement &operator=(const MyMatrixReplacement& rhs){return *this;}

  MyMatrixReplacement(const MyMatrixReplacement &rhs){}

  ~MyMatrixReplacement(){ }
};

namespace Eigen {
  namespace internal {
    template<>
    struct evaluator<MyMatrixReplacement>
      : evaluator_base<MyMatrixReplacement>
    {
      typedef MyMatrixReplacement XprType;
      typedef typename XprType::CoeffReturnType CoeffReturnType;
      enum { 
        CoeffReadCost = 1,
        Flags = Eigen::PacketAccessBit | Eigen::LinearAccessBit
      };
      
      evaluator(const MyMatrixReplacement& x)
      {}

      CoeffReturnType coeff(Index row, Index col) const
      {
        double entry = 0;
        if(row == col){
          entry = 1;
        }
        return entry;
      }

    };
  }
}


namespace Eigen {
namespace internal {
 
  template<>
  struct generic_product_impl<MyMatrixReplacement, Eigen::VectorXd>
  : generic_product_impl_base<MyMatrixReplacement,Eigen::VectorXd,generic_product_impl<MyMatrixReplacement,Eigen::VectorXd> >
  {
    typedef typename Product<MyMatrixReplacement,Eigen::VectorXd>::Scalar Scalar;
 
    template<typename Dest>
    static void scaleAndAddTo(Dest& dst, const MyMatrixReplacement& lhs, const Eigen::VectorXd& rhs, const Scalar& alpha)
    {
      assert(alpha==Scalar(1) && "scaling is not implemented");
      EIGEN_ONLY_USED_FOR_DEBUG(alpha);
      dst = rhs;
    }
  };
 
}
}


int main()
{
  Eigen::VectorXd vec(4);
  vec << 1, 2, 4, 8;
  MyMatrixReplacement mat;
  std::cout << mat << std::endl;
  Eigen::VectorXd x(4);
  Eigen::VectorXd vec2(4);
  vec2 = mat*vec;
  std::cout << vec2 << std::endl;
  MyMatrixReplacement mat2;
  mat = mat2;
  
  {
    Eigen::BiCGSTAB<MyMatrixReplacement, Eigen::IdentityPreconditioner> bicg;
    bicg.compute(mat);
    x = bicg.solve(vec);
    std::cout << "BiCGSTAB: #iterations: " << bicg.iterations() << ", estimated error: " << bicg.error() << std::endl;
  }
  
}

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