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
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.
OPAM installation:
opam install slap
Manual build (requiring Lacaml and cppo):
./configure
make
make install
make doc
).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 = @[%[email protected]]@.b = @[%[email protected]]@." pp_fmat a pp_rfvec b;
Format.printf "x = @[%[email protected]]@." pp_rfvec x;
Format.printf "a*x = @[%[email protected]]@." 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!
The following modules provide useful linear algebraic operations including BLAS and LAPACK functions.
Slap.S
: Single-precision (32-bit) real numbersSlap.D
: Double-precision (64-bit) real numbersSlap.C
: Single-precision (32-bit) complex numbersSlap.Z
: Double-precision (64-bit) complex numbersSlap.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.
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
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
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
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
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
.