Haskell 实现行列式、秩和逆矩阵计算 - 输入矩阵大小限制

Question

我是 Haskell 新手。 作为学术类(class)的一部分，我被要求在 Haskell 中实现一个函数，用于计算给定矩阵的行列式、秩和逆矩阵。

我使用高斯消除(对原始矩阵和另一个初始化为单位矩阵的矩阵执行相同的行运算)。我使用一次性方法“即时”累积排名和决定因素。

对于高达 600x600 的输入矩阵，该函数的行为符合预期(本例中的执行时间为 63 秒)。 任何超过该大小的输入都会导致过多的内存使用和不切实际的计算时间。

输入是 3 个不同的矩阵，大小分别为 800x800、800x800 和 1000x1000。

值得一提的是，我不允许使用任何外部库(例如 HMatrix 或 Data.Matrix)

我尽了最大努力来优化代码，但由于我是新手，所以我未能成功地使我的代码适用于大于 600x600 的尺寸。

import Data.Time.Clock
import System.IO

type Matrix = [[Double]]

-- Row operations

swapRows :: Int -> Int -> Matrix -> Matrix
swapRows i j mat = take i mat ++ [mat !! j] ++ (drop (i+1) . take j $ mat) ++ [mat !! i] ++ drop (j+1) mat

scaleRow :: Int -> Double -> Matrix -> Matrix
scaleRow i k mat = take i mat ++ [map (* k) (mat !! i)] ++ drop (i+1) mat

addRows :: Int -> Int -> Double -> Matrix -> Matrix
addRows i j k mat = take i mat ++ [zipWith (+) (mat !! i) (map (* k) (mat !! j))] ++ drop (i+1) mat

-- Gaussian elimination

eliminate :: Matrix -> (Matrix, Matrix, Double, Int)
eliminate mat = foldr eliminateRow (mat, identityMatrix (length mat), 1.0, rank) [0..length mat-1]
  where
    rank = length mat  -- Initialize rank as the number of rows

    eliminateRow row (mat, invMat, detAccum, rank) = foldl eliminateEntry (mat, invMat, detAccum, rank) [0..length mat-1]
      where
        eliminateEntry (mat, invMat, detAccum, rank) col
          | row == col = (scaleRow row (1 / pivot) mat, scaleRow row (1 / pivot) invMat, detAccum * pivot, rank)
          | mat !! col !! row /= 0 = (addRows col row (-ratio) mat, addRows col row (-ratio) invMat, detAccum, rank)
          | otherwise = (mat, invMat, detAccum, rank - 1)
          where
            pivot = mat !! row !! row
            ratio = (mat !! col !! row) / pivot

-- Matrix operations

identityMatrix :: Int -> Matrix
identityMatrix n = [[fromIntegral (fromEnum (i == j)) | j <- [1..n]] | i <- [1..n]]

-- create sub matrix n x n from matrix mat
subMatrix :: [[Double]] -> Int -> [[Double]]
subMatrix mat n = take n (map (take n) mat) 

-- Testing

main :: IO ()
main = do
  --let submat = [[1, 2], [3, 4]]  -- 3x3 invertible matrix
  
  contenm_headMulScal <- readFile "det_matrix(800 x 800).txt"
  let mat = map (map read . words) (lines contenm_headMulScal) :: [[Double]]
  
  let piece_dim = 600
  let submat = subMatrix mat piece_dim
  
  tt_start <- getCurrentTime
  let (refMat, inverse, determinant, rank) = eliminate submat  -- Calculate the ref matrix, determinant, and rank
  t_end <- getCurrentTime
  --putStrLn "Original Matrix:"
  --printMatrix submat
  putStrLn "\nDeterminant:"
  print determinant
  putStrLn "\nInverse Matrix:"
  --printMatrix inverse
  putStrLn $ "First element in the inverse matrix: " ++ show (head (head inverse))
  putStrLn "\nRank:"
  print rank
  tt_end <- getCurrentTime
  
  print "Runtime:"
  print (diffUTCTime tt_end tt_start)

