diff --git a/README.md b/README.md index c80a8e4..f158da4 100644 --- a/README.md +++ b/README.md @@ -43,6 +43,7 @@ This package comes with a set of Jupyter notebooks designed as a structured tuto - Examples with `numethods.orthogonal`. 4. [Tutorial 4: Root-Finding Methods](tutorials/tutorial4_root_finding.ipynb) + - Bisection method. - Fixed-point iteration. - Newton’s method. @@ -52,10 +53,10 @@ This package comes with a set of Jupyter notebooks designed as a structured tuto - Examples with `numethods.roots`. - [Polynomial Regression Demo](tutorials/polynomial_regression.ipynb) - - Step-by-step example of polynomial regression. - - Shows how to fit polynomials of different degrees to data. - - Visualizes fitted curves against the original data. - - Uses `PolyFit` from `numethods.fitting`. + +- Step-by-step example of polynomial regression. +- Shows how to fit polynomials of different degrees to data. +- Visualizes fitted curves against the original data. --- @@ -98,6 +99,13 @@ This package comes with a set of Jupyter notebooks designed as a structured tuto - **Second derivative**: `SecondDerivative` - **Richardson extrapolation**: `RichardsonExtrap` +### Eigenvalue methods + +- **Power Iteration** (dominant eigenvalue/vector): `PowerIteration` +- **Inverse Power Iteration** (optionally shifted): `InversePowerIteration` +- **Rayleigh Quotient Iteration**: `RayleighQuotientIteration` +- **QR eigenvalue iteration** (unshifted, educational): `QREigenvalues` + ### Matrix & Vector utilities - Minimal `Matrix` / `Vector` classes diff --git a/examples/demo.py b/examples/demo.py index c3596aa..a73b99b 100644 --- a/examples/demo.py +++ b/examples/demo.py @@ -1,6 +1,6 @@ import sys -sys.path.append("../") +# sys.path.append("../") from numethods import * @@ -41,6 +41,9 @@ def demo_qr(): def demo_roots(): + f = lambda x: x**2 - 2 + df = lambda x: 2 * x + # Newton steps = NewtonRoot(f, df, x0=1.0).trace() print("Newton Method Trace (x^2 - 2):") @@ -85,9 +88,37 @@ def demo_differentiation(): print("Second derivative:", SecondDerivative(f, x0)) +def demo_eigen(): + # Not sure about eigenvalue methods but let's include some demos + + A = Matrix([[4, 1, 1], [1, 3, 0], [1, 0, 2]]) + print("\n=== Power Iteration ===") + solver_pi = PowerIteration(A, tol=1e-12, max_iter=100) + lam, x = solver_pi.solve() + solver_pi.trace() + print(f"Dominant eigenvalue ≈ {lam:.6f}, eigenvector ≈ {x}\n") + + print("\n=== Inverse Power Iteration (shift=0) ===") + solver_ip = InversePowerIteration(A, shift=0.0, tol=1e-12, max_iter=100) + mu, x = solver_ip.solve() + solver_ip.trace() + print(f"Smallest eigenvalue ≈ {mu:.6f}, eigenvector ≈ {x}\n") + + print("\n=== Rayleigh Quotient Iteration ===") + solver_rqi = RayleighQuotientIteration(A, tol=1e-12, max_iter=20) + mu, x = solver_rqi.solve() + solver_rqi.trace() + print(f"Eigenvalue ≈ {mu:.6f}, eigenvector ≈ {x}\n") + + M = Matrix([[3, 1, 1], [-1, 3, 1], [1, 1, 3], [0, 2, 1]]) + U, S, V = SVD(M).solve() + print("Singular values:", S) + + if __name__ == "__main__": demo_linear_solvers() demo_roots() demo_interpolation() demo_qr() demo_differentiation() + demo_eigen() diff --git a/numethods/__init__.py b/numethods/__init__.py index 40799ca..c5758e4 100644 --- a/numethods/__init__.py +++ b/numethods/__init__.py @@ -14,6 +14,13 @@ from .differentiation import ( SecondDerivative, RichardsonExtrap, ) +from .eigen import ( + PowerIteration, + InversePowerIteration, + RayleighQuotientIteration, + QREigenvalues, + SVD, +) from .solvers import LUDecomposition, GaussJordan, Jacobi, GaussSeidel, Cholesky from .roots import Bisection, FixedPoint, Secant, NewtonRoot, print_trace from .interpolation import NewtonInterpolation, LagrangeInterpolation diff --git a/numethods/eigen.py b/numethods/eigen.py new file mode 100644 index 0000000..16da6e5 --- /dev/null +++ b/numethods/eigen.py @@ -0,0 +1,219 @@ +from __future__ import annotations +from .linalg import Matrix, Vector +from .orthogonal import QRHouseholder +from .solvers import LUDecomposition +from .exceptions import NonSquareMatrixError, ConvergenceError +from .linalg import Matrix, Vector +import math + + +def solve_linear(M: Matrix, b: Vector) -> Vector: + """Solve Mx = b using LU decomposition.""" + solver = LUDecomposition(M) + return solver.solve(b) + + +class PowerIteration: + def __init__(self, A: Matrix, tol: float = 1e-10, max_iter: int = 5000): + if not A.is_square(): + raise NonSquareMatrixError("A must be square") + self.A, self.tol, self.max_iter = A, tol, max_iter + self.history = [] + + def solve(self, x0: Vector | None = None) -> tuple[float, Vector]: + n = self.A.n + x = Vector([1.0] * n) if x0 is None else x0.copy() + lam_old = 0.0 + self.history.clear() + + for k in range(self.max_iter): + y = self.A @ x + nrm = y.norm2() + if nrm == 0.0: + raise ConvergenceError("Zero vector encountered") + x = (1.0 / nrm) * y + lam = (x.dot(self.A @ x)) / (x.dot(x)) + err = abs(lam - lam_old) + + self.history.append({"iter": k, "lambda": lam, "error": err}) + + if err <= self.tol * (1.0 + abs(lam)): + return lam, x + lam_old = lam + + raise ConvergenceError("Power iteration did not converge") + + def trace(self): + if not self.history: + print("No iterations stored. Run .solve() first.") + return + print("Power Iteration Trace") + print(f"{'iter':>6} | {'lambda':>12} | {'error':>12}") + print("-" * 40) + for row in self.history: + print(f"{row['iter']:6d} | {row['lambda']:12.6e} | {row['error']:12.6e}") + + +class InversePowerIteration: + def __init__( + self, A: Matrix, shift: float = 0.0, tol: float = 1e-10, max_iter: int = 5000 + ): + if not A.is_square(): + raise NonSquareMatrixError("A must be square") + self.A, self.shift, self.tol, self.max_iter = A, shift, tol, max_iter + self.history = [] + + def solve(self, x0: Vector | None = None) -> tuple[float, Vector]: + n = self.A.n + x = Vector([1.0] * n) if x0 is None else x0.copy() + mu_old = None + self.history.clear() + + for k in range(self.max_iter): + M = Matrix( + [ + [ + self.A.data[i][j] - (self.shift if i == j else 0.0) + for j in range(n) + ] + for i in range(n) + ] + ) + y = solve_linear(M, x) + nrm = y.norm2() + if nrm == 0.0: + raise ConvergenceError("Zero vector") + x = (1.0 / nrm) * y + mu = (x.dot(self.A @ x)) / (x.dot(x)) + err = abs(mu - mu_old) if mu_old is not None else float("inf") + + self.history.append({"iter": k, "mu": mu, "error": err}) + + if (mu_old is not None) and err <= self.tol * (1.0 + abs(mu)): + return mu, x + mu_old = mu + + raise ConvergenceError("Inverse/shifted power iteration did not converge") + + def trace(self): + if not self.history: + print("No iterations stored. Run .solve() first.") + return + print("Inverse/Shifted Power Iteration Trace") + print(f"{'iter':>6} | {'mu':>12} | {'error':>12}") + print("-" * 40) + for row in self.history: + print(f"{row['iter']:6d} | {row['mu']:12.6e} | {row['error']:12.6e}") + + +class RayleighQuotientIteration: + def __init__(self, A: Matrix, tol: float = 1e-12, max_iter: int = 1000): + if not A.is_square(): + raise NonSquareMatrixError("A must be square") + self.A, self.tol, self.max_iter = A, tol, max_iter + self.history = [] + + def solve(self, x0: Vector | None = None) -> tuple[float, Vector]: + n = self.A.n + x = Vector([1.0] * n) if x0 is None else x0.copy() + x = (1.0 / x.norm2()) * x + mu = (x.dot(self.A @ x)) / (x.dot(x)) + self.history.clear() + + for k in range(self.max_iter): + M = Matrix( + [ + [self.A.data[i][j] - (mu if i == j else 0.0) for j in range(n)] + for i in range(n) + ] + ) + y = solve_linear(M, x) + x = (1.0 / y.norm2()) * y + mu_new = (x.dot(self.A @ x)) / (x.dot(x)) + err = abs(mu_new - mu) + + self.history.append({"iter": k, "mu": mu_new, "error": err}) + + if err <= self.tol * (1.0 + abs(mu_new)): + return mu_new, x + mu = mu_new + + raise ConvergenceError("Rayleigh quotient iteration did not converge") + + def trace(self): + if not self.history: + print("No iterations stored. Run .solve() first.") + return + print("Rayleigh Quotient Iteration Trace") + print(f"{'iter':>6} | {'mu':>12} | {'error':>12}") + print("-" * 40) + for row in self.history: + print(f"{row['iter']:6d} | {row['mu']:12.6e} | {row['error']:12.6e}") + + +class QREigenvalues: + def __init__(self, A: Matrix, tol: float = 1e-10, max_iter: int = 10000): + if not A.is_square(): + raise NonSquareMatrixError("A must be square") + self.A0, self.tol, self.max_iter = A.copy(), tol, max_iter + + def solve(self) -> Matrix: + A = self.A0.copy() + n = A.n + for _ in range(self.max_iter): + qr = QRHouseholder(A) + Q, R = qr.Q, qr.R + A = R @ Q + off = 0.0 + for i in range(1, n): + off += sum(abs(A.data[i][j]) for j in range(0, i)) + if off <= self.tol: + return A + raise ConvergenceError("QR did not converge") + + +class SVD: + def __init__(self, A: Matrix, tol: float = 1e-10, max_iter: int = 10000): + self.A, A = A, A + self.tol, self.max_iter = tol, max_iter + + def _eig_sym(self, S: Matrix): + n = S.n + V = Matrix.identity(n) + A = S.copy() + for _ in range(self.max_iter): + qr = QRHouseholder(A) + Q, R = qr.Q, qr.R + A = R @ Q + V = V @ Q + off = 0.0 + for i in range(1, n): + off += sum(abs(A.data[i][j]) for j in range(0, i)) + if off <= self.tol: + break + return [A.data[i][i] for i in range(n)], V + + def solve(self) -> tuple[Matrix, Vector, Matrix]: + At = self.A.transpose() + S = At @ self.A + eigvals, V = self._eig_sym(S) + idx = sorted(range(len(eigvals)), key=lambda i: eigvals[i], reverse=True) + eigvals = [eigvals[i] for i in idx] + V = Matrix([[V.data[r][i] for i in idx] for r in range(V.m)]) + sing = [math.sqrt(ev) if ev > 0 else 0.0 for ev in eigvals] + Ucols = [] + for j, sv in enumerate(sing): + vj = V.col(j) + Av = self.A @ vj + if sv > 1e-14: + uj = (1.0 / sv) * Av + else: + nrm = Av.norm2() + uj = (1.0 / nrm) * Av if nrm > 0 else Vector([0.0] * self.A.m) + nrm = uj.norm2() + uj = (1.0 / nrm) * uj if nrm > 0 else uj + Ucols.append(uj.data) + U = Matrix([[Ucols[j][i] for j in range(len(Ucols))] for i in range(self.A.m)]) + + Sigma = Vector(sing) + return U, Sigma, V diff --git a/numethods/linalg.py b/numethods/linalg.py index 3dac7fc..3ba5634 100644 --- a/numethods/linalg.py +++ b/numethods/linalg.py @@ -1,6 +1,7 @@ from __future__ import annotations from typing import Iterable, Tuple, List, Union from .exceptions import NonSquareMatrixError, SingularMatrixError +import math Number = float # We'll use float throughout @@ -97,19 +98,26 @@ class Matrix: def col(self, j: int) -> Vector: return Vector(self.data[i][j] for i in range(self.m)) - def norm(self) -> Number: - return ( - max(sum(abs(self.data[i][j]) for i in range(self.m)) for j in range(self.n)) - if self.n > 0 - else 0.0 + def norm(self) -> float: + """Matrix 1-norm: max column sum.""" + return max( + sum(abs(self.data[i][j]) for i in range(self.m)) for j in range(self.n) ) - def norm_inf(self) -> Number: - return ( - max(sum(abs(self.data[i][j]) for j in range(self.n)) for i in range(self.m)) - if self.m > 0 - else 0.0 - ) + def norm_inf(self) -> float: + """Matrix infinity norm: max row sum.""" + return max(sum(abs(v) for v in row) for row in self.data) + + def norm2(self, tol: float = 1e-10, max_iter: int = 5000) -> float: + """Spectral norm via power iteration on A^T A.""" + # lazy import, avoids circular import + from .eigen import PowerIteration + + if not self.is_square(): + raise NonSquareMatrixError("Spectral norm requires square matrix") + AtA = self.T @ self + lam, _ = PowerIteration(AtA, tol=tol, max_iter=max_iter).solve() + return math.sqrt(lam) def norm_fro(self) -> Number: return ( @@ -117,6 +125,44 @@ class Matrix: ** 0.5 ) + def inverse(self) -> "Matrix": + # lazy import, avoids circular import + from .solvers import LUDecomposition + + """Compute A^{-1} using LU decomposition.""" + if not self.is_square(): + raise NonSquareMatrixError("Inverse requires square matrix") + n = self.n + solver = LUDecomposition(self) + cols = [] + for j in range(n): + e = Vector([1.0 if i == j else 0.0 for i in range(n)]) + x = solver.solve(e) + cols.append([x[i] for i in range(n)]) + return Matrix([[cols[j][i] for j in range(n)] for i in range(n)]) + + def condition_number(self, norm: str = "2") -> float: + """ + Compute condition number κ(A) = ||A|| * ||A^{-1}||. + norm: "1", "inf", or "2" + """ + if not self.is_square(): + raise NonSquareMatrixError("Condition number requires square matrix") + + if norm == "1": + normA = self.norm() + normAinv = self.inverse().norm() + elif norm == "inf": + normA = self.norm_inf() + normAinv = self.inverse().norm_inf() + elif norm == "2": + normA = self.norm2() + normAinv = self.inverse().norm2() + else: + raise ValueError("norm must be '1', 'inf', or '2'") + + return normA * normAinv + def __add__(self, other: "Matrix") -> "Matrix": if not isinstance(other, Matrix): raise TypeError("Can only add Matrix with Matrix")