diff --git a/src/Magnum/Math/Algorithms/Qr.h b/src/Magnum/Math/Algorithms/Qr.h
index f6c68f268..9531d810c 100644
--- a/src/Magnum/Math/Algorithms/Qr.h
+++ b/src/Magnum/Math/Algorithms/Qr.h
@@ -29,26 +29,44 @@
  * @brief Function @ref Magnum::Math::Algorithms::qr()
  */
 
-#include "Magnum/Math/Matrix.h"
+#include "Magnum/Math/Algorithms/GramSchmidt.h"
 
 namespace Magnum { namespace Math { namespace Algorithms {
 
 /**
 @brief QR decomposition
 
-Calculated using the classic [Gram-Schmidt process](https://en.wikipedia.org/wiki/QR_decomposition#Using_the_Gram–Schmidt_process).
+Calculated using the [Gram-Schmidt process](https://en.wikipedia.org/wiki/QR_decomposition#Using_the_Gram–Schmidt_process),
+in particular the [modifed Gram-Schmidt](https://en.wikipedia.org/wiki/Gram–Schmidt_process#Numerical_stability)
+from @ref gramSchmidtOrthonormalize(). Given the input matrix
+@f$ \boldsymbol{A} = (\boldsymbol{a}_1, \cdots, \boldsymbol{a}_n) @f$ and the set of
+column vectors @f$ \boldsymbol{e}_k @f$ coming from the Gram-Schmidt process,
+the resulting @f$ \boldsymbol{Q} @f$ and @f$ \boldsymbol{R} @f$ matrices are as
+follows: @f[
+    \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}
+@f]
+
+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](https://github.com/mosra/magnum/blob/master/src/Magnum/Math/Algorithms/Test/QrTest.cpp)
+for an example.
+@see @ref svd(), @ref Matrix3::rotationShear(), @ref Matrix4::rotationShear()
 */
 template<std::size_t size, class T> std::pair<Matrix<size, T>, Matrix<size, T>> qr(const Matrix<size, T>& matrix) {
-    Matrix<size, T> q, r;
-
+    const Matrix<size, T> q = gramSchmidtOrthonormalize(matrix);
+    Matrix<size, T> r{ZeroInit};
     for(std::size_t k = 0; k != size; ++k) {
-        Vector<size, T> p = matrix[k];
-        for(std::size_t j = 0; j != k; ++j) {
-            r[k][j] = Math::dot(p, q[j]);
-            p -= q[j]*r[k][j];
+        for(std::size_t j = 0; j <= k; ++j) {
+            r[k][j] = Math::dot(q[j], matrix[k]);
         }
-        r[k][k] = p.length();
-        q[k] = p/r[k][k];
     }
 
     return {q, r};
diff --git a/src/Magnum/Math/Algorithms/Test/QrTest.cpp b/src/Magnum/Math/Algorithms/Test/QrTest.cpp
index 9a28b7326..33a47b92e 100644
--- a/src/Magnum/Math/Algorithms/Test/QrTest.cpp
+++ b/src/Magnum/Math/Algorithms/Test/QrTest.cpp
@@ -27,6 +27,7 @@
 #include <Corrade/TestSuite/Tester.h>
 
 #include "Magnum/Math/Matrix3.h"
+#include "Magnum/Math/Matrix4.h"
 #include "Magnum/Math/Algorithms/Qr.h"
 
 namespace Magnum { namespace Math { namespace Algorithms { namespace Test {
@@ -35,35 +36,54 @@ struct QrTest: Corrade::TestSuite::Tester {
     explicit QrTest();
 
     void test();
+    void decomposeRotationShear();
 };
 
-typedef Matrix3<Float> Matrix3;
+using namespace Math::Literals;
+
+typedef Matrix3x3<Float> Matrix3x3;
+typedef Vector3<Float> Vector3;
+typedef Matrix4<Float> Matrix4;
 
 QrTest::QrTest() {
-    addTests({&QrTest::test});
+    addTests({&QrTest::test,
+              &QrTest::decomposeRotationShear});
 }
 
 void QrTest::test() {
-    Matrix3 a{{  0.0f,   3.0f,   4.0f},
-              {-20.0f,  27.0f,  11.0f},
-              {-14.0f,  -4.0f,  -2.0f}};
+    Matrix3x3 a{Vector3{  0.0f,   3.0f,   4.0f},
+                Vector3{-20.0f,  27.0f,  11.0f},
+                Vector3{-14.0f,  -4.0f,  -2.0f}};
 
-    Matrix3 q, r;
+    Matrix3x3 q, r;
     std::tie(q, r) = Algorithms::qr(a);
 
-    CORRADE_COMPARE(q*r, a);
-
-    Matrix3 qExpected = Matrix3{{  0.0f,  15.0f, 20.0f},
-                                {-20.0f,  12.0f, -9.0f},
-                                {-15.0f, -16.0f, 12.0f}}/25.0f;
+    auto qExpected = Matrix3x3{Vector3{  0.0f,  15.0f, 20.0f},
+                               Vector3{-20.0f,  12.0f, -9.0f},
+                               Vector3{-15.0f, -16.0f, 12.0f}}/25.0f;
     CORRADE_COMPARE(q, qExpected);
 
-    Matrix3 rExpected{{ 5.0f,  0.0f,  0.0f},
-                      {25.0f, 25.0f,  0.0f},
-                      {-4.0f, 10.0f, 10.0f}};
+    Matrix3x3 rExpected{Vector3{ 5.0f,  0.0f,  0.0f},
+                        Vector3{25.0f, 25.0f,  0.0f},
+                        Vector3{-4.0f, 10.0f, 10.0f}};
     CORRADE_COMPARE(r, rExpected);
 }
 
+void QrTest::decomposeRotationShear() {
+    Matrix4 a = Matrix4::scaling({1.5f, 2.0f, 1.0f})*Matrix4::rotationZ(35.0_degf);
+
+    Matrix3x3 q, r;
+    std::tie(q, r) = Algorithms::qr(a.rotationScaling());
+    CORRADE_COMPARE(q*r, a.rotationScaling());
+
+    auto q4 = Matrix4::from(q, {});
+    auto r4 = Matrix4::from(r, {});
+
+    CORRADE_COMPARE(q4, Matrix4::rotationZ(43.03357_degf));
+    CORRADE_COMPARE(r4.scaling(), (Vector3{1.68099f, 1.85048f, 1.0f}));
+    CORRADE_COMPARE(r4.rotationShear(), Matrix4::shearingXZ(0.274077f, 0.0f).rotationShear());
+}
+
 }}}}
 
 CORRADE_TEST_MAIN(Magnum::Math::Algorithms::Test::QrTest)