diff --git a/code/eigenvalue_methods.py b/code/eigenvalue_methods.py new file mode 100644 index 00000000..0d97defa --- /dev/null +++ b/code/eigenvalue_methods.py @@ -0,0 +1,27 @@ +from linear_algebra_utils import LinearAlgebraUtils + +class EigenvalueMethods: + @staticmethod + def power_method(A, max_iter, 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.vector_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, 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 = [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 new file mode 100644 index 00000000..5ff15f83 --- /dev/null +++ b/code/linear_algebra_utils.py @@ -0,0 +1,71 @@ +import math + +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_norm(v): + return sum(x*x for x in v)**0.5 + + @staticmethod + def vector_scalar_divide(x, scalar): + return [xi / scalar for xi in x] + + @staticmethod + def matrix_scalar_multiply(A, w): + return [[w * A[i][j] for j in range(len(A[0]))] for i in range(len(A))] + + @staticmethod + def vector_vector_subtraction(v1, v2): + return [x-y for x, y in zip(v1, v2)] + + @staticmethod + def vector_vector_addition(v1, v2): + return [x+y for x, y in zip(v1, v2)] + + @staticmethod + def scalar_matrix_multiply(omega, vector): + return [omega * x for x in vector] + + + @staticmethod + def matrix_norm(A): + return math.sqrt(sum(sum(element ** 2 for element in row) for row in A)) + + @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 diff --git a/code/main.py b/code/main.py index ee207c61..b468d54c 100644 --- a/code/main.py +++ b/code/main.py @@ -1,87 +1,5 @@ -import unittest -import numpy as np # For testing ONLY -from richardson import modified_richardson +import pytest -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}") - 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__": + # Run pytest and exit with the appropriate status code + pytest.main(["-v", "tests.py"]) 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..7451062a --- /dev/null +++ b/code/matrix_generator.py @@ -0,0 +1,26 @@ +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 = MatrixGenerator.generate_spd_matrix(size) + b = np.random.uniform(-1, 1, size) + return A, b + + def generate_identity_matrix(size): + return np.eye(size) 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..2d963d11 --- /dev/null +++ b/code/richardson_method.py @@ -0,0 +1,43 @@ +from linear_algebra_utils import LinearAlgebraUtils +from eigenvalue_methods import EigenvalueMethods +from matrix_generator import MatrixGenerator + +class RichardsonMethod: + 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, 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.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) + + @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 + + def solve(self): + x = self.x0[:] + if RichardsonMethod.convergence_norm(self.A, self.omega, self.I) >= 1: + 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, 0 diff --git a/code/tests.py b/code/tests.py new file mode 100644 index 00000000..1bdf7961 --- /dev/null +++ b/code/tests.py @@ -0,0 +1,80 @@ +import pytest +import numpy as np +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 + +def calcualte_norm_from_matrix_numpy(A, n, max_iterations): + lambda_min, lambda_max = calculate_eigenvalues(A) + omega = 2 / (lambda_min + lambda_max) + I = np.eye(n) + 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): + 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, info_richardson = richardson_solver.solve() + + solution_cg, info = cg(A, b) + + if info == 0: # SciPy CG converged + 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, info_richardson, A, b) + +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") + 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) + print(f"Difference between Richardson and CG solutions: {difference:.8f}") + 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}" + +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) + 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