{
 "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
}