Quick start to Lyla
Introduction
This document provides a very short introduction to the basics of the Lyla linear algebra library. For a more detailed introduction, see the reference documentation.
Matrix and Vector classes
Lyla is a template class library for matrices and vectors. Each matrix or vector is represented by a specific class with different properties. There are classes for dense and sparse matrices holding (almost) arbitrary elements as doubles or floats. There are some predefined and frequently used classes:
- Dense rectangular matrices with double values:
lyla.Matrix.DMatrix, lyla.Matrix.DRowMatrix, lyla.Matrix.DColMatrix, - Dense vectors with double values:
lyla.Vector.DVector, lyla.Vector.DRowVector, lyla.Vector.DColVector, - Sparse rectangular matrices with double values:
lyla.Matrix.DRowSparse, lyla.Matrix.DColSparse, - Sparse vectors with double values:
lyla.Vector.DRowSparseVector, lyla.Vector.DColSparseVector.
All matrix classes derive from the base class lyla.Matrix, all vector classes derive from lyla.Vector. Those two base classes provide a common interface for all matrices and vectors. The interface can be divided into six main parts:
We will give a short introduction to all parts.
Status information
Each matrix has the following properties:
- Returns the number of rows.
size_t n_rows()
- Returns the number of columns.
size_t n_cols()
- Returns the orientation of the matrix, either Orientation.RowMajor or Orientation.ColMajor.
Orientation orientation()
Each vector has the following additional properties:
- Returns the size of the vector.
size_t length()
Initialization
Each matrix and vector class contains many initialization functions. Those are used to create a matrix with a specific dimension and initial values of the elements. The following functions are provided:
- Initializes the matrix to be empty.
void init()
- Initializes the matrix to have n rows, m columns and all values equal to x.
void init( size_t n, size_t m, value_t x )
- Initializes the matrix to have n rows, m columns and element (i,j) has value init_func(i,j).
void init( size_t n, size_t m, value_t delegate(size_t, size_t) init_func )
- Initializes the matrix to have n rows, m columns and element (row_indices[i], col_indices[i]) has value value[i] for all i. If some element has no value given in those array, it's initialized to 0.
void init( size_t n, size_t m, size_t[] row_indices, size_t[] col_indices, value_t[] values )
- Initializes the matrix to be a copy of A.
void init( Matrix!(value_t) A )
Vectors provide the same interface but they do not accept arbitrary dimensions (e.g. a column-vector requires the column-dimension to be 1). In addition, vectors provide the following functions.
- Initializes the vector to have n elements and all values equal to x.
void init( size_t n, value_t x )
- Initializes the vector to have n elements and element i has value init_func(i).
void init( size_t n, value_t delegate(size_t) init_func )
- Initializes the vector to have n elements and element indices[i] has value value[i] for all i. If some element has no value given in those array, it's initialized to 0.
void init( size_t n, size_t[] indices, value_t[] values )
- Initializes the vector to be a copy of values.
void init( value_t[] values )
Furthermore each matrix and vector class provides constructors with the same signatures as the init-functions above. They can be used to create a matrix/vector and initialize it in one step.
Examples:
1. Create an empty matrix.
auto A = new DMatrix;
2. Create a 2x5 matrix and initialize all elements to 42.
auto A = new DMatrix(2, 5, 42);
3. Initialize an existing vector of class DVector to contain the elements 1, 2, ..., 23.
auto x = new DVector; ... x.init( 23, (size_t i) { return cast(double)(i + 1); });
Element access
This are the common operations one expects from a matrix:
- Returns the element with row-index i and column-index j.
value_t opIndex( size_t i, size_t j )
- Sets the value of the element with row-index i and column-index j to x.
value_t opIndexAssign( value_t x, size_t i, size_t j )
Vectors support one dimensional-access:
- Returns the element with index i.
value_t opIndex( size_t i )
- Sets the value of the element with index i to x.
value_t opIndexAssign( value_t x, size_t i )
Examples
1. Increase all elements by 1.
auto A = new DMatrix; ... for ( size_t i = 0; i < A.n_rows; ++i ) { for ( size_t j = 0; j < A.n_cols; ++j ) { A[i,j] = A[i,j] + 1; } }
2. Set vector-elements by matrix-access, get them by vector-access.
auto x = new DColVector( 42 ); for ( size_t i = 0; i < x.n_rows; ++i ) { for ( size_t j = 0; j < x.n_cols; ++j ) { x[i,j] = i + j; } } for ( size_t i = 0; i < x.length; ++i ) { assert( x[i] == i ); }
Iteration
Each matrix provides two kinds of iterators: Dense iterators, iterating over each logical matrix element, and non-zero iterators, iterating only over those elements A[i,j] != 0. The latter ones are mainly useful with sparse matrices.
Iterators are available for each row and each column as well as for all matrix elements. The following functions provide iterators:
- Returns an iterator over all matrix elements.
ModifyMatrixIterator!(value_t) elements()
- Returns an iterator over all non-zero matrix elements.
MatrixIterator!(value_t) nz_elements()
- Returns an iterator over all elements of row i.
ModifyIndexIterator!(value_t) row_elements( size_t i )
- Returns an iterator over all non-zero elements of row i.
IndexIterator!(value_t) nz_row_elements( size_t i )
- Returns an iterator over all elements of column j.
ModifyIndexIterator!(value_t) col_elements( size_t j)
- Returns an iterator over all non-zero elements of column j.
IndexIterator!(value_t) nz_col_elements( size_t j )
All iterators are inspired by tango collections. The first two functions returns matrix iterators which have the following interface:
- Returns the value of the current element.
value_t value()
- Returns the row-index of the current element.
size_t row()
- Returns the column-index of the current element.
size_t col()
- Returns true if there's another element and moves to that element, returns false otherwise.
bool more()
- Iterates over the (remaining) elements of this iterator.
int opApply(int delegate(ref size_t i, ref size_t j, ref value_t x) )
The other functions return index iterators with the following interface:
- Returns the value of the current element.
value_t value()
- Returns the index of the current element in its part (row/column).
size_t index()
- Returns true if there's another element and moves to that element, returns false otherwise.
bool more()
- Iterates over the (remaining) elements of this iterator.
int opApply(int delegate(ref size_t i, ref value_t x) )
The iterators returns by the nz_* functions are read-only. In contrast the (Modify*Iterator) iterators returned by the other functions provide write-access to the element the point to, i.e. they have the following additional function:
- Sets the value of the current element to x.
value_t value( value_t x )
The two functions elements and nz_elements for vectors have a slightly extended interface:
- Returns an iterator over all vector elements.
ModifyVectorIterator!(value_t) elements
- Returns an iterator over all non-zero vector elements.
VectorIterator!(value_t) nz_elements()
Vector iterators provide the following additional functions which deal with the one-dimensional index of the current element:
- Returns the one-dimensional vector index of the current element.
size_t index()
- Iterates over the (remaining) elements of this iterator, where i is the element's vector index.
int opApply(int delegate(ref size_t i, ref value_t x) )
Examples
1. Iterate over all elements of the 4th column:
A = new DMatrix; ... // using an explicit iterator for ( auto it = A.col_elements(4); it.more; ) { it.value = it.value + 1; } // using foreach foreach ( size_t i, ref double x; A.col_elements(4) ) { x -= 1; }
2. Iterator over non-zero elements:
for ( auto it = A.nz_elements; it.more; ) { writefln("%d %d %g", it.row, it.col, it.value); } foreach ( size_t i, size_t j, double x; A.nz_elements ) { writefln("%d %d %g", i, j, x); }
3. Iterate over all elements with |.| > 5 using a tolerance-iterator:
for ( auto it = tol_it(A.nz_elements, 5); it.more; ) { writefln("%d %d %g", it.row, it.col, it.value); }
4. Display a matrix row-wise
for ( size_t i = 0; i < A.n_rows; ++i ) { foreach( size_t j, double x; A.row_elements(i) ) { writef("%g ", x); } writefln(); }
Modification
Rectangular matrices and vectors may be modified by adding and removing elements. Matrices provide the following modification functions:
- Insert the rows of A before row i
void insert_rows( size_t i, Matrix!(value_t) A ) void prepend_rows( Matrix!(value_t) A ) // same as insert_rows( 0, A ) void append_rows( Matrix!(value_t) A ) // same as insert_rows( n_rows, A );
- Insert the columns of A before column j
void insert_cols( size_t j, Matrix!(value_t A) ) void prepend_cols( Matrix!(value_t) A ) // same as insert_cols( 0, A ) void append_cols( Matrix!(value_t) A ) // same as insert_cols( n_cols, A );
- Delete rows i, i+1, ..., i+n-1
void delete_rows( size_t i, size_t n )
- Delete rows del_indices[0], del_indices[1], ...
void delete_rows( size_t[] del_indices )
- Delete columns j, j+1, ..., j+n-1
void delete_cols( size_t j, size_t n )
- Delete columns del_indices[0], del_indices[1], ...
void delete_cols( size_t[] del_indices )
Examples:
1. Duplicate the rows of A
A.append_rows( A );
2. Delete columns 0,2,4,...
auto del_indices = new size_t[A.n_cols / 2]; foreach(i, x; del_indices) x = i*2; A.delete_cols( del_indices );
Operations
LYLA support many matrix operations. For example, adding two matrices B and C and storing the result in a matrix A is done by
// A = B + C add(A, B, C);
As a general rule, the result of an operation is always stored in the first argument of the function. Furthermore, for many operations it is possible to call the function as method call of the result matrix. Therefore, the above example is equivalent to
// A = B + C A.add(B, C);
Many functions may store the result in one of the operands:
// A = A + B add(A, B); A.add(B);
The same is true for unary matrix functions, i.e. functions that work on a single matrix:
// A[i,j] = sqrt(B[i,j]) sqrt(A, B); A.sqrt(B); // A[i,j] = sqrt(A[i,j]); sqrt(A); A.sqrt();
This table lists some of the binary operations for matrices:
add(A,B,C) | A = B + C |
sub(A,B,C) | A = B - C |
div(A,B,C) | A[i,j] = B[i,j] / C[i,j] |
pmul(A,B,C) | A[i,j] = B[i,j] * C[i,j] |
mul(A,B,C) | A = B * C |
This table lists some of the binary operations for matrices with scalars:
add(A,x,C) | A[i,j] = x + C[i,j] |
add(A,B,x) | A[i,j] = B[i,j] + x |
sub(A,x,C) | A[i,j] = x - C[i,j] |
sub(A,B,x) | A[i,j] = B[i,j] - x |
div(A,x,C) | A[i,j] = x / C[i,j] |
div(A,B,x) | A[i,j] = B[i,j] / x |
pmul(A,x,C) | A[i,j] = x * C[i,j] |
pmul(A,B,x) | A[i,j] = B[i,j] * x |
Unary operations:
abs(A,B) | A[i,j] = |B[i,j]| |
neg(A,B) | A[i,j] = -B[i,j] |
sqrt(A,B) | A[i,j] = sqrt(B[i,j]) |
sin(A,B) | A[i,j] = sin(B[i,j]) |
cos(A,B) | A[i,j] = cos(B[i,j]) |
Operators
Common matrix operators are overloaded: {{{#!d auto A = new DMatrix, B = new DMatrix; double x; ... A += B; // same as A.add(B); A += x; // same as A.add(x); A -= B; // same as A.sub(B); A -= x; // same as A.sub(x); A *= B; // same as A.mul(B); A *= x; // same as A.mul(x); A /= B; // same as A.div(B); A /= x; // same as A.div(x); }}}
Binary operators are overloaded, too, but there is one important difference to the explicit functions above. Since operators return a value instead of modifying their arguments, calling for example A + B always creates a new matrix:
auto A = new DMatrix, B = new DMatrix; ... auto C = A + B; // creates a new matrix C
The resulting matrix has usually the same orientation as the first operand and is sparse if and only if both operands are sparse. This means, if you use only one matrix type, the result of the operations will always have the same type as both operands:
auto A = new DMatrix, B = new DMatrix; DMatrix C; ... C = A * B; // works
BLAS like operations
LYLA support some BLAS-like functions. If BLAS-support is compiled in, the arithmetic for dense matrices call the BLAS-functions to perform the operations. Since the rule in LYLA is to store the result in the first argument (in contrast to BLAS where the result is usually stored in the last argument), must functions have another order of parameters than their BLAS counterparts:
- compute the inner product of two vectors x and y
value_t dot(Vector!(value_t) x, Vector!(value_t) y)
- scale a vector y by factor alpha: x = x * alpha
void scal( Vector!(value_t) x, value_t alpha )
- compute x = x + alpha * y
void xepay( Vector!(value_t) x, value_t alpha, Vector!(value_t) y)
- Multiply a matrix and a vector: x = x * alpha + A * y * beta
void vmul( Vector!(value_t) x, value_t alpha, Matrix!(value_t) A, Vector!(value_t) y, value_t beta )
- Multiply two matrices: A = A * alpha + B * C * beta
void mmul( Matrix!(value_t) A, value_t alpha, Matrix!(value_t) B, Matrix!(value_t) C, value_t beta )
Remark: BLAS functions may be called as methods of the destination matrix, too:
/// A = A * alpha + B * C * beta A.mmul( alpha, B, C, beta )