feat: made max iterations obligatory argument

code/eigenvalue_methods.py

@@ -2,7 +2,7 @@ from linear_algebra_utils import LinearAlgebraUtils
 
 class EigenvalueMethods:
     @staticmethod
-    def power_method(A, max_iter=1000, tol=1e-6):
+    def power_method(A, max_iter, tol=1e-6):
         n = len(A)
         x = [1] * n
         lambda_old = 0
@@ -18,7 +18,7 @@ class EigenvalueMethods:
         return lambda_new
 
     @staticmethod
-    def inverse_power_method(A, max_iter=1000, tol=1e-6):
+    def inverse_power_method(A, max_iter, tol=1e-6):
         n = len(A)
         I = [[1 if i == j else 0 for j in range(n)] for i in range(n)]
         A_inv = [[I[i][j] - A[i][j] for j in range(n)] for i in range(n)]

code/linear_algebra_utils.py

@@ -9,6 +9,14 @@ class LinearAlgebraUtils:
     def matrix_vector_multiply(A, x):
         return [LinearAlgebraUtils.dot_product(row, x) for row in A]
     
+    @staticmethod
+    def vector_norm(v):
+        return sum(x*x for x in v)**0.5
+
+    @staticmethod
+    def vector_scalar_divide(x, scalar):
+        return [xi / scalar for xi in x]
+
     @staticmethod
     def matrix_scalar_multiply(A, w):
         return [[w * A[i][j] for j in range(len(A[0]))] for i in range(len(A))]
@@ -25,18 +33,11 @@ class LinearAlgebraUtils:
     def scalar_matrix_multiply(omega, vector):
         return [omega * x for x in vector]
 
-    @staticmethod
-    def vector_norm(v):
-        return sum(x*x for x in v)**0.5
 
     @staticmethod
     def matrix_norm(A):
         return math.sqrt(sum(sum(element ** 2 for element in row) for row in A))
 
-    @staticmethod
-    def vector_scalar_divide(x, scalar):
-        return [xi / scalar for xi in x]
-
     @staticmethod
     def matrix_matrix_subtraction(A, B):
         return [[A[i][j] - B[i][j] for j in range(len(A[0]))] for i in range(len(A))]

code/matrix_generator.py

@@ -1,9 +1,24 @@
 import numpy as np
 
 class MatrixGenerator:
+    @staticmethod
+    def generate_spd_matrix(n: int) -> np.ndarray:
+        """
+        Generates a random symmetric positive definite matrix of size n x n.
+
+        Parameters:
+        n (int): The size of the matrix to generate.
+
+        Returns:
+        np.ndarray: A symmetric positive definite matrix of size n x n.
+        """
+        A = np.random.rand(n, n)
+        spd_matrix = np.dot(A, A.T) + n * np.eye(n)  # Adding n*I ensures positive definiteness
+        return spd_matrix
+
     @staticmethod
     def generate_random_matrix_and_vector(size):
-        A = np.random.uniform(-1, 1, (size, size))
+        A = MatrixGenerator.generate_spd_matrix(size)
         b = np.random.uniform(-1, 1, size)
         return A, b
 

code/richardson_method.py

@@ -3,29 +3,36 @@ from eigenvalue_methods import EigenvalueMethods
 from matrix_generator import MatrixGenerator
 
 class RichardsonMethod:
-    def __init__(self, A, b, size: int, x0=None, max_iterations=1000, tol=1e-5):
+    def __init__(self, A, b, max_iterations, size: int, x0=None, tol=1e-5):
         self.A = A
         self.b = b
         self.x0 = x0 if x0 is not None else [0.0] * len(b)
         self.max_iterations = max_iterations
         self.tol = tol
         self.I = MatrixGenerator.generate_identity_matrix(size)
-        self.lambda_min = EigenvalueMethods.inverse_power_method(self.A)
+        self.lambda_min, self.lambda_max = RichardsonMethod.calculate_eigenvalues(self.A, max_iterations)
         if self.lambda_min < 0:
             raise ValueError("Matrix A is not positive semi-definite.")
-        self.lambda_max = EigenvalueMethods.power_method(self.A)
+        self.omega = RichardsonMethod.calculate_omega(self.lambda_min, self.lambda_max)
-        self.omega = 2 / (self.lambda_min + self.lambda_max)
 
+    @staticmethod
+    def calculate_eigenvalues(A, max_iterations):
+        return EigenvalueMethods.inverse_power_method(A, max_iterations), EigenvalueMethods.power_method(A, max_iterations)
 
-    def will_converge(self) -> bool:
+    @staticmethod
-        wA = LinearAlgebraUtils.matrix_scalar_multiply(self.A, self.omega)
+    def calculate_omega(lambda_min, lambda_max):
-        IMinuswA = LinearAlgebraUtils.matrix_matrix_subtraction(self.I, wA)
+        return 2 / (lambda_min + lambda_max)
+    
+    @staticmethod
+    def convergence_norm(A, omega, I) -> bool:
+        wA = LinearAlgebraUtils.matrix_scalar_multiply(A, omega)
+        IMinuswA = LinearAlgebraUtils.matrix_matrix_subtraction(I, wA)
         norm = LinearAlgebraUtils.matrix_norm(IMinuswA)
-        return norm < 1
+        return norm
 
     def solve(self):
         x = self.x0[:]
-        if not self.will_converge():
+        if RichardsonMethod.convergence_norm(self.A, self.omega, self.I) >= 1:
             return "Richardson method for those values will NOT converge"
 
         for iteration in range(self.max_iterations):

code/tests.py

@@ -4,12 +4,39 @@ from scipy.sparse.linalg import cg
 from matrix_generator import MatrixGenerator
 from richardson_method import RichardsonMethod
 
+def calculate_norm_numpy(I, w, A):
+    # Calculate the difference between I and w * A
+    difference = I - w * A
+    
+    # Calculate the Euclidean norm of the difference
+    norm = np.linalg.norm(difference)
+    
+    return norm
+
+def calculate_eigenvalues(A):
+    # Calculate the eigenvalues of matrix A
+    eigenvalues = np.linalg.eigvals(A)
+    
+    # Find the minimum and maximum eigenvalues
+    lambda_min = np.min(eigenvalues)
+    lambda_max = np.max(eigenvalues)
+    
+    return lambda_min, lambda_max
+
 @pytest.mark.parametrize("n", [2, 3, 4, 5, 10, 20, 50, 100])
 def test_richardson_vs_cg(n: int):
+    print("matrix size: ", n)
     tolerance = 1e-5
+    max_iterations=1000
     A, b = MatrixGenerator.generate_random_matrix_and_vector(n)
+    lambda_min, lambda_max = calculate_eigenvalues(A)
+    print("eigenvalues: ", lambda_min, lambda_max, RichardsonMethod.calculate_eigenvalues(A, max_iterations))
+    omega = 2 / (lambda_min + lambda_max)
+    print("omega: ", omega, RichardsonMethod.calculate_omega(lambda_min, lambda_max))
+    I = np.eye(n)
+    print("norms: ", calculate_norm_numpy(I, omega, A), RichardsonMethod.convergence_norm(A, omega, I))
     
-    richardson_solver = RichardsonMethod(A, b, size=n, max_iterations=1000, tol=1e-7)
+    richardson_solver = RichardsonMethod(A, b, max_iterations, size=n, tol=1e-7)
     solution_richardson = richardson_solver.solve()
     
     solution_cg, info = cg(A, b)