Haskell行列の決定要因を計算します

決定要因を計算するためのコードを投稿することにしました。コードは機能していますが、名人のふりをしていません。Haskellでこの問題を解決するのは興味深いことでした。問題を解決するための2つのアプローチが考えられます:単純な再帰とガウス法です。



少し理論



ご存知のように、正方形行列n * nの決定要因は、nの合計です用語。それぞれは、各列から1つ、各行から1つだけのマトリックス要素を含む製品です。次の作品のサイン:

a11a22..。..。..。ann



置換のパリティによって決定され

ます\ begin {pmatrix} 1&2&…&n \\ {i} _ {1}&{i} _ {2}&…&{i} _ {n} \ end {pmatrix}
((12..。n12..。n




決定要因を計算する直接的な方法は、行または列の要素によってそれを代数的補数による任意の行または列の要素の積の合計に分解することです。次に、行列要素の代数的補数

aj

有る

((-1+jMj

ここで

Mj

-マイナーな要素(i、j)があります。 i番目の行とj番目の列を削除することによって元の決定子から取得された決定子。



このメソッドは、任意の決定要因を計算できるようにする再帰プロセスを生成します。しかし、このアルゴリズムのパフォーマンスには、まだ多くの要望があります-O(n!)。したがって、シンボリック計算を除いて(そして高すぎない決定要因を使用して)、そのような直接計算が使用されます。



ガウス法の方がはるかに生産的であることがわかります。その本質は、次の規定に基づいています。



1.上三角行列の決定要因\ begin {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.マトリックス内で2つの行(または2つの列)が交換されている場合、決定要因の値はその符号を反対に変更します。の対角要素の積に等しいです。この事実は、最初の行または最初の列の要素に関する決定要因の分解の直後に発生します。
((a11a12..。a1n0a22..。a2n00..。..。..。..。00..。ann












係数を選択することにより、マトリックスの最初の行を他のすべての行と一緒に追加し、最初の列を除くすべての位置の最初の列でゼロを取得できます。2行目でゼロを取得するには、最初の行に最初の行を追加して、を掛ける必要があります。

-a21/a11

3行目でゼロを取得するには、最初の行にを掛けたものを追加する必要があります

-a31/a11

最終的に、マトリックスはすべての要素が含まれる形式に縮小されます

an1

n> 1の場合、ゼロに等しくなります。



マトリックス内の場合、要素

a11

ゼロに等しいことが判明した場合、最初の列でゼロ以外の要素を見つけ(k番目の場所にあると仮定)、最初の行とk番目の行を入れ替えることができます。この変換により、決定要因はその符号を変更するだけであり、これを考慮に入れることができます。最初の列にゼロ以外の要素がない場合、決定要因はゼロに等しくなります。



さらに、同様の方法で動作すると、2番目の列でゼロを取得し、次に3番目の列でゼロを取得できます。文字列が追加されたときに、以前に取得されたゼロが変更されないことが重要です。どの行でも分母のゼロ以外の要素を見つけることができない場合、決定子はゼロに等しく、プロセスを停止できます。ガウスプロセスが正常に完了すると、主対角線の下にあるすべての要素がゼロに等しい行列が生成されます。上記のように、そのようなマトリックスの決定要因は、対角要素の積に等しくなります。



プログラミングに移りましょう。



浮動小数点データを処理します。行列を文字列のリストとして表します。まず、2つのタイプを定義しましょう。



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関数を記述しましょう。



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


この決定要因の値は54であり、検証できます。



ガウス法



いくつかのユーティリティ関数が必要になります(他の場所で使用できます)。1つ目は、マトリックス内の2つの行の交換です。



--    
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


上記のコードからわかるように、関数は1行ずつ実行されます。この場合、番号n1の行が検出されると、行n2が強制的に置換されます(またはその逆)。残りの行はそのまま残ります。



次の関数は、文字列r2に要素ごとに数値fを掛けた文字列r1を計算します。



--   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


次の段階で先頭の要素を見つけることができなかった場合、関数は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