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 f

Slap

Install / Use

README