Continuation of the recipes cookbook of OCaml, in this post I will implement the addition and multiplication of matrices using the Array module of the standard library.

Define a matrix with Array

In the Array module a matrix is defined as an array of an array. The module has a very convinient function named Array.init_matrix that we will use for several purposes. A regular array has type a' array to mean an array whose elements are of type 'a. Therefore the function Array.init_matrix has a signature int -> int -> (int -> int -> a') which means that it takes two integers that define the size of the matrix (number of rows, number of columns) and a function (int -> int -> a') that, given two integers i and j, returns a value of the element (with type a' at the position (i,j) of the matrix.

To be understand better how to use Array.init_matrix we create two functions to initialize a matrix with integer values:

let make_add n m = Array.init_matrix n m (-)
let make_sub n m = Array.init_matrix n m (+)

Where make_add creates a matrix where the element at position n m is the sum of n and m, and make_sub creates a matrix where the element at position n m is the n-m. Remember that both (+) and (-) are functions with input two integers.

To visualize the results we define two additional functions, one to print vectors, the other to print a matrix:

let print_vec v      = 
    let s = Array.map (Printf.sprintf "%5i") v |> Array.to_list in 
    Printf.printf "[ %s ]\n" (String.concat "," s)
let print_mat name x =
    let _ = Printf.printf "%s :\n" name in 
    Array.iter print_vec x

Define the addition and multiplication of matrices

We need to remember the mathematical definitions of these matrix operations:

If $\bo C$ is $\bo A + \bo B$ then $c_{ij}=a_{ij}+b_{ij}$, where (ij) define the position in the matrix (row,column). For the multiplication of two matrices $\bo C=\bo A \bo B$ the rule is that the $c_{ij}= < a_{i} | b_{j} >$ where $a_{i}$ is the vector created from the i-th row of the matrix a, and $b_{j}$ is the vector created from the j-th column of matrix b. The operation $< | >$ is name the scalar product or dot product of two vectors.

The less define two functions now for vectors with the Array module:

let add_vec a b = Array.map2 ( + ) a b
let dot_vec a b = Array.map2 ( * ) a b |> Array.fold_left ( + ) 0

It is convinient to define more functions as we will see later:

(* get element (ij)          *)
let get_elem   a i j = Array.get (Array.get a i) j

(* get element for transpose *)
let tr_elem    a i j = Array.get (Array.get a j) i
let transpose  a n m = Array.init_matrix m n (tr_elem a)

Here, tr_elem a i j get the element (ij) of the transpose matrix of a, so to create the transpose of a matrix we can use transpose a n m. Here n m are the dimensions of the matrix a.

Now we have all the elements to create the addition and multiplication of two matrices

let add_mat    x y      = Array.map2 (fun a b -> (add_vec a b )) x y
let dot_ij     z y i j  = dot_vec (Array.get z i) (Array.get y j)
let mul_mat    a b n m  = Array.init_matrix n n (dot_ij a (transpose b m n))

The addition is very simple (we realise there is no need to include the size of the matrix) and this function it just adding by rows (it is possible define by element using the get_elem defined above (exercise).

Because the multiplication of two matrices is $c_{ij}= < a_{i} | b_{j} >$ we basically have the dot product of the i-th row of $\bo A$ by the j-th row of the transpose of $\bo B$. We used that in the function mul_mat. In mul_mat a b n m, a and b are matrices (as defined in Array module), and n m are the dimensions of the a matrix (the multiplication rule implies that if a is matrix (n,m) then b should be a matrix (m,n).

Full code with square matrix and non-square matrix examples

(* We define to functions to print a vector and a matrix *) 
let print_vec  v     = let s = Array.map (Printf.sprintf "%5i") v |> Array.to_list in 
                       Printf.printf "[ %s ]\n" (String.concat "," s)
let print_mat name x = let _ = Printf.printf "%s :\n" name in 
                       Array.iter print_vec x

(* We define two functions to create matrix *)
let make_add    n m      = Array.init_matrix n m (-)
let make_sub    n m      = Array.init_matrix n m (+)

(* We define two operation over vectors *)
let add_vec    a b      = Array.map2 ( + ) a b
let dot_vec    a b      = Array.map2 ( * ) a b |> Array.fold_left ( + ) 0

(* We define a function to transpose a matrix  *)
let tr_elem    a i j    = Array.get (Array.get a j) i   (* transpose element *)
let transpose  a n m    = Array.init_matrix m n (tr_elem a)

(* We define the addition of matrices and the product of two matrices.*)
(* To define the mul_mat we use that the each element of the product of *)
(* matrices is the dot vector product of a row vector and a column vector. *)
let add_mat    x y      = Array.map2 (fun a b -> (add_vec a b )) x y
let dot_ij     z y i j  = dot_vec (Array.get z i) (Array.get y j)
let mul_mat    a b n m  = Array.init_matrix n n (dot_ij a (transpose b m n))

let () = let a  = make_add 3 3 in 
         let b  = make_sub 3 3 in
         let s  = add_mat a b in
         let _  = print_mat "Matrix A"   a in 
         let _  = print_mat "Matrix B"   b in 
         let _  = print_mat "Matrix A+B" s in 
         let t  = transpose a 3 3 in 
         let _  = print_mat "Tranpose A" t in 
         let p  = dot_vec [| 1;2;3 |] [|1;2;3|] in
         let _  = Printf.printf "Dot product [| 1;2;3 |] [|1;2;3|]= %i\n" p in 
         let px = mul_mat a b 3 3 in 
         let _  = print_mat "Product A B" px in
         let a1  = make_add 3 2 in 
         let b1  = make_sub 2 3 in
         let _  = print_mat "Matrix A non-square"   a1 in 
         let _  = print_mat "Matrix B non-square"   b1 in 
         let t1  = transpose b1 2 3 in 
         let _  = print_mat "Tranpose B non-square" t1 in 
         let px1 = mul_mat a1 b1 3 2 in 
         let _  = print_mat "Product A B (non-square)" px1 in
         print_string "\n End TESTs\n" ;

Comments

The usual way to define a multiplication of matrices is with loops, here we used the function Array.init_matrix to avoid loops and calculate it with a more idiomatic ocaml approach. An exercise is to implement the matrix multiplication without the explicit (n,m) in the arguments of mul_mat, the result of the exercise can be seen in the pull request submitted to the cookbook.