CODEWITHSUNDEEP
functional-programmingocaml
alexgolec
alexgolec

Reputation: 28262

New to OCaml: How would I go about implementing Gaussian Elimination?

I'm new to OCaml, and I'd like to implement Gaussian Elimination as an exercise. I can easily do it with a stateful algorithm, meaning keep a matrix in memory and recursively operating on it by passing around a reference to it.

This statefulness, however, smacks of imperative programming. I know there are capabilities in OCaml to do this, but I'd like to ask if there is some clever functional way I haven't thought of first.

Upvotes: 7

Views: 829

Answers (3)

Toby Kelsey
Toby Kelsey

Reputation: 166

The answers so far are using/emulating mutable data-types, but what does a functional approach look like?

To see, let's decompose the problem into some functional components:

Gaussian elimination involves a sequence of row operations, so it is useful first to define a function taking 2 rows and scaling factors, and returning the resultant row operation result.

The row operations we want should eliminate a variable (column) from a particular row, so lets define a function which takes a pair of rows and a column index and uses the previously defined row operation to return the modified row with that column entry zero.

Then we define two functions, one to convert a matrix into triangular form, and another to back-substitute a triangular matrix to the diagonal form (using the previously defined functions) by eliminating each column in turn. We could iterate or recurse over the columns, and the matrix could be defined as a list, vector or array of lists, vectors or arrays. The input is not changed, but a modified matrix is returned, so we can finally do:

let out_matrix = to_diagonal (to_triangular in_matrix);

What makes it functional is not whether the data-types (array or list) are mutable, but how they they are used. This approach may not be particularly 'clever' or be the most efficient way to do Gaussian eliminations in OCaml, but using pure functions lets you express the algorithm cleanly.

Upvotes: 3

nlucaroni
nlucaroni

Reputation: 47944

You can use a Map to emulate a matrix. The key would be a pair of integers referencing the row and column. You'll want to use your own get x y function to ensure x < n and y < n though, instead of accessing the Map directly. (edit) You can use the compare function in Pervasives directly.

module OrderedPairs = struct
    type t = int * int
    let compare = Pervasives.compare
end                     
module Pairs = Map.Make (OrderedPairs)

let get_ n set x y =
    assert( x < n && y < n ); 
    Pairs.find (x,y) set

let set_ n set x y v = 
    assert( x < n && y < n ); 
    Pairs.add (x,y) set v

Actually, having a general set of functions (get x y and set x y at a minimum), without specifying the implementation, would be an even better option. The functions then can be passed to the function, or be implemented in a module through a functor (a better solution, but having a set of functions just doing what you need would be a first step since you're new to OCaml). In this way you can use a Map, Array, Hashtbl, or a set of functions to access a file on the hard-drive to implement the matrix if you wanted. This is the really important aspect of functional programming; that you trust the interface over exploiting the side-effects, and not worry about the underlying implementation --since it's presumed to be pure.

Upvotes: 4

Jeffrey Scofield
Jeffrey Scofield

Reputation: 66823

OCaml arrays are mutable, and it's hard to avoid treating them just like arrays in an imperative language.

Haskell has immutable arrays, but from my (limited) experience with Haskell, you end up switching to monadic, mutable arrays in most cases. Immutable arrays are probably amazing for certain specific purposes. I've always imagined you could write a beautiful implementation of dynamic programming in Haskell, where the dependencies among array entries are defined entirely by the expressions in them. The key is that you really only need to specify the contents of each array entry one time. I don't think Gaussian elimination follows this pattern, and so it seems it might not be a good fit for immutable arrays. It would be interesting to see how it works out, however.

Upvotes: 4

Related Questions