Matrico matrix Functor

`matrico`’s Matrix Functor Module `matrix`

The matrix functor is the abstract definition of matrico’s core matrix type, which is a back-end module and particularly not user-facing, even though a subset of functions will be passed through to the user interface. This setup has two reasons:

Generic functional procedures on the matrix type simplify higher-level function definition,
and this allows to hide the low-level matrix type definition from the user and thus ensure consistency of the matrix data.

Furthermore, the matrix type is encapsuled in a functor, instead of a standard module, to keep the underlying entry-type of the matrix flexible. So while, matrico will by default use double-precision flonums, using single-precision flonums, rationals or integers instead remains possible. Beyond entry type, also the column type is kept exchangeable this way. This means instead of the standard dense matrix format, also a sparse matrix format could be implemented. The source file src/matrix.scm is here.

Matrix Type

Internally a matrix in matrico is designed as a record or structure via the (chicken base) module. The matrix itself is encoded as list of vectors, specifically a list of columns - a column-major matrix data model. This is in opposition to the typically chosen plain vector with strided rows or columns, or vector of vectors. The reasons for this uncommon choice are again twofold:

Using a list plays on the strength of Scheme while keeping each column (a vector) accessible in constant time.
For larger problems it is not necessarily preferable to use one large memory block as matrix. The list-of-columns approach allows spreading the columns throughout available memory and thus making allocation easier and improve balancing memory use.

It will be important for higher level functions using this matrix type to keep this underlying list-of-columns layout in mind. This means minimizing list traversal and maximizing per-column work.

As the matrix dimensions are constant after binding, the number of columns remains constant, hence it can be stored and no list traversal is needed to obtain this number. Thus, the matrix record is defined with two slots,

the list-of-columns storing the matrix entries,
and the number of columns,

yielding the implicitly defined functions:

(make-matrix cols lst) - Returns a matrix with list-of-columns argument and number of columns argument.
(matrix-cols mat) - Returns the number columns of the matrix argument.
(matrix-data mat) - Returns the list-of-columns payload of the matrix argument.
(matrix? any) - Returns a boolean answering if argument is a matrix.

The make-matrix function illustrates why this function cannot be user facing. Would it be, it is not guaranteed that a list of columns of equal length or even a list at all is passed as first slot, nor that the correct number of columns or a number at all is passed as second slot. The matrix functor takes care of assigning these correctly.

Local Functions

A certain composition of procedures appears in quite a few matrix functions, hence this was extracted and made into two local functions. These functions are incomplete map operations of matrix arguments.

(matrix-map* fun x) - Returns matrix resulting from applying function argument to each column of the matrix argument.
(matrix-map** fun x y) - Returns matrix resulting from applying function argument to each column of the two matrix arguments.

Note that these are not complete map operations on matrices, as only the matrix functor knows about the map operation on a column. Also, for performance reasons, two separate functions were made instead of one with variadic arguments.

Generators

To create matrices, generators or constructors are used. First, a basic function to specify dimensions in terms of rows and columns and the value of all entries, is available. Sometimes one wants to prescribe a complete matrix manually entry-by-entry, for this purpose a generator from a list-of-list, in row-major format, which is more natural for terminal input or reading from text files, is provided. To procedurally assemble a matrix, a function is included which calculates an entry based on its row and column index.

(make-matrix* rows cols val) - Returns matrix of row and column argument dimension with all entries set to value argument.
(make-matrix** lst) - Returns matrix from list-of-lists argument.
(matrix-generate fun rows cols) - Returns matrix resulting from computing entries by applying function argument to entries’ row and column indices.

I use the stars to signify variants of the standard generator make-matrix which is automatically provided through creating the matrix record.

Combiners

While the generator functions create matrices from non-matrix arguments, combiners create matrices from matrices, particularly by concatenation. This is implemented by horizontal concatenation, meaning joining lists of all matrix arguments (having the same number of rows), and vertical concatenation, meaning joining all respective columns of all matrix arguments (having the same number of columns).

(matrix-horcat . mat) - Returns matrix resulting from concatenating list-of-matrixes argument horizontally.
(matrix-vercat . mat) - Returns matrix resulting from concatenating list-of-matrixes argument vertically.

Dimensions

