Merge pull request #654 from unordinateur/master

QR factorisation of matrices #654

.gitignore

@@ -52,3 +52,5 @@ Makefile
 # local build(s)
 build*
 
+/.vs
+/CMakeSettings.json

glm/gtx/matrix_factorisation.hpp

@@ -0,0 +1,64 @@
+/// @ref gtx_matrix_factorisation
+/// @file glm/gtx/matrix_factorisation.hpp
+///
+/// @see core (dependence)
+///
+/// @defgroup gtx_matrix_factorisation GLM_GTX_matrix_factorisation
+/// @ingroup gtx
+///
+/// @brief Functions to factor matrices in various forms
+///
+/// <glm/gtx/matrix_factorisation.hpp> need to be included to use these functionalities.
+
+#pragma once
+
+// Dependency:
+#include "../glm.hpp"
+
+#ifndef GLM_ENABLE_EXPERIMENTAL
+#	error "GLM: GLM_GTX_matrix_factorisation is an experimental extension and may change in the future. Use #define GLM_ENABLE_EXPERIMENTAL before including it, if you really want to use it."
+#endif
+
+#if GLM_MESSAGES == GLM_MESSAGES_ENABLED && !defined(GLM_EXT_INCLUDED)
+#	pragma message("GLM: GLM_GTX_matrix_factorisation extension included")
+#endif
+
+/*
+Suggestions:
+ - Move helper functions flipud and fliplr to another file: They may be helpful in more general circumstances.
+ - Implement other types of matrix factorisation, such as: QL and LQ, L(D)U, eigendecompositions, etc...
+*/
+
+namespace glm{
+	/// @addtogroup gtx_matrix_factorisation
+	/// @{
+
+	/// Flips the matrix rows up and down.
+	/// From GLM_GTX_matrix_factorisation extension.
+	template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType>
+	GLM_FUNC_DECL matType<C, R, T, P> flipud(const matType<C, R, T, P>& in);
+
+	/// Flips the matrix columns right and left.
+	/// From GLM_GTX_matrix_factorisation extension.
+	template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType>
+	GLM_FUNC_DECL matType<C, R, T, P> fliplr(const matType<C, R, T, P>& in);
+
+	/// Performs QR factorisation of a matrix.
+	/// Returns 2 matrices, q and r, such that the columns of q are orthonormal and span the same subspace than those of the input matrix, r is an upper triangular matrix, and q*r=in.
+	/// Given an n-by-m input matrix, q has dimensions min(n,m)-by-m, and r has dimensions n-by-min(n,m).
+	/// From GLM_GTX_matrix_factorisation extension.
+	template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType>
+	GLM_FUNC_DECL void qr_decompose(matType<(C < R ? C : R), R, T, P>& q, matType<C, (C < R ? C : R), T, P>& r, const matType<C, R, T, P>& in);
+
+	/// Performs RQ factorisation of a matrix.
+	/// Returns 2 matrices, r and q, such that r is an upper triangular matrix, the rows of q are orthonormal and span the same subspace than those of the input matrix, and r*q=in.
+	/// Note that in the context of RQ factorisation, the diagonal is seen as starting in the lower-right corner of the matrix, instead of the usual upper-left.
+	/// Given an n-by-m input matrix, r has dimensions min(n,m)-by-m, and q has dimensions n-by-min(n,m).
+	/// From GLM_GTX_matrix_factorisation extension.
+	template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType>
+	GLM_FUNC_DECL void rq_decompose(matType<(C < R ? C : R), R, T, P>& r, matType<C, (C < R ? C : R), T, P>& q, const matType<C, R, T, P>& in);
+
+	/// @}
+}
+
+#include "matrix_factorisation.inl"

glm/gtx/matrix_factorisation.inl

