Skip to content

Instantly share code, notes, and snippets.

@tspspi

tspspi/krylovlinear01.ipynb

Last active March 21, 2022 10:56
Show Gist options
  • Save tspspi/0d7aa7e7b119dbca55c08ba601cf9dcf to your computer and use it in GitHub Desktop.
Save tspspi/0d7aa7e7b119dbca55c08ba601cf9dcf to your computer and use it in GitHub Desktop.
Download ZIP
Steepest descent, Conjugate gradinet and Generalized minimal residual
Display the source blob
Display the rendered blob
Raw
krylovlinear01.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "2679a027",
"metadata": {},
"source": [
"# Steepest descent, Conjugate gradinet and Generalized minimal residual\n",
"\n",
"## Introduction\n",
"\n",
"This notebook contains a mini implementation of the three linear equation solving algorithms described in the [accompanying blog article](https://www.tspi.at/2022/01/02/krylovsubspacelinearsolver.html). These are:\n",
"\n",
"* Steepest descent (only for symmetric positive definit matrices; slow but guaranteed to converge - because of it's slowness it's not often used in real world applications)\n",
"* Conjugate gradient (CG; only for symmetric positive definit matrices)\n",
"* Generalized minimal residual (GMRES; works for general systems of equations, requires more memory than cojugate gradient and converges a little bit slower, requires a little bit more ressources due to generalized minimization instead of line search)\n",
"\n",
"GMRES will use the [QR decomposition using Givens rotations](https://www.tspi.at/2021/12/08/qrgivens.html) as well as the [least squares solver using QR decomposition with Givens rotations](https://www.tspi.at/2021/12/18/leastsquaresqr.html). As of today the CG method as well as GMRES and their variants are the most commonly used algorithms for computational approaches such as finite difference and finite element methods.\n",
"\n",
"In addition it contains a simple implementation of a Gaussian elimination based solver to provide a reference solution for the example problems."
]
},
{
"cell_type": "markdown",
"id": "990fd06a",
"metadata": {},
"source": [
"## Basic functions\n",
"\n",
"The following block will just implement basic matrix functions and contains the QR decomposition functions as well as least squares implementations out of the linked blog articles. In contrast to the previous ones this ```Matrix``` class allows one to change the effective size of the matrix while still keeping a constant storage allocation as maximum dimension. This will be used only for GMRES at the end of this notebook."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "298f7b88",
"metadata": {},
"outputs": [],
"source": [
"import math\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "f2786145",
"metadata": {},
"outputs": [],
"source": [
"class Matrix:\n",
" def __init__(self, rows, cols, *, storageCols = None, storageRows = None, rounding = 8, initList = None, name = None, description = None):\n",
" if not storageCols:\n",
" storageCols = cols\n",
" else:\n",
" if storageCols < cols:\n",
" raise ValueError(\"Cannot store less than used columns\")\n",
" if not storageRows:\n",
" storageRows = rows\n",
" else:\n",
" if storageRows < rows:\n",
" raise ValueError(\"Cannot store less than used rows\")\n",
"\n",
" self.name = name\n",
" self.description = description\n",
" self.rows = rows\n",
" self.cols = cols\n",
" self.storageRows = storageRows\n",
" self.storageCols = storageCols\n",
" self.data = [ 0 ] * storageRows\n",
" self.rounding = rounding\n",
" for row in range(storageRows):\n",
" self.data[row] = [ 0.0 ] * storageCols\n",
"\n",
" if initList:\n",
" if not isinstance(initList, list):\n",
" raise ValueError(\"List initializer has to be a list ...\")\n",
" # Two ways: Either one list with rows * cols entries just ordered in row major\n",
" # order (first row, second row, etc.) or nested list\n",
" if len(initList) == (rows * cols):\n",
" for row in range(rows):\n",
" for col in range(cols):\n",
" try:\n",
" k = float(initList[col + row * cols])\n",
" self.data[row][col] = round(k, self.rounding)\n",
" except ValueError:\n",
" raise ValueError(\"Initializer list has to contain numbers\")\n",
" else:\n",
" if len(initList) != rows:\n",
" raise ValueError(\"Row mismatch in initializer list\")\n",
" for row in range(rows):\n",
" if not isinstance(initList[row], list):\n",
" raise ValueError(\"Columns have to be lists in initializer\")\n",
" if len(initList[row]) != cols:\n",
" raise ValueError(\"Column mismatch in initializer list\")\n",
" for col in range(cols):\n",
" try:\n",
" k = float(initList[row][col])\n",
" self.data[row][col] = round(k, self.rounding)\n",
" except ValueError:\n",
" raise ValueError(\"Initializer list has to contain numbers\")\n",
"\n",
" def setSize(self, rows, cols):\n",
" if cols > self.storageCols:\n",
" raise ValueError(\"Cannot extend over allocated storage\")\n",
" if rows > self.storageRows:\n",
" raise ValueError(\"Cannot extend over allocated storage\")\n",
" self.rows = rows\n",
" self.cols = cols\n",
" def get(self, row, col):\n",
" if (row < 0) or (row >= self.rows):\n",
" raise IndexError(\"Row {} invalid (matrix is {} x {})\".format(row, self.rows, self.cols))\n",
" if (col < 0) or (col >= self.cols):\n",
" raise IndexError(\"Column {} invalid (matrix is {} x {})\".format(col, self.rows, self.cols))\n",
" return self.data[row][col]\n",
" def set(self, row, col, value):\n",
" if (row < 0) or (row >= self.rows):\n",
" raise IndexError(\"Row {} invalid (matrix is {} x {})\".format(row, self.rows, self.cols))\n",
" if (col < 0) or (col >= self.cols):\n",
" raise IndexError(\"Column {} invalid (matrix is {} x {})\".format(col, self.rows, self.cols))\n",
" self.data[row][col] = round(value, self.rounding)\n",
" def copy(self, *, name = None, description = None):\n",
" result = Matrix(self.rows, self.cols, rounding = self.rounding, name = name, description = description)\n",
" for row in range(self.rows):\n",
" for col in range(self.cols):\n",
" result.set(row, col, self.get(row, col))\n",
" return result\n",
" def getName(self):\n",
" return self.name\n",
" def getDescription(self):\n",
" return self.description\n",
" def __str__(self):\n",
" resString = \"\"\n",
" if self.name:\n",
" resString = resString + \"{} = \".format(self.name)\n",
" for row in range(self.rows):\n",
" for col in range(self.cols):\n",
" resString = resString + \"{}\".format(self.data[row][col]) + \"\\t\"\n",
" resString = resString + \"\\n\"\n",
" return resString\n",
" def transpose(self, *, name = None, description = None):\n",
" result = Matrix(self.cols, self.rows, rounding = self.rounding, name = name, description = description)\n",
" for col in range(self.cols):\n",
" for row in range(self.rows):\n",
" result.set(col, row, self.get(row, col))\n",
" return result\n",
" def add(self, other, *, name = None, description = None):\n",
" if not isinstance(other, Matrix):\n",
" raise ValueError(\"Adding matrix with non matrix is not supported\")\n",
" if (other.cols != self.cols) or (other.rows != self.rows):\n",
" raise ValueError(\"Addition requires same dimensions\")\n",
" result = Matrix(self.rows, self.cols, rounding = self.rounding, name = name, description = description)\n",
" for col in range(self.cols):\n",
" for row in range(self.rows):\n",
" result.set(row, col, round(self.get(row, col) + other.get(row, col), self.rounding))\n",
" return result\n",
" def sub(self, other, *, name = None, description = None):\n",
" if not isinstance(other, Matrix):\n",
" raise ValueError(\"Adding matrix with non matrix is not supported\")\n",
" if (other.cols != self.cols) or (other.rows != self.rows):\n",
" raise ValueError(\"Addition requires same dimensions\")\n",
" result = Matrix(self.rows, self.cols, rounding = self.rounding, name = name, description = description)\n",
" for col in range(self.cols):\n",
" for row in range(self.rows):\n",
" result.set(row, col, round(self.get(row, col) - other.get(row, col), self.rounding))\n",
" return result\n",
" def multiply(self, other, *, name = None, description = None):\n",
" # Calculate self * other\n",
" if isinstance(other, Matrix):\n",
" if not self.cols == other.rows:\n",
" raise ValueError(\"Colums of left side {} have to match rows of right side {}\".format(self.cols, other.rows))\n",
" result = Matrix(self.rows, other.cols, rounding = self.rounding, name = name, description = description)\n",
" for row in range(self.rows):\n",
" for col in range(other.cols):\n",
" t = 0\n",
" for k in range(self.cols):\n",
" t = t + self.get(row, k) * other.get(k, col)\n",
" t = round(t, 4)\n",
" result.set(row, col, round(t, self.rounding))\n",
" return result\n",
" else:\n",
" try:\n",
" otherFloat = float(other)\n",
" result = Matrix(self.rows, self.cols, rounding = self.rounding, name = name, description = description)\n",
" for r in range(self.rows):\n",
" for c in range(self.cols):\n",
" result.set(r,c, round(self.get(r,c,)*otherFloat, self.rounding))\n",
" return result\n",
" except ValueError:\n",
" raise ValueError(\"Matrix has to be multiplied by matrix or scalar\")\n",
" def norm(self):\n",
" if self.cols != 1:\n",
" raise NotImplementedError(\"Currently norm is not implemented for matrices, only for vectors\")\n",
" norm = 0\n",
" for i in range(self.rows):\n",
" norm = norm + self.get(i, 0)*self.get(i,0)\n",
" return math.sqrt(norm)\n",
" def isZero(self):\n",
" for row in range(self.rows):\n",
" for col in range(self.cols):\n",
" if self.get(row, col) != 0:\n",
" return False\n",
" return True\n",
" def unitVector(self):\n",
" if self.cols != 1:\n",
" raise ValueError(\"Unit vector is only supported on vectory (single column) quantities\")\n",
" norm = 0\n",
" for i in range(self.rows):\n",
" norm = norm + self.get(i, 0)*self.get(i,0)\n",
" norm = math.sqrt(norm)\n",
" if norm != 0:\n",
" return self.multiply(1.0 / norm)\n",
" else:\n",
" return self.copy()\n",
" def getColumnVector(self, col):\n",
" res = Matrix(self.rows, 1, rounding = self.rounding)\n",
" for i in range(self.rows):\n",
" res.set(i,0,self.get(i, col))\n",
" return res\n",
" def setColumnVector(self, col, v):\n",
" if v.rows != self.rows:\n",
" raise ValueError(\"Row count has to match\")\n",
" if v.cols != 1:\n",
" raise ValueError(\"Vector is not a vector ...\")\n",
" if (col >= self.cols) or (col < 0):\n",
" raise ValueError(\"Column index is invalid\")\n",
" for i in range(self.rows):\n",
" self.set(i, col, v.get(i, 0))\n",
" def toListVector(self):\n",
" if self.cols != 1:\n",
" raise ValueError(\"Vector is not a vector ...\")\n",
" res = []\n",
" for i in range(self.rows):\n",
" res.append(self.get(i, 0))\n",
" return res\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "d872dda7",
"metadata": {},
"outputs": [],
"source": [
"class IdentityMatrix(Matrix):\n",
" def __init__(self, size, *, name = None, description = None):\n",
" super().__init__(size, size, name = name, description = description)\n",
" for k in range(size):\n",
" self.set(k,k,1)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "086a6ef4",
"metadata": {},
"outputs": [],
"source": [
"class GivensRotationMatrix(IdentityMatrix):\n",
" def __init__(self, size, i, j, x1, x2, *, name = None, description = None):\n",
" super().__init__(size, name = name, description = description)\n",
" if (i >= size) or (i < 0) or (j >= size) or (j < 0):\n",
" raise ValueError(\"Invalid first index i or j\")\n",
" # Calculate s and c\n",
" s = 1.0 / math.sqrt(x1 * x1 / (x2 * x2) + 1)\n",
" c = s * x1 / x2\n",
" self.set(i,i,round(c, self.rounding))\n",
" self.set(j,i,round(-1 * s, self.rounding))\n",
" self.set(i,j,round(s, self.rounding))\n",
" self.set(j,j,round(c, self.rounding))\n",
"class GivensRotationMatrixFromCS(IdentityMatrix):\n",
" def __init__(self, size, i, j, c, s, *, name = None, description = None):\n",
" super().__init__(size, name = name, description = description)\n",
" if (i >= size) or (i < 0) or (j >= size) or (j < 0):\n",
" raise ValueError(\"Invalid first index i or j\")\n",
" self.set(i, i, round(c, self.rounding))\n",
" self.set(j, i, round(-1 * s, self.rounding))\n",
" self.set(i, j, round(s, self.rounding))\n",
" self.set(j, j, round(c, self.rounding))"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "b72ac9ba",
"metadata": {},
"outputs": [],
"source": [
"def QRDecomposeWithGivens(matIn):\n",
" matR = matIn.copy()\n",
" matQ = IdentityMatrix(matIn.rows)\n",
" for col in range(matIn.cols):\n",
" for row in range(matIn.rows-1, col, -1):\n",
" if matR.get(row,col) == 0:\n",
" continue\n",
" rotationMatrix = GivensRotationMatrix(matIn.rows, row-1, row, matR.get(row-1, col), matR.get(row,col))\n",
" matR = rotationMatrix.multiply(matR)\n",
" matQ = matQ.multiply(rotationMatrix.transpose())\n",
"\n",
" return (matQ, matR)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "950a9725",
"metadata": {},
"outputs": [],
"source": [
"def leastSquares_CheckRowEmpty(matA, row, *, epsilon = 1e-5):\n",
" for k in range(matA.cols):\n",
" v = matA.get(row, k)\n",
" if v < (-1.0 * epsilon) or v > epsilon:\n",
" return False\n",
" return True\n",
"\n",
"\n",
"def leastSquares(matA, vecB):\n",
" if not isinstance(matA, Matrix):\n",
" raise ValueError(\"Coefficient matrix has to be a matrix type\")\n",
" if not isinstance(vecB, Matrix):\n",
" raise ValueError(\"Constant matrix has to be a vector (matrix) type\")\n",
" if vecB.cols != 1:\n",
" raise ValueError(\"Constant vector has to have only one row\")\n",
" if vecB.rows != matA.rows:\n",
" raise ValueError(\"Columns in coefficient matrix have to match the number of constants in constant vector\")\n",
"\n",
" q,r = QRDecomposeWithGivens(matA)\n",
" bTransformed = q.transpose().multiply(vecB)\n",
"\n",
" # Initialize our result vector that will contain our best fit for our unknowns\n",
" xResult = [ 0.0 ] * vecB.rows\n",
"\n",
" # Check how many rows we have to drop (i.e. how many variables we're going to ignore)\n",
" populatedRows = matA.rows\n",
" while leastSquares_CheckRowEmpty(r, populatedRows - 1) and (populatedRows > 0):\n",
" populatedRows = populatedRows - 1\n",
" if populatedRows == 0:\n",
" return None\n",
"\n",
" # Now perform backsubstitution ...\n",
" for row in range(populatedRows-1, -1, -1):\n",
" t = bTransformed.get(row, 0)\n",
" for i in range(row+1, populatedRows):\n",
" t = t - xResult[i] * r.get(row, i)\n",
" t = t / r.get(row, row)\n",
" xResult[row] = t\n",
"\n",
" return xResult"
]
},
{
"cell_type": "markdown",
"id": "37b4085b",
"metadata": {},
"source": [
"## The test input\n",
"\n",
"The following sets of equations will be used as test input - they have a known solution so they can be used to check if the methods are working. They are small scale of course so there is no benefit in using iterative numeric solvers for them. Also not all solvers are applicable to all problems as mentioned later.\n",
"\n",
"### Example 1, general non symmetric non positiv definit\n",
"\n",
"$$\n",
"a - b + c -d + e = 0 \\\\ \n",
"-4a + 3b - 2c + d = 0 \\\\ \n",
"16a + 8b + 4c + 2d + e = 6.75 \\\\ \n",
"24a + 12b + 2c = 0 \\\\ \n",
"32a + 12b + 4c +d = 0\n",
"$$\n",
"\n",
"The known algebraic solutions are:\n",
"\n",
"$$\n",
"a = -0.75 \\\\ \n",
"b = 1 \\\\ \n",
"c = 3 \\\\ \n",
"d = 0 \\\\ \n",
"e = -1.25\n",
"$$\n",
"\n",
"### Example 2, general non symmetric non positiv definit\n",
"\n",
"$$\n",
"1a + 0b = 3 \\\\ \n",
"0a + 2b = 4\n",
"$$\n",
"\n",
"The solutions should be obvious:\n",
"\n",
"$$\n",
"a = 3 \\\\\n",
"b = 2\n",
"$$\n",
"\n",
"### Example 3, symmetric, positiv definit\n",
"\n",
"This system of equations has been crafted by first defining a positiv definit 4 by 4 matrix\n",
"\n",
"$$\n",
"\\begin{pmatrix} 9 & -3 & 3 & 9 \\\\ -3 & 17 & -1 & -7 \\\\ 3 & -1 & 17 & 15 \\\\ 9 & -7 & 15 & 44 \\end{pmatrix}\n",
"$$\n",
"\n",
"Then it has just been assumed that the solutions should be $1, 2, 3, 4$ and the constants have been selected to satisfy this solutions:\n",
"\n",
"$$\n",
"9a - 3b + 3c + 9d = 48 \\\\ \n",
"-3a + 17b -1c -7d = 0 \\\\ \n",
"3a -1b + 17c + 15d = 112 \\\\\n",
"9a -7b +15c + 44d = 216\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "eae82b4b",
"metadata": {},
"outputs": [],
"source": [
"ex1_A = Matrix(5, 5, rounding = 8, initList = [ [ -4, 3, -2, 1, 0 ], [ 1, -1, 1, -1, 1 ], [ 16, 8, 4, 2, 1 ], [ 24, 12, 2, 0, 0 ], [ 32, 12, 4, 1, 0 ] ])\n",
"ex1_B = Matrix(5, 1, rounding = 8, initList = [ 0, 0, 6.76, 0, 0 ])\n",
"\n",
"ex2_A = Matrix(2, 2, rounding = 8, initList = [ [ 1, 0 ], [ 0, 2 ] ])\n",
"ex2_B = Matrix(2, 1, rounding = 8, initList = [ 3, 4 ])\n",
"\n",
"ex3_A = Matrix(4, 4, rounding = 8, initList = [ [ 9, -3, 3, 9 ], [ -3, 17, -1, -7 ], [ 3, -1, 17, 15 ], [ 9, -7, 15, 44 ] ])\n",
"ex3_B = Matrix(4, 1, rounding = 8, initList = [ 48, 0, 112, 216 ])"
]
},
{
"cell_type": "markdown",
"id": "8c094f53",
"metadata": {},
"source": [
"## For reference: Direct Gaussian elimination\n",
"\n",
"To get the algebraic result first (and verify the systems have been crafted properly) lets apply Gaussian elimination (with pivoting). The pivoting works by first searching for the largest numerical value in the matrix and then swapping rows and columns so it gets assigned to the next diagonal element. Then the diagonal element is brought to $1$ by dividing the whole row (i.e. dividing the equation in the given row by a constant). In the next step the equation is subtracted from all other equations below and above to produce a column of zeros except for the $1$ on the diagonal.\n",
"\n",
"After this process has been finished the results for each of the variables are available in the constant vector."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "2dae391a",
"metadata": {},
"outputs": [],
"source": [
"def linearGaussianElimination(matA_in, vecB_in):\n",
" matA = matA_in.copy()\n",
" vecB = vecB_in.copy()\n",
" \n",
" if matA.rows != matA.cols:\n",
" raise ValueError(\"Simple Gaussian solver requires n equations for n variables (not over or underdetermined)\")\n",
"\n",
" variablePositions = [*range(matA.cols)]\n",
"\n",
" for nCol in range(matA.cols):\n",
" # Find pivot value for our diagonal element (by searching for the largest\n",
" # numerical value in the remainign submatrix)\n",
" cmax = nCol\n",
" rmax = nCol\n",
" vmax = matA.get(rmax, cmax)\n",
" for c in range(nCol, matA.cols):\n",
" for r in range(nCol, matA.rows):\n",
" if matA.get(r,c) > vmax:\n",
" vmax = matA.get(r,c)\n",
" rmax = r\n",
" cmax = c\n",
"\n",
" # Swap in the largest element by row and column swaps into our diagonal\n",
" # if required. First swap columns\n",
" if nCol != cmax:\n",
" for r in range(matA.rows):\n",
" t = matA.get(r, nCol)\n",
" matA.set(r, nCol, matA.get(r, cmax))\n",
" matA.set(r, cmax, t)\n",
" t = variablePositions[nCol]\n",
" variablePositions[nCol] = variablePositions[cmax]\n",
" variablePositions[cmax] = t\n",
"\n",
" # Swap rows if required\n",
" if nCol != rmax:\n",
" for c in range(matA.cols):\n",
" t = matA.get(nCol, c)\n",
" matA.set(nCol, c, matA.get(rmax, c))\n",
" matA.set(rmax, c, t)\n",
" t = vecB.get(nCol, 0)\n",
" vecB.set(nCol, 0, vecB.get(rmax, 0))\n",
" vecB.set(rmax, 0, t)\n",
"\n",
" # Now divide current row to get a 1 on the diagonal\n",
" for i in range(matA.cols):\n",
" if i != nCol:\n",
" matA.set(nCol, i, matA.get(nCol, i) / matA.get(nCol, nCol))\n",
" vecB.set(nCol, 0, vecB.get(nCol, 0) / matA.get(nCol, nCol))\n",
" matA.set(nCol, nCol, 1.0)\n",
"\n",
" # And eliminate elements to the top and bottom\n",
" for r in range(matA.rows):\n",
" if r != nCol:\n",
" factor = matA.get(r, nCol)\n",
" for c in range(matA.cols):\n",
" matA.set(r, c, matA.get(r,c) - factor * matA.get(nCol, c))\n",
" vecB.set(r, 0, vecB.get(r,0) - factor * vecB.get(nCol, 0))\n",
" # Now create our result vector by swapping back variables into the correct position\n",
" res = [ 0.0 ] * matA.cols\n",
" for i in range(matA.cols):\n",
" res[variablePositions[i]] = vecB.get(i, 0)\n",
" return res"
]
},
{
"cell_type": "markdown",
"id": "3d99e2a8",
"metadata": {},
"source": [
"Now lets apply this to our example system of equations:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "fda9a3a8",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[-0.75, 1.0, 3.0, -0.0, -1.25]\n"
]
}
],
"source": [
"eliminationResults = linearGaussianElimination(ex1_A, ex1_B)\n",
"for i in range(len(eliminationResults)):\n",
" eliminationResults[i] = round(eliminationResults[i], 2)\n",
"print(eliminationResults)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "e06fecf2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[3.0, 2.0]\n"
]
}
],
"source": [
"eliminationResults2 = linearGaussianElimination(ex2_A, ex2_B)\n",
"for i in range(len(eliminationResults2)):\n",
" eliminationResults2[i] = round(eliminationResults2[i], 2)\n",
"print(eliminationResults2)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "e1195785",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1.0, 2.0, 3.0, 4.0]\n"
]
}
],
"source": [
"eliminationResults3 = linearGaussianElimination(ex3_A, ex3_B)\n",
"for i in range(len(eliminationResults3)):\n",
" eliminationResults3[i] = round(eliminationResults3[i], 2)\n",
"print(eliminationResults3)"
]
},
{
"cell_type": "markdown",
"id": "5a6652f5",
"metadata": {},
"source": [
"As one can see this matches the results that have also been derived by hand calculation above"
]
},
{
"cell_type": "markdown",
"id": "2efafb1c",
"metadata": {},
"source": [
"## Steepest descent\n",
"\n",
"The first algorithm that will be implemented ist steepest descent. This is usually also the slowest one - but it's the only one that is guaranteed to converge. Recall that steepest descent only works for __positive definit symmetric__ matrices.\n",
"\n",
"Steepest descent is based on the fact that solving $A\\vec{x} = \\vec{b}$ is equivalent to solving the minimization problem $\\phi(x) = \\vec{x}^T A \\vec{x} - \\vec{x}^T \\vec{b} \\to \\text{min}$. As shown in the [blog post](https://www.tspi.at/2022/01/02/krylovsubspacelinearsolver.html) the gradient $-\\nabla \\phi(\\vec{x}_k)$ is equivalent to the residuum $\\vec{r}_k$. After having discovered the direction of steepest descent of the residuum the second ingredient in steepest descent is the exact line search along that gradient - the solution of the exact line search is a second order polynomial solution so it's possible to be calculated exact. It turns out the distance that one has to step is determined by the spectral coefficients:\n",
"\n",
"$$\n",
"\\alpha_k = \\frac{\\vec{r}_k^T \\vec{r}_k}{\\vec{r}_k^T A \\vec{r}_k}\n",
"$$\n",
"\n",
"To take a step from a previously guessed or calculated estimate $\\vec{x}_k$ one just has to move into the direction of the gradient $\\vec{r}_k$ scaled by $\\alpha_k$ which yields the iteration rule\n",
"\n",
"$$\n",
"\\vec{x}_{k+1} = \\vec{x}_k + \\alpha_k \\vec{r}_k\n",
"$$\n",
"\n",
"The residuum can be calculated after every step:\n",
"\n",
"$$\n",
"\\vec{r}_k = \\vec{b} - A \\vec{x}_k\n",
"$$\n",
"\n",
"Then one repeats the calculation of the residuum and the movement along the gradient until the residuum is either $\\vec{0}$ or at least small enough (note that you should not compare floating point numbers with zero but check if their absolute value is within some kind of margin - usually called $\\epsilon$ - around zero)."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "a3216c2c",
"metadata": {},
"outputs": [],
"source": [
"def steepestDescent_Terminate(r):\n",
" stepSize = r.transpose().multiply(r).get(0,0)\n",
" if stepSize < 1e-4:\n",
" return True\n",
" else:\n",
" return False\n",
"\n",
"def steepestDescent(matA_in, vecB_in, rounding = 3, initialGuess = None, trace = False, convergencePlot = False):\n",
" # If we don't have an initial guess (which we should have) generate one ...\n",
" if not initialGuess:\n",
" # initialGuess = Matrix(vecB_in.rows, 1, rounding = vecB_in.rounding, initList = [ 0.0 ] * vecB_in.rows)\n",
" initialGuess = Matrix(vecB_in.rows, 1, rounding = vecB_in.rounding)\n",
" for i in range(vecB_in.rows):\n",
" initialGuess.set(i, 0, matA_in.get(i,i,))\n",
" else:\n",
" initialGuess = initialGuess.copy()\n",
"\n",
" if convergencePlot:\n",
" convergencePlotData = { 'iteration' : [], 'r' : [] }\n",
"\n",
" # Calculate initial residuum\n",
" xProjected = matA_in.multiply(initialGuess)\n",
" r = vecB_in.sub(xProjected)\n",
" x = initialGuess\n",
"\n",
" # Perform our iterations\n",
" iterationCount = 0\n",
" while not steepestDescent_Terminate(r):\n",
" if convergencePlot:\n",
" convergencePlotData['iteration'].append(len(convergencePlotData['iteration']))\n",
" convergencePlotData['r'].append(r.transpose().multiply(r).get(0,0))\n",
" \n",
" # Calculate spectral coefficient (distance to jump)\n",
" spectralCoefficient = r.transpose().multiply(r).get(0,0) / r.transpose().multiply(matA_in.multiply(r)).get(0,0)\n",
"\n",
" # Perform step into residuum direction\n",
" x = x.add(r.multiply(spectralCoefficient))\n",
"\n",
" # Update our residuum\n",
" xProjected = matA_in.multiply(x)\n",
" r = vecB_in.sub(xProjected)\n",
" iterationCount = iterationCount + 1\n",
"\n",
" # Round our result as specified\n",
" res = [ 0.0 ] * x.rows\n",
" for i in range(x.rows):\n",
" res[i] = round(x.get(i, 0), rounding)\n",
"\n",
" if convergencePlot:\n",
" plt.plot(convergencePlotData['iteration'], convergencePlotData['r'])\n",
" plt.title(\"Convergence plot (Residuum)\")\n",
" plt.grid()\n",
" plt.xlabel(\"Iteration\")\n",
" plt.ylabel(\"Residuum b - Ax\")\n",
" return (res, r)"
]
},
{
"cell_type": "markdown",
"id": "8af0f7ec",
"metadata": {},
"source": [
"Applying this method to our symmetric positiv definit problem:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "0c51393e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1.001, 2.0, 3.0, 4.0]\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAi60lEQVR4nO3deZxkZX3v8c+3V6a7hlm6hxEQZgJoiJqoNO4EGWMUCC4hcY+aSEDuDS7XFS++jFdNbhKjF7fEFTGIjEtEyQii0UaCgsDgMIJoHDZFBmYVpmft5Xf/OKdmaoqu7pqeOrWc832/XvXqU+c8dc6vTlf/6unnPM9zFBGYmVn+dLU6ADMzy4YTvJlZTjnBm5nllBO8mVlOOcGbmeWUE7yZWU45wZsdAEl/Kem6Bu/zh5Ke3Mh9TnOMMUnH1NjW8PdUZ0z9kn4u6bBmH7sonOBzRNIrJd2c/jGvl3SVpJNaHVdRSbpG0l/PUuYFwLaI+En6/L2SxtPf4W8l/UjSMw42logoRcRdB7ufRoqI3cBFwDtbHUteOcHnhKS3ABcCfw8sBY4G/gV4UQvD2o+knlbH0IbOBS6pWvfliCgBw8Ao8NWmR9U8XwJeK6m/1YHkkRN8DkhaALwP+JuI+HpEbI+I8Yj4j4h4e1qmX9KFku5PHxeW/6gknSLpPklvlbQhrf3/Vbrt6ZIekNRdcbw/lbQ2Xe6SdL6kOyVtlvQVSYvTbcslhaSzJP0K+L6kbkkfkrRJ0t2SzkvL9JTfi6TPpTH8RtIHyscuNyVI+mdJW9PXn1YR12JJn0/f31ZJ36jYdoakNRW14j+Y4XyGpDdKuiuN84OSpv1bkfRMSTdJeij9+cx0/d8Bfwh8PK2Nf3ya1/YBzwF+MN2+I2ICuBQ4UtKSOs7PcZJ+kMaySdKXq97TcenykKQrJD0s6Ubg2Ipy5d9ZT8W6vf+JpP9hfLFW+bTsB9JzPCbpP9LjXZoe7yZJyyve433AVuDptX4fNndO8PnwDOAQ4PIZylxA8kf0JOCJwFOBd1dsfxSwADgSOAv4hKRFEXEDsJ0kEZW9kqTmBfBG4MXAs4EjSP5YP1F17GcDvwc8HzgbOC2N44T0tZW+AEwAxwFPBp4HVDZzPA34BUnt9p+Az0lSuu0SYAB4PHAY8P8AJJ1A0hTwemAI+BRwxSy1xj8FTkxjfBHwuuoC6RfZt4CPpvv9MPAtSUMRcQHwX8B5afPIedMc4zHAVJrkHiH9AngNsJnkvM52ft4PfAdYBDwa+FiN9/YJYBdwePq+HvHeDtLLgVeTfJaOBa4HPg8sBu4A/raq/B0kn0lrtIhoqwfJH+IG4LY6y78U+BlwO/ClVsffonP2KuCBWcrcCZxe8fz5wD3p8inATqCnYvsG4Onp8geAi9Ll+SQJf1n6/A7gjypedzgwDvQAy4EAjqnY/n3g9RXPn5uW6SFpWtoNzKvY/gpgNF3+S2BdxbaB9LWPSo87BSya5r3/K/D+qnW/AJ5d41wFcGrF8/8JfK8ihuvS5VcDN1a99nrgL9Pla4C/nuF38qzq3xvwXmAP8FtgkiS5n5Jum+38/BvwaeDRNd7TcUB3+vs5vmLb31e8p/LvrPKzsPd9pPF9sWLbfuXTshdUbP8QcFXF8xcAa6piuxR4T6v/jvL4aMca/MXAqfUUlPQY4F3AsyLi8cCbswurrW0Ghmdp4z4CuLfi+b3pur37iKRJoGwHUEqXvwScmdZ4zwRuiYjyvpYBl6dNH78lSfiTJMmo7NdVcfy6xrZlQC+wvmJ/nyKpjZc9UF6IiB3pYgk4CtgSEVt5pGXAW8v7TPd7VNX7r1YZV/W5qnwv91atu5ek5lqPrSRfmNW+EhELSc7hbcBIun628/MOQMCNkm6XNF3NfAnJl2n1+2ukByuWd07zvLR/ceaTfKFZg7Vdgo+Ia4EtleskHSvp25JWS/ovScenm84GPlH+o46IDU0Ot11cT/Iv94tnKHM/SYIoOzpdN6uI+BlJEjiN/ZtnIEkUp0XEworHIRHxm8pdVCyvJ2k+KDuqal+7geGKfR2afnnP5tfAYkkLa2z7u6oYByLishn2VxlXrXNVfU7LZcvvfbapWn8JSNK0XwgRsYmkWem9kg5nlvMTEQ9ExNkRcUT6un8pt7tX2EjSxFP9/sq2pz8HKtY9qmp7rW1z9XvArQ3Yj1VpuwRfw6eBN0TECPA2kt4hAI8FHqukH/ENkuqq+edNRDwEvIek3fzFkgYk9Uo6TdI/pcUuA94taYmk4bT8F2vtcxpfImlvP5n9e3V8Evg7ScsA0v3P1HPnK8CbJB2ZJuO9XeQiYj1JG/KHJB2q5ALusZKePVtw6WuvIklqi9L3f3K6+TPAuZKepsSgpD+RNF3tuezt6X6OAt4EfHmaMleSfP5eKalH0suAxwGr0u0PAtP2PU9jHgf+k+QaRa0yPweuBt4x2/mR9BJJ5S/PrSRfMJNV+5sEvk7ypTEg6XHAayu2byT5gvoLJRfEX0fFRVhgDXCypKOVXNx/V63Y65F+uS0GbjiY/dj02j7BSyoBzwS+KmkNyb+kh6ebe0guVJ1C0hb52Ro1uNyLiA8DbyG5cLqRpLZ3HvCNtMgHgJuBtcBPgVvSdfW6jOQ8fz+tWZZ9BLgC+I6kbSR/qE+bYT+fIUlSa4GfkCTJCfYlotcAfSTXVbYCX2Pf73s2ryZpX/45yTWENwNExM0k/+19PN3nOpK29Jl8E1hNktC+BXyuukBEbAbOAN5K0kz2DuCMivPzEeDPlfTo+WiN43wqjXsmHwTOUTIgaKbz8xTgx5LGSH4nb4qIu6fZ33kkzSQPkDSJfr5q+9nA29P39HjgRxXv+bskX3ZrSc7PKg7OK4EvRNIn3hpMEe13w4+0G9WqiHiCpEOBX0TEI/7IJX0SuCEiLk6ffw84PyJuama8NndKujl+MiKqmzpaRlIAj4mIdU063nUk/6H+pBnHaxfpNZ1bgZML3LyaqbavwUfEw8Ddkl4CSYOlpHKXqm8AK9L1wyRNNm01Ws/2J2mepNPTJo0jSbrMzdS9M/ci4qSiJXdIRrJGxPFO7tlpuwQv6TKSi4a/q2TwzVkk3QDPknQrSXfIchvv1cBmST8jGfH39vTfZmtfAv4PSfPCT0h63bynpRGZ5VRbNtGYmdnBa7savJmZNUZbTf40PDwcy5cvn9Nrt2/fzuDgYGMDyoDjbLxOidVxNlanxAnZxrp69epNEbFk2o2tHkpb+RgZGYm5Gh0dnfNrm8lxNl6nxOo4G6tT4ozINlbg5uigqQrMzKwBnODNzHLKCd7MLKec4M3McsoJ3swsp5zgzcxyygnezCynOj7BRwQf/d4v+enGidkLm5kVSMcneEl85tq7WLtpcvbCZmYF0vEJHmCo1MdDuz1pmplZpVwk+OFSP9v2OMGbmVXKRYIfKvXxsBO8mdl+cpLg+53gzcyq5CLBDw/2MbYHJqec5M3MynKR4IdK/QSwZfueVodiZtY2cpHgh0v9AGzevrvFkZiZtY9cJPihUh8Am8dcgzczK8tFgh9OE/ymMdfgzczKcpHghwbTJhrX4M3M9spFgl8wr5duuQZvZlYpFwm+q0vM75Nr8GZmFXKR4IEkwbsXjZnZXrlJ8If2wSbX4M3M9spRgncN3sysUn4SfL/YtM01eDOzsvwk+D6xc3ySHXt8ZyczM8hZggf3hTczK8tNgp+fJnj3hTczS+QmwR/a7xq8mVml3CT4Ba7Bm5ntJzcJvtxEs9lzwpuZATlK8H3dotTf4xq8mVkqNwkeknnh3QZvZpbIPMFL6pb0E0mrsj7W0GCfa/BmZqlm1ODfBNzRhOMwXOp3Dd7MLJVpgpf0aOBPgM9meZyyoVK/56MxM0spIrLbufQ14P8C84G3RcQZ05Q5BzgHYOnSpSMrV66c07HGxsa4en0fq+4c53PPH6BLOojIszM2NkapVGp1GLPqlDihc2J1nI3VKXFCtrGuWLFidUScOO3GiMjkAZwB/Eu6fAqwarbXjIyMxFyNjo7G56+7K5a9c1Vs2rZrzvvJ2ujoaKtDqEunxBnRObE6zsbqlDgjso0VuDlq5NQsm2ieBbxQ0j3ASuA5kr6Y4fEYnp/cm9XzwpuZZdgGHxHviohHR8Ry4OXA9yPiL7I6HlTefNvt8GZmueoHP1zqA2CTR7OamdHTjINExDXANVkfZ6jkGryZWVmuavAL5/XS3SX3hTczI2cJvqtLLPZoVjMzIGcJHsrTFbgGb2aWuwQ/7NGsZmZADhO8Z5Q0M0vkL8EP9rsXjZkZOUzww/P72L5nkp17JlsdiplZS+UvwQ+WpytwLd7Mii13CX4oHc3qe7OaWdHlMMF7NKuZGeQxwQ+mNXj3pDGzgstdgh9Oa/AbXYM3s4LLXYKf19fNYF+3a/BmVni5S/Dge7OamUFuE7xHs5qZ5TPBD/a7H7yZFV4uE/yS+Z5R0swslwl+aLCfLdt3MzUVrQ7FzKxl8pngS31MBfx253irQzEza5mcJniPZjUzy2WCH07no3E7vJkVWU4TvGeUNDPLZYLfNx+NE7yZFVcuE/zCgT665CmDzazYcpngu7vE4kH3hTezYstlgoekHd5t8GZWZLlN8Ml8NE7wZlZcsyZ4SZdIWlDxfJmk72Ub1sEbGux3G7yZFVo9NfjrgB9LOl3S2cB3gQszjaoBPKOkmRVdz2wFIuJTkm4HRoFNwJMj4oHMIztIw6V+xnZPsGt8kkN6u1sdjplZ09XTRPNq4CLgNcDFwJWSnphxXAdt32hWt8ObWTHNWoMH/gw4KSI2AJdJupwk0T85y8AO1tBgeT6aPTx60UCLozEza756mmheXPX8RkmvzyyiBhlKa/C+dZ+ZFVU9NXgAJD0OeDnwCuAh4MRZyh8CXAv0p8f5WkT87dxDPTD75qPxhVYzK6YZE7ykZSQJ/RXABLAMODEi7qlj37uB50TEmKRe4DpJV0XEDQcZc1321uCd4M2soGpeZJX0I+BKoBf484gYAbbVmdyJxFj6tDd9NO0WSwN9PQz0dfsiq5kV1ky9aDYC84GlwJJ03QElaEndktYAG4DvRsSP5xLkXHk0q5kVmSJq5+x0BOufkTTRHAcsBJ4fETce0EGkhcDlwBsi4raqbecA5wAsXbp0ZOXKlQey673GxsYolUr7rXvf9TuZ1wNvf8q8Oe0zC9PF2Y46JU7onFgdZ2N1SpyQbawrVqxYHRHTXxONiLoewGHAG4AfAb+u93UVr/9b4G0zlRkZGYm5Gh0dfcS6sy6+MU698No57zML08XZjjolzojOidVxNlanxBmRbazAzVEjp9Y92VhEbIiIj0XEM4GTZisvaUlac0fSPOC5wM/rPV4jDJf63URjZoVVdzfJShFxbx3FDge+IKmbpK3/KxGxai7Hm6uhUh+bt+9hairo6lIzD21m1nJzSvD1iIi1tHi069BgP5NTwUM7x1mU3sbPzKwocjsfPHg0q5kV2wEleEm3ZBVIFjya1cyK7EBr8B3VkF1O8B7NamZFdKAJ/luZRJGRIU8ZbGYFdkAJPiLenVUgWVg00IeEu0qaWSHl+iJrd5dYPNDHJt+b1cwKKNcJHjwfjZkVV1394CX1AceTTDb2i4jomCpxMpq1Y8I1M2uYeu7J+ifAncBHgY8D6ySdlnVgjTJU6vdFVjMrpHpq8B8CVkTEOgBJx5L0prkqy8AaZWiwzzV4MyuketrgN5STe+oukvndO8JwqY9tuyfYNT7Z6lDMzJqqZg1e0pnp4u2SrgS+QtIG/xLgpibE1hBD6WCnLdv3cMTC9pkX3swsazM10bygYvlB4Nnp8kZgUWYRNVjlaFYneDMrkpoJPiL+qpmBZMWjWc2sqHLfD354sDzhmBO8mRVL7hP8vimD3ZPGzIol9wl+sL+Heb3dHs1qZoUzaz/49L6qrwGWV5aPiDdmFlWDJdMVuAZvZsVSz0CnK4EbgJ8CU9mGk42hUj8bXYM3s4KpJ8EfEhFvyTySDA0P9rH+oV2tDsPMrKnqaYO/RNLZkg6XtLj8yDyyBhoq9fm+rGZWOPXU4PcAHwQuIBnJSvrzmKyCarTyjJIRgdRRdx00M5uzehL8W4DjImJT1sFkZajUz8RU8NDOcRYO9LU6HDOzpqinieZ2YEfWgWRpeO9oVvekMbPiqKcGPwmskTQK7G3I7qhukoPl+Wh2c9xhpRZHY2bWHPUk+G+kj47l0axmVkSzJviI+EIzAsnSvhkl3ZPGzIqjnpGsd7Ov98xeEdExvWgWDfQiwUa3wZtZgdTTRHNixfIhJDf86Kh+8D3dXSwa6HMN3swKZdZeNBGxueLxm4i4EHhO9qE1lu/NamZFU08TzQkVT7tIavTzM4soIx7NamZFU08TzYcqlieAe4CXZhJNhoZL/fzs/odbHYaZWdPU04tmRTMCydqwZ5Q0s4Kpp4nmPdOtj4j3NT6c7AwN9rFt1wS7Jybp7+ludThmZpmrZ6qC7RWPSeA0kpt/zEjSUZJGJd0h6XZJbzqoSA/SUNoXfosHO5lZQdTTRFPZBo+kfwauqGPfE8BbI+IWSfOB1ZK+GxE/m1uoB6c8H83msT0cvmBeK0IwM2uqudyTdYA6pgqOiPURcUu6vA24AzhyDsdriHINfpPb4c2sIBTxiEGq+xeQfsq+kazdwBLgfRHx8boPIi0HrgWeEBEPV207BzgHYOnSpSMrV66sO/hKY2NjlEq1JxLbsGOKd1y7k7/+/T5OOrJ3TsdohNnibBedEid0TqyOs7E6JU7INtYVK1asjogTp90YETM+gGUVjyOBntleU/X6ErAaOHO2siMjIzFXo6OjM27ftms8lr1zVXzymnVzPkYjzBZnu+iUOCM6J1bH2VidEmdEtrECN0eNnFqzDb7itnzbqjYdKomI2DLbN4ukXuDfgUsj4uuzlc/SYF83/T1dnlHSzApjpousq0maZgQcDWxNlxcCvwJ+Z6YdK7k33ueAOyLiw40I9mBIYrjU7zZ4MyuMmhdZI+J3Ipkx8mrgBRExHBFDwBlAPbXxZwGvBp4jaU36OL0hUc/RcMnz0ZhZcdQzVcFTIuLc8pOIuErS+2d7UURcR1LjbxtDpX4efHhXq8MwM2uKerpJbpL0bknLJS2TdAGwOevAsuAZJc2sSOpJ8K8g6Rp5Ocmt+w5L13WcoVI/m7fvLvfuMTPLtXpGsm4BWjrNQKMMl/oYnwwe3jXBgnmt6wtvZtYMM3WTvDAi3izpP5j+ln0vzDSyDFTem9UJ3szybqYa/CXpz39uRiDNMJTOR7NpbA/HLGlxMGZmGauZ4CNidfrzB+V1khYBR0XE2ibE1nBDg/tq8GZmeTfrRVZJ10g6NB3ZeivweUktH7g0F+UZJTd5NKuZFUA9vWgWRDJB2JnA5yNiBHhutmFlY/Fgecpg1+DNLP/qSfA9kg4nuQ/rqozjyVRPdxeLBnrdF97MCqGeBP8+kukK7oyImyQdA/wy27CyM+T5aMysIOrpB/9V4KsVz+8C/izLoLLk0axmVhT1XGR9rKTvSbotff4Hkt6dfWjZGC71s2m7a/Bmln/1NNF8BngXMA6QdpF8eZZBZckzSppZUdST4Aci4saqdRNZBNMMQ6V+Hto5zp6JqVaHYmaWqXpnkzyWdLoCSX8OrM80qgyVR7NucV94M8u5euaD/xvg08Dxkn4D3A28KtOoMlQezbppbDePWnBIi6MxM8tOPb1o7gKeK2mQpMa/E3gZcG/GsWViyfx0sJNr8GaWczWbaNLpCd4l6eOS/hjYAbwWWEcy6KkjeT4aMyuK2WaT3ApcD5wNvAPoA14cEWuyDy0b+2aUdII3s3ybKcEfExG/DyDps8Am4OiI2NaUyDJS6u+hr6fLXSXNLPdm6kUzXl6IiEng7k5P7gCSGB7sY5MTvJnl3Ew1+CdKejhdFjAvfS4gIuLQzKPLyPD85N6sZmZ5NtMNP7qbGUgzDbkGb2YFUM9Ap9zxjJJmVgQFTfDJfDQRj7iXuJlZbhQywQ8P9rNncoptuzt2Sh0zs1kVM8GXR7O6Hd7McqyQCd6jWc2sCIqZ4D2a1cwKoJAJfrhUnlHSTTRmll+FTPCLB90Gb2b5V8gE39vdxcKBXo9mNbNcK2SCh2Q0q2vwZpZnmSV4SRdJ2iDptqyOcTCGSv1s9EVWM8uxLGvwFwOnZrj/gzJc6nM3STPLtcwSfERcC2zJav8Ha7jU79v2mVmuKcv5WCQtB1ZFxBNmKHMOcA7A0qVLR1auXDmnY42NjVEqleou/811e7h83Tiffd4APV2a0zHn4kDjbJVOiRM6J1bH2VidEidkG+uKFStWR8SJ026MiMwewHLgtnrLj4yMxFyNjo4eUPlLrr8nlr1zVTz40M45H3MuDjTOVumUOCM6J1bH2VidEmdEtrECN0eNnFrYXjTD6WhWX2g1s7wqbIIfKpXno3E7vJnlU5bdJC8Drgd+V9J9ks7K6lhzUZ6uwIOdzCyvZron60GJiFdkte9GKE845hq8meVVYZto5vf30Nfd5QnHzCy3CpvgJTFU6vOUwWaWW4VN8FC+N6sTvJnlU6ETvEezmlmeFTrBDw32+yKrmeVWoRP8cNoGHxlO12Bm1iqFTvBDpT52T0wxtnui1aGYmTVcsRP8oEezmll+FTrBD8/3aFYzy69CJ/ih9ObbHuxkZnlU6ARfno/Gg53MLI8KneAXD3o+GjPLr0In+L6eLhbM6/VoVjPLpUIneEi6Sm7yaFYzy6HCJ/jhwX7X4M0slwqf4JMZJV2DN7P8cYL3jJJmllOFT/DDpX627hhnYnKq1aGYmTVU4RN8+ebbW3a4mcbM8qXwCX7YfeHNLKcKn+CHPJrVzHLKCb7kGryZ5VPhE7znozGzvCp8gj/0kB56u+V7s5pZ7hQ+wUtK783qGryZ5UvhEzx4NKuZ5ZMTPElPGtfgzSxvnOCBYdfgzSyHnOBJetJs3r6biGh1KGZmDeMET3Jv1l3jU+zYM9nqUMzMGsYJHo9mNbN8coJn32hWt8ObWZ44wQNL0hq8e9KYWZ70tDqAdlCuwV9x6/3smZxi2eJBjh4aYMG83hZHZmY2d5kmeEmnAh8BuoHPRsQ/ZHm8uVpS6uf4R81n1dr1rFq7fu/6hQO9LFs8wNFDg+nPAZYtHmD58CCHze9HUgujNjObWWYJXlI38Angj4H7gJskXRERP8vqmHPV093Ft998MmO7J/jV5h38ast27t28g3u37OBXm3ew5tdb+dba+5mq6EV5SG8XRy8e4OjFgywbGmDZ0ABHLR5gsC+Z26a3u4u+ni56uvYt93Z3sXMi2DU+SW93F91d/oIws+xkWYN/KrAuIu4CkLQSeBHQdgm+rNTfw+OOOJTHHXHoI7aNT07xm60706S/7wvg3s3buW7dRnaNH8At//7z2wB0iST5d3fRW/Fl0Ih/DMr7EKp6Xt6u/Z5XPhGwY8cOBm/5wcEH0gCznY7tbRTrTBxnY3VKnDB7rIsG+vjKuc9o+HGzTPBHAr+ueH4f8LTqQpLOAc4BWLp0Kddcc82cDjY2Njbn1x6oo4CjSnBSCTgaIg7ht7uDTTuD8SmYmAompmAyYCJ9Xl7evms33b19TE7BRJD8nAomYip5zdTB98WPvT9jvxXVw7gqx3VVbzt03hQ9XTsPOpaDVc/Ys1L/FN1qfayzcZyN1SlxwuyxHjKxK5v8FRGZPICXkLS7l5+/GvjYTK8ZGRmJuRodHZ3za5vJcTZep8TqOBurU+KMyDZW4OaokVOz7CZ5H0llt+zRwP0ZHs/MzCpkmeBvAh4j6Xck9QEvB67I8HhmZlYhszb4iJiQdB5wNUk3yYsi4vasjmdmZvvLtB98RFwJXJnlMczMbHqeqsDMLKec4M3McsoJ3swsp5zgzcxyStFGt6mTtBG4d44vHwY2NTCcrDjOxuuUWB1nY3VKnJBtrMsiYsl0G9oqwR8MSTdHxImtjmM2jrPxOiVWx9lYnRIntC5WN9GYmeWUE7yZWU7lKcF/utUB1MlxNl6nxOo4G6tT4oQWxZqbNngzM9tfnmrwZmZWwQnezCynOirBSzpV0i8krZN0/jTbJemj6fa1kk5oUZxHSRqVdIek2yW9aZoyp0h6SNKa9PGeFsV6j6SfpjHcPM32lp9TSb9bcZ7WSHpY0puryrTsfEq6SNIGSbdVrFss6buSfpn+XFTjtTN+ppsQ5wcl/Tz93V4uaWGN1874OWlCnO+V9JuK3+/pNV7btPM5Q6xfrojzHklrarw2+3Na604g7fYgmXL4TuAYoA+4FXhcVZnTgatIbuP5dODHLYr1cOCEdHk+8N/TxHoKsKoNzus9wPAM29vinFZ9Dh4gGdzRFucTOBk4AbitYt0/Aeeny+cD/1jjvcz4mW5CnM8DetLlf5wuzno+J02I873A2+r4bDTtfNaKtWr7h4D3tOqcdlINfu9NvCNiD1C+iXelFwH/FokbgIWSDm92oBGxPiJuSZe3AXeQ3KO2E7XFOa3wR8CdETHXEc8NFxHXAluqVr8I+EK6/AXgxdO8tJ7PdKZxRsR3ImIifXoDyZ3XWqrG+axHU88nzByrkjvbvxS4LMsYZtJJCX66m3hXJ816yjSVpOXAk4EfT7P5GZJulXSVpMc3N7K9AviOpNXpDdCrtds5fTm1/2Da4XyWLY2I9ZB84QOHTVOm3c7t60j+W5vObJ+TZjgvbUq6qEaTV7udzz8EHoyIX9bYnvk57aQEr2nWVffxrKdM00gqAf8OvDkiHq7afAtJM8MTgY8B32hyeGXPiogTgNOAv5F0ctX2tjmnSm79+ELgq9NsbpfzeSDa6dxeAEwAl9YoMtvnJGv/ChwLPAlYT9L0Ua1tzmfqFcxce8/8nHZSgq/nJt5tc6NvSb0kyf3SiPh69faIeDgixtLlK4FeScNNDpOIuD/9uQG4nOTf3Eptc05J/hBuiYgHqze0y/ms8GC5KSv9uWGaMm1xbiW9FjgDeFWkjcPV6vicZCoiHoyIyYiYAj5T4/htcT4BJPUAZwJfrlWmGee0kxJ8PTfxvgJ4Tdrz4+nAQ+V/k5spbXv7HHBHRHy4RplHpeWQ9FSS38Xm5kUJkgYlzS8vk1xwu62qWFuc01TNGlE7nM8qVwCvTZdfC3xzmjItvzG9pFOBdwIvjIgdNcrU8znJVNV1nz+tcfyWn88KzwV+HhH3Tbexaec0yyu4jX6Q9Oj4b5Ir5Rek684Fzk2XBXwi3f5T4MQWxXkSyb+Ga4E16eP0qljPA24nudJ/A/DMFsR5THr8W9NY2vmcDpAk7AUV69rifJJ86awHxklqkWcBQ8D3gF+mPxenZY8ArpzpM93kONeRtFuXP6efrI6z1uekyXFekn7+1pIk7cNbfT5rxZquv7j82awo2/Rz6qkKzMxyqpOaaMzM7AA4wZuZ5ZQTvJlZTjnBm5nllBO8mVlOOcFbLkkaS38ul/TKBu/7f1c9/1Ej92/WKE7wlnfLgQNK8JK6ZymyX4KPiGceYExmTeEEb3n3D8AfpnNu/y9J3ekc6DelE1e9HvbOJz8q6UskA2qQ9I10Iqjby5NBSfoHYF66v0vTdeX/FpTu+7Z0nu+XVez7GklfUzL3+qXlUbdmWeppdQBmGTufZB7xMwDSRP1QRDxFUj/wQ0nfScs+FXhCRNydPn9dRGyRNA+4SdK/R8T5ks6LiCdNc6wzSSbDeiIwnL7m2nTbk4HHk8yN8kPgWcB1jX6zZpVcg7eieR7J3DprSKZwHgIek267sSK5A7xRUnnqg6MqytVyEnBZJJNiPQj8AHhKxb7vi2SyrDUkTUdmmXIN3opGwBsi4ur9VkqnANurnj8XeEZE7JB0DXBIHfuuZXfF8iT+27MmcA3e8m4byW0Ty64G/kc6nTOSHpvO5ldtAbA1Te7Hk9yusGy8/Poq1wIvS9v5l5Dczu3GhrwLszlwLcLybi0wkTa1XAx8hKR55Jb0QudGpr+d3reBcyWtBX5B0kxT9mlgraRbIuJVFesvB55BMkNgAO+IiAfSLwizpvNskmZmOeUmGjOznHKCNzPLKSd4M7OccoI3M8spJ3gzs5xygjczyykneDOznPr/z6rB4nyDzncAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"x,r = steepestDescent(ex3_A, ex3_B, convergencePlot = True)\n",
"print(x)"
]
},
{
"cell_type": "markdown",
"id": "1bb5869a",
"metadata": {},
"source": [
"As one can see this matches - up to rounding error - our result obtained by Gaussian elimination:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "d0e971b3",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1.0, 2.0, 3.0, 4.0]\n"
]
}
],
"source": [
"print(eliminationResults3)"
]
},
{
"cell_type": "markdown",
"id": "93d3f129",
"metadata": {},
"source": [
"Applying the same function to our other two toy problems won't work since they are _not_ symmetric positive definit"
]
},
{
"cell_type": "markdown",
"id": "2f33345f",
"metadata": {},
"source": [
"## Conjugate gradient (CG)\n",
"\n",
"The conjugate gradient algorithm is as of today the most used algorithm when implementing finite element calculations - it works on all symmetric positive definit coefficient matrices. Basically it's an extension to the steepest descent method - but in contrast to direct steepest descent it limits the line search for the minimum along the gradient to a Krylov subspace of the problem - into which it project the residuum vector at every step. This will be the vector $\\vec{p}$ during the iteration. Since the projective subspaces are choosen to be orthogonal to each other one calls this _conjugate_ gradient method. This speeds up the calculation significantly and might terminate way before expanding to the full space of the coefficient matrix.\n",
"\n",
"Since one projects the residuum into the direction of the approximate Eigenvectors of the system - those are formed by the Arnoldi iteration usually since a power iteration approaches the Eigenvectors and orthogonalization procedure using a modified Gram-Schmidt algorithm then tries to project onto each Eigenvector descending from the one with the larges Eigenvalue - thus the most important one - downwards. Thus the algorithm first minimizes into the direction that affects the residuum the most, then into the second most important direction and so on. Since one is usually only interested into approximate solutions one can often terminate after a few iterations."
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "1c33b73e",
"metadata": {},
"outputs": [],
"source": [
"def conjugateGradient_Terminate(r):\n",
" stepSize = r.transpose().multiply(r).get(0,0)\n",
" if stepSize < 1e-4:\n",
" return True\n",
" else:\n",
" return False\n",
"\n",
"def conjugateGradient(matA_in, vecB_in, rounding = 3, initialGuess = None, trace = False, convergencePlot = False):\n",
" # If we don't have an initial guess (which we should have) generate one ...\n",
" if not initialGuess:\n",
" # initialGuess = Matrix(vecB_in.rows, 1, rounding = vecB_in.rounding, initList = [ 0.0 ] * vecB_in.rows)\n",
" initialGuess = Matrix(vecB_in.rows, 1, rounding = vecB_in.rounding)\n",
" for i in range(vecB_in.rows):\n",
" initialGuess.set(i, 0, matA_in.get(i,i,))\n",
" else:\n",
" initialGuess = initialGuess.copy()\n",
"\n",
" if convergencePlot:\n",
" convergencePlotData = { 'iteration' : [], 'r' : [] }\n",
"\n",
" # Calculate initial residuum and rename \"initialGuess\" to x_k\n",
" # k will be the iteration counter\n",
" x = initialGuess\n",
" r = vecB_in.sub(matA_in.multiply(x))\n",
" k = 0\n",
"\n",
" while not conjugateGradient_Terminate(r):\n",
" if convergencePlot:\n",
" convergencePlotData['iteration'].append(len(convergencePlotData['iteration']))\n",
" convergencePlotData['r'].append(r.transpose().multiply(r).get(0,0))\n",
"\n",
" if k == 0:\n",
" pold = None\n",
" p = r.copy()\n",
" else:\n",
" s = -1.0 * pold.transpose().multiply(matA_in.multiply(r)).get(0,0) / pold.transpose().multiply(matA_in.multiply(pold)).get(0,0)\n",
" p = r.add(pold.multiply(s))\n",
" spectralCoefficient = p.transpose().multiply(r).get(0,0) / p.transpose().multiply(matA_in.multiply(p)).get(0,0)\n",
" x = x.add(p.multiply(spectralCoefficient))\n",
" r = vecB_in.sub(matA_in.multiply(x))\n",
" k = k + 1\n",
" pold = p\n",
"\n",
" # Round our result as specified\n",
" res = [ 0.0 ] * x.rows\n",
" for i in range(x.rows):\n",
" res[i] = round(x.get(i, 0), rounding)\n",
"\n",
" if convergencePlot:\n",
" plt.plot(convergencePlotData['iteration'], convergencePlotData['r'])\n",
" plt.title(\"Convergence plot (Residuum)\")\n",
" plt.grid()\n",
" plt.xlabel(\"Iteration\")\n",
" plt.ylabel(\"Residuum b - Ax\")\n",
" return (res, r)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "c4ed539f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1.0, 2.0, 3.0, 4.0]\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAApJUlEQVR4nO3dZ5gc5Znu8f8zSWkUUUB5lEYkkyRAREmgABgvGGe82GAbgW2SQfiYtY93l7V39yyIZLKJxoCcwIsBWxJCARElgRAIMcoICaGcRmnScz5UDbSGCT0z3VPdNffvuvpSd8X77VI/9U51dZW5OyIiEj85UQcQEZH0UIEXEYkpFXgRkZhSgRcRiSkVeBGRmFKBFxGJKRV4kUYws0vMbF6Kl/mKmR2XymXWso5SMxtcx7iUtynJTG3M7AMz69nS624tVOBjxMwuMrMF4Yd5g5n93cxOizpXa2Vms83sBw1M8yVgt7u/Hb7+NzMrD7fhDjN71cxObm4Wdy9091XNXU4qufsB4GHg/0SdJa5U4GPCzK4Dbgf+E+gFDADuAc6PMNZBzCwv6gwZ6Arg8RrD/uDuhUB3YBbwpxZP1XKeBL5rZm2iDhJHKvAxYGadgZuAH7v70+6+x93L3f1v7n5DOE0bM7vdzD4OH7dXf6jMbIyZrTOz681sU9j7vzQcN8rMPjGz3IT1fdnMFofPc8zsZ2a20sy2mtkfzaxbOK7IzNzMvm9ma4GXzCzXzKaY2RYzW21mV4bT5FW3xcweCjOsN7NfVa+7+lCCmd1iZtvD+c9JyNXNzB4J27fdzP6aMO48M1uU0Cs+up73083sajNbFea82cxq/ayY2SlmNt/Mdob/nhIO/zVwOnBX2Bu/q5Z5C4AzgTm1LdvdK4AngL5m1iOJ92eomc0Js2wxsz/UaNPQ8PkhZvasme0yszeBIQnTVW+zvIRhn/4lEv6F8fu6pg+n/VX4Hpea2d/C9T0Rrm++mRUltHEdsB0YVdf2kKZTgY+Hk4G2wDP1TPNzgg/RscAxwInALxLGHwp0BvoC3wfuNrOu7v46sIegEFW7iKDnBXA1cAEwGuhD8GG9u8a6RwOHAxOBy4BzwhzHh/MmegyoAIYCxwETgMTDHCcBJQS92/8BHjIzC8c9DrQHjgR6ArcBmNnxBIcCLgcOAe4Hnm2g1/hlYGSY8XzgezUnCHdkzwN3hsu9FXjezA5x958DLwNXhodHrqxlHcOAqrDIfU64A/gOsJXgfW3o/fkPYDrQFegH/KaOtt0N7Ad6h+36XNua6ZvAxQT/l4YArwGPAN2ApcC/1ph+KcH/SUk1d8+oB8EHcRPwXpLTfx14H1gCPBl1/ojes28DnzQwzUrg3ITXE4E14fMxwD4gL2H8JmBU+PxXwMPh844EBX9g+HopcFbCfL2BciAPKAIcGJww/iXg8oTX48Jp8ggOLR0A2iWM/xYwK3x+CbAiYVz7cN5Dw/VWAV1rafu9wH/UGFYCjK7jvXLg7ITXPwJmJmSYFz6/GHizxryvAZeEz2cDP6hnm5xac7sB/waUATuASoLiPiYc19D78zvgAaBfHW0aCuSG2+ewhHH/mdCm6m2W+H/h03aE+X6fMO6g6cNpf54wfgrw94TXXwIW1cj2BPDLqD9HcXxkYg/+UeDsZCY0s2HAjcCp7n4kcG36YmW0rUD3Bo5x9wE+THj9YTjs02V4cEig2l6gMHz+JHBh2OO9EHjL3auXNRB4Jjz0sYOg4FcSFKNqH9XI8VEd4wYC+cCGhOXdT9Abr/ZJ9RN33xs+LQT6A9vcfTufNxC4vnqZ4XL712h/TYm5ar5XiW35sMawDwl6rsnYTrDDrOmP7t6F4D18DxgRDm/o/fkpYMCbZrbEzGrrmfcg2JnWbF8qbUx4vq+W14UHT05Hgh2apFjGFXh3nwtsSxxmZkPM7B9mttDMXjazw8JRlwF3V3+o3X1TC8fNFK8R/Ml9QT3TfExQIKoNCIc1yN3fJygC53Dw4RkICsU57t4l4dHW3dcnLiLh+QaCwwfV+tdY1gGge8KyOoU774Z8BHQzsy51jPt1jYzt3f2pepaXmKuu96rme1o9bXXbG7pU63LAzKzWHYK7byE4rPRvZtabBt4fd//E3S9z9z7hfPdUH3dPsJngEE/N9lXbE/7bPmHYoTXG1zWuqQ4H3knBcqSGjCvwdXgAuMrdRwCTCc4OASgGii04j/h1M0uq5x837r4T+CXBcfMLzKy9meWb2Tlm9j/hZE8BvzCzHmbWPZz+93UtsxZPEhxvP4ODz+q4D/i1mQ0ECJdf35k7fwSuMbO+YTH+9BQ5d99AcAx5ipl1suAL3CFmNrqhcOG8fycoal3D9p8Rjv4tcIWZnWSBDmb2RTOrrfdc7YZwOf2Ba4A/1DLNCwT//y4yszwz+wZwBPBcOH4jUOu552HmcuBFgu8o6prmA2Aa8NOG3h8z+5qZVe88txPsYCprLK8SeJpgp9HezI4AvpswfjPBDuqfLfhC/HskfAkLLALOMLMBFny5f2Nd2ZMR7ty6Aa83ZzlSu4wv8GZWCJwC/MnMFhH8Sdo7HJ1H8EXVGIJjkQ/W0YOLPXe/FbiO4IvTzQS9vSuBv4aT/ApYACwG3gXeCocl6ymC9/mlsGdZ7Q7gWWC6me0m+KCeVM9yfktQpBYDbxMUyQo+K0TfAQoIvlfZDvyZz7Z3Qy4mOL78AcF3CNcCuPsCgr/27gqXuYLgWHp9/hdYSFDQngceqjmBu28FzgOuJzhM9lPgvIT35w7gqxac0XNnHeu5P8xdn5uBSRb8IKi+9+cE4A0zKyXYJte4++palnclwWGSTwgOiT5SY/xlwA1hm44EXk1o8wyCnd1igvfnOZrnIuAxD86JlxQz98y74Ud4GtVz7n6UmXUCStz9cx9yM7sPeN3dHw1fzwR+5u7zWzKvNJ0Fpzne5+41D3VExswcGObuK1poffMI/kJ9uyXWlynC73TeAc5oxYdX0yrje/DuvgtYbWZfg+CApZlVn1L1V2BsOLw7wSGbjPq1nhzMzNqZ2bnhIY2+BKfM1Xd6Z+y5+2mtrbhD8EtWdz9MxT19Mq7Am9lTBF8aDrfgxzffJzgN8Ptm9g7B6ZDVx3inAVvN7H2CX/zdEP7ZLJnLgH8nOLzwNsFZN7+MNJFITGXkIRoREWm+jOvBi4hIamTUxZ+6d+/uRUVFTZp3z549dOjQIbWBIhKXtsSlHaC2ZKK4tAOa15aFCxducfcetY3LqAJfVFTEggULmjTv7NmzGTNmTGoDRSQubYlLO0BtyURxaQc0ry1mVucvkXWIRkQkplTgRURiSgVeRCSmVOBFRGJKBV5EJKZU4EVEYkoFXkQkprK+wO8vr+SBuSsp2VbZ8MQiIq1IRv3QqakenreGwpwyJrnz2f2XRURat6zvwbfNz+Wqs4ayYkcVs0p01VERkWpZX+ABvj6yPz3aGTdPW0ZVla6OKSICMSnw+bk5XDA0n6UbdvHCexuijiMikhFiUeABTu6Tx7Cehdw6YxkVlVVRxxERiVxsCnyOGddPKGbV5j088/b6qOOIiEQuNgUeYOKRh/KFvp25/cXlHKjQaZMi0rrFqsBb2Itfv2Mff5j/UdRxREQiFasCDzC6uAcnFnXjNy+tYF+ZevEi0nrFrsCbGZMnDmfz7gP87rU1UccREYlM7Ao8wImDunFGcQ/unbOS3fvLo44jIhKJWBZ4gMkTitmxt5wHX14ddRQRkUjEtsAf3a8LE4/sxUPzVrN9T1nUcUREWlxsCzzA9ROGs6esgvvmrIw6iohIi4t1gS/u1ZELju3LY6+tYeOu/VHHERFpUbEu8ADXjhtGRaVz10sroo4iItKiYl/gBx7Sga+N7M/U+Wv5aNveqOOIiLSY2Bd4gKvPGoqZccfM5VFHERFpMa2iwPfu3I6LRw3k6bfWsWJTadRxRERaRKso8AA/HDOEtvm53DZjWdRRRERaRKsp8N0L2/C9Uwfx/LsbeG/9zqjjiIikXasp8ACXnTGYTm3zuFW9eBFpBVpVge/cLp/LRw/hpQ82sfDDbVHHERFJq1ZV4AEuPbWI7oUF3DytBHfdoFtE4qvVFfj2BXn8aMxQXl+1jVdWbI06johI2qS9wJtZrpm9bWbPpXtdybropAH06dyWm6erFy8i8dUSPfhrgKUtsJ6ktc3P5eqzhvHORzt4cemmqOOIiKRFWgu8mfUDvgg8mM71NMVXRvSj6JD2TJleQlWVevEiEj+WzkMUZvZn4L+AjsBkdz+vlmkmAZMAevXqNWLq1KlNWldpaSmFhYWNmue1jyu4f/EBrji6DaP65DVpvenQlLZkori0A9SWTBSXdkDz2jJ27NiF7j6y1pHunpYHcB5wT/h8DPBcQ/OMGDHCm2rWrFmNnqeysson3jbHx9w8y8srKpu87lRrSlsyUVza4a62ZKK4tMO9eW0BFngdNTWdh2hOBf7JzNYAU4Ezzez3aVxfo+XkGNeNL2b1lj385a11UccREUmptBV4d7/R3fu5exHwTeAld//ndK2vqcYf0Ytj+nfhjheXc6CiMuo4IiIp0+rOg6/JzJg8oZiPd+7nyTfWRh1HRCRlWqTAu/tsr+UL1kxx2tDunDSoG3fPWsHesoqo44iIpESr78FD0Iu/YeJwtpSW8eira6KOIyKSEirwoZFF3Rg7vAf3z1nFzn3lUccREWk2FfgE108Yzs595Tz48qqoo4iINJsKfIKj+nbm3C8cysPzVrO19EDUcUREmkUFvobrxhezr7ySe2evjDqKiEizqMDXMLRnR758XD9+9/qHbNi5L+o4IiJNpgJfi2vHDcPd+c1LK6KOIiLSZCrwtejfrT3fOKE/f5z/EWu37o06johIk6jA1+GqM4eRm2Pc/qJu0C0i2UkFvg69OrXlu6cU8cyi9SzbuDvqOCIijaYCX48rRg+hQ0Eet05XL15Eso8KfD26dSjge6cN4h9LPuHddTujjiMi0igq8A34wemD6NI+n1uml0QdRUSkUVTgG9CpbT5XjB7CnGWbmb9mW9RxRESSpgKfhO+eXESPjm24+R8l1bcjFBHJeCrwSWhXkMuVY4fy5pptzF2+Jeo4IiJJUYFP0jdP7E/fLu2YMl29eBHJDirwSWqTl8s144axeN1Opi3ZGHUcEZEGqcA3woXH9WVwjw7cOqOEyir14kUks6nAN0Jebg4/GVfMso2lPPvO+qjjiIjUSwW+kb74hd4c3rsTt81YTnllVdRxRETqpALfSDk5xuQJxazdtpc/LVgXdRwRkTqpwDfBmYf15LgBXbhz5nL2l1dGHUdEpFYq8E1gZtwwYTif7NrP71//MOo4IiK1UoFvolOGdueUIYdw7+yV7DlQEXUcEZHPUYFvhskTh7N1TxmPvLI66igiIp+jAt8Mxw/oyrjDe3L/3FXs3FsedRwRkYM0WODN7HEz65zweqCZzUxvrOxx3fjh7N5fwf1zV0YdRUTkIMn04OcBb5jZuWZ2GTADuD2tqbLIEX06cd7RvXnklTVs3n0g6jgiIp9qsMC7+/3AD4D/BW4CznD3v6U7WDa5bnwxZZVV3DN7RdRRREQ+lcwhmouBh4HvAI8CL5jZMWnOlVUG9yjkK8f35YnX17J+x76o44iIAMkdovkKcJq7P+XuNwJXEBR6SXD1WcNwnN/MXB51FBERILlDNBe4+6aE128Cl6c1VRbq17U9F504gD8tXMfqLXuijiMikvxpkmZ2hJndZGbLgXuSmL6tmb1pZu+Y2RIz+/dmJc0CPz5zKPm5xu0vLos6iohI/QU+PCXyZ2b2DvA48CNgvLuPTGLZB4Az3f0Y4FjgbDMb1dzAmaxnx7Zccsognn3nY0o+2R11HBFp5eos8Gb2KvACkA981d1HALvdfU0yC/ZAafgyP3zE/i4ZV4weTGFBHlOml0QdRURaufp68JuBjkAvoEc4rFEF2sxyzWwRsAmY4e5vNCVkNunSvoAfnD6Y6e9vZNFHO6KOIyKtmNV3A+nwF6xfAb4FDAW6ABPDL1qTX4lZF+AZ4Cp3f6/GuEnAJIBevXqNmDp1amMW/anS0lIKCwubNG+q7atwfjpnLwM65XDDCe0aPX8mtaU54tIOUFsyUVzaAc1ry9ixYxfWedjc3ZN6AD2Bq4BXgY+SnS9h/n8FJtc3zYgRI7ypZs2a1eR50+GBOSt94P95zl9buaXR82ZaW5oqLu1wV1syUVza4d68tgALvI6amvRZNO6+yd1/4+6nAKc1NL2Z9Qh77phZO2Ac8EGy68t2F588kF6d2nDLtJLqHZyISItq0tUk3T2Zu1z0BmaZ2WJgPsEx+Oeasr5s1DY/lyvPHMaCD7czu2Rz1HFEpBVK2+WC3X2xux/n7ke7+1HuflO61pWpvjGyP/26tuOW6SVUVakXLyItS9eDT6OCvBx+Mq6YJR/v4h9LPok6joi0Mo0q8Gb2VrqCxNUFx/VlaM9CpkwvoVK9eBFpQY3twVtaUsRYbo5x3fhiVm7ewzNvr486joi0Io0t8M+nJUXMnX3koRzZpxO3v7iMsoqqqOOISCvRqALv7r9IV5A4y8kxJk8czrrt+/jDgo+ijiMirYS+ZG0hY4p7MHJgV34zczn7yyujjiMirYAKfAsxC3rxm3Yf4HevrYk6joi0AkkVeDMrMLOjzewLZlaQ7lBxNWrwIZw+rDv3zl7J7v3lUccRkZhL5p6sXwRWAncCdwErzOycdAeLq8kThrN9bzkPz1sTdRQRiblkevBTgLHuPsbdRwNjgdvSGyu+junfhQlH9OLBl1exfU9Z1HFEJMaSKfCb3H1FwutVBNd3lya6fsJwSssquG/uyqijiEiM1XdHpwvN7EJgiZm9YGaXmNl3gb8RXDxMmmj4oR35p2P68Nira9i0a3/UcUQkpurrwX8pfLQFNgKjgTEEd3rqmvZkMfeTccWUVzp3z1rR8MQiIk2QV9cId7+0JYO0NkXdO/D1kf148s21XHbGYPp1bR91JBGJGZ0HH6GrzhyGYdw5c3nUUUQkhlTgI9SnSzu+PWoAf164jpWbS6OOIyIxowIfsR+NGUrb/Fxum7Es6igiEjN1HoOvFt5X9TtAUeL07n512lK1Ij06tuHSU4u4e9ZKfjRmF0f06RR1JBGJiWR68C8QFPd3gYUJD0mRSacPoWPbPG6dURJ1FBGJkQZ78EBbd78u7Ulasc7t87n8jMHcMn0Zb63dHnUcEYmJZHrwj5vZZWbW28y6VT/SnqyVufTUQRzSoYBbpqkXLyKpkUyBLwNuBl7js8MzC9IZqjXq0CaPH40dyqsrt/L+Vl0vXkSaL5kCfx0w1N2L3H1Q+Bic7mCt0bdPGkDvzm35y7Iy3HWDbhFpnmQK/BJgb7qDCLTNz+WqM4excmcVM5fqem4i0jzJFPhKYJGZ3W9md1Y/0h2stfrayH70bG/cMr2Eqir14kWk6ZI5i+av4UNaQH5uDhcMLeCBxbt5/t0NfOmYPlFHEpEs1WCBd/fHWiKIfGZU71xmbyzkthnLOOeoQ8nL1Q+ORaTxkrll32ozW1Xz0RLhWqscM64bP5xVW/bw9Fvro44jIlkqmUM0IxOetwW+Bug8+DSbeGQvju7XmTtmLuf84/rQJi836kgikmUa7MG7+9aEx3p3vx04M/3RWjczY/KE4azfsY+pb34UdRwRyULJXGzs+ISXOQQ9+o5pSySfOn1Yd04c1I27Zq3g6yP7065AvXgRSV4y395NSXj8FzAC+Ho6Q0nAzLhh4nA27z7AY6+tiTqOiGSZZM6iGdsSQaR2JxR1Y3RxD+6dvZKLThpAp7b5UUcSkSyRzCGaX9Y23N1vSn0cqc3kCcP50l3zePDl1Vw3vjjqOCKSJZI5RLMn4VEJnENwffh6mVl/M5tlZkvNbImZXdOspK3YF/p15pyjDuWhl1exbU9Z1HFEJEskcxbNlITHr4ExQN8kll0BXO/uhwOjgB+b2RHNStuKXTe+mL3lldw3Z2XUUUQkSzTlJ5LtgQavJunuG9z9rfD5bmApye0YpBbDenXky8f25bFX17Bx1/6o44hIFrCGLktrZu8C1RPlAj2Am9z9rqRXYlYEzAWOcvddNcZNAiYB9OrVa8TUqVOTDp+otLSUwsLCJs2baepqy6a9Vdz48j5G98vjO0e2iSBZ47SGbZKN4tKWuLQDmteWsWPHLnT3kbWOdPd6H8DAhEdfIK+heWrMX0hwk5ALG5p2xIgR3lSzZs1q8ryZpr62/MvTi33Ijc/72q17Wi5QE7WWbZJt4tKWuLTDvXltARZ4HTW1zkM0Cbfm253w2Ad0SvaWfWaWD/wFeMLdn05qdyT1uurMYeTmGLe/uDzqKCKS4eo7Bl99a76FwGZgGbA8fL6woQWbmQEPAUvd/dbmRxWAQzu35eJRA3nm7XWs2LQ76jgiksHqLPD+2a35pgFfcvfu7n4IcB6QTG/8VOBi4EwzWxQ+zk1J6lbuh2OG0C4/l1tnLIs6iohksGTOojnB3V+ofuHufwdGNzSTu89zd3P3o9392PDxQkPzScMOKWzD908bxAvvfsJ763dGHUdEMlQyBX6Lmf3CzIrMbKCZ/RzYmu5gUr8fnDGYzu3ymTK9JOooIpKhkinw3yI4NfIZglv39QyHSYQ6tc3n8tGDmVWymQVrtkUdR0QyUDK/ZN3m7te4+3Hh4xp3V0XJAJecUkT3wjbcPK2k+pRUEZFP1Xea5O3hv38zs2drPlosodSpfUEeV44dwhurtzFvxZao44hIhqnvapKPh//e0hJBpGm+ddIAfvvyam6ZVsJpQ7sTnJ0qIlL/aZILw3/nVD+AxcD28LlkgDZ5uVx91lDeWbeT6e9vjDqOiGSQBo/Bm9lsM6v+9eo7wCNmph8uZZCvHN+PQd07cOv0ZVRW6Vi8iASSOYumswcXCLsQeMTdRwDj0htLGiMvN4efjC+mZONunlv8cdRxRCRDJFPg88ysN8F9WJ9Lcx5povO+0JvDDu3IbTOWUV5ZFXUcEckAyRT4mwguV7DS3eeb2WCCa9JIBsnJMa6fMJw1W/fyl4Xroo4jIhkgmfPg/xRebuCH4etV7v6V9EeTxhp3eE+O6d+FO2YuZ395ZdRxRCRiyXzJWmxmM83svfD10Wb2i/RHk8YyM346cTgbdu7nyTfWRh1HRCKWzCGa3wI3AuUA7r4Y+GY6Q0nTnTq0OycPPoR7Zq9gb1lF1HFEJELJFPj27v5mjWGqHBls8sThbCkt45FX1kQdRUQilOzVJIcQ3pfVzL4KbEhrKmmWEQO7cuZhPbl/zkp27iuPOo6IRCSZAv9j4H7gMDNbD1wLXJHOUNJ8108oZtf+Cn47d1XUUUQkIsmcRbPK3ccRXDL4MGAMcFqac0kzHdmnM188ujcPv7KaLaUHoo4jIhGo72qSnczsRjO7y8zGA3uB7wIrCH70JBnuJ+OK2V9eyb2zV0YdRUQiUF8P/nFgOPAucBkwHfgacIG7n98C2aSZhvYs5MLj+/H46x+yYee+qOOISAurr8APdvdL3P1+gjs4jQTOc/dFLZJMUuKas4bh7tw5c0XUUUSkhdVX4D89/cLdK4HV7r47/ZEklfp3a8+3ThzAnxZ8xIdb90QdR0RaUH0F/hgz2xU+dgNHVz83s10tFVCa78qxQ8nLNW5/UZcQEmlN6rvhR667dwofHd09L+F5p5YMKc3Ts1NbvntyEX9dtJ6ST/RHmEhrkcx58BIDV4weQoeCPG6dURJ1FBFpISrwrUTXDgX84PRBTFuykcXrdkQdR0RagAp8K/L90wbRtX0+t0xfFnUUEWkBKvCtSMe2+Vwxeghzl23mjVVbo44jImmmAt/KfOfkInp0bMMt00tw1w26ReJMBb6VaVeQy9VnDmX+mu3MWbY56jgikkYq8K3QN04YQL+u7ZgyfZl68SIxpgLfChXk5XDNWcN4d/1Opi35JOo4IpImKvCt1JeP68vgHh24ZfoyKqvUixeJIxX4ViovN4frxw9nxaZS/nfR+qjjiEgapK3Am9nDZrbJzN5L1zqkec456lCO6N2J219cTnllVdRxRCTF0tmDfxQ4O43Ll2bKyTEmTyxm7ba9/HHBR1HHEZEUS1uBd/e5wLZ0LV9SY+zwnhw/oAt3zlzO/vLKqOOISApZOk+TM7Mi4Dl3P6qeaSYBkwB69eo1YurUqU1aV2lpKYWFhU2aN9O0dFuWbq3k/83fzzeHF3D2oPyULVfbJDPFpS1xaQc0ry1jx45d6O4jax3p7ml7AEXAe8lOP2LECG+qWbNmNXneTBNFW77929f9uJum++795SlbprZJZopLW+LSDvfmtQVY4HXUVJ1FIwBMnjicbXvKeGTe6qijiEiKqMALAMf278K4w3vxwNxV7NhbFnUcEUmBdJ4m+RTwGjDczNaZ2ffTtS5JjesnFFNaVsH9c1dFHUVEUiCdZ9F8y917u3u+u/dz94fStS5JjcN7d+JLR/fh0VfWsGn3/qjjiEgz6RCNHOQn44spq6zinlkro44iIs2kAi8HGdS9A189vh9PvrGW9Tv2RR1HRJpBBV4+5+pxwwC488XlEScRkeZQgZfP6dulHRedNIA/v7WO1Vv2RB1HRJpIBV5q9eOxQynIzeG2GbpBt0i2UoGXWvXo2IZLTi3i2Xc+ZumGXVHHEZEmUIGXOl1+xmA6tsljynT14kWykQq81KlL+wImnTGYF5du5O2126OOIyKNpAIv9br0tEF061CgXrxIFlKBl3oVtsnjR2OGMG/FFl5duSXqOCLSCCrw0qB/HjWQXp3acMu0kurLQItIFlCBlwa1zc/l6rOG8dbaHcwq2RR1HBFJkgq8JOXrI/szoFt7bpm2jKoq9eJFsoEKvCQlPzeHa8cN4/0Nu/j7e59EHUdEkqACL0k7/9i+DOtZyJQZJVRUVkUdR0QaoAIvScvNMa4bX8yqzXt45u31UccRkQaowEujnH3UoXyhb2fumLmcsgr14kUymQq8NIqZcf2EYtZt38cf5q+NOo6I1EMFXhptdHEPTijqyp0vrWBfWWXUcUSkDirw0mhmxuQJw9m8+wC/e21N1HFEpA4q8NIkJw0+hDOKe3DvnJXs3l8edRwRqYUKvDTZ5AnF7NhbzkPzVkcdRURqoQIvTXZ0vy5MPLIXD768mu17yqKOIyI1qMBLs1w/YTh7yiq4b87KqKOISA0q8NIsxb06csGxfXnstTVs2rU/6jgikkAFXprt2nHDqKh07pq1IuooIpJABV6abeAhHfjayP489eZaPtq2N+o4IhJSgZeUuPqsoZgZd8xcHnUUEQmpwEtK9O7cjotHDeTpt9axYlNp1HFEBBV4SaEfjhlC2/xcbntRN+gWyQQq8JIy3Qvb8L1TB/H84g0s+Xhn1HFEWj0VeEmpy84YTKe2eUyZrl68SNRU4CWlOrfL5/LRQ3jpg02s2K4rTYpEKS/qABI/l55axCOvrObRJQfYWPA+BXk5FOTmBv/m5VCQawnPg+H54bA2Nab9dHjCsNwci7qJIlkhrQXezM4G7gBygQfd/b/TuT7JDO0L8vi/5x3BL/6yiMdf/5DySqeyylO2/Bwj3DnkUJCXS5uEHcFnw3PIzw13GDWGVU/XppZh1dMV1Bi2ckcl3dfv/Gx54fD86mlzc8jRjkcyTNoKvJnlAncD44F1wHwze9bd30/XOiVznH9sXzrvWM6YMWMAqKxyyiqqKKuo4kBlJeWVn70uq6iirLKSsgqnrPLgYeUVzoHEYRVVlFdWfTrdgcRhFVUHzb+7vIKtNYZVT1e9zEZ5fV69o/Nz7eCi38BOIz9hJ/O5HUbeZzun/Nwa89azzJrD8nIMM+14Wqt09uBPBFa4+yoAM5sKnA+owLdCuTlGu4Jc2hXkAvlRxwHA3Smv9M/tHA7UssNY8PYiDjviqISdz2c7jAO1DPtsZ/b5YaUHKj63M0qct7wydX/tmBH8JZNQ9MsOHKDt6y81e7nNncZoeCH1LWPfvn20nz8riWXUv56kdn8NtiWJRdSTI6d8H2FfKKXSWeD7Ah8lvF4HnFRzIjObBEwC6NWrF7Nnz27SykpLS5s8b6aJS1vi0g6AIe3202bzB7SpbaQR7LOatN/Koea5DlXuVFTx2SN8XV4FFVVOeRVU1nhdET6vnqe8xjLKq5yKqioqqqooa1NFfn7dN2nxFOxfGlqENzhFwwspz6kiL+9A0plqXUXzYyS3ngbG5xdUpuWzks4CX9vu6nPtdPcHgAcARo4c6WOauBubPXs2TZ0308SlLXFpB6gtmSgu7YD0tSWdp0muA/onvO4HfJzG9YmISIJ0Fvj5wDAzG2RmBcA3gWfTuD4REUmQtkM07l5hZlcC0whOk3zY3Zeka30iInKwtJ4H7+4vAC+kcx0iIlI7XapARCSmVOBFRGJKBV5EJKZU4EVEYso8FT9bSxEz2wx82MTZuwNbUhgnSnFpS1zaAWpLJopLO6B5bRno7j1qG5FRBb45zGyBu4+MOkcqxKUtcWkHqC2ZKC7tgPS1RYdoRERiSgVeRCSm4lTgH4g6QArFpS1xaQeoLZkoLu2ANLUlNsfgRUTkYHHqwYuISAIVeBGRmMqqAm9mZ5tZiZmtMLOf1TLezOzOcPxiMzs+ipzJSKItY8xsp5ktCh+/jCJnQ8zsYTPbZGbv1TE+m7ZJQ23Jlm3S38xmmdlSM1tiZtfUMk1WbJck25It26Wtmb1pZu+Ebfn3WqZJ7XZx96x4EFxyeCUwGCgA3gGOqDHNucDfCe4mNQp4I+rczWjLGOC5qLMm0ZYzgOOB9+oYnxXbJMm2ZMs26Q0cHz7vCCzL4s9KMm3Jlu1iQGH4PB94AxiVzu2STT34T2/i7e5lQPVNvBOdD/zOA68DXcysd0sHTUIybckK7j4X2FbPJNmyTZJpS1Zw9w3u/lb4fDewlOAeyYmyYrsk2ZasEL7XpeHL6rv41jzLJaXbJZsKfG038a65oZOZJhMkm/Pk8M+5v5vZkS0TLeWyZZskK6u2iZkVAccR9BYTZd12qactkCXbxcxyzWwRsAmY4e5p3S5pveFHiiVzE++kbvSdAZLJ+RbBNSZKzexc4K/AsHQHS4Ns2SbJyKptYmaFwF+Aa919V83RtcySsdulgbZkzXZx90rgWDPrAjxjZke5e+J3PindLtnUg0/mJt7ZcqPvBnO6+67qP+c8uDNWvpl1b7mIKZMt26RB2bRNzCyfoCA+4e5P1zJJ1myXhtqSTdulmrvvAGYDZ9cYldLtkk0FPpmbeD8LfCf8JnoUsNPdN7R00CQ02BYzO9TMLHx+IsG22triSZsvW7ZJg7Jlm4QZHwKWuvutdUyWFdslmbZk0XbpEfbcMbN2wDjggxqTpXS7ZM0hGq/jJt5mdkU4/j6C+7+eC6wA9gKXRpW3Pkm25avAD82sAtgHfNPDr9kziZk9RXAWQ3czWwf8K8GXR1m1TSCptmTFNgFOBS4G3g2P9wL8CzAAsm67JNOWbNkuvYHHzCyXYCf0R3d/Lp01TJcqEBGJqWw6RCMiIo2gAi8iElMq8CIiMaUCLyISUyrwIiIxpQIvsWRmpeG/RWZ2UYqX/S81Xr+ayuWLpIoKvMRdEdCoAh+ep1yfgwq8u5/SyEwiLUIFXuLuv4HTw+uE/yS82NPNZjY/vN725fDpNcVnmdmTwLvhsL+a2cLw2t2TwmH/DbQLl/dEOKz6rwULl/2emb1rZt9IWPZsM/uzmX1gZk9U//JSJJ2y5pesIk30M2Cyu58HEBbqne5+gpm1AV4xs+nhtCcCR7n76vD199x9W/iz8vlm9hd3/5mZXenux9ayrguBY4FjgO7hPHPDcccBRxJcV+QVgl9ozkt1Y0USqQcvrc0Egmt9LCK47OwhfHblwTcTijvA1Wb2DvA6wQWgGrpC4WnAU+5e6e4bgTnACQnLXufuVcAigkNHImmlHry0NgZc5e7TDhpoNgbYU+P1OOBkd99rZrOBtkksuy4HEp5Xos+etAD14CXudhPc6q3aNIILU+UDmFmxmXWoZb7OwPawuB9GcPu0auXV89cwF/hGeJy/B8EtAN9MSStEmkC9CIm7xUBFeKjlUeAOgsMjb4VfdG4GLqhlvn8AV5jZYqCE4DBNtQeAxWb2lrt/O2H4M8DJBPfYdeCn7v5JuIMQaXG6mqSISEzpEI2ISEypwIuIxJQKvIhITKnAi4jElAq8iEhMqcCLiMSUCryISEz9f9aYPWf3reGAAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"x,r = conjugateGradient(ex3_A, ex3_B, convergencePlot = True)\n",
"print(x)"
]
},
{
"cell_type": "markdown",
"id": "73a5773c",
"metadata": {},
"source": [
"Again this of course matches our expected output well."
]
},
{
"cell_type": "markdown",
"id": "218a681a",
"metadata": {},
"source": [
"## Generalized minimal residual (GMRES)\n",
"\n",
"GMRES and it's variants are the currently most flexible and generic algorithms in use to solve massive amounts of linear equations. This is mainly also due to it's capability to work with matrices that are not symmetric positiv definit. In contrast to steepest decent and conjugate gradient methods GMRES on the other hand performs a general minimum search in the given Krylov subspace at the specific iteration by solving a local least squares problem. The result is refined at each iteration in each new subspace that grows at most one dimension per iteration. On the other hand it's not a typical iterative method since it has to perform a least squares minimization on every iteration - but it can keep intermediate results while performing QR decomposition. One can implement this particular efficient when using QR decomposition with Givens rotations.\n",
"\n",
"Since each subspace is a subspace of the previous ones again the residuum always is guaranteed to decrease with every iteration. When one has expanded the subspace over the whole space spanned by one's coefficient matrix $A$ the method yields the exact solution - the hope usually is that the approximation is good enough even after a few steps though. This algorithm works pretty well when also performing a continuous QR decomposition by iteratively applying Givens rotations to the problem. In contrast to solving the least squares problem directly to determine the solutions this saves much computational power in many cases by solving the problem in a reduced (Krylov) subspace. See [the blog post](https://www.tspi.at/2022/01/02/krylovsubspacelinearsolver.html#generalized-minimal-residual-method-gmres) for details.\n",
"\n",
"Using [Givens rotations for the least squares problem](https://www.tspi.at/2021/12/18/leastsquaresqr.html) is of particular interest for this method since the Arnoldi iteration yields a upper Hessenberg matrix on each iteration and thus an efficient implementation is possible.\n",
"\n",
"$$\n",
"A \\vec{q}_k = \\sum_{i=1}^{k+1} h_{i,k} \\vec{q}_i\n",
"$$\n",
"\n",
"$$\n",
"A Q^k = Q^{k+1} H^k\n",
"$$\n",
"\n",
"In this case $H \\in \\mathfrak{R}^{k+1 \\times k}$ is an upper Hessenberg matrix - i.e. all elements below the first sub diagonal are zero. The upper Hessenberg structure makes it pretty simple - at the $k$ step one only needs to apply $k$ givens rotations to the _last column_.\n",
"\n",
"$$\n",
"G(1,2,\\theta_1)^T * G(2,3,\\theta_2)^T, \\ldots, G(k-1, k, \\theta_{k-1})^T\n",
"$$\n",
"\n",
"In the end the transform $G(k, k+1, \\theta_k)^T$ has to be applied only to the result since this has no effect on the upper previous columns due to upper hessenberg structure. At step $k$ the cost is only $O(k)$ to update the QR factorization.\n",
"\n",
"### Elements used\n",
"\n",
"Lets sum up which elements are going to be used for GMRES. Note that - in contrast to CG and steepest descent - the size of the matrices (at least the used ones) change from iteration to iteration. Usually one selects a maximum number of iteration that one's ready to take - this then also determines the maximum size of the matrices so one can preallocate them for all allowed iterations. Especially for large scale applications it's usually not practicable to iterate over the whole number of possible iterations - since then in the end one could have done a direct QR decomposition for the whole system which is usually prohibited by the size of the problems. In this case one might use restarting - and restart processing of the problem with a new better initial guess for the solution that one can reuse from the previous iteration.\n",
"\n",
"$$\n",
"A * \\vec{x} = \\vec{b}\n",
"$$\n",
"\n",
"This is our initial problem with the coefficient matrix $A \\in \\mathfrak{R}^{n \\times n}$ and the constant vector $\\vec{b} \\in \\mathfrak{R}^{n}$. The vector $\\vec{x} \\in \\mathfrak{R}^n$ will contain our solutions (result of the algorithm).\n",
"\n",
"The matrix $A$ will iteratively be decomposed like during a $QR$ decomposition - into an orthogonal and an upper triangular part. The ortohonal part will be formed by the Krylov subspace vectors that are orthonormalized by using a Gram-Schmidt like procedure (this is also called Arnoldi iteration). These vectors will be $\\vec{q}_i \\in \\mathfrak{R}^n$. The matrix $Q$ is formed by appending one column after each other - at the $k$'th iteration $Q^{(k)} = [ \\vec{q}_1 \\mid \\vec{q}_2 \\mid \\ldots \\vec{q}_k ]$. This means that $Q \\in \\mathfrak{R}^{n \\times (m+1)}$ where $m$ is the maximum number of iterations we are taking with $m \\leq n$. $Q^{(k)} \\in \\mathfrak{R}^{n \\times k}$\n",
"\n",
"For least squares we have to store our approximated solution inside our subspace $\\vec{\\lambda}^{(k)} \\in \\vec{R}^{k}$. This is related to our real solution via $\\vec{x}^{(k)} = Q^{(k)} * \\vec{\\lambda}^{(k)}$\n",
"\n",
"We're also going to build and store our Hessenberg matrix incrementally. $H^{(k)} \\in \\mathfrak{R}^{(k+1) \\times k}$. This matrix can reach a maximum dimension of $H \\in \\mathfrak{R}^{(m+1) \\times m}$.\n",
"\n",
"Since we're going to apply Givens rotations to transform our Hessenberg matrix to upper diagonal form we also have to store the parameters for our Givens rotations since we want to re-apply them incrementally to the new columns. For this we only need to store the evaluated cosine ($c$) and sine ($s$) parts - with $\\vec{c} \\in \\mathfrak{R}^{m}$ and $\\vec{s} \\in \\mathfrak{R}^{m}$.\n",
"\n",
"The building of our Hessenberg matrix can be done by calculating vectors $\\vec{h}^{(k)} \\in \\mathfrak{(k+1)}$ that will be zero extended after each iteration. Note that $h_{(k+1), k}$ is also equal to our residuum after each iteration so we don't have to calculate that explicitly again."
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "1af50771",
"metadata": {},
"outputs": [],
"source": [
"def generalizezMinimalResidual(matA_in, vecB_in, x0_in = None, errorThreshold = 1e-6, maximumIterations = None, rounding = 3, trace = False, convergencePlot = False):\n",
" iterationCounter = 0\n",
"\n",
" if convergencePlot:\n",
" convergencePlotData = { 'iteration' : [], 'r' : [] }\n",
"\n",
" if maximumIterations:\n",
" if (maximumIterations > matA_in.cols) or (maximumIterations < 1):\n",
" maxIter = matA_in.cols\n",
" else:\n",
" maxIter = maximumIterations\n",
" else:\n",
" maxIter = matA_in.cols\n",
"\n",
" # Allocate our storage\n",
"\n",
" # Q is our orthogonal matrix that is built by concatenating the columns of our\n",
" # Krylov subspace base vectors\n",
" Q = Matrix(matA_in.rows, 1, storageRows = matA_in.rows, storageCols = maxIter+1)\n",
"\n",
" # Hessenberg matrix storage\n",
" H = Matrix(1, 0, storageRows = maxIter+1, storageCols = maxIter)\n",
"\n",
" # Storage for Givens rotation parameters (store sin(theta) and cos(theta) instead of angle)\n",
" s = [ 0.0 ] * maxIter\n",
" c = [ 0.0 ] * maxIter\n",
"\n",
" if x0_in:\n",
" # Calculate first residuum\n",
" r = vecB_in.sub(matA_in.multiply(x0_in))\n",
" else:\n",
" # If initial guess is zero the initial residuum is b\n",
" r = vecB_in.copy()\n",
"\n",
" # Our first base vector in our normalized Krylov subspace is just the initial residuum\n",
" # This will allow us to use the vector e_1 * |b| as our starting \"solution\"\n",
"\n",
" q = r.unitVector()\n",
" Q.setColumnVector(0, q)\n",
"\n",
" rNorm = r.norm()\n",
" error = rNorm\n",
"\n",
" if convergencePlot:\n",
" convergencePlotData['iteration'].append(len(convergencePlotData['iteration']))\n",
" convergencePlotData['r'].append(error)\n",
"\n",
"\n",
" # Build or local result vector. In our first Krylov space dimension this is just a vector\n",
" # with the same length as the original b but with direction along the first axis.\n",
" # We zero extend the vector to match the Q dimensions\n",
" b = Matrix(vecB_in.rows+1, vecB_in.cols)\n",
" b.set(0, 0, rNorm)\n",
"\n",
" if trace:\n",
" print(\"Before iteration:\\n------------\\nQ:\\n{}\\nH:\\n{}\\nb:\\n{}\\nError: {}\\n\".format(Q, H, b, error))\n",
"\n",
" for k in range(maxIter):\n",
" iterationCounter = iterationCounter + 1\n",
"\n",
" # Extend our Hesse matrix storage size\n",
" H.setSize(H.rows+1, H.cols+1)\n",
"\n",
" # Perform Arnoldi iteration. First build the next base vector by projection\n",
" q = matA_in.multiply(q)\n",
"\n",
" # And orthogonalize\n",
" for i in range(k+1):\n",
" # Project new basis vector onto previous one ...\n",
" qi = Q.getColumnVector(i)\n",
" hik = qi.transpose().multiply(q).get(0,0)\n",
" H.set(i, k, hik)\n",
" # ... and subtract the scaled part (this removes any components along the previous basis vector)\n",
" # which is the essence of Gram-Schmidt orthogonalization\n",
" q = q.sub(qi.multiply(hik))\n",
"\n",
" # Now we normalize our new vector and append to the Q matrix\n",
" h_knext_k = q.norm()\n",
"\n",
" error = h_knext_k\n",
" if error < 0:\n",
" error = -1.0 * error\n",
" if convergencePlot:\n",
" convergencePlotData['iteration'].append(len(convergencePlotData['iteration']))\n",
" convergencePlotData['r'].append(error)\n",
"\n",
" H.set(k+1, k, h_knext_k)\n",
" q = q.unitVector()\n",
" Q.setSize(Q.rows, Q.cols+1)\n",
" Q.setColumnVector(k+1, q)\n",
" \n",
" # Solve least squares using Givens rotations\n",
" # First apply the previous rotations\n",
"\n",
" hk = H.getColumnVector(k)\n",
" for i in range(k):\n",
" rotationMatrix = GivensRotationMatrixFromCS(hk.rows, i, i+1, c[i], s[i])\n",
" hk = rotationMatrix.multiply(hk)\n",
"\n",
" # And generate our new Givens rotation matrix for our new dimension ...\n",
" if not q.isZero():\n",
" x1 = hk.get(k, 0)\n",
" x2 = hk.get(k+1, 0)\n",
" s[k] = 1.0 / math.sqrt(x1 * x1 / (x2 * x2) + 1)\n",
" c[k] = s[k] * x1 / x2\n",
" rotationMatrix = GivensRotationMatrixFromCS(hk.rows, k, k+1, c[k], s[k])\n",
"\n",
" # Rotate our expected result into our subspace. In our implementation this means we\n",
" # require a second matrix - note that this might be solved by pre-allocating the whole\n",
" # matrix and then filling with zero in a real-world implementation - or by just applying\n",
" # the equations and not building a matrix since it's only a Givens rotation ...\n",
" rotationMatrix2 = GivensRotationMatrixFromCS(b.rows, k, k+1, c[k], s[k])\n",
" b = rotationMatrix2.multiply(b)\n",
"\n",
" # Last rotation of our new dimension to maintain upper triagonal structure\n",
" hk = rotationMatrix.multiply(hk)\n",
" H.setColumnVector(k, hk)\n",
"\n",
" if trace:\n",
" print(\"Iteration {}:\\n------------\\nQ:\\n{}\\nH:\\n{}\\nb:\\n{}\\nError: {}\\n\\nSines: {}\\nCosines: {}\".format(k, Q, H, b, error, s,c))\n",
"\n",
" # For trace also generate intermediate results so one can watch how the\n",
" # results converge to the real values. Usually one is not required to perform\n",
" # this matrix-vector product which saves much time. We can avoid this since\n",
" # q.norm() is a measure for the error of the next step\n",
"\n",
" xResult = [ 0.0 ] * (H.rows - 1)\n",
" for row in range(H.rows - 2, -1, -1):\n",
" t = b.get(row, 0)\n",
" for i in range(row + 1, H.rows-1):\n",
" t = t - xResult[i] * H.get(row, i)\n",
" t = t / H.get(row, row)\n",
" xResult[row] = t\n",
"\n",
" Q.setSize(Q.rows, Q.cols - 1)\n",
" xResult = Matrix(len(xResult), 1, initList = xResult)\n",
" xResult = Q.multiply(xResult)\n",
" Q.setSize(Q.rows, Q.cols + 1)\n",
"\n",
" print(\"Estimated result at iteration {}:\\n{}\".format(k, xResult))\n",
"\n",
" # If our error is smaller than a specified threshold we're done. Else\n",
" # we iterate up to the specified maximum number of iterations\n",
" if error < errorThreshold:\n",
" break\n",
"\n",
" # Perform least squares using QR decomposition. Note that we already have\n",
" # an upper triangular form of H so we solve H*lambda = b (b inside our projection\n",
" # space obviously) by backsubstitution ...\n",
"\n",
" xResult = [ 0.0 ] * (H.rows - 1)\n",
" for row in range(H.rows - 2, -1, -1):\n",
" t = b.get(row, 0)\n",
" for i in range(row + 1, H.rows-1):\n",
" t = t - xResult[i] * H.get(row, i)\n",
" t = t / H.get(row, row)\n",
" xResult[row] = t\n",
"\n",
" # Then project our local solution lambda back through our orthogonal projection\n",
" # into the space the problem had been defined in\n",
" Q.setSize(Q.rows, Q.cols - 1)\n",
" xResult = Matrix(len(xResult), 1, initList = xResult)\n",
" xResult = Q.multiply(xResult)\n",
" # Restore size\n",
" Q.setSize(Q.rows, Q.cols + 1)\n",
"\n",
" if convergencePlot:\n",
" plt.plot(convergencePlotData['iteration'], convergencePlotData['r'])\n",
" plt.title(\"Convergence plot (Error)\")\n",
" plt.grid()\n",
" plt.xlabel(\"Iteration\")\n",
" plt.ylabel(\"Error h_(k+1,k)\")\n",
" plt.show()\n",
"\n",
" return (xResult, error, iterationCounter)\n"
]
},
{
"cell_type": "markdown",
"id": "eed8a4d4",
"metadata": {},
"source": [
"To test the algorithm we're going to run it on the same system of equations as shown before. For example 3 the whole trace output is printed to allow to compare numerical values with other implementations and to see how the results evolve."
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "900e7502",
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Before iteration:\n",
"------------\n",
"Q:\n",
"0.19354839\t\n",
"0.0\t\n",
"0.4516129\t\n",
"0.87096774\t\n",
"\n",
"H:\n",
"\n",
"\n",
"b:\n",
"248.0\t\n",
"0.0\t\n",
"0.0\t\n",
"0.0\t\n",
"0.0\t\n",
"\n",
"Error: 248.0\n",
"\n",
"Iteration 0:\n",
"------------\n",
"Q:\n",
"0.19354839\t0.10030592\t\n",
"0.0\t-0.93321574\t\n",
"0.4516129\t-0.31490809\t\n",
"0.87096774\t0.14099646\t\n",
"\n",
"H:\n",
"53.0935\t\n",
"-0.0\t\n",
"\n",
"b:\n",
"245.4195\t\n",
"-35.6826\t\n",
"0.0\t\n",
"0.0\t\n",
"0.0\t\n",
"\n",
"Error: 7.639176770044404\n",
"\n",
"Sines: [0.14388146788245634, 0.0, 0.0, 0.0]\n",
"Cosines: [0.9895949288471468, 0.0, 0.0, 0.0]\n",
"Estimated result at iteration 0:\n",
"0.8947\t\n",
"0.0\t\n",
"2.0875\t\n",
"4.026\t\n",
"\n",
"Iteration 1:\n",
"------------\n",
"Q:\n",
"0.19354839\t0.10030592\t0.90324337\t\n",
"0.0\t-0.93321574\t-0.04244795\t\n",
"0.4516129\t-0.31490809\t0.26251434\t\n",
"0.87096774\t0.14099646\t-0.33680232\t\n",
"\n",
"H:\n",
"53.0935\t10.1503\t\n",
"-0.0\t16.739\t\n",
"0.0\t0.0\t\n",
"\n",
"b:\n",
"245.4195\t\n",
"-35.6396\t\n",
"1.7512\t\n",
"0.0\t\n",
"0.0\t\n",
"\n",
"Error: 0.8214917293418378\n",
"\n",
"Sines: [0.14388146788245634, 0.04907709208536911, 0.0, 0.0]\n",
"Cosines: [0.9895949288471468, 0.9987949934958845, 0.0, 0.0]\n",
"Estimated result at iteration 1:\n",
"0.7599\t\n",
"1.9869\t\n",
"2.9418\t\n",
"4.0803\t\n",
"\n",
"Iteration 2:\n",
"------------\n",
"Q:\n",
"0.19354839\t0.10030592\t0.90324337\t-0.36974265\t\n",
"0.0\t-0.93321574\t-0.04244795\t-0.35631504\t\n",
"0.4516129\t-0.31490809\t0.26251434\t0.79248289\t\n",
"0.87096774\t0.14099646\t-0.33680232\t-0.32909092\t\n",
"\n",
"H:\n",
"53.0935\t10.1503\t0.1205\t\n",
"-0.0\t16.739\t1.1507\t\n",
"0.0\t0.0\t6.8781\t\n",
"0.0\t0.0\t-0.0\t\n",
"\n",
"b:\n",
"245.4195\t\n",
"-35.6396\t\n",
"1.74\t\n",
"-0.1974\t\n",
"0.0\t\n",
"\n",
"Error: 0.7751625296736495\n",
"\n",
"Sines: [0.14388146788245634, 0.04907709208536911, 0.11270514830690037, 0.0]\n",
"Cosines: [0.9895949288471468, 0.9987949934958845, 0.9936284766174525, 0.0]\n",
"Estimated result at iteration 2:\n",
"0.9872\t\n",
"1.9924\t\n",
"3.015\t\n",
"3.995\t\n",
"\n",
"Iteration 3:\n",
"------------\n",
"Q:\n",
"0.19354839\t0.10030592\t0.90324337\t-0.36974265\t0.12451571\t\n",
"0.0\t-0.93321574\t-0.04244795\t-0.35631504\t-0.81133475\t\n",
"0.4516129\t-0.31490809\t0.26251434\t0.79248289\t-0.07641849\t\n",
"0.87096774\t0.14099646\t-0.33680232\t-0.32909092\t0.56603178\t\n",
"\n",
"H:\n",
"53.0935\t10.1503\t0.1205\t-0.021\t\n",
"-0.0\t16.739\t1.1507\t0.0292\t\n",
"0.0\t0.0\t6.8781\t1.8483\t\n",
"0.0\t0.0\t-0.0\t9.4229\t\n",
"0.0\t0.0\t0.0\t-0.0\t\n",
"\n",
"b:\n",
"245.4195\t\n",
"-35.6396\t\n",
"1.74\t\n",
"-0.1974\t\n",
"0.0001\t\n",
"\n",
"Error: 0.005742327708394915\n",
"\n",
"Sines: [0.14388146788245634, 0.04907709208536911, 0.11270514830690037, 0.0006049092059919701]\n",
"Cosines: [0.9895949288471468, 0.9987949934958845, 0.9936284766174525, 0.9999998170424097]\n",
"Estimated result at iteration 3:\n",
"1.0\t\n",
"2.0\t\n",
"3.0\t\n",
"4.0\t\n",
"\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEWCAYAAABi5jCmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAo/UlEQVR4nO3deXxddZ3/8dcna5d0X9IlSVukbGVr0wZQllZEEBmLoJBUHRz1h86Ao6MzCo4/fzojDo9xdNRRZkYHRmaUxrJXFgExFZGlG6W0lNLSQpvSfU+3JM3n98c5wZvbJHfJvTk3ue/n43EeOcv3nPvOSXI/Oed7zznm7oiIiLQriDqAiIjkFhUGERHpQIVBREQ6UGEQEZEOVBhERKQDFQYREelAhUEkx5nZJ83s2Qxv849mNj2T20zx9c82s+eien3pngqD9IiZzTOzpWbWZGZbzexxM7sw6lz5yswWmdlnErT5M+Cgu78UTn/TzFrCn2H7sC+bOd19JbAvzCI5RoVB0mZmXwJ+AHwHKAeqgDuAuRHG6sDMiqLOkIM+B/xv3LxfuXtZzDC8sxU725+p7uOY9r8EPpvKutI7VBgkLWY2DPgH4CZ3f8DdD7l7i7v/2t3/LmxTamY/MLO3w+EHZlYaLpttZo1m9mUz2xEebfxFuOx8M9tmZoUxr/dhM1sZjheY2S1m9oaZ7TazBWY2Mlw22czczD5tZpuA35lZoZl9z8x2mdlGM7s5bFPU/r2Y2Z1hhi1m9u32124/jWNm/2Jme8P1PxCTa6SZ/Xf4/e01s4dill1lZivMbJ+ZPWdmZ3ezP93M/trMNoQ5v2tmnf59mtm7zWyJme0Pv747nH8bcBHw4/C//h93sm4J8F7g94l/yh2y3WRm64B1MT+7r5rZNuC/k/xZv9M+3PQi4NL2dpI7VBgkXRcAA4AHu2nz98D5wLnAOUAN8PWY5eOAYcBE4NPAT8xshLu/ABwieANrNw+4Jxz/a+Bq4BJgArAX+Enca18CnA5cDvwf4ANhjhnhurHuBlqBk4HpwPuB2NMx5wFrgdHAPwN3mpmFy/4XGARMA8YC/wpgZjOAuwj+Ix4F/CewMMGb4IeBmWHGucCn4huEBfBR4Efhdr8PPGpmo9z974E/ADeH//Xf3MlrTAXa3L2xmxyduZpgP5wRTo8DRgKTgBtJ7mcd2x533wK0AKemmEWyzd01aEh5AD4GbEvQ5g3gypjpy4E3w/HZwBGgKGb5DuD8cPzbwF3h+BCCQjEpnF4DXBqz3niCN5giYDLgwEkxy38HfDZm+n1hmyKCU2DHgIExy+uAhnD8k8D6mGWDwnXHha/bBozo5Hv/d+Af4+atBS7pYl85cEXM9F8BT8dkeDYc/wSwOG7d54FPhuOLgM908zN5T/zPDfgm0Azsixka4rK9N2Z6dth+QAo/6w7tY9ptAS6O+vdZQ8dB518lXbuB0WZW5O6tXbSZALwVM/1WOO+dbcStexgoC8fvAZ4zs78ErgGWu3v7tiYBD5pZW8y6xwne5NttjsuxuYtlk4BiYOufDgIoiGuzrX3E3Q+H7coI/gPe4+57OdEk4AYz+3zMvBI6fv/xYl8zfl/Ffi9vxc17i+CoKxl7CQptvAXu/vEkswHsdPej3eSKzx/fvt0QgkIkOUSnkiRdzwNHOfG0TKy3Cd4g21WF8xJy91cJ3lw+QMfTSBC8SX3A3YfHDAM8ODXxziZixrcCFTHTlXHbOgaMjtnWUHeflkTMzcBIMxvexbLb4jIOcvf53WwvNldX+yp+n7a3bf/eE90ueR1gZpZsIWkXv9346UQ/6xNymdkEgmK5NsUskmUqDJIWd98PfIOgX+BqMxtkZsVm9gEz++ew2Xzg62Y2xsxGh+1/kcLL3EPQn3AxcG/M/P8AbjOzSQDh9rv7JNQC4AtmNjF8E/9qzPexFXgS+J6ZDQ07tt9lZpckCheu+zhwh5mNCL//i8PFPwM+Z2bnWWCwmX3QzDr7b73d34XbqQS+APyqkzaPAadY8DHhIjO7nuC8/yPh8u3ASd1kbgF+S9AHk0np/KxnA79z92MZziI9pMIgaXP37wNfIuhk3EnwX/LNwENhk28DS4GVwCvA8nBesubzpzePXTHzfwgsBJ40s4PACwQdo135GcGb/0rgJYI311aC008Af07wn+urBKda7iPoP0jGJwj6N14j6CP5IoC7LyXo9P5xuM31BH0F3XkYWAasIOhgvjO+gbvvBq4CvkxwOu8rwFUx++eHwEfCT0j9qIvX+c8wd6zrreN1DE1mNjZB3ljp/Kw/RlDkJceYux7UI/kl/Ljpf7h7/CmZyJiZA1PdfX0vvd6zwOc9vMitt5nZWcBP3f2CKF5fuqfCIP2emQ0E5hAcNZQD9wMvuPsXo8wVq7cLg0h3dCpJ8oEB3yI4pfMSwcddvxFpIpEcpiMGERHpQEcMIiLSQZ+/wG306NE+efLktNY9dOgQgwcPzmygDMjVXJC72ZQrNcqVmv6Ya9myZbvcfUynC6O+9LqnQ3V1taeroaEh7XWzKVdzueduNuVKjXKlpj/mApZ6F++rOpUkIiIdqDCIiEgHKgwiItJBVguDmVWaWYOZrTGz1Wb2hXD+N8MHoqwIhytj1rnVzNab2Vozuzyb+URE5ETZ/lRSK/Bld18e3jxsmZk9FS77V3f/l9jGZnYGUEvw0JMJwG/N7BR3P46IiPSKrB4xuPtWd18ejh8kuOK0u9v9zgXq3f2Yu28kuPFYTTYziohIR73Wx2Bmkwkem/hiOOtmM1tpZneZ2Yhw3kQ6PhCkkeQfQCIiIhnQK7fEMLMygoeP3+buD5hZObCL4OEd/wiMd/dPmdlPgOfd/RfhencCj7n7/XHbu5HwubHl5eXV9fX1KWfadqiN37xxmHnTBlNSaIlX6EVNTU2UlZUlbhiBXM2mXKlRrtT0x1xz5sxZ5u4zO13Y1QUOmRoIHpv4BPClLpZPBlaF47cCt8YsewK4oLvtp3uB2x/X7fRJX33EH1i+Oa31sylXL6Zxz91sypUa5UpNf8xFVBe4WfBw3DuBNR481KV9fuxDUD4MrArHFwK1ZlZqZlOAqcDibGQ7/6RRlA8y5i+Of5StiEh+y/ankt5D8KSoV8xsRTjva0CdmZ1LcCrpTeCzAO6+2swWEDxJqxW4ybP0iaSCAuPiiiLufX0P63c0cfLY3DtMFBGJQlYLg7s/S3Av/HiPdbPObcBtWQsV48KJxTy4vpX6xZv4+lVn9MZLiojkvLy+8nlYqfH+aeXcv7yRY626VEJEBPK8MADUzqpi7+EWnli9PeooIiI5Ie8Lw4Unj6ZixEDmv7gp6igiIjkh7wtDQYFRV1PF8xt2s3HXoajjiIhELu8LA8BHqysoLDDql+ioQUREhQEYO3QAl542lvuWNtLc2hZ1HBGRSKkwhOrOq2L3oWZ+u0ad0CKS31QYQhdPHcPE4QOZv1ink0Qkv6kwhAoLjOtmVvKHdbvYtPtw1HFERCKjwhDjulkVFBj8aqmOGkQkf6kwxBg/bCBzTh3LgqWNtBxXJ7SI5CcVhjh1NVXsPHiMp9fsiDqKiEgkVBjizD51DOOGDtA1DSKSt1QY4hQVFnDdzAp+//pOGveqE1pE8o8KQyeum1UJwIIleoiPiOQfFYZOVIwYxCWnjGHB0kZa1QktInlGhaELtbOq2HbgKIvW7ow6iohIr1Jh6MKlp49lzJBSXQktInlHhaELxWEndMPaHWzdfyTqOCIivUaFoRvXz6yizWHBksaoo4iI9BoVhm5UjRrERVNH86slmzje5lHHERHpFSoMCdTVVPH2/qM8s06d0CKSH1QYEnjf6eWMLivRM6FFJG+oMCRQUlTAtdUVPP3aDrYfOBp1HBGRrFNhSELtrCqOtzn3LtWV0CLS/6kwJGHK6MG8+12jqF+ymTZ1QotIP6fCkKTamioa9x7h2fW7oo4iIpJVKgxJunxaOSMGFetKaBHp91QYklRaVMhHqit46tXt7Dx4LOo4IiJZo8KQgutnVdHa5ty3TFdCi0j/pcKQgpPHllEzZST1SzapE1pE+q2sFgYzqzSzBjNbY2arzewL4fyRZvaUma0Lv46IWedWM1tvZmvN7PJs5kvHvJoq3tp9mBc27I46iohIVmT7iKEV+LK7nw6cD9xkZmcAtwBPu/tU4OlwmnBZLTANuAK4w8wKs5wxJVecOY5hA4u5R53QItJPZbUwuPtWd18ejh8E1gATgbnA3WGzu4Grw/G5QL27H3P3jcB6oCabGVM1oLiQa2ZM5InV29jdpE5oEel/eq2PwcwmA9OBF4Fyd98KQfEAxobNJgKxlxc3hvNySl1NFS3HnQeWb4k6iohIxpl79jtRzawM+D1wm7s/YGb73H14zPK97j7CzH4CPO/uvwjn3wk85u73x23vRuBGgPLy8ur6+vq0cjU1NVFWVpbWut9+4QhNzc4/XTQQM0trG9nIlW25mk25UqNcqemPuebMmbPM3Wd2utDdszoAxcATwJdi5q0Fxofj44G14fitwK0x7Z4ALuhu+9XV1Z6uhoaGtNe9d+lmn/TVR/z5N3alvY2u9CRXtuVqNuVKjXKlpj/mApZ6F++r2f5UkgF3Amvc/fsxixYCN4TjNwAPx8yvNbNSM5sCTAUWZzNjuj541niGDCiiXp3QItLPZLuP4T3AJ4D3mtmKcLgSuB24zMzWAZeF07j7amAB8CrwG+Amdz+e5YxpGVhSyIenT+SxVdvYe6g56jgiIhlTlM2Nu/uzQFcn4C/tYp3bgNuyFiqDamdV8T/Pv8UDL23h0xdOiTqOiEhG6MrnHjhjwlDOrRxO/eJN7X0iIiJ9ngpDD9XVVLJuRxPL3tobdRQRkYxQYeihq86eQFlpka6EFpF+Q4WhhwaXFjH33Ak8unIr+w+3RB1HRKTHVBgyoK6mimOtbTy0QldCi0jfp8KQAWdOHMZZE4cxX53QItIPqDBkSF1NFa9tO8iKzfuijiIi0iMqDBnyoXMnMKikUM+EFpE+T4UhQ8pKi/jQORP49ctbOXhUndAi0nepMGRQXU0VR1qO8/CKt6OOIiKSNhWGDDq7Yhinjx/KPS+qE1pE+i4VhgwyM+bVVPLq1gO8smV/1HFERNKiwpBhc6dPZEBxAfMXb07cWEQkB6kwZNjQAcVcdfYEFq7YQtOx1qjjiIikTIUhC+pqqjjUfJxfv6xOaBHpe1QYsmBG1XBOLR+ip7uJSJ+kwpAFZkZtTSUvN+5nlTqhRaSPUWHIkg9Pn0hpUQH1S3TUICJ9iwpDlgwfVMIHzxrPwy+9zeFmdUKLSN+RUmEws8FmVpitMP1NbU0VB4+18sjKrVFHERFJWreFwcwKzGyemT1qZjuA14CtZrbazL5rZlN7J2bfNGvyCN41ZrBurCcifUqiI4YG4F3ArcA4d69097HARcALwO1m9vEsZ+yzzIy6mipe2rSP17YdiDqOiEhSEhWG97n7P7r7Sndva5/p7nvc/X53vxb4VXYj9m3XzKigpLCAel0JLSJ9RLeFwd1bAMzs0/HLzOz22DbSuZGDS7jizHE8sLyRoy3Ho44jIpJQsp3PHzGzj7VPmNkdwJjsROp/6mqqOHC0lcdeUSe0iOS+ZAvDNcAnzazOzP4HaHb3E44ipHPnnzSSKaPVCS0ifUOiTyWNNLORwEDgM8BXgAPAP4TzJQlmRu2sSpa8uZd12w9GHUdEpFuJjhiWAUvDrw3AcOCDMfMlSddWV1BcaNQvUSe0iOS2ou4WuvuU3grS340uK+X9Z4zj/uWN/N3lpzKgWNcJikhuSuuWGGY208wmZjpMf1dXU8W+wy08sXpb1FFERLqU7r2SPg88Yma6hiEF737XKKpGDlIntIjktLQKg7vf4O7TCTqku2Rmd5nZDjNbFTPvm2a2xcxWhMOVMctuNbP1ZrbWzC5PJ1suKygwrp9VyQsb9rBhZ1PUcUREOpX23VXN7DR3T/QRm58DV3Qy/1/d/dxweCzc3hlALTAtXOeO/njDvo/OrKCoQJ3QIpK7enLb7ScTNXD3Z4A9SW5vLlDv7sfcfSOwHqjpQb6cNHbIAN53ejn3LWvkWKuuhBaR3GPu3vVCsx91tQi4wd2HJnwBs8nAI+5+Zjj9TeCTBNdDLAW+7O57zezHwAvu/ouw3Z3A4+5+XyfbvBG4EaC8vLy6vr4+UYxONTU1UVZWlta6PbFyZyvfX3aMvzqnlJrxJ34wLKpcycjVbMqVGuVKTX/MNWfOnGXuPrPThe7e5QAcJHgDvqGTYVd368ZsYzKwKma6HCgkOFq5DbgrnP8T4OMx7e4Erk20/erqak9XQ0ND2uv2ROvxNn/3Pz3t8372fKfLo8qVjFzNplypUa7U9MdcwFLv4n212+sYgCXhm/pz8QvC//xT5u7bY7bxM+CRcLIRqIxpWgG8nc5r5LrCguBK6O899Tpv7T7EpFGDo44kIvKORH0MHwFWdLbA07z4zczGx0x+GGj/xNJCoNbMSs1sCjAVWJzOa/QFH51ZSYGhTmgRyTmJbru9x90Pp7txM5sPPA+camaN4e27/9nMXjGzlcAc4G/C11oNLABeBX4D3OTu/bZ3dtywAbz3tHLuXdpIy/G2xCuIiPSSRKeSOjCzf3P3zyfb3t3rOpl9ZzftbyPod8gL886r5LdrtvP0mu1cceb4xCuIiPSCVD+u+p6spMhTl5wylvHDBnCPnu4mIjmkJ9cxSA8VFhjXzazkD+t2snlP2mfsREQyKmFhMLONZrbBzDYCZ7SPm9mGXsjX7103qxIDFizVUYOI5IaEhcHdp7j7SeGnkNa0j7v7Sb2Qr9+bOHwgs08dy6+WbKZVndAikgN0KikH1M6qZMfBY/zutR1RRxERSbkw3JuVFHnuvaeNZeyQUl3TICI5IanCYGalAO7+nZh5euZzhhQVFnD9rEoWrd3Bln1Hoo4jInku2SOGB8ysuH0ivHr5qexEyk/XzazEgQU6ahCRiCVbGB4C7jWzwvBuqU8At2YrVD6qHDmIi6aOYcHSzbR1c8dbEZFsS6owuPvPCI4QHgJ+DXzO3RM+j0FSM6+mkq37j7JyZ7+9E4iI9AHd3hLDzL4UO0lw99MVwPlmdr67fz+L2fLOpaeXM7qslN83tvLFqMOISN5KdMQwJGYoAx4keLJa+zzJoOLCAj46s4KXdx5n2/6jUccRkTzV7RGDu3+rt4JIoHZWJf++6A3uXbqZz186Neo4IpKHuj1iMLOfmtmZXSwbbGafMrOPZSdafpo0ajBnjCqgfslm2trUCS0ivS/RqaQ7gG+Y2Rozu9fM7jCzu8zsD8BzBKeTTngms/TM7Ipituw7wh/W74o6iojkoUSnklYA15lZGTATGA8cIbhn0trsx8tPM8oLGTW4hPkvbuKSU8ZEHUdE8kxSD+px9yZgUVfLzex+d782U6HyXVGBcW11BXc9u5EdB48ydsiAqCOJSB7J1E30dKfVDKudVUlrm3Pfssaoo4hInslUYVAvaYadNKaM808aSf1idUKLSO/SbbdzWF1NFZv2HOa5N3ZHHUVE8kimCoNlaDsS4/Jp4xg+qJj5SzZFHUVE8kimCsNXM7QdiTGguJBrZ1Tw5Opt7Go6FnUcEckTyT6P4T1m9pSZvd7ZM591Q73sqauppOW4c786oUWklyR7xHAn8H3gQmAWwTUNs7IVSv7k5LFDmDV5BPVLNuO6HbeI9IJkC8N+d3/c3Xe4++72IavJ5B11NVVs3HWIFzbsiTqKiOSBRPdKmmFmM4AGM/uumV3QPi+cL73gyrPGM3RAEfMXqxNaRLIv0ZXP34ubnhkz7sB7MxtHOjOguJBrZlRwz4ub2HuomRGDS6KOJCL9WKJ7Jc3prSDSvdqaSn7+3Jvcv7yRz1ykC81FJHt0gVsfcdq4oUyvGs78xZvUCS0iWaXC0IfU1VTxxs5DLH1rb9RRRKQfS1gYzKzAzN7dG2Gke1edPZ4hpUXMf1Gd0CKSPQkLg7u3cWIndFLCh/rsMLNVMfNGhhfLrQu/johZdquZrTeztWZ2eTqv2Z8NKili7vQJPPrKVvYfbok6joj0U8meSnrSzK41s1TvifRz4Iq4ebcAT7v7VODpcBozOwOoBaaF69xhZoUpvl6/V1dTxbHWNh58SVdCi0h2JFsYvgTcCzSb2QEzO2hmBxKt5O7PAPFXZc0F7g7H7waujplf7+7H3H0jsB6oSTJf3pg2YRjnVAxj/mJdCS0i2WHZfnMxs8nAI+5+Zji9z92Hxyzf6+4jzOzHwAvu/otw/p3A4+5+wjOlzexG4EaA8vLy6vr6+rSyNTU1UVZWlta62ZQo16LNLfx8dTNfP28AJ4/o3YOqvrrPoqJcqVGu1PQk15w5c5a5+8xOF7p7UgPwIeBfwuGqFNabDKyKmd4Xt3xv+PUnwMdj5t8JXJto+9XV1Z6uhoaGtNfNpkS5Dh5t8TP+7+P+twtW9E6gGH11n0VFuVKjXKnpSS5gqXfxvprs3VVvB74AvBoOXwjnpWO7mY0Ptzse2BHObwQqY9pVAG+n+Rr9WllpER86dyK/Xvk2B46qE1pEMivZPoYrgcvc/S53v4ugc/jKNF9zIXBDOH4D8HDM/FozKzWzKcBUYHGar9Hv1dVUcrSljYdf2hJ1FBHpZ1K5wG14zPiwZFYws/nA88CpZtZoZp8GbgcuM7N1wGXhNO6+GlhAcETyG+Amdz+eQr68ctbEYUybMJR71AktIhmW6CZ67b4DvGRmDQSP8bwYuDXRSu5e18WiS7tofxtwW5KZ8pqZUVdTxdcfWsXKxv2cUzk86kgi0k8kdeUz0AacDzwQDhe4e3ofBZKMmXvuBAYWF+p23CKSUcle+Xyzu29194Xu/rC7b+uFbJLAkAHF/Nk541n48ts0HWuNOo6I9BPJ9jE8ZWZ/a2aV4S0tRprZyKwmk6TU1VRxuPk4C1foA1wikhnJ9jF8Kvx6U8w8B/RggIidWzmc08YNYf7iTcw7ryrqOCLSDyTbx3CLu0+JG1QUckB7J/QrW/azasv+qOOISD+QbB/DTYnaSXSunj6R0qICdUKLSEaoj6EfGDawmA+ePZ6HV7zNIXVCi0gPJVsYPkVw1PAMsCwclmYrlKRuXk0VTcdaeXTl1qijiEgfl1Tns7tPyXYQ6ZnqSSOYOraMexZv4rpZlYlXEBHpQrdHDGb2lZjxj8Yt+062QknqzIzamipWbN7Hmq0JH5UhItKlRKeSamPG42+BEf9kNonYNdMnUlJUQL06oUWkBxIVButivLNpidiIwSVceeY4HnhpC0eadf9BEUlPosLgXYx3Ni05oLamioNHW3n0FXVCi0h6EhWGc9qf8QycHY63T5/VC/kkRedNGclJowfrdJKIpK3bwuDuhe4+1N2HuHtRON4+XdxbISV57VdCL31rL69vPxh1HBHpg1J5UI/0EdfMmEhxoelKaBFJiwpDPzSqrJTLp43jgeVbONqiTmgRSY0KQz81r6aK/Uda+M0qPTpDRFKjwtBPnX/SKCaNGsQ9Op0kIilSYeinCgqM2llVLN64hzd2NkUdR0T6EBWGfuwj1RUUFZg+uioiKVFh6MfGDCnlsjPKuW9ZI8da1QktIslRYejn6mqq2Hu4hSdXb486ioj0ESoM/dyFJ4+mYsRAXdMgIklTYejngk7oSp57Yzdv7joUdRwR6QNUGPLAR2dWUlhg1C/ZHHUUEekDVBjyQPnQAVx62ljuW7aZ5ta2qOOISI5TYcgTdedVsaupmd+uUSe0iHRPhSFPXDx1DBOHqxNaRBJTYcgThQXGdTMr+cO6XWzeczjqOCKSw1QY8sh1syooMKhfoqMGEelaZIXBzN40s1fMbIWZLQ3njTSzp8xsXfh1RFT5+qPxwwYy59Sx3Lu0kZbj6oQWkc5FfcQwx93PdfeZ4fQtwNPuPhV4OpyWDKqrqWLHwWP87rUdUUcRkRwVdWGINxe4Oxy/G7g6uij90+xTxzBu6AB1QotIl8zdo3lhs43AXsCB/3T3n5rZPncfHtNmr7ufcDrJzG4EbgQoLy+vrq+vTytDU1MTZWVlaa2bTdnO9cC6Zn79RgvfvWQgowem9r9Bvu6zdClXapQrNT3JNWfOnGUxZ2s6cvdIBmBC+HUs8DJwMbAvrs3eRNuprq72dDU0NKS9bjZlO9fmPYd88i2P+PeeXJvyuvm6z9KlXKlRrtT0JBew1Lt4X43sVJK7vx1+3QE8CNQA281sPED4VSfCs6BixCAuOWUMC5ZsplWd0CISJ5LCYGaDzWxI+zjwfmAVsBC4IWx2A/BwFPnyQe2sKrYdOMqitTujjiIiOSaqI4Zy4FkzexlYDDzq7r8BbgcuM7N1wGXhtGTBpaePZcyQUl3TICInKIriRd19A3BOJ/N3A5f2fqL8U1xYwHUzK/j3RW+wdf8Rxg8bGHUkEckRufZxVelF18+sos1hwZLGqKOISA5RYchjVaMGcdHU0SxYupnjbdF8bFlEco8KQ56rq6liy74jPLNOndAiElBhyHPvO72cUYNLmP+iOqFFJKDCkOdKigr4yMwKnn5tBzsOHI06jojkABUGoXZWFcfbnHuXqRNaRFQYBJgyejAXnDSK+iWbaFMntEjeU2EQIHgm9OY9R/jjG7uijiIiEVNhEAAun1bOiEHFuh23iKgwSKC0qJBrZ1Tw5Ort7Dx4LOo4IhIhFQZ5R21NFa1tzv3L1Qktks9UGOQdJ48to2bKSOoXqxNaJJ+pMEgHdTWVvLn7MC9s2B11FBGJiAqDdPCBM8czbGAx85dsjjqKiEREhUE6GFBcyDUzJvLEqm3sblIntEg+UmGQE9TVVNF8vI0Hlm+JOoqIRECFQU5wSvkQqieNYP6STQTPDBeRfKLCIJ2qq6liw85DLN64J+ooItLLVBikUx88azxDBhTpSmiRPKTCIJ0aWFLIh6dP5LFV29h3uDnqOCLSi1QYpEu1s6poblUntEi+UWGQLp0xYSjnVA5n/mJ1QovkExUG6da8mkrW7Whi+aa9UUcRkV6iwiDduursCZSVFnHPi7oSWiRfqDBItwaXFvGhcyfw6Ctvs/9IS9RxRKQXqDBIQvNqqjja0sbDK9QJLZIPVBgkoTMnDuOsicO450V1Qovkg6KoA0jfUFtTyd8/uIqvHjAmvPpHhg8qZvjAYoYPKjlxfFAJwwcWM2JQCUMGFFFQYFHHF5EUqDBIUq6ZXsFbuw+zct1blAwsZs+hZjbsPMS+w80cONra5XpmMGxgN0Wkk4IyfFAxQwcUq6CIRESFQZIysKSQr115OosWbWf27JoOy1qPt3HgaCv7Djez93AL+480s+9wSzg0s+9IML73cHPWCsr2Q23sO9ysgiKSASoM0mNFhQWMHFzCyMElKa2X6YLy1T88pSMUkQzIucJgZlcAPwQKgf9y99sjjiRZksmC8vyylYyfdHJGjlCGDSphRBcFZdigoO9EBUX6s5wqDGZWCPwEuAxoBJaY2UJ3fzXaZJJLOisoBduKmX3hlG7Xy/Ypr/iCMmxgMZs3NfPK8XXdrt/1sp4VnO5W37ihmTW8kda6iVIlim3dbGHDmy1seHYjBQYFBYaZBeMdvhoFBcF4/HJrX27ty2PWeWc5FJ6w7Zi2BSdub8fhNjbvOXzC9goKTtz2CduLydLTn2lvyanCANQA6919A4CZ1QNzARUG6bHITnmtez0L300GvP5a1Ak691qO/rk/09DjTdgJRSyYLmwvJCkWmunDm5k9u+ffWrxcKwwTgdh7LzQC58U3MrMbgRsBysvLWbRoUVov1tTUlPa62ZSruSB3s/VmrgJgZDhQEg7DYlsYUEqbl3Cw6RBlgwd3up2eXBGSaN1El5s0NR2irCyNXAm2mzBXguVNTYcYNHgw7tBG8H24Ow60eTB4+3w6m/Z3xjtso306dr2Y7Xa5zXAbR48epaS0tMM2O2v3p9fyLl879nWCed4xf/t0Z/ljljkwgNbs/N67e84MwEcJ+hXapz8B/Ft361RXV3u6Ghoa0l43m3I1l3vuZlOu1ChXavpjLmCpd/G+mmtXPjcClTHTFcDbEWUREclLuVYYlgBTzWyKmZUAtcDCiDOJiOSVnOpjcPdWM7sZeILg46p3ufvqiGOJiOSVnCoMAO7+GPBY1DlERPJVrp1KEhGRiKkwiIhIByoMIiLSgQqDiIh0YN7Hn8hlZjuBt9JcfTSwK4NxMiVXc0HuZlOu1ChXavpjrknuPqazBX2+MPSEmS1195lR54iXq7kgd7MpV2qUKzX5lkunkkREpAMVBhER6SDfC8NPow7QhVzNBbmbTblSo1ypyatced3HICIiJ8r3IwYREYmjwiAiIh3kRWEwsyvMbK2ZrTezWzpZbmb2o3D5SjObkSO5ZpvZfjNbEQ7f6KVcd5nZDjNb1cXyqPZXoly9vr/MrNLMGsxsjZmtNrMvdNKm1/dXkrmi2F8DzGyxmb0c5vpWJ22i+v1KJltUf5OFZvaSmT3SybLM76+unuDTXwaC23e/AZxE8BDGl4Ez4tpcCTxO8EzG84EXcyTXbOCRCPbZxcAMYFUXy3t9fyWZq9f3FzAemBGODwFez5Hfr2RyRbG/DCgLx4uBF4Hzo95fKWSL6m/yS8A9nb12NvZXPhwx1ADr3X2DuzcD9cDcuDZzgf/xwAvAcDMbnwO5IuHuzwB7umkSxf5KJlevc/et7r48HD8IrCF4dnmsXt9fSebqdeE+aAoni8Mh/hMwUf1+JZOt15lZBfBB4L+6aJLx/ZUPhWEisDlmupET/0CSaRNFLoALwkPbx81sWpYzJSuK/ZWsyPaXmU0GphP8pxkr0v3VTS6IYH+Fp0VWADuAp9w9Z/ZXEtmg9/fZD4CvAG1dLM/4/sqHwmCdzIv/LyCZNpmWzGsuJ7ifyTnAvwEPZTlTsqLYX8mIbH+ZWRlwP/BFdz8Qv7iTVXplfyXIFcn+cvfj7n4uwTPda8zszLgmke2vJLL16j4zs6uAHe6+rLtmnczr0f7Kh8LQCFTGTFcAb6fRptdzufuB9kNbD55sV2xmo7OcKxlR7K+EotpfZlZM8Ob7S3d/oJMmkeyvRLmi/v1y933AIuCKuEWR/351lS2CffYe4ENm9ibB6eb3mtkv4tpkfH/lQ2FYAkw1sylmVgLUAgvj2iwE/jzs3T8f2O/uW6POZWbjzMzC8RqCn9fuLOdKRhT7K6Eo9lf4encCa9z9+1006/X9lUyuiPbXGDMbHo4PBN4HvBbXLJLfr2Sy9fY+c/db3b3C3ScTvEf8zt0/Htcs4/sr5575nGnu3mpmNwNPEHwS6C53X21mnwuX/wfBM6avBNYDh4G/yJFcHwH+0sxagSNArYcfQ8gmM5tP8OmL0WbWCPw/go64yPZXkrmi2F/vAT4BvBKemwb4GlAVkyuK/ZVMrij213jgbjMrJHhTXeDuj0T995hCtkj+JuNle3/plhgiItJBPpxKEhGRFKgwiIhIByoMIiLSgQqDiIh0oMIgIiIdqDCIxDCzpvDrZDObl+Ftfy1u+rlMbl8kU1QYRDo3GUipMISff+9Oh8Lg7u9OMZNIr1BhEOnc7cBFFtxz/2/Cm6t918yWWHDP+8/CO/fnbzCze4BXwnkPmdkyC+7pf2M473ZgYLi9X4bz2o9OLNz2KjN7xcyuj9n2IjO7z8xeM7Nftl91K5JN/f7KZ5E03QL8rbtfBRC+we9391lmVgr80cyeDNvWAGe6+8Zw+lPuvie8rcISM7vf3W8xs5vDG7TFuwY4FzgHGB2u80y4bDowjeDeN38kuKL52Ux/syKxdMQgkpz3E9yPZgXB7atHAVPDZYtjigLAX5vZy8ALBDc3m0r3LgTmh3f23A78HpgVs+1Gd28DVhCc4hLJKh0xiCTHgM+7+xMdZprNBg7FTb8PuMDdD5vZImBAEtvuyrGY8ePob1Z6gY4YRDp3kOCRmO2eILh5WjGAmZ1iZoM7WW8YsDcsCqcRPGqxXUv7+nGeAa4P+zHGEDzCdHFGvguRNOi/D5HOrQRaw1NCPwd+SHAaZ3nYAbwTuLqT9X4DfM7MVgJrCU4ntfspsNLMlrv7x2LmPwhcQPDcbwe+4u7bwsIi0ut0d1UREelAp5JERKQDFQYREelAhUFERDpQYRARkQ5UGEREpAMVBhER6UCFQUREOvj/6P8K9kXkef8AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Took 4 iterations\n",
"1.0\t\n",
"2.0\t\n",
"3.0\t\n",
"4.0\t\n",
"\n",
"Returned error: 0.005742327708394915\n"
]
}
],
"source": [
"result, error, iterations = generalizezMinimalResidual(ex3_A, ex3_B, trace = True, convergencePlot = True)\n",
"print(\"Took {} iterations\".format(iterations))\n",
"print(result)\n",
"print(\"Returned error: {}\".format(error))"
]
},
{
"cell_type": "markdown",
"id": "dbc5f0b5",
"metadata": {},
"source": [
"Now let's run the algorithm on our other system:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "670fa9ad",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEWCAYAAACDoeeyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAt7klEQVR4nO3dd3xV9f348dc7i5BFyCCMBEIIQ0AIEGQpJmr9Yt3WWbVaqzhwfdE6fvVrtbXW4ihad5UOtUZFoW5FTVDAxd57g7JX2Enevz/uCb2kIfOee5J738/H4z64Z37en1x458PnnvM+oqoYY4wJTRFeB2CMMcY9luSNMSaEWZI3xpgQZkneGGNCmCV5Y4wJYZbkjTEmhFmSNyaIRORqEZkS4HNOFZF+gTxnPdvvIyLTvGrf1MySvDlCRH4uItNFpFREfhCRj0TkRK/jClciUiIi19ayz9nAHlWd5Sw/ICKHnc+w8rXTzThVdS6w04nFNDGW5A0AIjIaGAs8DGQAHYFngXM9DOsoIhLldQxN0A3AK1XWvaGqCX6v5OoOrO7nWd+fsd/+rwHX1+dYExyW5A0i0gr4HTBKVd9R1b2qelhV31PVXzv7tBCRsSKy0XmNFZEWzrYCEVkvIneIyGbnfwG/dLYNFpEfRSTSr73zRWSu8z5CRO4RkRUisk1E3hSRFGdbtoioiPxKRNYCX4hIpIg8LiJbRWSViNzs7BNV2RcRedmJYYOIPFTZduVUiYg8JiI7nOPP8IsrRUT+5vRvh4hM9Nt2lojMFpGdIjJNRPrU8PNUEblVRFY6cT4qItX+WxORoSLyvYjscv4c6qz/A3AS8LQzGn+6mmNjgFOAybV/ykfFNkpElgHL/D67u0XkR+Bvdfysj+zvnLoEOLVyP9N0WJI3AEOAWGBCDfv8BhgM5AF9gROA+/y2twVaAR2AXwHPiEhrVf0G2IsvGVX6OfAv5/2twHnAyUB7YAfwTJW2TwaOA/4HuA44w4mjv3Osv38AZUAu0A84HfCf8hgELAHSgDHAyyIizrZXgDigF9AG+DOAiPQHxuEbqaYCLwDv1pLQzgfynRjPBa6puoPzy+wD4CnnvE8AH4hIqqr+BvgKuNkZjd9cTRtdgQpVXV9DHNU5D9/Poaez3BZIAToBI6nbZ+2/P6q6ATgMdK9nLMZtqmqvMH8BlwM/1rLPCuCnfsv/A6x23hcA+4Eov+2bgcHO+4eAcc77RHxJv5OzvAg41e+4dviSRRSQDSiQ47f9C+B6v+XTnH2i8E0zHQRa+m2/DCh23l8NLPfbFucc29ZptwJoXU3fnwN+X2XdEuDkY/ysFBjht3wT8LlfDFOc91cC31U59mvgaud9CXBtDZ/JsKqfG/AAcAjY6fcqrhLbKX7LBc7+sfX4rI/a32+/DcBwr/8+2+vol81xGoBtQJqIRKlq2TH2aQ+s8Vte46w7co4qx+4DEpz3/wKmiciNwAXATFWtPFcnYIKIVPgdW44vYVdaVyWOdcfY1gmIBn74z+CciCr7/Fj5RlX3Ofsl4BuZblfVHfy3TsBVInKL37oYju5/Vf5tVv1Z+fdlTZV1a/D9b6guduD7pVnVm6p6RR1jA9iiqgdqiKtq/FX3r5SI75eKaUJsusaAb/R4gP+e+vC3EV+yq9TRWVcrVV2IL1GcwdFTNeBLOGeoarLfK1Z9//0/cgq/9z8AmX7LWVXOdRBI8ztXkqr2qkOY64AUEUk+xrY/VIkxTlVfr+F8/nEd62dV9WdauW9l32srEbsMEBGp6y+FSlXPW3W5ts/6v+ISkfb4fvEtqWcsxmWW5A2qugu4H988+nkiEici0SJyhoiMcXZ7HbhPRNJFJM3Z/9V6NPMvfPPvw4G3/NY/D/xBRDoBOOev6YqeN4HbRKSDk5Dv9uvHD8CnwOMikuR8qdtFRE6uLTjn2I+AZ0WktdP/4c7mvwI3iMgg8YkXkTNFpLpRdKVfO+fJAm4D3qhmnw+BbuK7dDVKRC7BN0/+vrN9E5BTQ8yHgc/wfWcRSA35rAuAL1T1YIBjMY1kSd4AoKpPAKPxfcG2Bd/o9WZgorPLQ8B0YC4wD5jprKur1/lPItjqt/5J4F3gUxHZA3yD70vBY/krvkQ+F5iFL1GW4ZviAfgFvhHlQnzTGePxzbfXxZX4vg9YjO87hdsBVHU6vi98n3bOuRzf3HpN/g3MAGbj+3L15ao7qOo24CzgDnxTZncBZ/n9fJ4ELnSu9HnqGO284MTt7xI5+jr5UhFpU0u8/hryWV+O7xe2aWJE1R4aYpov5xLI51W16rSHZ0REga6qujxI7U0BblHnhqhgE5HjgRdVdYgX7ZuaWZI3zYqItAQK8Y3mM4C3gW9U9XYv4/IX7CRvTE1susY0NwI8iG/aZBa+SzDv9zQiY5owG8kbY0wIs5G8McaEMFdvhhKR7hx96VgOcL+qjq1u/7S0NM3Ozm5we3v37iU+Pr7BxzdH4dbncOsvWJ/DRWP6PGPGjK2qml7dNleTvKouwVf/AqdI1AZqqI+SnZ3N9OnTG9xeSUkJBQUFDT6+OQq3Podbf8H6HC4a02cRqXrn9BHBnK45FVjhdzu7McYYlwXti1cRGYevZsnTVdaPxKlkl5GRMaCoqKjBbZSWlpKQkFD7jiEk3Pocbv0F63O4aEyfCwsLZ6hqfnXbgpLknbrXG4FeqrrpWPvl5+erTdfUT7j1Odz6C9bncNHI6ZpjJvlgTdecgW8Uf8wEb4wxJvCCleQvw1e7xBhjTBC5nuRFJA74CfCO220ZY4w5musPDVHVffgebWaMMSbIQuLJUBUVyphPltC5oqL2nY0xJoyERJLfsHM/b3y/lkOHD5N93DYG5dh/HIwxBkKkdk1WShwTbhpGYoxwxcvfMn5GfR9eb4wxoSkkkjxAdlo8/ze4JSd0TuHOt+Yw5uPFVFRYhU1jTHgLmSQPEB8t/P2XJ3DZCR15tmQFN702k/2Hyms/0BhjQlRIJXmA6MgIHj6/N/edeRyfLPyRi1/4mk27D3gdljHGeCLkkjyAiHDtSTn89cp8Vmwp5bxnprJg4y6vwzLGmKALySRf6bSeGYy/YSgCXPT810xaaFUVjDHhJaSTPEDP9klMHDWMrm0SGPnKdF78cgX2yENjTLgI+SQP0CYplqKRQzijd1se/nAx974zj0NlduOUMSb0hUWSB2gZE8nTl/XnllNyKfp+HVeN+45d+w57HZYxxrgqbJI8QESEcMfp3Xni4r7MWLOD85+dyqqte70OyxhjXBNWSb7SBf0zee26QezYd4jzn53KNyu3eR2SMca4IiyTPMDA7BQmjhpGanwMV778LW9OX+d1SMYYE3Bhm+QBOqXG885Nwxick8pd4+fyx48WWSkEY0xICeskD9CqZTTjrh7I5YM68sLkldz42gz2HSrzOixjjAmIsE/y4CuF8NB5vfnt2T2ZtHATF7/wNT/uslIIxpjmz5K8Q0T45bDOvHRVPqu27OXcZ6Ywf4OVQjDGNG+W5Ks4pUcG428cSlREBBc9/zWfLPjR65CMMabBLMlX47h2SUwYNZTubRO54dUZPD/ZSiEYY5on15O8iCSLyHgRWSwii0RkiNttBkKbxFiKRg7mzOPb8chHi7n77blWCsEY0+wE4xmvTwIfq+qFIhIDxAWhzYCIjY7kqUv7kZOewFOfL2Pt9n08f8UAkuNivA7NGGPqxNWRvIgkAcOBlwFU9ZCq7nSzzUCLiBBG/6QbYy/JY+aanZz/7DRWbin1OixjjKkTcXOuWUTygBeBhUBfYAZwm6ru9dtnJDASICMjY0BRUVGD2ystLSUhIaExIddo2Y5ynpp1gAqFm/NiOS410rW26srtPjc14dZfsD6Hi8b0ubCwcIaq5le7UVVdewH5QBkwyFl+Evj9sfYfMGCANkZxcXGjjq+Ltdv26mmPl2iXez/Qou/WuN5ebYLR56Yk3Pqran0OF43pMzBdj5FX3f7idT2wXlW/dZbHA/1dbtNVWSlxvH3TUIbmpnH32/P444eLKLdSCMaYJsrVJK+qPwLrRKS7s+pUfFM3zVpSbDTjrsrnF0M68cKXK7nh1RnsPWilEIwxTU8wrpO/BXhNROYCecDDQWjTdVGREfzu3N48cHZPPl+0iYue/5ofdu33OixjjDmK60leVWerar6q9lHV81R1h9ttBtPVwzrz8tUDWbt9H+c+PZW563d6HZIxxhxhd7wGQGH3Nrx941CiIyO4+IWv+Xj+D16HZIwxgCX5gOneNpGJo4ZxXLskbnh1Js+WLLdSCMYYz1mSD6D0xBa8ft1gzunbnjEfL+HX460UgjHGW8EoaxBWYqMjefLSPHLS4xn72X9KIaTEWykEY0zw2UjeBSLC7ad148lL85i9bifnPzuVFVYKwRjjAUvyLjo3rwOvXzeYvQfLOP+ZqUxdvtXrkIwxYcaSvMsGdGrNhJuG0a5VS64a9x2vf7fW65CMMWHEknwQZKXEMf7GIQzLTePed+bx0PsLrRSCMSYoLMkHSWJsNC9flc/VQ7N5acoqrn9lupVCMMa4zpJ8EEVFRvDAOb343bm9KF6yhQuf/5qNO60UgjHGPZbkPfCLIdmMu3og67fv49xnpjJn3U6vQzLGhChL8h45uVs6b980lBZRvlIIH86zUgjGmMCzJO+hbhmJ/HvUMHp3aMVNr83kmWIrhWCMCSxL8h5LTWjBa9cO4ry89jz6yRLueHMOB8vKvQ7LGBMirKxBExAbHcmfL8kjJz2BJyYtZd2OfbxwZb6VQjDGNJqN5JsIEeHWU7vyl8v6MXf9Ls57ZirLN+/xOixjTDNnSb6JObtve4pGDmbfoXLOf3YaXy3b4nVIxphmzJJ8E9SvY2smjhpKh+SWXP2373n1mzVeh2SMaaYsyTdRma3jeOuGIQzvmsZ9E+fzu/esFIIxpv4syTdhibHRvHTVQK4Z1plxU1dx3T+nU2qlEIwx9WBJvomLjBDuP7snD53Xm8lLt3Dhc9PYYKUQjDF15HqSF5HVIjJPRGaLyHS32wtVVwzuxN9/OZANO/dz7tNTmbV2h9chGWOagWCN5AtVNU9V84PUXkg6qWs6E24aSsuYCC598Rven7vR65CMMU2cTdc0M7ltEpl40zD6ZLbi5n/N4pnZB1i6ya6nN8ZUT9yulSIiq4AdgAIvqOqLVbaPBEYCZGRkDCgqKmpwW6WlpSQkJDQi2ubjcIXy7orDfLr6EIfKhRPaRnJubgztE0L793Y4fcaVrM/hoTF9LiwsnHGsmZJgJPn2qrpRRNoAk4BbVPXL6vbNz8/X6dMbPm1fUlJCQUFBg49vjt77tJiFFe34x7TV7D9czjl923PrqV3pkh6a/0DC8TO2PoeHxvRZRI6Z5F0f9qnqRufPzcAE4AS32wwniTHC3SN68NVdhYwcnsOnCzbxkycmM/qN2azautfr8IwxHnM1yYtIvIgkVr4HTgfmu9lmuEpNaMG9ZxzHV3cX8qsTO/Ph/B847YnJ3PHmHNZss2RvTLhyuwplBjBBRCrb+peqfuxym2EtLaEFvzmzJ9cNz+GFySt59Zs1TJy9gQv6deCWU7rSMTXO6xCNMUHkapJX1ZVAXzfbMNVrkxjL/53Vk+uH5/Dc5BW89u1aJszawIUDMhlVmEtWiiV7Y8JBaF+KYWiTFMtvz+7Fl78u5PJBHXln5gYKHyvh3nfm2Z2zxoQBS/Jhom2rWB48tzeT7yrgshM6Mn7GOgoeLea+ifP4YZcle2NClSX5MNOuVUt+f15vSn5dyEX5Wbzx/TpOHlPC/f+ez4+7DngdnjEmwCzJh6kOyS15+PzjKb6zgJ8N6MC/vl3L8EeLeeDdBWzabcnemFBhST7MZbaO448X9KH4zgLOz+vAK9+sYfiYYn733kI277Fkb0xzV68k71z3HulWMMY7WSlx/OnCPnxxx8mc3bc9//h6NcPHFPPQ+wvZsueg1+EZYxqoxiQvIhEi8nMR+UBENgOLgR9EZIGIPCoiXYMTpgmWTqnxPHZRXz4ffTI/Pb4d46auYviYYh7+cBHbSi3ZG9Pc1DaSLwa6APcCbVU1S1XbACcB3wCPiMgVLsdoPJCdFs8TF+fx2eiTGdG7LS99tZKTxhTzyEeL2b73kNfhGWPqqLaboU5T1cNVV6rqduBt4G0RiXYlMtMk5KQn8OdL8hhVmMtTny/jhS9X8MrXq7lqaDbXnZRD6/gYr0M0xtSgxpF8ZYIXkV9V3SYij/jvY0JbbpsEnrqsH5/ePpzCHm14bvIKThpTzGOfLGHnPhvZG9NU1fWL1wtF5PLKBRF5Fkh3JyTTlHXNSOTpn/fn49uGM7xbGk8XL+ekPxXzxKSl7Npvv++NaWrqWrvmAuBdEakAzgC2q+pN7oVlmrrubRN59vIBLPphN09+toynPl/G36au4lcnduaaEzuTFGuzeMY0BbVdXZMiIilAS+Ba4C5gN/A7Z70Jc8e1S+L5Kwfwwa0nMjgnlbGfLePER77gL58vY88BG9kb47XaRvIz8D22T/z+PNN5KZDjanSm2ejVvhV//UU+8zfsYuxnS3l80lJenrqK607K4aqh2SS0cLuqtTGmOjX+y1PVzsEKxISG3h1a8dJVA5m7fidjP1vGo58s4aWvVnLd8ByuGpJNvCV7Y4KqQWUNRCRfRDoEOhgTOvpkJjPu6oFMHDWMPpnJjPl4CSeNKeaFySvYd6jM6/CMCRsNrV1zC/C+iLwRyGBM6MnLSuYf15zA2zcOpVf7JP740WKGjynmpa9Wsv9QudfhGRPyGpTkVfUqVe2H78tYY2o1oFNrXvnVIMbfMIQebZN46INFnDSmmJenrOLAYUv2xrilwVUoRaSHqu4JZDAm9OVnp/DqtYN48/ohdG2TwO/fX8jwMcX8baole2Pc0JhSw58GLAoTdk7onMLrIwdTNHIw2WnxPPjeQk5+tJh/fr2ag2WW7I0JlBovdRCRp461CUiuayNOeeLpwAZVPavO0ZmQNzgnlTdGDubrFdv482dLuf/fC3iuZAU3FeZycX4mLaKssrUxjVHb9Wy/BO4Aqqsxe1k92rkNWAQk1eMYEyZEhKG5aQzpksrU5b5k/38T5/Nc8XJGnZLLRQOyiImy59sY0xC1JfnvgfmqOq3qBhF5oC4NiEgmvpun/gCMrm+AJnyICCd2TWNYbipfLdvKnz9bym8mzOfZ4hXcckouPxuQ6XWIxjQ7oqrH3ugrXXBAVfc1uAGR8cAfgUTgzqrTNSIyEhgJkJGRMaCoqKihTVFaWkpCQkKDj2+OQrnPqsq8reVMWH6YVbsqSG8p/E9mBafmxCMiXocXNKH8GR+L9bl+CgsLZ6hqfnXbarvjdXuDWnSIyFnAZlWdISIFx2jjReBFgPz8fC0oqHa3OikpKaExxzdHod7nQuAWVYqXbObPk5bx6rJddO/ZmSuHZHsdWtCE+mdcHetz4NT3Ga9/qef5hwHniMhqoAg4RURerec5TJgTEU7pkcG/Rw2jd1okD32wiOWb7epdY+qivt9mDavPzqp6r6pmqmo2cCnwhara4wJNg0RECNf2jiG+RRS3vj7bLrU0pg7skgXTrCTHRvCnn/Vh4Q+7eeLTpV6HY0yTV2tJQBFZxX/KDLcTkZXOe1XVOpcaVtUSoKRhYRrzHz/pmcHPB3Xkxa9WcnK3dIbmpnkdkjFNVq0jeVXtrKo5TtnhRZXv65PgjQm0+848js6p8Yx+c449Y9aYGth0jWmW4mKiePLSfmwtPcj/mzCPmi4FNiac1TfJv+VKFMY0wPGZrRh9ejc+nPcj42es9zocY5qkOiV5EWkBoKoP+62zZ7waz10/vAuDOqfwwLsLWLNtr9fhGNPk1HUk/46IRFcuiEg7YJI7IRlTd5ERwhOX5BERIdz+xmzKyiu8DsmYJqWuSX4i8JaIRIpINvAJcK9bQRlTHx2SW/Lw+ccza+1O/vLFcq/DMaZJqdNTlVX1ryISgy/ZZwPXV1e0zBivnN23PcWLN/OXL5YxvFsaAzrZbKIxUMtIXkRGV76AWCALmA0MdtYZ02Q8eG4v2ie35PY3ZrPnwGGvwzGmSahtuibR75UATACW+60zpslIjI1m7CV5bNixnwfeXeh1OMY0CbVVoXwwWIEYEwj52SncXJjLU18sp7BHOmf1ae91SMZ4qrbpmhdFpPcxtsWLyDUicrk7oRnTMLec2pW8rGT+3zvz2Lhzv9fhGOOp2qZrngXuF5FFIvKWiDwrIuNE5CtgGr4pm/GuR2lMPURHRjD2kjzKKpTRb86mvMLuhjXhq7bpmtnAxSKSAOQD7YD9+GrYLHE/PGMaJjstngfO6cVd4+fy169WcsPJXbwOyRhP1PUSylJqqCApIm+r6s8CFZQxgXDRgEyKF2/m8U+XcGJuGr07tPI6JGOCLlAFyqwipWlyRIQ/XnA8qfEtuLVoFvsP2UNGTPgJVJK3SU/TJCXHxfD4xX1ZuWUvD31gl1Wa8GOlhk3IG5abxsjhObz27Vo+W7jJ63CMCapAJXkJ0HmMccUdp3ejZ7sk7np7Lpv3HPA6HGOCJlBJ/u4AnccYV7SIiuTJS/PYe7CMX7811x4yYsJGXevJDxORSSKyVERWisgq51mvAKjqp+6FaExgdM1I5DdnHsfkpVv459drvA7HmKCo0yWUwMvA/wIzgDpfoiAiscCXQAunrfGq+tv6BmlMoFw5uBPFizfzhw8XMaRLKt0yrASTCW11na7ZpaofqepmVd1W+arDcQeBU1S1L5AHjBCRwQ0N1pjGEhHGXNiXxBZR3Pr6LA6W2WWVJrTVVrumv4j0B4pF5FERGVK5zllfI/UpdRajnZdNhhpPpSe2YMyFfVj84x4e/dhu3DahTWr6AkpEims4VlX1lFobEInEN82TCzyjqndX2T4SGAmQkZExoKioqC5xV6u0tJSEhIQGH98chVufA9nffy48yBdry/h1fiy90iIDck43hNtnDNbn+iosLJyhqvnVblTVoLyAZKAY6H2sfQYMGKCNUVxc3Kjjm6Nw63Mg+7vvYJme8lixnvCHSbq99GDAzhto4fYZq1qf6wuYrsfIq0G7GUpVd+KrfzMiWG0aU5OWMZE8eWk/tu89xD3v2GWVJjS5muRFJF1Ekp33LYHTgMVutmlMffTu0Io7T+/OJws28eb0dV6HY0zA1ZrkRSRCRIY28Pzt8H1pOxf4Hpikqu838FzGuOK6k3IY2iWVB99byKqte70Ox5iAqjXJq2oF8HhDTq6qc1W1n6r2UdXeqvq7hpzHGDdFRAiPX9yX6MgIbi+axeHyCq9DMiZg6jpd86mI/ExErEaNCUntWrXk4fOPZ876XTz1+TKvwzEmYOp6x+toIB4oF5H9+AqSqaomuRaZMUF2Zp92FC/J5Jni5Qzvls7A7BSvQzKm0eo0klfVRFWNUNVoVU1yli3Bm5DzwDm9yGwdx+1Fs9l94LDX4RjTaHW+ukZEzhGRx5zXWW4GZYxXElpEMfbSPH7cfYDf/nuB1+EY02h1rUL5CHAbsNB53easMybk9O/YmltOyWXCrA38e/YGr8MxplHqOpL/KfATVR2nquPw3dD0U/fCMsZbNxfm0r9jMvdNnM/6Hfu8DseYBqvPzVDJfu/tsfcmpEVFRjD2kn6owug35lBeYXfDmuaprkn+YWCWiPxdRP6Br+DYw+6FZYz3OqbG8eA5vfhu9Xaen7zC63CMaZBaL6EUkQigAhgMDMR3+eTdqvqjy7EZ47kL+nfgiyWb+fOkpZzUNY0+mcleh2RMvdT1jtebVfUHVX1XVf9tCd6ECxHh4fOOJz2xBbcVzWbfoTKvQzKmXuo6XTNJRO4UkSwRSal8uRqZMU1Eq7honrg4j9Xb9vL79xd6HY4x9VLXO16vcf4c5bdOgZzAhmNM0zSkSyrXD+/C85NXUNC9Df/Tq63XIRlTJ3WqQgnco6qdq7wswZuwMvon3ejdIYl73p7Lpt0HvA7HmDqp65z8qNr2MybUxUT5Lqvcf7icO9+aQ4VdVmmaAZuTN6YectskcN+ZPflq2Vb+Pm211+EYUyubkzemni4f1JGSJZt55OPFDM1NpUdbq9Vnmq66VqGsOh9vc/ImbIkIf/pZH5Jio7nt9dkcOFzudUjGHFONSV5E7vJ7f1GVbXbHqwlbqQktePSiPizZtIc/fWyPLTZNV20j+Uv93t9bZduIAMdiTLNS2L0NVw/N5m9TVzN56RavwzGmWrUleTnG++qWjQk795zRg24ZCdz51hy2lR70Ohxj/kttSV6P8b665f/iXI1TLCKLRGSBiNxW7wiNacJioyMZe0k/du07zD3vzEPVLqs0TUttSb6viOwWkT1AH+d95fLxdTh/GXCHqh6Hr8DZKBHp2ciYjWlSerZP4q4R3Zm0cBOvf7fO63CMOUqNSV5VI/2e6RrlvK9cjq7t5E5Rs5nO+z3AIqBDYEI3pum4ZlhnTsxN4/fvL2TFllKvwzHmCAnWfy9FJBv4Euitqrv91o8ERgJkZGQMKCoqanAbpaWlJCQkNDLS5iXc+tyU+7vjQAX3Td1PessI7hscS1REYL62asp9dov1uX4KCwtnqGp+tRtV1fUXkIDvQSMX1LTfgAEDtDGKi4sbdXxzFG59bur9/WjeD9rp7vf1Tx8tCtg5m3qf3WB9rh9guh4jr9bn8X8NIiLRwNvAa6r6jtvtGeOlEb3bckl+Fs9NXsE3K7d5HY4x7iZ5ERHgZWCRqj7hZlvGNBX3n92TTilxjH5jNrv2H/Y6HBPm3B7JDwOuBE4RkdnO66cut2mMp+JbRDH20n5s2nOQ+ybOt8sqjafqWqCsQVR1CnbTlAlDeVnJ/O9pXXns06Wc0iOd8/tleh2SCVOuz8kbE65uLMhlYHZr7p+4gHXb93kdjglTluSNcUlkhPDExXkA/O8bsykrr/A2IBOWLMkb46KslDh+f15vpq/ZwXMlK7wOx4QhS/LGuOy8fh04p297xn6+jFlrd3gdjgkzluSNCYLfn9ebtkmx3P7GbPYeLPM6HBNGLMkbEwStWkbzxMV9Wbt9Hw++t8DrcEwYsSRvTJAMyknlpoIuvDl9PR/N+8HrcEyYsCRvTBDdflo3+mS24p535vHjrgNeh2PCgCV5Y4IoOjKCsZfkcaisgjvemk1Fhd0Na9xlSd6YIMtJT+D+s3sydfk2xk1d5XU4JsRZkjfGA5cOzOL0nhmM+XgJCzfurv0AYxrIkrwxHhARHvlZH5LjormtaBYHDpd7HZIJUZbkjfFISnwMj13Ul2WbS/njh4u8DseEKEvyxnhoeLd0rhnWmX98vYbixZu9DseEIEvyxnjsrhHd6dE2kV+Pn8PW0oNeh2NCjCV5YzwWGx3J2Evz2H2gjLvHz7WHjJiAsiRvTBPQo20S94zoweeLN/Pqt2u9DseEEEvyxjQRVw/NZni3dP7wwUKWb97jdTgmRFiSN6aJiIgQHruwD3ExUdxWNJtDZfaQEdN4luSNaULaJMXyyAXHs2Djbh6ftMTrcEwIsCRvTBNzeq+2XHZCR178ciXTVmz1OhzTzLma5EVknIhsFpH5brZjTKj5v7OOo3NqPHe8OYe9h+1qG9Nwbo/k/w6McLkNY0JOXEwUT17ajy17DvLSvIOs2FJql1aaBoly8+Sq+qWIZLvZhjGh6vjMVtw9ogd/+HARpz4+mXatYhmWm8aJuWkMzU2lTWKs1yGaZkDcHh04Sf59Ve19jO0jgZEAGRkZA4qKihrcVmlpKQkJCQ0+vjkKtz6HW38BVm0pZfX+FizYVs6i7eXsPexbn5kg9EyNpGdqJN1TImkZJd4GGkDh+Dk3ps+FhYUzVDW/um2eJ3l/+fn5On369Aa3VVJSQkFBQYOPb47Crc/h1l84us/lFcrCjbuZsnwrU5dv5bvV2zlUVkFUhJCXlewb6XdNIy8rmejI5ntdRbh/zvUlIsdM8q5O1xhjAisyQjg+sxXHZ7bixoIuHDhczow1O5iyfCvTlm/lqS+W8eTny4iPiWRQTuqR6Z1uGQmIhM5I39SdJXljmrHY6EiG5aYxLDcNgJ37DvHNym3OSH8bXziVLdMTWzCsS+qRfdsnt/QybBNEriZ5EXkdKADSRGQ98FtVfdnNNo0JZ8lxMYzo3Y4RvdsBsH7HPqYt9yX9r5ZtZeLsjQDkpMdzopPwB+ek0qpltJdhGxe5fXXNZW6e3xhTs8zWcVw8MI6LB2ZRUaEs2bSHqcu3MmX5Vt6avp5/fr2GCIE+mclHkn7/Tsm0iIr0OnQTIDZdY0yYiIgQjmuXxHHtkrj2pBwOlVUwa+0Opi7fytQV23hu8gqeLl5ObHQEA7NTjiT9nu2SiIiw+fzmypK8MWEqJiqCQTmpDMpJZTSw58Bhvl25/ciVO3/8aDEAreOiGep8gXtibhpZKXHeBm7qxZK8MQaAxNhoTuuZwWk9MwDYtPvAkamdqcu38sHcHwDomBJ35KqdIV1SSYmP8TJsUwtL8saYamUkxXJB/0wu6J+JqrJiSylTlm1lyvJtvDdnI69/txYR6NU+6UjSH5idQmy0zec3JZbkjTG1EhFy2ySS2yaRq4d1pqy8gjnrdzHNGemPm7KKFyavJCYqgvxOrY9cqnl8h1ZE2ny+pyzJG2PqLSoyggGdWjOgU2tuObUr+w6V8d2q7c70zjYe/WQJj36yhKTYKIZ0ST3yJW7ntHi7KSvILMkbYxotLiaKgu5tKOjeBoCtpQeZtmIbU5f5RvqfLNgEQPvKImtd0xjaJY30xBZehh0WLMkbYwIuLaEF5/Rtzzl926OqrNm2z1d6YcVWPl24ibdmrAegR9vEI/P5J3ROIb6FpaRAs5+oMcZVIkJ2WjzZafFcMbjTfxVZe+WbNbw8ZRVREUL/jr75/Njd5Qw+XG5f4gaAJXljTFDVVGRt6vKtjP18Karw6PRP6Nk+if4dW9OvYzL9O7Yms3VLm9OvJ0vyxhhPVVdk7W/vf8nhpExmrt3Bm9PX8fdpqwHfNFC/jslHkn6fzFbExVgaq4n9dIwxTUpyXAz92kRRUNADgLLyCpZs2sOstTuZuXYHs9fuZNJC3xe5kRFCj7aJR432O6XG2WjfjyV5Y0yTFhUZQa/2rejVvhVXDO4EwI69h5i9zpf0Z67dwYRZG3jlmzUApMTH0C8rmf6dWtMvK5k+WckkhPEXuuHbc2NMs9U6PobCHm0o7OG7ZLO8Qlm22Rntr9nBrHU7+dyppR8h0C0j8UjS79+pNTlhdL2+JXljTLPnm7ZJokfbJC47oSMAu/YdZta6HUemed6bs5F/fbsWgFYto49M7/TrmEzfrGSSYkOzpr4leWNMSGoVF33UDVoVFb76O5VJf9banfx5qe9KHhHo2ibhqLn9LukJIVFi2ZK8MSYsREQIXTMS6ZqRyMUDswDYfeAwc9ftOjK3/9H8Hyn6fh0AibFR5GX9Z7TfL6s1reKa32jfkrwxJmwlxUZzYldfmQUAVWXl1r1Hjfb/8sUyKtS3f5f0ePp3bO2b3++YTNc2iU2+AJsleWOMcYgIXdIT6JKewIUDMgEoPVjG3PU7j3yp+/nizUfKMiS0iKJvVqujRvutm1h9fUvyxhhTg4QWUQzt4iuoBhypxTNr3Q5mrtnJrHU7eLZkBeXOcD8nLZ48vy91u2ckEhUZ4Vn8rid5ERkBPAlEAi+p6iNut2mMMW7xr8Vzfj/faH/foTLmrt91ZJrny6VbeGfmBgDiYiLpk1k52m9N/47JpCYEr/qmq0leRCKBZ4CfAOuB70XkXVVd6Ga7xhgTTHExUQzOSWVwTirgG+2v37H/yLz+rLU7ePHLlZQ5o/1OqXF+N2y1pke7RNdic3skfwKwXFVXAohIEXAuYEneGBOyRISslDiyUuI4N68DAAcOlzN/g3Mlz5qdTFuxjYmzNwIQGx1BrxShoMCFWFQ18GetPLnIhcAIVb3WWb4SGKSqN/vtMxIYCZCRkTGgqKiowe2VlpaSkJDQuKCbmXDrc7j1F6zPoUpV2X5AWbGzguU7y6koO8wVxzesz4WFhTNUNb+6bW6P5Ku7tuio3yqq+iLwIkB+fr4WNOJXWUlJCY05vjkKtz6HW3/B+hwu3Oqz21/5rgey/JYzgY0ut2mMMcbhdpL/HugqIp1FJAa4FHjX5TaNMcY4XJ2uUdUyEbkZ+ATfJZTjVHWBm20aY4z5D9evk1fVD4EP3W7HGGPMf/PuNixjjDGusyRvjDEhzJK8McaEMEvyxhgTwly947W+RGQLsKYRp0gDtgYonOYi3Pocbv0F63O4aEyfO6lqenUbmlSSbywRmX6sW3tDVbj1Odz6C9bncOFWn226xhhjQpgleWOMCWGhluRf9DoAD4Rbn8Otv2B9Dheu9Dmk5uSNMcYcLdRG8sYYY/xYkjfGmBAWEkleREaIyBIRWS4i93gdTzCIyDgR2Swi872OJRhEJEtEikVkkYgsEJHbvI7JbSISKyLficgcp88Peh1TMIhIpIjMEpH3vY4lWERktYjME5HZIjI9oOdu7nPyzsPCl+L3sHDgslB/WLiIDAdKgX+qam+v43GbiLQD2qnqTBFJBGYA54Xy5ywiAsSraqmIRANTgNtU9RuPQ3OViIwG8oEkVT3L63iCQURWA/mqGvAbwEJhJH/kYeGqegiofFh4SFPVL4HtXscRLKr6g6rOdN7vARYBHbyNyl3qU+osRjuv5j0qq4WIZAJnAi95HUuoCIUk3wFY57e8nhD/xx/uRCQb6Ad863EornOmLmYDm4FJqhrqfR4L3AVUeBxHsCnwqYjMEJGRgTxxKCT5Wh8WbkKHiCQAbwO3q+pur+Nxm6qWq2oevucjnyAiITs1JyJnAZtVdYbXsXhgmKr2B84ARjnTsQERCkneHhYeJpx56beB11T1Ha/jCSZV3QmUACO8jcRVw4BznPnpIuAUEXnV25CCQ1U3On9uBibgm4YOiFBI8vaw8DDgfAn5MrBIVZ/wOp5gEJF0EUl23rcETgMWexqUi1T1XlXNVNVsfP+Ov1DVKzwOy3UiEu9cTICIxAOnAwG7aq7ZJ3lVLQMqHxa+CHgzHB4WLiKvA18D3UVkvYj8yuuYXDYMuBLf6G628/qp10G5rB1QLCJz8Q1mJqlq2FxWGEYygCkiMgf4DvhAVT8O1Mmb/SWUxhhjjq3Zj+SNMcYcmyV5Y4wJYZbkjTEmhFmSN8aYEGZJ3hhjQpgleROyRKTU+TNbRH4e4HP/vyrL0wJ5fmMCxZK8CQfZQL2SvFPdtCZHJXlVHVrPmIwJCkvyJhw8Apzk3ED1v07Rr0dF5HsRmSsi1wOISIFTs/5fwDxn3USnaNSCysJRIvII0NI532vOusr/NYhz7vlOffBL/M5dIiLjRWSxiLzm3MVrjKuivA7AmCC4B7izsja5k6x3qepAEWkBTBWRT519TwB6q+oqZ/kaVd3ulBX4XkTeVtV7RORmp3BYVRcAeUBfIM055ktnWz+gF77aSlPx3cU7JdCdNcafjeRNODod+IVTwvdbIBXo6mz7zi/BA9zq3G7+Db5CeF2p2YnA6071yE3AZGCg37nXq2oFMBvfNJIxrrKRvAlHAtyiqp8ctVKkANhbZfk0YIiq7hOREiC2Duc+loN+78uxf38mCGwkb8LBHiDRb/kT4EandDEi0s2p/ldVK2CHk+B7AIP9th2uPL6KL4FLnHn/dGA4vqJTxnjCRhImHMwFypxpl78DT+KbKpnpfPm5BTivmuM+Bm5wqkAuwTdlU+lFYK6IzFTVy/3WTwCGAHPwPbzmLlX90fklYUzQWRVKY4wJYTZdY4wxIcySvDHGhDBL8sYYE8IsyRtjTAizJG+MMSHMkrwxxoQwS/LGGBPC/j+Pv3opVUheOQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Took 5 iterations\n",
"-0.7512\t\n",
"1.0016\t\n",
"3.0046\t\n",
"-0.0002\t\n",
"-1.252\t\n",
"\n",
"Returned error: 0.0004358158644427713\n"
]
}
],
"source": [
"result, error, iterations = generalizezMinimalResidual(ex1_A, ex1_B, convergencePlot = True)\n",
"print(\"Took {} iterations\".format(iterations))\n",
"print(result)\n",
"print(\"Returned error: {}\".format(error))"
]
}
],
"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.8.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment