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)