Add sequential version of Richardson's algorithm and first tests

code/Richardson_numpy.py

@@ -0,0 +1,70 @@
+import numpy as np
+from scipy.sparse.linalg import cg
+
+def generate_random_matrix_and_vector(size):
+    A = np.random.uniform(-1, 1, (size, size))
+    A = np.dot(A.T, A) + np.eye(size) * size  # dodanie `size` do przekątnej zwiększa wartości własne -> obejście problemu z overflow
+    b = np.random.uniform(-1, 1, size)
+    return A, b
+
+def richardson(A, b, x0=None, max_iterations=1000, tol=1e-5):
+
+    if x0 is None:
+        x0 = np.zeros_like(b, dtype=float)
+    
+    x = np.array(x0, dtype=float)
+    n = len(b)
+
+    eigenvalues = np.linalg.eigvals(A)
+    lambda_min = np.min(eigenvalues)
+    lambda_max = np.max(eigenvalues)
+
+    if lambda_min < 0:
+        raise ValueError("Matrix A is not positive semi-definite.")
+    
+    omega = 2/(lambda_min + lambda_max)
+
+    for iteration in range(max_iterations):
+        Ax = np.dot(A, x)
+        residual = b - Ax
+
+        if np.linalg.norm(residual) < tol:
+            print('Convergence achieved.')
+            return x
+        
+        x = x + omega * residual
+
+    print('Maximum number of iterations reached without convergence.')
+    return x
+
+def run_tests():
+    test_sizes = [2, 3, 4, 5, 10, 20, 50, 100]
+    tolerance = 1e-5
+    
+    for n in test_sizes:
+        print(f"\nRunning test for n = {n}")
+        
+        A, b = generate_random_matrix_and_vector(n)
+        # print("\nMatrix A:")
+        # print(A)
+        # print("\nVector b:")
+        # print(b)
+
+        solution_richardson = richardson(A, b, max_iterations=10000000, tol=1e-7)
+        print("Richardson Method Solution:", solution_richardson)
+
+        solution_cg, info = cg(A, b)
+        if info == 0:
+            print("SciPy Conjugate Gradient solution:", solution_cg)
+        else:
+            print("SciPy Conjugate Gradient did not converge.")
+
+        difference = np.linalg.norm(solution_richardson - solution_cg)
+        print(f"Difference between Richardson and CG solutions: {difference:.8f}")
+        
+        if difference < tolerance:
+            print("The solutions are effectively the same.")
+        else:
+            print("The solutions are different!")
+
+run_tests()

code/eigenvalue_methods.py

@@ -0,0 +1,26 @@
+from linear_algebra_utils import LinearAlgebraUtils
+
+class EigenvalueMethods:
+    @staticmethod
+    def power_method(A, max_iter=1000, tol=1e-6):
+        n = len(A)
+        x = [1] * n
+        lambda_old = 0
+        
+        for _ in range(max_iter):
+            x = LinearAlgebraUtils.matrix_vector_multiply(A, x)
+            lambda_new = LinearAlgebraUtils.vector_norm(x)
+            x = LinearAlgebraUtils.scalar_divide(x, lambda_new)
+            if abs(lambda_new - lambda_old) < tol:
+                break
+            lambda_old = lambda_new
+            
+        return lambda_new
+
+    @staticmethod
+    def inverse_power_method(A, max_iter=1000, 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)]
+        
+        return 1 / EigenvalueMethods.power_method(A_inv, max_iter, tol)

code/linear_algebra_utils.py

@@ -0,0 +1,28 @@
+class LinearAlgebraUtils:
+    @staticmethod
+    def dot_product(v1, v2):
+        return sum(x*y for x, y in zip(v1, v2))
+
+    @staticmethod
+    def matrix_vector_multiply(A, x):
+        return [LinearAlgebraUtils.dot_product(row, x) for row in A]
+
+    @staticmethod
+    def vector_subtraction(v1, v2):
+        return [x-y for x, y in zip(v1, v2)]
+
+    @staticmethod
+    def vector_addition(v1, v2):
+        return [x+y for x, y in zip(v1, v2)]
+
+    @staticmethod
+    def scalar_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 scalar_divide(x, scalar):
+        return [xi / scalar for xi in x]

code/main.py

@@ -1,87 +1,7 @@
-import unittest
+from tests import run_tests
-import numpy as np # For testing ONLY
-from richardson import modified_richardson
 
