计算矩阵的Haskell行列式

我决定发布用于计算行列式的代码。该代码可以正常工作,尽管它不会伪装成专家。在Haskell上解决这个问题很有趣。考虑了两种解决问题的方法:简单递归和高斯方法。



一点理论



如您所知,平方矩阵n * n的行列式n!项,每个项都是一个乘积,其中每列恰好包含一个矩阵元素,而每行恰好包含一个矩阵元素。下一块的标志:

一种1个一世1个一种2一世2.........一种ñ一世ñ



由替换的奇偶性决定:

\ begin {pmatrix} 1&2&…&n \\ {i} _ {1}&{i} _ {2}&…&{i} _ {n} \ end {pmatrix}
1个2...ñ一世1个一世2...一世ñ




计算行列式的直接方法包括通过行或列的元素将其分解为任何行或列的元素的代数补码之和。反过来,矩阵元素的代数补码

一种一世Ĵ

--1个一世+Ĵ中号一世Ĵ

其中

中号一世Ĵ

-有一个次要元素(i,j),即通过删除第i行和第j列从原始行列式获得的行列式。



这种方法会生成一个递归过程,从而可以计算任何行列式。但是该算法的性能还有很多需要改进的地方-O(n!)。因此,除了符号计算(以及行列式不太高)外,都使用这种直接计算。



高斯方法的生产力更高。其实质是基于以下规定:



1.上三角矩阵的行列式\开始{pmatrix} {a} _ {1,1}&{a} _ {1,2}&…&{a} _ {1,n} \\ 0&{a} _ {2,2}&…&{a} _ {2,n} \\ 0&0&…&... \\ 0&0&…&{a} _ {n,n} \\\ end { pmatrix}等于其对角元素的乘积。根据行列式在第一行或第一列的元素上的分解,可以立即得出这一事实。 2.如果在矩阵中将一行的元素与另一行的元素相乘并乘以相同的数字,则行列式的值将保持不变。 3.如果矩阵中两行(或两列)互换,则行列式的值将变为相反的符号。等于其对角元素的乘积。根据行列在第一行或第一列的元素的分解,可以立即得出这一事实。
一种1个1个一种1个2...一种1个ñ0一种22...一种2ñ00............00...一种ññ












通过选择系数,我们可以将矩阵的第一行与所有其他行相加,并在除第一行以外的所有位置的第一列中获得零。要在第二行中获得零,您需要将第一行与第二行相乘

--一种21个/一种1个1个

要在第三行中获得零,请将第一行乘以

--一种31个/一种1个1个

等等 最终,矩阵将简化为所有元素的形式

一种ñ1个

对于n> 1将等于零。



如果在矩阵中

一种1个1个

结果等于0,那么您可以在第一列中找到一个非零元素(假设它在第k位)并交换第一行和第k行。通过这种转换,行列式可以简单地更改其符号,可以将其考虑在内。如果第一列中没有非零元素,则行列式等于零。



此外,以类似的方式操作,您可以在第二列中获取零,然后在第三列中获取,以此类推。重要的是,添加字符串时,较早获得的零不会改变。如果对于任何行都找不到分母的非零元素,则行列式等于零,然后可以停止该过程。高斯过程的正常完成会生成一个矩阵,其中位于主对角线以下的所有元素都等于零。如上所述,这种矩阵的行列式等于对角元素的乘积。



让我们继续编程。



我们使用浮点数据。我们将矩阵表示为字符串列表。首先,让我们定义两种类型:



type Row    = [Double]
type Matrix = [Row]


简单递归



毫不犹豫的是,我们将通过第一行(即零行)的元素来扩展行列式。我们需要一个用于构造次要元素的程序,该程序可通过删除第一行和第k列获得。



--  k-     
deln :: Matrix -> Int -> Matrix
deln matrix k = map (\ r -> (take (k) r)++(drop (k+1) r)) matrix


这是未成年人:



--  k-   
minor :: Matrix -> Int -> Double
minor matrix k = det $ deln (drop 1 matrix) k


请注意:未成年人是决定因素。我们正在调用尚未实现的det函数。要实现det,我们将必须形成第一行下一个元素的行列式与第一行下一个元素的乘积之和。为了避免麻烦的表达式,让我们创建一个单独的函数来形成和号:



sgn :: Int -> Double
sgn n = if n `rem` 2 == 0 then 1.0 else (-1.0)


现在我们可以计算行列式:



