A bit of theory
As you know, the determinant of a square matrix n * n is the sum of n! terms, each of which is a product containing exactly one matrix element from each column and exactly one from each row. The sign of the next piece:
determined by the parity of the substitution:
\ begin {pmatrix} 1 & 2 &… & n \\ {i} _ {1} & {i} _ {2} &… & {i} _ {n} \ end {pmatrix}
A direct method for calculating the determinant consists in decomposing it by the elements of a row or column into the sum of the products of elements of a row or column by their algebraic complements. In turn, the algebraic complement of the matrix element
there is
wherein
- there is a minor element (i, j), i.e. determinant obtained from the original determinant by deleting the i-th row and j-th column.
Such a method generates a recursive process that allows any determinant to be computed. But the performance of this algorithm leaves much to be desired - O (n!). Therefore, such a direct computation is used, except for symbolic calculations (and with determinants of not too high order).
The Gauss method turns out to be much more productive. Its essence is based on the following provisions:
1. Determinant of the upper triangular matrix \ 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} is equal to the product of its diagonal elements. This fact immediately follows from the decomposition of the determinant into the elements of the first row or the first column. 2. If in the matrix to the elements of one row add the elements of another row, multiplied by the same number, then the value of the determinant will not change. 3. If two rows (or two columns) are interchanged in the matrix, then the value of the determinant will change its sign to the opposite. equal to the product of its diagonal elements. This fact immediately follows from the decomposition of the determinant into the elements of the first row or the first column.
By choosing the coefficients, we can add the first row of the matrix with all the others and get zeros in the first column in all positions except the first. To get zero in the second line, you need to add to the second line the first, multiplied by
To get zero in the third line, you need to add the first line multiplied by
etc. Ultimately, the matrix will be reduced to a form in which all elements
for n> 1 will be equal to zero.
If in the matrix the element
turns out to be equal to zero, then you can find a nonzero element in the first column (suppose it is in the kth place) and swap the first and kth rows. With this transformation, the determinant simply changes its sign, which can be taken into account. If there are no nonzero elements in the first column, then the determinant is equal to zero.
Further, acting in a similar way, you can get zeros in the second column, then in the third, etc. It is important that the zeros obtained earlier will not change when the strings are added. If for any line it is not possible to find a nonzero element for the denominator, then the determinant is equal to zero and the process can be stopped. The normal completion of the Gaussian process generates a matrix in which all the elements located below the main diagonal are equal to zero. As mentioned above, the determinant of such a matrix is equal to the product of the diagonal elements.
Let's move on to programming.
We work with floating point data. We represent matrices as lists of strings. First, let's define two types:
type Row = [Double]
type Matrix = [Row]
Simple recursion
No more hesitation, we will expand the determinant by the elements of the first (i.e. zero) row. We need a program for constructing the minor, obtained by deleting the first row and the kth column.
-- k-
deln :: Matrix -> Int -> Matrix
deln matrix k = map (\ r -> (take (k) r)++(drop (k+1) r)) matrix
And here is the minor:
-- k-
minor :: Matrix -> Int -> Double
minor matrix k = det $ deln (drop 1 matrix) k
Please note: minor is a determinant. We are calling the det function, which we have not yet implemented. To implement det, we will have to form the alternating sum of the products of the next element of the first row by the determinant of the next minor. To avoid cumbersome expressions, let's create a separate function to form the sum sign:
sgn :: Int -> Double
sgn n = if n `rem` 2 == 0 then 1.0 else (-1.0)
Now we can calculate the determinant:
--
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
The code is very simple and doesn't require any special comments. To test the performance of our functions, let's write the main function:
main = print $ det [[1,2,3],[4,5,6],[7,8,(-9)]]The value of this determinant is 54, which can be verified.
Gauss method
We'll need a few utility functions (which can be used elsewhere). The first is the interchange of two rows in the matrix:
--
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
As you can see from the above code, the function goes through line by line. In this case, if a line with number n1 is encountered, the line n2 is forcibly inserted (and vice versa). The rest of the lines remain in place.
The following function calculates the string r1 added with the string r2 multiplied element by element by the number f:
-- r1+f*r2
comb :: Row -> Row -> Double -> Row
comb r1 r2 f = zipWith (\ x y -> x+f*y) r1 r2
Everything here is extremely transparent: actions are performed on the rows of the matrix (ie, on [Double] lists). But the following function performs this transformation on the matrix (and, naturally, gets a new matrix):
-- 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
The getNz function looks for the number of the first nonzero element in the list. It is needed when the next diagonal element is equal to zero.
--
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..]
If all elements of the list are zero, the function will return -1.
The search function checks if the matrix is suitable for the next transformation (it must have a nonzero next diagonal element). If it is not, the matrix is transformed by permutation of rows.
--
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
If it is impossible to find the leading (nonzero) element (the matrix is degenerate), then the function will return it unchanged. The mkzero function forms zeros in the next column of the matrix:
--
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)
The triangle function forms the upper triangular shape of the matrix:
--
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
If at the next stage it was not possible to find the leading element, the function returns a zero matrix of the 1st order. Now you can compose the parade function of reducing the matrix to the upper triangular form:
--
gauss :: Matrix -> Matrix
gauss matrix = triangle matrix 0
To calculate the determinant, we need to multiply the diagonal elements. To do this, let's create a separate function:
--
proddiag :: Matrix -> Double
proddiag matrix = product $ map (\ (r,k) -> r !!k) $ zip matrix [0,1..]
Well, and the "bow" is the actual calculation of the determinant:
--
det :: Matrix -> Double
det matrix = proddiag $ triangle matrix 0
Let's check how this function works:
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
Thanks to those who read to the end!
The code can be downloaded here