-class TestModifiedRichardson(unittest.TestCase):
+def main():
-    def setUp(self):
+    run_tests()
-        self.A_2x2 = np.random.rand(2, 2).tolist()
-        self.b_2x2 = np.random.rand(2).tolist()
-        self.x0_2x2 = np.random.rand(2).tolist()
-        self.alpha_2x2 = 0.1
 
-        self.A_3x3 = np.random.rand(3, 3).tolist()
+if __name__ == "__main__":
-        self.b_3x3 = np.random.rand(3).tolist()
+    main()
-        self.x0_3x3 = np.random.rand(3).tolist()
-        self.alpha_3x3 = 0.15
-
-    def test_convergence_2x2(self):
-        print("Testing 2x2 Convergence")
-        print(f"A: {self.A_2x2}")
-        print(f"b: {self.b_2x2}")
-        print(f"x0: {self.x0_2x2}")
-        result = modified_richardson(self.A_2x2, self.b_2x2, self.x0_2x2, self.alpha_2x2)
-        expected_solution = np.linalg.solve(np.array(self.A_2x2), np.array(self.b_2x2))
-        print(f"Result: {result}")
-        print(f"Expected: {expected_solution}")
-        for r, e in zip(result, expected_solution):
-            self.assertAlmostEqual(r, e, places=4)
-
-    def test_convergence_3x3(self):
-        print("Testing 3x3 Convergence")
-        print(f"A: {self.A_3x3}")
-        print(f"b: {self.b_3x3}")
-        print(f"x0: {self.x0_3x3}")
-        result = modified_richardson(self.A_3x3, self.b_3x3, self.x0_3x3, self.alpha_3x3)
-        expected_solution = np.linalg.solve(np.array(self.A_3x3), np.array(self.b_3x3))
-        print(f"Result: {result}")
-        print(f"Expected: {expected_solution}")
-        for r, e in zip(result, expected_solution):
-            self.assertAlmostEqual(r, e, places=2)
-
-    def test_invalid_alpha(self):
-        with self.assertRaises(ValueError):
-            modified_richardson(self.A_2x2, self.b_2x2, self.x0_2x2, alpha=-0.1)
-
-    def test_non_square_matrix(self):
-        A = [[1, 2, 3], [4, 5, 6]]  # Not a square matrix
-        b = [7, 8]
-        with self.assertRaises(ValueError):
-            modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
-
-    def test_dimension_mismatch(self):
-        b = [1, 2, 3]  # Length mismatch with A_2x2
-        with self.assertRaises(ValueError):
-            modified_richardson(self.A_2x2, b, self.x0_2x2, self.alpha_2x2)
-
-    def test_zero_matrix(self):
-        A = [[0, 0], [0, 0]]
-        b = [0, 0]
-        result = modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
-        # Solution should be [0, 0]
-        print("Testing Zero Matrix")
-        print(f"A: {A}")
-        print(f"b: {b}")
-        print(f"Result: {result}")
-        self.assertEqual(result, [0, 0])
-
-    def test_large_system(self):
-        # A large test case designed to take a long time to converge
-        size = 1000
-        A = np.random.rand(size, size) + size * np.eye(size)  # Large diagonally dominant matrix
-        b = np.random.rand(size)
-        x0 = np.random.rand(size)
-        alpha = 0.01 / size  # Small alpha to ensure convergence
-
-        print("Testing Large System")
-        #print(f"A: {A}")
-        #print(f"b: {b}")
-        #print(f"x0: {x0}")
-        result = modified_richardson(A.tolist(), b.tolist(), x0.tolist(), alpha, tol=1e-6, max_iter=500000)
-        expected_solution = np.linalg.solve(A, b)
-        print(f"Result: {result}")
-        print(f"Expected: {expected_solution}")
-        for r, e in zip(result, expected_solution):
-            self.assertAlmostEqual(r, e, places=2)
-
-if __name__ == '__main__':
-    unittest.main()

code/main_abstract.py

