Repository files navigation

SLAP is a linear algebra library inOCaml with type-basedstatic 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 vectorscauses runtime errors like memory corruption or wrong computation.

SLAP helps your debug by detecting inconsistency of dimensions

• atcompile time (instead of runtime) and
• athigher level (i.e., at a caller site rather than somewhere deep insideof a call stack).

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

This library is a wrapper ofLacaml, abinding of two widely used linear algebra libraries BLAS (Basic Linear AlgebraSubprograms) and LAPACK (Linear Algebra PACKage) for FORTRAN.This provides many useful and high-performance linear algebraic operations withtype-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 (requiringLacaml andcppo):

./configuremakemake install

• 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 bymake doc).
• This library interface was announced atML Family Workshop 2014 in Gothenburg,Sweden: A Simple and Practical Linear Algebra Library Interface with StaticSize Checking, by Akinori Abe and Eijiro Sumii (Tohoku University).PDF Abstract,PDF Slides,PDF Supplement.(The talk was accepted byOCaml Workshop 2014, but it waspresented at ML Workshop.)

The following code (examples/linsys/jacobi.ml)is an implementation ofJacobi method (to solvesystem of linear equations).

openSlap.IoopenSlap.CommonopenSlap.SizeopenSlap.Dletjacobiab=let d_inv=Vec.reci (Mat.diag a)in(* reciprocal diagonal elements*)let r=Mat.mapi (funijaij ->if i= jthen0.0else aij) ainlet y=Vec.create (Vec.dim b)in(* temporary memory*)letrecloopzx=    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*)ifVec.ssqr_diff x z<1e-10then zelse loop x z(* Check convergence*)inlet x0=Vec.make (Vec.dim b)1.0in(* the initial values of `x'*)let z=Vec.create (Vec.dim b)in(* temporary memory*)  loop z x0let()=let a= [%mat [5.0,1.0,0.0;1.0,3.0,1.0;0.0,1.0,4.0]]inlet b= [%vec [7.0;10.0;14.0]]inlet x= jacobi a binlet x= jacobi a binFormat.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 equationsa * x = b wherea isa n-by-n matrix, andx andb is a n-dimensional vectors. This code canbe compiled by

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

anda.out outputs:

a=510131014b=71014x=123a*x=71014

OK. Vectorx is computed. Try to modify any one ofthe dimensions ofa,b andx in the above code, e.g.,

...let()=let a= ...inlet 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 runtimeerror like an exception):

File"jacobi.ml", line31, characters19-20:Error:This expression hastype         (two, 'a) vec= (two,float, rprec, 'a)Slap_vec.t       but an expression was expected oftype         (three, 'b) vec= (three,float, rprec, 'b)Slap_vec.tType two= z s s isnot compatiblewithtype three= z s s sType z isnot compatiblewithtype z s

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

The following modules provide useful linear algebraic operations includingBLAS 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

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

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

Typesz 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 asingleton type on natural numbers; i.e., evaluation of a term (i.e.,expression) of type'n talways results in the natural numbercorresponding to'n. Thereforez s s s s t is the type of terms alwaysevaluated to four.

Creation of a four-dimensional vector:

#openSlap.D;;#let x=Vec.init four (funi -> float_of_int i);;valx : (zssss,'a)vec=R1R2R3R41234

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

Vectors of different dimensionsalways have different types:

#let y=Vec.init five (funi -> float_of_int i);;valy : (zsssss,'a)vec=R1R2R3R4R512345

The addition of four-dimensional vectorx and five-dimensional vectorycauses a type error (at compile time):

#Vec.add x y;;Error:This expression hastype  (z s s s s s, 'a) vec= (z s s s s s,float, rprec, 'a)Slap.Vec.tbut an expression was expected oftype  (z s s s s, 'b) vec= (z s s s s,float, rprec, 'b)Slap.Vec.tType z s isnot compatiblewithtype z

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

#let z=Vec.map (func -> c*.2.0) x;;valz : (zssss,'a)vec=R1R2R3R42468#Vec.addxz;;- : (zssss,'a)vec=R1R2R3R436912

Creation of a 3-by-5 matrix:

#let a=Mat.init three five (funij -> float_of_int (10* i+ j));;vala : (zsss,zsssss,'a)mat=C1C2C3C4C5R11112131415R22122232425R33132333435

('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 ofthe third type parameter is omitted.)

BLAS functiongemm multiplies two general matrices:

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

executesc := alpha * a * b + beta * c with matricesa,b andc, andscalar valuesalpha andbeta. The parameterstransa andtransb specifyno transpose (Slap.Common.normal), transpose (Slap.Common.trans) orconjugate transpose (Slap.Common.conjtr) of matricesa andb,respectively. (conjtr can be used only for complex operations inSlap.CandSlap.Z.) For example, iftransa=normal andtransb=trans, thengemm executesc := alpha * a * b^T + beta * c (whereb^T is the transposeofb). When you computea * a^T bygemm, a 3-by-3 matrix is returned sincea is a 3-by-5 matrix:

#openSlap.Common;;# gemm~transa:normal~transb:trans a a;;- : (z s s s, z s s s, 'a) mat=C1C2C3R185515052155R2150526553805R3215538055455

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

# gemm~transa:normal~transb:normal a a;;Error:This expression hastype  (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.tbut an expression was expected oftype  (z s s s s s, 'b, 'c) mat=  (z s s s s s, 'b,float, rprec, 'c)Slap.Mat.tType z isnot compatiblewithtype 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 avector 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 fileand call a SLAP function to convert it to a vector or a matrix.

Conversion of a list into a vector:

#moduleX= (valVec.of_list [1.;2.;3.] : Vec.CNTVEC);;moduleX :Slap.D.Vec.CNTVEC

The returned moduleX has the following signature:

module typeSlap.D.Vec.CNTVEC=sigtypen (* a type to represent the dimension of a vector *)valvalue : (n,'cnt)vec(* the instance of a vector*)end

The instance of a vector isX.value:

#X.value;;- : (X.n, 'cnt) vec=R1R2R3123

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

You can also convert a list into a matrix:

#moduleA= (valMat.of_list [[1.;2.;3.];                               [4.;5.;6.]] :Mat.CNTMAT);;#A.value;;- : (A.m,A.n, 'cnt) mat=C1C2C3R1123R2456

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

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

The third line should be ill-typed because the dimensions ofx andy areprobably different. (Even if"file1" and"file2" were the same path, theaddition should be ill-typed because the file might change between the twoloads.) Thus, the return type ofloadvec should be different every time it iscalled (regardless of the specific values of the argument). We call such areturn typegenerative because the function returns a value of a fresh typefor each call. The vector type with generative size information essentiallycorresponds to an existentially quantified sized type likeexists n. n vec.

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

Actually, typeX.n (returned byVec.of_list) is different for each call ofthe function, i.e., a generative phantom type:

#moduleY= (valVec.of_list [4.;5.] : Vec.CNTVEC);;#Vec.addX.valueY.value;;Error:This expression hastype  (Y.n, 'a) vec= (Y.n,float, rprec, 'a)Slap.Vec.tbut an expression was expected oftype  (X.n, 'b) vec= (X.n,float, rprec, 'b)Slap.Vec.tTypeY.n isnot compatiblewithtypeX.n

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

valVec.of_list_dyn : 'nSize.t ->floatlist -> ('n, 'cnt) vec

It also converts a list into a vector, but differs fromVec.of_list: You needto give the length of a list to the first parameter as a size. For example, ifyou consider that two listslst1 andlst2 (loaded from different files) havethe same length, you can add them as follows:

#let lst1= [1.;2.;3.;4.;5.];;(* loaded from a file*)vallst1 :floatlist= [1.;2.;3.;4.;5.]#letlst2= [6.;7.;8.;9.;10.];;(* loaded from another file*)vallst2 :floatlist= [6.;7.;8.;9.;10.]#moduleX= (valVec.of_list lst1 : Vec.CNTVEC);;moduleX :Slap.D.Vec.CNTVEC#let y=Vec.of_list_dyn (Vec.dimX.value) lst2;;valy : (X.n,'a)vec=R1R2R3R4R5678910#Vec.addX.valuey;;- : (X.n,'a)vec=R1R2R3R4R579111315

Vec.of_list_dyn raises anexception (at runtime) if the given size is notequal to the length, i.e., the lengths oflst1 andlst2are different in the above case. This dynamic check is unavoidable because theequality of sizes of two vectors loaded from different files cannot bestatically guaranteed. We gave functions containing dynamic checks the suffix_dyn.