@@ -0,0 +1,76 @@
+/// @ref gtx_matrix_factorisation
+/// @file glm/gtx/matrix_factorisation.inl
+
+namespace glm {
+	template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType>
+	GLM_FUNC_QUALIFIER matType<C, R, T, P> flipud(const matType<C, R, T, P>& in) {
+		matType<R, C, T, P> tin = transpose(in);
+		tin = fliplr(tin);
+		matType<C, R, T, P> out = transpose(tin);
+
+		return out;
+	}
+
+	template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType>
+	GLM_FUNC_QUALIFIER matType<C, R, T, P> fliplr(const matType<C, R, T, P>& in) {
+		matType<C, R, T, P> out;
+		for (length_t i = 0; i < C; i++) {
+			out[i] = in[(C - i) - 1];
+		}
+
+		return out;
+	}
+
+	template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType>
+	GLM_FUNC_QUALIFIER void qr_decompose(matType<(C < R ? C : R), R, T, P>& q, matType<C, (C < R ? C : R), T, P>& r, const matType<C, R, T, P>& in) {
+		// Uses modified Gram-Schmidt method
+		// Source: https://en.wikipedia.org/wiki/Gram–Schmidt_process
+		// And https://en.wikipedia.org/wiki/QR_decomposition
+
+		//For all the linearly independs columns of the input...
+		// (there can be no more linearly independents columns than there are rows.)
+		for (length_t i = 0; i < (C < R ? C : R); i++) {
+			//Copy in Q the input's i-th column.
+			q[i] = in[i];
+
+			//j = [0,i[
+			// Make that column orthogonal to all the previous ones by substracting to it the non-orthogonal projection of all the previous columns.
+			// Also: Fill the zero elements of R
+			for (length_t j = 0; j < i; j++) {
+				q[i] -= dot(q[i], q[j])*q[j];
+				r[j][i] = 0;
+			}
+
+			//Now, Q i-th column is orthogonal to all the previous columns. Normalize it.
+			q[i] = normalize(q[i]);
+
+			//j = [i,C[
+			//Finally, compute the corresponding coefficients of R by computing the projection of the resulting column on the other columns of the input.
+			for (length_t j = i; j < C; j++) {
+				r[j][i] = dot(in[j], q[i]);
+			}
+		}
+	}
+
+	template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType>
+	GLM_FUNC_QUALIFIER void rq_decompose(matType<(C < R ? C : R), R, T, P>& r, matType<C, (C < R ? C : R), T, P>& q, const matType<C, R, T, P>& in) {
+		// From https://en.wikipedia.org/wiki/QR_decomposition:
+		// The RQ decomposition transforms a matrix A into the product of an upper triangular matrix R (also known as right-triangular) and an orthogonal matrix Q. The only difference from QR decomposition is the order of these matrices.
+		// QR decomposition is Gram–Schmidt orthogonalization of columns of A, started from the first column.
+		// RQ decomposition is Gram–Schmidt orthogonalization of rows of A, started from the last row.
+
+		matType<R, C, T, P> tin = transpose(in);
+		tin = fliplr(tin);
+
+		matType<R, (C < R ? C : R), T, P> tr;
+		matType<(C < R ? C : R), C, T, P> tq;
+		qr_decompose(tq, tr, tin);
+
+		tr = fliplr(tr);
+		r = transpose(tr);
+		r = fliplr(r);
+
+		tq = fliplr(tq);
+		q = transpose(tq);
+	}
+} //namespace glm

test/gtx/CMakeLists.txt

@@ -21,6 +21,7 @@ glmCreateTestGTC(gtx_io)
 glmCreateTestGTC(gtx_log_base)
 glmCreateTestGTC(gtx_matrix_cross_product)
 glmCreateTestGTC(gtx_matrix_decompose)
+glmCreateTestGTC(gtx_matrix_factorisation)
 glmCreateTestGTC(gtx_matrix_interpolation)
 glmCreateTestGTC(gtx_matrix_major_storage)
 glmCreateTestGTC(gtx_matrix_operation)

test/gtx/gtx_matrix_factorisation.cpp

