Working on getting the QR decomposition to compile

.gitignore

@@ -1 +1,2 @@
-build/
+build/
+venv/

.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

Matrix.cpp

@@ -17,6 +17,16 @@ Matrix<rows, columns>::Matrix(const std::array<float, rows * columns> &array) {
   this->setMatrixToArray(array);
 }
 
+template <uint8_t rows, uint8_t columns>
+Matrix<rows, columns>::Matrix(const Matrix<rows, columns> &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 <uint8_t rows, uint8_t columns>
 template <typename... Args>
 Matrix<rows, columns>::Matrix(Args... args) {
@@ -30,16 +40,6 @@ Matrix<rows, columns>::Matrix(Args... args) {
   memcpy(this->matrix.begin(), initList.begin(), minSize * sizeof(float));
 }
 
-template <uint8_t rows, uint8_t columns>
-Matrix<rows, columns>::Matrix(const Matrix<rows, columns> &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 <uint8_t rows, uint8_t columns>
 void Matrix<rows, columns>::setMatrixToArray(
     const std::array<float, rows * columns> &array) {
@@ -86,7 +86,7 @@ Matrix<rows, columns>::Sub(const Matrix<rows, columns> &other,
 
 template <uint8_t rows, uint8_t columns>
 template <uint8_t other_columns>
-Matrix<rows, columns> &
+Matrix<rows, other_columns> &
 Matrix<rows, columns>::Mult(const Matrix<columns, other_columns> &other,
                             Matrix<rows, other_columns> &result) const {
   // allocate some buffers for all of our dot products
@@ -353,11 +353,7 @@ float Matrix<rows, columns>::dotProduct(const Matrix<vector_size, 1> &vec1,
 
 template <uint8_t rows, uint8_t columns>
 void Matrix<rows, columns>::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 <uint8_t rows, uint8_t columns>
@@ -440,4 +436,73 @@ Matrix<rows, columns>::Normalize(Matrix<rows, columns> &result) const {
   return result;
 }
 
+template <uint8_t rows, uint8_t columns>
+Matrix<rows, rows> Matrix<rows, columns>::Eye() {
+  Matrix<rows, rows> i_matrix;
+  i_matrix.Fill(0);
+  for (uint8_t i{0}; i < rows; i++) {
+    i_matrix[i][i] = 1;
+  }
+  return i_matrix;
+}
+
+template <uint8_t rows, uint8_t columns>
+void Matrix<rows, columns>::QR_Decomposition(Matrix<rows, columns> &Q,
+                                             Matrix<rows, columns> &R) const {
+  Q = Matrix<rows, columns>::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<houseHoldVectorSize, 1> x{};
+    this->SubMatrix(row, row, x);
+
+    Matrix<houseHoldVectorSize, 1> e1{};
+    e1.Fill(0);
+    if (x[0][0] >= 0) {
+      e1[0][0] = x.Norm();
+    } else {
+      e1[0][0] = -x.Norm();
+    }
+
+    Matrix<houseHoldVectorSize, 1> v = x + e1;
+    v = v * (1 / v.Norm()); // normalize V
+
+    // ************************************
+    // Apply the reflection to the R matrix
+    // ************************************
+    // initialize R's submatrix
+    Matrix<houseHoldVectorSize, subMatrixSize> 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<houseHoldVectorSize, subMatrixSize> vR_outer{};
+    // calculate the reflection
+    R_subMatrix =
+        R_subMatrix - 2 * Matrix<rows, columns>::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<rows, houseHoldVectorSize> Q_subMatrix{};
+    Q.SubMatrix(0, row, Q_subMatrix);
+    // create some temporary buffers
+    Matrix<rows, 1> Qv{};
+    Matrix<rows, houseHoldVectorSize> Qv_outer{};
+
+    Q_subMatrix = Q_subMatrix - 2 * Matrix<rows, columns>::OuterProduct(
+                                        Q_subMatrix.Mult(v, Qv), v, Qv_outer);
+    Q.CopySubMatrixInto(0, row, Q_subMatrix);
+  }
+}
+
 #endif // MATRIX_H_

Matrix.hpp

@@ -35,6 +35,7 @@ public:
    * @brief Initialize a matrix directly with any number of arguments
    */
   template <typename... Args> 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 <uint8_t other_columns>
-  Matrix<rows, columns> &Mult(const Matrix<columns, other_columns> &other,
-                              Matrix<rows, other_columns> &result) const;
+  Matrix<rows, other_columns> &Mult(const Matrix<columns, other_columns> &other,
+                                    Matrix<rows, other_columns> &result) const;
 
   /**
    * @brief Multiply the matrix by a scalar
@@ -125,6 +126,66 @@ public:
    */
   Matrix<rows, columns> &Normalize(Matrix<rows, columns> &result) const;
 
+  /**
+   * @brief return an identity matrix of the specified size
+   */
+  static Matrix<rows, rows> 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 <uint8_t subRows, uint8_t subColumns>
+  Matrix<subRows, subColumns> &
+  SubMatrix(uint8_t rowIndex, uint8_t columnIndex,
+            Matrix<subRows, subColumns> &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 <uint8_t subRows, uint8_t subColumns>
+  void CopySubMatrixInto(uint8_t rowIndex, uint8_t columnIndex,
+                         const Matrix<subRows, subColumns> &subMatrix) {}
+
+  /**
+   * @brief Returns the norm of the matrix
+   */
+  float Norm() { return 0; }
+
+  template <uint8_t vec1Length, uint8_t vec2Length>
+  static Matrix<vec1Length, vec2Length> &
+  OuterProduct(const Matrix<1, vec1Length> &vec1,
+               const Matrix<1, vec2Length> &vec2,
+               Matrix<vec1Length, vec2Length> &result) {
+    return result;
+  }
+
+  template <uint8_t vec1Length, uint8_t vec2Length>
+  static Matrix<vec1Length, vec2Length> &
+  OuterProduct(const Matrix<vec1Length, 1> &vec1,
+               const Matrix<vec2Length, 1> &vec2,
+               Matrix<vec1Length, vec2Length> &result) {
+    return result;
+  }
+  /**
+   * @brief Calulcate the QR decomposition of a matrix
+   * @param Q the
+   */
+  void QR_Decomposition(Matrix<rows, columns> &Q,
+                        Matrix<rows, columns> &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

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)

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();