Last days I was checking the cookbook in the new webpage of OCAML. This cookbook webpage is open to contributions, for that it is only need to fork the repository ocaml.org and follow the steps to submmit the new recipe to the cookbook as a regular pull-request.

Inverse matrix

One of the recipes that needs help is the calculation of the inverse of a matrix. This recipe is quite easy using the Owl library as it has an specific function implemented for this. However, it is interesting to include a method that ensure that the input matrix is invertible (like assess that the determinant is different than zero), and perform a checking of the result.

Remember that if $\bo M$ is a square matrix (n-by-n) then it is invertible if, and only if, $det(\bo M)\=0$. If $\bo M$ is invertible and it's inverse matrix is named $\bo X $ then the product $\bo M \bo X=\bo I$, where $\bo I $ represents the identity matrix (n-by-n0. The identity matrix has the following key property regarding the matrix product: $\bo C \bo I=\bo I \bo C=\bo C$, for any $\bo C$ square matrix (n-by-n). Also $\bo I$ it is a diagonal matrix with ones in the diagonal (and all other elements zero).

Setting up a project with dune

> dune init project invertmatrix

and we edit the file bin/dune to be:

(executable
 (public_name invm)
 (name main)
 (libraries owl))

so we will edit the code at bin/main.ml, and with the bin/dune above defined we can use dune exec invm just to run the exectuable.

Define an invertible matrix

We can define an invertible random matrix using the Owl function Owl.Mat.uniform n m that needs two arguments for the matrix dimensions (number of rows n and number of columns m). In our case it will be an square matrix, therefore n=m. This random generated matrix would be very likely invertible but we can ensure this by calculating the determinant of the matrix. If it is 0. then we need to recompute, otherwise it is a good input matrix candidate.

let rec create_invertible m =
        let input = Owl.Mat.uniform m m in
        match Owl.Linalg.D.det input with
        | 0. -> create_invertible m 
        | _  -> input

Calculate and check the matrix inverse

Now given the input matrix we can use the function Owl.Mat.inv to invert the matrix and we can check that the input matrix multipled by the inverse matrix is the identity matrix (knowing that the non-diagonal elements will be zero, but within the precision given by the float numbers).

(* 1. Define an input matrix that is invertible      *)
let input  = create_invertible 5

(* 2. Calculates the inverse using `Owl.Mat.inv`     *)
let inverse = Owl.Mat.inv input

(* 3. Calculates the matrix product: input * inverse *)
(*    it uses the operator *@                        *)
let id_test = Owl.Mat.(input *@ inverse)

(* 4. Result should be an identity matrix.           *)
(*    Non-diagonal elements should be reasonably     *) 
(*    close to zero (round-off effect). We can test  *)
(*    this using the =~ also available in Owl.Mat    *)
let test_result input inverse = 
        Owl.Mat.(input *@ inverse =~ eye 5);;

Complete code for cookbook

let rec create_invertible m =
        let input = Owl.Mat.uniform m m in
        match Owl.Linalg.D.det input with
        | 0. -> create_invertible m 
        | _  -> input

let input  = create_invertible 5

let inverse = Owl.Mat.inv input

let id_test = Owl.Mat.(input *@ inverse)

let test_result input inverse = 
        Owl.Mat.(input *@ inverse =~ eye 5);;

let () = 
    print_string "** Input matrix (random numbers)" ;
    Owl.Mat.(input  |> print) ; flush_all() ;
    print_string "\n** Inverse of input matrix" ;
    Owl.Mat.(inverse |> print) ; flush_all() ;
    print_string "\n** Product input and inverse =~ identity" ;
    Owl.Mat.(id_test |> print) ; flush_all() ;
    print_string "\n** Is the inverse correct?  " ;
    test_result input inverse |> string_of_bool |> print_string ; 
    print_endline ""

Suggesting the recipe

To suggest this recipe it is only need to fork the official repository, and create in the corresponding directory the file with the proposal for the recipe. I created a file 01_owl.ml, see the repository of ocaml.org. Note that it is also possible to propose additional recipes. The final step is to create a full request.

Further exercises

There are other properties related to inverse of a matrix in particular about the named invertible matrix theorem. Further exercises could test these other methods to assess if a matrix is invertible. This would give us better knowledge of the Owl library.