diff --git a/.gitignore b/.gitignore index d163863..f188146 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ -build/ \ No newline at end of file +build/ +venv/ \ No newline at end of file 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/qr-decom.py b/qr-decom.py new file mode 100644 index 0000000..fa3b2cb --- /dev/null +++ b/qr-decom.py @@ -0,0 +1,159 @@ +import numpy as np + +# QR decomposition using the householder reflection method +def householder_reflection(A): + """ + Perform QR decomposition using Householder reflection. + + Arguments: + A -- A matrix to be decomposed (m x n). + + Returns: + Q -- Orthogonal matrix (m x m). + R -- Upper triangular matrix (m x n). + """ + A = A.astype(float) # Ensure the matrix is of type float + m, n = A.shape + Q = np.eye(m) # Initialize Q as an identity matrix + R = A.copy() # R starts as a copy of A + + # Apply Householder reflections for each column + for k in range(n): + # Step 1: Compute the Householder vector + x = R[k:m, k] + e1 = np.zeros_like(x) + e1[0] = np.linalg.norm(x) if x[0] >= 0 else -np.linalg.norm(x) + v = x + e1 + v = v / np.linalg.norm(v) # Normalize v + + # Step 2: Apply the reflection to the matrix + R[k:m, k:n] = R[k:m, k:n] - 2 * np.outer(v, v.T @ R[k:m, k:n]) + + # Step 3: Apply the reflection to Q + Q[:, k:m] = Q[:, k:m] - 2 * np.outer(Q[:, k:m] @ v, v) + + # The resulting Q and R are the QR decomposition + return Q, R + +# Example usage +A = np.array([[12, -51, 4], + [6, 167, -68], + [-4, 24, -41]]) + +Q, R = householder_reflection(A) +print("Q matrix:") +print(Q) +print("\nR matrix:") +print(R) +print("Multiplied Together:") +print(Q@R) + +def svd_decomposition(A): + """ + Perform Singular Value Decomposition (SVD) from scratch. + + Arguments: + A -- The matrix to be decomposed (m x n). + + Returns: + U -- Orthogonal matrix of left singular vectors (m x m). + Sigma -- Diagonal matrix of singular values (m x n). + Vt -- Orthogonal matrix of right singular vectors (n x n). + """ + # Step 1: Compute A^T A + AtA = np.dot(A.T, A) # A transpose multiplied by A + + # Step 2: Compute the eigenvalues and eigenvectors of A^T A + eigenvalues, V = np.linalg.eig(AtA) + + # Step 3: Sort eigenvalues in descending order and sort V accordingly + sorted_indices = np.argsort(eigenvalues)[::-1] # Indices to sort eigenvalues in descending order + eigenvalues = eigenvalues[sorted_indices] + V = V[:, sorted_indices] + + # Step 4: Compute the singular values (sqrt of eigenvalues) + singular_values = np.sqrt(eigenvalues) + + # Step 5: Construct the Sigma matrix + m, n = A.shape + Sigma = np.zeros((m, n)) # Initialize Sigma as a zero matrix + for i in range(min(m, n)): + Sigma[i, i] = singular_values[i] # Place the singular values on the diagonal + + # Step 6: Compute the U matrix using A * V = U * Sigma + U = np.dot(A, V) # A * V gives us the unnormalized U + # Normalize the columns of U + for i in range(U.shape[1]): + U[:, i] = U[:, i] / singular_values[i] # Normalize each column by the corresponding singular value + + # Step 7: Return U, Sigma, Vt + return U, Sigma, V.T # V.T is the transpose of V + +# Example usage +A = np.array([[12, -51, 4], + [6, 167, -68], + [-4, 24, -41]]) + +U, Sigma, Vt = svd_decomposition(A) + +print("\nSVD DECOMPOSITION\nU matrix:") +print(U) +print("\nSigma matrix:") +print(Sigma) +print("\nVt matrix:") +print(Vt) +print("Multiplied together:") +print(U@Sigma@Vt) + +def eigen_decomposition_qr(A, max_iter=1000, tol=1e-9): + """ + Compute the eigenvalues and eigenvectors of a matrix A using the QR algorithm + with QR decomposition. + + Arguments: + A -- A square matrix (n x n). + max_iter -- Maximum number of iterations for convergence (default 1000). + tol -- Tolerance for convergence (default 1e-9). + + Returns: + eigenvalues -- List of eigenvalues. + eigenvectors -- Matrix of eigenvectors. + """ + # Make a copy of A to perform the iteration + A_copy = A.copy() + n = A_copy.shape[0] + + # Initialize the matrix for eigenvectors (this will accumulate the Q matrices) + eigenvectors = np.eye(n) + + # Perform QR iterations + for _ in range(max_iter): + # Perform QR decomposition on A_copy + Q, R = householder_reflection(A_copy) + + # Update A_copy to be R * Q (QR algorithm step) + A_copy = R @ Q + + # Accumulate the eigenvectors + eigenvectors = eigenvectors @ Q + + # Check for convergence: if the off-diagonal elements are small enough, we stop + off_diagonal_norm = np.linalg.norm(np.tril(A_copy, -1)) # Norm of the lower triangle (off-diagonal) + if off_diagonal_norm < tol: + break + + # The eigenvalues are the diagonal elements of the matrix A_copy + eigenvalues = np.diag(A_copy) + + return eigenvalues, eigenvectors + +# Example usage +A = np.array([[12, -51, 4], + [6, 167, -68], + [-4, 24, -41]]) + +eigenvalues, eigenvectors = eigen_decomposition_qr(A) + + +print("\n\nEigenvalues:", eigenvalues) +print("Eigenvectors:\n", eigenvectors) \ No newline at end of file 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();