@@ -1,101 +0,0 @@
-import unittest
-import numpy as np # For testing ONLY
-from richardson_abstract import modified_richardson
-
-class TestModifiedRichardson(unittest.TestCase):
-    def setUp(self):
-        self.A_2x2 = np.random.rand(2, 2).tolist()
-        self.b_2x2 = np.random.rand(2).tolist()
-        self.x0_2x2 = np.random.rand(2).tolist()
-        self.alpha_2x2 = 0.1
-
-        self.A_3x3 = np.random.rand(3, 3).tolist()
-        self.b_3x3 = np.random.rand(3).tolist()
-        self.x0_3x3 = np.random.rand(3).tolist()
-        self.alpha_3x3 = 0.15
-
-    def test_convergence_2x2(self):
-        print("Testing 2x2 Convergence")
-        print(f"A: {self.A_2x2}")
-        print(f"b: {self.b_2x2}")
-        print(f"x0: {self.x0_2x2}")
-        richardson = modified_richardson(self.A_2x2, self.b_2x2, self.x0_2x2, self.alpha_2x2)
-        result = richardson()
-        # result = modified_richardson(self.A_2x2, self.b_2x2, self.x0_2x2, self.alpha_2x2)
-        expected_solution = np.linalg.solve(np.array(self.A_2x2), np.array(self.b_2x2))
-        print(f"Result: {result}")
-        print(f"Expected: {expected_solution}")
-        for r, e in zip(result, expected_solution):
-            self.assertAlmostEqual(r, e, places=4)
-
-    def test_convergence_3x3(self):
-        print("Testing 3x3 Convergence")
-        print(f"A: {self.A_3x3}")
-        print(f"b: {self.b_3x3}")
-        print(f"x0: {self.x0_3x3}")
-        richardson = modified_richardson(self.A_3x3, self.b_3x3, self.x0_3x3, self.alpha_3x3)
-        result = richardson()
-        # result = modified_richardson(self.A_3x3, self.b_3x3, self.x0_3x3, self.alpha_3x3)
-        expected_solution = np.linalg.solve(np.array(self.A_3x3), np.array(self.b_3x3))
-        print(f"Result: {result}")
-        print(f"Expected: {expected_solution}")
-        for r, e in zip(result, expected_solution):
-            self.assertAlmostEqual(r, e, places=2)
-
-    def test_invalid_alpha(self):
-        richardson = modified_richardson(self.A_2x2, self.b_2x2, self.x0_2x2, alpha=-0.1)
-        with self.assertRaises(ValueError):
-            richardson()
-            # modified_richardson(self.A_2x2, self.b_2x2, self.x0_2x2, alpha=-0.1)
-
-    def test_non_square_matrix(self):
-        A = [[1, 2, 3], [4, 5, 6]]  # Not a square matrix
-        b = [7, 8]
-        richardson = modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
-        with self.assertRaises(ValueError):
-            richardson()
-            # modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
-
-    def test_dimension_mismatch(self):
-        b = [1, 2, 3]  # Length mismatch with A_2x2
-        richardson = modified_richardson(self.A_2x2, b, self.x0_2x2, self.alpha_2x2)
-        with self.assertRaises(ValueError):
-            richardson()
-            # modified_richardson(self.A_2x2, b, self.x0_2x2, self.alpha_2x2)
-
-    def test_zero_matrix(self):
-        A = [[0, 0], [0, 0]]
-        b = [0, 0]
-        richardson = modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
-        result = richardson()
-        # result = modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
-        # Solution should be [0, 0]
-        print("Testing Zero Matrix")
-        print(f"A: {A}")
-        print(f"b: {b}")
-        print(f"Result: {result}")
-        self.assertEqual(result, [0, 0])
-
-    def test_large_system(self):
-        # A large test case designed to take a long time to converge
-        size = 10 #1000
-        A = np.random.rand(size, size) + size * np.eye(size)  # Large diagonally dominant matrix
-        b = np.random.rand(size)
-        x0 = np.random.rand(size)
-        alpha = 0.01 / size  # Small alpha to ensure convergence
-
-        print("Testing Large System")
-        #print(f"A: {A}")
-        #print(f"b: {b}")
-        #print(f"x0: {x0}")
-        richardson = modified_richardson(A.tolist(), b.tolist(), x0.tolist(), alpha, tol=1e-6, max_iter=500000)
-        result = richardson()
-        # result = modified_richardson(A.tolist(), b.tolist(), x0.tolist(), alpha, tol=1e-6, max_iter=500000)
-        expected_solution = np.linalg.solve(A, b)
-        print(f"Result: {result}")
-        print(f"Expected: {expected_solution}")
-        for r, e in zip(result, expected_solution):
-            self.assertAlmostEqual(r, e, places=2)
-
-if __name__ == '__main__':
-    unittest.main()

code/matrix_generator.py

@@ -0,0 +1,9 @@
+import numpy as np
+
+class MatrixGenerator:
+    @staticmethod
+    def generate_random_matrix_and_vector(size):
+        A = np.random.uniform(-1, 1, (size, size))
+        A = np.dot(A.T, A) + np.eye(size) * size # dodanie `size` do przekątnej zwiększa wartości własne -> obejście problemu z overflow
+        b = np.random.uniform(-1, 1, size)
+        return A, b

