update eigen.py

fixed -> matrix and vector imported twice from linalg

.gitignore

@@ -1,38 +1,3 @@
-Logo
-Issues
-Pull Requests
-Milestones
-Explore
-PhiloMath
-/
-data-fitting-models
-Private
-generated from PhiloMath/default-project-template
-Code
-Issues
-Pull Requests
-Actions
-Packages
-Projects
-Releases
-Wiki
-Activity
-Settings
-Files
-.gitignore
-LICENSE
-README.md
-data-fitting-models
-/
-.gitignore
-
-Ömür Uğur
-f9860d7f48
-Initial commit
-2 days ago
-357 lines
-6.4 KiB
-Plaintext
 # ---> JupyterNotebooks
 # gitignore template for Jupyter Notebooks
 # website: http://jupyter.org/

README.md

@@ -10,53 +10,9 @@ A lightweight, from-scratch, object-oriented Python package implementing classic
 - Lightweight, no dependencies.
 - Consistent object-oriented API (.solve() etc).
 
----
-
 ## Tutorial Series
 
-This package comes with a set of Jupyter notebooks designed as a structured tutorial series in **numerical methods**, both mathematically rigorous and hands-on with code.
+This package comes with a set of Jupyter notebooks designed as a structured tutorial series in **numerical methods**, both mathematically rigorous and hands-on with code. See [Tutorials](./tutorials/README.md).
-
-### Core Tutorials
-
-1. [Tutorial 1: Vectors and Matrices](tutorials/tutorial1_vectors.ipynb)
-
-   - Definitions of vectors and matrices.
-   - Vector operations: addition, scalar multiplication, dot product, norms.
-   - Matrix operations: addition, multiplication, transpose, inverse.
-   - Matrix and vector norms.
-   - Examples with `numethods.linalg`.
-
-2. [Tutorial 2: Linear Systems of Equations](tutorials/tutorial2_linear_systems.ipynb)
-
-   - Gaussian elimination and Gauss–Jordan.
-   - LU decomposition.
-   - Cholesky decomposition.
-   - Iterative methods: Jacobi and Gauss-Seidel.
-   - Examples with `numethods.solvers`.
-
-3. [Tutorial 3: Orthogonalization and QR Factorization](tutorials/tutorial3_orthogonalization.ipynb)
-
-   - Inner products and orthogonality.
-   - Gram–Schmidt process (classical and modified).
-   - Householder reflections.
-   - QR decomposition and applications.
-   - Examples with `numethods.orthogonal`.
-
-4. [Tutorial 4: Root-Finding Methods](tutorials/tutorial4_root_finding.ipynb)
-   - Bisection method.
-   - Fixed-point iteration.
-   - Newton’s method.
-   - Secant method.
-   - Convergence analysis and error behavior.
-   - Trace outputs for iteration history.
-   - 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.
-
----
 
 ## Features
 
@@ -97,6 +53,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

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()

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

numethods/eigen.py

