update eigen.py

fixed -> matrix and vector imported twice from linalg

.gitignore

@@ -0,0 +1,366 @@
+# ---> JupyterNotebooks
+# gitignore template for Jupyter Notebooks
+# website: http://jupyter.org/
+
+.ipynb_checkpoints
+*/.ipynb_checkpoints/*
+
+# IPython
+profile_default/
+ipython_config.py
+
+# Remove previous ipynb_checkpoints
+#   git rm -r .ipynb_checkpoints/
+
+# ---> VisualStudioCode
+.vscode/*
+!.vscode/settings.json
+!.vscode/tasks.json
+!.vscode/launch.json
+!.vscode/extensions.json
+!.vscode/*.code-snippets
+
+# Local History for Visual Studio Code
+.history/
+
+# Built Visual Studio Code Extensions
+*.vsix
+
+# ---> Python
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+cover/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+.pybuilder/
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+#   For a library or package, you might want to ignore these files since the code is
+#   intended to run in multiple environments; otherwise, check them in:
+# .python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# UV
+#   Similar to Pipfile.lock, it is generally recommended to include uv.lock in version control.
+#   This is especially recommended for binary packages to ensure reproducibility, and is more
+#   commonly ignored for libraries.
+#uv.lock
+
+# poetry
+#   Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control.
+#   This is especially recommended for binary packages to ensure reproducibility, and is more
+#   commonly ignored for libraries.
+#   https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control
+#poetry.lock
+
+# pdm
+#   Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control.
+#pdm.lock
+#   pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it
+#   in version control.
+#   https://pdm.fming.dev/latest/usage/project/#working-with-version-control
+.pdm.toml
+.pdm-python
+.pdm-build/
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+# pytype static type analyzer
+.pytype/
+
+# Cython debug symbols
+cython_debug/
+
+# PyCharm
+#  JetBrains specific template is maintained in a separate JetBrains.gitignore that can
+#  be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore
+#  and can be added to the global gitignore or merged into this file.  For a more nuclear
+#  option (not recommended) you can uncomment the following to ignore the entire idea folder.
+#.idea/
+
+# Ruff stuff:
+.ruff_cache/
+
+# PyPI configuration file
+.pypirc
+
+# ---> Julia
+# Files generated by invoking Julia with --code-coverage
+*.jl.cov
+*.jl.*.cov
+
+# Files generated by invoking Julia with --track-allocation
+*.jl.mem
+
+# System-specific files and directories generated by the BinaryProvider and BinDeps packages
+# They contain absolute paths specific to the host computer, and so should not be committed
+deps/deps.jl
+deps/build.log
+deps/downloads/
+deps/usr/
+deps/src/
+
+# Build artifacts for creating documentation generated by the Documenter package
+docs/build/
+docs/site/
+
+# File generated by Pkg, the package manager, based on a corresponding Project.toml
+# It records a fixed state of all packages used by the project. As such, it should not be
+# committed for packages, but should be committed for applications that require a static
+# environment.
+Manifest.toml
+
+# ---> R
+# History files
+.Rhistory
+.Rapp.history
+
+# Session Data files
+.RData
+.RDataTmp
+
+# User-specific files
+.Ruserdata
+
+# Example code in package build process
+*-Ex.R
+
+# Output files from R CMD build
+/*.tar.gz
+
+# Output files from R CMD check
+/*.Rcheck/
+
+# RStudio files
+.Rproj.user/
+
+# produced vignettes
+vignettes/*.html
+vignettes/*.pdf
+
+# OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3
+.httr-oauth
+
+# knitr and R markdown default cache directories
+*_cache/
+/cache/
+
+# Temporary files created by R markdown
+*.utf8.md
+*.knit.md
+
+# R Environment Variables
+.Renviron
+
+# pkgdown site
+docs/
+
+# translation temp files
+po/*~
+
+# RStudio Connect folder
+rsconnect/
+
+# ---> MATLAB
+# Windows default autosave extension
+*.asv
+
+# OSX / *nix default autosave extension
+*.m~
+
+# Compiled MEX binaries (all platforms)
+*.mex*
+
+# Packaged app and toolbox files
+*.mlappinstall
+*.mltbx
+
+# Generated helpsearch folders
+helpsearch*/
+
+# Simulink code generation folders
+slprj/
+sccprj/
+
+# Matlab code generation folders
+codegen/
+
+# Simulink autosave extension
+*.autosave
+
+# Simulink cache files
+*.slxc
+
+# Octave session info
+octave-workspace
+
+# ---> VirtualEnv
+# Virtualenv
+# http://iamzed.com/2009/05/07/a-primer-on-virtualenv/
+.Python
+[Bb]in
+[Ii]nclude
+[Ll]ib
+[Ll]ib64
+[Ll]ocal
+[Ss]cripts
+pyvenv.cfg
+.venv
+pip-selfcheck.json
+
+# ---> macOS
+# General
+.DS_Store
+.AppleDouble
+.LSOverride
+
+# Icon must end with two \r
+Icon
+
+
+# Thumbnails
+._*
+
+# Files that might appear in the root of a volume
+.DocumentRevisions-V100
+.fseventsd
+.Spotlight-V100
+.TemporaryItems
+.Trashes
+.VolumeIcon.icns
+.com.apple.timemachine.donotpresent
+
+# Directories potentially created on remote AFP share
+.AppleDB
+.AppleDesktop
+Network Trash Folder
+Temporary Items
+.apdisk
+
+Powered by Gitea
+Version:
+1.24.5
+Page:
+152ms
+Template:
+8ms
+Licenses
+API

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,6 +1,23 @@
+import sys
+
+# sys.path.append("../")
 from numethods import *
 
