Sized Linear Algebra Package (SLAP)

SLAP is a linear algebra library in OCaml with type-based static size checking for matrix operations.

Many programming languages for numerical analysis (e.g., MatLab, GNU Octave, SciLab, etc.) and linear algebra libraries (e.g., BLAS, LAPACK, NumPy, etc.) do not statically (i.e., at compile time) guarantee consistency of dimensions of vectors and matrices. Dimensional inconsistency, e.g., addition of two- and three-dimensional vectors causes runtime errors like memory corruption or wrong computation.

SLAP helps your debug by detecting inconsistency of dimensions

• at compile time (instead of runtime) and
• at higher level (i.e., at a caller site rather than somewhere deep inside of a call stack).

For example, addition of vectors of different sizes causes a type error at compile time, and dynamic errors such as exceptions do not happen. For most high-level matrix operations, the consistency of sizes is verified statically. (Certain low-level operations, like accesses to elements by indices, need dynamic checks.)

This library is a wrapper of Lacaml, a binding of two widely used linear algebra libraries BLAS (Basic Linear Algebra Subprograms) and LAPACK (Linear Algebra PACKage) for FORTRAN. This provides many useful and high-performance linear algebraic operations with type-based static size checking, e.g., least squares problems, linear equations, Cholesky, QR-factorization, eigenvalue problems and singular value decomposition (SVD). Most of their interfaces are compatible with Lacaml functions.

Install

OPAM installation:

opam install slap

Manual build (requiring Lacaml and cppo):

./configure
make
make install

Documentation

• Web page: http://akabe.github.io/slap/
  • Tutorial: http://akabe.github.io/slap/usage.html
  • Online API documentation: http://akabe.github.io/slap/api/ (generated by make doc).
• This library interface was announced at ML Family Workshop 2014 in Gothenburg, Sweden: A Simple and Practical Linear Algebra Library Interface with Static Size Checking, by Akinori Abe and Eijiro Sumii (Tohoku University). PDF Abstract, PDF Slides, PDF Supplement. (The talk was accepted by OCaml Workshop 2014, but it was presented at ML Workshop.)

Demo

The following code (examples/linsys/jacobi.ml) is an implementation of Jacobi method (to solve system of linear equations).

open Slap.Io
open Slap.Common
open Slap.Size
open Slap.D

let jacobi a b =
  let d_inv = Vec.reci (Mat.diag a) in (* reciprocal diagonal elements *)
  let r = Mat.mapi (fun i j aij -> if i = j then 0.0 else aij) a in
  let y = Vec.create (Vec.dim b) in (* temporary memory *)
  let rec loop z x =
    ignore (copy ~y b); (* y := b *)
    ignore (gemv ~y ~trans:normal ~alpha:(-1.0) ~beta:1.0 r x); (* y := y-r*x *)
    ignore (Vec.mul ~z d_inv y); (* z := element-wise mult. of d_inv and y *)
    if Vec.ssqr_diff x z < 1e-10 then z else loop x z (* Check convergence *)
  in
  let x0 = Vec.make (Vec.dim b) 1.0 in (* the initial values of `x' *)
  let z = Vec.create (Vec.dim b) in (* temporary memory *)
  loop z x0

let () =
  let a = [%mat [5.0, 1.0, 0.0;
                 1.0, 3.0, 1.0;
                 0.0, 1.0, 4.0]] in
  let b = [%vec [7.0; 10.0; 14.0]] in
  let x = jacobi a b in
  let x = jacobi a b in
  Format.printf "a = @[%a@]@.b = @[%a@]@." pp_fmat a pp_rfvec b;
  Format.printf "x = @[%a@]@." pp_rfvec x;
  Format.printf "a*x = @[%a@]@." pp_rfvec (gemv ~trans:normal a x)

jacobi a b solves a system of linear equations a * x = b where a is a n-by-n matrix, and x and b is a n-dimensional vectors. This code can be compiled by

ocamlfind ocamlopt -package slap,slap.ppx -linkpkg -short-paths jacobi.ml

and a.out outputs:

a = 5 1 0
    1 3 1
    0 1 4
b = 7 10 14
x = 1 2 3
a*x = 7 10 14

OK. Vector x is computed. Try to modify any one of the dimensions of a, b and x in the above code, e.g.,

...

let () =
  let a = ... in
  let b = [%vec [7.0; 10.0]] in (* remove the last element `14.0' *)
  ...

and compile the changed code. Then OCaml reports a type error (not a runtime error like an exception):

File "jacobi.ml", line 31, characters 19-20:
Error: This expression has type
         (two, 'a) vec = (two, float, rprec, 'a) Slap_vec.t
       but an expression was expected of type
         (three, 'b) vec = (three, float, rprec, 'b) Slap_vec.t
       Type two = z s s is not compatible with type three = z s s s
       Type z is not compatible with type z s

By using SLAP, your mistake (i.e., a bug) is captured at compile time!

Usage

The following modules provide useful linear algebraic operations including BLAS and LAPACK functions.

• Slap.S: Single-precision (32-bit) real numbers
• Slap.D: Double-precision (64-bit) real numbers
• Slap.C: Single-precision (32-bit) complex numbers
• Slap.Z: Double-precision (64-bit) complex numbers

Simple computation

Dimensions (sizes)

Slap.Size provides operations on sizes (i.e., natural numbers as dimensions of vectors and matrices) of curious types. Look at the type of Slap.Size.four:

# open Slap.Size;;
# four;;
- : z s s s s t = 4

Types z and 'n s correspond to zero and 'n + 1, respectively. Thus, z s s s s represents 0+1+1+1+1 = 4. 'n t (= 'n Slap.Size.t) is a singleton type on natural numbers; i.e., evaluation of a term (i.e., expression) of type 'n t always results in the natural number corresponding to 'n. Therefore z s s s s t is the type of terms always evaluated to four.