code/richardson.py

@@ -1,59 +0,0 @@
-def modified_richardson(A, b, x0, alpha, tol=1e-6, max_iter=1000):
-    """
-    Solves the system of linear equations Ax = b using the Modified Richardson iteration method.
-
-    Parameters:
-    A : list of lists
-        Coefficient matrix (n x n).
-    b : list
-        Right-hand side vector (n).
-    x0 : list
-        Initial guess for the solution (n).
-    alpha : float
-        Relaxation parameter (0 < alpha < 2 / max(eigenvalue(A))).
-    tol : float, optional
-        Tolerance for the stopping criterion (default is 1e-6).
-    max_iter : int, optional
-        Maximum number of iterations (default is 1000).
-
-    Returns:
-    x : list
-        Approximate solution to the system of equations.
-    """
-    n = len(A)
-    x = x0[:]
-
-    if len(A) != len(A[0]):
-        raise ValueError("Matrix A must be square.")
-    if len(b) != n:
-        raise ValueError("Dimension mismatch between A and b.")
-    if alpha <= 0:
-        raise ValueError("Alpha must be greater than 0.")
-
-    def vector_norm(v):
-        return sum(vi ** 2 for vi in v) ** 0.5
-
-    def mat_vec_mult(mat, vec):
-        return [sum(mat[i][j] * vec[j] for j in range(len(vec))) for i in range(len(mat))]
-
-    def vec_sub(v1, v2):
-        return [v1[i] - v2[i] for i in range(len(v1))]
-
-    def vec_add(v1, v2):
-        return [v1[i] + v2[i] for i in range(len(v1))]
-
-    def vec_scalar_mult(scalar, vec):
-        return [scalar * vi for vi in vec]
-
-    r = vec_sub(b, mat_vec_mult(A, x))
-    iteration = 0
-
-    while vector_norm(r) > tol and iteration < max_iter:
-        x = vec_add(x, vec_scalar_mult(alpha, r))
-        r = vec_sub(b, mat_vec_mult(A, x))
-        iteration += 1
-
-    if iteration == max_iter:
-        raise ValueError("Maximum number of iterations reached before convergence")
-
-    return x

code/richardson_abstract.py

@@ -1,128 +0,0 @@
-from abc import ABC, abstractmethod
-import threading
-
-class modified_richardson_base(ABC):
-    """
-    Solves the system of linear equations Ax = b using the Modified Richardson iteration method.
-
-    Parameters:
-    A : list of lists
-        Coefficient matrix (n x n).
-    b : list
-        Right-hand side vector (n).
-    x0 : list
-        Initial guess for the solution (n).
-    alpha : float
-        Relaxation parameter (0 < alpha < 2 / max(eigenvalue(A))).
-    tol : float, optional
-        Tolerance for the stopping criterion (default is 1e-6).
-    max_iter : int, optional
-        Maximum number of iterations (default is 1000).
-
-    Returns:
-    x : list
-        Approximate solution to the system of equations.
-    """
-    def __init__(self, A, b, x0, alpha, tol=1e-6, max_iter=1000):
-        self.A = A
-        self.b = b
-        self.x0 = x0
-        self.alpha = alpha
-        self.tol = tol
-        self.max_iter = max_iter
-        self.n = len(A)
-        self.x = self.x0[:]
-
-    def check_input_data(self):
-        if len(self.A) != len(self.A[0]):
-            raise ValueError("Matrix A must be square.")
-        if len(self.b) != self.n:
-            raise ValueError("Dimension mismatch between A and b.")
-        if self.alpha <= 0:
-            raise ValueError("Alpha must be greater than 0.")
-
-    @abstractmethod
-    def vector_norm(self, v):
-        pass
-
-    @abstractmethod
-    def mat_vec_mult(self, mat, vec):
-        pass
-
-    @abstractmethod
-    def vec_sub(self, v1, v2):
-        pass
-
-    @abstractmethod
-    def vec_add(self, v1, v2):
-        pass
-
-    @abstractmethod
-    def vec_scalar_mult(self, scalar, vec):
-        pass
-
-    def __call__(self):
-        self.check_input_data()
-        x = self.x
-        
-        r = self.vec_sub(self.b, self.mat_vec_mult(self.A, x))
-        iteration = 0
-
-        while self.vector_norm(r) > self.tol and iteration < self.max_iter:
-            x = self.vec_add(x, self.vec_scalar_mult(self.alpha, r))
-            r = self.vec_sub(self.b, self.mat_vec_mult(self.A, x))
-            iteration += 1
-
-        if iteration == self.max_iter:
-            raise ValueError("Maximum number of iterations reached before convergence")
-
-        return x
-
-
-class modified_richardson(modified_richardson_base):
-    def vector_norm(self, v):
-        return sum(vi ** 2 for vi in v) ** 0.5
-
-    def mat_vec_mult(self, mat, vec):
-        return [sum(mat[i][j] * vec[j] for j in range(len(vec))) for i in range(len(mat))]
-
-    def vec_sub(self, v1, v2):
-        return [v1[i] - v2[i] for i in range(len(v1))]
-
-    def vec_add(self, v1, v2):
-        return [v1[i] + v2[i] for i in range(len(v1))]
-
-    def vec_scalar_mult(self, scalar, vec):
-        return [scalar * vi for vi in vec]
-
-    
-class modified_richardson_with_threads(modified_richardson_base):
-    def compute_with_threads(self, v1, v2, function: function):
-        """ Executes the provided function on many threads """
-        result = [0] * len(v1)
-        threads = []
-
-        for i in range(len(v1)):
-            t = threading.Thread(target=function, args=(v1, v2, result, i))
-            threads.append(t)
-            t.start()
-
-        for t in threads:
-            t.join()
-
-        return result
-    
-    def subtract_elements(self, v1, v2, result, index):
-        """ This function is executed by single thread. It calculates one row in matrix """
-        result[index] = v1[index] - v2[index]
-
-    def add_elements(self, v1, v2, result, index):
-        """ As above """
-        result[index] = v1[index] - v2[index]
-
-    def vec_sub(self, v1, v2):
-        return self.compute_with_threads(v1, v2, self.subtract_elements)
-
-    def vec_add(self, v1, v2):
-        return self.compute_with_threads(v1, v2, self.add_elements)
-    