To determine the dimensions of a matrix, a set of functions is provided. These include the number of rows, the number of entries, and if a matrix has the shape of a mathematical scalar, vector or matrix; the latter is sometimes known as rank of a multi-dimensional array, but this is only correct in the sense of the tensor rank, and not as “rank” typically used in linear algebra.

(matrix-rows mat) - Returns the number of rows of the matrix argument.
(matrix-numel mat) - Returns the total number of entries of the matrix argument.
(matrix-dims mat) - Returns 0 for 1x1, 1 for Nx1 or 1xN, or 2 for NxM dimensions of matrix argument.

I note again that the matrix-cols function is automatically provided through creating the matrix record.

Accessors

There are various ways matrices can be accessed. Most basic is access to a single element; I also provide a specific accessor to the entry at the first row and first column, as this entry is always available and used for converting between a scalar matrix and its entry type. Also a function to return an element form the first column is included, not for performance as one may suspect, but readability, as it highlights that using this call no list traversal is needed. To manipulate single entries there is a function which returns a copy of the argument matrix, yet with a changed entry, and a function mutating the entry of the argument matrix. Beyond single entries, a single column, a single row or a sub-matrix may need to be extracted, for which functionality is included, the latter needs the sublist function form the utils module. Another important accessor for matrices returns the matrix’s diagonal - all entries with the same row and column index of a square matrix.

(matrix-ref00 mat) - Returns the entry of the first column and the first row of the matrix argument.
(matrix-ref*0 mat row) - Returns the entry of the first column at the row argument of the matrix argument.
(matrix-ref mat row col) - Returns the entry at the row argument and column argument of the matrix argument.
(matrix-set mat row col val) - Returns matrix with entry at the row argument and column argument of the matrix argument to value argument.
(matrix-set! mat row col val) - Returns void, mutate entry at the row argument and column argument of the matrix argument to value argument.
(matrix-col mat col) - Returns (column) matrix being at the column index argument of the matrix argument.
(matrix-row mat row) - Returns (row) matrix being at the row index argument of the matrix argument.
(matrix-diag mat) - Returns (column) matrix being the diagonal of the square matrix argument.
(matrix-submatrix mat row1 row2 col1 col2) - Returns matrix being the entries between row index arguments and column index arguments of matrix argument.

Predicates

Given a matrix, one may want to test its shape: Is it a one-column matrix or one-row matrix, is it a vector so either a column matrix or row matrix, or a scalar, meaning a column matrix and row matrix. Furthermore, I may want to know if a matrix has the same number of rows and columns. Given two matrices, I may want to know if those have they same number columns, rows or both. Lastly, it may be interesting to test for a provided predicate if any or all entries of a matrix fulfill it.

(matrix-col? mat) - Returns boolean answering if matrix argument has only one column.
(matrix-row? mat) - Returns boolean answering if matrix argument has only one row.
(matrix-scalar? mat) - Returns boolean answering if matrix argument has only one column and one row.
(matrix-vector? mat) - Returns boolean answering if matrix argument has only either one column or one row.
(matrix-square? mat) - Returns boolean answering if matrix argument has the same number of rows and columns.
(matrix-samecols? x y) - Returns boolean answering if matrix arguments have same number of columns.
(matrix-samerows? x y) - Returns boolean answering if matrix arguments have same number of rows.
(matrix-samedims? x y) - Returns boolean answering if matrix arguments have same number of columns and same number of rows.
(matrix-any? pred mat) - Returns boolean answering if any entry of matrix argument fulfills predicate argument.
(matrix-all? pred mat) - Returns boolean answering if all entries of matrix argument fulfill predicate argument.

I note that the matrix? predicate is automatically provided through creating the matrix record.

Reducers

There are three reducers or folders required for a two-dimensional array, such as a matrix:

folding all columns to a row matrix,
folding all rows to a column matrix,
and folding all columns and rows to a scalar,

which are implemented as three functions, as opposed to a single function with an extra argument, due to the very different internals of these functions and their non-uniform return type.

(matrix-colfold fun ini mat) - Returns (row) matrix resulting from folding all columns of matrix argument pairwise by function argument starting with initial argument.
(matrix-rowfold fun ini mat) - Returns (column) matrix resulting from folding all rows of matrix argument pairwise by function argument starting with initial argument.
(matrix-allfold fun ini mat) - Returns entry-type resulting from folding row of all folded column of matrix argument by function argument starting with initial argument.

