CODEWITHSUNDEEP
mathhaskelllinear-algebra
Alexei
Alexei

Reputation: 347

Implement Gauss-Jordan elimination in Haskell

We want to program the gauss-elimination to calculate a basis (linear algebra) as exercise for ourselves. It is not homework.

I thought first of [[Int]] as structure for our matrix. I thought then that we can sort the lists lexicographically. But then we must calculate with the matrix. And there is the problem. Can someone give us some hints.

Upvotes: 2

Views: 2657

Answers (2)

Jan
Jan

Reputation: 11726

Consider using matrices from the hmatrix package. Among its modules you can find both a fast implementation of a matrix and a lot of linear algebra algorithms. Browsing their sources might help you with your doubts.

Here's a simple example of adding one row to another by splitting the matrix into rows.

import Numeric.Container
import Data.Packed.Matrix

addRow :: Container Vector t => Int -> Int -> Matrix t -> Matrix t
addRow from to m = let rows = toRows m in
  fromRows $ take to rows ++
             [(rows !! from) `add` (rows !! to)] ++
             drop (to + 1) rows

Another example, this time by using matrix multiplication.

addRow :: (Product e, Container Vector e) =>
          Int -> Int -> Matrix e -> Matrix e
addRow from to m = m `add` (e <> m)
  where
    nrows = rows m
    e = buildMatrix nrows nrows
        (\(r,c) -> if (r,c) /= (to,from) then 0 else 1)

Cf. Container, Vector, Product.

Upvotes: 4

redxaxder
redxaxder

Reputation: 600

It will be easier if you use [[Rational]] instead of [[Int]] since you get nice division.

You probably want to start by implementing the elementary row operations.

swap :: Int -> Int -> [[Rational]] -> [[Rational]
swap r1 r2 m = --a matrix with r1 and r2 swapped

scale :: Int -> Rational -> [[Rational]] -> [[Rational]]
scale r c m = --a matrix with row r multiplied by c

addrow :: Int -> Int -> Rational -> [[Rational]] -> [[Rational]]
addrow r1 r2 c m = --a matrix with (c * r1) added to r2

In order to actually do guassian elimination, you need a way to decide what multiple of one row to add to another to get a zero. So given two rows..

5 4 3 2 1
7 6 5 4 3

We want to add c times row 1 to row 2 so that the 7 becomes a zero. So 7 + c * 5 = 0 and c = -7/5. So in order to solve for c all we need are the first elements of each row. Here's a function that finds c:

whatc :: Rational -> Rational -> Rational
whatc _ 0 = 0
whatc a b = - a / b

Also, as others have said, using lists to represent your matrix will give you worse performance. But if you're just trying to understand the algorithm, lists should be fine.

Upvotes: 1

Related Questions