--   
det :: Matrix -> Double
det [[a,b],[c,d]] = a*d-b*c
det matrix = sum $ map (\c -> ((matrix !! 0)!!c)*(sgn c)*(minor matrix c))  [0..n]
             where n = length matrix - 1


该代码非常简单,不需要任何特殊注释。为了测试我们函数的性能,让我们编写主要函数:



main = print $ det [[1,2,3],[4,5,6],[7,8,(-9)]]


该行列式的值是54,可以验证。



高斯法



我们将需要一些实用程序功能(可在其他地方使用)。首先是矩阵中两行的互换:



--    
swap :: Matrix -> Int -> Int -> Matrix
swap matrix n1 n2 = map row [0..n]
                    where n=length matrix - 1
                          row k | k==n1 = matrix !! n2
                                | k==n2 = matrix !! n1
                                | otherwise = matrix !! k


从上面的代码可以看出,该函数逐行进行。在这种情况下,如果遇到数字为n1的字符串,则将强制插入字符串n2(反之亦然)。其余的行保持不变。



以下函数计算字符串r1加上字符串r2并逐个元素乘以数字f:



--   r1+f*r2
comb :: Row -> Row -> Double -> Row
comb r1 r2 f = zipWith (\ x y -> x+f*y) r1 r2


此处的所有内容都是非常透明的:在矩阵的行(即[Double]列表)上执行操作。但是下面的函数在矩阵上执行此转换(自然会得到一个新的矩阵):



--    r1  r2,   f
trans :: Matrix -> Int -> Int -> Double -> Matrix
trans matrix n1 n2 f = map row [0..n]
                       where n=length matrix - 1
                             row k | k==n1 = comb (matrix !! n1) (matrix !! n2) f
                                   | otherwise = matrix !! k


getNz函数在列表中查找第一个非零元素的编号。当下一个对角线元素等于零时,将需要它。



--     
getNz :: Row -> Int
getNz xs = if length tmp == 0 then (-1) else snd $ head tmp
           where tmp=dropWhile (\ (x,k) -> (abs x) <= 1.0e-10) $ zip xs [0..]


如果列表的所有元素均为零,则该函数将返回-1。



搜索功能检查矩阵是否适合下一个变换(它的下一个对角线元素必须为非零)。如果不是,则通过行的排列来变换矩阵。



--        
search :: Matrix -> Int -> Matrix
search matrix k | (abs ((matrix !! k) !! k)) > 1.0e-10 = matrix
                | nz < 0 = matrix  --      
                | otherwise = swap matrix k p 
                           where n   = length matrix
                                 lst = map (\ r -> r !! k) $ drop k matrix
                                 nz  = getNz lst
                                 p   = k + nz


如果找不到前导(非零)元素(矩阵退化),则该函数将其返回不变。mkzero函数在矩阵的下一列中形成零:



--     
mkzero :: Matrix -> Int -> Int -> Matrix
mkzero matrix k p | p>n-1 = matrix
                  | otherwise = mkzero (trans matrix p k (-f)) k (p+1)
                    where n = length matrix
                          f = ((matrix !! p) !! k)/((matrix !! k) !! k)


三角函数形成矩阵的上三角形状:



--     
triangle :: Matrix -> Int -> Matrix
triangle matrix k | k>=n = matrix
                  | (abs v) <= 1.0e-10 = [[0.0]] --  
                  | otherwise = triangle (mkzero tmp k k1) k1 
                    where n   = length matrix
                          tmp = search matrix k
                          v   = (tmp !! k) !! k --  
                          k1  = k+1


如果在下一阶段找不到前导元素,则该函数返回一阶零矩阵。现在,您可以构成将矩阵简化为上三角形式的游行函数:



--  
gauss :: Matrix -> Matrix
gauss matrix = triangle matrix 0 


要计算行列式,我们需要将对角线元素相乘。为此,我们创建一个单独的函数:



--   
proddiag :: Matrix -> Double
proddiag matrix = product $ map (\ (r,k) -> r !!k) $ zip matrix [0,1..]


好吧,“弓”是行列式的实际计算:



--  
det :: Matrix -> Double
det matrix = proddiag $ triangle matrix 0


让我们检查一下此函数的工作方式:



main = print $ det  [[1,2,3],[4,5,6],[7,8,-9]]

[1 of 1] Compiling Main             ( main.hs, main.o ) 
Linking a.out ...                                                                                 
54.0     


感谢那些读到最后的人!



可以在这里下载代码



All Articles