[Matrix][decompositions] Add Eigenvalue Decomposition for Symmetric Matrices

numojo/core/matrix.mojo

@@ -18,6 +18,8 @@ from numojo.core.ndarray import NDArray
 from numojo.core.own_data import OwnData
 from numojo.core.utility import _get_offset
 from numojo.routines.manipulation import broadcast_to, reorder_layout
+from numojo.routines.linalg.misc import issymmetric
+
 
 # ===----------------------------------------------------------------------===#
 # Matrix struct
@@ -1209,10 +1211,16 @@ struct Matrix[dtype: DType = DType.float64](
 
     fn trace(self) raises -> Scalar[dtype]:
         """
-        Transpose of matrix.
+        Trace of matrix.
         """
         return numojo.linalg.trace(self)
 
+    fn issymmetric(self) -> Bool:
+        """
+        Transpose of matrix.
+        """
+        return issymmetric(self)
+
     fn transpose(self) -> Self:
         """
         Transpose of matrix.

numojo/routines/linalg/__init__.mojo

@@ -1,4 +1,4 @@
-from .decompositions import lu_decomposition, qr
+from .decompositions import lu_decomposition, qr, eig
 from .norms import det, trace
 from .products import cross, dot, matmul
 from .solving import inv, solve, lstsq

numojo/routines/linalg/decompositions.mojo

@@ -7,10 +7,11 @@ from algorithm import parallelize, vectorize
 import math as builtin_math
 
 from numojo.core.ndarray import NDArray
-from numojo.core.matrix import Matrix
+from numojo.core.matrix import Matrix, issymmetric
 from numojo.routines.creation import zeros, eye, full
 
 
+@always_inline
 fn _compute_householder[
     dtype: DType
 ](mut H: Matrix[dtype], mut R: Matrix[dtype], work_index: Int) raises -> None:
@@ -69,6 +70,7 @@ fn _compute_householder[
     vectorize[scaling_factor_increment_vec, simd_width](rRows)
 
 
+@always_inline
 fn _apply_householder[
     dtype: DType
 ](
@@ -396,3 +398,80 @@ fn qr[
             R = R[:inner, :]
 
     return Q^, R^
+
+
+# ===----------------------------------------------------------------------=== #
+# Eigenvalue Decomposition (symmetric) via the QR algorithm
+# ===----------------------------------------------------------------------=== #
+
+
+fn eig[
+    dtype: DType
+](
+    A: Matrix[dtype],
+    tol: Scalar[dtype] = 1.0e-12,
+    max_iter: Int = 10000,
+) raises -> Tuple[Matrix[dtype], Matrix[dtype]]:
+    """
+    Computes the eigenvalue decomposition for symmetric matrices using the QR algorithm.
+    For best performance, the input matrix should be in column-major order.
+
+    Args:
+        A: The input matrix. Must be square and symmetric.
+        tol: Convergence tolerance for off-diagonal elements.
+        max_iter: Maximum number of iterations for the QR algorithm.
+
+    Returns:
+        A tuple `(Q, D)` where:
+            - Q: A matrix whose columns are the eigenvectors
+            - D: A diagonal matrix containing the eigenvalues
+
+    Raises:
+        Error: If the matrix is not square or symmetric
+        Error: If the algorithm does not converge within max_iter iterations
+    """
+    if A.shape[0] != A.shape[1]:
+        raise Error("Matrix is not square.")
+
+    var n = A.shape[0]
+    if not issymmetric(A):
+        raise Error("Matrix is not symmetric.")
+
+    var T: Matrix[dtype]
+    if A.flags.C_CONTIGUOUS:
+        T = A.reorder_layout()
+    else:
+        T = A
+
+    var Q_total = Matrix.identity[dtype](n)
+
+    for k in range(max_iter):
+        var Qk: Matrix[dtype]
+        var Rk: Matrix[dtype]
+        Qk, Rk = qr(T, mode="complete")
+
+        T = Rk @ Qk
+        Q_total = Q_total @ Qk
+
+        var offdiag_norm: Scalar[dtype] = 0
+        for i in range(n):
+            for j in range(i + 1, n):
+                var v = T._load(i, j)
+                offdiag_norm += v * v
+        if builtin_math.sqrt(offdiag_norm) < tol:
+            break
+    else:
+        raise Error(
+            String("QR algorithm did not converge in {} iterations.").format(
+                max_iter
+            )
+        )
+
+    var D = Matrix.zeros[dtype](shape=(n, n), order=A.order())
+    for i in range(n):
+        D._store(i, i, T._load(i, i))
+
+    if A.flags.C_CONTIGUOUS:
+        Q_total = Q_total.reorder_layout()
+
+    return Q_total^, D^

numojo/routines/linalg/misc.mojo

@@ -9,6 +9,9 @@
 # Miscellaneous Linear Algebra Routines
 # ===----------------------------------------------------------------------=== #
 
+from sys import simdwidthof
+from algorithm import parallelize, vectorize
+
 from numojo.core.ndarray import NDArray
 
 
@@ -59,3 +62,40 @@ fn diagonal[
             res.item(i) = a.item(i - offset, i)
 
     return res
+
+
+fn issymmetric[
+    dtype: DType
+](
+    A: Matrix[dtype], rtol: Scalar[dtype] = 1e-5, atol: Scalar[dtype] = 1e-8
+) -> Bool:
+    """
+    Returns True if A is symmetric, False otherwise.
+
+    Parameters:
+        dtype: Data type of the Matrix Elements.
+
+    Args:
+        A: A Matrix.
+        rtol: Relative tolerance for comparison.
+        atol: Absolute tolerance for comparison.
+
+    Returns:
+        True if the array is symmetric, False otherwise.
+    """
+
+    if A.shape[0] != A.shape[1]:
+        return False
+
+    var n = A.shape[0]
+
+    for i in range(n):
+        for j in range(i + 1, n):
+            var a_ij = A._load(i, j)
+            var a_ji = A._load(j, i)
+            var diff = abs(a_ij - a_ji)
+            var allowed_error = atol + rtol * max(abs(a_ij), abs(a_ji))
+            if diff > allowed_error:
+                return False
+
+    return True

tests/core/test_matrix.mojo

@@ -319,6 +319,79 @@ def test_qr_decomposition_asym_complete2():
     assert_true(np.allclose(A_test.to_numpy(), A.to_numpy(), atol=1e-14))
 
 
+def test_eigen_decomposition():
+    var np = Python.import_module("numpy")
+
+    # Create a symmetric matrix by adding a matrix to its transpose
+    var A_random = Matrix.rand[f64]((10, 10), order=order)
+    var A = A_random + A_random.transpose()
+    var Anp = A.to_numpy()
+
+    # Compute eigendecomposition
+    Q, Lambda = nm.linalg.eig(A)
+
+    # Use NumPy for comparison
+    namedtuple = np.linalg.eig(Anp)
+
+    np_eigenvalues = namedtuple.eigenvalues
+    print(np_eigenvalues)
+    print(Lambda.to_numpy())
+    print(np.diag(Lambda.to_numpy()))
+
+    # Sort eigenvalues and eigenvectors for comparison (numpy doesn't guarantee order)
+    var np_sorted_eigenvalues = np.sort(np_eigenvalues)
+    var eigenvalues = np.diag(Lambda.to_numpy())
+    var sorted_eigenvalues = np.sort(eigenvalues)
+
+    assert_true(
+        np.allclose(sorted_eigenvalues, np_sorted_eigenvalues, atol=1e-10),
+        "Eigenvalues don't match expected values",
+    )
+
+    # Check that eigenvectors are orthogonal (Q^T Q = I)
+    var id = Q.transpose() @ Q
+    assert_true(
+        np.allclose(id.to_numpy(), np.eye(Q.shape[0]), atol=1e-10),
+        "Eigenvectors are not orthogonal",
+    )
+
+    # Check that A = Q * Lambda * Q^T (eigendecomposition property)
+    var A_reconstructed = Q @ Lambda @ Q.transpose()
+    print(A_reconstructed - A)
+    assert_true(
+        np.allclose(A_reconstructed.to_numpy(), Anp, atol=1e-10),
+        "A ≠ Q * Lambda * Q^T",
+    )
+
+    # Verify A*v = λ*v for each eigenvector and eigenvalue
+    for i in range(A.shape[0]):
+        var eigenvector = Matrix.zeros[f64]((A.shape[0], 1), order=order)
+        for j in range(A.shape[0]):
+            eigenvector[j, 0] = Q[j, i]
+
+        var Av = A @ eigenvector
+        var lambda_times_v = eigenvector * Lambda[i, i]
+
+        assert_true(
+            np.allclose(Av.to_numpy(), lambda_times_v.to_numpy(), atol=1e-10),
+            "Eigenvector verification failed: A*v ≠ λ*v",
+        )
+
+    # Verify A*v = λ*v for each eigenvector and eigenvalue
+    for i in range(A.shape[0]):
+        var eigenvector = Matrix.zeros[f64]((A.shape[0], 1), order=order)
+        for j in range(A.shape[0]):
+            eigenvector[j, 0] = Q[j, i]
+
+        var Av = A @ eigenvector
+        var lambda_times_v = eigenvector * Lambda[i, i]
+
+        assert_true(
+            np.allclose(Av.to_numpy(), lambda_times_v.to_numpy(), atol=1e-10),
+            "Eigenvector verification failed: A*v ≠ λ*v",
+        )
+
+
 # ===-----------------------------------------------------------------------===#
 # Mathematics
 # ===-----------------------------------------------------------------------===#