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

问：

我是 Haskell 的新手。 作为一门学术课程的一部分，我被要求在 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 的“一块”。它有效。拿一个更大的部分（700x700，或所有输入），它变得不切实际 - 大约一个小时的计算，由于内存使用过多，计算机完全卡住了。

感谢丹尼尔·瓦格纳（Daniel Wagner）指出： 对于想在家玩的人，请尝试：

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

代替从磁盘加载的矩阵。

我将不胜感激任何建议。