Вычисляем определитель матрицы на Хаскелле

Решил выложить код вычисления определителей. Код рабочий, хотя и не претендует на виртуозность. Просто было интересно решить эту задачу именно на Хаскелле. Рассмотрены два подхода к решению задачи: простая рекурсия и метод Гаусса.

Немного теории

Как известно, определитель квадратной матрицы n*n — это сумма n! слагаемых, каждое из которых есть произведение, содержащее ровно по одному элементу матрицы из каждого столбца и ровно по одному из каждой строки. Знак очередного произведения:

определяется чётностью подстановки:

Прямой метод вычисления определителя состоит в разложении его по элементам строки или столбца в сумму произведений элементов какой-либо строки или столбца на их алгебраические дополнения. В свою очередь, алгебраическое дополнение элемента матрицы

есть

при этом

— есть минор элемента (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. Если в матрице поменять местами две строки (или два столбца), то значение определителя изменит знак на противоположный.

Мы можем, подбирая коэффициенты, складывать первую строку матрицы со всеми остальными и получать в первом столбце нули во всех позициях, кроме первой. Для получения нуля во второй строке, нужно прибавить ко второй строке первую, умноженную на

Для получения нуля в третьей строке, нужно к третьей строке прибавить первую строку, умноженную на

и т.д. В конечном итоге, матрица приведется к виду, в котором все элементы

при n>1 будут равны нулю.

Если же в матрице элемент

оказался равным нулю, то можно найти в первом столбце ненулевой элемент (предположим, он оказался на 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:

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 проверяет, подходит ли матрица для очередного преобразования (у нее должен быть ненулевым очередной диагональный элемент). Если это не так, матрица преобразовывается перестановкой строк.

-- Поиск ведущего элемента и перестановка строк при необходимости
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 формирует верхнюю треугольную форму матрицы:

-- Получение верхней треугольной формы матрицы
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     

Спасибо тем, кто дочитал до конца!

Код можно скачать здесь