From 3609f67ffc42cb0ceeeab6a0a919269fa257fb59 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 3 Feb 2026 14:48:04 -0600 Subject: [PATCH] Add preconditioning demo --- demos/linear_systems/Preconditioning.ipynb | 265 +++++++++++++++++++++ 1 file changed, 265 insertions(+) create mode 100644 demos/linear_systems/Preconditioning.ipynb diff --git a/demos/linear_systems/Preconditioning.ipynb b/demos/linear_systems/Preconditioning.ipynb new file mode 100644 index 0000000..8b6e058 --- /dev/null +++ b/demos/linear_systems/Preconditioning.ipynb @@ -0,0 +1,265 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "7e4b7662-ec8a-46d1-b150-2da8e0a4a77e", + "metadata": {}, + "source": [ + "# Preconditioning\n", + "\n", + "Copyright (C) 2026 Andreas Kloeckner\n", + "\n", + "
\n", + "MIT License\n", + "Permission is hereby granted, free of charge, to any person obtaining a copy\n", + "of this software and associated documentation files (the \"Software\"), to deal\n", + "in the Software without restriction, including without limitation the rights\n", + "to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n", + "copies of the Software, and to permit persons to whom the Software is\n", + "furnished to do so, subject to the following conditions:\n", + "\n", + "The above copyright notice and this permission notice shall be included in\n", + "all copies or substantial portions of the Software.\n", + "\n", + "THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n", + "IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n", + "FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n", + "AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n", + "LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n", + "OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN\n", + "THE SOFTWARE.\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": 126, + "id": "b883645e-2eea-41f5-9e24-0f33697659c9", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import numpy.linalg as la\n", + "\n", + "np.set_printoptions(precision=3, linewidth=100, threshold=30)\n", + "n = 200" + ] + }, + { + "cell_type": "code", + "execution_count": 152, + "id": "42a0c17a-4bef-4c27-9d0c-9610f9a288fc", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[-2.000e+00, 1.000e+00, 0.000e+00, ..., 0.000e+00, 0.000e+00, 0.000e+00],\n", + " [ 8.704e-01, -1.741e+00, 8.704e-01, ..., 0.000e+00, 0.000e+00, 0.000e+00],\n", + " [ 0.000e+00, 7.575e-01, -1.515e+00, ..., 0.000e+00, 0.000e+00, 0.000e+00],\n", + " ...,\n", + " [ 0.000e+00, 0.000e+00, 0.000e+00, ..., -2.640e-12, 1.320e-12, 0.000e+00],\n", + " [ 0.000e+00, 0.000e+00, 0.000e+00, ..., 1.149e-12, -2.298e-12, 1.149e-12],\n", + " [ 0.000e+00, 0.000e+00, 0.000e+00, ..., 0.000e+00, 1.000e-12, -2.000e-12]],\n", + " shape=(200, 200))" + ] + }, + "execution_count": 152, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "scale = 10**np.linspace(0, -12, n)\n", + "D = -2*np.diag(scale) + np.diag(scale[:-1], 1) + np.diag(scale[1:], -1)\n", + "D" + ] + }, + { + "cell_type": "code", + "execution_count": 153, + "id": "1c07e587-a800-4069-86d4-fd0aabe22241", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "np.float64(263950477606337.56)" + ] + }, + "execution_count": 153, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "la.cond(D, 2)" + ] + }, + { + "cell_type": "code", + "execution_count": 154, + "id": "d00da4cd-43d6-4d29-bd93-7fb8c61c31f1", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([-1.777e+00, 9.652e-01, -1.031e+00, ..., -4.409e-12, 1.772e-12, -2.863e-13], shape=(200,))" + ] + }, + "execution_count": 154, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "xtrue = np.random.randn(n)\n", + "b = D @ xtrue\n", + "b" + ] + }, + { + "cell_type": "code", + "execution_count": 155, + "id": "7b6839eb-f17f-482a-844a-111466a54715", + "metadata": {}, + "outputs": [], + "source": [ + "xhat = la.solve(D, b)" + ] + }, + { + "cell_type": "markdown", + "id": "57ba8ac8-53fa-442a-bfe8-a465cac73861", + "metadata": {}, + "source": [ + "Now fix the row scaling in $A$:" + ] + }, + { + "cell_type": "code", + "execution_count": 156, + "id": "05ae5301-7aa7-45d0-9f08-0142c3e80342", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[-2. 1. 0. ... 0. 0. 0.]\n", + " [ 1. -2. 1. ... 0. 0. 0.]\n", + " [ 0. 1. -2. ... 0. 0. 0.]\n", + " ...\n", + " [ 0. 0. 0. ... -2. 1. 0.]\n", + " [ 0. 0. 0. ... 1. -2. 1.]\n", + " [ 0. 0. 0. ... 0. 1. -2.]]\n" + ] + } + ], + "source": [ + "#clear\n", + "M = np.diag(1/scale)\n", + "print(M@D)" + ] + }, + { + "cell_type": "markdown", + "id": "1742aeb2-16d1-49dd-be87-2ccb48c93c4b", + "metadata": {}, + "source": [ + "And solve the modified linear system:" + ] + }, + { + "cell_type": "code", + "execution_count": 157, + "id": "8bf53b00-db75-4b40-b411-90108bc715e5", + "metadata": {}, + "outputs": [], + "source": [ + "#clear\n", + "x = la.solve(M@D, M@b)" + ] + }, + { + "cell_type": "code", + "execution_count": 158, + "id": "e0d29d27-b7a0-4af5-bea4-3c2921a0116f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(np.float64(9.847440878158656e-14), np.float64(6.403057875184715e-14))" + ] + }, + "execution_count": 158, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "la.norm(xtrue-x, 2)/la.norm(x, 2), la.norm(xtrue-xhat, 2)/la.norm(x, 2)" + ] + }, + { + "cell_type": "markdown", + "id": "b15bdf84-6534-4425-9468-667e19df9f13", + "metadata": {}, + "source": [ + "What do you think the residuals will be in each case?" + ] + }, + { + "cell_type": "code", + "execution_count": 159, + "id": "08bd15c6-2104-4dbc-9c74-d90ee1c9701b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(np.float64(8.510308541902259e-17), np.float64(8.74628147040319e-17))" + ] + }, + "execution_count": 159, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "la.norm(D@x - b, 2)/la.norm(b, 2), la.norm(D@xhat - b, 2)/la.norm(b, 2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eaa16442-5bdd-4ae5-a76b-616ef9a51bc4", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.14.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}