code/richardson_method.py

@@ -0,0 +1,33 @@
+from linear_algebra_utils import LinearAlgebraUtils
+from eigenvalue_methods import EigenvalueMethods
+
+class RichardsonMethod:
+    def __init__(self, A, b, x0=None, max_iterations=1000, 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
+
+    def solve(self):
+        x = self.x0[:]
+        lambda_min = EigenvalueMethods.inverse_power_method(self.A)
+        lambda_max = EigenvalueMethods.power_method(self.A)
+
+        if lambda_min < 0:
+            raise ValueError("Matrix A is not positive semi-definite.")
+        
+        omega = 2 / (lambda_min + lambda_max)
+
+        for iteration in range(self.max_iterations):
+            Ax = LinearAlgebraUtils.matrix_vector_multiply(self.A, x)
+            residual = LinearAlgebraUtils.vector_subtraction(self.b, Ax)
+
+            if LinearAlgebraUtils.vector_norm(residual) < self.tol:
+                print('Convergence achieved.')
+                return x
+            
+            x = LinearAlgebraUtils.vector_addition(x, LinearAlgebraUtils.scalar_multiply(omega, residual))
+
+        print('Maximum number of iterations reached without convergence.')
+        return x

code/tests.py

@@ -0,0 +1,31 @@
+import numpy as np
+from scipy.sparse.linalg import cg
+from matrix_generator import MatrixGenerator
+from richardson_method import RichardsonMethod
+
+def run_tests():
+    test_sizes = [2, 3, 4, 5, 10, 20, 50, 100]
+    tolerance = 1e-5
+    
+    for n in test_sizes:
+        print(f"\nRunning test for n = {n}")
+        
+        A, b = MatrixGenerator.generate_random_matrix_and_vector(n)
+
+        richardson_solver = RichardsonMethod(A, b, max_iterations=10000000, tol=1e-7)
+        solution_richardson = richardson_solver.solve()
+        print("Richardson Method Solution:", solution_richardson)
+
+        solution_cg, info = cg(A, b)
+        if info == 0:
+            print("SciPy Conjugate Gradient solution:", solution_cg)
+        else:
+            print("SciPy Conjugate Gradient did not converge.")
+
+        difference = np.linalg.norm(solution_richardson - solution_cg)
+        print(f"Difference between Richardson and CG solutions: {difference:.8f}")
+        
+        if difference < tolerance:
+            print("The solutions are effectively the same.")
+        else:
+            print("The solutions are different!")