Magnum::Math::Algorithms namespace

Algorithms.

Various matrix and vector algorithms.

This library is built as part of Magnum by default. To use this library with CMake, find the Magnum package and link to the Magnum::Magnum target:

find_package(Magnum REQUIRED)

# ...
target_link_libraries(your-app PRIVATE Magnum::Magnum)

See Downloading and building and Usage with CMake for more information.

Functions

template<std::size_t size, std::size_t rows, class T>
auto gaussJordanInPlaceTransposed(RectangularMatrix<size, size, T>& a, RectangularMatrix<size, rows, T>& t) -> bool
In-place Gauss-Jordan elimination on transposed matrices.
template<std::size_t size, std::size_t cols, class T>
auto gaussJordanInPlace(RectangularMatrix<size, size, T>& a, RectangularMatrix<cols, size, T>& t) -> bool
In-place Gauss-Jordan elimination.
template<std::size_t size, class T>
auto gaussJordanInverted(Matrix<size, T> matrix) -> Matrix<size, T>
Gauss-Jordan matrix inversion.
template<std::size_t cols, std::size_t rows, class T>
void gramSchmidtOrthogonalizeInPlace(RectangularMatrix<cols, rows, T>& matrix)
In-place Gram-Schmidt matrix orthogonalization.
template<std::size_t cols, std::size_t rows, class T>
auto gramSchmidtOrthogonalize(RectangularMatrix<cols, rows, T> matrix) -> RectangularMatrix<cols, rows, T>
Gram-Schmidt matrix orthogonalization.
template<std::size_t cols, std::size_t rows, class T>
void gramSchmidtOrthonormalizeInPlace(RectangularMatrix<cols, rows, T>& matrix)
In-place Gram-Schmidt matrix orthonormalization.
template<std::size_t cols, std::size_t rows, class T>
auto gramSchmidtOrthonormalize(RectangularMatrix<cols, rows, T> matrix) -> RectangularMatrix<cols, rows, T>
Gram-Schmidt matrix orthonormalization.
template<class Iterator, class T = typename std::decay<decltype(*std::declval<Iterator>())>::type>
auto kahanSum(Iterator begin, Iterator end, T sum = T(0), T* compensation = nullptr) -> T
Kahan summation algorithm.
template<std::size_t size, class T>
auto qr(const Matrix<size, T>& matrix) -> Containers::Pair<Matrix<size, T>, Matrix<size, T>>
QR decomposition.
template<std::size_t cols, std::size_t rows, class T>
auto svd(RectangularMatrix<cols, rows, T> m) -> Containers::Optional<Containers::Triple<RectangularMatrix<cols, rows, T>, Vector<cols, T>, Matrix<cols, T>>>
Singular Value Decomposition.

Function documentation

#include <Magnum/Math/Algorithms/GaussJordan.h>
template<std::size_t size, std::size_t rows, class T>
bool Magnum::Math::Algorithms::gaussJordanInPlaceTransposed(RectangularMatrix<size, size, T>& a, RectangularMatrix<size, rows, T>& t)

In-place Gauss-Jordan elimination on transposed matrices.

Parameters
a Transposed left side of augmented matrix
t Transposed right side of augmented matrix
Returns True if a is regular, false if a is singular (and thus the system cannot be solved).

As Gauss-Jordan elimination works on rows and matrices in OpenGL are column-major, it is more efficient to operate on transposed matrices and treat columns as rows. See also gaussJordanInPlace() which works with non-transposed matrices.

The function eliminates matrix a and solves t in place. For efficiency reasons, only pure Gaussian elimination is done on a and the final backsubstitution is done only on t, as a would always end with identity matrix anyway.

Based on an ultra-compact Python code by Jarno Elonen, http://elonen.iki.fi/code/misc-notes/python-gaussj/index.html.

#include <Magnum/Math/Algorithms/GaussJordan.h>
template<std::size_t size, std::size_t cols, class T>
bool Magnum::Math::Algorithms::gaussJordanInPlace(RectangularMatrix<size, size, T>& a, RectangularMatrix<cols, size, T>& t)

In-place Gauss-Jordan elimination.

Transposes the matrices, calls gaussJordanInPlaceTransposed() on them and then transposes them back.

#include <Magnum/Math/Algorithms/GaussJordan.h>
template<std::size_t size, class T>
Matrix<size, T> Magnum::Math::Algorithms::gaussJordanInverted(Matrix<size, T> matrix)

Gauss-Jordan matrix inversion.

