C# 中的高效矩阵编程

Emanuele Ruffaldi

4.89/5 (44投票s)

2002年12月12日

4分钟阅读

460393

10949

快速矩阵表达式求值，基于动态代码生成和部分求值

引言

本文介绍了一种在 C# 语言中编写清晰高效矩阵数学代码的技术。矩阵运算通过操作符重定义来表示，但代码是动态生成的，遵循部分求值的原则。这提供了一种像 C++ 中表达式模板那样高效而优雅的解决方案。本文是我在 C# 动态代码生成和操作符重载方面研究的延续（此处）。

问题所在

像 C# 这样的面向对象语言，可以定义非常优雅的矩阵类，它们封装了元素访问，并且通过操作符重载机制，你可以轻松地编写

Matrix m = new Matrix(10,5);
Matrix n = new Matrix(10,5);
Matrix r;
double d = 2;	
	
r = m + n;
r = n*2;

这些类的经典实现为每个子表达式创建一个临时的 Matrix 对象，并且在下面的代码中，当生成两个临时对象时，问题就更加明显了。

Matrix a,b,c,r; // with m and n initialized

r = a+3*(b+c);

这种方法给基于对象的矩阵代码带来了很大的开销，通常数学程序员会使用数组和索引，从而失去了所有的面向对象特性。这类代码可以称为 C#Tran（来自 C#/Fortran），就像 Veldhuizen 提出的 JavaTran 一样。

double [] a, b, c, r;

for(int i = 0; i < a.length; i++)
    r[i] = a[i]+3*(b[i]+c[i]);

解决方案

解决这个问题的一个可能方法是延迟求值：当我们编写 b+c 时，我们生成一个 MatrixAdd 对象来封装该操作。它也重载了操作符，我们可以编写 3*(b+c) 来生成一个 MatrixAddMul 对象，该对象高效地求值该操作。这次的开销与所有可能的组合操作有关，并且如果我们的特定表达式没有被延迟操作覆盖，我们就会回退到低效的基于临时对象的代码。

我的提议是使用一种“生成式编程”方法：为特定操作生成专门的代码。我利用 .NET 框架的动态代码能力，在运行时直接在 IL 中生成基于 for 循环的代码。编译后的表达式可以被求值一次，或者存储到类型为 Evaluator 的委托中以供重用。

以下是使用示例

Vector a = new Vector(3);
Vector b = new Vector(3);
Vector r = new Vector(3);
	
(a+b*3).AssignTo(r)();

生成的 IL 代码等同于以下内容。您应该注意到 n 的大小在生成时已知，因此 for 循环的值是常数：这就是部分求值，因为代码是针对特定值专门生成的，JIT 可以生成优化的代码。

int n = a.rows();
for(int i = 0; i < <<n>>; i++)
{
    r[i] = a[i]+b[i]*3;
}

这是另一个使用 Matrix 对象的示例

Matrix a, b, c, r;
	
// evaluate directly
(a + Expr.Cos(b*c*3)).AssignTo(r)();

// store and reuse the compiled expression ...
Evaluator ev = (a + Expr.Cos(b*c*3)).AssignTo(r);
	
// and then evaluate
ev();

这是等效代码

int n = a.rows();
int m = a.cols();

for(int i = 0; i < <<n>>; i++)
    for (int j = 0; j < <<m>>; j++)
    {
        double d = 0;
        for(int k = 0; k < <<q>>; k++)
            d += b[i,k]*q[k,j];
        r[i,j] = a[i,j]+cos(d*3);
    }

实现

在内部，每个表达式都存储为派生自 Expr 类的对象的树，并分配给派生自 LeftExpr 的对象。我们可以缓慢地求值表达式，或者动态生成求值表达式的 IL 代码。我们考虑了向量和二维矩阵的操作，但该系统可以扩展到更多维度。在生成代码之前，我们检查操作中涉及的矩阵的大小，当它们不正确时，抛出 SizeMismatchException 异常。

类为 Expr 的表达式具有以下虚方法

// evaluates the expression against the two indices
float Eval(int i, int j)	
// compile the expression using the context
void  Compile(ILGenerator g, CompilerContext cc)
// a property to find the size of the expression
OpSize Size { get }

LeftExpr 是一个专门的 Expr，用于处理赋值

// assign a value v to the Expression
void Assign(int i,  int j, float v);
// compile the assignment (two phase compilation)
void CompileAssign(ILGenerator g, CompilerContext cc, bool post);

我定义的类的继承结构如下图所示

动态生成的类比上一篇文章中描述的类更复杂，这次我们需要存储所有数组对象、参数值和索引 i、j。要分析生成的代码，请在生成第一个表达式之前在 Compiler 类中设置 SafeMode 变量，然后在最后使用 Compiler.Save() 生成程序集。

子矩阵访问

使用此系统获得的功能之一是快速访问子矩阵。您可以对矩阵的各个部分（如行、列或方阵的主对角线）进行表达式求值。由于动态编译，对矩阵不同部分的访问非常高效。我还使用了易于阅读的表示法（我认为），基于 Tag 类方法。

现在的代码比标准的数组操作更清晰。以下代码显示了一些子矩阵访问及其 MATLAB 等效项（请记住 MATLAB 使用 1 基索引）。

Matrix m;

// a row, like MATLAB m(3,:)
m[2][Matrix.All]
	
// a range of rows like m(3:4, :)
m[Matrix.Range(2,3)][Matrix.All]
	
// a column, like m(:,5)
m[Matrix.All][4]
	
// columns, like m(:,3:6)
m[Matrix.Range(2,5)[Matrix.All]
	
// a submatrix m(2:3,4:5)
m[Matrix.Range(1,2)][Matrix.Range(3,4)]
	
// a single element m(3,4)
m[2][3]
	
// the third element taking the matrix as a row vector
// there is no MATLAB equivalent because MATLAB 
//uses column ordered matrices
m[2]
	
// some elements taking the matrix as a row vector
m[Matrix.Range(2,3)]
	
// the diagonal
m.Diagonal()

这些操作基于 MatrixRowRef 结构，该结构存储第一个数组访问并生成一个 SubMatrix 对象。它是一个派生自 LeftExpr 的类，因此可以在赋值中使用。

Vector v = new Vector(3);
Matrix m = new Matrix(3,5);
(v + 22).AssignTo(m[Matrix.All][2])();

速度考量

此系统生成的代码与基于 for 的代码一样高效，可以提供高速数学求值，但在生成阶段会有一点开销，因此当表达式被多次求值时，这种方法才有效。特别是第一个动态求值的表达式比其他表达式稍慢，因为创建了动态程序集。

结论

.NET 的动态代码生成，结合 C# 的自定义操作符，让我们能够构建优雅而高效的数学表达式，从而可以将 C# 用作数学语言。引用 Jim Blinn 的话，“我们必须在类库中尽可能多地融入智能，以便让这些类的用户尽可能轻松”。此代码的最终用户只需要知道它快速且易于使用。

参考文献

表达式模板，Todd Veldhuizen，1994（使用 C++ 模板的原始文章）
Java 中的表达式模板，Todd Veldhuizen（Java 的可能解决方案，使用完整的 Java 编译器）
优化 C++ 向量表达式，Jim Blinn，2000，IEEE Cg（伟大 Blinn 的精彩解释）