Ole*_*828 19 performance haskell functional-programming
我知道Haskell一点点,我想知道是否有可能在Haskell中编写类似于矩阵矩阵产品的东西,它具有以下所有特性:
IO或Statemonad(我不关心函数体中发生了什么.也就是说,我不关心函数体是否使用monad,只要整个函数是纯的).我可能想在纯函数中使用这个矩阵矩阵乘积.float[]非float[][])和大小为1000x1000的矩阵,以及一个单核CPU.(如果你每个周期得到0.5-2个浮点运算,你可能会在球场.)(我不希望这听起来像是一个挑战,但请注意Java可以轻松满足以上所有要求.)
我已经知道了
致意见:
我有三个原因:首先,"没有malloc或指针"的要求尚未定义(我要求你编写任何不使用指针的Haskell代码);
我看到很多Haskell程序没有使用Ptr.也许它指的是在机器指令级别使用指针的事实？那不是我要表达的意思.我指的是Haskell源代码的抽象级别.
第二,对CS研究的攻击是不合适的(而且我无法想象比使用别人已经为你编写的代码更简单的事情); 第三,Hackage上有很多矩阵包(要求这个问题的准备工作应该包括审查和拒绝每个).
你的#2和#3似乎是相同的("使用现有的库").我对矩阵产品感兴趣,只是测试Haskell可以自己做什么,以及它是否允许你"简单易懂".我可以很容易地想出一个没有任何现成库的数值问题,但是我必须解释这个问题,而每个人都已经知道矩阵产品是什么.
Java如何可能满足1.？任何Java方法本质上都是
:: IORef Arg -> ... -> IORef This -> IO Ret
这是我问题的根源,实际上是(+1).虽然Java并没有声称追踪纯度,但Haskell确实如此.在Java中,注释中指出了函数是否纯粹.我可以声称矩阵产品是纯的,即使我在函数体中做了突变.问题在于Haskell的方法(在类型系统中编码的纯度)是否与效率,内存安全性和简单性兼容.
Dan*_*her 19
和Java一样高效.具体来说,我们假设我在谈论一个简单的三重循环,单精度,连续的列主要布局(float [],而不是float [] [])和大小为1000x1000的矩阵,以及一个单核CPU.(如果你每个周期得到0.5-2个浮点运算,你可能在球场)
所以像
public class MatrixProd {
static float[] matProd(float[] a, int ra, int ca, float[] b, int rb, int cb) {
if (ca != rb) {
throw new IllegalArgumentException("Matrices not fitting");
}
float[] c = new float[ra*cb];
for(int i = 0; i < ra; ++i) {
for(int j = 0; j < cb; ++j) {
float sum = 0;
for(int k = 0; k < ca; ++k) {
sum += a[i*ca+k]*b[k*cb+j];
}
c[i*cb+j] = sum;
}
}
return c;
}
static float[] mkMat(int rs, int cs, float x, float d) {
float[] arr = new float[rs*cs];
for(int i = 0; i < rs; ++i) {
for(int j = 0; j < cs; ++j) {
arr[i*cs+j] = x;
x += d;
}
}
return arr;
}
public static void main(String[] args) {
int sz = 100;
float strt = -32, del = 0.0625f;
if (args.length > 0) {
sz = Integer.parseInt(args[0]);
}
if (args.length > 1) {
strt = Float.parseFloat(args[1]);
}
if (args.length > 2) {
del = Float.parseFloat(args[2]);
}
float[] a = mkMat(sz,sz,strt,del);
float[] b = mkMat(sz,sz,strt-16,del);
System.out.println(a[sz*sz-1]);
System.out.println(b[sz*sz-1]);
long t0 = System.currentTimeMillis();
float[] c = matProd(a,sz,sz,b,sz,sz);
System.out.println(c[sz*sz-1]);
long t1 = System.currentTimeMillis();
double dur = (t1-t0)*1e-3;
System.out.println(dur);
}
}
我想？(我没有在编码之前正确阅读规范,所以布局是行主要的,但由于访问模式是相同的,因为混合布局没有区别,所以我认为没关系.)
我没有花时间去思考一个聪明的算法或低级优化技巧(无论如何我在Java中都没有取得多少成就).我刚写了简单的循环,因为
我不希望这听起来像是一个挑战,但请注意Java可以轻松满足上述所有要求
这就是Java 轻松提供的东西,所以我会接受它.
(如果你每个周期得到0.5-2个浮点运算,你可能在球场)
无论如何,我都害怕,无论是在Java还是在Haskell.使用简单的三重循环,太多缓存未命中以达到该吞吐量.
在Haskell中做同样的事情,再一次不考虑聪明,简单明了的三重循环:
{-# LANGUAGE BangPatterns #-}
module MatProd where
import Data.Array.ST
import Data.Array.Unboxed
matProd :: UArray Int Float -> Int -> Int -> UArray Int Float -> Int -> Int -> UArray Int Float
matProd a ra ca b rb cb =
let (al,ah) = bounds a
(bl,bh) = bounds b
{-# INLINE getA #-}
getA i j = a!(i*ca + j)
{-# INLINE getB #-}
getB i j = b!(i*cb + j)
{-# INLINE idx #-}
idx i j = i*cb + j
in if al /= 0 || ah+1 /= ra*ca || bl /= 0 || bh+1 /= rb*cb || ca /= rb
then error $ "Matrices not fitting: " ++ show (ra,ca,al,ah,rb,cb,bl,bh)
else runSTUArray $ do
arr <- newArray (0,ra*cb-1) 0
let outer i j
| ra <= i = return arr
| cb <= j = outer (i+1) 0
| otherwise = do
!x <- inner i j 0 0
writeArray arr (idx i j) x
outer i (j+1)
inner i j k !y
| ca <= k = return y
| otherwise = inner i j (k+1) (y + getA i k * getB k j)
outer 0 0
mkMat :: Int -> Int -> Float -> Float -> UArray Int Float
mkMat rs cs x d = runSTUArray $ do
let !r = rs - 1
!c = cs - 1
{-# INLINE idx #-}
idx i j = cs*i + j
arr <- newArray (0,rs*cs-1) 0
let outer i j y
| r < i = return arr
| c < j = outer (i+1) 0 y
| otherwise = do
writeArray arr (idx i j) y
outer i (j+1) (y + d)
outer 0 0 x
和调用模块
module Main (main) where
import System.Environment (getArgs)
import Data.Array.Unboxed
import System.CPUTime
import Text.Printf
import MatProd
main :: IO ()
main = do
args <- getArgs
let (sz, strt, del) = case args of
(a:b:c:_) -> (read a, read b, read c)
(a:b:_) -> (read a, read b, 0.0625)
(a:_) -> (read a, -32, 0.0625)
_ -> (100, -32, 0.0625)
a = mkMat sz sz strt del
b = mkMat sz sz (strt - 16) del
print (a!(sz*sz-1))
print (b!(sz*sz-1))
t0 <- getCPUTime
let c = matProd a sz sz b sz sz
print $ c!(sz*sz-1)
t1 <- getCPUTime
printf "%.6f\n" (fromInteger (t1-t0)*1e-12 :: Double)
所以我们在两种语言中都做了几乎完全相同的事情.-O2使用javac 编译带有Java 的Haskell
$ java MatrixProd 1000 "-13.7" 0.013
12915.623
12899.999
8.3592897E10
8.193
$ ./vmmult 1000 "-13.7" 0.013
12915.623
12899.999
8.35929e10
8.558699
由此产生的时间非常接近.
如果我们将Java代码编译为native,with gcj -O3 -Wall -Wextra --main=MatrixProd -fno-bounds-check -fno-store-check -o jmatProd MatrixProd.java,
$ ./jmatProd 1000 "-13.7" 0.013
12915.623
12899.999
8.3592896512E10
8.215
仍然没有太大的区别.
作为特殊奖励,C(gcc -O3)中的算法相同:
$ ./cmatProd 1000 "-13.7" 0.013
12915.623047
12899.999023
8.35929e+10
8.079759
因此,当涉及使用浮点数的计算密集型任务时,这表明简单的Java和简单的Haskell之间没有根本的区别(当处理中到大数字的整数运算时,GHC使用GMP使得Haskell的表现远远超过Java的BigInteger对于许多任务,但这当然是库问题,而不是语言问题),并且两者都与此算法接近于C.
但公平地说,这是因为访问模式会导致每隔一纳秒出现一次缓存错误,因此在所有三种语言中,这种计算都受内存限制.
如果我们通过将行主矩阵与列主矩阵相乘来改进访问模式,所有都变得更快,gcc编译的C完成它1.18s,java需要1.23s而ghc编译的Haskell需要大约5.8s,通过使用llvm后端可以减少到3秒.
在这里,数组库的范围检查真的很痛.使用未经检查的数组访问(在检查错误之后,由于检查已经在控制循环的代码中完成),GHC的本机后端在2.4s内完成,通过llvm后端让计算在1.55s完成,这是不错的,虽然明显慢于C和Java.使用来自GHC.Prim而不是数组库的原语,llvm后端产生在1.16s中运行的代码(同样,没有对每次访问进行边界检查,但是在计算过程中只生成有效的索引,在这种情况下可以很容易地证明之前,所以在这里,没有牺牲内存安全性¹;检查每次访问会使时间达到1.96秒,仍然明显优于阵列库的边界检查).
底线:GHC需要(更多)更快的分支进行边界检查,并且优化器有改进的余地,但原则上,"Haskell的方法(在类型系统中编码的纯度)与效率,内存安全性和简单性兼容",我们还没有.目前,人们必须决定一个人愿意牺牲多少.
¹是的,这是一个特例,通常省略边界检查会牺牲记忆安全性,或者至少更难证明它没有.
Die*_*Epp 11
攻击这个问题有两个角度.
沿着这些方向进行的研究正在进行中.现在,有很多Haskell程序员比我更聪明; 事实上,我不断被提醒和谦卑.其中一个可能会来并纠正我,但我不知道有任何简单的方法将安全的Haskell原语组合成一个顶级的矩阵乘法例程.你谈到的那些论文听起来像是一个好的开始.
但是,我不是计算机科学研究员.我想知道在Haskell中是否可以保持简单的事情.
如果你引用那些论文,也许我们可以帮助破译它们.
按照这些思路,软件工程很容易理解,简单,甚至容易.一个精明的Haskell编码器会在BLAS周围使用一个薄的包装器,或者在Hackage中寻找这样的包装器.
解读尖端研究是一个持续的过程,将知识从研究人员转移到工程师.计算机科学研究员CAR Hoare首先发现了quicksort并发表了一篇关于它的论文.今天,它是一个罕见的计算机科学毕业生,无法亲自实现记忆中的快速排序(至少,最近毕业的那些).
几乎这个问题在历史上曾被问过几次.
是否有可能在Fortran中编写与汇编一样快的矩阵运算？
是否有可能在C中编写与Fortran一样快的矩阵运算？
是否有可能在Java中编写与C一样快的矩阵运算？
是否有可能在Haskell中编写与Java一样快的矩阵运算？
到目前为止,答案一直是"尚未",其次是"足够接近".使这成为可能的进步来自编写代码的改进,编译器的改进以及编程语言本身的改进.
作为一个具体的例子,在C99编译器在过去十年中普及之前,C在许多实际应用中无法超越Fortran.在Fortran中,假设不同的数组彼此具有不同的存储,而在C中,通常不是这种情况.因此,允许Fortran编译器对C编译器进行优化.好吧,直到C99出来并且您可以将restrict限定符添加到您的代码中.
Fortran编译器等待.最终处理器变得非常复杂,以至于良好的汇编编写变得更加困难,并且编译器变得足够复杂以至于Fortran很快.
然后C程序员一直等到2000年才能编写与Fortran匹配的代码.在那之前,他们使用Fortran或汇编程序(或两者)编写的库,或者降低速度.
同样,Java程序员必须等待JIT编译器,并且必须等待特定的优化才能出现.JIT编译器最初是一个深奥的研究概念,直到它们成为日常生活的一部分.为了避免每个阵列访问的测试和分支,边界检查优化也是必要的.
因此,很明显Haskell程序员正在"等待",就像之前的Java,C和Fortran程序员一样.我们还在等什么？
也许我们只是在等待某人编写代码,并向我们展示它是如何完成的.
也许我们正在等待编译器变得更好.
也许我们正在等待Haskell语言本身的更新.
也许我们正在等待上述的一些组合.
纯度和单子在Haskell中得到很多混淆.原因是因为在Haskell中,不纯的函数总是使用IOmonad.例如,Statemonad是100%纯.所以当你说"纯粹"和"类型签名不使用Statemonad"时,那些实际上是完全独立和单独的要求.
但是,你也可以IO在实现纯函数时使用monad,事实上,它很容易:
addSix :: Int -> Int
addSix n = unsafePerformIO $ return (n + 6)
好的,是的,这是一个愚蠢的功能,但它是纯粹的.它甚至显然是纯粹的.纯度测试有两个方面:
它是否为相同的输入提供相同的结果？是.
它是否会产生任何语义上重要的副作用？没有.
我们喜欢纯度的原因是因为纯函数比不纯函数更容易编写和操作.它们的实施方式并不重要.我不知道你是否意识到这一点,但Integer和ByteString周围不纯的C函数都基本包装,尽管界面是纯粹的.(有关新实现的工作Integer,我不知道它有多远.)
问题在于Haskell的方法(在类型系统中编码的纯度)是否与效率,内存安全性和简单性兼容.
该部分的答案是"是",因为我们可以从BLAS中获取简单的函数并将它们放在一个纯粹的,类型安全的包装器中.包装器的类型编码函数的安全性,即使Haskell编译器无法证明函数的实现是纯粹的.我们unsafePerformIO在其实现中的使用既承认我们已经证明了函数的纯洁性,也是一种让步,我们无法想出在Haskell类型系统中表达该证明的方法.
但答案也"尚未",因为我不知道如何在Haskell中完全实现该函数.
该领域的研究正在进行中.人们正在研究像Coq这样的证明系统和像Agda这样的新语言,以及GHC本身的发展.为了看看我们需要什么样的类型系统来证明可以安全地使用高性能BLAS例程.这些工具也可以与Java等其他语言一起使用.例如,您可以在Coq中编写一个证明您的Java实现是纯粹的.
我为"是和否"答案道歉,但没有其他答案会认可工程师(关心"是")和研究人员(关心"尚未")的贡献.
PS请引用论文.
像Java一样,Haskell不是编写数字代码的最佳语言.
Haskell的数字密码代码是......平均值.英特尔和海湾合作委员会都没有经过多年的研究.
相反,Haskell为您提供了一种方法,可以将"快速"代码与应用程序的其余部分完美地连接起来.请记住,3%的代码负责97%的应用程序运行时间.1
使用Haskell,您可以通过非常好的C外部函数接口以一种非常好地与其他代码进行交互的方式调用这些高度优化的函数.实际上,如果您愿意,可以使用体系结构的汇编语言编写数字代码,从而获得更高的性能!在应用程序的性能较重的部分插入C不是一个错误 - 这是一个功能.
但我离题了.
通过隔离这些高度优化的函数,并使用与其余Haskell代码类似的接口,您可以使用Haskell非常强大的重写规则执行高级优化,这允许您编写规则,例如reverse . reverse == id在编译时自动减少复杂表达式2.这导致了极其快速,纯粹功能和易于使用的库,如Data.Text 3和Data.Vector [4].
通过结合高级和低级优化,我们最终得到了更优化的实现,每一半("C/asm"和"Haskell")相对容易阅读.低级优化以其原生语言(C或汇编语言)完成,高级优化获得特殊的DSL(Haskell重写规则),其余代码完全无视它.
总而言之,是的,Haskell可以比Java更快.但它通过C为原始FLOPS作弊.这在Java中要难得多(而且Java的FFI的开销要高得多),因此可以避免.在Haskell中,它很自然.如果您的应用程序花费了大量时间进行数值计算,那么可能不是查看Haskell或Java,而是查看Fortran以满足您的需求.如果您的应用程序将大部分时间花在性能敏感代码的一小部分上,那么Haskell FFI是您最好的选择.如果您的应用程序没有花费任何时间在数字代码...然后使用你喜欢的任何.=)
Haskell(也不是Java,就此而言)不是Fortran.
1这些数字已经弥补,但你明白我的观点.
2 http://www.cse.unsw.edu.au/~dons/papers/CLS07.html
3 http://hackage.haskell.org/package/text
[4] http://hackage.haskell.org/package/vector
现在已经不在了,回答你的实际问题:
不,在Haskell中编写矩阵乘法目前并不聪明.目前,REPA是规范的方式[5].实现部分破坏了内存安全性(他们使用unsafeSlice),但"破坏内存安全"与该功能隔离,实际上非常安全(但编译器不易验证),并且如果出现问题则易于删除(替换" unsafeSlice"with"slice").
但这是Haskell!很少有功能的性能特征是孤立的.这可能是一件坏事(在空间泄漏的情况下),或者非常非常好的事情.
尽管使用的矩阵乘法算法是天真的,但它在原始基准测试中表现更差.但我们的代码很少看起来像基准测试.
如果你是一个拥有数百万个数据点的科学家想要繁殖庞大的矩阵怎么办？[7]
对于那些人,我们有mmultP [6].这执行矩阵乘法,但是是数据并行的,并且受REPA嵌套数据并行性的影响.另请注意,代码与顺序版本基本不变.
对于那些不会使大矩阵相乘,而是乘以大量小矩阵的人,往往会有其他代码与所述矩阵交互.可能将其切割成列向量并找到它们的点积,可能找到它的特征值,也许完全是其他东西.与C不同,Haskell知道虽然你喜欢孤立地解决问题,但通常在那里找不到最有效的解决方案.
像ByteString,Text和Vector一样,REPA数组也需要融合.2顺便说一下你应该读2 - 这是一篇写得很好的论文.这与相关代码的积极内联和REPA高度并行的特性相结合,使我们能够在幕后用非常先进的高级优化来表达这些高级数学概念.
虽然目前还不知道在纯函数式语言中编写有效矩阵乘法的方法,但我们可以稍微接近(没有自动矢量化,一些过多的解引用来获取实际数据等),但没有什么比IFORT或GCC可以做到.但是在岛上不存在程序,并且使得整个岛屿在Haskell上比Java更容易表现得更好.
[7]实际上,最好的方法是使用GPU.有一些可用于Haskell的GPU DSL可以实现本地化.他们真的很整洁!