{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(section:linear-equations-pivoting)=\n",
    "# Partial Pivoting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Last revised on August 25, 2025**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**References:**\n",
    "\n",
    "- {cite}`Chasnov` Section 3.3, *Partial pivoting*.\n",
    "\n",
    "- {cite}`Sauer` Section 2.4.1, *Partial Pivoting*.\n",
    "\n",
    "- {cite}`Burden-Faires` Section 6.2, *Pivoting Stratgies*.\n",
    "\n",
    "- {cite}`Dionne` Section 4.1, *Gaussian Elimination with Backward Substitution*.\n",
    "\n",
    "- {cite}`Chenney-Kincaid` Section 2.2, *Gaussian Elimination with Scaled Partial Pivoting*.\n",
    "\n",
    "- {cite}`Trefethen-Bau` Lecture 21, *Pivoting*. (However, this again is done directly in the context of the LU factorization, as in section {ref}`section:linear-equations-plu-factorization`).\n",
    "\n",
    "- {cite}`Dahlquist-Bjorck-v2` Section 7.2.3, *Pivoting Strategies*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "The basic row reduction method can fail due to divisoion by zero (and to have very large rouding errors when a denominator is extremely close to zero.\n",
    "A more robust modification is to swap the order of the equations to avaid these problems: *partial pivotng*.\n",
    "\n",
    "Here we look at a particularly robust version of this strategy, *Maximal Element Partial Pivoting*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{prf:remark}\n",
    ":label: no-scaled-partial-pivoting\n",
    "\n",
    "Some references describe the method of *scaled* partial pivoting, but here we present instead a version without the \"scaling\", because not only is it simpler, but modern research shows that it is esentially always as good, once the problem is set up in a \"sane\" way.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Maximum Norm of a Vector\n",
    "\n",
    "To discuss errors in linear algebra, it is useful to have a measure of the \"size\" of a vector, such as in describing the differerence between an exact solution vector and an approximation. \n",
    "For this, the usual *Euclidean norm*\n",
    "\n",
    "$$\\| x \\| = \\| x \\|_2 = \\sqrt{\\sum_{i=1}^n |x_i|^2}$$\n",
    "\n",
    "is not always well-suited.\n",
    "Instead we will mostly use the **maximum norm** or **infinity norm**\n",
    "\n",
    "$$\n",
    "\\| x \\|_{max} = \\| x \\|_\\infty : = \\max_{i=1}^n |x_i|\n",
    "$$\n",
    "\n",
    "For more details — including an extension to measuring the size of a matrix and an explanation of that mysterious name *infinity norm* — see section {ref}`section:linear-equations-error-bounds`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What can go wrong with naive row reduction?\n",
    "\n",
    "We have noted two problems with the naive algorithm for row reduction: total failure due the division be zero, and loss of precision due to dividing by very small values — or more preciselt calculations the lead to intermediate values far larger than the final results.\n",
    "The culprits in all cases are the same: the denominators are first the *pivot elements* $a_{k,k}^{(k-1)}$ in evaluation of $l_{i,k}$ during row reduction and then the $u_{k,k}$ in back substitution.\n",
    "Further, those $a_{k,k}^{(k-1)}$ are the final updated values at indices $(k,k)$, so are the same as $u_{k,k}$.\n",
    "\n",
    "Thus it is exactly these *main diagonal elements* that we must deal with."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The basic fix: Partial Pivoting\n",
    "\n",
    "The basic strategy is that at step $k$, we can swap equation $k$ with any equation $i$, $i > k$.\n",
    "Note that this involves swapping those rows of array `A` and also those elements of the array `b` for the right-hand side: $b_k \\leftrightarrow b_i$.\n",
    "\n",
    "This approach of swapping equations (swapping rows in arrays `A` and `b`) is called **pivoting**, or more specifically **partial pivoting**, to distinguish from the more elaborate strategy where to columns of `A` are also reordered (which is equivalent to reordeting the unknowns in the equations).\n",
    "The row that is swapped with row $k$ is sometimes called the **pivot row**, and the new denominator is the corresponding **pivot element**.\n",
    "\n",
    "This approach is robust so long as one is using exact arithmetic: it works for any well-posed system because so long as the $Ax = b$ has a unique solution — so that the original matrix $A$ is non-singular — at least one of the $a_{i,k}^{(k-1)}, i \\geq k$ will be non-zero, and thus the swap will give a new element in position $(k,k)$ that is non-zero.\n",
    "(I will stop caring about superscripts to distinguish updates, but if you wish to, the elements of the new row $k$ could be called either $a_{k,j}^{(k)}$ or even $u_{k,j}$, since those values are in their final state.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Handling rounding error: Maximal Element Partial Pivoting\n",
    "\n",
    "The final refinement is to seek the smallest possible magnitudes for intermediate values, and thus the smallest absolute errors in them, by making the multipliers $l_{i,k}$ small, in turn by making the denominator $a_{k,k}^{(k-1)} = u_{k,k}$ as large as possible in magnitude:\n",
    "\n",
    "*At step $k$, choose the pivot row $p_k \\geq k$ so that $|a_{p_k,k}^{(k-1)}| \\geq |a_{i,k}^{(k-1)}|$ for all $i \\geq k$. If there is more that one such element of largest magnitude, use the lowest value:\n",
    "in particular, if $p_k = k$ works, use it and do not swap!*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A consequence of this is that the mutipliers used in the row operations all have $|l_{i,k}| = \\left| \\frac{a_{p_i,k}^{(k-1)}}{a_{p_k,k}^{(k-1)}} \\right| < 1.$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Swapping values and rows in Python\n",
    "\n",
    "In Julia (as in Python) the value of two variables `a` and `b` can be swapped via tuples with `(a, b) = (b, a)`,\n",
    "and combined with array slicing, this can also be used to swap rows (or columns)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 1)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = 1\n",
    "b = 2\n",
    "(a, b) = (b, a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       " 11  12  13\n",
       " 21  22  23\n",
       " 31  32  33"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [11 12 13 ; 21 22 23 ; 31 32 33]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       " 21  22  23\n",
       " 11  12  13\n",
       " 31  32  33"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(A[1,:], A[2,:]) = (A[2,:], A[1,:])\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "See [Exercise 1](#exercise-1)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some demonstrations\n",
    "\n",
    "No row reduction is done here, so entire rows are swapped rather than just the elements from column $k$ onward:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, get the matrix pretty-printer seen in {ref}`section:linear-equations-row-reduction`;\n",
    "this time from {ref}`module:NumericalMethods`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"NumericalMethods.jl\")\n",
    "using .NumericalMethods: printmatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initially A is\n",
      "[ 1 -6 2 \n",
      "  3 5 -6 \n",
      "  4 2 7 ]\n"
     ]
    }
   ],
   "source": [
    "A = [1 -6 2 ; 3 5 -6 ; 4 2 7]\n",
    "n = 3\n",
    "println(\"Initially A is\")\n",
    "printmatrix(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After swapping rows 1 <-> 3 using slicing and a temporary row, A is\n",
      "[ 4 2 7 \n",
      "  3 5 -6 \n",
      "  1 -6 2 ]\n"
     ]
    }
   ],
   "source": [
    "k = 1\n",
    "p = 3\n",
    "temp = copy(A[k,:])\n",
    "A[k,:] = A[p,:]\n",
    "A[p,:] = temp\n",
    "println(\"After swapping rows 1 <-> 3 using slicing and a temporary row, A is\")\n",
    "printmatrix(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After swapping rows 2 <-> 3 using a loop and tuples of elements (no temp) A is:\n",
      "[ 4 2 7 \n",
      "  1 -6 2 \n",
      "  3 5 -6 ]\n"
     ]
    }
   ],
   "source": [
    "k = 2\n",
    "p = 3\n",
    "for j in 1:n\n",
    "    ( A[k,j] , A[p,j] ) = ( A[p,j] , A[k,j] )\n",
    "end\n",
    "println(\"After swapping rows 2 <-> 3 using a loop and tuples of elements (no temp) A is:\")\n",
    "printmatrix(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After swapping rows 1 <-> 2 using tuples of slices (no loop or temp) A is:\n",
      "[ 1 -6 2 \n",
      "  4 2 7 \n",
      "  3 5 -6 ]\n"
     ]
    }
   ],
   "source": [
    "k = 1\n",
    "p = 2\n",
    "( A[k,:] , A[p,:] ) = ( copy(A[p,:]) , copy(A[k,:]) )\n",
    "println(\"After swapping rows 1 <-> 2 using tuples of slices (no loop or temp) A is:\")\n",
    "printmatrix(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## When is it safe to do without pivoting?\n",
    "\n",
    "{prf:ref}`theorem-row-reduction-preserves-sdd` shows that diagonal dominance guarantees that pivoting is not necessary because the diagonal elements are never zero.\n",
    "\n",
    "In the case of the matrix $A$ being\n",
    "column-wise SDD as in {prf:ref}`definition-columnwise-strictly-diagonally-dominant`,\n",
    "the situation is even better; there is no reason for pivoting:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{prf:theorem}\n",
    ":label: theorem-no-pivoting-columnwise-sdd\n",
    "\n",
    "If matrix $A$ is column-wise SDD, maximal element partial pivoting in fact does no row-swaps;\n",
    "it does the same thing as naive row reduction.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Being row-wise SDD is more \"natural\" and common than being column-wise SDD, because the former is a property \"within\" each of the equations that go into the matrix.\n",
    "\n",
    "This might seem unfortunate, but there is a way to get the benefits of the above nice result also for row-wise SDD matrices, which we will see in the section {ref}`section:linear-equations-lu-factorization` with the\n",
    "[Crout decomposition](crout-decomposition)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{exercise-start} On Julia coding\n",
    ":label: linear-equations-pivoting-exercise-1\n",
    "```\n",
    "\n",
    "Explain why we cannot just swap the relevant elements of rows $k$ and $p$ with:\n",
    "\n",
    "    for j in 1:n\n",
    "        A[k,j] = A[p,j]\n",
    "        A[p,j] = A[k,j]\n",
    "    end\n",
    "\n",
    "or in vectorized form:\n",
    "\n",
    "    A[k,:] = A[p,:]\n",
    "    A[p,:] = A[k,:]\n",
    "\n",
    "Describe what happens instead.\n",
    "\n",
    "```{exercise-end}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{exercise-start}\n",
    ":label: linear-equations-pivoting-exercise-2\n",
    "```\n",
    "Solve $A \\vec{x} = \\vec b$\n",
    "with\n",
    "\n",
    "$$\n",
    "A = \\left[ \\begin{array}{ccc} 4. & 2. & 1. \\\\ 9. & 3. & 1. \\\\ 25. & 5. & 1. \\end{array} \\right],\n",
    "\\quad\n",
    "\\vec b = \\left[ \\begin{array}{c} 0.693147  \\\\ 1.098612 \\\\ 1.609438 \\end{array} \\right]\n",
    "$$\n",
    "as seen in {ref}`linear-equations-row-reduction-exercise-1` of {ref}`section:linear-equations-row-reduction`,\n",
    "now using *maximal element partial pivoting,* computing all intermediate values as decimal approximations rounded to four significant digits -- no fractions!\n",
    "Call this approximate solution $\\vec{x}_2$.\n",
    "\n",
    "Then compute the residual $\\vec{r}_2 = \\vec{b} - A \\vec{x}_2$ and its norm $\\| \\vec{r}_2 \\|_\\infty$.\n",
    "This should be done with high accuracy (not rounding to four decimal places) and could be done using a Python notebook as a \"calculator\".\n",
    "\n",
    "Next, use Python software to compute the condition number of $A$, $K_\\infty(A)$,\n",
    "and use this to get a bound on the absolute error in each of the above approximations.\n",
    "\n",
    "Finally some \"cheating\": use Python software to compute the \"exact\" solution and use this to compute the actual absolute error;\n",
    "compare this to the bound computed above.\n",
    "\n",
    "```{exercise-end}\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.11.4",
   "language": "julia",
   "name": "julia-1.11"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}