printMatrix :: Matrix -> IO ()
printMatrix mat = mapM_ (putStrLn . unwords . map show) mat

例如，如您所见，我从 800x800 输入中取出了 600x600 的“一 block ”。有用。采取更大的一 block (700x700，或所有输入)，它变得不切实际 - 大约一个小时的计算，计算机由于内存使用过多而完全卡住。

感谢 Daniel Wagner 指出: 对于想在家玩的人，请尝试:

countingMatrix :: Int -> Matrix
countingMatrix n = [[fromIntegral (j*n+i) | j <- [1..n]] | i <- 
[1..n]] 

代替从磁盘加载的矩阵。

如果有任何建议，我将不胜感激。

Answer 1

似乎存在空间泄漏/惰性问题。我不怪你:我发现获得我想要的评估行为非常挑剔，即使我已经很确定问题是什么!

您要确保的是，在创建新矩阵时，您不会保留对旧矩阵的任何引用。如果这样做，那么 GC 将无法收集旧矩阵。特别是，引用包括潜在计算，例如 take i oldMatrix 或 oldMatrix !!我或类似的。所以!让我们讨论一下为了实现这一目标我必须做出的改变。

首先:当我们创建一个新矩阵时，它主要是旧矩阵中行的副本，但有一两个新行，我们需要一种方式来表示“遍历整个矩阵，并强制进行足够的计算”您将指针复制到特定行，而不是稍后在旧矩阵中查找该行的计算”。为了实现这一点，我们将创建一个函数来强制列表中的每个元素，但仅限于 WHNF。请注意，由于我们传递的是一个矩阵，因此这不是完整的评估!我们不会一直强制到矩阵元素，只强制到矩阵行的级别。在我第一次尝试时，我就犯了这个错误，一直强制到元素级别，并无意中引入了非常严重的时间性能回归。

seqElements :: [a] -> [a]
seqElements [] = []
seqElements as@(a:at) = a `seq` seqElements at `seq` as

我们将在每个行操作开始时调用它。

swapRows i j mat = seqElements $ {- old implementation -}
scaleRows i k mat = seqElements $ {- old implementation -}
addRows i j k mat = seqElements $ {- old implementatiotn -}

这构成了我们手动强制注释的“基本情况”。不幸的是，我们现在需要将其一直传播到表达式树——每个调用者需要确保它使用字段中的其中一个数据结构创建的任何数据结构也将该数据结构的评估与该字段的评估联系起来。在eliminateEntry中，这意味着我们需要更严格的四元组版本。

quadruple a b c d = a `seq` b `seq` c `seq` d `seq` (a, b, c, d)
-- then, later, we replace (,,,) with quadruple:
        eliminateEntry (mat, invMat, detAccum, rank) col
          | row == col = quadruple (scaleRow row (1 / pivot) mat) (scaleRow row (1 / pivot) invMat) (detAccum * pivot) rank
          | mat !! col !! row /= 0 = quadruple (addRows col row (-ratio) mat) (addRows col row (-ratio) invMat) detAccum rank
          | otherwise = quadruple mat invMat detAccum (rank - 1)

最后，在 eliminateRow 中，我建议从 foldl 切换到 foldl'。在我的测试中，这似乎没有什么区别，但这是一个值得养成的好习惯； foldl 提供的额外惰性几乎是不需要的，而且常常是有害的。

    eliminateRow row (mat, invMat, detAccum, rank) = foldl' eliminateEntry (mat, invMat, detAccum, rank) [0..length mat-1]

通过这些更改，我现在发现内存使用量要低得多:600x600 矩阵约为 50 秒、96MiB，800x800 矩阵约为 128 秒、167MiB。

如果需要的话，可能有机会进行进一步的重大优化；例如我预计切换到矩阵的可变数组表示将是另一个巨大的提升。