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';