Made my own equally wrong QR factorization

.vscode/settings.json

@@ -76,5 +76,8 @@
     "clangd.enable": true,
     "C_Cpp.dimInactiveRegions": false,
     "editor.defaultFormatter": "xaver.clang-format",
-    "clangd.inactiveRegions.useBackgroundHighlight": true
+    "clangd.inactiveRegions.useBackgroundHighlight": true,
+    "clangd.arguments": [
+        "--compile-commands-dir=${workspaceFolder}/build"
+    ],
 }

src/Matrix.cpp

@@ -451,16 +451,6 @@ float Matrix<rows, columns>::EuclideanNorm() const {
 
   return sqrt(sum);
 }
-template <uint8_t rows, uint8_t columns>
-float Matrix<rows, columns>::L2Norm() const {
-  float sum{0};
-  Matrix<rows, 1> columnMatrix{};
-  for (uint8_t column = 0; column < columns; column++) {
-    this.GetColumn(column, columnMatrix);
-    sum += columnMatrix.EuclideanNorm();
-  }
-  return sqrt(sum);
-}
 
 template <uint8_t rows, uint8_t columns>
 template <uint8_t sub_rows, uint8_t sub_columns, uint8_t row_offset,
@@ -485,20 +475,37 @@ Matrix<sub_rows, sub_columns> Matrix<rows, columns>::SubMatrix() const {
 }
 
 template <uint8_t rows, uint8_t columns>
