From 52d6012092c4cb1390656df87637d65ffce6c71b Mon Sep 17 00:00:00 2001 From: aleksandrasob Date: Sun, 20 Oct 2024 16:30:42 +0200 Subject: [PATCH] Add sequential version of Richardson's algorithm and first tests --- code/Richardson_numpy.py | 70 +++++++++++++++++++ code/eigenvalue_methods.py | 26 +++++++ code/linear_algebra_utils.py | 28 ++++++++ code/main.py | 90 ++---------------------- code/main_abstract.py | 101 --------------------------- code/matrix_generator.py | 9 +++ code/richardson.py | 59 ---------------- code/richardson_abstract.py | 128 ----------------------------------- code/richardson_method.py | 33 +++++++++ code/tests.py | 31 +++++++++ 10 files changed, 202 insertions(+), 373 deletions(-) create mode 100644 code/Richardson_numpy.py create mode 100644 code/eigenvalue_methods.py create mode 100644 code/linear_algebra_utils.py delete mode 100644 code/main_abstract.py create mode 100644 code/matrix_generator.py delete mode 100644 code/richardson.py delete mode 100644 code/richardson_abstract.py create mode 100644 code/richardson_method.py create mode 100644 code/tests.py diff --git a/code/Richardson_numpy.py b/code/Richardson_numpy.py new file mode 100644 index 00000000..9de5a394 --- /dev/null +++ b/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() diff --git a/code/eigenvalue_methods.py b/code/eigenvalue_methods.py new file mode 100644 index 00000000..5f218d7d --- /dev/null +++ b/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) diff --git a/code/linear_algebra_utils.py b/code/linear_algebra_utils.py new file mode 100644 index 00000000..9b3647fa --- /dev/null +++ b/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] diff --git a/code/main.py b/code/main.py index ee207c61..16977fa9 100644 --- a/code/main.py +++ b/code/main.py @@ -1,87 +1,7 @@ -import unittest -import numpy as np # For testing ONLY -from richardson import modified_richardson +from tests import run_tests -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 +def main(): + run_tests() - 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}") - 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() +if __name__ == "__main__": + main() diff --git a/code/main_abstract.py b/code/main_abstract.py deleted file mode 100644 index 3398ee28..00000000 --- a/code/main_abstract.py +++ /dev/null @@ -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() diff --git a/code/matrix_generator.py b/code/matrix_generator.py new file mode 100644 index 00000000..ab4addd2 --- /dev/null +++ b/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 diff --git a/code/richardson.py b/code/richardson.py deleted file mode 100644 index a94b48e2..00000000 --- a/code/richardson.py +++ /dev/null @@ -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 \ No newline at end of file diff --git a/code/richardson_abstract.py b/code/richardson_abstract.py deleted file mode 100644 index 0cbaa03c..00000000 --- a/code/richardson_abstract.py +++ /dev/null @@ -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) - \ No newline at end of file diff --git a/code/richardson_method.py b/code/richardson_method.py new file mode 100644 index 00000000..86d07927 --- /dev/null +++ b/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 diff --git a/code/tests.py b/code/tests.py new file mode 100644 index 00000000..dcc6b654 --- /dev/null +++ b/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!")