Uses the Gauss-Jordan elimination to perform a matrix inversion. Since $ (\boldsymbol{A}^{-1})^T = (\boldsymbol{A}^T)^{-1} $ , passes matrix and an identity matrix to gaussJordanInPlaceTransposed(); returning the inverted matrix. Expects that the matrix is invertible.

#include <Magnum/Math/Algorithms/GramSchmidt.h>
template<std::size_t cols, std::size_t rows, class T>
void Magnum::Math::Algorithms::gramSchmidtOrthogonalizeInPlace(RectangularMatrix<cols, rows, T>& matrix)

In-place Gram-Schmidt matrix orthogonalization.

Parameters
matrix in/out Matrix to perform orthogonalization on

Performs the Gram-Schmidt process. With a projection operator defined as

\[ \operatorname{proj}_{\boldsymbol{u}}\,(\boldsymbol{v}) = \frac{\boldsymbol{u} \cdot \boldsymbol{v}}{\boldsymbol{u} \cdot \boldsymbol{u}} \boldsymbol{u} \]

the process works as follows, with $ \boldsymbol{v}_k $ being columns of matrix and $ \boldsymbol{u}_k $ columns of the output:

\[ \boldsymbol{u}_k = \boldsymbol{v}_k - \sum_{j = 1}^{k - 1} \operatorname{proj}_{\boldsymbol{u}_j}\,(\boldsymbol{v}_k) \]

Note that the above is not performed directly due to numerical instability, the stable modified Gram-Schmidt algorithm is used instead.

#include <Magnum/Math/Algorithms/GramSchmidt.h>
template<std::size_t cols, std::size_t rows, class T>
RectangularMatrix<cols, rows, T> Magnum::Math::Algorithms::gramSchmidtOrthogonalize(RectangularMatrix<cols, rows, T> matrix)

Gram-Schmidt matrix orthogonalization.

Unlike gramSchmidtOrthogonalizeInPlace() returns the modified matrix instead of performing the orthogonalization in-place.

#include <Magnum/Math/Algorithms/GramSchmidt.h>
template<std::size_t cols, std::size_t rows, class T>
void Magnum::Math::Algorithms::gramSchmidtOrthonormalizeInPlace(RectangularMatrix<cols, rows, T>& matrix)

In-place Gram-Schmidt matrix orthonormalization.

Parameters
matrix in/out Matrix to perform orthonormalization on

Performs the Gram-Schmidt process. With a projection operator defined as

\[ \operatorname{proj}_{\boldsymbol{u}}\,(\boldsymbol{v}) = \frac{\boldsymbol{u} \cdot \boldsymbol{v}}{\boldsymbol{u} \cdot \boldsymbol{u}} \boldsymbol{u} \]

the Gram-Schmidt process works as follows, with $ \boldsymbol{v}_k $ being columns of matrix and $ \boldsymbol{e}_k $ columns of the output:

\[ \boldsymbol{u}_k = \boldsymbol{v}_k - \sum_{j = 1}^{k - 1} \operatorname{proj}_{\boldsymbol{u}_j}\,(\boldsymbol{v}_k), ~~~~~~ \boldsymbol{e}_k = \frac{\boldsymbol{u}_k}{|\boldsymbol{u}_k|} \]

In particular, this adds the normalization step on top of gramSchmidtOrthogonalizeInPlace(). Note that the above is not performed directly due to numerical instability, the stable modified Gram-Schmidt algorithm is used instead.

#include <Magnum/Math/Algorithms/GramSchmidt.h>
template<std::size_t cols, std::size_t rows, class T>
RectangularMatrix<cols, rows, T> Magnum::Math::Algorithms::gramSchmidtOrthonormalize(RectangularMatrix<cols, rows, T> matrix)

Gram-Schmidt matrix orthonormalization.

Unlike gramSchmidtOrthonormalizeInPlace() returns the modified matrix instead of performing the orthonormalization in-place.

#include <Magnum/Math/Algorithms/KahanSum.h>
template<class Iterator, class T = typename std::decay<decltype(*std::declval<Iterator>())>::type>
T Magnum::Math::Algorithms::kahanSum(Iterator begin, Iterator end, T sum = T(0), T* compensation = nullptr)

Kahan summation algorithm.

Parameters
begin in Range begin
end in Range end
sum in Initial value for the sum
compensation in/out Floating-point roundoff error compensation value. If non- nullptr, value behind the pointer is used as initial compensation value and the resulting compensation value is saved back there.

