{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Partial Pivoting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**References:**\n", "\n", "- Section 2.4.1 *Partial Pivoting* of {cite}`Sauer`.\n", "- Section 6.2 *Pivoting Stratgies* of {cite}`Burden-Faires`.\n", "- Section 7.1 of {cite}`Chenney-Kincaid`.\n", "\n", "```{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": [ "## 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", "Here we look at a particularly robust version of this strategy, *Maximal Element Partial Pivoting*." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# As in recent sections, we import some items from modules individually, so they can be used by \"first name only\".\n", "from numpy import array" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What can go wrong with Naive Gaussian Elimination?\n", "\n", "We have noted two problems with the naive algorithm for Gaussian elimination: 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", "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": [ "## Swapping rows in Python\n", "\n", "I will not give a detailed algorithm for this, since we will soon implement an even better variant.\n", "\n", "However, here are some notes on swapping values and how to avoid a possible pitfall." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1, on Python coding\n", "\n", "a) Explain why we cannot just swap the relevant elements of rows $k$ and $p$ with:\n", "\n", " for j in range(k,n):\n", " A[k,j] = A[p,j]\n", " A[p,j] = A[k,j]\n", "\n", "or with vectorized \"slicing\":\n", "\n", " A[k,k:] = A[p,k:]\n", " A[p,k:] = A[k,k:]\n", "\n", "Describe what happens instead.\n", " \n", "b) A common strategy to avoid this problem uses an intermediate temporary copy of each value being \"moved\".\n", "This can be combined with slicing, but be careful: arrays (including slices of arrays) must be copied with the method `.copy()`:\n", "\n", " temp = A[k,k:].copy()\n", " A[k,k:] = A[p,k:]\n", " A[p,k:] = temp\n", "\n", "3) Python also has an elegant alternative for swapping a pair of values:\n", "\n", " for j in range(k,n):\n", " ( A[k,j] , A[p,j] ) = ( A[p,j] , A[k,j] )\n", "\n", "This can also be done with slicing, with care to copy those slices:\n", "\n", " ( A[k,k:] , A[p,k:] ) = ( A[p,k:].copy() , A[k,k:].copy() )" ] }, { "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": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Initially,\n", "A =\n", " [[ 1. -6. 2.]\n", " [ 3. 5. -6.]\n", " [ 4. 2. 7.]]\n" ] } ], "source": [ "A = array([[1. , -6. , 2.],[3. , 5. , -6.],[4. , 2. , 7.]])\n", "n = 3\n", "print(f\"Initially,\\nA =\\n {A}\")" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "After swapping rows 1 <-> 3 (row indices 0 <-> 2) using slicing and a temporary row,\n", "A =\n", " [[ 4. 2. 7.]\n", " [ 3. 5. -6.]\n", " [ 1. -6. 2.]]\n" ] } ], "source": [ "k = 0\n", "p = 2\n", "temp = A[k,k:].copy()\n", "A[k,k:] = A[p,k:]\n", "A[p,k:] = temp\n", "print(\"After swapping rows 1 <-> 3 (row indices 0 <-> 2) using slicing and a temporary row,\")\n", "print(f\"A =\\n {A}\")" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "After swapping rows 2 <-> 3 using a loop and tuples of elements, no temp,\n", "A =\n", " [[ 4. 2. 7.]\n", " [ 1. -6. 2.]\n", " [ 3. 5. -6.]]\n" ] } ], "source": [ "k = 1\n", "p = 2\n", "for j in range(n):\n", " ( A[k,j] , A[p,j] ) = ( A[p,j] , A[k,j] )\n", "print(f\"After swapping rows 2 <-> 3 using a loop and tuples of elements, no temp,\")\n", "print(f\"A =\\n {A}\")" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "After swapping rows 1 <-> 2 using tuples of slices, no loop or temp,\n", "A =\n", " [[ 1. -6. 2.]\n", " [ 4. 2. 7.]\n", " [ 3. 5. -6.]]\n" ] } ], "source": [ "k = 0\n", "p = 1\n", "( A[k,k:] , A[p,k:] ) = ( A[p,k:].copy() , A[k,k:].copy() )\n", "print(f\"After swapping rows 1 <-> 2 using tuples of slices, no loop or temp,\")\n", "print(f\"A =\\n {A}\")" ] } ], "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.9.18" } }, "nbformat": 4, "nbformat_minor": 4 }