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/

DISCLAIMER.md

@@ -0,0 +1,16 @@
+Disclaimer
+
+This package is developed for learning and research purposes only.
+
+Most of the methods here are inspired by or extended from homework assignments I worked on during my PhD coursework over the past two years. What you see are more object-oriented and "fancier" versions of those exercises.
+
+The algorithms are based primarily on course notes and Scientific Computing: An Introductory Survey by Michael T. Heath. Since I will be taking my PhD qualifying exam in a couple of months (wish me luck!), I've been preparing these codes and notes as part of my study, and I wanted to share them.
+
+I'm sure many of you reading this already have a deeper understanding of these topics than I do. If you find some implementations or explanations overly simple, that's intentional: my goal here is to understand and learn, step by step.
+
+There are mistakes, and there will be mistakes. As our dear Ömür Hoca always says:
+
+“Do not trust me. Do not trust someone else. Only trust yourself.”
+
+Please use these materials as a study companion, not as a definitive reference.
+Good luck, and have fun exploring numerical methods!

README.md

@@ -1,8 +1,19 @@
 # numethods
 
-A small, from-scratch, object-oriented Python package implementing classic numerical methods.  
+A lightweight, from-scratch, object-oriented Python package implementing classic numerical methods.  
 **No NumPy / SciPy solvers used**, algorithms are implemented transparently for learning and research.
 
+## Why this might be useful
+
+- Great for teaching/learning numerical methods step by step.
+- Good reference for people writing their own solvers in C/Fortran/Julia.
+- 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. See [Tutorials](./tutorials/README.md).
+
 ## Features
 
 ### Linear system solvers
@@ -25,7 +36,7 @@ A small, from-scratch, object-oriented Python package implementing classic numer
 - **Newton** (divided differences): `NewtonInterpolation`
 - **Lagrange** polynomials: `LagrangeInterpolation`
 
-### Orthogonalization, QR, and Least Squares (NEW)
+### Orthogonalization, QR, and Least Squares
 
 - **Classical Gram-Schmidt**: `QRGramSchmidt`
 - **Modified Gram-Schmidt**: `QRModifiedGramSchmidt`
@@ -33,13 +44,30 @@ A small, from-scratch, object-oriented Python package implementing classic numer
 - **QR-based linear solver** (square systems): `QRSolver`
 - **Least Squares** for overdetermined systems (via QR): `LeastSquaresSolver`
 
+### Numerical Differentiation
+
+- **Forward difference**: `ForwardDiff`
+- **Backward difference**: `BackwardDiff`
+- **Central difference (2nd order)**: `CentralDiff`
+- **Central difference (4th order)**: `CentralDiff4th`
+- **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
-- `@` operator for **matrix multiplication** (NEW)
+- `@` operator for **matrix multiplication**
 - `*` for **scalar**-matrix multiplication
 - `.T` for transpose
 - Forward / backward substitution helpers
+- Norms, dot products, row/column access
 
 ---
 

examples/demo.py

@@ -1,28 +1,8 @@
 import sys
-sys.path.append('../')
+
+# sys.path.append("../")
 from numethods import *
 
-A = Matrix([[2, -1], [1, 2], [1, 1]])
-b = Vector([1, 2, 3])
-
-# Factorization
-qr = QRHouseholder(A)
-Q, R = qr.Q, qr.R
-print("Q =", Q)
-print("R =", R)
-
-qrm = QRModifiedGramSchmidt(A)
-Qm, Rm = qrm.Q, qrm.R
-print("Qm =", Qm)
-print("Rm =", Rm)
-print("Q^T Q =", Q.T @ Q)
-print("Qm^T Qm =", Qm.T @ Qm)
-print("A=Qm Rm =", Qm @ Rm)
-
-# Solve Ax = b (least squares, since A is tall)
-x_ls = LeastSquaresSolver(A, b).solve()
-print("Least squares solution:", x_ls)
-
 
 def demo_linear_solvers():
     A = Matrix([[4, -1, 0], [-1, 4, -1], [0, -1, 3]])
@@ -35,15 +15,55 @@ def demo_linear_solvers():
     print("Gauss-Seidel:", GaussSeidel(A, b, tol=1e-12).solve())
 
 
+def demo_qr():
+    A = Matrix([[2, -1], [1, 2], [1, 1]])
+
+    b = Vector([1, 2, 3])
+
+    # Factorization
+    qr = QRHouseholder(A)
+    Q, R = qr.Q, qr.R
+    print("Q =", Q)
+    print("R =", R)
+
+    qrm = QRModifiedGramSchmidt(A)
+    Qm, Rm = qrm.Q, qrm.R
+    print("Qm =", Qm)
+    print("Rm =", Rm)
+    print("Q^T Q =", Q.T @ Q)
+    print("Qm^T Qm =", Qm.T @ Qm)
+    print("A=Qm Rm =", Qm @ Rm)
+    print("A=Q R =", Q @ R)
+
+    # Solve Ax = b (least squares, since A is tall)
+    x_ls = LeastSquaresSolver(A, b).solve()
+    print("Least squares solution:", x_ls)
+
+
 def demo_roots():
-    f = lambda x: x**3 - x - 2
-    df = lambda x: 3 * x**2 - 1
-    print("Bisection:", Bisection(f, 1, 2).solve())
-    print("Secant:", Secant(f, 1.0, 2.0).solve())
-    print("Newton root:", NewtonRoot(f, df, 1.5).solve())
-    # A simple contraction for demonstration; not general-purpose
-    g = lambda x: (x + 2 / x**2) / 2
-    print("Fixed point (demo):", FixedPoint(g, 1.5).solve())
+    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):")
+    print_trace(steps)
+
+    # Secant
+    steps = Secant(f, 0, 2).trace()
+    print("\nSecant Method Trace (x^2 - 2):")
+    print_trace(steps)
+
+    # Bisection
+    steps = Bisection(f, 0, 2).trace()
+    print("\nBisection Method Trace (x^2 - 2):")
+    print_trace(steps)
+
+    # Fixed-point: solve
+    g = lambda x: 0.5 * (x + 2 / x)
+    steps = FixedPoint(g, 1.0).trace()
+    print("\nFixed-Point Iteration Trace (x^2 - 2):")
+    print_trace(steps)
 
 
 def demo_interpolation():
@@ -56,7 +76,49 @@ def demo_interpolation():
     print("Lagrange interpolation at", t, "=", lagr.evaluate(t))
 
 
+def demo_differentiation():
+    f = lambda x: x**3  # f'(x) = 3x^2, f''(x) = 6x
+    x0 = 2.0
+
+    print("Forward  :", ForwardDiff(f, x0))
+    print("Backward :", BackwardDiff(f, x0))
+    print("Central  :", CentralDiff(f, x0))
+    print("4th order:", CentralDiff4th(f, x0))
+    print("Richardson:", RichardsonExtrap(f, x0))
+    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

@@ -6,7 +6,22 @@ from .orthogonal import (
     QRSolver,
     LeastSquaresSolver,
 )