+
+def demo_linear_solvers():
+    A = Matrix([[4, -1, 0], [-1, 4, -1], [0, -1, 3]])
+    b = Vector([15, 10, 10])
+
+    print("LU:", LUDecomposition(A).solve(b))
+    print("Gauss-Jordan:", GaussJordan(A).solve(b))
+    print("Cholesky:", Cholesky(A).solve(b))
+    print("Jacobi:", Jacobi(A, b, tol=1e-12).solve())
+    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
@@ -16,32 +33,37 @@ 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_linear_solvers():
-    A = Matrix([[4, -1, 0], [-1, 4, -1], [0, -1, 3]])
-    b = Vector([15, 10, 10])
-
-    print("LU:", LUDecomposition(A).solve(b))
-    print("Gauss-Jordan:", GaussJordan(A).solve(b))
-    print("Cholesky:", Cholesky(A).solve(b))
-    print("Jacobi:", Jacobi(A, b, tol=1e-12).solve())
-    print("Gauss-Seidel:", GaussSeidel(A, b, tol=1e-12).solve())
-
-
 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():
@@ -54,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,15 +98,104 @@ 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."""
+    def __matmul__(self, other: Union["Matrix", "Vector"]):
+        if isinstance(other, Matrix):
             if self.n != other.m:
-            raise ValueError("Matrix dimensions do not align for multiplication")
+                raise ValueError("dims")
             return Matrix(
                 [
                     [
@@ -107,12 +205,22 @@ class Matrix:
                     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)
@@ -26,9 +34,42 @@ class Bisection:
                 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
+        raise ConvergenceError("Bisection did not converge")
+
 
 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,16 +81,34 @@ 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
@@ -59,9 +118,42 @@ class Secant:
             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:
@@ -75,3 +167,45 @@ class NewtonRoot:
                 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]
+                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,6 +115,7 @@ 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
@@ -84,14 +129,28 @@ class Jacobi:
                 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,6 +159,7 @@ 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
@@ -113,13 +173,24 @@ class GaussSeidel:
                 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")
@@ -130,6 +201,7 @@ class Cholesky:
                     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:
@@ -140,12 +212,22 @@ class Cholesky:
                 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
+}