计算向量表达式的范数||aW+bX+cY||

Question

我是博士生。在我的论文介绍中，我对线性代数工具的表达性和性能之间的折衷很感兴趣。

作为一个简单的例子，我使用向量表达式范数的计算。我的例子的C代码是:

float normExpression3(float a, float *W, float b, float *X, float c, float*Y){
double norm = 0;
for (int i=0; i<n; ++i) // n in [3e6; 2e8]
{
    float tmp = a*W[i]+b*X[i]+c*Y[i];
    norm+=tmp*tmp;
}
return sqrtf(norm);

}

我比较了用不同的技术取得的成绩。由于向量很大(数百万个元素)，性能受到内存带宽的限制。然而，不同的方法之间存在着巨大的差异。

我写的优化的C版本没有表现力(一个新的函数必须被写为第四个向量)，而且非常难看(线程化和矢量化)，但实现了6.4 GFlops。另一方面，MATLAB代码非常好:

result = norm(a*W+b*X+c*Y)

但只能达到0.28 GFlops。

c++模板表达式 la Blitz++为用户提供了表现力和性能(6.5 GFlops)。

作为我分析的一部分，我想知道函数式语言与这些方法相比如何。我想在Haskell或OCaml中展示一个例子(我敢说，两者都被认为非常适合这种操作)。

这些语言我都不懂。我可以学习其中的一个来提供我的例子，但这不是一个公平的比较:我不确定是否能够提供一个同时允许表现力和性能的实现。

我的两个问题是:哪种语言最适合?2)如何有效地计算向量表达式的范数而不影响实现的通用性?

提前，谢谢!

k .

右舵

---

编辑:修正float的norm蓄电池类型为double

Answer 1

值得注意的是，下面是您的函数的OCaml版本:

let normExpression3 a w b x c y =
    let n = Array.length w in
    if not (n = Array.length x && n = Array.length y)
        then invalid_arg "normExpression3";
    let (@) = (Array.unsafe_get : float array -> int -> float) in
    let rec accum a w b x c y n i norm =
        if i = n then sqrt norm else
        let t = a *. (w @ i) +. b *. (x @ i) +. c *. (y @ i) in
        accum a w b x c y n (i + 1) (norm +. t)
    in accum a w b x c y n 0 0.

对性能做了一些考虑，即:

• 未检查的数组访问(或者更确切地说，手动提升到循环之外的数组边界检查)
• 单态数组访问
• 避免浮点累加器装箱和拆箱的递归内循环
• 内循环的lambda提升以避免引用闭合值

最后一个优化应该检查一个封闭的内部循环，因为有这么多参数，寄存器溢出可能会超过引用封闭参数的成本。

请注意，除非在基准;-)中竞争，否则通常不会为这种优化而费心，还请注意，您将需要使用64位OCaml进行测试，因为数组被限制为4兆元素。

Answer 2

1)哪种语言最适合?

都可用于此类任务。主要关注的是所需库的可用性(例如，用于向量或矩阵的库)，以及是否需要并行性。

Haskell中的vector和repa等库非常适合。在repa的情况下，你也得到了并行性。

2)如何有效地计算向量表达式的范数而不影响实现的通用性?

一种方法是使用元编程技术从高级描述生成计算内核的专门实现。在函数式语言中，这是一种相对常见的技术，基于小型特定于领域的语言，使用自定义代码生成器。

参见《高性能蒙特卡罗方法专用模拟器生成器》或《FFTW上的OCaml工作》。

Answer 3

不是函数式语言本身的答案，但请注意，计算norm的实现(在C中)与matlab实际计算的方式不同。

是的，确实有很好的理由。很可能你对norm的近似对于任何实际用例来说都是毫无用处的(因为它目前是如何实现的)。请不要低估计算norm s的数值近似的"困难"。

Answer 4

正如don所说，考虑Repa。下面是一些简单的代码来帮助你开始。

import Data.Array.Repa as R
len :: Int
len = 50000
main = do
    let ws = R.fromList (Z :. len) [0..len-1]
        xs = R.fromList (Z :. len) [10498..10498 + len - 1]
        ys = R.fromList (Z :. len) [8422..8422 + len - 1]
    print (multSum 52 73 81 ws xs ys)
multSum a b c ws xs ys = R.map (a*) ws +^ R.map (b*) xs +^ R.map (c*) ys

这仍然让您找到一种从磁盘获取数据并将其放入Repa数组的好方法。我认为将其全部读取为Llazy ByteString并使用Repa.fromFunction应该做的，也许有人会以更聪明的方式加入