+from .differentiation import (
+    ForwardDiff,
+    BackwardDiff,
+    CentralDiff,
+    CentralDiff4th,
+    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
+from .roots import Bisection, FixedPoint, Secant, NewtonRoot, print_trace
 from .interpolation import NewtonInterpolation, LagrangeInterpolation
 from .exceptions import *

numethods/differentiation.py

@@ -0,0 +1,51 @@
+from __future__ import annotations
+from typing import Callable
+
+
+# ----------------------------
+# First derivative approximations
+# ----------------------------
+
+
+def ForwardDiff(f: Callable[[float], float], x: float, h: float = 1e-5) -> float:
+    """Forward finite difference approximation of f'(x)."""
+    return (f(x + h) - f(x)) / h
+
+
+def BackwardDiff(f: Callable[[float], float], x: float, h: float = 1e-5) -> float:
+    """Backward finite difference approximation of f'(x)."""
+    return (f(x) - f(x - h)) / h
+
+
+def CentralDiff(f: Callable[[float], float], x: float, h: float = 1e-5) -> float:
+    """Central finite difference approximation of f'(x) (2nd-order accurate)."""
+    return (f(x + h) - f(x - h)) / (2 * h)
+
+
+def CentralDiff4th(f: Callable[[float], float], x: float, h: float = 1e-5) -> float:
+    """Fourth-order accurate central difference approximation of f'(x)."""
+    return (-f(x + 2 * h) + 8 * f(x + h) - 8 * f(x - h) + f(x - 2 * h)) / (12 * h)
+
+
+# ----------------------------
+# Second derivative
+# ----------------------------
+
+
+def SecondDerivative(f: Callable[[float], float], x: float, h: float = 1e-5) -> float:
+    """Central difference approximation of second derivative f''(x)."""
+    return (f(x + h) - 2 * f(x) + f(x - h)) / (h**2)
+
+
+# ----------------------------
+# Richardson Extrapolation
+# ----------------------------
+
+
+def RichardsonExtrap(f: Callable[[float], float], x: float, h: float = 1e-2) -> float:
+    """Richardson extrapolation to improve derivative accuracy.
+    Combines estimates with step h and h/2.
+    """
+    D_h = CentralDiff(f, x, h)
+    D_h2 = CentralDiff(f, x, h / 2)
+    return (4 * D_h2 - D_h) / 3

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
+from typing import Iterable, Tuple, List, Union
 from .exceptions import NonSquareMatrixError, SingularMatrixError
+import math
 
 Number = float  # We'll use float throughout
 
@@ -21,16 +22,24 @@ class Vector:
     def copy(self) -> "Vector":
         return Vector(self.data[:])
 
+    def norm(self) -> Number:
+        return sum(abs(x) for x in self.data)
+
     def norm_inf(self) -> Number:
         return max(abs(x) for x in self.data) if self.data else 0.0
 
+    def norm2(self) -> Number:
+        return sum(x * x for x in self.data) ** 0.5
+
     def __add__(self, other: "Vector") -> "Vector":
-        assert len(self) == len(other)
-        return Vector(a + b for a, b in zip(self.data, other.data))
+        if len(self) != len(other):
+            raise ValueError("Vector dimensions must match for addition")
+        return Vector([a + b for a, b in zip(self.data, other.data)])
 
     def __sub__(self, other: "Vector") -> "Vector":
-        assert len(self) == len(other)
-        return Vector(a - b for a, b in zip(self.data, other.data))
+        if len(self) != len(other):
+            raise ValueError("Vector dimensions must match for subtraction")
+        return Vector([a - b for a, b in zip(self.data, other.data)])
 
     def __mul__(self, scalar: Number) -> "Vector":
         return Vector(scalar * x for x in self.data)
@@ -89,30 +98,129 @@ class Matrix:
     def col(self, j: int) -> Vector:
         return Vector(self.data[i][j] for i in range(self.m))
 
+    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) -> 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 (
+            sum(self.data[i][j] ** 2 for i in range(self.m) for j in range(self.n))
+            ** 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")
+        if self.m != other.m or self.n != other.n:
+            raise ValueError("Matrix dimensions must match for addition")
+        return Matrix(
+            [
+                [self.data[i][j] + other.data[i][j] for j in range(self.n)]
+                for i in range(self.m)
+            ]
+        )
+
+    def __sub__(self, other: "Matrix") -> "Matrix":
+        if not isinstance(other, Matrix):
+            raise TypeError("Can only subtract Matrix with Matrix")
+        if self.m != other.m or self.n != other.n:
+            raise ValueError("Matrix dimensions must match for subtraction")
+        return Matrix(
+            [
+                [self.data[i][j] - other.data[i][j] for j in range(self.n)]
+                for i in range(self.m)
+            ]
+        )
+
     def transpose(self) -> "Matrix":
         return Matrix([[self.data[i][j] for i in range(self.m)] for j in range(self.n)])
 
     T = property(transpose)
 
-    def __matmul__(self, other: "Matrix") -> "Matrix":
-        """Matrix multiplication with @ operator."""
-        if self.n != other.m:
-            raise ValueError("Matrix dimensions do not align for multiplication")
-        return Matrix(
-            [
+    def __matmul__(self, other: Union["Matrix", "Vector"]):
+        if isinstance(other, Matrix):
+            if self.n != other.m:
+                raise ValueError("dims")
+            return Matrix(
                 [
-                    sum(self.data[i][k] * other.data[k][j] for k in range(self.n))
-                    for j in range(other.n)
+                    [
+                        sum(self.data[i][k] * other.data[k][j] for k in range(self.n))
+                        for j in range(other.n)
+                    ]
+                    for i in range(self.m)
                 ]
-                for i in range(self.m)
-            ]
-        )
+            )
+        elif isinstance(other, Vector):
+            if self.n != len(other):
+                raise ValueError("dims")
+            return Vector(
+                [
+                    sum(self.data[i][k] * other[k] for k in range(self.n))
+                    for i in range(self.m)
+                ]
+            )
+        else:
+            raise TypeError("Unsupported @")
 
-    def __mul__(self, other):
-        """Overload * for scalar multiplication."""
-        if isinstance(other, (int, float)):
-            return Matrix([[val * other for val in row] for row in self.data])
-        raise TypeError("Use @ for matrix multiplication, * only supports scalars")
+    def __mul__(self, s):
+        if isinstance(s, (int, float)):
+            return Matrix([[v * s for v in row] for row in self.data])
+        raise TypeError("Use @ for matrix multiply; * is scalar")
 
     __rmul__ = __mul__
 
@@ -142,6 +250,7 @@ class Matrix:
 
 
 def forward_substitution(L: Matrix, b: Vector) -> Vector:
+    """Solve Lx = b for x using forward substitution"""
     if not L.is_square():
         raise NonSquareMatrixError("L must be square")
     n = L.n
@@ -155,6 +264,7 @@ def forward_substitution(L: Matrix, b: Vector) -> Vector:
 
 
 def backward_substitution(U: Matrix, b: Vector) -> Vector:
+    """Solve Ux = b for x using backward substitution"""
     if not U.is_square():
         raise NonSquareMatrixError("U must be square")
     n = U.n

numethods/orthogonal.py

@@ -79,7 +79,7 @@ class QRHouseholder:
             sign = 1.0 if x[0] >= 0 else -1.0
             u1 = x[0] + sign * normx
             v = [xi / u1 if i > 0 else 1.0 for i, xi in enumerate(x)]
-            normv = sum(vi * vi for vi in v) ** 0.5
+            normv = Vector(v).norm2()
             v = [vi / normv for vi in v]
             for j in range(k, n):
                 s = sum(v[i] * self.R.data[k + i][j] for i in range(len(v)))
@@ -89,7 +89,7 @@ class QRHouseholder:
                 s = sum(v[i] * self.Q.data[j][k + i] for i in range(len(v)))
                 for i in range(len(v)):
                     self.Q.data[j][k + i] -= 2 * s * v[i]
-#        self.Q = self.Q.transpose()
+        # self.Q = self.Q.transpose()
 
 
 class QRSolver:

numethods/roots.py

@@ -1,9 +1,17 @@
 from __future__ import annotations
-from typing import Callable
+from typing import Callable, List, Dict, Any
 from .exceptions import ConvergenceError, DomainError
 
+
 class Bisection:
-    def __init__(self, f: Callable[[float], float], a: float, b: float, tol: float = 1e-10, max_iter: int = 10_000):
+    def __init__(
+        self,
+        f: Callable[[float], float],
+        a: float,
+        b: float,
+        tol: float = 1e-10,
+        max_iter: int = 10_000,
+    ):
         if a >= b:
             raise ValueError("Require a < b")
         fa, fb = f(a), f(b)
@@ -16,11 +24,38 @@ class Bisection:
         a, b, f = self.a, self.b, self.f
         fa, fb = f(a), f(b)
         for _ in range(self.max_iter):
-            c = 0.5*(a+b)
+            c = 0.5 * (a + b)
             fc = f(c)
-            if abs(fc) <= self.tol or 0.5*(b-a) <= self.tol:
+            if abs(fc) <= self.tol or 0.5 * (b - a) <= self.tol:
                 return c
-            if fa*fc < 0:
+            if fa * fc < 0:
+                b, fb = c, fc
+            else:
+                a, fa = c, fc
+        raise ConvergenceError("Bisection did not converge")
+
+    def trace(self) -> List[Dict[str, Any]]:
+        steps = []
+        a, b, f = self.a, self.b, self.f
+        fa, fb = f(a), f(b)
+        for k in range(self.max_iter):
+            c = 0.5 * (a + b)
+            fc = f(c)
+            steps.append(
+                {
+                    "iter": k,
+                    "a": a,
+                    "b": b,
+                    "c": c,
+                    "f(a)": fa,
+                    "f(b)": fb,
+                    "f(c)": fc,
+                    "interval": b - a,
+                }
+            )
+            if abs(fc) <= self.tol or 0.5 * (b - a) <= self.tol:
+                return steps
+            if fa * fc < 0:
                 b, fb = c, fc
             else:
                 a, fa = c, fc
@@ -28,7 +63,13 @@ class Bisection:
 
 
 class FixedPoint:
-    def __init__(self, g: Callable[[float], float], x0: float, tol: float = 1e-10, max_iter: int = 10_000):
+    def __init__(
+        self,
+        g: Callable[[float], float],
+        x0: float,
+        tol: float = 1e-10,
+        max_iter: int = 10_000,
+    ):
         self.g, self.x0, self.tol, self.max_iter = g, x0, tol, max_iter
 
     def solve(self) -> float:
@@ -40,28 +81,79 @@ class FixedPoint:
             x = x_new
         raise ConvergenceError("Fixed-point iteration did not converge")
 
+    def trace(self) -> List[Dict[str, Any]]:
+        steps = []
+        x = self.x0
+        for k in range(self.max_iter):
+            x_new = self.g(x)
+            steps.append({"iter": k, "x": x, "x_new": x_new, "error": abs(x_new - x)})
+            if abs(x_new - x) <= self.tol * (1.0 + abs(x_new)):
+                return steps
+            x = x_new
+        raise ConvergenceError("Fixed-point iteration did not converge")
+
 
 class Secant:
-    def __init__(self, f: Callable[[float], float], x0: float, x1: float, tol: float = 1e-10, max_iter: int = 10_000):
+    def __init__(
+        self,
+        f: Callable[[float], float],
+        x0: float,
+        x1: float,
+        tol: float = 1e-10,
+        max_iter: int = 10_000,
+    ):
         self.f, self.x0, self.x1, self.tol, self.max_iter = f, x0, x1, tol, max_iter
 
     def solve(self) -> float:
         x0, x1, f = self.x0, self.x1, self.f
         f0, f1 = f(x0), f(x1)
         for _ in range(self.max_iter):
-            denom = (f1 - f0)
+            denom = f1 - f0
             if abs(denom) < 1e-20:
                 raise ConvergenceError("Secant encountered nearly zero denominator")
-            x2 = x1 - f1*(x1 - x0)/denom
+            x2 = x1 - f1 * (x1 - x0) / denom
             if abs(x2 - x1) <= self.tol * (1.0 + abs(x2)):
                 return x2
             x0, x1 = x1, x2
             f0, f1 = f1, f(x1)
         raise ConvergenceError("Secant did not converge")
 
+    def trace(self) -> List[Dict[str, Any]]:
+        steps = []
+        x0, x1, f = self.x0, self.x1, self.f
+        f0, f1 = f(x0), f(x1)
+        for k in range(self.max_iter):
+            denom = f1 - f0
+            if abs(denom) < 1e-20:
+                raise ConvergenceError("Secant encountered nearly zero denominator")
+            x2 = x1 - f1 * (x1 - x0) / denom
+            steps.append(
+                {
+                    "iter": k,
+                    "x0": x0,
+                    "x1": x1,
+                    "x2": x2,
+                    "f(x0)": f0,
+                    "f(x1)": f1,
+                    "error": abs(x2 - x1),
+                }
+            )
+            if abs(x2 - x1) <= self.tol * (1.0 + abs(x2)):
+                return steps
+            x0, x1 = x1, x2
+            f0, f1 = f1, f(x1)
+        raise ConvergenceError("Secant did not converge")
+
 
 class NewtonRoot:
-    def __init__(self, f: Callable[[float], float], df: Callable[[float], float], x0: float, tol: float = 1e-10, max_iter: int = 10_000):
+    def __init__(
+        self,
+        f: Callable[[float], float],
+        df: Callable[[float], float],
+        x0: float,
+        tol: float = 1e-10,
+        max_iter: int = 10_000,
+    ):
         self.f, self.df, self.x0, self.tol, self.max_iter = f, df, x0, tol, max_iter
 
     def solve(self) -> float:
@@ -70,8 +162,50 @@ class NewtonRoot:
             dfx = self.df(x)
             if abs(dfx) < 1e-20:
                 raise ConvergenceError("Derivative near zero in Newton method")
-            x_new = x - self.f(x)/dfx
+            x_new = x - self.f(x) / dfx
             if abs(x_new - x) <= self.tol * (1.0 + abs(x_new)):
                 return x_new
             x = x_new
         raise ConvergenceError("Newton method did not converge")
+
+    def trace(self) -> List[Dict[str, Any]]:
+        steps = []
+        x = self.x0
+        for k in range(self.max_iter):
+            dfx = self.df(x)
+            if abs(dfx) < 1e-20:
+                raise ConvergenceError("Derivative near zero in Newton method")
+            x_new = x - self.f(x) / dfx
+            steps.append(
+                {
+                    "iter": k,
+                    "x": x,
+                    "f(x)": self.f(x),
+                    "df(x)": dfx,
+                    "x_new": x_new,
+                    "error": abs(x_new - x),
+                }
+            )
+            if abs(x_new - x) <= self.tol * (1.0 + abs(x_new)):
+                return steps
+            x = x_new
+        raise ConvergenceError("Newton method did not converge")
+
+
+def print_trace(steps: List[Dict[str, Any]]):
+    if not steps:
+        print("No steps recorded.")
+        return
+    # Get headers from dict keys
+    headers = list(steps[0].keys())
+    # Print header
+    print(" | ".join(f"{h:>10}" for h in headers))
+    print("-" * (13 * len(headers)))
+    # Print rows
+    for row in steps:
+        print(
+            " | ".join(
+                f"{row[h]:>10.6g}" if isinstance(row[h], (int, float)) else str(row[h])
+                for h in headers
+            )
+        )

numethods/solvers.py

@@ -1,9 +1,17 @@
 from __future__ import annotations
 from .linalg import Matrix, Vector, forward_substitution, backward_substitution
-from .exceptions import NonSquareMatrixError, SingularMatrixError, NotSymmetricError, NotPositiveDefiniteError, ConvergenceError
+from .exceptions import (
+    NonSquareMatrixError,
+    SingularMatrixError,
+    NotSymmetricError,
+    NotPositiveDefiniteError,
+    ConvergenceError,
+)
+
 
 class LUDecomposition:
-    """LU with partial pivoting: PA = LU"""
+    """LU decomposition with partial pivoting: PA = LU"""
+
     def __init__(self, A: Matrix):
         if not A.is_square():
             raise NonSquareMatrixError("A must be square")
@@ -11,6 +19,7 @@ class LUDecomposition:
         self.L = Matrix.identity(self.n)
         self.U = A.copy()
         self.P = Matrix.identity(self.n)
+        self.steps: list[tuple[int, Matrix, Matrix, Matrix]] = []  # store pivot steps
         self._decompose()
 
     def _decompose(self) -> None:
@@ -22,26 +31,47 @@ class LUDecomposition:
             self.U.swap_rows(k, pivot_row)
             self.P.swap_rows(k, pivot_row)
             if k > 0:
-                self.L.data[k][:k], self.L.data[pivot_row][:k] = self.L.data[pivot_row][:k], self.L.data[k][:k]
-            for i in range(k+1, n):
+                self.L.data[k][:k], self.L.data[pivot_row][:k] = (
+                    self.L.data[pivot_row][:k],
+                    self.L.data[k][:k],
+                )
+            for i in range(k + 1, n):
                 m = self.U.data[i][k] / self.U.data[k][k]
                 self.L.data[i][k] = m
                 for j in range(k, n):
                     self.U.data[i][j] -= m * self.U.data[k][j]
+            # record step
+            self.steps.append((k, self.L.copy(), self.U.copy(), self.P.copy()))
 
     def solve(self, b: Vector) -> Vector:
-        Pb = Vector([sum(self.P.data[i][j]*b[j] for j in range(self.n)) for i in range(self.n)])
+        Pb = Vector(
+            [
+                sum(self.P.data[i][j] * b[j] for j in range(self.n))
+                for i in range(self.n)
+            ]
+        )
         y = forward_substitution(self.L, Pb)
         x = backward_substitution(self.U, y)
         return x
 
+    def trace(self):
+        print("LU Decomposition Trace (steps of elimination)")
+        for k, L, U, P in self.steps:
+            print(f"\nStep {k}:")
+            print(f"L = {L}")
+            print(f"U = {U}")
+            print(f"P = {P}")
+
 
 class GaussJordan:
+    """Gauss-Jordan elimination."""
+
     def __init__(self, A: Matrix):
         if not A.is_square():
             raise NonSquareMatrixError("A must be square")
         self.n = A.n
         self.A = A.copy()
+        self.steps: list[tuple[int, Matrix]] = []
 
     def solve(self, b: Vector) -> Vector:
         n = self.n
@@ -57,12 +87,26 @@ class GaussJordan:
                 if r == col:
                     continue
                 factor = Ab.data[r][col]
-                Ab.data[r] = [rv - factor*cv for rv, cv in zip(Ab.data[r], Ab.data[col])]
+                Ab.data[r] = [
+                    rv - factor * cv for rv, cv in zip(Ab.data[r], Ab.data[col])
+                ]
+            # record step
+            self.steps.append((col, Ab.copy()))
         return Vector(row[-1] for row in Ab.data)
 
+    def trace(self):
+        print("Gauss-Jordan Trace (row reduction steps)")
+        for step, Ab in self.steps:
+            print(f"\nColumn {step}:")
+            print(f"Augmented matrix = {Ab}")
+
 
 class Jacobi:
-    def __init__(self, A: Matrix, b: Vector, tol: float = 1e-10, max_iter: int = 10_000):
+    """Jacobi iterative method for Ax = b."""
+
+    def __init__(
+        self, A: Matrix, b: Vector, tol: float = 1e-10, max_iter: int = 10_000
+    ):
         if not A.is_square():
             raise NonSquareMatrixError("A must be square")
         if A.n != len(b):
@@ -71,27 +115,42 @@ class Jacobi:
         self.b = b.copy()
         self.tol = tol
         self.max_iter = max_iter
+        self.history: list[float] = []
 
     def solve(self, x0: Vector | None = None) -> Vector:
         n = self.A.n
-        x = Vector([0.0]*n) if x0 is None else x0.copy()
+        x = Vector([0.0] * n) if x0 is None else x0.copy()
         for _ in range(self.max_iter):
-            x_new = [0.0]*n
+            x_new = [0.0] * n
             for i in range(n):
                 diag = self.A.data[i][i]
                 if abs(diag) < 1e-15:
                     raise SingularMatrixError("Zero diagonal entry in Jacobi")
-                s = sum(self.A.data[i][j]*x[j] for j in range(n) if j != i)
+                s = sum(self.A.data[i][j] * x[j] for j in range(n) if j != i)
                 x_new[i] = (self.b[i] - s) / diag
             x_new = Vector(x_new)
-            if (x_new - x).norm_inf() <= self.tol * (1.0 + x_new.norm_inf()):
+            r = (self.A @ x_new) - self.b
+            res_norm = r.norm2()
+            self.history.append(res_norm)
+            if res_norm <= self.tol * (1.0 + x_new.norm2()):
                 return x_new
             x = x_new
         raise ConvergenceError("Jacobi did not converge within max_iter")
 
+    def trace(self):
+        print("Jacobi Iteration Trace")
+        print(f"{'iter':>6} | {'residual norm':>14}")
+        print("-" * 26)
+        for k, res in enumerate(self.history):
+            print(f"{k:6d} | {res:14.6e}")
+
 
 class GaussSeidel:
-    def __init__(self, A: Matrix, b: Vector, tol: float = 1e-10, max_iter: int = 10_000):
+    """Gauss-Seidel iterative method for Ax = b."""
+
+    def __init__(
+        self, A: Matrix, b: Vector, tol: float = 1e-10, max_iter: int = 10_000
+    ):
         if not A.is_square():
             raise NonSquareMatrixError("A must be square")
         if A.n != len(b):
@@ -100,52 +159,75 @@ class GaussSeidel:
         self.b = b.copy()
         self.tol = tol
         self.max_iter = max_iter
+        self.history: list[float] = []
 
     def solve(self, x0: Vector | None = None) -> Vector:
         n = self.A.n
-        x = Vector([0.0]*n) if x0 is None else x0.copy()
+        x = Vector([0.0] * n) if x0 is None else x0.copy()
         for _ in range(self.max_iter):
             x_old = x.copy()
             for i in range(n):
                 diag = self.A.data[i][i]
                 if abs(diag) < 1e-15:
                     raise SingularMatrixError("Zero diagonal entry in Gauss-Seidel")
-                s1 = sum(self.A.data[i][j]*x[j] for j in range(i))
-                s2 = sum(self.A.data[i][j]*x_old[j] for j in range(i+1, n))
+                s1 = sum(self.A.data[i][j] * x[j] for j in range(i))
+                s2 = sum(self.A.data[i][j] * x_old[j] for j in range(i + 1, n))
                 x[i] = (self.b[i] - s1 - s2) / diag
-            if (x - x_old).norm_inf() <= self.tol * (1.0 + x.norm_inf()):
+            r = (self.A @ x) - self.b
+            res_norm = r.norm2()
+            self.history.append(res_norm)
+            if res_norm <= self.tol * (1.0 + x.norm2()):
                 return x
         raise ConvergenceError("Gauss-Seidel did not converge within max_iter")
 
+    def trace(self):
+        print("Gauss-Seidel Iteration Trace")
+        print(f"{'iter':>6} | {'residual norm':>14}")
+        print("-" * 26)
+        for k, res in enumerate(self.history):
+            print(f"{k:6d} | {res:14.6e}")
+
 
 class Cholesky:
-    """A = L L^T for SPD matrices."""
+    """Cholesky factorization A = L L^T for SPD matrices."""
+
     def __init__(self, A: Matrix):
         if not A.is_square():
             raise NonSquareMatrixError("A must be square")
         n = A.n
         for i in range(n):
-            for j in range(i+1, n):
+            for j in range(i + 1, n):
                 if abs(A.data[i][j] - A.data[j][i]) > 1e-12:
                     raise NotSymmetricError("Matrix is not symmetric")
         self.n = n
         self.L = Matrix.zeros(n, n)
+        self.steps: list[tuple[int, Matrix]] = []
         self._decompose(A.copy())
 
     def _decompose(self, A: Matrix) -> None:
         n = self.n
         for i in range(n):
-            for j in range(i+1):
-                s = sum(self.L.data[i][k]*self.L.data[j][k] for k in range(j))
+            for j in range(i + 1):
+                s = sum(self.L.data[i][k] * self.L.data[j][k] for k in range(j))
                 if i == j:
                     val = A.data[i][i] - s
                     if val <= 0.0:
-                        raise NotPositiveDefiniteError("Matrix is not positive definite")
+                        raise NotPositiveDefiniteError(
+                            "Matrix is not positive definite"
+                        )
                     self.L.data[i][j] = val**0.5
                 else:
                     self.L.data[i][j] = (A.data[i][j] - s) / self.L.data[j][j]
+            # record after each row i
+            self.steps.append((i, self.L.copy()))
 
     def solve(self, b: Vector) -> Vector:
         y = forward_substitution(self.L, b)
         x = backward_substitution(self.L.transpose(), y)
         return x
+
+    def trace(self):
+        print("Cholesky Decomposition Trace")
+        for i, L in self.steps:
+            print(f"\nRow {i}:")
+            print(f"L = {L}")

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/tutorial1_vectors_matrices.ipynb

@@ -0,0 +1,476 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "a9d0ab82",
+   "metadata": {},
+   "source": [
+    "# Tutorial 1 - Vectors & Matrices in `numethods`\n",
+    "\n",
+    "---\n",
+    "\n",
+    "## 1. Introduction\n",
+    "\n",
+    "In numerical computing, **vectors** and **matrices** are the basic objects.  \n",
+    "Almost every algorithm (linear solvers, eigenvalue problems, optimization, curve fitting, etc.) is built upon them.  \n",
+    "\n",
+    "In this tutorial, we will:\n",
+    "- Define what vectors and matrices are.\n",
+    "- Review their basic operations (addition, scalar multiplication, multiplication).\n",
+    "- Define and compute vector and matrix **norms**.\n",
+    "- Show how these operations work in the `numethods` package."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b91edc5c",
+   "metadata": {},
+   "source": [
+    "## 2. Importing our package"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "2107103a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from numethods.linalg import Vector, Matrix"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f64a4d9a",
+   "metadata": {},
+   "source": [
+    "## 3. What is a vector?\n",
+    "\n",
+    "A **vector** is an ordered collection of numbers (scalars).  \n",
+    "We can think of a vector as a column:\n",
+    "\n",
+    "$$\n",
+    "v =\n",
+    "\\begin{bmatrix}\n",
+    "v_1 \\\\\n",
+    "v_2 \\\\\n",
+    "\\vdots \\\\\n",
+    "v_n\n",
+    "\\end{bmatrix}, \\quad v \\in \\mathbb{R}^n\n",
+    "$$\n",
+    "\n",
+    "or as a row:\n",
+    "\n",
+    "$$\n",
+    "v^T = \\begin{bmatrix} v_1 & v_2 & \\cdots & v_n \\end{bmatrix}.\n",
+    "$$\n",
+    "\n",
+    "### Example\n",
+    "A vector in $\\mathbb{R}^3$:\n",
+    "\n",
+    "$$\n",
+    "v = \\begin{bmatrix} 1 \\\\ 2 \\\\ 3 \\end{bmatrix}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "46c96e42",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "v = Vector([1.0, 2.0, 3.0])\n"
+     ]
+    }
+   ],
+   "source": [
+    "v = Vector([1, 2, 3])\n",
+    "print(\"v =\", v)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5e951a1b",
+   "metadata": {},
+   "source": [
+    "## 4. What is a matrix?\n",
+    "\n",
+    "A **matrix** is a rectangular array of numbers with rows and columns:\n",
+    "\n",
+    "$$\n",
+    "A =\n",
+    "\\begin{bmatrix}\n",
+    "a_{11} & a_{12} & \\cdots & a_{1n} \\\\\n",
+    "a_{21} & a_{22} & \\cdots & a_{2n} \\\\\n",
+    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
+    "a_{m1} & a_{m2} & \\cdots & a_{mn}\n",
+    "\\end{bmatrix}, \\quad A \\in \\mathbb{R}^{m \\times n}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "a092da2f",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "A =\n",
+      " Matrix([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]])\n"
+     ]
+    }
+   ],
+   "source": [
+    "A = Matrix([[1, 2, 3],\n",
+    "            [4, 5, 6],\n",
+    "            [7, 8, 9]])\n",
+    "print(\"A =\\n\", A)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4f0d9ec4",
+   "metadata": {},
+   "source": [
+    "## 5. Basic operations\n",
+    "\n",
+    "### 5.1 Vector addition and subtraction\n",
+    "For $u, v \\in \\mathbb{R}^n$:\n",
+    "\n",
+    "$$\n",
+    "u + v = \\begin{bmatrix} u_1 + v_1 \\\\ u_2 + v_2 \\\\ \\vdots \\\\ u_n + v_n \\end{bmatrix},\n",
+    "\\quad\n",
+    "u - v = \\begin{bmatrix} u_1 - v_1 \\\\ u_2 - v_2 \\\\ \\vdots \\\\ u_n - v_n \\end{bmatrix}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "b7571e71",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "u + v = Vector([4.0, 4.0, 4.0])\n",
+      "u - v = Vector([2.0, 0.0, -2.0])\n"
+     ]
+    }
+   ],
+   "source": [
+    "u = Vector([3, 2, 1])\n",
+    "v = Vector([1, 2, 3])\n",
+    "\n",
+    "print(\"u + v =\", u + v)\n",
+    "print(\"u - v =\", u - v)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1aa8f396",
+   "metadata": {},
+   "source": [
+    "### 5.2 Scalar multiplication\n",
+    "For $\\alpha \\in \\mathbb{R}, v \\in \\mathbb{R}^n$:\n",
+    "\n",
+    "$$\n",
+    "\\alpha v = \\begin{bmatrix} \\alpha v_1 \\\\ \\alpha v_2 \\\\ \\vdots \\\\ \\alpha v_n \\end{bmatrix}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "b4168dfe",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "2 * v = Vector([2.0, 4.0, 6.0])\n",
+      "v * 2 = Vector([2.0, 4.0, 6.0])\n"
+     ]
+    }
+   ],
+   "source": [
+    "v = Vector([1, 2, 3])\n",
+    "\n",
+    "print(\"2 * v =\", 2 * v)\n",
+    "print(\"v * 2 =\", v * 2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e4919b59",
+   "metadata": {},
+   "source": [
+    "### 5.3 Matrix addition and subtraction\n",
+    "For $A, B \\in \\mathbb{R}^{m \\times n}$:\n",
+    "\n",
+    "$$\n",
+    "A + B = [ a_{ij} + b_{ij} ], \\quad\n",
+    "A - B = [ a_{ij} - b_{ij} ].\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "bd544880",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "A + B =\n",
+      " Matrix([[6.0, 8.0], [10.0, 12.0]])\n",
+      "A - B =\n",
+      " Matrix([[-4.0, -4.0], [-4.0, -4.0]])\n"
+     ]
+    }
+   ],
+   "source": [
+    "A = Matrix([[1, 2], [3, 4]])\n",
+    "B = Matrix([[5, 6], [7, 8]])\n",
+    "\n",
+    "print(\"A + B =\\n\", A + B)\n",
+    "print(\"A - B =\\n\", A - B)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0a1fb7f1",
+   "metadata": {},
+   "source": [
+    "### 5.4 Matrix-Vector multiplication\n",
+    "For $A \\in \\mathbb{R}^{m \\times n}, v \\in \\mathbb{R}^n$:\n",
+    "\n",
+    "$$\n",
+    "(Av)_i = \\sum_{j=1}^n a_{ij} v_j.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "d65fd20e",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "A @ v = Vector([-2.0, -2.0, -2.0])\n"
+     ]
+    }
+   ],
+   "source": [
+    "A = Matrix([[1, 2, 3],\n",
+    "            [4, 5, 6],\n",
+    "            [7, 8, 9]])\n",
+    "v = Vector([1, 0, -1])\n",
+    "\n",
+    "print(\"A @ v =\", A @ v)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8c2d30ec",
+   "metadata": {},
+   "source": [
+    "### 5.5 Matrix-Matrix multiplication\n",
+    "For $A \\in \\mathbb{R}^{m \\times n}, B \\in \\mathbb{R}^{n \\times p}$:\n",
+    "\n",
+    "$$\n",
+    "(AB)_{ij} = \\sum_{k=1}^n a_{ik} b_{kj}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "f46a96ec",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "A @ B =\n",
+      " Matrix([[4.0, 4.0], [10.0, 8.0]])\n"
+     ]
+    }
+   ],
+   "source": [
+    "A = Matrix([[1, 2],\n",
+    "            [3, 4]])\n",
+    "B = Matrix([[2, 0],\n",
+    "            [1, 2]])\n",
+    "\n",
+    "print(\"A @ B =\\n\", A @ B)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "470e38cb",
+   "metadata": {},
+   "source": [
+    "### 5.6 Transpose\n",
+    "For $A \\in \\mathbb{R}^{m \\times n}$:\n",
+    "\n",
+    "$$\n",
+    "A^T_{ij} = A_{ji}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "42168bd4",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "A^T =\n",
+      " Matrix([[1.0, 4.0], [2.0, 5.0], [3.0, 6.0]])\n"
+     ]
+    }
+   ],
+   "source": [
+    "A = Matrix([[1, 2, 3],\n",
+    "            [4, 5, 6]])\n",
+    "\n",
+    "print(\"A^T =\\n\", A.T)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "226efb8f",
+   "metadata": {},
+   "source": [
+    "## 6. Norms\n",
+    "\n",
+    "Norms measure the **size** or **length** of vectors and matrices."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "120cf212",
+   "metadata": {},
+   "source": [
+    "### 6.1 Vector norms\n",
+    "- **1-norm**:  $\\|v\\|_1 = \\sum |v_i|$  \n",
+    "- **2-norm (Euclidean)**: $\\|v\\|_2 = \\sqrt{\\sum v_i^2}$  \n",
+    "- **∞-norm**: $\\|v\\|_\\infty = \\max |v_i|$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "id": "3b187412",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "‖v‖₁ = 12.0\n",
+      "‖v‖₂ = 7.0710678118654755\n",
+      "‖v‖∞ = 5.0\n"
+     ]
+    }
+   ],
+   "source": [
+    "v = Vector([3, -4, 5])\n",
+    "\n",
+    "print(\"‖v‖₁ =\", v.norm())\n",
+    "print(\"‖v‖₂ =\", v.norm2())\n",
+    "print(\"‖v‖∞ =\", v.norm_inf())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3ad8369e",
+   "metadata": {},
+   "source": [
+    "### 6.2 Matrix norms\n",
+    "- **Frobenius norm**: $\\|A\\|_F = \\sqrt{\\sum_{i,j} a_{ij}^2}$  \n",
+    "- **1-norm** (maximum column sum): $\\|A\\|_1 = \\max_j \\sum_i |a_{ij}|$  \n",
+    "- **∞-norm** (maximum row sum): $\\|A\\|_\\infty = \\max_i \\sum_j |a_{ij}|$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "72a8ae7c",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "‖A‖_F = 5.477225575051661\n",
+      "‖A‖₁ = 6.0\n",
+      "‖A‖∞ = 7.0\n"
+     ]
+    }
+   ],
+   "source": [
+    "A = Matrix([[1, -2],\n",
+    "            [3,  4]])\n",
+    "\n",
+    "print(\"‖A‖_F =\", A.norm_fro())\n",
+    "print(\"‖A‖₁ =\", A.norm())\n",
+    "print(\"‖A‖∞ =\", A.norm_inf())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "eab94be3",
+   "metadata": {},
+   "source": [
+    "## 7. Summary\n",
+    "- A **vector** is an element of $\\mathbb{R}^n$, a list of numbers.  \n",
+    "- A **matrix** is a rectangular array of numbers $\\mathbb{R}^{m \\times n}$.  \n",
+    "- We can add, subtract, and scale vectors/matrices.  \n",
+    "- Multiplication extends naturally: matrix-vector and matrix-matrix.  \n",
+    "- **Transpose** flips rows and columns.  \n",
+    "- **Norms** measure the size of vectors and matrices."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

tutorials/tutorial3_orthogonalization.ipynb

@@ -0,0 +1,309 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "b2277afb",
+   "metadata": {},
+   "source": [
+    "# Tutorial 3 - Orthogonalization & QR Decomposition"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1e4a68da",
+   "metadata": {},
+   "source": [
+    "\n",
+    "In this tutorial, we study **orthogonalization** methods and the **QR decomposition**, which are central to numerical linear algebra.\n",
+    "\n",
+    "We will cover:\n",
+    "\n",
+    "- Orthogonal and orthonormal vectors\n",
+    "- QR decomposition\n",
+    "- Classical Gram-Schmidt\n",
+    "- Modified Gram-Schmidt\n",
+    "- Householder transformations\n",
+    "- Applications: least squares\n",
+    "- Examples with the `numethods` package\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "515dc61f",
+   "metadata": {},
+   "source": [
+    "\n",
+    "## 1. Why Orthogonalization?\n",
+    "\n",
+    "- Orthogonal vectors are easier to work with numerically.\n",
+    "- Many algorithms are more **stable** when using orthogonal bases.\n",
+    "- Key applications:\n",
+    "  - Solving **least squares problems**\n",
+    "  - Computing **eigenvalues**\n",
+    "  - Ensuring numerical stability in projections\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "af61a18f",
+   "metadata": {},
+   "source": [
+    "\n",
+    "## 2. Definitions\n",
+    "\n",
+    "### Orthogonal and Orthonormal vectors\n",
+    "\n",
+    "Two vectors $u, v \\in \\mathbb{R}^n$ are **orthogonal** if\n",
+    "\n",
+    "$$ u \\cdot v = 0. $$\n",
+    "\n",
+    "A set of vectors $\\{q_1, \\dots, q_m\\}$ is **orthonormal** if\n",
+    "\n",
+    "$$ q_i \\cdot q_j = \\begin{cases} 1 & i = j, \\\\ 0 & i \\neq j. \\end{cases} $$\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4f2330f1",
+   "metadata": {},
+   "source": [
+    "\n",
+    "### QR Decomposition\n",
+    "\n",
+    "For any $A \\in \\mathbb{R}^{m \\times n}$ with linearly independent columns, we can write\n",
+    "\n",
+    "$$ A = QR, $$\n",
+    "\n",
+    "- $Q \\in \\mathbb{R}^{m \\times n}$ has orthonormal columns ($Q^T Q = I$)\n",
+    "- $R \\in \\mathbb{R}^{n \\times n}$ is upper triangular\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "66567b83",
+   "metadata": {},
+   "source": [
+    "\n",
+    "## 3. Gram-Schmidt Orthogonalization\n",
+    "\n",
+    "### Classical Gram-Schmidt (CGS)\n",
+    "\n",
+    "Given linearly independent vectors $a_1, \\dots, a_n$:\n",
+    "\n",
+    "$$\n",
+    "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\\|}, \\qquad k = 2, \\ldots, n\n",
+    "$$\n",
+    "\n",
+    "Matrix form:\n",
+    "\n",
+    "$$ A = QR, \\quad R_{jk} = q_j^T a_k. $$\n",
+    "\n",
+    "⚠️ CGS can lose orthogonality in finite precision arithmetic.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "3fe97ce3",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Q (CGS): Matrix([[0.8164965809277261, -0.5520524474738834], [0.4082482904638631, 0.7590721152765896], [0.4082482904638631, 0.34503277967117707]])\n",
+      "R (CGS): Matrix([[2.449489742783178, 0.4082482904638631], [0.0, 2.41522945769824]])\n"
+     ]
+    }
+   ],
+   "source": [
+    "from numethods import Matrix, Vector\n",
+    "from numethods import QRGramSchmidt, QRModifiedGramSchmidt, QRHouseholder, LeastSquaresSolver\n",
+    "\n",
+    "# Example matrix\n",
+    "A = Matrix([[2, -1], [1, 2], [1, 1]])\n",
+    "\n",
+    "# Classical Gram-Schmidt\n",
+    "qrg = QRGramSchmidt(A)\n",
+    "print(\"Q (CGS):\", qrg.Q)\n",
+    "print(\"R (CGS):\", qrg.R)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ba84b59a",
+   "metadata": {},
+   "source": [
+    "\n",
+    "### Modified Gram-Schmidt (MGS)\n",
+    "\n",
+    "Same idea, but orthogonalization is done step by step:\n",
+    "\n",
+    "```\n",
+    "for k = 1..n:\n",
+    "    q_k = a_k\n",
+    "    for j = 1..k-1:\n",
+    "        r_jk = q_j^T q_k\n",
+    "        q_k = q_k - r_jk q_j\n",
+    "    r_kk = ||q_k||\n",
+    "    q_k = q_k / r_kk\n",
+    "```\n",
+    "MGS is more stable than CGS.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "e01e25ff",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Q (MGS): Matrix([[0.8164965809277261, -0.5520524474738834], [0.4082482904638631, 0.7590721152765896], [0.4082482904638631, 0.34503277967117707]])\n",
+      "R (MGS): Matrix([[2.449489742783178, 0.4082482904638631], [0.0, 2.41522945769824]])\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Modified Gram-Schmidt\n",
+    "qrm = QRModifiedGramSchmidt(A)\n",
+    "print(\"Q (MGS):\", qrm.Q)\n",
+    "print(\"R (MGS):\", qrm.R)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d893d189",
+   "metadata": {},
+   "source": [
+    "\n",
+    "## 4. Householder Reflections\n",
+    "\n",
+    "A more stable method uses **Householder matrices**.\n",
+    "\n",
+    "For a vector $x \\in \\mathbb{R}^m$:\n",
+    "\n",
+    "$$\n",
+    "v = x \\pm \\|x\\| e_1, \\quad H = I - 2 \\frac{vv^T}{v^T v}.\n",
+    "$$\n",
+    "\n",
+    "- $H$ is orthogonal ($H^T H = I$).\n",
+    "- Applying $H$ zeros out all but the first component of $x$.\n",
+    "\n",
+    "QR via Householder:\n",
+    "\n",
+    "$$\n",
+    "R = H_n H_{n-1} \\cdots H_1 A, \\quad Q = H_1^T H_2^T \\cdots H_n^T.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "15dfc35c",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Q (Householder): Matrix([[-0.8164965809277258, 0.552052447473883, -0.16903085094570333], [-0.40824829046386296, -0.7590721152765892, -0.5070925528371099], [-0.40824829046386296, -0.34503277967117707, 0.8451542547285166]])\n",
+      "R (Householder): Matrix([[-2.449489742783177, -0.408248290463863], [2.220446049250313e-16, -2.415229457698238], [2.220446049250313e-16, 2.220446049250313e-16]])\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Householder QR\n",
+    "qrh = QRHouseholder(A)\n",
+    "print(\"Q (Householder):\", qrh.Q)\n",
+    "print(\"R (Householder):\", qrh.R)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2b6c612f",
+   "metadata": {},
+   "source": [
+    "\n",
+    "## 5. Applications of QR\n",
+    "\n",
+    "### Least Squares\n",
+    "\n",
+    "We want to solve\n",
+    "\n",
+    "$$ \\min_x \\Vert Ax - b \\Vert_2^2. $$\n",
+    "\n",
+    "If $A = QR$, then\n",
+    "\n",
+    "$$ \\min_x \\Vert Ax - b \\Vert_2^2 = \\min_x \\Vert QRx - b \\Vert_2^2. $$\n",
+    "\n",
+    "Since $Q$ has orthonormal columns, and the normal equations boils down to\n",
+    "\n",
+    "$$ R x = Q^T b, $$\n",
+    "\n",
+    "we can therefore solve for $x$ by using back-substitution.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "25b399b7",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Least squares solution: Vector([1.0285714285714287, 0.828571428571429])\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Least squares example\n",
+    "b = Vector([1, 2, 3])\n",
+    "x_ls = LeastSquaresSolver(A, b).solve()\n",
+    "print(\"Least squares solution:\", x_ls)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a94c88c8",
+   "metadata": {},
+   "source": [
+    "\n",
+    "## 6. Key Takeaways\n",
+    "\n",
+    "- CGS is simple but numerically unstable.\n",
+    "- MGS is more stable and preferred if using Gram-Schmidt.\n",
+    "- Householder QR is the standard in practice (stable and efficient).\n",
+    "- QR decomposition underlies least squares, eigenvalue methods, and more.\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.13.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

tutorials/tutorial4_root_finding.ipynb

@@ -0,0 +1,366 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "3c51341a",
+   "metadata": {},
+   "source": [
+    "# Tutorial 4: Root-Finding Methods\n",
+    "---\n",
+    "In this tutorial, we study classical **root-finding algorithms** for nonlinear equations. We will:\n",
+    "\n",
+    "- Define the root-finding problem mathematically\n",
+    "- Derive several algorithms (bisection, fixed-point, Newton, secant)\n",
+    "- Discuss convergence conditions and error behavior\n",
+    "- Compare methods with worked examples using the `numethods` package.\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2999aa4f",
+   "metadata": {},
+   "source": [
+    "## 1. Problem Setup and Notation\n",
+    "\n",
+    "We seek to solve a nonlinear scalar equation:\n",
+    "$$\n",
+    "f(x) = 0, \\quad f: \\mathbb{R} \\to \\mathbb{R}.\n",
+    "$$\n",
+    "with a continuously differentiable function.\n",
+    "\n",
+    "\n",
+    "### Root, residual, and error\n",
+    "- A **root** $r$ satisfies $f(r)=0$.\n",
+    "- **Absolute error:** $(e_k = |x_k - x^\\star|)$.\n",
+    "- **Residual:** $(r_k = |f(x_k)|)$. \n",
+    "Note that small residual does not always imply small error.\n",
+    "\n",
+    "### Multiplicity\n",
+    "A root $r$ has **multiplicity** $m$ if\n",
+    "$$\n",
+    "f(r) = f'(r) = \\dots = f^{(m-1)}(r) = 0, \\quad f^{(m)}(r) \\neq 0.\n",
+    "$$\n",
+    "If $x^\\star$ satisfies $f(x^\\star)=0$ and $f'(x^\\star)\\ne 0$, we say the root is **simple** (multiplicity 1).\n",
+    "\n",
+    "If $f'(x^\\star)=\\cdots=f^{(m-1)}(x^\\star)=0$ and $f^{(m)}(x^\\star)\\ne 0$, we say the root has **multiplicity** (m).\n",
+    "\n",
+    "- **Simple roots** ($m=1$): most methods converge rapidly.\n",
+    "- **Multiple roots** ($m>1$): convergence often slows.\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b610cea2",
+   "metadata": {},
+   "source": [
+    "## 2. Bisection Method\n",
+    "\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",
+    "- Repeatedly bisect interval and select subinterval containing the root.\n",
+    "\n",
+    "**Iteration:**\n",
+    "$$\n",
+    "c_k = \\frac{a_k+b_k}{2}, \\quad f(c_k).\n",
+    "$$\n",
+    "\n",
+    "**Error bound:** interval length halves each step:\n",
+    "$$\n",
+    "|c_k-r| \\le \\frac{b-a}{2^k}.\n",
+    "$$\n",
+    "- Convergence: **linear**, guaranteed.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d26f75fb",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Bisection root: 1.4142135605216026\n",
+      "\n",
+      "Bisection Method Trace (x^2 - 2):\n",
+      "      iter |          a |          b |          c |       f(a) |       f(b) |       f(c) |   interval\n",
+      "--------------------------------------------------------------------------------------------------------\n",
+      "         0 |          0 |          2 |          1 |         -2 |          2 |         -1 |          2\n",
+      "         1 |          1 |          2 |        1.5 |         -1 |          2 |       0.25 |          1\n",
+      "         2 |          1 |        1.5 |       1.25 |         -1 |       0.25 |    -0.4375 |        0.5\n",
+      "         3 |       1.25 |        1.5 |      1.375 |    -0.4375 |       0.25 |  -0.109375 |       0.25\n",
+      "         4 |      1.375 |        1.5 |     1.4375 |  -0.109375 |       0.25 |  0.0664062 |      0.125\n",
+      "         5 |      1.375 |     1.4375 |    1.40625 |  -0.109375 |  0.0664062 | -0.0224609 |     0.0625\n",
+      "         6 |    1.40625 |     1.4375 |    1.42188 | -0.0224609 |  0.0664062 |  0.0217285 |    0.03125\n",
+      "         7 |    1.40625 |    1.42188 |    1.41406 | -0.0224609 |  0.0217285 | -0.000427246 |   0.015625\n",
+      "         8 |    1.41406 |    1.42188 |    1.41797 | -0.000427246 |  0.0217285 |  0.0106354 |  0.0078125\n",
+      "         9 |    1.41406 |    1.41797 |    1.41602 | -0.000427246 |  0.0106354 | 0.00510025 | 0.00390625\n",
+      "        10 |    1.41406 |    1.41602 |    1.41504 | -0.000427246 | 0.00510025 | 0.00233555 | 0.00195312\n",
+      "        11 |    1.41406 |    1.41504 |    1.41455 | -0.000427246 | 0.00233555 | 0.000953913 | 0.000976562\n",
+      "        12 |    1.41406 |    1.41455 |    1.41431 | -0.000427246 | 0.000953913 | 0.000263274 | 0.000488281\n",
+      "        13 |    1.41406 |    1.41431 |    1.41418 | -0.000427246 | 0.000263274 | -8.20011e-05 | 0.000244141\n",
+      "        14 |    1.41418 |    1.41431 |    1.41425 | -8.20011e-05 | 0.000263274 | 9.06326e-05 | 0.00012207\n",
+      "        15 |    1.41418 |    1.41425 |    1.41422 | -8.20011e-05 | 9.06326e-05 | 4.31482e-06 | 6.10352e-05\n",
+      "        16 |    1.41418 |    1.41422 |     1.4142 | -8.20011e-05 | 4.31482e-06 | -3.88434e-05 | 3.05176e-05\n",
+      "        17 |     1.4142 |    1.41422 |    1.41421 | -3.88434e-05 | 4.31482e-06 | -1.72643e-05 | 1.52588e-05\n",
+      "        18 |    1.41421 |    1.41422 |    1.41421 | -1.72643e-05 | 4.31482e-06 | -6.47477e-06 | 7.62939e-06\n",
+      "        19 |    1.41421 |    1.41422 |    1.41421 | -6.47477e-06 | 4.31482e-06 | -1.07998e-06 | 3.8147e-06\n",
+      "        20 |    1.41421 |    1.41422 |    1.41421 | -1.07998e-06 | 4.31482e-06 | 1.61742e-06 | 1.90735e-06\n",
+      "        21 |    1.41421 |    1.41421 |    1.41421 | -1.07998e-06 | 1.61742e-06 | 2.68718e-07 | 9.53674e-07\n",
+      "        22 |    1.41421 |    1.41421 |    1.41421 | -1.07998e-06 | 2.68718e-07 | -4.05632e-07 | 4.76837e-07\n",
+      "        23 |    1.41421 |    1.41421 |    1.41421 | -4.05632e-07 | 2.68718e-07 | -6.84571e-08 | 2.38419e-07\n",
+      "        24 |    1.41421 |    1.41421 |    1.41421 | -6.84571e-08 | 2.68718e-07 | 1.0013e-07 | 1.19209e-07\n",
+      "        25 |    1.41421 |    1.41421 |    1.41421 | -6.84571e-08 | 1.0013e-07 | 1.58366e-08 | 5.96046e-08\n",
+      "        26 |    1.41421 |    1.41421 |    1.41421 | -6.84571e-08 | 1.58366e-08 | -2.63102e-08 | 2.98023e-08\n",
+      "        27 |    1.41421 |    1.41421 |    1.41421 | -2.63102e-08 | 1.58366e-08 | -5.23681e-09 | 1.49012e-08\n"
+     ]
+    }
+   ],
+   "source": [
+    "from numethods import Bisection, FixedPoint, NewtonRoot, Secant, print_trace\n",
+    "import math\n",
+    "\n",
+    "# Example function: f(x) = x^2 - 2\n",
+    "f = lambda x: x**2 - 2\n",
+    "df = lambda x: 2*x\n",
+    "\n",
+    "# Bisection\n",
+    "bisect = Bisection(f, 0, 2, tol=1e-8)\n",
+    "root_b = bisect.solve()\n",
+    "print('Bisection root:', root_b)\n",
+    "\n",
+    "steps = bisect.trace()\n",
+    "print(\"\\nBisection Method Trace (x^2 - 2):\")\n",
+    "print_trace(steps)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3be4a510",
+   "metadata": {},
+   "source": [
+    "## 3. Fixed-Point Iteration\n",
+    "- Rewrite equation as $x=g(x)$.\n",
+    "- Iterate:\n",
+    "$$\n",
+    "x_{k+1} = g(x_k).\n",
+    "$$\n",
+    "\n",
+    "**Convergence theorem (Banach fixed-point):** If g is continuously differentiable near $(x^\\star)$ and\n",
+    "$$\n",
+    "|g'(x_\\star)| < 1,\n",
+    "$$\n",
+    "then for initial guesses $x_0$ sufficiently close to $x^\\star$, the iterates converge **linearly** to $x^\\star$ with asymptotic rate $|g'(x^\\star)|$.\n",
+    "\n",
+    "**Choice of g.** Different rearrangements yield different g's with different convergence properties.\n",
+    "A poor choice (with $(|g'|\\ge 1))$ can diverge.\n",
+    "\n",
+    "- Rate: linear with factor $|g'(r)|$.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "436ce6f6",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Fixed-point failed: Fixed-point iteration did not converge\n",
+      "\n",
+      "Fixed-Point Iteration Trace (x^2 - 2):\n",
+      "      iter |          x |      x_new |      error\n",
+      "----------------------------------------------------\n",
+      "         0 |          1 |        1.5 |        0.5\n",
+      "         1 |        1.5 |    1.41667 |  0.0833333\n",
+      "         2 |    1.41667 |    1.41422 | 0.00245098\n",
+      "         3 |    1.41422 |    1.41421 | 2.1239e-06\n",
+      "         4 |    1.41421 |    1.41421 | 1.59495e-12\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Fixed point: g(x)=sqrt(2) ~ rewriting\n",
+    "g = lambda x: (2/x)  # not always convergent, but demonstrates\n",
+    "try:\n",
+    "    fp = FixedPoint(g, 1.0, tol=1e-8)\n",
+    "    root_fp = fp.solve()\n",
+    "    print('Fixed-point root:', root_fp)\n",
+    "except Exception as e:\n",
+    "    print('Fixed-point failed:', e)\n",
+    "\n",
+    "g = lambda x: 0.5 * (x + 2 / x)\n",
+    "steps = FixedPoint(g, 1.0).trace()\n",
+    "print(\"\\nFixed-Point Iteration Trace (x^2 - 2):\")\n",
+    "print_trace(steps)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "40e66a30",
+   "metadata": {},
+   "source": [
+    "## 4. Newton’s Method\n",
+    "From Taylor expansion:\n",
+    "$$\n",
+    "f(x) \\approx f(x_k) + f'(x_k)(x-x_k).\n",
+    "$$\n",
+    "Set $f(x)=0$ to solve for next iterate:\n",
+    "$$\n",
+    "x_{k+1} = x_k - \\frac{f(x_k)}{f'(x_k)}.\n",
+    "$$\n",
+    "\n",
+    "- **Quadratic convergence** for simple roots.\n",
+    "- For multiple roots: drops to linear.\n",
+    "- Requires derivative, sensitive to initial guess.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "7ebf9068",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Newton root: 1.4142135623730951\n",
+      "Newton Method Trace (x^2 - 2):\n",
+      "      iter |          x |       f(x) |      df(x) |      x_new |      error\n",
+      "------------------------------------------------------------------------------\n",
+      "         0 |          1 |         -1 |          2 |        1.5 |        0.5\n",
+      "         1 |        1.5 |       0.25 |          3 |    1.41667 |  0.0833333\n",
+      "         2 |    1.41667 | 0.00694444 |    2.83333 |    1.41422 | 0.00245098\n",
+      "         3 |    1.41422 | 6.0073e-06 |    2.82843 |    1.41421 | 2.1239e-06\n",
+      "         4 |    1.41421 | 4.51061e-12 |    2.82843 |    1.41421 | 1.59472e-12\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Newton\n",
+    "newton = NewtonRoot(f, df, 1.0, tol=1e-12)\n",
+    "root_n = newton.solve()\n",
+    "print('Newton root:', root_n)\n",
+    "\n",
+    "steps = newton.trace()\n",
+    "print(\"Newton Method Trace (x^2 - 2):\")\n",
+    "print_trace(steps)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "84888305",
+   "metadata": {},
+   "source": [
+    "## 5. Secant Method\n",
+    "- Avoids derivative by approximating slope with finite difference:\n",
+    "$$\n",
+    "x_{k+1} = x_k - f(x_k) \\frac{x_k - x_{k-1}}{f(x_k)-f(x_{k-1})}.\n",
+    "$$\n",
+    "\n",
+    "- Convergence order: $\\approx 1.618$ (superlinear).\n",
+    "- More efficient than Newton if derivative expensive.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "f2318bf3",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Secant root: 1.414213562373095\n",
+      "\n",
+      "Secant Method Trace (x^2 - 2):\n",
+      "      iter |         x0 |         x1 |         x2 |      f(x0) |      f(x1) |      error\n",
+      "-------------------------------------------------------------------------------------------\n",
+      "         0 |          0 |          2 |          1 |         -2 |          2 |          1\n",
+      "         1 |          2 |          1 |    1.33333 |          2 |         -1 |   0.333333\n",
+      "         2 |          1 |    1.33333 |    1.42857 |         -1 |  -0.222222 |  0.0952381\n",
+      "         3 |    1.33333 |    1.42857 |    1.41379 |  -0.222222 |  0.0408163 |  0.0147783\n",
+      "         4 |    1.42857 |    1.41379 |    1.41421 |  0.0408163 | -0.00118906 | 0.000418335\n",
+      "         5 |    1.41379 |    1.41421 |    1.41421 | -0.00118906 | -6.00729e-06 | 2.12421e-06\n",
+      "         6 |    1.41421 |    1.41421 |    1.41421 | -6.00729e-06 | 8.93146e-10 | 3.15775e-10\n",
+      "         7 |    1.41421 |    1.41421 |    1.41421 | 8.93146e-10 | -8.88178e-16 | 2.22045e-16\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Secant\n",
+    "sec = Secant(f, 0, 2, tol=1e-12)\n",
+    "root_s = sec.solve()\n",
+    "print('Secant root:', root_s)\n",
+    "\n",
+    "steps = sec.trace()\n",
+    "print(\"\\nSecant Method Trace (x^2 - 2):\")\n",
+    "print_trace(steps)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7ed8a0b5",
+   "metadata": {},
+   "source": [
+    "## 6. Stopping Criteria\n",
+    "We stop iteration when:\n",
+    "- $|f(x_k)| \\le \\varepsilon$ (residual small), or\n",
+    "- $|x_{k+1}-x_k| \\le \\varepsilon(1+|x_{k+1}|)$ (relative step small).\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "833b538c",
+   "metadata": {},
+   "source": [
+    "## 7. Comparison of Methods\n",
+    "| Method | Requires derivative | Convergence rate | Guarantee? |\n",
+    "|--------|---------------------|------------------|------------|\n",
+    "| Bisection | No | Linear | Yes (if sign change) |\n",
+    "| Fixed-Point | No | Linear | Not always |\n",
+    "| Newton | Yes | Quadratic | Locally (good guess) |\n",
+    "| Secant | No | ~1.618 (superlinear) | Locally (good guess) |\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0d372aa0",
+   "metadata": {},
+   "source": [
+    "## 8. Exercises\n",
+    "1. Apply all four methods to $f(x)=\\cos x - x$.\n",
+    "2. Try Newton’s method on $f(x)=(x-1)^2$ and compare convergence rate with $f(x)=x^2-2$.\n",
+    "3. Modify Secant to stop if denominator too small.\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.13.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}