Sci.Alg: Vector and Matrix Algebra

A vector and matrix algebra library with a focus on:

• simplicity: a simple API built around a single kind of objects: 2-D arrays
• generality: dense homogeneous arrays of Lua numbers (doubles) or arbitrary value-type FFI ctype
• performance: carefully designed for efficiency with use of BLAS and LAPACK kernels
• safety: out of bounds access and dimensional compatibility are checked, algorithms are unaffected by aliasing
• consistency: 1-based indexing as for Lua tables

This library implements 2-dimensional arrays which can be used at the same time as vectors and as matrices:

local alg = require 'sci.alg'

local x = alg.mat(2, 3) -- A 2-rows x 3-columns array, i.e. a matrix.

-- Can be used as a matrix:
for r=1,x:nrow() do 
  for c=1,x:ncol() do
    x[{r, c}] = r*10 + c
  end
end

print(x)
--> +11.00000,+12.00000,+13.00000
--> +21.00000,+22.00000,+23.00000

-- Can be used as a vector too, access is left-to-right, top-to-bottom:
for i=1,#x do
  x[i] = i
end

print(x)
--> +1.000000,+2.000000,+3.000000
--> +4.000000,+5.000000,+6.000000

local y = alg.vec(6) -- A 6-rows x 1-column array, i.e. a (column) vector.

print(y)
--> +0.000000
--> +0.000000
--> +0.000000
--> +0.000000
--> +0.000000
--> +0.000000

In practice, most code will interface with each single array as if it were only a vector or only a matrix. No performance penalty is involved in this dual representation.

The length of an array is defined as its number of rows times its number of columns. Arrays with zero rows, columns or length are allowed.

API

alg = require 'sci.alg'

Returns the algebra module, no global is set.

A = alg.vec(n)

Returns a zero-initialized array of Lua numbers of length n with n rows and 1 column.

See alg.typeof() for arrays of different element type.

A = alg.tovec(t)

Returns an array of Lua numbers with 1 column initialized with the elements of t which must be a Lua table sequence:

local x = alg.tovec{1, 2, 3}

print(x)
--> +1.000000
--> +2.000000
--> +3.000000

See alg.typeof() for arrays of different element type.

A = alg.mat(nrow, ncol)

Returns a zero-initialized array of Lua numbers of length nrow*ncol with nrow rows and ncol columns.

See alg.typeof() for arrays of different element type.

A = alg.tomat(tt)

Returns an array of Lua numbers initialized with the elements of tt which must be a Lua table sequence of Lua table sequences of constant length:

local x = alg.tomat{{1, 2, 3},
                    {4, 5, 6}}

print(x)
--> +1.000000,+2.000000,+3.000000
--> +4.000000,+5.000000,+6.000000

See alg.typeof() for arrays of different element type.

A = alg.join(...)

Returns an array which is obtained by stacking the passed arrays vertically. Horizontal stacking is obtained via the .. concatenation operator:

local x1, x2 = alg.tovec{1, 2, 3}, alg.tovec{4, 5}
local x1x2 = alg.join(x1, x2)

print(x1x2)
--> +1.000000
--> +2.000000
--> +3.000000
--> +4.000000
--> +5.000000

local x, y = alg.tomat{{1, 2}, {3, 4}}, alg.tomat{{5, 6}, {7, 8}}
local xy = alg.join(x..y..x, y..x..y)
print(xy)
--> +1.000000,+2.000000,+5.000000,+6.000000,+1.000000,+2.000000
--> +3.000000,+4.000000,+7.000000,+8.000000,+3.000000,+4.000000
--> +5.000000,+6.000000,+1.000000,+2.000000,+5.000000,+6.000000
--> +7.000000,+8.000000,+3.000000,+4.000000,+7.000000,+8.000000

Please notice that .. does not create a new array and its use is limited to the alg.join() efunction.

A2 = A:new()

Returns a zero-filled array with the same dimensions of A.

A2 = A:copy()

Returns an array which is an independent copy of A.

n = #A

Returns the length of A.

n = A:nrow()

Returns the number of rows of A.

n = A:ncol()

Returns the number of columns of A.

el = A[i]; A[i] = el

Gets / sets the element at position i, access is left-to-right, top-to-bottom and indexing starts at 1.

el = A[{r,c}]; A[{r,c}] = el

Gets / sets the element at row r and column c, indexing starts at 1.

t = A:totable()

Returns a Lua table sequence of Lua table sequences which are initialized with the elements of A.

s = tostring(A)

Returns a string representation of the array.

alg.mul(C, A, B, At, Bt)

The matrix C is overwritten with the result of the matrix-matrix multiplication A times B. The remaining arguments are booleans specifying weather the multiplication is of transposed matrices or not. Note that matrix-vector multiplication is a special case of matrix-matrix multiplication and hence already covered by this function. The element type of all matrices involved must be constant and one of double, float, complex, complex float .

alg.pow(C, A, s)

The square matrix A is exponentiated to the s-th power, where s is a positive integer, and the matrix C is overwritten with the result. The element type of all matrices involved must be constant and one of double, float, complex, complex float.

algct = alg.typeof(ctype)

Takes as input a FFI ctype or a FFI ctype string declaration and returns a table algct with keys algct.vec(), algct.tovec(), algct.mat(), algct.tomat() where these constructors return arrays of element type ctype.

The argument ctype must be one of the following value types:

bool
char             
int8_t           
int16_t          
int32_t          
int64_t          
uint8_t          
uint16_t         
uint32_t         
uint64_t         
float            
double           
complex float
complex  

or the dual number diff.dn from the automatic differentiation module sci.diff.

The constructors for arrays of doubles (i.e. Lua numbers) have already been defined in the alg module:

-- These are already available in alg:
alg.vec   = alg.typeof('double').vec
alg.tovec = alg.typeof('double').tovec
alg.mat   = alg.typeof('double').mat
alg.tomat = alg.typeof('double').tomat