@@ -0,0 +1,218 @@
+from __future__ import annotations
+from .linalg import Matrix, Vector
+from .orthogonal import QRHouseholder
+from .solvers import LUDecomposition
+from .exceptions import NonSquareMatrixError, ConvergenceError
+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

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:
+    def norm(self) -> float:
-        return (
+        """Matrix 1-norm: max column sum."""
-            max(sum(abs(self.data[i][j]) for i in range(self.m)) for j in range(self.n))
+        return max(
-            if self.n > 0
+            sum(abs(self.data[i][j]) for i in range(self.m)) for j in range(self.n)
-            else 0.0
         )
 
-    def norm_inf(self) -> Number:
+    def norm_inf(self) -> float:
-        return (
+        """Matrix infinity norm: max row sum."""
-            max(sum(abs(self.data[i][j]) for j in range(self.n)) for i in range(self.m))
+        return max(sum(abs(v) for v in row) for row in self.data)
-            if self.m > 0
+
-            else 0.0
+    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")

tutorials/README.md

@@ -0,0 +1,47 @@
+# Tutorial Series
+
+This package comes with a set of Jupyter notebooks designed as a structured tutorial series in **numerical methods**, both mathematically rigorous and hands-on with code.
+
+## Core Tutorials
+
+1. [Tutorial 1: Vectors and Matrices](./tutorial1_vectors_matrices.ipynb)
+
+   - Definitions of vectors and matrices.
+   - Vector operations: addition, scalar multiplication, dot product, norms.
+   - Matrix operations: addition, multiplication, transpose, inverse.
+   - Matrix and vector norms.
+   - Examples with `numethods.linalg`.
+
+2. [Tutorial 2: Linear Systems of Equations](./tutorial2_linear_systems.ipynb)
+
+   - Gaussian elimination and Gauss–Jordan.
+   - LU decomposition.
+   - Cholesky decomposition.
+   - Iterative methods: Jacobi and Gauss-Seidel.
+   - Examples with `numethods.solvers`.
+
+3. [Tutorial 3: Orthogonalization and QR Factorization](./tutorial3_orthogonalization.ipynb)
+
+   - Inner products and orthogonality.
+   - Gram–Schmidt process (classical and modified).
+   - Householder reflections.
+   - QR decomposition and applications.
+   - Examples with `numethods.orthogonal`.
+
+4. [Tutorial 4: Root-Finding Methods](./tutorial4_root_finding.ipynb)
+
+   - Bisection method.
+   - Fixed-point iteration.
+   - Newton’s method.
+   - Secant method.
+   - Convergence analysis and error behavior.
+   - Trace outputs for iteration history.
+   - Examples with `numethods.roots`.
+
+- [Polynomial Regression Demo](./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.
+
+---

tutorials/tutorial3_orthogonalization.ipynb

@@ -94,7 +94,7 @@
     "q_1 = \\frac{a_1}{\\|a_1\\|}\n",
     "$$\n",
     "$$\n",
-    "q_k = \\frac{a_k - \\sum_{j=1}^{k-1} (q_j \\cdot a_k) q_j}{\\left\\|a_k - \\sum_{j=1}^{k-1} (q_j \\cdot a_k) q_j\\right\\|}\n",
+    "q_k = \\frac{a_k - \\sum_{j=1}^{k-1} (q_j \\cdot a_k) q_j}{\\left\\|a_k - \\sum_{j=1}^{k-1} (q_j \\cdot a_k) q_j\\right\\|}, \\qquad k = 2, \\ldots, n\n",
     "$$\n",
     "\n",
     "Matrix form:\n",
@@ -236,17 +236,17 @@
     "\n",
     "We want to solve\n",
     "\n",
-    "$$ \\min_x \\|Ax - b\\|_2. $$\n",
+    "$$ \\min_x \\Vert Ax - b \\Vert_2^2. $$\n",
     "\n",
     "If $A = QR$, then\n",
     "\n",
-    "$$ \\min_x \\|Ax - b\\|_2 = \\min_x \\|QRx - b\\|_2. $$\n",
+    "$$ \\min_x \\Vert Ax - b \\Vert_2^2 = \\min_x \\Vert QRx - b \\Vert_2^2. $$\n",
     "\n",
-    "Since $Q$ has orthonormal columns:\n",
+    "Since $Q$ has orthonormal columns, and the normal equations boils down to\n",
     "\n",
-    "$$ R x = Q^T b. $$\n",
+    "$$ R x = Q^T b, $$\n",
     "\n",
-    "So we can solve using back-substitution.\n"
+    "we can therefore solve for $x$ by using back-substitution.\n"
    ]
   },
   {

tutorials/tutorial4_root_finding.ipynb

@@ -55,7 +55,7 @@
    "source": [
     "## 2. Bisection Method\n",
     "\n",
-    "**Assumption (Intermediate Value Theorem):** If f is continuous on ([a,b]) and (f(a),f(b) < 0),\n",
+    "**Assumption (Intermediate Value Theorem):** If f is continuous on $[a,b]$ and $f(a)f(b) < 0$,\n",
     "then there exists $x^\\star$ in (a,b) with $f(x^\\star)=0$.\n",
     "\n",
     "- Assumes $f$ is continuous on $[a,b]$ with $f(a)f(b)<0$.\n",