Matrico Module IIII
matrico’s matrico Module (Part IIII)
Series:
For the fourth part in this series of posts about the main matrico module of the matrico project,
I summarize the core numerical functionality.
These are again contained in the src/mx.scm Scheme source code file,
that is included into the module’s main file matrico.scm.
Linear Algebra
Usual basic linear algebra functions that matrico provides are concatenation of matrices,
vectorization and transposition of matrices, symmetric and anti-symmetric (skew-symmetric) parts of a matrix,
and extracting the diagonal of a matrix:
(mx-horcat A B)- Returns matrix consisting of a copy of horizontally concatenated argument matrixes.(mx-vercat A B)- Returns matrix consisting of a copy of vertically concatenated argument matrixes.(mx-vec A)- Returns column-matrix with entries of argument matrix column-wise stacked.(mx-transpose A)- Returns matrix with entries of argument matrix with swapped row and column indices.(mx-sympart A)- Returns matrix being the symmetric part of argument matrix.(mx-skewpart A)- Returns matrix being the skew-symmetric part of argument matrix.(mx-diagonal A)- Returns diagonal matrix argument column matrix.
Linear Problems
A core competency of numerical linear algebra is the solution of linear problems,
meaning solving a system of linear equations, typically presented in vectorized form
as matrix times unknown vector equals vector: A x = b.
A class of algorithms to solve such linear problems utilize matrix decompositions.
The idea behind these approaches is to factorize the matrix A such that the factors
facilitate a solution.
Among the available matrix decomposition algorithms, I selected the QR decomposition (or QR factorization).
This algorithm factors the matrix A = QR into an orthogonal matrix Q and a right triangular matrix R.
While the former (Q) can be inverted by transposition,
the latter (R) corresponds to back substitution;
this means A = QR => QR x = b => R x = Q^T b.
Advantages of the QR decomposition compared to other decompositions are on the one hand, that the matrix A
does not need to be symmetric, invertible, or even square (in the singular or rectangular case a least-squares solution is obtained),
on the other hand, the QR decomposition is the basis for more complex algorithms,
such as the Singular Value Decomposition (SVD), or eigenvalue computation via the Francis-QR algorithm.
Among the various variants of practical QR factorization algorithms, I specifically chose the Gram-Schmidt-based Schwarz-Rutishauser algorithm,
also known as column-wise modified Gram Schmidt;
because first, this column-wise algorithm fits matrico’s matrix data structure,
second, its lower complexity, third, its numerical stability, and lastly its simple implementation.
matrico provides functions computing QR decompositions for different applications,
all using the same core algorithm:
(mx-qr A)- Returns a pair of matrixes (Q and R) for an argument matrix(mx-solver A)- Returns a function (expecting a column matrix argument) that solves the linear problem formed with the argument matrix.(mx-solve A b)- Returns a column matrix solving the linear problem formed by argument matrix and column matrix.(mx-orth A)- Returns matrix representing an orthogonalization of argument matrix.
Use mx-qr if no linear problem is to be solved, but the factor matrices are required.
Use mx-solver if multiple linear problems need to be solved with the same matrix A but different right-hand-sides b, as the returned function caches the QR decomposition as closure and thus prevents recomputation.
Use mx-solve if only a single solution to a linear problem is required.
Determinants
The QR decomposition also allows to compute a matrix’s determinant up to its sign. Thus, the absolute value of the determinant is computable. Given a positive definite matrix, also the logarithm of the determinant can be computed, which is useful for badly conditioned matrices.
(mx-absdet A)- Returns the absolute value of the determinant of the argument matrix.(mx-logdet A)- Returns the logarithm of the determinant of the argument matrix.
Traces
Functions computing accumulated quantites from the matrix diagonal are subsumed under “Traces”. Usually, a matrix trace refers to the sum of its diagonal values; here, also a the product of diagonal values is provided, as well as functions calculating the (additive) trace of a product of matrices, since the trace of a product can be computed more efficiently than computing first the product and then its trace.
(mx-trace A)- Returns the sum of diagonal entries of the matrix argument.(mx-multrace A)- Returns the product of diagonal entries of the matrix argument.(mx-prodtrace A B)- Returns the sum of diagonal entries of the product of matrix arguments.(mx-prodtrace* AT B)- Returns the sum of diagonal entries of the product of matrix arguments, with first factor transposed.
Matrix multiplication
Matrix multiplication is also an important pillar of (numerical) linear algebra.
It has interpretations, among others, as compositions of linear transformations,
or all possible inner products of two sets vectors (column matrices).
matrico provides variants of the generic matrix multiplication
mainly for convenience:
(mx-scalar AT B)- Returns scalar product (inner product) of column-matrix arguments.(mx-dyadic A B)- Returns dyadic product (outer product) of column-matrix argument and row-matrix argument.(mx-dot* AT B)- Returns product of matrix arguments with first factor transposed.(mx-dot A B)- Returns product of matrix arguments.(mx-gram A)- Returns product of transposed matrix argument with itself.(mx-gram* A)- Returns product of matrix argument with its transposed self.(mx-square A)- Returns product of matrix argument with itself.
The mx-dot* functions should always be preferred for non-column-matrices, as it is the most performant.
Multivariate Statistics
Some statistical functions are available in matrico, too:
cross-covariance,
cross-correlation,
variance, standard deviation,
angle, and coherence.
For these functions the argument matrices are assumed to be arranged as
each row corresponding to all observations of one variable,
while each column represents all variables of one observation:
(mx-xcov A B)- Returns cross-covariance matrix of matrix arguments.(mx-cov A)- Returns covariance of matrix argument.(mx-var A)- Returns variance of matrix argument.(mx-std A)- Returns standard deviation of matrix argument.(mx-xcor A B)- Returns cross-correlation of matrix arguments.(mx-cor A)- Returns correlation of matrix argument.(mx-angle A B)- Returns angles between matrix argument columns.(mx-coher A B)- Returns coherence of matrix arguments.
Analysis and Calculus
Analysis deals with derivatives and integrals,
so numerical analysis approximates derivation and integration.
matrico provides a function to compute differences of consecutive columns of a matrix,
which can be used for a finite difference approximation of a difference quotient.
To approximate integrals by Riemann sums, a column-wise trapezoid rule is implemented.
(mx-diff A)- Returns matrix of differences of consecutive columns of argument matrix.(mx-trapz A)- Returns column-matrix of sums of consecutive columns of argument matrix.
Ordinary Differential Equations
A central subject of numerical analysis is the approximation of solutions to differential equations.
matrico is focusing on ordinary differential equations (ODE), differential equations in a one dimensional variable,
as opposed to partial differential equations in multiple dimensional variables.
Ordinary differential equations are prominently used for dynamical systems,
where the variable is (one-dimensional) time, and describes the temporal evolution of a system.
Integrators, such as Runge-Kutta methods,
compute a sequence of point-wise approximations to the solution function’s values, called trajectory for evolution systems.
In matrico trajectories are obtained by a single step methods, meaning the next (time) step approximation depends only the current approximation step.
Opposed to, for example, MATLAB, matrico uses fixed step methods,
and not adaptive time step methods.
The two selected solvers implemented in matrico are explicit Runge-Kutta method of second (nonlinear) order,
and were chosen as they feature low complexity in terms of compute and memory.
Particularly, both methods are of the “two register” class (“2R*”);
this means that internally an iteration requires merely two vector (column-matrix) registers,
of which one is constant for each sub-step iteration.
This is the minimum memory consumption possible.
Furthermore, both time steppers induce a family of methods as either allows to increase the number of sub-steps or stages, to increase stability,
and hence approximate more complicated problems:
First, the family of implemented hyperbolic Runge-Kutta methods is aimed at transport, advection, or non-normal problems and require a stability region around the imaginary axis.
These hyperbolic Runge-Kutta methods have the maximum hyperbolic stability limit for an explicit second order method of chosen number of stages.
Second, the family of implemented strong stability preserving Runge-Kutta methods is aimed at general and parabolic problems, which require a stability region in the negative-real half-plane.
Strong stability preserving (SSP) means that if a first order explicit Euler method adheres to some constraint for the given ODE, so will this second-order method, which is optimal among second-order methods in this regard.
(mx-ode2-hyp num sys tim x0)- Returns matrix holding trajectory solved by the hyperbolic Runge-Kutta method with fixnum argument many stages, procedure argument vector field, pair argument variable discretization (step, horizon) and column-matrix initial value.(mx-ode2-ssp num sys tim x0)- Returns matrix holding trajectory solved by the strong stability preserving Runge-Kutta method with fixnum argument many stages, procedure argument vector field, pair argument variable discretization (step, horizon) and column-matrix initial value.
Instead of just a vector field also a pair of a vector field procedure and an output function procedure can be passed, which allows to simulate input-output systems, typical in system theory and control theory.
This sums up the highest-level functionality of matrico.
The next post will mop up utility and private functions for completeness,
and thereby conclude this series.