Mappers

The unary map of a matrix is just a composition of the outer list map and the inner column map, applying the argument function to each entry. The binary map is more involved, as not only matrices of same dimension should be mappable, but also two matrices that fall in one of the cases:

one argument is scalar, the other has any dimension,
both arguments share either column or row dimension,
or both arguments are vectors.

This functionality is called “automatic brodcasting” (in Octave, NumPy, Julia) or “implicit expansion” (in MATLAB).

(matrix-map fun mat) - Returns matrix resulting from applying function argument to all matrix argument’s entries.
(matrix-broadcast fun x y) - Returns matrix resulting from applying function to matrix arguments under broadcasting rules.

Specific

There a some operations specific to matrices, or rather two-dimensional arrays, and independent of the type of the entries that are generically provided in this functor module. There is (mathematical) vectorization, meaning stacking all columns of a matrix, and transposition of a matrix, which swaps row and column indices of each entry. Whereas the transposition is implemented elegantly using a similar approach to list-of-list matrix transposition.

The vector “axpy” operation is here generalized to matrices. This is neither a broadcasted operation due to performance reasons nor a plain mapper due to the non-matrix argument.

Given two column matrices with the same number of rows, the scalar product is an essential operation in (vector) linear algebra, also known as inner product. The outer product counter part is not defined here, as it can be more easily achieved by a broadcasting / expansion operation. In a sense, matrix multiplication is then a generalization if the scalar product for two matrices (where the number or columns of the first operand is equal to the number of rows of the second operand), as it computes all possible scalar products of all rows of the first operand and all columns of the second operand. I call the matrix multiplication dot-product as NumPy does, since it is mathematically more expressive. The implemented function is called dot* (dot-star, a special variant of dot) though, because the first argument is assumed to be transposed: This allows an efficient implementation for the list-of-column matrix data model. This may seem a little odd as there is no standard version of dot, but the intention here is to use an efficient algorithm by default, and the frontend can then offer a standard dot matrix maultiplication and transpose the arguments as needed.

Lastly, for performance and convenience, twin functions to “explode” a matrix into a list of column matrices, and “implode” a list of column matrices (with the same number of rows) into one matrix, where the latter wraps the matrix-horcat function are included.

(matrix-vec mat) - Returns (column) matrix resulting form concatenating columns of matrix argument vertically.
(matrix-transpose mat) - Returns matrix with entries of argument matrix, where the entries’ row and column indices are exchanged.
(matrix-axpy a x y) - Returns matrix resulting from applying “axpy” function to scalar and matrix arguments.
(matrix-scalar xt y) - Returns entry-type resulting from computing scalar product of column matrix arguments with same number of rows.
(matrix-dot* xt y) - Returns matrix resulting from matrix multiplication of the matrix arguments, where the first argument is assumed transposed.
(matrix-explode mat) - Returns list-of-column-matrixes from matrix argument.
(matrix-implode lst) - Returns matrix from list-of-column-matrixes argument.

Utilities

Lastly, some matrix utility functions are needed, these have, by nature, side-effects. Most importantly is printing a matrix to the standard output, so typically the terminal. This functions needs a function argument that processes the internal entry-to-string conversion, so that each printed entry has the same character width, which is especially necessary for flonums, see fptaper from the fpmath module. This print function utilizes the “square bracket” symbols from the miscellaneous technical unicode block.

Also, matrices need to be stored and restored from disk; this can be done in a native format using S-expressions. Additionally, matrices need to be exported, so other numerical packages, such as MATLAB, Octave, NumPy or Julia can load them. The common format for all of these is the CSV (comma-separated-value) format.

(matrix-print fun mat) - Returns void, prints argument matrix on the standard output, preprocessing matrix entries by function argument.
(matrix-export str mat) - Returns void, output argument matrix to CSV file at string argument path.
(matrix-save str mat) - Returns void, output argument matrix to S-expression formatted file at string argument path.
(matrix-load str) - Returns matrix loaded from string argument path.

In the next post, I will focus on the (chicken memory representation) which provides an essential function for the matrix functions matrix-horcat, matrix-vercat, matrix-col, matrix-submatrix, matrix-set, matrix-vec, and matrix-explode, not discussed here.

matrico’s Matrix Functor Module matrix