The mlina C++ library is a lightweight and easy starting point for applications that require computations with small matrices and vectors of fixed size. Apart from the fundamental arithmetic operations, mlina supports solving dense linear systems and features dedicated types and algorithms for geometric computations in three dimensions. The goal is to provide a value-centric, minimalist, and easy-to-learn API that does not require any configuration or specialized algorithmic knowledge.
mlina has been developed by Mario Laux.
mlina is a single-file header-only library with no external dependencies for C++11 or higher. You can simply #include "mlina.h" and you are good to go. Everything is defined within the namespace mlina, so you can either use using namespace mlina; or specify the namespace explicitly where needed. Documentation of the public API can be found in the header file, but the following tutorial is probably the better starting point.
The logo, this tutorial, and the header file are all published under the MIT license. Internal benchmarks and tests are not part of this public repository.
All matrices and vectors in mlina have a fixed size, which is also part of their type. All sizes and indices are of type std::size_t. Vectors are always treated as column vectors, i.e. they are matrices with exactly one column. Matrices store their entries directly as part of the object and no dynamic memory allocations are performed.
The entries of both matrices and vectors can be initialized upon construction by specifying all elements in a single initializer list row by row.
// construct a 2x3 matrix with entries 0, 1, 2 in the first
// row and entries 3, 4, 5 in the second row (the size is
// fixed, the formatting of the initializer list is just a
// visual aid)
Matrix<2, 3> M { 0, 1, 2,
3, 4, 5 };
// make a 3D (column-) vector with entries 6, 7, 8
// (Vector<n> is an alias for Matrix<n, 1>)
Vector<3> v { 6,
7,
8 };
// have a look at the contents of M and v
std::cout << "M = " << M << '\n'
<< "v = " << v << '\n';
// access (read & write) individual elements (all indices
// are zero-based)
std::cout << "M(1, 1) = " << M(1, 1) << '\n';
// set elements individually
M(1, 2) = 55;
std::cout << "after update: M = " << M << '\n';
// vectors can alternatively be accessed via the subscript
// operator
std::cout << "v[0] = " << v[0] << '\n';
// there are dedicated functions to construct zero- and
// identity matrices
auto Z = zero<2, 4>();
auto Id3 = identity<3>();
std::cout << "Z = " << Z << '\n'
<< "Id3 = " << Id3 << '\n';mlina uses double precision for all computations and does not perform any additional rounding internally. However, values with a very small absolute value are replaced with 0 for printing via operator<< (as they are most likely unwanted artifacts of floating point arithmetic):
Vector<6> v { 1.0, 1e-13, 1e-14, 1e-15, 1e-16, 1e-17 };
std::cout << v << '\n';Functions and mathematical operations in mlina never modify their arguments.
// construct example matrices
Matrix<2, 3> A { 2, -2, 1,
3, 6, -4 };
Matrix<2, 3> B { 1, 0, -6,
2, -1, 3 };
// construct example vectors (= matrices with one column)
Vector<3> x { 2, 0, 9 };
Vector<3> y { -3, 1, 8 };
// the usual arithmetic operations are supported (provided
// that all the dimensions match)
auto r1 = (0.5 * A + B / 2) * (x - y);
auto r2 = (A + B) * x/2 - 0.5 * (A * y + B * y);
// compare the two results
std::cout << "r1 = " << r1 << '\n'
<< "r2 = " << r2 << '\n';
// transposition works as usual (AT is a Matrix<3, 2> and
// xT is a Matrix<1, 3>)
auto AT = transpose(A);
auto xT = transpose(x);
std::cout << "A = " << A << '\n'
<< "x = " << x << '\n'
<< "A * x = " << A * x << '\n'
<< "x^T * A^T = " << xT * AT << '\n';mlina provides functions to compute the Euclidean norm of vectors, the distance between two points, and the inner product (these work in any dimension):
// vectors of matching dimension support the standard inner
// product (dot product), which yields a scalar
Vector<2> x { 5, 1 }, y { 2, -10 };
std::cout << "x = " << x << '\n'
<< "y = " << y << '\n'
<< "x * y = " << dot(x, y) << '\n';
// the inner product can be expressed as a matrix product,
// but will then yield a 1x1-matrix instead of a scalar
std::cout << "x^T * y = " << transpose(x) * y << '\n';
// compute the Euclidean norm (= "length") of a vector
Vector<2> v { 4, 5 };
std::cout << "v = " << v << '\n'
<< "|v| = " << norm(v) << '\n';
// compute the distance between two points
Vector<2> p { 3, -1 }, q { -2, 11 };
std::cout << "p = " << p << '\n'
<< "q = " << q << '\n'
<< "dist(p, q) = " << dist(p, q) << '\n'
<< "|p-q| = " << norm(p - q) << '\n';For two three-dimensional vectors, one can also take the cross product:
// there is a type alias Vec3 for three-dimensional vectors
Vec3 x { 2, -3, 5 }, y { 3, 1, 8 };
// "cross product" or "vector product"
auto z = cross(x, y);
// verify orthogonality
std::cout << "x = " << x << '\n'
<< "y = " << y << '\n'
<< "z = " << z << '\n'
<< "x * z = " << dot(x, z) << '\n'
<< "y * z = " << dot(y, z) << '\n';mlina provides a versatile solve function that can handle a variety of linear systems. Let's start with the most basic case:
// solve the linear system A * x = b given by
//
// 2 * x0 - 3 * x1 = -4
// 5 * x0 + 2 * x1 = 9
//
Matrix<2, 2> A { 2, -3,
5, 2 };
Vector<2> b { -4,
9 };
auto x = solve(A, b);
// verify the solution
std::cout << "A = " << A << '\n'
<< "x = " << x << '\n'
<< "b = " << b << '\n'
<< "A * x = " << A * x << '\n';The solve function can also handle matrix equations, i.e. it can solve linear systems for an arbitrary number of right-hand sides simultaneously. This is more efficient than multiple individual calls to solve as it has to analyze the coefficient matrix only once.
// solve A * X = B for X
Matrix<3, 3> A { 2, -1, 4,
6, 0, -1,
7, 1, 3 };
Matrix<3, 2> B { 13, 12,
24, 3,
23, 18 };
// the solution X is a Matrix<3, 2> (solve automatically
// rejects inputs that don't match in size)
auto X = solve(A, B);
std::cout << "A = " << A << '\n'
<< "X = " << X << '\n'
<< "B = " << B << '\n'
<< "A * X = " << A * X << '\n';If a square matrix is invertible, its inverse can also be used to obtain solutions of linear systems. This technique is recommended if you know in advance that the matrix A is not (close to being) singular:
// set up a linear system A * x = b
Matrix<3, 3> A { 3, -1, -4,
-1, 3, 2,
-1, -2, 1 };
Vec3 b { 2, 6, 1 };
// compute the inverse A^-1 (if A were singular, this would
// throw an exception)
auto AI = inverse(A);
// obtain the solution as x = A^-1 * b
std::cout << "A = " << A << '\n'
<< "b = " << b << '\n'
<< "A^-1 = " << AI << '\n'
<< "x = A^-1 * b = " << AI * b << '\n';
// solve can obtain the same solution, but is usually slower
std::cout << "x = " << solve(A, b) << '\n';If a linear system has multiple (infinitely many) solutions, solve will obtain the one with the smallest Euclidean norm.
// compute the point within the plane 3x + 5y - 2z = 19 that
// is closest to the origin (i.e. has smallest norm)
Matrix<1, 3> M { 3, 5, -2 };
Matrix<1, 1> b { 19 };
// the vector pointing to p will be orthogonal to the plane,
// i.e. parallel to its normal vector (3, 5, -2)
auto p = solve(M, b);
// verify that p actually lies within the plane
std::cout << "M = " << M << '\n'
<< "p = " << p << '\n'
<< "b = " << b << '\n'
<< "M * p = " << M * p << '\n';When presented with an inconsistent or overconstrained linear system AX=B without an exact solution, solve will find the X that minimizes the error AX-B in the least-squares sense:
// find the values a and b such that the line ax+by=1 fits
// the given data points as well as possible
Matrix<5, 2> data { 0, 3.056,
1, 2.424,
2, 2.004,
3, 1.425,
4, 0.946 };
Vector<5> rhs { 1,
1,
1,
1,
1 };
// compute the "best" solution of data * line = rhs
auto line = solve(data, rhs);
// inspect the solution and check its quality
std::cout << "data = " << data << '\n'
<< "line = " << line << '\n'
<< "error = " << data * line - rhs << '\n';
// convert line to form y=mx+c
std::cout << "m = " << -line[0] / line[1] << '\n'
<< "c = " << 1. / line[1] << '\n';The solve function can also be used to compute the pseudo-inverse of an arbitrary matrix:
// a non-invertible matrix (could even be rectangular)
Matrix<2, 2> M { 3, 4,
6, 8 };
// solving for the identity as the right-hand side computes
// the pseudo-inverse
auto P = solve(M, identity<2>());
// verify pseudo-inverse properties
std::cout << "M = " << M << '\n'
<< "P = " << P << '\n'
<< "P * M * P = " << P * M * P << '\n'
<< "M * P * M = " << M * P * M << '\n';Although they can be represented as matrices, rotations have a lot more structure that can be used to simplify certain operations and to rule out others. For example, inverting a rotation matrix is computationally trivial, while adding two rotation matrices or setting a single component of it individually does not make sense. Due to such semantic differences, mlina provides the dedicated type SO3 for rotations in three dimensions.
// construct a rotation from an axis of rotation (does not
// have to be normalized) and an angle (in radians)
SO3 R { Vec3 { 1, 3, -2 }, .25 * M_PI };
// rotations are printed just like matrices (but internally,
// they are represented more efficiently)
std::cout << "R = " << R << '\n';
// rotations can be applied to vectors (as if they were
// matrices, but internally, specialized algorithms are
// used)
Vec3 x { 1, 0, 0 };
Vec3 a { 1, 3, -2 };
std::cout << "x = " << x << '\n'
<< "a = " << a << '\n'
<< "R * x = " << R * x << '\n'
<< "R * a = " << R * a << '\n';
// for example, the norm is invariant under rotations
std::cout << "|a| = " << norm(a) << '\n'
<< "|R * a| = " << norm(R * a) << '\n';
// as another example, rotations distribute over the outer
// product
Vec3 p { 1, -1, 5 }, q { -2, 0, 3 };
std::cout << "p = " << p << '\n'
<< "q = " << q << '\n'
<< "Rp x Rq = " << cross(R * p, R * q) << '\n'
<< "R(p x q) = " << R * cross(p, q) << '\n';
// rotations can also be multiplied with rotations (but not
// with general matrices), which results in a new rotation
auto R2 = R * R;
// rotations can be queried for their axis and their angle
std::cout << "R^2 axis = " << R2.axis() << '\n'
<< "R^2 angle = " << R2.angle() / M_PI << "pi\n";
// a 360 degree rotation is the identity
std::cout << "R^8 = " << R2 * R2 * R2 * R2 << '\n';It is possible to convert rotations into plain matrices and vice versa:
// construct 90 degree rotation about the z-axis
SO3 R { Vec3 { 0, 0, 1 }, M_PI / 2 };
// convert the rotation into a plain matrix
Matrix<3, 3> M = R.matrix();
// computations on M are mathematically equivalent but
// usually less efficient
std::cout << "R = " << R << '\n'
<< "M = " << M << '\n'
<< "R * R = " << R * R << '\n'
<< "M * M = " << M * M << '\n'
<< "R^-1 = " << inverse(R) << '\n'
<< "M^-1 = " << inverse(M) << '\n';
// if (an approximation of) a rotation matrix is given, it
// can be converted into a rotation for further processing
SO3 R2 { M * M };
std::cout << "R^2 axis = " << R2.axis() << '\n'
<< "R^2 angle = " << R2.angle() / M_PI << "pi\n";Homogeneous transformations combine a rotation and a translation forming the group of all isometries of Euclidean space that preserve handedness. Homogeneous transformations are sometimes also referred to as "poses".
// a homogeneous transformation can be constructed from a
// rotation and a translation
SO3 R { Vec3 { 0, 0, 1 }, M_PI / 2 };
Vec3 t { 4, -1, 2 };
SE3 T { R, t };
// homogeneous transformations are printed in their 4x4
// representation as a homogeneous transformation matrix
std::cout << "R = " << R << '\n'
<< "t = " << t << '\n'
<< "T = " << T << '\n';
// when applied to a vector, T first performs the rotation
// followed by the translation, i.e. Tv = Rv + t
Vec3 v { 5, 7, -2 };
std::cout << "v = " << v << '\n'
<< "Tv = " << T * v << '\n'
<< "Rv + t = " << R * v + t << '\n';
// homogeneous transformations can be queried for their
// rotation (SO3) and their translation (Vec3) part
std::cout << "T rotation = " << T.rotation() << '\n'
<< "T translation = " << T.translation() << '\n';
// homogeneous transformations can be composed
auto T2 = T * T;
std::cout << "T^2 = " << T2 << '\n';
// the analogous computations in matrix form yield the same
// results, but are way less efficient as they cannot make
// use of the special structure of SE3
Matrix<4, 4> M = T.matrix();
std::cout << "M^2 = " << M * M << '\n';
// inverting homogeneous transformations is computationally
// cheap (while the analogous computation on their matrix
// representation is expensive)
std::cout << "T^-1 = " << inverse(T) << '\n'
<< "M^-1 = " << inverse(M) << '\n';mlina supports linear interpolation between points in space (Vector<n>), orientations (SO3) and also poses (SE3) by means of the Interpolation types. Instances can be constructed from an initial value A and a target value B. Such an interpolation is a functor representing "linear" interpolation between A and B. For example, f(0) will yield A, f(0.5) will result in the midpoint between A and B, and f(1) will yield B.
The most basic case is interpolation between points:
// two points in space
Vec3 p { 2, 5, 1 }, q { 6, 1, -1 };
std::cout << "p = " << p << '\n'
<< "q = " << q << '\n';
// linear interpolation between two points in space (works
// in any dimension)
Interpolation<Vec3> f { p, q };
// inspect some points along the line connecting p and q
std::cout << "f(0) = " << f(0) << '\n'
<< "f(0.25) = " << f(0.25) << '\n'
<< "f(0.5) = " << f(0.5) << '\n'
<< "f(1) = " << f(1) << '\n';
// an interpolation may be evaluated for progress values
// outside of [0,1]
std::cout << "f(-1) = " << f(-1) << '\n'
<< "f(2) = " << f(2) << '\n';
// as the progress value increases from 0 to 1, the distance
// to the target point decreases linearly
for (int i = 0; i <= 10; ++i) {
double x = 0.1 * i;
std::cout << "d(f(" << x << "), q) = "
<< dist(f(x), q) << '\n';
}It is conceptually clear that an orientation (SO3) can be smoothly changed into any given target orientation: visually, this corresponds to rotating a rigid body until the desired orientation is reached. However, it might not be obvious how exactly this can be done in the most direct way, or how such a trajectory through SO3 could be implemented. mlina can help:
// two orientations (rotations default to the identity)
SO3 R1;
SO3 R2 { Vec3 { 1, 1, 1 }, 0.6 * M_PI };
std::cout << "R1 = " << R1 << '\n'
<< "R2 = " << R2 << '\n';
// interpolate between orientations
Interpolation<SO3> f { R1, R2 };
// every "point" along f is a rotation (note that
// component-wise linear interpolation of the matrix
// representations would not make any sense)
std::cout << "f(0) = " << f(0) << '\n'
<< "f(0.5) = " << f(0.5) << '\n'
<< "f(1) = " << f(1) << '\n';
// there is also a distance measure for rotations (the
// result is the "difference angle" between 0 and pi)
std::cout << "d(R1, R2) = "
<< dist(R1, R2) / M_PI << "pi\n";
// the distance of two rotations is zero if and only if they
// are identical (dist defines a metric in the usual
// mathematical sense)
std::cout << "d(R2, R2) = " << dist(R2, R2) << '\n';
// as the progress value increases, the distance to the
// target orientation decreases linearly
for (int i = 0; i <= 10; ++i) {
double x = .1 * i;
std::cout << "d(f(" << x << "), R2) = "
<< dist(f(x), R2) / M_PI
<< "pi\n";
}For convenience, mlina also provides interpolation between poses (SE3). These are simply a combination of the interpolations between positions and orientations. Since there is no natural way of comparing distances of positions and distances of orientations, there is no distance measure for poses in mlina.
// two poses (e.g. of a rigid body)
SE3 T1 {
SO3 {},
Vec3 { 2, 5, 1 }
};
SE3 T2 {
SO3 { Vec3 { 1, 1, 1 }, 0.6 * M_PI },
Vec3 { 6, 1, -1 }
};
std::cout << "T1 = " << T1 << '\n'
<< "T2 = " << T2 << '\n';
// interpolation between poses
Interpolation<SE3> f { T1, T2 };
// inspect some poses along the way
std::cout << "f(0) = " << f(0) << '\n'
<< "f(0.5) = " << f(0.5) << '\n'
<< "f(1) = " << f(1) << '\n';