haskell矩阵实现性能

Question

我听说过很多关于用Haskell编写的程序的惊人性能,并希望进行一些测试.因此,我为矩阵运算编写了一个"库",只是为了将它的性能与用纯C编写的相同内容进行比较.首先,我测试了500000个矩阵的乘法性能,并注意到它是......永无止境的(即结束于所以10分钟后出现内存异常)!在研究了haskell后,我设法摆脱了懒惰,我设法获得的最佳结果比C中的等效物慢20倍.所以,问题是:你能否查看下面的代码并告诉它的性能是否可以改善了一点？20次仍然令我失望.

import Prelude hiding (foldr, foldl, product)
import Data.Monoid
import Data.Foldable

import Text.Printf
import System.CPUTime

import System.Environment

data Vector a = Vec3 a a a
              | Vec4 a a a a
                deriving Show

instance Foldable Vector where
  foldMap f (Vec3 a b c) = f a `mappend` f b `mappend` f c
  foldMap f (Vec4 a b c d) = f a `mappend` f b `mappend` f c `mappend` f d

data Matr a = Matr !a !a !a !a
                   !a !a !a !a
                   !a !a !a !a
                   !a !a !a !a

instance Show a => Show (Matr a) where
  show m = foldr f [] $ matrRows m
            where f a b = show a ++ "\n" ++ b

matrCols (Matr a0 b0 c0 d0 a1 b1 c1 d1 a2 b2 c2 d2 a3 b3 c3 d3)
              = [Vec4 a0 a1 a2 a3, Vec4 b0 b1 b2 b3, Vec4 c0 c1 c2 c3, Vec4 d0 d1 d2 d3]

matrRows (Matr a0 b0 c0 d0 a1 b1 c1 d1 a2 b2 c2 d2 a3 b3 c3 d3)
              = [Vec4 a0 b0 c0 d0, Vec4 a1 b1 c1 d1, Vec4 a2 b2 c2 d2, Vec4 a3 b3 c3 d3]

matrFromList [a0, b0, c0, d0, a1, b1, c1, d1, a2, b2, c2, d2, a3, b3, c3, d3]
        = Matr a0 b0 c0 d0
               a1 b1 c1 d1
               a2 b2 c2 d2
               a3 b3 c3 d3

matrId :: Matr Double
matrId  = Matr 1 0 0 0
               0 1 0 0
               0 0 1 0
               0 0 0 1

normalise (Vec4 x y z w) = Vec4 (x/w) (y/w) (z/w) 1

mult a b = matrFromList [f r c | r <- matrRows a, c <- matrCols b] where 
            f a b = foldr (+) 0 $ zipWith (*) (toList a) (toList b)

Answer 1

首先,我怀疑你是否会通过这种实现获得出色的表现.不同表示之间的转换太多.你最好将代码放在像vector vector这样的东西上.此外,您不提供所有测试代码,因此可能还有其他问题我们无法在此处.这是因为生产对消费的管道对Haskell性能有很大影响,而且你没有提供任何一个结束.

现在,有两个具体问题:

1)您的向量被定义为3或4元素向量.这意味着对于每个向量,都需要额外检查以查看正在使用的元素数量.在C中,我想你的实现可能更接近

struct vec {
  double *vec;
  int length;
}

你应该在Haskell做类似的事情; 这是怎样vector和bytestring例如实现.

即使您不更改Vector定义,也要严格使用字段.您还应该添加UNPACK编译指示(到Vector和Matrix)或使用编译-funbox-strict-fields.

2)mult改为

mult a b = matrFromList [f r c | r <- matrRows a, c <- matrCols b] where 
            f a b = Data.List.foldl' (+) 0 $ zipWith (*) (toList a) (toList b)

foldl'在这种情况下,额外的严格性将提供更好的性能foldr.

仅这一变化可能会产生很大的不同,但是如果没有看到其余的代码,就很难说.

Answer 2

回答我自己的问题只是为了分享我昨天得到的新结果：

1. 我将 ghc 升级到最新版本，性能确实没有那么糟糕（只差了大约 7 倍）。

2. 我还尝试以一种愚蠢而简单的方式实现矩阵（请参见下面的列表），并获得了真正可接受的性能 - 仅比 C 语言慢大约 2 倍。

   data Matr a = Matr ( a, a, a, a
                      , a, a, a, a
                      , a, a, a, a
                      , a, a, a, a) 
   
   mult (Matr (!a0,  !b0,  !c0,  !d0,  
               !a1,  !b1,  !c1,  !d1,  
               !a2,  !b2,  !c2,  !d2,  
               !a3,  !b3,  !c3,  !d3))
        (Matr (!a0', !b0', !c0', !d0', 
               !a1', !b1', !c1', !d1', 
               !a2', !b2', !c2', !d2', 
               !a3', !b3', !c3', !d3'))
        = Matr ( a0'', b0'', c0'', d0''
               , a1'', b1'', c1'', d1''
               , a2'', b2'', c2'', d2''
               , a3'', b3'', c3'', d3'')
           where a0'' = a0 * a0' + b0 * a1' + c0 * a2' + d0 * a3'
                 b0'' = a0 * b0' + b0 * b1' + c0 * b2' + d0 * b3'
                 c0'' = a0 * c0' + b0 * c1' + c0 * c2' + d0 * c3'
                 d0'' = a0 * d0' + b0 * d1' + c0 * d2' + d0 * d3'
   
                 a1'' = a1 * a0' + b1 * a1' + c1 * a2' + d1 * a3'
                 b1'' = a1 * b0' + b1 * b1' + c1 * b2' + d1 * b3'
                 c1'' = a1 * c0' + b1 * c1' + c1 * c2' + d1 * c3'
                 d1'' = a1 * d0' + b1 * d1' + c1 * d2' + d1 * d3'
   
                 a2'' = a2 * a0' + b2 * a1' + c2 * a2' + d2 * a3'
                 b2'' = a2 * b0' + b2 * b1' + c2 * b2' + d2 * b3'
                 c2'' = a2 * c0' + b2 * c1' + c2 * c2' + d2 * c3'
                 d2'' = a2 * d0' + b2 * d1' + c2 * d2' + d2 * d3'
   
                 a3'' = a3 * a0' + b3 * a1' + c3 * a2' + d3 * a3'
                 b3'' = a3 * b0' + b3 * b1' + c3 * b2' + d3 * b3'
                 c3'' = a3 * c0' + b3 * c1' + c3 * c2' + d3 * c3'
                 d3'' = a3 * d0' + b3 * d1' + c3 * d2' + d3 * d3'