Calculates a sum of a large range of floating-point numbers with roundoff error compensation. Compared to for example std::accumulate() the algorithm significantly reduces numerical error in the total. See the Kahan summation algorithm article on Wikipedia for an in-depth explanation. Example with summation of a hundred million ones:

std::vector<Float> data(100000000, 1.0f);
Float a = std::accumulate(data.begin(), data.end(), 0.0f);      // 1.667e7f
Float b = Math::Algorithms::kahanSum(data.begin(), data.end()); // 1.000e8f

If required, it is also possible to use this algorithm on non-contiguous ranges or single values (for example when calculating sum of pixel values in an image with some row padding or when the inputs are generated / converted from other values):

Containers::ArrayView<UnsignedByte> pixels;
Float sum = 0.0f, c = 0.0f;
for(UnsignedByte pixel: pixels) {
    Float value = Math::unpack<Float>(pixel);
    sum = Math::Algorithms::kahanSum(&value, &value + 1, sum, &c);
}

#include <Magnum/Math/Algorithms/Qr.h>
template<std::size_t size, class T>
Containers::Pair<Matrix<size, T>, Matrix<size, T>> Magnum::Math::Algorithms::qr(const Matrix<size, T>& matrix)

QR decomposition.

Calculated using the Gram-Schmidt process, in particular the modifed Gram-Schmidt from gramSchmidtOrthonormalize(). Given the input matrix $ \boldsymbol{A} = (\boldsymbol{a}_1, \cdots, \boldsymbol{a}_n) $ and the set of column vectors $ \boldsymbol{e}_k $ coming from the Gram-Schmidt process, the resulting $ \boldsymbol{Q} $ and $ \boldsymbol{R} $ matrices are as follows:

\[ \begin{array}{rcl} \boldsymbol{Q} & = & \left( \boldsymbol{e}_1, \cdots, \boldsymbol{e}_n \right) \\[10pt] \boldsymbol{R} & = & \begin{pmatrix} \boldsymbol{e}_1 \cdot \boldsymbol{a}_1 & \boldsymbol{e}_1 \cdot \boldsymbol{a}_2 & \boldsymbol{e}_1 \cdot \boldsymbol{a}_3 & \ldots \\ 0 & \boldsymbol{e}_2 \cdot \boldsymbol{a}_2 & \boldsymbol{e}_2 \cdot \boldsymbol{a}_3 & \ldots \\ 0 & 0 & \boldsymbol{e}_3 \cdot \boldsymbol{a}_3 & \ldots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix} \end{array} \]

One possible use is to decompose a transformation matrix into separate rotation and scaling/shear parts. Note, however, that the decomposition is not unique. See the associated test case for an example.

#include <Magnum/Math/Algorithms/Svd.h>
template<std::size_t cols, std::size_t rows, class T>
Containers::Optional<Containers::Triple<RectangularMatrix<cols, rows, T>, Vector<cols, T>, Matrix<cols, T>>> Magnum::Math::Algorithms::svd(RectangularMatrix<cols, rows, T> m)

Singular Value Decomposition.

Performs Thin SVD on given matrix where rows >= cols:

\[ \boldsymbol{M} = \boldsymbol{U} \boldsymbol{\Sigma} \boldsymbol{V}^* \]

Returns first cols column vectors of $ \boldsymbol{U} $ , diagonal of $ \boldsymbol{\Sigma} $ and non-transposed $ \boldsymbol{V} $ . If the solution doesn't converge, returns Containers::NullOpt.

Full $ \boldsymbol{U} $ , $ \boldsymbol{\Sigma} $ matrices and original $ \boldsymbol{M} $ matrix can be reconstructed from the values as following:

Math::RectangularMatrix<cols, rows, Double> m;

auto svd = Math::Algorithms::svd(m);
Math::RectangularMatrix<cols, rows, Double> uPart = svd->first();
Math::Vector<cols, Double> wDiagonal = svd->second();

// Extend U
Math::Matrix<rows, Double> u{Math::ZeroInit};
for(std::size_t i = 0; i != rows; ++i)
    u[i] = uPart[i];

// Diagonal W
Math::RectangularMatrix<cols, rows, Double> w =
    Math::RectangularMatrix<cols, rows, Double>::fromDiagonal(wDiagonal);

// u*w*v.transposed() == m

One possible use is to decompose a transformation matrix into separate rotation and scaling parts. Note, however, that the decomposition is not unique. See the associated test case for an example. Implementation based on Golub, G. H.; Reinsch, C. (1970). "Singular value decomposition and least squares solutions".