diff --git a/glm/gtx/matrix_factorisation.hpp b/glm/gtx/matrix_factorisation.hpp index dcdd6989..c553c121 100644 --- a/glm/gtx/matrix_factorisation.hpp +++ b/glm/gtx/matrix_factorisation.hpp @@ -27,7 +27,6 @@ /* Suggestions: - Move helper functions flipud and flip lr to another file: They may be helpful in more general circumstances. - - When rq_decompose is fed a matrix that has more rows than columns, the resulting r matrix is NOT upper triangular. Is that a bug? - Implement other types of matrix factorisation, such as: QL and LQ, L(D)U, eigendecompositions, etc... */ @@ -46,15 +45,16 @@ namespace glm{ GLM_FUNC_DECL matType fliplr(const matType& in); /// Performs QR factorisation of a matrix. - /// Returns 2 matrices, q and r, such that q columns are orthonormal, r is an upper triangular matrix, and q*r=in. - /// r is a square matrix whose dimensions are the same than the width of the input matrix, and q has the same dimensions than the input 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 class matType> GLM_FUNC_DECL void qr_decompose(matType& q, matType& r, const matType& in); /// Performs RQ factorisation of a matrix. - /// Returns 2 matrices, r and q, such that r is an upper triangular matrix, q rows are orthonormal, and r*q=in. - /// q has the same dimensions than the input matrix, and r is a square matrix whose dimensions are the same than the height of the input 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 class matType> GLM_FUNC_DECL void rq_decompose(matType& r, matType& q, const matType& in); diff --git a/glm/gtx/matrix_factorisation.inl b/glm/gtx/matrix_factorisation.inl index ec9ae0f1..f53d8280 100644 --- a/glm/gtx/matrix_factorisation.inl +++ b/glm/gtx/matrix_factorisation.inl @@ -29,21 +29,25 @@ namespace glm { // 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 < std::min(R, C); 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]; - } - - q[i] = normalize(q[i]); - } - - for (length_t i = 0; i < std::min(R, C); i++) { - for (length_t j = 0; j < i; 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]); } diff --git a/test/gtx/gtx_matrix_factorisation.cpp b/test/gtx/gtx_matrix_factorisation.cpp index 9ea46874..d45c352e 100644 --- a/test/gtx/gtx_matrix_factorisation.cpp +++ b/test/gtx/gtx_matrix_factorisation.cpp @@ -22,17 +22,17 @@ int test_qr(matType m) { //Test if the columns of q are orthonormal for (glm::length_t i = 0; i < std::min(C, R); i++) { - if ((length(q[i]) - 1) > epsilon) return 1; + if ((length(q[i]) - 1) > epsilon) return 2; for (glm::length_t j = 0; j epsilon) return 1; + if (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 < std::min(C, R); j++) { - if (r[i][j] != 0) return 1; + if (r[i][j] != 0) return 4; } } @@ -61,17 +61,17 @@ int test_rq(matType m) { matType tq = transpose(q); for (glm::length_t i = 0; i < std::min(C, R); i++) { - if ((length(tq[i]) - 1) > epsilon) return 1; + if ((length(tq[i]) - 1) > epsilon) return 2; for (glm::length_t j = 0; j epsilon) return 1; + if (abs(dot(tq[i], tq[j])) > epsilon) return 3; } } //Test if the matrix r is upper triangular for (glm::length_t i = 0; i < std::min(C, R); i++) { - for (glm::length_t j = i + 1; j < R; j++) { - if (r[i][j] != 0) return 1; + for (glm::length_t j = R - std::min(C, R) + i + 1; j < R; j++) { + if (r[i][j] != 0) return 4; } } @@ -85,19 +85,19 @@ int main() 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 1; + 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 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 1; + 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 1; + 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 1; + if (test_rq(glm::dmat4x3(12, 6, -4, -51, 167, 24, 4, -68, -41, 7, 2, 15))) return 6; return 0; }