From 0b55d293767eb5652ee0e4cafef46cfc25c96cc0 Mon Sep 17 00:00:00 2001 From: Quinn Henthorne Date: Wed, 18 Dec 2024 16:30:53 -0500 Subject: [PATCH] Working on getting the QR decomposition to compile --- .vscode/settings.json | 3 +- Matrix.cpp | 97 +++++++++++++++++++++++++++++++------ Matrix.hpp | 69 +++++++++++++++++++++++++- unit-tests/matrix-tests.cpp | 12 +++++ 4 files changed, 162 insertions(+), 19 deletions(-) diff --git a/.vscode/settings.json b/.vscode/settings.json index 71e5c54..da785c2 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -70,7 +70,8 @@ "thread": "cpp", "typeinfo": "cpp", "variant": "cpp", - "shared_mutex": "cpp" + "shared_mutex": "cpp", + "complex": "cpp" }, "clangd.enable": false, "C_Cpp.dimInactiveRegions": false diff --git a/Matrix.cpp b/Matrix.cpp index 11954b7..736982b 100644 --- a/Matrix.cpp +++ b/Matrix.cpp @@ -17,6 +17,16 @@ Matrix::Matrix(const std::array &array) { this->setMatrixToArray(array); } +template +Matrix::Matrix(const Matrix &other) { + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + this->matrix[row_idx * columns + column_idx] = + other.Get(row_idx, column_idx); + } + } +} + template template Matrix::Matrix(Args... args) { @@ -30,16 +40,6 @@ Matrix::Matrix(Args... args) { memcpy(this->matrix.begin(), initList.begin(), minSize * sizeof(float)); } -template -Matrix::Matrix(const Matrix &other) { - for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { - for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { - this->matrix[row_idx * columns + column_idx] = - other.Get(row_idx, column_idx); - } - } -} - template void Matrix::setMatrixToArray( const std::array &array) { @@ -86,7 +86,7 @@ Matrix::Sub(const Matrix &other, template template -Matrix & +Matrix & Matrix::Mult(const Matrix &other, Matrix &result) const { // allocate some buffers for all of our dot products @@ -353,11 +353,7 @@ float Matrix::dotProduct(const Matrix &vec1, template void Matrix::Fill(float value) { - for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { - for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { - this->matrix[row_idx * columns + column_idx] = value; - } - } + this->matrix.fill(value); } template @@ -440,4 +436,73 @@ Matrix::Normalize(Matrix &result) const { return result; } +template +Matrix Matrix::Eye() { + Matrix i_matrix; + i_matrix.Fill(0); + for (uint8_t i{0}; i < rows; i++) { + i_matrix[i][i] = 1; + } + return i_matrix; +} + +template +void Matrix::QR_Decomposition(Matrix &Q, + Matrix &R) const { + Q = Matrix::Eye(); // Q starts as the identity matrix + R = *this; // R starts as a copy of this matrix (For this algorithm we'll call + // this matrix A) + + for (uint8_t row{0}; row < rows; row++) { + // compute the householder vector + const uint8_t houseHoldVectorSize{rows - row}; + const uint8_t subMatrixSize{columns - row}; + Matrix x{}; + this->SubMatrix(row, row, x); + + Matrix e1{}; + e1.Fill(0); + if (x[0][0] >= 0) { + e1[0][0] = x.Norm(); + } else { + e1[0][0] = -x.Norm(); + } + + Matrix v = x + e1; + v = v * (1 / v.Norm()); // normalize V + + // ************************************ + // Apply the reflection to the R matrix + // ************************************ + // initialize R's submatrix + Matrix R_subMatrix{}; + R.SubMatrix(row, row, R_subMatrix); + // create some temporary buffers + Matrix<1, subMatrixSize> vR{}; + Matrix<1, houseHoldVectorSize> v_T{}; + v.Transpose(v_T); + Matrix vR_outer{}; + // calculate the reflection + R_subMatrix = + R_subMatrix - 2 * Matrix::OuterProduct( + v_T, v_T.Mult(R_subMatrix, vR), vR_outer); + // save the reflection back to R + R.CopySubMatrixInto(row, row, R_subMatrix); + + // ************************************ + // Apply the reflection to the Q matrix + // ************************************ + // initialize Q's submatrix + Matrix Q_subMatrix{}; + Q.SubMatrix(0, row, Q_subMatrix); + // create some temporary buffers + Matrix Qv{}; + Matrix Qv_outer{}; + + Q_subMatrix = Q_subMatrix - 2 * Matrix::OuterProduct( + Q_subMatrix.Mult(v, Qv), v, Qv_outer); + Q.CopySubMatrixInto(0, row, Q_subMatrix); + } +} + #endif // MATRIX_H_ \ No newline at end of file diff --git a/Matrix.hpp b/Matrix.hpp index 4b5479b..b52eadc 100644 --- a/Matrix.hpp +++ b/Matrix.hpp @@ -35,6 +35,7 @@ public: * @brief Initialize a matrix directly with any number of arguments */ template Matrix(Args... args); + /** * @brief Set all elements in this to value */ @@ -64,8 +65,8 @@ public: * @param result A buffer to store the result into */ template - Matrix &Mult(const Matrix &other, - Matrix &result) const; + Matrix &Mult(const Matrix &other, + Matrix &result) const; /** * @brief Multiply the matrix by a scalar @@ -125,6 +126,66 @@ public: */ Matrix &Normalize(Matrix &result) const; + /** + * @brief return an identity matrix of the specified size + */ + static Matrix Eye(); + + /** + * @brief write a copy of a sub matrix into the given result matrix. + * @param rowIndex The row index to start the copy from + * @param columnIndex the column index to start the copy from + * @param result the matrix buffer to write the sub matrix into. The size of + * the matrix buffer allows the function to determine the end indices of the + * sub matrix + */ + template + Matrix & + SubMatrix(uint8_t rowIndex, uint8_t columnIndex, + Matrix &result) const { + return result; + } + + /** + * @brief write a copy of a sub matrix into this matrix starting at the given + * idnex. + * @param rowIndex The row index to start the copy from + * @param columnIndex the column index to start the copy from + * @param subMatrix The submatrix to copy into this matrix. The size of + * the matrix buffer allows the function to determine the end indices of the + * sub matrix + */ + template + void CopySubMatrixInto(uint8_t rowIndex, uint8_t columnIndex, + const Matrix &subMatrix) {} + + /** + * @brief Returns the norm of the matrix + */ + float Norm() { return 0; } + + template + static Matrix & + OuterProduct(const Matrix<1, vec1Length> &vec1, + const Matrix<1, vec2Length> &vec2, + Matrix &result) { + return result; + } + + template + static Matrix & + OuterProduct(const Matrix &vec1, + const Matrix &vec2, + Matrix &result) { + return result; + } + /** + * @brief Calulcate the QR decomposition of a matrix + * @param Q the + */ + void QR_Decomposition(Matrix &Q, + Matrix &R) const; + /** * @brief Get a row from the matrix * @param row_index the row index to get @@ -160,6 +221,10 @@ public: */ float Get(uint8_t row_index, uint8_t column_index) const; + // ******************************************************* + // ************** OPERATOR OVERRIDES ********************* + // ******************************************************* + /** * @brief get the specified row of the matrix returned as a reference to the * internal array diff --git a/unit-tests/matrix-tests.cpp b/unit-tests/matrix-tests.cpp index 09a3da9..a1d5ec9 100644 --- a/unit-tests/matrix-tests.cpp +++ b/unit-tests/matrix-tests.cpp @@ -182,6 +182,18 @@ TEST_CASE("Elementary Matrix Operations", "Matrix") { REQUIRE(minorMat4.Get(1, 0) == 4); REQUIRE(minorMat4.Get(1, 1) == 5); } + SECTION("Identity Matrix") { + Matrix<3, 3> mat4{Matrix<3, 3>::Eye()}; + REQUIRE(mat4.Get(0, 0) == 1); + REQUIRE(mat4.Get(0, 1) == 0); + REQUIRE(mat4.Get(0, 2) == 0); + REQUIRE(mat4.Get(1, 0) == 0); + REQUIRE(mat4.Get(1, 1) == 1); + REQUIRE(mat4.Get(1, 2) == 0); + REQUIRE(mat4.Get(2, 0) == 0); + REQUIRE(mat4.Get(2, 1) == 0); + REQUIRE(mat4.Get(2, 2) == 1); + } SECTION("Determinant") { float det1 = mat1.Det();