From 52d6012092c4cb1390656df87637d65ffce6c71b Mon Sep 17 00:00:00 2001 From: aleksandrasob Date: Sun, 20 Oct 2024 16:30:42 +0200 Subject: [PATCH 01/10] 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!") From d3b0e8e0cc420f93241711ab0c3dd9a804df041e Mon Sep 17 00:00:00 2001 From: Krzysztof Rudnicki Date: Sun, 20 Oct 2024 17:47:37 +0200 Subject: [PATCH 02/10] feat: convergence method suggested by wikipedia --- code/eigenvalue_methods.py | 2 +- code/linear_algebra_utils.py | 22 ++++++++++++++++++---- code/matrix_generator.py | 3 +++ code/richardson_method.py | 13 ++++++++----- code/tests.py | 2 +- 5 files changed, 31 insertions(+), 11 deletions(-) diff --git a/code/eigenvalue_methods.py b/code/eigenvalue_methods.py index 5f218d7d..8c7dc2c6 100644 --- a/code/eigenvalue_methods.py +++ b/code/eigenvalue_methods.py @@ -10,7 +10,7 @@ class EigenvalueMethods: 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) + x = LinearAlgebraUtils.vector_scalar_divide(x, lambda_new) if abs(lambda_new - lambda_old) < tol: break lambda_old = lambda_new diff --git a/code/linear_algebra_utils.py b/code/linear_algebra_utils.py index 9b3647fa..c6175668 100644 --- a/code/linear_algebra_utils.py +++ b/code/linear_algebra_utils.py @@ -1,3 +1,5 @@ +import math + class LinearAlgebraUtils: @staticmethod def dot_product(v1, v2): @@ -8,15 +10,19 @@ class LinearAlgebraUtils: return [LinearAlgebraUtils.dot_product(row, x) for row in A] @staticmethod - def vector_subtraction(v1, v2): + def matrix_scalar_multiply(A, w): + return [[w * A[i][j] for j in range(len(A[0]))] for i in range(len(A))] + + @staticmethod + def vector_vector_subtraction(v1, v2): return [x-y for x, y in zip(v1, v2)] @staticmethod - def vector_addition(v1, v2): + def vector_vector_addition(v1, v2): return [x+y for x, y in zip(v1, v2)] @staticmethod - def scalar_multiply(omega, vector): + def scalar_matrix_multiply(omega, vector): return [omega * x for x in vector] @staticmethod @@ -24,5 +30,13 @@ class LinearAlgebraUtils: return sum(x*x for x in v)**0.5 @staticmethod - def scalar_divide(x, scalar): + 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))] diff --git a/code/matrix_generator.py b/code/matrix_generator.py index ab4addd2..0b99273b 100644 --- a/code/matrix_generator.py +++ b/code/matrix_generator.py @@ -7,3 +7,6 @@ class MatrixGenerator: 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 generate_identity_matrix(size): + return np.eye(size) diff --git a/code/richardson_method.py b/code/richardson_method.py index 86d07927..0dfc753b 100644 --- a/code/richardson_method.py +++ b/code/richardson_method.py @@ -1,13 +1,15 @@ from linear_algebra_utils import LinearAlgebraUtils from eigenvalue_methods import EigenvalueMethods +from matrix_generator import MatrixGenerator class RichardsonMethod: - def __init__(self, A, b, x0=None, max_iterations=1000, tol=1e-5): + def __init__(self, A, b, size: int, 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 + self.I = MatrixGenerator.generate_identity_matrix(size) def solve(self): x = self.x0[:] @@ -21,13 +23,14 @@ class RichardsonMethod: 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: + residual = LinearAlgebraUtils.vector_vector_subtraction(self.b, Ax) + wA = LinearAlgebraUtils.matrix_scalar_multiply(self.A, omega) + IMinuswA = LinearAlgebraUtils.matrix_matrix_subtraction(self.I, wA) + if LinearAlgebraUtils.matrix_norm(IMinuswA) < 1: print('Convergence achieved.') return x - x = LinearAlgebraUtils.vector_addition(x, LinearAlgebraUtils.scalar_multiply(omega, residual)) + x = LinearAlgebraUtils.vector_vector_addition(x, LinearAlgebraUtils.scalar_matrix_multiply(omega, residual)) print('Maximum number of iterations reached without convergence.') return x diff --git a/code/tests.py b/code/tests.py index dcc6b654..0dc23ccc 100644 --- a/code/tests.py +++ b/code/tests.py @@ -12,7 +12,7 @@ def run_tests(): A, b = MatrixGenerator.generate_random_matrix_and_vector(n) - richardson_solver = RichardsonMethod(A, b, max_iterations=10000000, tol=1e-7) + richardson_solver = RichardsonMethod(A, b, size=n, max_iterations=10000000, tol=1e-7) solution_richardson = richardson_solver.solve() print("Richardson Method Solution:", solution_richardson) From 76ace2b76c8fae65be8df4f91bc21d657eb5431c Mon Sep 17 00:00:00 2001 From: Krzysztof Rudnicki Date: Sun, 20 Oct 2024 18:09:06 +0200 Subject: [PATCH 03/10] feat: reworked convergence calculations --- code/Richardson_numpy.py | 70 --------------------------------------- code/richardson_method.py | 28 ++++++++-------- code/tests.py | 5 +-- 3 files changed, 16 insertions(+), 87 deletions(-) delete mode 100644 code/Richardson_numpy.py diff --git a/code/Richardson_numpy.py b/code/Richardson_numpy.py deleted file mode 100644 index 9de5a394..00000000 --- a/code/Richardson_numpy.py +++ /dev/null @@ -1,70 +0,0 @@ -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/richardson_method.py b/code/richardson_method.py index 0dfc753b..59aa2964 100644 --- a/code/richardson_method.py +++ b/code/richardson_method.py @@ -10,27 +10,25 @@ class RichardsonMethod: self.max_iterations = max_iterations self.tol = tol self.I = MatrixGenerator.generate_identity_matrix(size) + self.lambda_min = EigenvalueMethods.inverse_power_method(self.A) + if self.lambda_min < 0: + raise ValueError("Matrix A is not positive semi-definite.") + self.lambda_max = EigenvalueMethods.power_method(self.A) + self.omega = 2 / (self.lambda_min + self.lambda_max) + + def will_converge(self) -> bool: + wA = LinearAlgebraUtils.matrix_scalar_multiply(self.A, self.omega) + IMinuswA = LinearAlgebraUtils.matrix_matrix_subtraction(self.I, wA) + return LinearAlgebraUtils.matrix_norm(IMinuswA) < 1 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) - + if not self.will_converge(): + print("Richardson method for those values will NOT converge") for iteration in range(self.max_iterations): Ax = LinearAlgebraUtils.matrix_vector_multiply(self.A, x) residual = LinearAlgebraUtils.vector_vector_subtraction(self.b, Ax) - wA = LinearAlgebraUtils.matrix_scalar_multiply(self.A, omega) - IMinuswA = LinearAlgebraUtils.matrix_matrix_subtraction(self.I, wA) - if LinearAlgebraUtils.matrix_norm(IMinuswA) < 1: - print('Convergence achieved.') - return x - - x = LinearAlgebraUtils.vector_vector_addition(x, LinearAlgebraUtils.scalar_matrix_multiply(omega, residual)) + x = LinearAlgebraUtils.vector_vector_addition(x, LinearAlgebraUtils.scalar_matrix_multiply(self.omega, residual)) print('Maximum number of iterations reached without convergence.') return x diff --git a/code/tests.py b/code/tests.py index 0dc23ccc..80852a32 100644 --- a/code/tests.py +++ b/code/tests.py @@ -11,8 +11,9 @@ def run_tests(): print(f"\nRunning test for n = {n}") A, b = MatrixGenerator.generate_random_matrix_and_vector(n) - - richardson_solver = RichardsonMethod(A, b, size=n, max_iterations=10000000, tol=1e-7) + print("A: ", A) + print("b: ", b) + richardson_solver = RichardsonMethod(A, b, size=n, max_iterations=1000, tol=1e-7) solution_richardson = richardson_solver.solve() print("Richardson Method Solution:", solution_richardson) From ea98dc9712ac7045cb0d98f6695701af54dd255a Mon Sep 17 00:00:00 2001 From: Krzysztof Rudnicki Date: Sun, 20 Oct 2024 18:27:51 +0200 Subject: [PATCH 04/10] feat: changed to proper pytests --- code/main.py | 8 +++---- code/richardson_method.py | 8 ++++--- code/tests.py | 48 +++++++++++++++++++++------------------ 3 files changed, 34 insertions(+), 30 deletions(-) diff --git a/code/main.py b/code/main.py index 16977fa9..b468d54c 100644 --- a/code/main.py +++ b/code/main.py @@ -1,7 +1,5 @@ -from tests import run_tests - -def main(): - run_tests() +import pytest if __name__ == "__main__": - main() + # Run pytest and exit with the appropriate status code + pytest.main(["-v", "tests.py"]) diff --git a/code/richardson_method.py b/code/richardson_method.py index 59aa2964..256c2a9c 100644 --- a/code/richardson_method.py +++ b/code/richardson_method.py @@ -15,20 +15,22 @@ class RichardsonMethod: raise ValueError("Matrix A is not positive semi-definite.") self.lambda_max = EigenvalueMethods.power_method(self.A) self.omega = 2 / (self.lambda_min + self.lambda_max) + def will_converge(self) -> bool: wA = LinearAlgebraUtils.matrix_scalar_multiply(self.A, self.omega) IMinuswA = LinearAlgebraUtils.matrix_matrix_subtraction(self.I, wA) - return LinearAlgebraUtils.matrix_norm(IMinuswA) < 1 + norm = LinearAlgebraUtils.matrix_norm(IMinuswA) + return norm < 1 def solve(self): x = self.x0[:] if not self.will_converge(): - print("Richardson method for those values will NOT converge") + return "Richardson method for those values will NOT converge" + for iteration in range(self.max_iterations): Ax = LinearAlgebraUtils.matrix_vector_multiply(self.A, x) residual = LinearAlgebraUtils.vector_vector_subtraction(self.b, Ax) x = LinearAlgebraUtils.vector_vector_addition(x, LinearAlgebraUtils.scalar_matrix_multiply(self.omega, residual)) - print('Maximum number of iterations reached without convergence.') return x diff --git a/code/tests.py b/code/tests.py index 80852a32..3b3a89b9 100644 --- a/code/tests.py +++ b/code/tests.py @@ -1,32 +1,36 @@ +import pytest 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] +@pytest.mark.parametrize("n", [2, 3, 4, 5, 10, 20, 50, 100]) +def test_richardson_vs_cg(n): tolerance = 1e-5 + A, b = MatrixGenerator.generate_random_matrix_and_vector(n) - for n in test_sizes: - print(f"\nRunning test for n = {n}") - - A, b = MatrixGenerator.generate_random_matrix_and_vector(n) - print("A: ", A) - print("b: ", b) - richardson_solver = RichardsonMethod(A, b, size=n, max_iterations=1000, 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.") + richardson_solver = RichardsonMethod(A, b, size=n, max_iterations=1000, tol=1e-7) + solution_richardson = richardson_solver.solve() + + solution_cg, info = cg(A, b) + + if info == 0: # SciPy CG converged + assert_scipy_converged(solution_richardson, solution_cg, tolerance) + else: # SciPy CG did not converge + assert_scipy_not_converged(solution_richardson) +def assert_scipy_converged(solution_richardson, solution_cg, tolerance): + if solution_richardson == "Richardson method for those values will NOT converge": + print("Richardson did not converge, while SciPy did") + assert False, "Richardson did not converge, while SciPy did" + else: 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!") + assert difference < tolerance, f"The solutions are different! Difference: {difference:.8f}" + +def assert_scipy_not_converged(solution_richardson): + if solution_richardson == "Richardson method for those values will NOT converge": + print("Richardson and SciPy did not converge") + else: + print("Richardson converged while SciPy did not:", solution_richardson) + assert False, "Richardson converged while SciPy did not" From 838d5ed56303a72eb00f5941ca2f82631cf078fb Mon Sep 17 00:00:00 2001 From: Krzysztof Rudnicki Date: Sun, 20 Oct 2024 18:32:12 +0200 Subject: [PATCH 05/10] fix: removed adding "size" to self values --- code/matrix_generator.py | 1 - code/tests.py | 15 +++++++++++---- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/code/matrix_generator.py b/code/matrix_generator.py index 0b99273b..24ad3355 100644 --- a/code/matrix_generator.py +++ b/code/matrix_generator.py @@ -4,7 +4,6 @@ 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/tests.py b/code/tests.py index 3b3a89b9..949a8d32 100644 --- a/code/tests.py +++ b/code/tests.py @@ -15,22 +15,29 @@ def test_richardson_vs_cg(n): solution_cg, info = cg(A, b) if info == 0: # SciPy CG converged - assert_scipy_converged(solution_richardson, solution_cg, tolerance) + assert_scipy_converged(solution_richardson, solution_cg, tolerance, A, b) else: # SciPy CG did not converge - assert_scipy_not_converged(solution_richardson) + assert_scipy_not_converged(solution_richardson, A, b) -def assert_scipy_converged(solution_richardson, solution_cg, tolerance): +def assert_scipy_converged(solution_richardson, solution_cg, tolerance, A, b): if solution_richardson == "Richardson method for those values will NOT converge": print("Richardson did not converge, while SciPy did") + print("Matrix A:\n", A) + print("Vector b:\n", b) assert False, "Richardson did not converge, while SciPy did" else: difference = np.linalg.norm(solution_richardson - solution_cg) print(f"Difference between Richardson and CG solutions: {difference:.8f}") + if difference >= tolerance: + print("Matrix A:\n", A) + print("Vector b:\n", b) assert difference < tolerance, f"The solutions are different! Difference: {difference:.8f}" -def assert_scipy_not_converged(solution_richardson): +def assert_scipy_not_converged(solution_richardson, A, b): if solution_richardson == "Richardson method for those values will NOT converge": print("Richardson and SciPy did not converge") else: print("Richardson converged while SciPy did not:", solution_richardson) + print("Matrix A:\n", A) + print("Vector b:\n", b) assert False, "Richardson converged while SciPy did not" From 9014c37064441494be6f54ece5ccb17abf672f59 Mon Sep 17 00:00:00 2001 From: Krzysztof Rudnicki Date: Sun, 20 Oct 2024 18:39:01 +0200 Subject: [PATCH 06/10] feat: can run tests independently --- code/tests.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/code/tests.py b/code/tests.py index 949a8d32..8b3fd978 100644 --- a/code/tests.py +++ b/code/tests.py @@ -5,7 +5,7 @@ from matrix_generator import MatrixGenerator from richardson_method import RichardsonMethod @pytest.mark.parametrize("n", [2, 3, 4, 5, 10, 20, 50, 100]) -def test_richardson_vs_cg(n): +def test_richardson_vs_cg(n: int): tolerance = 1e-5 A, b = MatrixGenerator.generate_random_matrix_and_vector(n) @@ -28,7 +28,11 @@ def assert_scipy_converged(solution_richardson, solution_cg, tolerance, A, b): else: difference = np.linalg.norm(solution_richardson - solution_cg) print(f"Difference between Richardson and CG solutions: {difference:.8f}") - if difference >= tolerance: + if difference < tolerance: + print("Both Richardson and CG converged and calculated correct values.") + print("Solution CG:\n", solution_cg) + print("Solution Richardson:\n", solution_richardson) + else: print("Matrix A:\n", A) print("Vector b:\n", b) assert difference < tolerance, f"The solutions are different! Difference: {difference:.8f}" @@ -41,3 +45,8 @@ def assert_scipy_not_converged(solution_richardson, A, b): print("Matrix A:\n", A) print("Vector b:\n", b) assert False, "Richardson converged while SciPy did not" + +if __name__ == "__main__": + # Run pytest and exit with the appropriate status code + for n in [2, 3, 4, 5, 10, 20, 50, 100]: + test_richardson_vs_cg(n) \ No newline at end of file From d74da05da87de4f0ecd36111ab064bc0b21c68c6 Mon Sep 17 00:00:00 2001 From: Krzysztof Rudnicki Date: Sun, 20 Oct 2024 19:06:15 +0200 Subject: [PATCH 07/10] feat: made max iterations obligatory argument --- code/eigenvalue_methods.py | 4 ++-- code/linear_algebra_utils.py | 15 ++++++++------- code/matrix_generator.py | 17 ++++++++++++++++- code/richardson_method.py | 25 ++++++++++++++++--------- code/tests.py | 29 ++++++++++++++++++++++++++++- 5 files changed, 70 insertions(+), 20 deletions(-) diff --git a/code/eigenvalue_methods.py b/code/eigenvalue_methods.py index 8c7dc2c6..00302d42 100644 --- a/code/eigenvalue_methods.py +++ b/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)] diff --git a/code/linear_algebra_utils.py b/code/linear_algebra_utils.py index c6175668..4e08ee6c 100644 --- a/code/linear_algebra_utils.py +++ b/code/linear_algebra_utils.py @@ -8,6 +8,14 @@ class LinearAlgebraUtils: @staticmethod 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): @@ -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))] diff --git a/code/matrix_generator.py b/code/matrix_generator.py index 24ad3355..7451062a 100644 --- a/code/matrix_generator.py +++ b/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 diff --git a/code/richardson_method.py b/code/richardson_method.py index 256c2a9c..88ae9890 100644 --- a/code/richardson_method.py +++ 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 = 2 / (self.lambda_min + self.lambda_max) + self.omega = RichardsonMethod.calculate_omega(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) + + @staticmethod + def calculate_omega(lambda_min, lambda_max): + return 2 / (lambda_min + lambda_max) - def will_converge(self) -> bool: - wA = LinearAlgebraUtils.matrix_scalar_multiply(self.A, self.omega) - IMinuswA = LinearAlgebraUtils.matrix_matrix_subtraction(self.I, wA) + @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): diff --git a/code/tests.py b/code/tests.py index 8b3fd978..e3defbbd 100644 --- a/code/tests.py +++ b/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) From 70c0da16c086171d358780849f767d6c6472160c Mon Sep 17 00:00:00 2001 From: Krzysztof Rudnicki Date: Sun, 20 Oct 2024 19:19:31 +0200 Subject: [PATCH 08/10] fix: calculations of inverse matrix --- code/eigenvalue_methods.py | 5 +++-- code/linear_algebra_utils.py | 28 ++++++++++++++++++++++++++++ 2 files changed, 31 insertions(+), 2 deletions(-) diff --git a/code/eigenvalue_methods.py b/code/eigenvalue_methods.py index 00302d42..0d97defa 100644 --- a/code/eigenvalue_methods.py +++ b/code/eigenvalue_methods.py @@ -21,6 +21,7 @@ class EigenvalueMethods: 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)] - + A_inv = [LinearAlgebraUtils.gaussian_elimination(A.tolist(), I_col) for I_col in I] + A_inv = list(map(list, zip(*A_inv))) + return 1 / EigenvalueMethods.power_method(A_inv, max_iter, tol) diff --git a/code/linear_algebra_utils.py b/code/linear_algebra_utils.py index 4e08ee6c..5ff15f83 100644 --- a/code/linear_algebra_utils.py +++ b/code/linear_algebra_utils.py @@ -41,3 +41,31 @@ class LinearAlgebraUtils: @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))] + + @staticmethod + def gaussian_elimination(A, b): + n = len(A) + M = [row[:] for row in A] + + for i in range(n): + M[i].append(b[i]) + + for k in range(n): + if M[k][k] == 0: + for i in range(k + 1, n): + if M[i][k] != 0: + M[k], M[i] = M[i], M[k] + break + + for i in range(k + 1, n): + factor = M[i][k] / M[k][k] + for j in range(k, n + 1): + M[i][j] -= factor * M[k][j] + + x = [0] * n + for i in range(n - 1, -1, -1): + x[i] = M[i][-1] / M[i][i] + for k in range(i - 1, -1, -1): + M[k][-1] -= M[k][i] * x[i] + + return x From cb179a76ce2b367f03180d03847457efeb32abf3 Mon Sep 17 00:00:00 2001 From: Krzysztof Rudnicki Date: Sun, 20 Oct 2024 19:20:38 +0200 Subject: [PATCH 09/10] chore: moved debug prints to separate function --- code/tests.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/code/tests.py b/code/tests.py index e3defbbd..6e16b4d4 100644 --- a/code/tests.py +++ b/code/tests.py @@ -23,19 +23,20 @@ def calculate_eigenvalues(A): 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) +def numpy_additional_richardson_calculations(): 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)) - + +@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) richardson_solver = RichardsonMethod(A, b, max_iterations, size=n, tol=1e-7) solution_richardson = richardson_solver.solve() From 6946035518a5cbaae21e990e5c6ff9355f7ef1cd Mon Sep 17 00:00:00 2001 From: Krzysztof Rudnicki Date: Sun, 20 Oct 2024 19:34:12 +0200 Subject: [PATCH 10/10] feat: compare norms calculated by richardson vs numpy --- code/richardson_method.py | 4 ++-- code/tests.py | 26 +++++++++++++------------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/code/richardson_method.py b/code/richardson_method.py index 88ae9890..2d963d11 100644 --- a/code/richardson_method.py +++ b/code/richardson_method.py @@ -33,11 +33,11 @@ class RichardsonMethod: def solve(self): x = self.x0[:] if RichardsonMethod.convergence_norm(self.A, self.omega, self.I) >= 1: - return "Richardson method for those values will NOT converge" + return RichardsonMethod.convergence_norm(self.A, self.omega, self.I), "Richardson method for those values will NOT converge", for iteration in range(self.max_iterations): Ax = LinearAlgebraUtils.matrix_vector_multiply(self.A, x) residual = LinearAlgebraUtils.vector_vector_subtraction(self.b, Ax) x = LinearAlgebraUtils.vector_vector_addition(x, LinearAlgebraUtils.scalar_matrix_multiply(self.omega, residual)) - return x + return x, 0 diff --git a/code/tests.py b/code/tests.py index 6e16b4d4..1bdf7961 100644 --- a/code/tests.py +++ b/code/tests.py @@ -23,13 +23,13 @@ def calculate_eigenvalues(A): return lambda_min, lambda_max -def numpy_additional_richardson_calculations(): +def calcualte_norm_from_matrix_numpy(A, n, max_iterations): 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)) + return calculate_norm_numpy(I, omega, A) + + @pytest.mark.parametrize("n", [2, 3, 4, 5, 10, 20, 50, 100]) def test_richardson_vs_cg(n: int): @@ -38,20 +38,20 @@ def test_richardson_vs_cg(n: int): max_iterations=1000 A, b = MatrixGenerator.generate_random_matrix_and_vector(n) richardson_solver = RichardsonMethod(A, b, max_iterations, size=n, tol=1e-7) - solution_richardson = richardson_solver.solve() + solution_richardson, info_richardson = richardson_solver.solve() solution_cg, info = cg(A, b) if info == 0: # SciPy CG converged - assert_scipy_converged(solution_richardson, solution_cg, tolerance, A, b) + assert_scipy_converged(solution_richardson, info_richardson, solution_cg, tolerance, A, b, max_iterations, n) else: # SciPy CG did not converge - assert_scipy_not_converged(solution_richardson, A, b) + assert_scipy_not_converged(solution_richardson, info_richardson, A, b) -def assert_scipy_converged(solution_richardson, solution_cg, tolerance, A, b): - if solution_richardson == "Richardson method for those values will NOT converge": +def assert_scipy_converged(solution_richardson, info_richardson, solution_cg, tolerance, A, b, max_iterations, n): + if info_richardson == "Richardson method for those values will NOT converge": print("Richardson did not converge, while SciPy did") - print("Matrix A:\n", A) - print("Vector b:\n", b) + numpy_norm = calcualte_norm_from_matrix_numpy(A, n, max_iterations) + print("Numpy norm: ", numpy_norm, " Richardson norm: ", solution_richardson) assert False, "Richardson did not converge, while SciPy did" else: difference = np.linalg.norm(solution_richardson - solution_cg) @@ -65,8 +65,8 @@ def assert_scipy_converged(solution_richardson, solution_cg, tolerance, A, b): print("Vector b:\n", b) assert difference < tolerance, f"The solutions are different! Difference: {difference:.8f}" -def assert_scipy_not_converged(solution_richardson, A, b): - if solution_richardson == "Richardson method for those values will NOT converge": +def assert_scipy_not_converged(solution_richardson, info_richardson, A, b): + if info_richardson == "Richardson method for those values will NOT converge": print("Richardson and SciPy did not converge") else: print("Richardson converged while SciPy did not:", solution_richardson)