From 75edad3d0adb2f9947748e01bf77a199bc4fbf28 Mon Sep 17 00:00:00 2001 From: Cynopolis Date: Mon, 2 Jun 2025 20:39:44 -0400 Subject: [PATCH] Made my own equally wrong QR factorization --- .vscode/settings.json | 5 +- src/Matrix.cpp | 93 +++++++++++++++++++++---------------- src/Matrix.hpp | 11 ++--- unit-tests/matrix-tests.cpp | 10 ++-- 4 files changed, 66 insertions(+), 53 deletions(-) diff --git a/.vscode/settings.json b/.vscode/settings.json index 8b34ef8..ebb6cf3 100644 --- a/.vscode/settings.json +++ b/.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" + ], } \ No newline at end of file diff --git a/src/Matrix.cpp b/src/Matrix.cpp index f253d89..a7fb1a2 100644 --- a/src/Matrix.cpp +++ b/src/Matrix.cpp @@ -451,16 +451,6 @@ float Matrix::EuclideanNorm() const { return sqrt(sum); } -template -float Matrix::L2Norm() const { - float sum{0}; - Matrix columnMatrix{}; - for (uint8_t column = 0; column < columns; column++) { - this.GetColumn(column, columnMatrix); - sum += columnMatrix.EuclideanNorm(); - } - return sqrt(sum); -} template template Matrix::SubMatrix() const { } template -template +template void Matrix::SetSubMatrix( + uint8_t rowOffset, uint8_t columnOffset, const Matrix &sub_matrix) { - static_assert(sub_rows + row_offset <= rows, - "The submatrix you're trying to set is out of bounds (rows)"); - static_assert( - sub_columns + column_offset <= columns, - "The submatrix you're trying to set is out of bounds (columns)"); + int16_t adjustedSubRows = sub_rows; + int16_t adjustedSubColumns = sub_columns; + int16_t adjustedRowOffset = rowOffset; + int16_t adjustedColumnOffset = columnOffset; - for (uint8_t row_idx{0}; row_idx < sub_rows; row_idx++) { - for (uint8_t column_idx{0}; column_idx < sub_columns; column_idx++) { - this->matrix[(row_idx + row_offset) * columns + column_idx + - column_offset] = sub_matrix.Get(row_idx, column_idx); + // a bunch of safety checks to make sure we don't overflow the matrix + if (sub_rows > rows) { + adjustedSubRows = rows; + } + if (sub_columns > columns) { + adjustedSubColumns = columns; + } + + if (adjustedSubRows + adjustedRowOffset >= rows) { + adjustedRowOffset = + std::max(0, static_cast(rows) - adjustedSubRows); + } + + if (adjustedSubColumns + adjustedColumnOffset >= columns) { + adjustedColumnOffset = + std::max(0, static_cast(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::QRDecomposition(Matrix &Q, Matrix &R) const { static_assert(columns <= rows, "QR decomposition requires columns <= rows"); - Matrix a_col, u, q_col, proj; Q.Fill(0); R.Fill(0); + Matrix a_col, e, u, Q_column_k{}; + Matrix<1, rows> a_T, e_T{}; - for (uint8_t k = 0; k < columns; ++k) { - this->GetColumn(k, a_col); + for (uint8_t column = 0; column < columns; column++) { + this->GetColumn(column, a_col); u = a_col; - - for (uint8_t j = 0; j < k; ++j) { - Q.GetColumn(j, q_col); - float r_jk = Matrix::DotProduct( - q_col, u); // FIXED: use u instead of a_col - R[j][k] = r_jk; - - proj = q_col * r_jk; - u = u - proj; + // ----------------------- + // ----- CALCULATE Q ----- + // ----------------------- + for (uint8_t k = 0; k <= column; k++) { + Q.GetColumn(k, Q_column_k); + Matrix<1, rows> Q_column_k_T = Q_column_k.Transpose(); + u = u - Q_column_k * (Q_column_k_T * a_col); } - - float norm = sqrt(Matrix::DotProduct(u, u)); - if (norm < 1e-12f) - norm = 1e-12f; // for stability - - for (uint8_t i = 0; i < rows; ++i) { - Q[i][k] = u[i][0] / norm; + float norm = u.EuclideanNorm(); + if (norm > 1e-4) { + u = u / norm; + } else { + u.Fill(0); } + 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); + } } } diff --git a/src/Matrix.hpp b/src/Matrix.hpp index 0abd9b1..47eb5cb 100644 --- a/src/Matrix.hpp +++ b/src/Matrix.hpp @@ -129,13 +129,12 @@ public: Matrix 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 SubMatrix() const; - template - void SetSubMatrix(const Matrix &sub_matrix); + template + void SetSubMatrix(uint8_t rowOffset, uint8_t columnOffset, + const Matrix &sub_matrix); /** * @brief take the dot product of the two vectors diff --git a/unit-tests/matrix-tests.cpp b/unit-tests/matrix-tests.cpp index f218b6b..896b34e 100644 --- a/unit-tests/matrix-tests.cpp +++ b/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 &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};