@@ -0,0 +1,103 @@
+#define GLM_ENABLE_EXPERIMENTAL
+#include <glm/gtx/matrix_factorisation.hpp>
+
+const double epsilon = 1e-10f;
+
+template <glm::length_t C, glm::length_t R, typename T, glm::precision P, template<glm::length_t, glm::length_t, typename, glm::precision> class matType>
+int test_qr(matType<C, R, T, P> m) {
+	matType<(C < R ? C : R), R, T, P> q(-999);
+	matType<C, (C < R ? C : R), T, P> r(-999);
+
+	glm::qr_decompose(q, r, m);
+
+	//Test if q*r really equals the input matrix
+	matType<C, R, T, P> tm = q*r;
+	matType<C, R, T, P> err = tm - m;
+
+	for (glm::length_t i = 0; i < C; i++) {
+		for (glm::length_t j = 0; j < R; j++) {
+			if (std::abs(err[i][j]) > epsilon) return 1;
+		}
+	}
+
+	//Test if the columns of q are orthonormal
+	for (glm::length_t i = 0; i < (C < R ? C : R); i++) {
+		if ((length(q[i]) - 1) > epsilon) return 2;
+
+		for (glm::length_t j = 0; j<i; j++) {
+			if (std::abs(dot(q[i], q[j])) > epsilon) return 3;
+		}
+	}
+
+	//Test if the matrix r is upper triangular
+	for (glm::length_t i = 0; i < C; i++) {
+		for (glm::length_t j = i + 1; j < (C < R ? C : R); j++) {
+			if (r[i][j] != 0) return 4;
+		}
+	}
+
+	return 0;
+}
+
+template <glm::length_t C, glm::length_t R, typename T, glm::precision P, template<glm::length_t, glm::length_t, typename, glm::precision> class matType>
+int test_rq(matType<C, R, T, P> m) {
+	matType<C, (C < R ? C : R), T, P> q(-999);
+	matType<(C < R ? C : R), R, T, P> r(-999);
+
+	glm::rq_decompose(r, q, m);
+
+	//Test if q*r really equals the input matrix
+	matType<C, R, T, P> tm = r*q;
+	matType<C, R, T, P> err = tm - m;
+
+	for (glm::length_t i = 0; i < C; i++) {
+		for (glm::length_t j = 0; j < R; j++) {
+			if (std::abs(err[i][j]) > epsilon) return 1;
+		}
+	}
+
+
+	//Test if the rows of q are orthonormal
+	matType<(C < R ? C : R), C, T, P> tq = transpose(q);
+
+	for (glm::length_t i = 0; i < (C < R ? C : R); i++) {
+		if ((length(tq[i]) - 1) > epsilon) return 2;
+
+		for (glm::length_t j = 0; j<i; j++) {
+			if (std::abs(dot(tq[i], tq[j])) > epsilon) return 3;
+		}
+	}
+
+	//Test if the matrix r is upper triangular
+	for (glm::length_t i = 0; i < (C < R ? C : R); i++) {
+		for (glm::length_t j = R - (C < R ? C : R) + i + 1; j < R; j++) {
+			if (r[i][j] != 0) return 4;
+		}
+	}
+
+	return 0;
+}
+
+int main()
+{
+
+	//Test QR square
+	if(test_qr(glm::dmat3(12, 6, -4, -51, 167, 24, 4, -68, -41))) return 1;
+
+	//Test RQ square
+	if (test_rq(glm::dmat3(12, 6, -4, -51, 167, 24, 4, -68, -41))) return 2;
+
+	//Test QR triangular 1
+	if (test_qr(glm::dmat3x4(12, 6, -4, -51, 167, 24, 4, -68, -41, 7, 2, 15))) return 3;
+
+	//Test QR triangular 2
+	if (test_qr(glm::dmat4x3(12, 6, -4, -51, 167, 24, 4, -68, -41, 7, 2, 15))) return 4;
+
+	//Test RQ triangular 1 : Fails at the triangular test
+	if (test_rq(glm::dmat3x4(12, 6, -4, -51, 167, 24, 4, -68, -41, 7, 2, 15))) return 5;
+
+	//Test QR triangular 2
+	if (test_rq(glm::dmat4x3(12, 6, -4, -51, 167, 24, 4, -68, -41, 7, 2, 15))) return 6;
+
+	return 0;
+}