Vectors

Creation of a four-dimensional vector:

# open Slap.D;;
# let x = Vec.init four (fun i -> float_of_int i);;
val x : (z s s s s, 'a) vec = R1 R2 R3 R4
                               1  2  3  4

Vec.init creates a vector initialized by the given function. ('n, 'a) vec is the type of 'n-dimensional vectors. So (z s s s s, 'a) vec is the type of four-dimensional vectors. (Description of the second type parameter is omitted.)

Vectors of different dimensions always have different types:

# let y = Vec.init five (fun i -> float_of_int i);;
val y : (z s s s s s, 'a) vec = R1 R2 R3 R4 R5
                                 1  2  3  4  5

The addition of four-dimensional vector x and five-dimensional vector y causes a type error (at compile time):

# Vec.add x y;;
Error: This expression has type
  (z s s s s s, 'a) vec = (z s s s s s, float, rprec, 'a) Slap.Vec.t
but an expression was expected of type
  (z s s s s, 'b) vec = (z s s s s, float, rprec, 'b) Slap.Vec.t
Type z s is not compatible with type z

Of course, addition of vectors of the same dimension succeeds:

# let z = Vec.map (fun c -> c *. 2.0) x;;
val z : (z s s s s, 'a) vec = R1 R2 R3 R4
                               2  4  6  8
# Vec.add x z;;
- : (z s s s s, 'a) vec = R1 R2 R3 R4
                           3  6  9 12

Matrices

Creation of a 3-by-5 matrix:

# let a = Mat.init three five (fun i j -> float_of_int (10 * i + j));;
val a : (z s s s, z s s s s s, 'a) mat =
     C1 C2 C3 C4 C5
  R1 11 12 13 14 15
  R2 21 22 23 24 25
  R3 31 32 33 34 35

('m, 'n, 'a) mat is the type of 'm-by-'n matrices. Thus (z s s s, z s s s s s, 'a) mat is the type of 3-by-5 matrices. (Description of the third type parameter is omitted.)

BLAS function gemm multiplies two general matrices:

gemm ?beta ?c ~transa ?alpha a ~transb b

executes c := alpha * a * b + beta * c with matrices a, b and c, and scalar values alpha and beta. The parameters transa and transb specify no transpose (Slap.Common.normal), transpose (Slap.Common.trans) or conjugate transpose (Slap.Common.conjtr) of matrices a and b, respectively. (conjtr can be used only for complex operations in Slap.C and Slap.Z.) For example, if transa=normal and transb=trans, then gemm executes c := alpha * a * b^T + beta * c (where b^T is the transpose of b). When you compute a * a^T by gemm, a 3-by-3 matrix is returned since a is a 3-by-5 matrix:

# open Slap.Common;;
# gemm ~transa:normal ~transb:trans a a;;
- : (z s s s, z s s s, 'a) mat =
     C1   C2   C3
R1  855 1505 2155
R2 1505 2655 3805
R3 2155 3805 5455

a * a causes a type error since the number of columns of the first matrix is not equal to the number of rows of the second matrix:

# gemm ~transa:normal ~transb:normal a a;;
Error: This expression has type
  (z s s s, z s s s s s, 'a) mat =
  (z s s s, z s s s s s, float, rprec, 'a) Slap.Mat.t
but an expression was expected of type
  (z s s s s s, 'b, 'c) mat =
  (z s s s s s, 'b, float, rprec, 'c) Slap.Mat.t
Type z is not compatible with type z s s

Sizes decided at runtime

SLAP can safely treat sizes that are unknown until runtime (e.g., the dimension of a vector loaded from a file)! Unfortunately, SLAP does not provide functions to load a vector or a matrix from a file. (Maybe such operations will be implemented.) You need to write a function to load a list or an array from a file and call a SLAP function to convert it to a vector or a matrix.

Conversion of a list into a vector:

# module X = (val Vec.of_list [1.; 2.; 3.] : Vec.CNTVEC);;
module X : Slap.D.Vec.CNTVEC

The returned module X has the following signature:

module type Slap.D.Vec.CNTVEC = sig
  type n (* a type to represent the dimension of a vector *)
  val value : (n, 'cnt) vec (* the instance of a vector *)
end

The instance of a vector is X.value:

# X.value;;
- : (X.n, 'cnt) vec = R1 R2 R3
                       1  2  3

It can be treated as stated above. It's very easy!

You can also convert a list into a matrix:

# module A = (val Mat.of_list [[1.; 2.; 3.];
                               [4.; 5.; 6.]] : Mat.CNTMAT);;
# A.value;;
- : (A.m, A.n, 'cnt) mat =    C1 C2 C3
                           R1  1  2  3
                           R2  4  5  6

Idea of generative phantom type

In this section, we explain our basic idea of static size checking. For example, let's think about the function loadvec : string -> (?, _) vec. It returns a vector of some dimension, loaded from the given path. The dimension is decided at runtime, but we need to type it at compile time. How do we represent the return type ?? Consider the following code for example:

let (x : (?1, _) vec) = loadvec "file1" in
let (y : (?2, _) vec) = loadvec "file2" in
Vec.add x y

The third line should be ill-typed because the dimensions of x and y are probably different. (Even if "file1" and "file2" were the same path, the addition should be ill-typed because the file might change between the two loads.) Thus, the return type of loadvec should be different every time it is called (regardless of the specific values of the argument). We call such a return type generative because the function returns a value of a fresh type for each call. The vector type with generative size information essentially corresponds to an existentially quantified sized type like exists n. n vec.

Type parameters 'm, 'n and 'a of types 'n Size.t, ('n, 'a) vec and ('m, 'n, 'a) mat are phantom, meaning that they do not appear on the right hand side of the type definition. A phantom type parameter is often instantiated with a type that has no value (i.e., no constructor) which we call a phantom type. Then we call the type ? a generative phantom type.

Actually, type X.n (returned by Vec.of_list) is different for each call of the function, i.e., a generative phantom type:

# module Y = (val Vec.of_list [4.; 5.] : Vec.CNTVEC);;
# Vec.add X.value Y.value;;
Error: This expression has type
  (Y.n, 'a) vec = (Y.n, float, rprec, 'a) Slap.Vec.t
but an expression was expected of type
  (X.n, 'b) vec = (X.n, float, rprec, 'b) Slap.Vec.t
Type Y.n is not compatible with type X.n

Addition of vectors loaded from different files

When you want to add vectors loaded from different files, you can use Vec.of_list_dyn:

val Vec.of_list_dyn : 'n Size.t -> float list -> ('n, 'cnt) vec

It also converts a list into a vector, but differs from Vec.of_list: You need to give the length of a list to the first parameter as a size. For example, if you consider that two lists lst1 and lst2 (loaded from different files) have the same length, you can add them as follows:

# let lst1 = [1.; 2.; 3.; 4.; 5.];; (* loaded from a file *)
val lst1 : float list = [1.; 2.; 3.; 4.; 5.]
# let lst2 = [6.; 7.; 8.; 9.; 10.];; (* loaded from another file *)
val lst2 : float list = [6.; 7.; 8.; 9.; 10.]
# module X = (val Vec.of_list lst1 : Vec.CNTVEC);;
module X : Slap.D.Vec.CNTVEC
# let y = Vec.of_list_dyn (Vec.dim X.value) lst2;;
val y : (X.n, 'a) vec = R1 R2 R3 R4 R5
                         6  7  8  9 10
# Vec.add X.value y;;
- : (X.n, 'a) vec = R1 R2 R3 R4 R5
                     7  9 11 13 15

Vec.of_list_dyn raises an exception (at runtime) if the given size is not equal to the length, i.e., the lengths of lst1 and lst2 are different in the above case. This dynamic check is unavoidable because the equality of sizes of two vectors loaded from different files cannot be statically guaranteed. We gave functions containing dynamic checks the suffix _dyn.