-template <uint8_t sub_rows, uint8_t sub_columns, uint8_t row_offset,
+template <uint8_t sub_rows, uint8_t sub_columns>
-          uint8_t column_offset>
 void Matrix<rows, columns>::SetSubMatrix(
+    uint8_t rowOffset, uint8_t columnOffset,
     const Matrix<sub_rows, sub_columns> &sub_matrix) {
-  static_assert(sub_rows + row_offset <= rows,
+  int16_t adjustedSubRows = sub_rows;
-                "The submatrix you're trying to set is out of bounds (rows)");
+  int16_t adjustedSubColumns = sub_columns;
-  static_assert(
+  int16_t adjustedRowOffset = rowOffset;
-      sub_columns + column_offset <= columns,
+  int16_t adjustedColumnOffset = columnOffset;
-      "The submatrix you're trying to set is out of bounds (columns)");
 
-  for (uint8_t row_idx{0}; row_idx < sub_rows; row_idx++) {
+  // a bunch of safety checks to make sure we don't overflow the matrix
-    for (uint8_t column_idx{0}; column_idx < sub_columns; column_idx++) {
+  if (sub_rows > rows) {
-      this->matrix[(row_idx + row_offset) * columns + column_idx +
+    adjustedSubRows = rows;
-                   column_offset] = sub_matrix.Get(row_idx, column_idx);
+  }
+  if (sub_columns > columns) {
+    adjustedSubColumns = columns;
+  }
+
+  if (adjustedSubRows + adjustedRowOffset >= rows) {
+    adjustedRowOffset =
+        std::max(0, static_cast<int16_t>(rows) - adjustedSubRows);
+  }
+
+  if (adjustedSubColumns + adjustedColumnOffset >= columns) {
+    adjustedColumnOffset =
+        std::max(0, static_cast<int16_t>(columns) - adjustedSubColumns);
+  }
+
+  for (uint8_t row_idx{0}; row_idx < adjustedSubRows; row_idx++) {
+    for (uint8_t column_idx{0}; column_idx < adjustedSubColumns; column_idx++) {
+      this->matrix[(row_idx + adjustedRowOffset) * columns + column_idx +
+                   adjustedColumnOffset] = sub_matrix.Get(row_idx, column_idx);
     }
   }
 }
@@ -510,33 +517,37 @@ void Matrix<rows, columns>::QRDecomposition(Matrix<rows, columns> &Q,
                                             Matrix<columns, columns> &R) const {
   static_assert(columns <= rows, "QR decomposition requires columns <= rows");
 
-  Matrix<rows, 1> a_col, u, q_col, proj;
   Q.Fill(0);
   R.Fill(0);
+  Matrix<rows, 1> a_col, e, u, Q_column_k{};
+  Matrix<1, rows> a_T, e_T{};
 
-  for (uint8_t k = 0; k < columns; ++k) {
+  for (uint8_t column = 0; column < columns; column++) {
-    this->GetColumn(k, a_col);
+    this->GetColumn(column, a_col);
     u = a_col;
-
+    // -----------------------
-    for (uint8_t j = 0; j < k; ++j) {
+    // ----- CALCULATE Q -----
-      Q.GetColumn(j, q_col);
+    // -----------------------
-      float r_jk = Matrix<rows, 1>::DotProduct(
+    for (uint8_t k = 0; k <= column; k++) {
-          q_col, u); // FIXED: use u instead of a_col
+      Q.GetColumn(k, Q_column_k);
-      R[j][k] = r_jk;
+      Matrix<1, rows> Q_column_k_T = Q_column_k.Transpose();
-
+      u = u - Q_column_k * (Q_column_k_T * a_col);
-      proj = q_col * r_jk;
-      u = u - proj;
     }
-
+    float norm = u.EuclideanNorm();
-    float norm = sqrt(Matrix<rows, 1>::DotProduct(u, u));
+    if (norm > 1e-4) {
-    if (norm < 1e-12f)
+      u = u / norm;
-      norm = 1e-12f; // for stability
+    } else {
-
+      u.Fill(0);
-    for (uint8_t i = 0; i < rows; ++i) {
-      Q[i][k] = u[i][0] / norm;
     }
+    Q.SetSubMatrix(0, column, u);
 
-    R[k][k] = norm;
+    // -----------------------
+    // ----- CALCULATE R -----
+    // -----------------------
+    for (uint8_t k = 0; k <= column; k++) {
+      Q.GetColumn(k, e);
+      R[k][column] = (a_col.Transpose() * e).Get(0, 0);
+    }
   }
 }
 

src/Matrix.hpp

@@ -129,13 +129,12 @@ public:
   Matrix<columns, rows> Transpose() const;
 
   /**
-   * @brief reduce the matrix so the sum of its elements equal 1
+   * @brief Returns the euclidean magnitude of the matrix. Also known as the L2
+   * norm
    * @param result a buffer to store the result into
    */
   float EuclideanNorm() const;
 
-  float L2Norm() const;
-
   /**
    * @brief Get a row from the matrix
    * @param row_index the row index to get
@@ -209,9 +208,9 @@ public:
             uint8_t column_offset>
   Matrix<sub_rows, sub_columns> SubMatrix() const;
 
-  template <uint8_t sub_rows, uint8_t sub_columns, uint8_t row_offset,
+  template <uint8_t sub_rows, uint8_t sub_columns>
-            uint8_t column_offset>
+  void SetSubMatrix(uint8_t rowOffset, uint8_t columnOffset,
-  void SetSubMatrix(const Matrix<sub_rows, sub_columns> &sub_matrix);
+                    const Matrix<sub_rows, sub_columns> &sub_matrix);
 
   /**
    * @brief take the dot product of the two vectors

unit-tests/matrix-tests.cpp

@@ -336,27 +336,27 @@ TEST_CASE("Elementary Matrix Operations", "Matrix") {
     Matrix<3, 3> mat4 = startMatrix;
 
     Matrix<2, 2> mat5{10, 11, 12, 13};
-    mat4.SetSubMatrix<2, 2, 0, 0>(mat5);
+    mat4.SetSubMatrix(0, 0, mat5);
     REQUIRE(mat4.Get(0, 0) == 10);
     REQUIRE(mat4.Get(0, 1) == 11);
     REQUIRE(mat4.Get(1, 0) == 12);
     REQUIRE(mat4.Get(1, 1) == 13);
 
     mat4 = startMatrix;
-    mat4.SetSubMatrix<2, 2, 1, 1>(mat5);
+    mat4.SetSubMatrix(1, 1, mat5);
     REQUIRE(mat4.Get(1, 1) == 10);
     REQUIRE(mat4.Get(1, 2) == 11);
     REQUIRE(mat4.Get(2, 1) == 12);
     REQUIRE(mat4.Get(2, 2) == 13);
 
     Matrix<3, 1> mat6{10, 11, 12};
-    mat4.SetSubMatrix<3, 1, 0, 0>(mat6);
+    mat4.SetSubMatrix(0, 0, mat6);
     REQUIRE(mat4.Get(0, 0) == 10);
     REQUIRE(mat4.Get(1, 0) == 11);
     REQUIRE(mat4.Get(2, 0) == 12);
 
     Matrix<1, 3> mat7{10, 11, 12};
-    mat4.SetSubMatrix<1, 3, 0, 0>(mat7);
+    mat4.SetSubMatrix(0, 0, mat7);
     REQUIRE(mat4.Get(0, 0) == 10);
     REQUIRE(mat4.Get(0, 1) == 11);
     REQUIRE(mat4.Get(0, 2) == 12);
@@ -375,7 +375,7 @@ float matrixSum(const Matrix<rows, columns> &matrix) {
 
 // TODO: Add test for scalar division
 
-TEST_CASE("Normalization", "Matrix") {
+TEST_CASE("Euclidean Norm", "Matrix") {
 
   SECTION("2x2 Normalize") {
     Matrix<2, 2> mat1{1, 2, 3, 4};