Calculer le déterminant Haskell d'une matrice

J'ai décidé d'afficher le code de calcul des déterminants. Le code fonctionne, même s'il ne prétend pas être virtuose. C'était juste intéressant de résoudre ce problème sur Haskell. Deux approches pour résoudre le problème sont envisagées: la récursion simple et la méthode gaussienne.

Un peu de théorie

Comme vous le savez, le déterminant d'une matrice carrée n * n est la somme de n! termes, chacun étant un produit contenant exactement un élément de matrice de chaque colonne et exactement un de chaque ligne. Le signe de la pièce suivante:

$$

$${a}_{1,i1}*{a}_{2,i2}*...{a}_{n,in}$$

est déterminé par la parité de la substitution:

$$\begin{pmatrix}1 & 2 &… & n \\ {i}_{1} & {i}_{2} &… & {i}_{n} \end{pmatrix}$$



La méthode directe de calcul du déterminant consiste à le décomposer par les éléments d'une ligne ou d'une colonne en somme des produits des éléments de n'importe quelle ligne ou colonne par leurs compléments algébriques. À son tour, le complément algébrique de l'élément de matrice

$$

$${a}_{i,j}$$

il y a

$$

$${(-1)}^{i+j}*{M}_{i,j}$$

où

$$

$${M}_{i,j}$$

- il y a un élément mineur (i, j), i.e. déterminant obtenu à partir du déterminant original en supprimant la i-ème ligne et la j-ème colonne.



Cette méthode génère un processus récursif qui permet de calculer n'importe quel déterminant. Mais les performances de cet algorithme laissent beaucoup à désirer - O (n!). Par conséquent, un tel calcul direct est utilisé, sauf pour les calculs symboliques (et avec des déterminants d'ordre pas trop élevé).



La méthode Gauss s'avère beaucoup plus productive. Son essence est basée sur les dispositions suivantes:



1. Déterminant de la matrice triangulaire supérieure $$\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}$$ est égal au produit de ses éléments diagonaux. Ce fait découle immédiatement de la décomposition du déterminant en termes d'éléments de la première ligne ou de la première colonne.



2. Si dans la matrice aux éléments d'une ligne, ajoutez les éléments d'une autre ligne, multipliés par le même nombre, alors la valeur du déterminant ne changera pas.



3. Si deux lignes (ou deux colonnes) sont interchangées dans la matrice, la valeur du déterminant changera son signe en l'inverse.



On peut, en choisissant les coefficients, ajouter la première ligne de la matrice avec toutes les autres et obtenir des zéros dans la première colonne dans toutes les positions sauf la première. Pour obtenir zéro dans la deuxième ligne, vous devez ajouter à la deuxième ligne la première, multipliée par

$$

$$-{a}_{2,1}/{a}_{1,1}$$

Pour obtenir zéro sur la troisième ligne, ajoutez la première ligne multipliée par

$$

$$-{a}_{3,1}/{a}_{1,1}$$

etc. En fin de compte, la matrice sera réduite à une forme dans laquelle tous les éléments

$$

$${a}_{n,1}$$

pour n> 1 sera égal à zéro.



Si dans la matrice l'élément

$$

$${a}_{1,1}$$

s'avère être égal à zéro, alors vous pouvez trouver un élément différent de zéro dans la première colonne (supposons qu'il soit à la kème place) et permuter les première et kème lignes. Avec cette transformation, le déterminant change simplement son signe, qui peut être pris en compte. S'il n'y a aucun élément différent de zéro dans la première colonne, le déterminant est égal à zéro.



De plus, agissant de la même manière, vous pouvez obtenir des zéros dans la deuxième colonne, puis dans la troisième, etc. Il est important que les zéros obtenus précédemment ne changent pas lorsque les chaînes sont ajoutées. Si pour une ligne quelconque, il n'est pas possible de trouver un élément différent de zéro pour le dénominateur, alors le déterminant est égal à zéro et le processus peut être arrêté. L'achèvement normal du processus gaussien génère une matrice dans laquelle tous les éléments situés sous la diagonale principale sont égaux à zéro. Comme mentionné ci-dessus, le déterminant d'une telle matrice est égal au produit des éléments diagonaux.

Passons à la programmation.

Nous travaillons avec des données en virgule flottante. Nous représentons les matrices comme des listes de chaînes. Tout d'abord, définissons deux types:

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

Récursivité simple

Rien d'hésitant, nous allons étendre le déterminant par les éléments de la première ligne (c'est-à-dire zéro). Nous avons besoin d'un programme pour construire le mineur, obtenu en supprimant la première ligne et la kème colonne.

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

Et voici le mineur:

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

Remarque: le mineur est un déterminant. Nous appelons la fonction det, que nous n'avons pas encore implémentée. Pour implémenter det, nous devrons former la somme alternée des produits de l'élément suivant de la première ligne par le déterminant du mineur suivant. Pour éviter les expressions encombrantes, créons une fonction distincte pour former le signe somme:

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

Maintenant, nous pouvons calculer le déterminant:

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

Le code est très simple et ne nécessite aucun commentaire particulier. Pour tester les performances de nos fonctions, écrivons la fonction principale:

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

La valeur de ce déterminant est 54, ce qui peut être vérifié.

Méthode Gauss

Nous aurons besoin de quelques fonctions utilitaires (qui peuvent être utilisées ailleurs). Le premier est l'échange de deux lignes dans la matrice:

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

Comme vous pouvez le voir dans le code ci-dessus, la fonction passe ligne par ligne. Dans ce cas, si une ligne de numéro n1 est rencontrée, la ligne n2 est insérée de force (et vice versa). Le reste des lignes reste en place.



La fonction suivante calcule la chaîne r1 ajoutée avec la chaîne r2 multipliée élément par élément par le nombre f:

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

Tout ici est extrêmement transparent: les actions sont effectuées sur les lignes de la matrice (c'est-à-dire sur les listes [Double]). Mais la fonction suivante effectue cette transformation sur la matrice (et, naturellement, obtient une nouvelle matrice):

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

La fonction getNz recherche le numéro du premier élément différent de zéro dans la liste. Il est nécessaire lorsque l'élément diagonal suivant est égal à zéro.

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

Si tous les éléments de la liste sont à zéro, la fonction renverra -1.



La fonction de recherche vérifie si la matrice est adaptée à la transformation suivante (elle doit avoir un élément diagonal suivant différent de zéro). Si ce n'est pas le cas, la matrice est transformée par permutation de lignes.

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

S'il est impossible de trouver l'élément principal (différent de zéro) (la matrice est dégénérée), alors la fonction le renverra inchangé. La fonction mkzero forme des zéros dans la colonne suivante de la matrice:

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

La fonction triangle forme la forme triangulaire supérieure de la matrice:

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

Si à l'étape suivante il n'a pas été possible de trouver l'élément de tête, la fonction renvoie une matrice nulle du 1er ordre. Maintenant, nous pouvons composer la fonction cérémonielle de réduction de la matrice à la forme triangulaire supérieure:

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

Pour calculer le déterminant, nous devons multiplier les éléments diagonaux. Pour ce faire, créons une fonction distincte:

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

Eh bien, et «l'arc» est le calcul réel du déterminant:

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

Voyons comment fonctionne cette fonction:

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     

Merci à ceux qui ont lu jusqu'au bout!



Le code peut être téléchargé ici