已知矩阵的模板矩阵表达式

Question

这是一道思维练习，不是一个具体的问题，但我很想听听你的意见。假设我有一些使用模板（例如Eigen，ublas等）的矩阵表达式DSL。

现在假设我有一些常量矩阵，例如：

Matrix2 sigma1 = {{0,1}, {1,0}};
Matrix2 sigma2 = {{0,i}, {-i,0}};
... etc

我有一些涉及运行时值的矩阵操作：

a*sigma1 + b*sigma2; // a,b runtime

您有什么想法可以实现常量矩阵，以便编译器能够最优化表达式？特别是如何将(i,j)操作符解析为一个常量？

- Anycorn

---

为了让它成为一个真实的场景，转换矩阵就是这样的，例如，某些给定角度的旋转或某些给定因数的缩放可以被优化掉。 - TC1

“表达式模板”是我想到的。 - dpiskyulev

你可以使用表达式模板，但是需要注意一些（大多数？）编译器在处理浮点数乘以0或1时无法进行优化。这是因为它们仍然应该能够处理代码中的a或b如果是Nan或inf。GCC有一个选项（-ffinite-math-only），可以解决这个问题。（我认为） - sbabbi

@sdabbi ffast-math会丢弃一些IEEE要求。我的目标是看看是否有一种通过模板手段来优化那些乘法的方法。 - Anycorn

@dpiskyulev 是的，我有表达式模板，但问题是如何使用。这有点棘手 :) - Anycorn

Answer 1

据我所了解的问题空间：给定一个领域特定语言，我们想要确定最小的转换，使得某些运算符（例如(i,j)）在操作数据时会查找常量，而不是计算数学公式（例如a*sigma1 + b*sigma2）。

让我们探索一些可能性：

• 直接执行数学运算

  零级实现。如果编译器没有进行明确的优化，那么直接执行指令会发生什么？答案是取决于情况。但是，在大多数处理器上，您的代码执行将落入CPU缓存中，其中您的汇编和分支执行将被优化到最佳状态。这个过程的核心非常酷炫, 但我们假设你想超越这些能力，直接在代码中解决操作的常量。

• 使用编译器编译器来限制和捕获空间

  第一阶段的优化是使用编译器编译器来限制和捕获可能的输入和输出空间。虽然这将有效地解决您的输入范围问题，但除此之外它对性能没有任何帮助。因此，我们必须继续前进。

• 字符串化和宏展开

  第二阶段的优化是直接对值空间进行字符串或宏展开。虽然这充满了角落情况和令人惊讶的实现级坑洼, 但如果需要该操作，您的编译器可以直接执行它。(另请参见：循环展开)

• 手动推导和堆栈边界可满足性闭合形式表达式
  (例如，使用查找表)

  最后，我们的第三阶段优化是直接限制您的空间。这要求您具有具有有界输入和输出空间的明确定义的闭合形式关系才能有效工作。如果无法确定此关系或不包含边界，则您就无法使用，并且应考虑保留当前实现，除非已知存在更好的实现。

这些优化技术中，最适用于线性代数操作的是后两种，考虑到您所描述的问题范围。由于大多数操作（如矩阵平移、旋转和缩放操作）本质上是确定性的，因此您确实可以有效地进行优化和限制空间。

对于更理论性的答案，我建议咨询http://cs.stackexchange.com、http://cstheory.stackexchange.com和http://math.stackexchange.com。两者都有许多帖子专门讨论方程整个类别的可决定性和闭合形式的多项式解的证明。

- MrGomez

Answer 2

好的，这个结构很丑陋，但与我在原始帖子上的评论有关。

使用这个结构，应该可以定义你需要的相关操作，但编写所有适当的特化将是很多工作。你可能也更喜欢交换行/列。

在最后定义矩阵肯定不如在原始帖子中那样优雅，但也许这可以改进，尤其是在C++11中使用'auto'。

//-----------------------------------------------------------------------------
struct Unused {};

struct Imaginary {
    Imaginary() {}
    Imaginary(Unused const& unused) {}
};
struct MinusImaginary {
    MinusImaginary() {}
    MinusImaginary(Unused const& unused) {}
};

//-----------------------------------------------------------------------------
template <int I, int F = 0>
struct Fixed {
    Fixed() {}
    Fixed(Unused const& unused) {}
};

//-----------------------------------------------------------------------------
struct Float
{
    Float(float value) : value_(value) {}
    const float value_;
};

//-----------------------------------------------------------------------------
template <typename COL0, typename COL1>
struct Vector2
{
    typedef COL0 col0_t;
    typedef COL1 col1_t;

    template <typename T0, typename T1>
    Vector2(T0 const& t0, T1 const& t1)
        : col0_(t0)
        , col1_(t1)
    {}

    COL0 col0_;
    COL1 col1_;
};

//-----------------------------------------------------------------------------
template <typename ROW0, typename ROW1>
struct Matrix2
{
    typedef ROW0 row0_t;
    typedef ROW1 row1_t;

    Matrix2()
        : row0_(Unused(), Unused())
        , row1_(Unused(), Unused())
    {
    }
    template <typename M00, typename M01, typename M10, typename M11>
    Matrix2(M00 const& m00, M01 const& m01, M10 const& m10, M11 const& m11)
        : row0_(m00, m01)
        , row1_(m10, m11)
    {
    }

    ROW0 row0_;
    ROW1 row1_;
};

//-----------------------------------------------------------------------------
Matrix2<
    Vector2< Fixed<0>, Fixed<1> >,
    Vector2< Fixed<1>, Fixed<0> > 
> sigma1;

const float f = 0.1f;

//-----------------------------------------------------------------------------
Matrix2<
    Vector2< Fixed<0>, Imaginary >,
    Vector2< MinusImaginary, Fixed<0> > 
> sigma2;

//-----------------------------------------------------------------------------
Matrix2<
    Vector2< Fixed<0>, Float >,
    Vector2< Float, Fixed<0> > 
> m3(Unused(), 0.2f,
     0.8f, Unused());


// EDIT: Nicer initialization syntax in c++11

//-----------------------------------------------------------------------------
template <typename M00, typename M01, typename M10, typename M11>
Matrix2< Vector2<M00, M01>, Vector2<M10, M11> >
MakeMatrix(M00 const& m00, M01 const& m01, M10 const& m10, M11 const& m11)
{
    return Matrix2< Vector2<M00, M01>, Vector2<M10, M11> >(m00,m01,m10,m11);
}

auto m4 = MakeMatrix(Fixed<0>(),  Float(0.2f), 
                     Float(0.8f), Fixed<0>()  );

- Pete