{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 10. Implementing QR Factorization" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We used QR factorization in computing eigenvalues and to compute least squares regression. It is an important building block in numerical linear algebra.\n", "\n", "\"One algorithm in numerical linear algebra is more important than all the others: QR factorization.\" --Trefethen, page 48\n", "\n", "Recall that for any matrix $A$, $A = QR$ where $Q$ is orthogonal and $R$ is upper-triangular." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Reminder**: The *QR algorithm*, which we looked at in the last lesson, uses the *QR decomposition*, but don't confuse the two." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### In Numpy ###" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "\n", "np.set_printoptions(suppress=True, precision=4)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "n = 5\n", "A = np.random.rand(n,n)\n", "npQ, npR = np.linalg.qr(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check that Q is orthogonal:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(True, True)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(np.eye(n), npQ @ npQ.T), np.allclose(np.eye(n), npQ.T @ npQ)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check that R is triangular" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-0.8524, -0.7872, -1.1163, -1.2248, -0.7587],\n", " [ 0. , -0.9363, -0.2958, -0.7666, -0.632 ],\n", " [ 0. , 0. , 0.4645, -0.1744, -0.3542],\n", " [ 0. , 0. , 0. , 0.4328, -0.2567],\n", " [ 0. , 0. , 0. , 0. , 0.1111]])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "npR" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Linear Algebra Review: Projections" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When vector $\\mathbf{b}$ is projected onto a line $\\mathbf{a}$, its projection $\\mathbf{p}$ is the part of $\\mathbf{b}$ along that line $\\mathbf{a}$.\n", "\n", "Let's look at interactive graphic (3.4) for [section 3.2.2: Projections](http://immersivemath.com/ila/ch03_dotproduct/ch03.html) of the [Immersive Linear Algebra online book](http://immersivemath.com/ila/index.html).\n", "\n", "\"projection\"\n", "(source: [Immersive Math](http://immersivemath.com/ila/ch03_dotproduct/ch03.html))\n", "\n", "And here is what it looks like to project a vector onto a plane:\n", "\n", "\"projection\"\n", "(source: [The Linear Algebra View of Least-Squares Regression](https://medium.com/@andrew.chamberlain/the-linear-algebra-view-of-least-squares-regression-f67044b7f39b))\n", "\n", "When vector $\\mathbf{b}$ is projected onto a line $\\mathbf{a}$, its projection $\\mathbf{p}$ is the part of $\\mathbf{b}$ along that line $\\mathbf{a}$. So $\\mathbf{p}$ is some multiple of $\\mathbf{a}$. Let $\\mathbf{p} = \\hat{x}\\mathbf{a}$ where $\\hat{x}$ is a scalar." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Orthogonality" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**The key to projection is orthogonality:** The line *from* $\\mathbf{b}$ to $\\mathbf{p}$ (which can be written $\\mathbf{b} - \\hat{x}\\mathbf{a}$) is perpendicular to $\\mathbf{a}$.\n", "\n", "This means that $$ \\mathbf{a} \\cdot (\\mathbf{b} - \\hat{x}\\mathbf{a}) = 0 $$\n", "\n", "and so $$\\hat{x} = \\frac{\\mathbf{a} \\cdot \\mathbf{b}}{\\mathbf{a} \\cdot \\mathbf{a}} $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gram-Schmidt ##" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Classical Gram-Schmidt (unstable) ###" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For each $j$, calculate a single projection $$v_j = P_ja_j$$ where $P_j$ projects onto the space orthogonal to the span of $q_1,\\ldots,q_{j-1}$." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def cgs(A):\n", " m, n = A.shape\n", " Q = np.zeros([m,n], dtype=np.float64)\n", " R = np.zeros([n,n], dtype=np.float64)\n", " for j in range(n):\n", " v = A[:,j]\n", " for i in range(j):\n", " R[i,j] = np.dot(Q[:,i], A[:,j])\n", " v = v - (R[i,j] * Q[:,i])\n", " R[j,j] = np.linalg.norm(v)\n", " Q[:, j] = v / R[j,j]\n", " return Q, R" ] }, { "cell_type": "code", "execution_count": 190, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Q, R = cgs(A)" ] }, { "cell_type": "code", "execution_count": 191, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 191, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(A, Q @ R)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check if Q is unitary:" ] }, { "cell_type": "code", "execution_count": 109, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 109, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(np.eye(len(Q)), Q.dot(Q.T))" ] }, { "cell_type": "code", "execution_count": 115, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 115, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(npQ, -Q)" ] }, { "cell_type": "code", "execution_count": 192, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0.02771, 0.02006, -0.0164 , ..., 0.00351, 0.00198, 0.00639],\n", " [ 0. , 0.10006, -0.00501, ..., 0.07689, -0.0379 , -0.03095],\n", " [ 0. , 0. , 0.01229, ..., 0.01635, 0.02988, 0.01442],\n", " ..., \n", " [ 0. , 0. , 0. , ..., 0. , -0. , -0. ],\n", " [ 0. , 0. , 0. , ..., 0. , 0. , -0. ],\n", " [ 0. , 0. , 0. , ..., 0. , 0. , 0. ]])" ] }, "execution_count": 192, "metadata": {}, "output_type": "execute_result" } ], "source": [ "R" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Gram-Schmidt should remind you a bit of the Arnoldi Iteration (used to transform a matrix to Hessenberg form) since it is also a structured orthogonalization." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Modified Gram-Schmidt ###" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Classical (unstable) Gram-Schmidt: for each $j$, calculate a single projection $$v_j = P_ja_j$$ where $P_j$ projects onto the space orthogonal to the span of $q_1,\\ldots,q_{j-1}$.\n", "\n", "Modified Gram-Schmidt: for each $j$, calculate $j-1$ projections $$P_j = P_{\\perp q_{j-1}\\cdots\\perp q_{2}\\perp q_{1}}$$" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "n = 3\n", "A = np.random.rand(n,n).astype(np.float64)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def cgs(A):\n", " m, n = A.shape\n", " Q = np.zeros([m,n], dtype=np.float64)\n", " R = np.zeros([n,n], dtype=np.float64)\n", " for j in range(n):\n", " v = A[:,j]\n", " for i in range(j):\n", " R[i,j] = np.dot(Q[:,i], A[:,j])\n", " v = v - (R[i,j] * Q[:,i])\n", " R[j,j] = np.linalg.norm(v)\n", " Q[:, j] = v / R[j,j]\n", " return Q, R" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def mgs(A):\n", " V = A.copy()\n", " m, n = A.shape\n", " Q = np.zeros([m,n], dtype=np.float64)\n", " R = np.zeros([n,n], dtype=np.float64)\n", " for i in range(n):\n", " R[i,i] = np.linalg.norm(V[:,i])\n", " Q[:,i] = V[:,i] / R[i,i]\n", " for j in range(i, n):\n", " R[i,j] = np.dot(Q[:,i],V[:,j])\n", " V[:,j] = V[:,j] - R[i,j]*Q[:,i]\n", " return Q, R" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Q, R = mgs(A)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(np.eye(len(Q)), Q.dot(Q.T.conj()))" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(A, np.matmul(Q,R))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Householder" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Intro" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{array}{ l | l | c }\n", "\\hline\n", "Gram-Schmidt & Triangular\\, Orthogonalization & A R_1 R_2 \\cdots R_n = Q \\\\\n", "Householder & Orthogonal\\, Triangularization & Q_n \\cdots Q_2 Q_1 A = R \\\\\n", "\\hline\n", "\\end{array}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Householder reflections lead to a more nearly orthogonal matrix Q with rounding errors\n", "\n", "Gram-Schmidt can be stopped part-way, leaving a reduced QR of 1st n columns of A" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Initialization" ] }, { "cell_type": "code", "execution_count": 113, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "n = 4\n", "A = np.random.rand(n,n).astype(np.float64)\n", "\n", "Q = np.zeros([n,n], dtype=np.float64)\n", "R = np.zeros([n,n], dtype=np.float64)" ] }, { "cell_type": "code", "execution_count": 114, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0.5435, 0.6379, 0.4011, 0.5773],\n", " [ 0.0054, 0.8049, 0.6804, 0.0821],\n", " [ 0.2832, 0.2416, 0.8656, 0.8099],\n", " [ 0.1139, 0.9621, 0.7623, 0.5648]])" ] }, "execution_count": 114, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A" ] }, { "cell_type": "code", "execution_count": 115, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from scipy.linalg import block_diag" ] }, { "cell_type": "code", "execution_count": 116, "metadata": { "collapsed": true }, "outputs": [], "source": [ "np.set_printoptions(5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Algorithm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I added more computations and more info than needed, since it's illustrative of how the algorithm works. This version returns the *Householder Reflectors* as well." ] }, { "cell_type": "code", "execution_count": 117, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def householder_lots(A):\n", " m, n = A.shape\n", " R = np.copy(A)\n", " V = []\n", " Fs = []\n", " for k in range(n):\n", " v = np.copy(R[k:,k])\n", " v = np.reshape(v, (n-k, 1))\n", " v[0] += np.sign(v[0]) * np.linalg.norm(v)\n", " v /= np.linalg.norm(v)\n", " R[k:,k:] = R[k:,k:] - 2*np.matmul(v, np.matmul(v.T, R[k:,k:]))\n", " V.append(v)\n", " F = np.eye(n-k) - 2 * np.matmul(v, v.T)/np.matmul(v.T, v)\n", " Fs.append(F)\n", " return R, V, Fs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check that R is upper triangular:" ] }, { "cell_type": "code", "execution_count": 121, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-0.62337, -0.84873, -0.88817, -0.97516],\n", " [ 0. , -1.14818, -0.86417, -0.30109],\n", " [ 0. , 0. , -0.64691, -0.45234],\n", " [-0. , 0. , 0. , -0.26191]])" ] }, "execution_count": 121, "metadata": {}, "output_type": "execute_result" } ], "source": [ "R" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a check, we will calculate $Q^T$ and $R$ using the block matrices $F$. The matrices $F$ are *the householder reflectors*.\n", "\n", "Note that this is not a computationally efficient way of working with $Q$. In most cases, you do not actually need $Q$. For instance, if you are using QR to solve least squares, you just need $Q^*b$.\n", "\n", "- See page 74 of Trefethen for techniques to implicitly calculate the product $Q^*b$ or $Qx$.\n", "- See [these lecture notes](http://www.cs.cornell.edu/~bindel/class/cs6210-f09/lec18.pdf) for a different implementation of Householder that calculates Q simultaneously as part of R. " ] }, { "cell_type": "code", "execution_count": 175, "metadata": { "collapsed": true }, "outputs": [], "source": [ "QT = np.matmul(block_diag(np.eye(3), F[3]), \n", " np.matmul(block_diag(np.eye(2), F[2]), \n", " np.matmul(block_diag(np.eye(1), F[1]), F[0])))" ] }, { "cell_type": "code", "execution_count": 178, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-0.69502, 0.10379, -0.71146],\n", " [ 0.10379, 0.99364, 0.04356],\n", " [-0.71146, 0.04356, 0.70138]])" ] }, "execution_count": 178, "metadata": {}, "output_type": "execute_result" } ], "source": [ "F[1]" ] }, { "cell_type": "code", "execution_count": 177, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 1. , 0. , 0. , 0. ],\n", " [ 0. , -0.69502, 0.10379, -0.71146],\n", " [ 0. , 0.10379, 0.99364, 0.04356],\n", " [ 0. , -0.71146, 0.04356, 0.70138]])" ] }, "execution_count": 177, "metadata": {}, "output_type": "execute_result" } ], "source": [ "block_diag(np.eye(1), F[1])" ] }, { "cell_type": "code", "execution_count": 193, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 1. , 0. , 0. , 0. ],\n", " [ 0. , 1. , 0. , 0. ],\n", " [ 0. , 0. , -0.99989, 0.01452],\n", " [ 0. , 0. , 0.01452, 0.99989]])" ] }, "execution_count": 193, "metadata": {}, "output_type": "execute_result" } ], "source": [ "block_diag(np.eye(2), F[2])" ] }, { "cell_type": "code", "execution_count": 194, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 1., 0., 0., 0.],\n", " [ 0., 1., 0., 0.],\n", " [ 0., 0., 1., 0.],\n", " [ 0., 0., 0., -1.]])" ] }, "execution_count": 194, "metadata": {}, "output_type": "execute_result" } ], "source": [ "block_diag(np.eye(3), F[3])" ] }, { "cell_type": "code", "execution_count": 176, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-0.87185, -0.00861, -0.45431, -0.18279],\n", " [ 0.08888, -0.69462, 0.12536, -0.70278],\n", " [-0.46028, 0.10167, 0.88193, -0.00138],\n", " [-0.14187, -0.71211, 0.00913, 0.68753]])" ] }, "execution_count": 176, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.matmul(block_diag(np.eye(1), F[1]), F[0])" ] }, { "cell_type": "code", "execution_count": 123, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-0.87185, -0.00861, -0.45431, -0.18279],\n", " [ 0.08888, -0.69462, 0.12536, -0.70278],\n", " [ 0.45817, -0.112 , -0.88171, 0.01136],\n", " [ 0.14854, 0.71056, -0.02193, -0.68743]])" ] }, "execution_count": 123, "metadata": {}, "output_type": "execute_result" } ], "source": [ "QT" ] }, { "cell_type": "code", "execution_count": 124, "metadata": { "collapsed": true }, "outputs": [], "source": [ "R2 = np.matmul(block_diag(np.eye(3), F[3]), \n", " np.matmul(block_diag(np.eye(2), F[2]),\n", " np.matmul(block_diag(np.eye(1), F[1]), \n", " np.matmul(F[0], A))))" ] }, { "cell_type": "code", "execution_count": 125, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 125, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(A, np.matmul(np.transpose(QT), R2))" ] }, { "cell_type": "code", "execution_count": 126, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 126, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(R, R2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's a concise version of Householder (although I do create a new R, instead of overwriting A and computing it in place)." ] }, { "cell_type": "code", "execution_count": 179, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def householder(A):\n", " m, n = A.shape\n", " R = np.copy(A)\n", " Q = np.eye(m)\n", " V = []\n", " for k in range(n):\n", " v = np.copy(R[k:,k])\n", " v = np.reshape(v, (n-k, 1))\n", " v[0] += np.sign(v[0]) * np.linalg.norm(v)\n", " v /= np.linalg.norm(v)\n", " R[k:,k:] = R[k:,k:] - 2 * v @ v.T @ R[k:,k:]\n", " V.append(v)\n", " return R, V" ] }, { "cell_type": "code", "execution_count": 157, "metadata": { "collapsed": true }, "outputs": [], "source": [ "RH, VH = householder(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check that R is diagonal:" ] }, { "cell_type": "code", "execution_count": 129, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-0.62337, -0.84873, -0.88817, -0.97516],\n", " [-0. , -1.14818, -0.86417, -0.30109],\n", " [-0. , -0. , -0.64691, -0.45234],\n", " [-0. , 0. , 0. , -0.26191]])" ] }, "execution_count": 129, "metadata": {}, "output_type": "execute_result" } ], "source": [ "RH" ] }, { "cell_type": "code", "execution_count": 130, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[array([[ 0.96743],\n", " [ 0.00445],\n", " [ 0.2348 ],\n", " [ 0.09447]]), array([[ 0.9206 ],\n", " [-0.05637],\n", " [ 0.38641]]), array([[ 0.99997],\n", " [-0.00726]]), array([[ 1.]])]" ] }, "execution_count": 130, "metadata": {}, "output_type": "execute_result" } ], "source": [ "VH" ] }, { "cell_type": "code", "execution_count": 131, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 131, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(R, RH)" ] }, { "cell_type": "code", "execution_count": 132, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def implicit_Qx(V,x):\n", " n = len(x)\n", " for k in range(n-1,-1,-1):\n", " x[k:n] -= 2*np.matmul(v[-k], np.matmul(v[-k], x[k:n]))" ] }, { "cell_type": "code", "execution_count": 133, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "array([[ 0.54348, 0.63791, 0.40114, 0.57728],\n", " [ 0.00537, 0.80485, 0.68037, 0.0821 ],\n", " [ 0.2832 , 0.24164, 0.86556, 0.80986],\n", " [ 0.11395, 0.96205, 0.76232, 0.56475]])" ] }, "execution_count": 133, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Both classical and modified Gram-Schmidt require $2mn^2$ flops." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gotchas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some things to be careful about:\n", "- when you've copied values vs. when you have two variables pointing to the same memory location\n", "- the difference between a vector of length n and a 1 x n matrix (np.matmul deals with them differently)" ] }, { "cell_type": "markdown", "metadata": { "heading_collapsed": true }, "source": [ "## Analogy" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "| | $A = QR$ | $A = QHQ^*$ |\n", "|------------------------------|--------------|-------------|\n", "| orthogonal structuring | Householder | Householder |\n", "| structured orthogonalization | Gram-Schmidt | Arnoldi |\n", "\n", "**Gram-Schmidt and Arnoldi**: succession of triangular operations, can be stopped part way and first $n$ columns are correct\n", "\n", "**Householder**: succession of orthogoal operations. Leads to more nearly orthogonal A in presence of rounding errors.\n", "\n", "Note that to compute a Hessenberg reduction $A = QHQ^*$, Householder reflectors are applied to two sides of A, rather than just one." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Examples" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following examples are taken from Lecture 9 of Trefethen and Bau, although translated from MATLAB into Python" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Ex: Classical vs Modified Gram-Schmidt ###" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This example is experiment 2 from section 9 of Trefethen. We want to construct a square matrix A with random singular vectors and widely varying singular values spaced by factors of 2 between $2^{-1}$ and $2^{-(n+1)}$" ] }, { "cell_type": "code", "execution_count": 68, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "from matplotlib import rcParams\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 183, "metadata": { "collapsed": true }, "outputs": [], "source": [ "n = 100\n", "U, X = np.linalg.qr(np.random.randn(n,n)) # set U to a random orthogonal matrix\n", "V, X = np.linalg.qr(np.random.randn(n,n)) # set V to a random orthogonal matrix\n", "S = np.diag(np.power(2,np.arange(-1,-(n+1),-1), dtype=float)) # Set S to a diagonal matrix w/ exp \n", " # values between 2^-1 and 2^-(n+1)" ] }, { "cell_type": "code", "execution_count": 184, "metadata": { "collapsed": true }, "outputs": [], "source": [ "A = np.matmul(U,np.matmul(S,V))" ] }, { "cell_type": "code", "execution_count": 188, "metadata": { "collapsed": true }, "outputs": [], "source": [ "QC, RC = cgs(A)\n", "QM, RM = mgs(A)" ] }, { "cell_type": "code", "execution_count": 189, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAJKCAYAAACh7AgvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlYlOX6wPHvCwrImopbKZu5kLikmJppqFSWWWpKFqUo\n5VaZ/TJLzf2IuVRkmUsmmmEnyXI7LS6BnsrTUUvrqGkpQpkLouKCO8/vj5GRYWZgBmaF+3Ndc73O\nM+9yz+LMzbNqSimEEEIIIYT783B2AEIIIYQQwjYksRNCCCGEqCAksRNCCCGEqCAksRNCCCGEqCAk\nsRNCCCGEqCDcOrHTNO19TdP+1DTtrKZpRzRNS9Y0zcvZcQkhhBBCOINbJ3bAe0BTpVQg0PLGbZxz\nQxJCCCGEcI4qzg6gPJRSe4vc1YACoJGTwhFCCCGEcCqXq7HTNK2xpmlTNU37j6ZpOZqmndM0bZem\naeM1TfMzsf9rmqadB06gq7FLdnjQQgghhBAuQHO1lSc0TXsDeA5YC/wHuAp0AeKAX4D2SqmLJo6L\nBOKBBUqpv0q6RnBwsAoLC7Nx5EIIIYQQtrdz586TSqlaluzrik2xnwEzlFJ5RcoWaJr2OzAeSETX\nt86AUmqfpmm7geXoEkGzwsLC2LFjhw1DFkIIIYSwD03Tsizd1+WaYpVSO4oldYU+vbGNKuHwqkBj\n20clhBBCCOH6XC6xK0H9G9vjAJqmBWmalqBp2i2aTgvgdeAbp0UohBBCCOFEbpHYaZrmCUwArgEr\nbhQr4CngEHAOWA18CbzgjBiFEEIIIZzNFfvYmZIMdADGKaX2AyilzgKxlp5A07QhwBCAkJAQe8Qo\nhBBCCOFULl9jp2naNOB5YJFSakZZz6OUWqSUilZKRdeqZdHAEiGEEEIIt+LSiZ2maZPR9ZtLAYY5\nNxohhBBCCNfmsondjaRuErAMeEa52oR7QgghhBAuxiX72GmaNhFdUrccGKyUKnBySEIIYbW8vDxO\nnjzJlStXnB2KEMLFeHl5ERwcTFBQkE3P63KJnaZpzwFTgGxgE/CkpmlFdzmulNrojNiEEMJSly5d\n4vjx49SvX59q1apR7HtMCFGJKaW4ePEif/31F97e3vj4+Njs3C6X2AFtb2xD0DXDFrcFkMROCOHS\ncnJyqFWrFr6+vs4ORQjhYjRNw9fXl+DgYHJycmjQoIHNzu1yfeyUUglKKa2EW4yzYxRCiNJcunQJ\nf39/Z4chhHBhAQEBXLp0yabndLnETgghKoJr165RpYorNooIIVxFlSpVuHbtmk3PKYmdEELYifSr\nE0KUxB7fEZLYCSGEEEJUEJLYCSGEEEJUEJLYCSGEEKW4dOkSmqYxbJhrLoLUvn17mjZt6uwwbCo6\nOpqoqChnh+F2JLETQghRJpqmWXw7fPiws8M1aefOncTFxdGwYUN8fHwIDg6mVatWDB8+nF9//dXZ\n4bmsnj17omkav/32m9l9rl+/zm233UZwcLBM0u1AMmRLCCFEmSxfvtzg/r///W8WLVrEkCFD6NSp\nk8FjtWrVcmRoFlm1ahVxcXHUrVuXAQMGEBERwenTp9m/fz9r1qyhWbNmNG/eHAAfHx8uXrwoI51v\nSExMZP369aSkpDBz5kyT+3zzzTf8/fffvPjii3h5eTk4wspLPqFCCCHK5KmnnjK4f+3aNRYtWkSH\nDh2MHjNHKUV+fj5+fn72CLFEr732GoGBgezcuZO6desaPHb9+nVOnTplUGbL1QHcQX5+Pt7e3nh6\neho99vDDD1OnTh2WL19OUlKSyX1SUlIAXRIoHEeaYoUQQjjE119/jaZpfPLJJ7zzzjs0bdoUb29v\n3n33XQDq1q1L9+7dzR73z3/+06D84sWLTJ06lTvuuAMfHx9q1KhBr169LGpCLSgo4ODBgzRr1swo\nqQPw9PQ0qGU01ceuaNnWrVu555578PX1pVatWgwbNoz8/Hyj827atIl27drh4+NDvXr1GD16NLt2\n7ULTNN544w39fgsWLEDTNP7zn/8YncPS/nQ//PADAwYMoFGjRvj6+hIYGEjnzp1Zv3690b79+/fH\nx8eHY8eOMWDAAGrXro2/vz85OTkmz12lShUGDBjA0aNH+eqrr4wez83NZe3atURHR+trPQE++ugj\nevToQYMGDfD29qZ27dr069ePffv2lfp8AIKDg3n44YeNytevX4+maXz22WcG5fn5+UyePJnIyEj9\nZ6R3797s2bPHYL9r164xa9YsoqKi8Pf3JygoiMjISIYMGWJRXK5EauyEEMIdbdsGGRkQEwMdOjg7\nGqvMnDmTvLw8Bg8eTO3atYmIiLD6HJcvXyY2NpadO3cycOBARo4cyalTp/Q1ht9//z0tW7Y0e7yH\nhwfh4eH8/PPPbN++nbZt25rdtzT//e9/SUtL45lnnuGpp55i8+bNLFy4EC8vL+bOnavfb/PmzTz4\n4IPUrl2bcePGERAQwD//+U+2bNlS5muXJC0tjUOHDtG/f39CQkLIyclh6dKl9OzZk1WrVtGnTx+D\n/QsKCujWrRvh4eFMnDiRc+fOUa1aNbPnHzx4MLNnzyYlJcUo2UpNTeXKlStGtXXvvPMO4eHhDB8+\nnFq1anHgwAE++OADNm3axO7duwkJCbHZ87906RJdu3Zl165dJCQk8OKLL5Kbm8vChQtp374927Zt\n0w/OGD9+PLNmzeKxxx7j+eefB+DgwYOsXr3aZvE4jFKq0t3atGmj7O6HH5RKStJthRCVzt69e+13\n8h9+UKpaNaU8PXVbF/meSUlJUYBKSUkx+fhXX32lAFWrVi2Vm5tr9HidOnXUAw88YPa4Tz75RF+W\nlJSkPDw81Lfffmuwb25urqpbt67J8xS3fPlyBShN01TLli3V8OHD1ZIlS1RWVpbRvhcvXlSAGjp0\nqFGZp6en+umnnwz279q1q/L29laXLl3Sl7Vo0UL5+vqq7Oxsfdnly5dVmzZtFKBmzJihL58/f74C\n1LZt24xiadeunWrSpEmpZefPnzc69ty5cyo8PFzdeeedBuWPP/64AlRiYqLRMSW5++67VdWqVVVO\nTo5BeatWrVS1atXUmTNnSo3pp59+Uh4eHuqVV14xKG/Tpo1q1qyZQVnNmjVVjx49jM6xbt06Bai0\ntDR92dSpU5Wnp6faunWrwb45OTmqdu3aBudp1KiRatu2bSnP1j4s+a4AdigLcxxpirWHbdugWzeY\nMEG33bbN2REJISqSjAy4cgWuX9dtMzKcHZFVBg8eTI0aNcp1jo8//pjmzZvTvHlzTp48qb8V1jql\np6eXulTTU089RXp6Or169eLw4cPMnz+fwYMHExYWRt++fY362Jlz7733cueddxqUde3alcuXL/Pn\nn38CkJWVxS+//ELfvn0NFnz38vJi5MiRVj57yxTtt5ifn09ubi6XLl3i3nvvZdeuXVy+fNnomNGj\nR1t1jcTERK5evUpqaqq+bNeuXezatYvHHnuMoKAgkzEppTh79iwnT56kQYMGhIWF8eOPP1p17dJ8\n/PHHtGrVisjISIPPCECXLl3YvHkzBQUFAAQFBZGZmcl///tfm8bgDNIUaw+mvnTdrKlECOHCYmLA\ny0v3/eLlpbvvRho3blyu45VS7N+/n+vXr5c42vb06dOljsaNiYkhJiZGf8709HTee+89Vq1axdWr\nV1mzZk2p8ZhqSq5Zsyag62t2++23k5mZCUCTJk2M9jVVZgtHjx5l/PjxrFu3Tp/QFJWXl0ft2rX1\n9z08PGjYsKFV14iLi+PFF18kJSWFF198EYAlS5YAugS+uB9//JGJEyfy3XffGfVBtOUAmuvXr/P7\n77+jlCrxM5CXl0f16tWZNWsWffv2pV27djRo0IAuXbrQo0cP+vTp43Yjod0rWnfh5l+6QggX16ED\nbN7stn3sfH19TZabWzezeM2butmtxmDAQXHFa4tKomkaTZs2pWnTpgwYMIDIyEjWrl3LyZMnCQ4O\nLvFYUyNCi8ZaFiWtIWrJovHXr1+nW7duZGZm8uKLL9KmTRuCgoLw8PBg4cKFfPbZZ/raqkJVq1al\natWqVsXp7+9PXFwcS5Ys4aeffiIqKooVK1YQERFBTLHfvgMHDhATE0OdOnWYMmUKt99+O35+fmia\nxtChQ43iMcXaz0jbtm1JSkoqMX7Q1eBlZmby9ddfk56ezrfffstHH31Eq1at2Lp1KwEBAaXG5iok\nsbMHN//SFUK4gQ4dKtx3S40aNUw2fx46dMjgfmHN0smTJ+nWrZvNF1L38/OjefPm/Pnnnxw5cqTU\nxM4SYWFhAOzfv9/oMVNlhU3VxV8PpRSZmZml1kTu2LGDffv2kZSUxNixYw0ee++996wJvVSJiYks\nWbKElJQU7r33XnJzcxk1apTR+5KWlsalS5dIS0szGKxSUFBATk6ORQMnLP2MVKlShfDwcHJzc4mN\njbXoeQQGBhIXF0dcXBwAs2fPZsyYMXz88ccMHz7conO4AuljZy8dOsDYsYZfvNu2wYwZ0udOCCFM\naNy4Mb/++isnTpzQl128eJH58+cb7TtgwACysrKYN2+eyXMdP3681OuZmqYDdE2YP/74I97e3mUa\nsWtKWFgYUVFRfPbZZ/p+dwBXrlwxGDlbqLC5etOmTQblS5cutajvX2EtYvEaw59++ol//etfVsdf\nkrvvvpumTZuyYsUKFixYgIeHBwkJCRbH9M4773Du3DmLrtW4cWN2795t0LScn5/PggULjPYdMGAA\nhw4dMvkYGH5GTDVVF/abtLSvpauQGjtHKRxQUdg8u3lzhftrWwghyuP5559n9erVdO3alSFDhnDx\n4kWWLVtmskn1lVdeYfPmzbzwwgt88803xMTE4O/vT3Z2Nhs3bqRmzZpmE7dCPXv2JCwsjIcffpim\nTZvi6enJwYMH+eijj8jNzSUpKcmmTXBvvfUWDz74IO3bt2fYsGEEBATwySef6B8vWsPVsmVL7rnn\nHt555x2uXr1KVFQUO3fuZP369frav5K0aNGCxo0b849//IMzZ87QqFEj9u3bxwcffECLFi346aef\nbPa8QFdrV/iedO/enfr16xvt07NnTyZNmkRcXBwjRozA39+frVu3kp6ebvE0J88//zzr16+na9eu\nPPvss+Tn55OSkmJyMM5rr73Gt99+y/Dhw/nqq6/o3Lkzfn5+ZGdns2HDBurVq8e6desACAkJ4YEH\nHiA6Opp69erx119/sXDhQnx9fenbt2/5XhxHs3T4bEW6OWS6k+KSknRTE4Bum5Tk+BiEEA5j1+lO\nXJSl050UnbakuA8++EDdfvvtqmrVqioiIkK99dZb6ssvvzR53JUrV9Sbb76pWrdurXx9fZWfn59q\n1KiRevrpp9XmzZtLjXfFihVq4MCBKjIyUgUFBakqVaqoOnXqqIceekitWbPGYN+SpjspWlbI3HQl\n33zzjYqOjlbe3t6qTp066sUXX1RbtmxRgHrnnXcM9v3rr79U7969lb+/v/L391c9evRQBw4csHi6\nkz/++EP17t1b1axZU/n6+qp27dqpdevWqVdffVUB6ujRo/p9H3/8ceXt7V3qa2bO8ePHVdWqVRWg\nVq5caXa/DRs2qHbt2ik/Pz9VvXp19eijj6r9+/ebnNrEVJlSute2YcOGqmrVqqphw4YqOTlZrV27\n1mi6E6V008nMnj1b3XnnnapatWr6z8iAAQMMpsqZPHmyuvvuu1VwcLDy8vJSDRo0UE888YT65Zdf\nyvyaWMrW051oqowdO91ZdHS02rFjh2MvKjV2QlQq+/btIzIy0tlhCDeQmprKU089xRdffEGvXr2c\nHY5wMEu+KzRN26mUirbkfNIU6ygyoEIIISq1goICrl27hpeXl77s8uXLJCcn4+3tTefOnZ0Ynago\nJLFzgFmzIK/mJlLPPUP2lWxCfgwhfu9ignJjGdPJfZcFEkIIYbmzZ88SGRlJfHw8jRs3Jicnh08+\n+YQ9e/YwadKkck/aLARIYucQeTU3kTSyJfQLh/AssnaFk5TWknFjlkK3EdI8K4QQlUC1atW4//77\n+fzzzzl27BgATZs2ZdGiRTz77LNOjk5UFJLYOUDquWd0SV3aSoieDzuGQ784UtnFdFmhQgghKgVv\nb2+WLVvm7DBEBSeJnQNk52VDeJYuqds6ETpPhfAMstHAy0dWqBBCCCGETUhi5wAhQSFk7QrX1dR1\nnqrbhqcT0ioTNn8ifeyEEEIIYROS2DlAfMBiktJaQr84CM+A8HRIW0l8p92mlwXaJgMqhBBCCGE9\nSewcICg3lnFzN5F6LpPsPI2QVpnEd9pNUK6J9etkvjshhBBClJEkdg4wZgxALNM5XPrOGRm6pE4G\nVAghhBDCSh7ODkAUExOjq6nz9JQBFUIIIYSwitTYuRpZoUIIIYQQZSSJnY0ZrDKRl01IUAjxATdW\nmRhj4UlkQIUQQgghykCaYm2scJWJrF3hKJRulYmRLcmruansJy0cUDFhgm67bZvtAhZCCDe3dOlS\nNE0jIyOjxDKAzMxMevXqRa1atdA0jYSEBACDf9s7NmGZsLAwYhzcHcmaa7rqeyuJnY3pVpmI060y\n8e0U3bZfnK68rEwNqBBCCBeQkZGBpmlomsbzzz9vcp8TJ07g5eWFpmkO/6EuLiEhgS1btvDqq6+y\nfPlyhg4d6tR4TDl48CCjRo0iKiqKwMBAvLy8uPXWW3nooYeYP38+Fy5ccHaIZXbs2DFGjx5NVFQU\nAQEBBAYG0qhRI/r378/nn3/u7PDKLSMjg8mTJ3PmzBmnxSBNsTZmdpWJPK3sJy0cUCErVAghXJSP\njw8rVqzgzTffxNvb2+Cx5cuXo5SiShXH/eQ8/fTT9O/fHy8vL33Z5cuX+fe//83zzz/P6NGjDfa/\nePEinp6eDovPnKVLlzJs2DCqVKlCXFwcw4YNw9fXl2PHjrF161aef/55Vq9ezTfffOPsUK2WlZXF\nXXfdxdmzZ4mPj2f48OEA/PHHH6Snp5OSkkKfPn2cGuP+/fvRtLL/XmdkZDBlyhQSEhK45ZZbbBiZ\n5SSxs7ESV5koKxlQIYRwcb179+aTTz5hzZo1xMXFGTyWkpLCQw89xObNmx0Wj6enp1Gidvz4cZRS\n1KhRw2h/Hx8fR4Vm1ubNm0lMTCQqKoovv/yS2267zeDxcePGcejQIT799NNSz3Xu3DkCAgLsFWqZ\nzJkzhxMnTrB69WoeffRRo8ePHTvmhKgMFf+jxB1JU6yNxQcs1je/0nWSvlk2PmBx+U7coQOMHWuY\n1G3bBjNmSJ87ISqZ1F9TCUsOw2OKB2HJYaT+murskGjdujUtWrQgJSXFoPy///0ve/bsYdCgQWaP\nXb16NR07dsTPzw9/f386duzImjVrTO77wQcf0LRpU7y9vbn99ttJTk5GKWW0X/H+TwkJCYSGhgIw\nZcoUffNx4ePm+tht2rSJ+++/n1tuuQUfHx9atGjBggULyhWbOWNujLBbuXKlUVJXKCIigrFjxxqU\nFfYL+/nnn3nggQcICgqiRYsWgC7Be/3112nXrh3BwcH62F577TXy8/MNzlPYrL506VLef/99mjRp\ngo+PD82bN2f9+vUA/Prrr3Tv3p3AwEBq1qzJyJEjuXr1qkXP7/fffwegW7duJh+vW7euyfLffvuN\nHj16EBAQQFBQEH379jVKAidPnoymaezdu5dRo0ZRr149fH196datG/v37wfg888/p3Xr1lSrVo2w\nsDAWLVpkdC1zfewseW8TEhKYMmUKAOHh4frP2OTJk0t8XWxNauxszNJVJso9elZWqBCiUkr9NZUh\n64aQf1X3o5yVl8WQdUMAiG8e78zQGDx4MP/3f//HkSNH9InJkiVLqF27Ng8//LDJY95//32ee+45\nmjZtysSJEwFdUtarVy8WLlzIkCFD9PsmJyfz0ksv0bJlS5KSksjPz2fOnDnUrl271NiGDh1Kq1at\neOmll+jdu7e+yS8yMtLsMYsWLWLYsGG0b9+e8ePH4+fnx8aNGxk+fDgHDx5k9uzZNokNdIM6fvrp\nJzp37kyTJk0sOqao7OxsunbtSr9+/Xjsscc4f/48AEeOHGHx4sU89thjPPnkk1SpUoUtW7Ywa9Ys\nfv75Z5NNuvPmzeP06dM888wz+Pj4MHfuXHr37k1aWhrPPvssTzzxBL169WLDhg28++671K5dm9df\nf73UGBs2bAjokqRRo0ZZ1OR55MgRYmJi6N27N7Nnz2b37t0sXLiQs2fPsmHDBqP9Bw4ciL+/P+PG\njSMnJ4c333yTBx54gGnTpjFmzBiGDx/O4MGD+fDDDxk6dCh33HEH99xzT4kxWPreDh06lLNnz/LF\nF1/w9ttvExwcDKBPsh1GKVXpbm3atFHONm7xRoXvCcXAGMVkdFvfE2rc4o2WnSApSSlPT6VAt01K\nsm/AQgir7N271y7nDX07VPedUewW+naoXa5XmvT0dAWo2bNnq5MnTyovLy81ffp0pZRS+fn5Kigo\nSL388stKKaX8/PzUvffeqz/21KlTys/PTzVs2FDl5eXpy/Py8lRERITy9/dXp0+fVkopdfr0aeXr\n66siIyPVhQsX9Pv++eefys/PTwEqPT1dX56SkmJUlpmZqQA1adIko+cBqIEDB+rv//3338rb21s9\n8cQTRvuOHDlSeXh4qIMHD5YpNlPWrl2rADVy5Eijxy5cuKBycnIMbgUFBfrHQ0NDFaA++OADo2Mv\nX76srly5YlT++uuvK0D9+OOP+rLC9/LWW29VZ86c0Zfv3r1bAUrTNLVq1SqD87Ru3VrVrVu3xOdW\n6ODBgyowMFABqkGDBurJJ59Ub7/9ttqxY4fJ/Quf16effmpQPmLECAWo3377TV82adIkBaiHH37Y\n4LV55513FKACAgJUdna2vvzEiRPK29tb9e/f3+iaRT+j1r63hXFkZmZa9JooZdl3BbBDWZjjSFOs\nk5R79KysUCFEpZSdl21VuSPVrFmTRx55hKVLlwK6pq+8vDwGDx5scv+NGzdy4cIFRo4cSWBgoL48\nMDCQkSNHcv78eTZt0k0VtWHDBvLz83nuuefw9fXV71u/fn3i421fU/nZZ59x+fJlEhMTOXnypMGt\nZ8+eFBQU2DS2s2fPAhi8DoUmTpxIrVq1DG65ubkG+9SoUcNkc7eXlxdVq1YF4Nq1a5w+fZqTJ08S\nG6trRfrxxx+NjklISCAoKEh/v0WLFgQGBnLrrbcaDW645557OHbsmL6GsCQRERHs3r2b5557DoAV\nK1bw0ksvER0dTYsWLdi5c6fRMbfeeqtRn82uXbsCN5t2ixo5cqRBTWCnTp0AeOSRR2jQoIG+vFat\nWjRp0sTkOYpy9OfOFqQp1knKPXpWBlQIUSmFBIWQlZdlstwVDBo0iB49evDdd9+xZMkS7rrrLu64\n4w6T+2Zm6gaVNWvWzOixwrJDhw4ZbJs2bWq0r7nzl8e+ffsA9AmQKcePH7dZbIUJXWGCV9TQoUPp\n3r07ALNnzzbZBNmwYUOzo3rff/99FixYwJ49eygoKDB47PTp00b7R0REGJVVr17dIDEqWg6Qm5uL\nv78/58+fN0ryatSooR+dHBYWxnvvvcd7773H0aNH+e6771i+fDnr1q3j4YcfZs+ePQaDW0zFUrNm\nTf01S4u9ML7w8HCTsWdlGf9fKsrRnztbkMTOSUyOnr1Qh4A2XxGWHGZZvztZoUKISmd6t+kGfewA\nfKv6Mr3bdCdGddMDDzzAbbfdxpQpU0hPT2f+/PnODqlM1I2O8R999BH16tUzuY+ppKOsoqKiANi1\na5fRY40aNaJRo0YAfPzxxyaPL1qbVNRbb73Fyy+/zP3338/IkSO59dZb8fLy4siRIyQkJBgleoDZ\nBLGk6WAKX685c+boBxAUSk9PNzkgoV69evTr149+/foRHx/PihUr+PLLL3nqqaesumZ5Yjd1Dncn\niZ2TxAcsJimtpa45NjwDwtPhn19w9n+Pc7Z/bwjP0q1akdaScXM3Aeb/atSTARVCVHiFAyTGbx6v\n/wNwerfpTh84UcjT05MBAwYwY8YMqlWrxhNPPGF238LEaM+ePUYjJffu3WuwT+H2t99+M7uvLRUm\nUsHBwSXW2tkqtvDwcFq3bs13333H/v37yzSAwpTly5cTFhbGV199hYfHzd5XX3/9tU3OX9yAAQOM\nBiO0bNmy1OPat2/PihUrOHLkiF3iKitr39vyzIFnK9LHzkl0o2d3E9oqEw2N0FaZBA5MgGaflr3f\nnaxQIUSlEN88nsOjDlMwqYDDow67TFJXaNiwYUyaNIkFCxaY7DNW6L777sPPz493332Xc+fO6cvP\nnTvHu+++i7+/P/fdd59+32rVqjFv3jyDaTr++usvVqxYYfPnEBcXh7e3N5MmTeLixYtGj+fl5XH5\n8mWbxjZz5kz9tf/++2+T+1hbw+Tp6YmmaQbHXbt2jTfeeMOq81gqIiKC2NhYg1thc2hGRobJ17Kg\noIB169YBrte8ae176+/vD8CpU6ccFmNxbl1jp2laODAXuBu4BiwBxiuljOuWXYyuaTWW6RzWl3lM\n8YBH1oD/8RL73ZmdKuWvJxnjNU1WqBBCOFVISIhFc3fdcsstzJo1i+eee4527drp55FbunQpf/zx\nBwsXLtR34q9evTrTpk1j9OjR3H333QwYMID8/HwWLFhAo0aN+Pnnn236HOrXr8/8+fN55plniIyM\n5OmnnyY0NJScnBx+/fVXVq9ezd69ewkLC7NZbLGxsfppOBo3bky/fv1o06YNvr6+HD9+nK1bt7Jh\nwwbq1atn8YTKffv2ZezYsTz44IP06dOHs2fPsmLFCv2ACkeaM2cO33//PT179qR169YEBQVx7Ngx\nVq1axc6dO+nSpQs9evRweFwlsfa9bd++PQCvvvoq8fHx+Pj4EBUVpW9qdwS3Tew0TfME1gHfAH2B\n2sB64Aww04mhlZmlq1bk1dxE0siW0C+8WJPtbnhKBlQIIdzHiBEjqFevHrNnz9b3zWrZsiVffPEF\nvXr1Mtj35Zdfxt/fn7feeouxY8fSoEEDRo8eTVBQkNmRt+UxaNAgGjduzJw5c1i4cCFnzpwhODiY\nJk2aMG1FiAqlAAAgAElEQVTaNIMJdW0VW0JCAp06dWLu3Lls2rSJtLQ0rly5QnBwMC1btmTevHk8\n/fTT+Pn5WXS+V155BaUUH374IS+++CJ169bl8ccfZ9CgQQ6vHXv99ddJS0tj69atfPPNN5w6dQo/\nPz8iIyN58803ee655wyai12FNe9tx44dmTlzJgsWLODZZ5/l2rVrTJo0yaGJneauHQc1TbsD+B/g\np5S6eKMsAZiklDIe/lJEdHS02rFjh/2DtNL4DwsTthv97jJjIG0l4+buZnrizT4eYclhugQwbaVu\nVO2O4dAvjtBWmRweddj0yWVQhRAOtW/fvhInvxVCCLDsu0LTtJ1KqWhLzueQGjtN0xoDTwH3Aw0B\nH+AgkAYkK6UulOW0xbaF/w7TNC1QKWU8ZtzFWbpqhdVTpcigCiGEEKJScFRT7GDgOWAtkApcBboA\n/wDiNE1rX6TW7Z/A4yWcq4tSKgPYjy45nK5p2jigLvDSjX0CAbdL7Ez1uyvsTxeWfLM/XcDfj3B2\nx4PwW58Sm2z1TA2qkMROCCGEqHAcldh9BsxQSuUVKVugadrvwHggEXjvRvmzwPMlnCsPQCl1TdO0\nnsDbQBZwCvgQXf864xkX3ZSp/nT8c6nuwf69b06VkraS+E67TZ+kcJUKGVQhhBBCVGgOSeyUUuY6\ntH2KLrGLKrLvOeCcmf2Ln/c34MHC+5qmPQdsL2PTrkvSLT1WrD9d1KcEtvmK6pElN9nqySoVQggh\nRKXg7FGx9W9sj5flYE3TWgCHgEvomnbHAwNtE5prMNmfruskzqGRN+rmrC6mmmwNVq2QVSqEEEKI\nCs9p44pvTFcyAd38c2WdXbIfumbYPOAN4Fml1EbbROgaQoJCdKNji06BkhljtC5kYZNt1q5wFEo3\nBcrIluTV3GT6xIUDKiZM0G23bbP/kxFCCCGEXTmzxi4Z6ACMU0rtL8sJlFIT0CWHpdI0bQgwBHST\nZ7oLk0uPmehPZ7LJtl8cqecyDQZj6BUbUDFr6kXy+pqY9NjcOrVCCCGEcDlOqbHTNG0augESi5RS\nMxxxTaXUIqVUtFIqulatWo64pE2YWnps3FxzU6Bk3GyyjZ5/YwqUbNMnLhxQ4ekJXl7ktc+2rsZP\nCCGEEC7H4TV2mqZNBl4HUoBhjr6+uzE1BYoplq5aAYX98S6QOimQ7EvHCfEJ5PSh1dD0suU1fkII\nIYRwOQ5N7G4kdZOAZcAzyl2XvXBBljbZQtEpVCIh/DhZ+yIh7QPdsaWsUyuEEEII1+WwxE7TtIno\nkrrlwGClVEEphwgrWLpqBZjvjwdYVOMnhBBCCNfkqCXFngOmANnAJuBJTTOoCTpe0UazOpqlTbZg\nZgoV0CV6FtT4CSGEEMI1OWrwRNsb2xB0zbDLi93GOygOgZkpVP7Xn8D+gwitsQVNQWiNLYwb8y/z\nkx4LIUQZxcTEEBYW5uwwANeKpSLKyMhA0zSWLl3qstesaJ8BhyR2SqkEpZRWwi3GEXEInfiAxTdr\n57pO0m339eH5u0dxuN/3FHgncbjf90yflCBTnQghLJKfn09ycjKdOnWiRo0aVK1alTp16vDQQw+x\ndOlSrl275uwQHW7dunX06dOH+vXr4+3tjb+/P5GRkSQmJvLtt986O7xy2bp1K4888ghhYWF4e3tT\nu3ZtoqOjGTlyJIcOHXJ2eOWWnJzs0GTUlpy98oRwghL743VAVqgQQljljz/+oEePHhw4cIDY2FjG\njh1LcHAwJ06cYNOmTQwaNIi9e/cya9YsZ4dqZMOGDdh6HN/Fixd58sknWb16NU2aNGHAgAFERERw\n/fp1Dhw4wPr161myZAkrVqzgiSeesOm1HWH+/PmMGDGCiIgIBg4cSIMGDcjJyWHfvn188skndO7c\nmYiICKfF17lzZy5evEjVqlXLfI7k5GTCwsJISEiwXWAOIoldBVe41JjBxMM1dRMPHx5zuPQTFK5Q\nceWKbt67zZsluRNC6F28eJGHH36YQ4cOsWrVKvr06WPw+Kuvvsr27dvZvn27kyIsmZeXl83POXz4\ncFavXs0rr7zCG2+8gYeHYePYnDlz+OKLL/D19S3xPEopLly4gL+/v81jLKtr164xbtw4QkJC+Pnn\nnwkMDDR4/MqVK5w/f95J0el4eHjg4+Pj1BicyWlLignHsHqpsSJmzYLxH2wgbOhFPF6/TtjQi4z/\nYAMu+Ee3EJXGrFmQnm5Ylp6O0/5fLl68mP379/Pyyy8bJXWF2rZty4gRI0o8z3//+18SEhJo3Lgx\nvr6+BAQE0LFjR7744gujff/8808GDx5MaGiovhnw7rvvZtmyZfp9CgoKSE5OpkWLFgQEBBAYGEiT\nJk1ITEzk6tWr+v3M9a/6448/GDRoEPXr18fLy4tbb72VRx99lJ07d5b4PH755ReWLVtGx44dmTlz\nplFSB6BpGn369KF79+76sqL9wubNm8cdd9yBj48Pc+bMsfr1SUhIQNM0cnNzSUhIIDg4mICAAHr1\n6sWxY8cAWLRoEZGRkfj4+NC0aVPWrFlT4vMqdPLkSc6cOUPbtm2NkjrQJco1atQweWxKSgrNmjXD\n29ub0NBQkzW4YWFhxMTEsHv3bmJjY/H396d27dq8/PLLXLt2jUuXLjF69Ghuu+02fHx86Ny5M/v2\n7TM4h7k+dqdPn+bZZ58lODgYPz8/YmJiTL6fmqaRlZXFli1b0DRNfzt8+LBFr5GzSY1dBWf1UmNF\n5NXcRNKnI6BfBtySQdbpGJLSRjBu7iZABlUI4Qxt20JcHKxcCV266JK6wvvO8NlnnwEwZMiQcp3n\niy++4LfffiMuLo7Q0FByc3NZtmwZffr0ITU1lSeffBLQ1Rjdd999HDlyhBEjRtC4cWPy8vL45Zdf\n+Pe//83AgQMBmD59OhMnTqRnz54MGzYMT09PMjMzWbt2LZcvXy6xmW7Hjh1069aNq1evkpiYSFRU\nFKdOnWLLli388MMPtGnTxuyxq1atAiAxMZFisz9YJDk5mdzcXJ599lnq1q1LgwYNrHp9iurevTv1\n69dn6tSp/PHHH8ydO5fevXvTp08fFi1aRGJiIj4+PsydO5e+ffty4MABwsPDS4yvTp06+Pv7s3Xr\nVvbv30+TJk0sel4LFizg+PHjJCYmcsstt/Dxxx/z6quvUr9+faPY//rrL+677z4ef/xx+vbty4YN\nG3jrrbeoUqUKe/bs4eLFi7z22mucPHmSOXPm0KtXL/bt22cyiS509epVHnjgAbZv387TTz9N+/bt\n2bVrF7GxsdSsWdNg3+XLl/PSSy8RHBzM+PE3x3a6zapVSqlKd2vTpo2qLLTJmmIyis5TFCjddjJK\nm6yVemzo26GKgTEK3xO643xPKAbGqNC3Q+0fuBBubu/evXY797ffKhUcrNSECbrtt9/a7VKlqlGj\nhgoMDLTqmHvvvVeFhoYalJ0/f95ovwsXLqjGjRuryMhIfdnu3bsVoGbOnFniNe68806D4yyNpaCg\nQDVr1kx5e3ur3bt3G+1//fr1Es/Xp08fBaiffvrJ6LHc3FyVk5Ojv+Xl5ekfS09PV4CqXr26On78\nuNGxlr4+Sik1cOBABagRI0YYlL/00ksKUA0aNDC4duFr+tprr5X43ArNmTNHAcrT01O1bdtWjRw5\nUn388cfq6NGjRvsWPq969eqpM2fOGMQeHBys2rdvb7B/aGioAtTKlSsNylu3bq00TVOPPPKIKigo\n0Je/8847ClBff/210TVTUlL0ZQsXLlSAmjhxosF53377bQUYfR5DQ0PVvffea9HrUV6WfFcAO5SF\nOY40xVZwJqc2yYzRlZeixPVnt22DGTN0WyGEQ3XpAsOHw7Rpum2XLs6L5ezZswQEBJT7PH5+fvp/\n5+fnk5ubS35+Pl27dmXfvn2cPXsWgKCgIADS09M5ceKE2fMFBQVx5MgRvvvuO6vi2LVrF3v27GHQ\noEG0aNHC6PGSaoUAfZymmikbN25MrVq19DdTtWwDBgygdu3aRuWWvj5FjRo1yuB+p06d9NcoGl+L\nFi0IDAzk999/L/G5FXr55ZdZu3Yt999/P3v37mXu3Lk89dRT1K9fn8TERPLz842OGTRokP69A/D1\n9aV9+/Ymr3nbbbfRr18/g7J77rkHpRQvvPCCQU1o4XMqLfbVq1fj6enJyy+/bFA+fPhwk++VO5Om\n2ArOmqXGijO7/mzkPhlQIYQTpafD/PkwYYJu26WL85K7wMBAzp07V+7znDhxgtdff501a9aYTNjO\nnDlDYGAgoaGhjB8/nhkzZlCvXj1atWpFt27d6NevH23bttXvn5SURK9evejUqRO33norMTEx9OjR\ng759+5Y4YKIwQbjzzjvL9DwKkwRTidbnn3/OlStXALjvvvtMHt+4cWOT5Za+PkUVH5lavXp1AJPN\nrdWrVyc3N1d/Pycnh+vXr+vve3p6GjRF9uzZk549e3L9+nX27t3L5s2beeedd1iyZAlVqlRh4cKF\nJcYCULNmTYNrFjIXn6nHCstNnaeoQ4cOUa9ePaPXyNvbm4iICE6fPl3i8e5EauwqON3UJrsJbZWJ\nhkZoq0zGzTW91FhxJue7S1tJ/P7huqTu+nXdNiPD/k9ECAEY9qmbOlW3jYszHlDhKFFRUZw9e7Zc\nc5cppbj//vtZtmwZAwcO5NNPP+Xrr79m48aN+lqtgoKbq1D+4x//4Pfffyc5OZmGDRuyePFi7rrr\nLl599VX9Ph06dODgwYN89tln9O7dm127dhEfH0+rVq04depU2Z9wKaKiogBdzV9xnTt3JjY2lthY\n89+/pkbKWvv6FPL09DR5DXPlqsi0L23btqVevXr6W9Gkufi5mjdvzqhRo9i+fTtBQUEsW7bMICks\n6ZrWxGdp7JWdJHYV3JgxMD0xlsOjDlMwqYDDow4zPTHWoomHzSaF1RJ0NXWenrptTIy9n4YQ4obt\n228OnADdduVKXbkzPPbYY4BudGxZ/fLLL+zevZvXXnuNWbNmERcXxwMPPEBsbKxRglAoIiKCF154\ngZUrV/L333/TuXNnZs2aZVCb5e/vz2OPPcZ7773Hnj17mDdvHvv27ePDDz80G0thjZmpxMwSha/H\nhx9+aLNkoyyvT3mlpqayceNG/S01NbXUY4KDg2nYsCGXL1/m5MmTdomrrCIiIjh69KhRTerly5dN\n/lFSloEvrkISO2GW2aRwXqiu+XXaNGmGFcLBxowxbnbt0gWnrRLzzDPP0KRJE+bMmWN2yoydO3fy\n/vvvmz1HYS1M8UTof//7n9F0Hnl5eQbTlQD4+PgQGRkJoG9SM5VYtG7dGqDEGruWLVvSrFkzlixZ\nwp49e4weLy1Za9GiBQMGDOD777/ntddeM1mTZm3CZ83rYysdO3bU1y7GxsbSsWNHQNe/b8uWLSaP\n+f3339m7dy/BwcEuN4L00Ucf5fr167z55psG5fPnzzfZbO7v72/Xml17kj52omw6dJAVKoQQ+Pr6\nsn79enr06EGvXr24//77ue+++6hZsyY5OTmkp6fzzTff8Morr5g9R2RkJM2aNWPWrFnk5+fTpEkT\nDhw4wMKFC2nevLnBXGPp6ekMGTKExx57jCZNmuDv78/OnTtZvHgx7dq100+/ERkZSfv27WnXrh23\n3norR48eZdGiRXh5edG/f3+zsWiaRkpKCt26deOuu+7ST3dy5swZtmzZQvfu3XnhhRdKfE0WLFhA\nXl4es2bNYs2aNfTp04eIiAiuXr1Kdna2foqY0qYWKcvrY2/5+fnExMQQFRVF9+7dadSoEUopfvvt\nNz766CMuXbrEvHnzSh1k4miDBg1i0aJFTJ06lczMTDp06MDPP/9MWloaDRs2NFryrn379nz44YdM\nmDCByMhIPDw86Nmzp8EgFlcliZ2wDVmhQohK6/bbb+fnn39m4cKFrFq1iunTp3P+/HmqV6/OnXfe\nSUpKCvHx8WaP9/T05F//+hejR49m2bJlXLhwgaioKJYtW8bu3bsNEpeWLVvSp08fMjIySE1N5fr1\n64SEhDBu3DiDEY8vv/wyX375JXPnziUvL4/atWvTvn17xo4dS8uWLUt8Pm3btmX79u1MmzaNlStX\nsmDBAoKDg7nrrrv0NVclqVatGl988QVr165l6dKlLFu2jJycHKpWrUqDBg3o1KkTixYtoouFI16s\neX3s7ZZbbmHJkiVs2LCBtWvXcvToUS5dukStWrW49957eeGFFyx+Xo7k5eXFxo0beeWVV1i9ejWr\nVq2ibdu2bNy4kdGjRxtNPjx9+nROnTrFvHnzOHPmDEopMjMz3SKx0ypjh8Po6Gi1Y8cOZ4fhlkwu\nURawmKDPPBiz8X7dgApPT10z7dixzg5XCKfZt2+fvnlQCCHMseS7QtO0nUqpaEvO51p1pcLlmV2i\nrH22DKgQQgghnEyaYoVVzC5RFpTJ9M2bpY+dEEII4USS2Amr6FajyLq5GkXnqTdWo9BkQIUQQgjh\nZJLYCauYXY2iVabxzjKgQgghhHAo6WMnrGJ2NYoAE5OTZmTIChVCCCGEA0liJ6xi1RJlMTEyoEII\nIYRwIGmKFVbRzW4fy3QOl75zhw665lfpYycqKaWUWy9NJISwL3tMOSeJnbAvGVAhKqkqVapw7do1\nqlat6uxQhBAu6tq1a1SpYttUTBI74VgyoEJUEj4+PvrVF4QQwpRz587h4+Nj03NKHzvhWDKgQlQS\ntWrVIicnh/z8fLs0twgh3JdSivz8fE6ePEmtWrVsem6psROOVTigorDGTgZUiArKx8eHOnXqcOzY\nMS5fvuzscIQQLsbb25s6derYvMZOEjthF2bXlM2NZYwMqBCVRFBQEEFBQc4OQwhRiUhiJ8rNVBJ3\n24E5/PDRg/BkOIRn6daUTWvJuLmboEOs6YROBlUIIYQQ5SKJnSi3vJqbSBrZUreG7I0kLivtXujy\nuvGasucyTU+VIoMqhBBCiHKTxE6UW+q5Z3RJXbEkjvAMuFTdeE1ZU0wNqpDETgghhLCKJHai3LLz\nsiE8S5fUFUniyIyxbE1ZsGhQRYn99sbY8QkKIYQQbkISO1FuIUEhZO0KN0zifE7Dd+Nu1tyFp+vW\nlO202/RJLFilwlSTr77fHiaWNBNCCCEqGUnsRLnFBywmKa2lYRK3Yh13D/iaI1GZZOdphLTKJOL8\nGbZ9FUrYuTDTNW6lrFJhrsnXbL89IYQQopKRCYpFuQXlxjJu7m5CW2WioRHaKpNx8//Do437cnjU\nYQomFXB41GE6PJhF+le3kLUrHIXS1biNbElezU2mT1w4oGLCBOjW7UaTb8bNJt/o+Tf67WU79PkK\nIVzDrFmQnm5Ylp6uK5d4HMeVnrcrxeIsktiJchszBqYnxhokcdMTjfu96Wrc4nQ1bt9O0W37xenK\nTSk2oCKEION+e5kxhASF2PspCuEU8iNVsrZtIS7u5muUnq6737atxONIlj5vR3yey/seVIj/c0qp\nSndr06aNEo6nTdYUk1F0nqJA6baTUdpkzfQBP/ygVLVqSnl6KlWtmho3OUXhe0IxMEZ3noExCt8T\natzijY59IkKU08yZSn37rWHZt9/qyouXBQff3Lf4fXHzNZkwwfxrY+nr7ah4KiJLnrejPs/leQ9c\n9f8csENZmOM4Pclyxk0SO+cIfTtUn4zReYo+SQt9O9T8QT/8oFRSklI//KBmzlRq3OKNKvTtUKVN\n1lTo26Fq3OKNdvlyFsKerPnxqKyJgjUmTND9mk2YYPpxR/9YlxZPRWXJ83bU59mSWMwl/EOGuN7/\nOUnsJLFzSeMWb7RNjVuRZE8Id2XND1xlTRQsYenr6KiEwpUScVetqSz6ebZHjNZ+Jkwl/K72f04S\nO0nsXJJNatyKNc9KcifcmSvVcLgja2vi7P1j7WrNeOWJx5qEqzw10G++advXzNrnbOr/lyv+n5PE\nThI7t2Iq4evy2AHV5bEDxklg9826pA5026QkZ4cvRJm4Up8ke3BEbVFZko/SfqzLE7epY4cM0d3K\ncj5bKGuSYs1nr7x9RguTO1skUmV5/4om/K76f04SO0ns3IrJJlrv07pb8WbbySlSYyccztZJiqU/\nHvZIjhzVPOdKP5BlqVFydA1Sed+Xko4va02lrWuu7BFjeRV/jkOGOK752hqS2Eli51bMDaowO9BC\n+tgJB7P1j70j+z4VZ+q5+Prqak1sHY8tEwNb16SVdKytExpH1M7aqzbMEQlXWV/v8v4/cqU/Pkoj\niZ0kdm7F3DQoVk+NIsmesCNX6ndjqx80e/VzKspWiYGjf4Rt3cHfEf0pbf2+OuIzX573tbyfCWf+\ngWUtSewksXMrVtfYFScDKtyGO32RmuIqI+Vs0bxX/LnYo0nKGTVf1rB0ugtHJkjl/YyVJSE1td+b\nbyrl52f/RNrWf6S4Ym2bLUhiJ4mdW7Gqj52pqVGSkmRAhZtwVtOHLRJKV/sBsSQRM/fjbK55zpad\nyO31XtsyuTYVY2CgUkFBtmvStOZ1sHWNXXlitFfzvD3Y+w8uV/iDtEIkdoA38AFwCDgHHABeKLZP\nHPAdcB44bOm5JbFzLVaNijX1H0lq7NyKMxKk8v5wuWpfHEsSseIJibnaJ1OJiyXvlbkfvQcftP2P\noT0+O9bUVJYlgXDUKiO2Ot5VVvCwlKs3F9tKRUns/IBpwO3o1rRtBRwH4orscx/QH3hRErtKTvrY\nuRVnNGkW/wGwpnnNXX7QzP3IldY8V1Kzm6us6mDpdco73UVp17dXAmHPUbGWsuV77Yj/M45MuJxd\nY18hEjuTwepq8OaaKO8liZ0wIsmeS3LmF2Rp/cqcXQNnqZJ+0MryHMu7tJIjXkd71XxZEru7N1Va\nwtL30Nr9bJV0ucIcgc7sY+vwxA5oDEwF/gPk3Gg63QWMB/xsdI2qwG/AMyYek8ROGJLmWatY+qVZ\n3i9SZzZpWFKbZQ+OnIvOlp3+rX2vbPk6OqpDfXlqAcszuMCa5+eKNV+WvtdlTfidOZjDHGf/EeiM\nxO6NG8lcKvACMAz4FFDAbqBakX3/eaPc3C3GzDUWAtsBLxOPSWInDMmACquY+mI31YncXMfy8nxh\nW/ojVZ5jzf1wlWeOL0f1nbKUrWuV7LGqg6Vs8ZpZknw4a0SmNc/PEbWF9nyvy5Lw2+P/a3k48w/S\nQs5I7KKBIBPl/7iRrD1fpCwACC7hVtXEed4CfgGCzVxfErtKyuzAi9jtKvQlTWmTUKEvaapL7HbL\nB2NUUqa+sC0tc2R8tqohKW8NQFl+nF2hqdLWyvujV95m4JJicsRn1BGrOhTftzw1sfb4A6m0GlFb\nvofOaA51hT62LtPHDmh+I7FbUI5zJAO/ArVK2EcSu0rK4qlSrJk+pRIz9aVpaVlZOatmyNFToLjK\nHHi2Zs/Z/21Z22OP5K68n0drnp+t+ofa+g+k0v6v2rLW1Vl/VJri6GTPlRK7B28kdlPKePxc4H/m\nkjrAE/AB+gFZN/7tXdp5JbGrOKyZ3NjiCY8rKWfV2Fn7A1CWCVjtyRVGVDqKvUZumqqdK+tr5i5r\n4Zanxq68tVfuUINsy1pKe3B086xLJHY3kq4fgKtAkzIcH3ojKbyEbp66wttXRfZJMNFHr9SaO0ns\nKg5rliOzaomySsbUl5Q9+tiVdn1rR+TZ48ve1jWIjv4BsCdbJTOl1c65w2vmqCZNc/u6wxqwZWXq\nObviKGRH/sHmKonduzcSrbH2uoaV8QwBdgA7QkJCyv0iC9cgNXa24ahRsSUp6xxatu5QbemPrqX7\nuUKtoi2V98fMklpgWyxn5srKOyrWVv1DXbUG2Z3+zzgqQXZ6YoduYmEFLLTH+ct7kxq7ikP62Lkm\na7+YLfmhsWbd0/IqbzwVXXlfb3ernXM1jhwAIcyrNDV2wOQbSd0SQLP1+W1xk8Su4rB0OTJ92Rt1\nlDYZFfpGnUo7KtbV5sWyVfOerb9cy5K8VIZkz9Y1dhW9ds7VVIbPqCM4OkF2WmJXJKlbCnjY8ty2\nvEliJwxUshUqHPWFZGkC4Iq1D2VNXip6bYg9+9gJ4U4qxahYYOKNpO4jV07qlCR24oaZM5UaNznF\nYL67cZNTKsVfro5qQnDGqg7OXh3D1fsvlYcrrGcqRFm4+2fPGRMUP3cjqcsCBgBPFbvdZ4vr2Oom\niZ1Qykz/PDfud2ftF5e9ky5nJTjlTcwcsZi6va4rhDDN3WuLnZHYLTUx7UjRW4YtrmOrmyR2Qinz\nI2rddaRsWfq12SvpcvaXqDNrzaQZVwjX5M616U4fFevqN0nshFLm58Bzpbnt7DG61NYJhD2aQ23B\nGfN0STOuEK7NlefvK4k1iZ0HQlRSIUEhkBkDO4ZD56m6bWaMrnzbNpgxQ7d1orZtIS4O0tN199PT\ndffbtjW9f5cuMHw4TJum23bpYrzP9u2wcuXNx7p00d3fvt12MX7+OfTvbxzbmDFlu4a10tNh/nyY\nMEG3LYzN3sr72lry/gkhysZZ3wsOZ2kGWJFuUmMnlCqhj93kFKWqVVPK01O3dfJoWWtqcZzdr81W\n13XFkbKOIDV2QtiHO38vKCU1dkJYJCg3lnFzdxPaKhMNjdBWmYybu5ug/4TAlStw/bpum5FRruvM\nmmX8l2F6uq7cEpbW4hTW5q1cCVOn6rZFa9LsydY1TdbWVBZl6xpJR3Hm+ydEReeu3wtlYmkGWJFu\nUmMnzDE3BcqDDzqvBskR88GVlz1qmipb7ZWMihVCmIMMnpDETpSNuebZJ0fvUAHVL6o6I/orbbKm\n6ozorwKqXyxzcmZqtn1TAw5MrQnpaoth27OJw107OgshhC1JYieJnSgjc1Og1JxZU3kndjco907s\nrj7+5WOLz13a+piBgUoFBZWexFmzALgjaoHsdY3KVmMnhBDmSGIniZ0oI3NToOhvxcotnfPOVJJi\naZml5ytpP3frMOyucQshhD1Yk9jJ4AkhijA3BQpgsjw7L7vUc5rrFA/GAw4sHYRgzX6F15s48WYc\nrnYpoacAACAASURBVD6NRqXq6CyEELZkaQZYkW5SYyfMMdfHzvfh8SbL64zoX+o5zTVVDhli/xq7\nQtJXTQgh3BdSYydE2ZibAqXRmeF4PzEAwjN0O4Zn4P3EALp4Gc64a2pqk7ZtTdc0ff65YS1er17Q\nu3fp011YOy1GpZmUUwghBJouEaxcoqOj1Y4dO5wdhnAzqb+mMn7zeLLzsgkJCqFD1hdEROWSevIp\nsi8dJ8SnDh1P/os189uwbp2u+bBoEla0+XPWLF3CV7Rs6FDdduHCm2Xp6bqksOiKDaaONbVfYXnR\n65uLRwghhOvSNG2nUiraon0lsROibMZ/uImkkS2hX5yuJi8zBtJW8uSIbDYsbcPw4boaMmcmUdYk\ngUIIIVyTJHalkMRO2EJYchhZu8IhbSVEz9cNqOgXR2irTAacOsy0abrmz6lTnR2pEEIId2ZNYid9\n7IQoo+y8bF1NXfR82DpRtw3PIGtXhPRpE0II4RRVnB2AEO4qJChEV2NXdAoUn9N4fP86K9fdnL5E\n+rQJIYRwFKmxE6KM4gMW65ph+8VB10m6bfo0+g8/TBefbTBjBl18tsn8a0IIIRxGEjshysjk1Cjz\n/0PBsVOEpXXE4/I4wtI68vepyTJQQQghhEPI4AkhbCj111SGfD6IfK7qy3ypyqI+KcQ3j3diZEII\nIdyVDJ4QwknGbx5vkNQB5HOV8ZvHOykiIYQQlYkkdkLYkLm1Y7PysghLDsNjigdhyWGk/prq4MiE\nEEJUBpLYCWFDIUEhJss1NLLyslAosvKyGLJuiCR3QgghbE4SOyFsaHq36fhW9TUo09BQGPZlzb+a\nL82zQgghbE4SOyFsKL55PIt6LiI0KFQ3UjYo1CipK2Su2VYIIYQoK5mgWAgbi28ebzACNiw5jKy8\nLKP9zDXbCiGEEGUlNXZC2Jmp5lnfqr5M7zbdSREJIYSoqCSxE8LOTDXPLuq5SFert023QgXbtjk7\nTCGEEBWATFAshLNs2wbdusGVK+DlBZs3Q4cOzo5KCCGEi5EJioVwBxkZuqTu+nXdNiPD2REJIYRw\nc5LYCeEsMTG6mjpPT902JsbZEQkhhHBzMipWCGfp0EHX/JqRoUvqpBlWCCFEOUliJ4QzdehgnNBt\n2ybJnhBCiDKRxE4IVyIDKoQQQpSD9LETwpXIgAohhBDlIImdEC4ktfFlwkYW4DEJwkYWkNr4srND\nEkII4UYksRPCRaT+msqQ/bPJClIoDbKCFEP2z2bEv0YQlhyGxxQPwpLDSP011dmhCiGEcFHSx04I\nFzF+83jyr+YblOVfzWfBjgUodBOJZ+VlMWTdEACD9WiFEEIIkBo7IVxGdl62yfLCpK5Q/tV8xm8e\n74iQhBBCuBlJ7IRwktRfUw2aWGtUq2HxseaSQCGEEJWbNMUK4QSpv6YyZN0QfdNrVl4WVT2q4uXp\nxZXrV/T7aWhGNXYAIUEhDotVCCGE+5AaOyGcwFR/uqsFVwnwCiA0KBQNjdCgUIZFD8O3qq/Bfr5V\nfZnebbojwxVCCOEmpMZOCCcw15R66uIpTo45aVDWMaQj4zePJzsvm5CgEKZ3m64bOCErVAghhChG\nEjshnCAkKISsvCyT5cXFN483GgGbumoy47+fSnagIiRNY/rfE4l/bLK9whVCCOEmpClWCCeY3m16\nmZtYU39NZcj/kgznu/tfksxvJ4QQwrUTO03T3tc07U9N085qmnZE07RkTdO8LH1cCFcV3zyeRT0X\nGfSnW9RzkUVz043fPJ58rhqU5XNVpkARQgiBppTxiDtXoWnaHUCWUuqCpmnBQBqwRSk12ZLHzYmO\njlY7duywb/BC2InHFA+TI2U1NAomFTghIiGEEPakadpOpVS0Jfu6dI2dUmqvUurCjbsaUAA0svRx\nISoic1Od6Mu3bYMZM3RbM4rPoSfNuEIIUTHYJLHTNK2xpmlTNU37j6ZpOZqmndM0bZemaeM1TfMr\n57lf0zTtPHACaAkkW/O4EBVNif3ztm2Dbt1gwgTd1kRyVziHXlZeFgqlX6ZMkjshhHB/tqqxGwy8\nBBwEpgKvAPuBfwA/aJpWrXBHTdP+qWmaKuEWU/TESqk3lFL+wB3AAuCoNY8LUdGU2D8vIwOuXIHr\n13XbjAyj482tSSt99IQQwv3ZarqTz4AZSqm8ImULNE37HRgPJALv3Sh/Fni+hHPlmSpUSu3TNG03\nsBzoYu3jQlQkpqZAAXRz2nl56ZI6Ly/d/WLMzaEny5QJIYT7s0mNnVJqR7GkrtCnN7ZRRfY9p5Q6\nWcLtqonzFKoKNC7H40JUaKn+hwibFIjHhOuETQok1f+Q0T6l9tETQgjhtuw9eKL+je1xaw/UNC1I\n07QETdNu0XRaAK8D31jyuBCVjb7v3KXjKCDr0vGbfeeKDKgozxx6QgghXJvdpjvRNM0T+DfQFohS\nSu238vhA4HOgNeCFbnDE58CkG9OblPh4SeeW6U5ERRSWHGZyNYtQnzocnnL2ZvPs5s2k+h8yvUyZ\nEEIIl2PNdCf2TOzeRdeXbpxSaoZdLmJdPEOAIQAhISFtsrKMfwCFcGfm5rcDCD0D2UEQkgfTa/Yj\n/vWVDo5OCCFEWTl9HjtN06ahS+oWuUJSB6CUWqSUilZKRdeqVcvZ4Qhhc+b6yGlA1i3olh+7BYao\ntTK1iRBCVFA2T+w0TZuMrq9bCjDM1ucXQphmqu+chmZUh5dfcFmmNhFCiArKpondjaRuErAMeEa5\n8nplQlQwpua3M9c0m52XbdEKFUIIIdyLreaxQ9O0ieiSuuXAYKWULFophIMVn9/O3ICKEJ/aupUp\nigyooEMHR4YqhBDCDmy1pNhzwBQgG9gEPKlp2lNFbvfZ4jpCCOuYndrkaudSV6gQQgjhfmxVY9f2\nxjYEXTNscVuAjTa6lhDCQoW1d0ZTm5yPgKT1Ja5QIYQQwv3YbboTVybz2AmBrm9dRoYuqZNmWCGE\ncFnWTHdisz52Qgg306GDcUInyZ4QQrg1SeyEEDrbtsmACiGEcHP2XitWCOEuMjJkQIUQQrg5SeyE\nEDoxMbqaOk9PGVAhhBBuSppihRA6HTroml+lj50Q/9/e/UfZXdd3Hn++MySQARyLWpXizPgDsF2C\nbh26nbXbhsa2iqW61h9bRmsKOLKCAttjWxggQBgj2kqsrIeOVeCwY6G0Fk2LP3M2VEs47ejSpqtG\nKpuJP9OIMgQmJCH57B/fO+HOnXtn7r1zf9/n45x7bub7/d7v/Q6fE3jx+Xzen4/Utuyxk/S04WG4\n4oqjoW5yxySDNz6PFdcFgzc+zz1mJanFGewkFTW5Y5LRey5g+sk9JGD6yT2M3nNB0XA3uWOSwc2D\nrLhuBYObBw2AktQkBjtJRY1tHWP2yIF5x2aPHGBs69i8Y5M7JhndMsr0zDSJxPTMNKNbRg13ktQE\nBjtJRe2e2V3W8bGtY8wemp13bPbQ7IIAKEmqP4OdpKL6+/rLOl5uAJQk1Z/BThKwcJ7cOaeeQ+/K\n3nnX9K7sZXzdeLaY8aZNsH172QFQklR/BjtJRefJ3f7Pt/P2l72dgb4BgmCgb4CJcycYefxF2Q4V\nV18N69Yx/oL1pQOgJKmhXMdOUsl5cvc+dC+7Lts1/+JNm+btUDHyrWPh3AnGto6xe2Y3/X39jK8b\nZ2TNSON+AUkSYLCTRGXz5CZPO8DYe46w+xnQ/9gRxk87wMiakQVBbnLHpGFPkhrMoVhJZc+Tm9wx\nyejODzLdl0gB032J0Z0fXLC0iUugSFJzGOwkMb5uvKx5cosubZJXUOESKJLUHA7FSjo6RLrU0Omi\nQ7br1mVz71atYvcfPln6OklS3RjsJAEUnSdXqL+vn+mZ6YXH6YOD+44WVPTTxzSPFv28JKl+HIqV\nVLaSQ7ZnXAqrVkFPD6xaxfgZl7oEiiQ1gT12ksq26JDtyb8B27bB2rWMDA/DjlOtipWkBouUUrOf\noeGGhobS1NRUsx9D6nzbtx8NewwPN/tpJKktRcRXU0pD5Vxrj52k+ti+fV5BBVu3Gu4kqc6cYyep\nPrZtm7dDBdu2NfuJJKnjGewk1cfatfMKKli7ttlPJEkdz6FYSfUxPJwNvzrHTpIaxmAnqX6Gh4sH\nOosqJKkuDHaSGsuiCkmqG+fYSWosiyokqW4MdpIay6IKSaobh2IlNZZFFZJUNwY7SY1XrKjCggpJ\nWjaDnaSGmtwxuXAP2cdfZEGFJNWAwU5Sw0zumGR0yyizh2YBmJ6ZZnTLKBx5LSOFBRUGO0mqmMUT\nkhpmbOvY0VA3Z/bQLGMr/96CCkmqAXvsJDXM7pndRY9PP7mHwQ3PZfeTe+g/7hmMn/AwI9hjJ0mV\nMthJapj+vn6mZ6YXHA+C6Sf3AFnIG90yCpDNvbOgQpLK5lCspIYZXzdO78reeceCIJHmHZs9NMul\nn3kXg3e/khUHrmTw7lcy+dfXNvBJJak9GewkNczImhEmzp1goG+AIBjoG1gQ6uY88tRjTPclUsB0\nX2L0X9/H5I7JBj+xJLWXSKn4v1Q72dDQUJqammr2Y0gCBjcPFh2eLWagb4Bdl+2q7wNJUouJiK+m\nlIbKudYeO0lNVWx4tpRSxReSpIzBTlJTFRuefdbqZxW9tr+vP9uhYtOm7F2SNI9VsZKabmTNCCNr\nRo7+XLiQMUDvyl7GX7DeHSokaRH22ElqOcV68SbOnYBvfJ3Bd+5nxVWHGXznfia33sTkjkkGNw+y\n4roVDG4etMBCUldri+KJiFgN7ACel1I6Ie/4bcB5wMG8y9+YUvrcYvezeEJqP5M7Jhm95wJmjxw4\nemwlPURPDwcPP/2vgN6VvUycOzGvB1CS2lknFk9cD5Qqm5tIKZ2Q91o01ElqT2Nbx+aFOoBDHJ4X\n6iC3RdnWsUY+miS1jJYPdhHxCuDVwI3NfhZJzVNJRazVs5K6VU2CXUScFhHXR8QDEbE3IvZFxIMR\nMRYRxy/jvscAHwMuZv5wa76RiPhxRHwj930WhEgdqL+vvy7XSlInqVWP3fnA5cC3yYZN3wvsBG4A\n7s/NkQMgIu6MiLTIa23efd8L/J+U0t+X+N4/BU4Hng28DVgPbKjR7ySphRRb727lipWsomfesV5W\nMr5uvJGPJkkto1bB7q+AU1JKIymlj6SUbkkpvQUYB84ELsi79h3AcxZ5/QNARLwEuIgs3BWVUvpa\nSunfU0pHUkpTZKHuv9Xod5LUQopVyt76+lv5xBlXMTATRIKBmWDijCstnJDUtWoybJkLVcXcBYwB\nZ+Rduw/YV8Ztfwl4LvCtiABYCRwfET8C3lCiFy8BUcGjS2ojhevdAbAGRk7+Ddi2Dd601nXtJHW1\nes9HOyX3vqeKz/4l8KW8n4eB24CXA3sBIuItwOeAx4A1wDXA3VU+q6R2NTy8MNBt356FvbVrDXuS\nukbdgl1E9ABXA08Bn6z08ymlWeDosvMRsTc7nL6bd9m7gFvIevN+ANwBbFrGY0vqBNu3u0OFpK5U\nz+VONpP1sl2TUtq53JullLblL06cO/YrKaWfyq1fd2pK6fqU0qFin4+I0YiYioipvXv3LvdxJLWy\nbduyUHf4cPa+bZs7VEjqCnUJdhGxEbiEbPHgluhBSylNpJSGUkpDz3nOc5r9OJLqae3arKeupwdW\nrWLytAOMbhllemaaRGJ6ZprRLaOGO0kdp+bBLiKuBa4CbiWrapWkxhoezoZfN26ErVsZ+85tzB6a\nnXeJO1RI6kQ1nWOXC3UbgNuBC1M7bEQrqTPlFVTs/kLxnSjcoUJSp6lZj11EXEMW6u4Azk8pHanV\nvSVpOUrtROEOFZI6Ta22FLsYuA7YTbZEyXkR8da816/V4nskqRrj68bpZeW8Y+5QIakT1arH7qzc\nez/ZMOwdBS8nskhqmMIKWICJM650hwpJHS+6cRrc0NBQmpoqtVmGpHY2uWOS0S2j84olelf2MnHu\nBCOPv8hFiyW1nYj4akppqJxr67mOnSQ13NjWsdIVsMPDcMUVR0Pd5I5JBm98HiuuCwZvfJ7Ln0hq\ne/XeUkySGqpUpWvh8ckdk4zecwGzRw4AMP3kHkbvueDo+bGtY+ye2U1/Xz/j68YdtpXUFgx2kjpK\nf18/0zPTRY/nG9s6djTUzZk9coBLP3sp+5/af7TXb24xY8BwJ6nlORQrqaOMrxund2XvvGO9K3sX\nVMCW6tl7ZP8jLmYsqW0Z7CR1lJE1I0ycO8FA3wBBMNA3kBVOFPS2VbqGnYsZS2oHDsVK6jgja0aW\nHDYdXzdetHp29TGreWT/IwuudzFjSe3AHjtJXalUz96HX/NhelccO+/a3hXHupixpLZgj52krlWy\nZ+/Tn2bs0bvZ3Qf9MzD+rN+ycEJSWzDYSVKBkXWXM7Lub+HgQVi1CrZe3uxHkqSyGOwkqdDwMGzd\n6i4VktqOc+wkqRh3qZDUhuyxk6QlLLZLhXPvJLUSe+wkaQmldqlw0WJJrcZgJ0lLKHf/WUlqNoOd\nJC2h1OLELlosqdUY7CRpCYvuP7t9O2zalL1LUpNZPCFJS5grkBjbOsbumd309/Uzvm6ckcdfBOvW\n5a13t9WlUSQ1lcFOkspQdJeKTZuyUHf4cPa+bZvBTlJTORQrSVWaPO0Ag+85wooNMPieI0yedmDp\nD0lSHdljJ0lVmNwxyejODzLblwCY7kuM7vwg7Dg1G6J11wpJTWCwk6QqjG0dY/bQ7Lxjs4dmGbv3\n9xm57jHn3UlqCoOdJFWh1Bp200/uYfCdsLsP+mf2M771Jjjh4YWFF+5YIakOnGMnSVUotYZdANPP\nhBTZ++8d/hTnf/p8pmemSSSmZ6YZ3TJadK/ZyR2TDG4eZMV1KxjcPOh+tJIqZrCTpCoUW9suCFLB\ndYc4zMHDB+cdmz00u2A7sskdk4xuGS0rAEpSKQY7SarCyJoRJs6dYKBvgCAY6BsgLYh1pRUO5Zac\ns+d+tJIq4Bw7SapS4dp2g5sHmZ6ZLuuzhUO57kcrqRbssZOkGik2PLtyxUpW0TPvWC8rs+3I8rgf\nraRaMNhJUo0UG5699fW38okzrmJgJogEAzPBxBlXAswrlDjn1HNK70crSWWKlMqfE9IphoaG0tTU\nVLMfQ1I32b796KLFkyc8zOiW0Xlz6npX9vL2l72dex+6d96yKFBkj1qXSpG6SkR8NaU0VNa1BjtJ\naqxSc/EG+gbYddmuoz/PVcoWBsCJcycMd1IXqSTYORQrSQ1WbqGElbKSKmWwk6QGK7dQotJKWRc4\nlmSwk6QGG183Tu+KY+cd611xbEWVsoUh7l1/9y4XOJZksJOkRhtZM8LE6z/OwHHPJYCB457LxOs/\nvmDeXLHlU3pX9nLOqecsCHG3TN1S82FbewCl9mPxhCS1krzqWYaHmdwxuaAqdmzrWNkLIUNWlFFp\nVW2pwo1ilbsWckj1ZVXsEgx2klrS9u2wbh0cPAirVsHWrTA8vOCyFdetKHv7smz/2qevLTeclarc\nLXY/q3Sl+rIqVpLa0bZtWag7fDh737at6GWl5t4FseDnwgA4e2iWW6ZuWXIuXqkCjWL3s0pXah0G\nO0lqFWvXZj11PT3Z+9q1RS8rNffuoqGL5u16UapXr1g4u/Szl86bT3fS6pPKfmz3s5Vah8FOklrF\n8HA2/LpxY8lhWCi+ddnEuRN89LUfZddluziy4Qi7LtvFQN9A2V/9yP5H5vXiPXbgMVb1rJp3TWGP\n4JyTVp9UdpHFcgoyWq2Yo9WeRwLn2ElS6ysoqChXsQKIYsOzpTxr9bM4YdUJR+finXPqOdz+z7fP\nu9/KFSuJCA4ePnj02Ny8O5i/HVqxz5e7lVoln23EfD93BVEjWTyxBIOdpLZRZkFFKYVVtcUCUilB\ncGTDkUXv9/jBx3lk/yMLPvus1c9i/1P7ywqVhceLhcVyP1tpuCpWdQxL789b7rZwUi0Y7JZgsJPU\nNjZtgquvzgoqenqyYdorrljWLcsNZ+WElEoqdBulsKexVC9esV63xXog8++x2O9dzfIy9VIsuNqj\n2H4Mdksw2ElqG8vssSvHcoYVS/VctZJSQ7aVrAdYGHIrWQ7G4WItV8cEu4i4DTgPOJh3+I0ppc/l\nzj9e8JFjgW+klM5c7L4GO0ltpco5dpWotmenVHhYfczqor2AhcGnkjl/tfxs78resoaj8+X3xBUb\n0q7XcHG1SoXPcns0W0039z52WrB7PKV0SZnX/wtwZ0rpfYtdZ7CT1BEaEPjKUWqeWjk7V1RSkFHO\nZyvREz0cTofLuracnrhKei6LhSsoPrev2kBT7jB5O+wo0u29j10Z7CLiF4D7gf6U0vcXu9ZgJ6nt\nNWCIdrnKDSTVFjAU+2yp+YKlFPbcVVK4Ue7wbDkWC7PFqoGLBZrl/LNopSHkYrq9WKXhwS4iTgPe\nCvw68GLgOODbwN3A5pTSE1Xe9zbgdUAC9gD/C7gxpfRUkWv/DDg5pXTuUvc12Elqe3UoqugElSzx\nMtA3cHSu3WKhslRYK6waXu7yMsWU6lWsdimaSiw37C1n6LTws+W2wXK/t1U1I9i9H7gY+AzwAHAI\nOBt4M/AvwC+mlPbnrr0TeMsitzs7pbQtd+3PA98FfgT8PPAXZEOtVxd8//HA94HfTSl9eqnnNdhJ\nantt0GPXLOUs8VLJMF4lvUXLWV5mOUoFyMIAWGmP5lLf06g1C8vtNe3UIdtmBLsh4KGU0kzB8RuA\nMeDdKaWbc8dOJCtyKGUmpXSoxPecB1yXUjq14Ph6YBPwgmK9eYUMdpI6QovMsWsHy+09Wk5YWM4Q\naSXzAItpRK9iK61ZWOmQbbv07rXMHLuIWEPWY/dnKaWLanC/3wE2ppReUnD8K8BXUkp/VM59DHaS\nOpZhry5qGQAqWT+v2By7SlTbq7jcsFcvhWsEQnXD5lBZYF/OPNBaaKVg9xrgXuD6lNKGKj7/FuBz\nwGPAGuAu4J6U0hV515wOfAM4PaX0UDn3NdhJ6kgOz7aNSoJCuT1+y1lWpR3CXjnDrpXMp1xs15Tl\nbKVXj3DXEsEuInqALwNnAWeklHZWcY/7gDOBlcAPgDuATflDtRHxAeA/pZR+ZYl7jQKjAP39/a+Y\nnm7tBTUlqWIWVHSFUj1Nta5iLXeuYiPWLKxk2LVU0Ue1PZ+VPHe9qnQrCXbH1Pzbn7YZGAaurCbU\nASwV1nLX/EGZ95oAJiDrsavmeSSppa1dm/XUzfXYrV3b7CdSHcyFm3oPA46sGVlwz1f2v7IpaxYW\n+/12z+wu+tyJtGDIdmzrWNXD2ZX0UpZ6pkaqS49dRGwErgImUkrvrPkXLJNDsZI6lnPs1ASNWLOw\nUCWFEo3a07gVeuxqHuwi4lpgA3ArcEFqwRWQDXaSuophTx2okuKHcrdXK3f+YivPsVtR4y++lizU\n3Q5c2IqhTpK6ylxBxdVXZ+/btzf7iaSaGFkzwsS5Ewz0DRAEA30DJYPV+Lpxelf2zjvWu7KXD7/m\nw+y6bBdHNhxh12W7+PBrPlz0uouGLpr3Pbe+/lY+8bpPlPXdjVazOXYRcQ1ZqLsDOD+ldGSJj0iS\n6m3btmzO3eHD2fu2bfbaqWMUmwdY6jpYesi30vmLrRDkCtVqgeKLgZuB3cDVQGGo25NS+uKyv6hG\nHIqV1DVcAkVqe82oij0r995PNgxb6D6gZYKdJHWN4eEszDnHTuoKNQl2KaX1wPpa3EuSVGPDwwsD\nnQUVUkeq5zp2kqRW5PCs1LFqWhUrSWoDxQoqJHUEg50kdZu5HSp6etyhQuowDsVKUrexoELqWAY7\nSepGFlRIHclgJ0myoELqEM6xkyRZUCF1CIOdJMmCCqlDOBQrSbKgQuoQBjtJUsaCCqntGewkScVZ\nUCG1HefYSZKKs6BCajsGO0lScRZUSG3HoVhJUnEWVEhtx2AnSSrNggqprRjsJEnls6BCamnOsZMk\nlc+CCqmlGewkSeWzoEJqaQ7FSpLKZ0GF1NIMdpKkylhQIbUsg50kaXksqJBahnPsJEnLY0GF1DIM\ndpKk5bGgQmoZDsVKkpbHggqpZRjsJEnLV6ygAiyqkBrMYCdJqg+LKqSGc46dJKk+LKqQGs5gJ0mq\nD4sqpIZzKFaSVB8WVUgNZ7CTJNWPu1RIDWWwkyQ1jgUVUl05x06S1DgWVEh1ZbCTJDWOBRVSXTkU\nK0lqHAsqpLoy2EmSGsuCCqluDHaSpOayoEKqGefYSZKay4IKqWYMdpKk5rKgQqoZh2IlSc1lQYVU\nMwY7SVLzWVAh1YTBTpLUeiyokKriHDtJUuuxoEKqisFOktR6LKiQqtLywS4iXhsRX4uIJyLihxHx\n3rxzL4yILRHxSETsiYhNEdHyv5MkaQlzBRUbNzoMK1WgpefYRcSvAxPA7wL3Ab1Af+5cD7AF+Dzw\nRuCngb8FHgVubMbzSpJqyIIKqWItHeyAjcDGlNLW3M+PAf+a+/PpwM8BZ6WUDgDfiYibgA0Y7CSp\n81hQIS2pJsOWEXFaRFwfEQ9ExN6I2BcRD0bEWEQcX+U9jwfOAp4XEd/MDbV+JiJeOHdJwfvcnwcj\n4hnV/zaSpJZkQYW0pFrNRzsfuBz4NnA98F5gJ3ADcH9ErJ67MCLujIi0yGtt7tKfIgtqvw28Gngh\n8EPgUxERuft/GxiPiNW5wHd57rMGO0nqNBZUSEuKlNLybxIxBDyUUpopOH4DMAa8O6V0c+7YicCx\ni9xuJqV0KCL6yObLvSOl9Oe5zz4b2AsMpJR2R8RLgZuAVwA/Bj5ONgx7YkrpiVJfMDQ0lKampqr8\nbSVJTeMcO3WhiPhqSmmonGtrMscupVQqJd1FFuzOyLt2H7CvjHvORMQ0UDJ5ppS+Cbxm7ueIuBj4\np8VCnSSpjVlQIS2q3sUTp+Te91T5+VuASyPiC2Q9dRuBr6aUdgNExJnAw8CTwNlkIfLty3piSVL7\nsKBCmqdua77lliO5GngK+GSVt/kA8Fnga8D3gJOBN+SdfxMwDcwA7ycbtv1iiecZjYipiJjaE06e\nqwAAFMNJREFUu3dvlY8jSWopFlRI89Rkjl3RG0d8BLgEuDKltKkuX1Il59hJUoewx05doOFz7Io8\nwEayUDfRaqFOktRB5naocI6dBNQh2EXEtcBVwK3ARbW+vyRJ81hQIR1V02CXC3UbgNuBC1O9xnkl\nSSrF4Vl1sZoVT0TENWSh7g7g/JTSkVrdW5KksllQoS5Wkx673Ppx1wG7gS8B52WbQxy1p1S1qiRJ\nNTW3Q8Vcj507VKiL1Goo9qzcez/ZMGyh+wCDnSSp/iyoUBer1c4T64H1tbiXJEnLZkGFulS9d56Q\nJKn5LKhQl6jbzhOSJLUMCyrUJQx2kqTON1dQ0dNjQYU6mkOxkqTOZ0GFuoTBTpLUHYoVVIBFFeoo\nBjtJUveyqEIdxjl2kqTuZVGFOozBTpLUvSyqUIdxKFaS1L0sqlCHMdhJkrqbu1SogxjsJEnKZ0GF\n2phz7CRJymdBhdqYwU6SpHwWVKiNORQrSVI+CyrUxgx2kiQVsqBCbcpgJ0nSUiyoUJtwjp0kSUux\noEJtwmAnSdJSLKhQm3AoVpKkpVhQoTZhsJMkqRwWVKgNGOwkSaqGBRVqQc6xkySpGhZUqAUZ7CRJ\nqoYFFWpBDsVKklQNCyrUggx2kiRVy4IKtRiDnSRJtWJBhZrMOXaSJNWKBRVqMoOdJEm1YkGFmsyh\nWEmSasWCCjWZwU6SpFqyoEJNZLCTJKmeLKhQAznHTpKkerKgQg1ksJMkqZ4sqFADORQrSVI9WVCh\nBjLYSZJUbxZUqEEMdpIkNZoFFaoT59hJktRoFlSoTgx2kiQ1mgUVqhOHYiVJajQLKlQnBjtJkprB\nggrVgcFOkqRWYEGFasA5dpIktQILKlQDBjtJklqBBRWqgZYOdhHx/Ij464j4UUQ8EhH3RMQpeeff\nHBFfiYjHI2JXEx9VkqTlmSuo2LjRYVhVrdXn2H2U7BlfCBwGPgZ8Avj13PmfADcDzwUub8YDSpJU\nM8UKKsCiCpWt1YPdi4E/TintA4iITwIfnzuZUvpi7vjrm/N4kiTVmUUVqkBNhmIj4rSIuD4iHoiI\nvRGxLyIejIixiDh+Gbf+EPDGiHhmRJwIvA3YUotnliSpLVhUoQrUao7d+WRDod8GrgfeC+wEbgDu\nj4jVcxdGxJ0RkRZ5rc2771eAZwI/Bh4FTgeurNEzS5LU+iyqUAVqNRT7V8CmlNJM3rFbIuIhYAy4\ngGwuHMA7gEsWudcMQESsAL4EfAo4h2yO3R8A2yLi5SmlQzV6dkmSWpe7VKgCNQl2KaWpEqfuIgt2\nZ+Rduw/YV8ZtTwIGgD9NKT0OEBEfAq4lm3v3zWU8siRJ7cNdKlSmehdPzC1NsqfSD6aUfhQR/wa8\nKyI2kPXYXUpWCbsLICJ6gJW5V0TEcdlH04EaPLskSa3JggqVULd17HKh62rgKeCTVd7mdcCZwHfJ\nwuFvAL+ZUnoyd/5twH7gL4H+3J93lnie0YiYioipvXv3Vvk4kiS1AAsqVEI9e+w2A8PAlSmlomFr\nKSmlrwOvXuT8bcBtZd5rApgAGBoaStU8jyRJLWGuoGKux86CCuXUJdhFxEayAomJlNKmenyHJEld\ny4IKlVDzYBcR1wJXAbcCF9X6/pIkCQsqVFRNg10u1G0AbgcuTCk55ClJUiNYUCFqWDwREdeQhbo7\ngPNTSkdqdW9JkrQECypEjXrsIuJi4DpgN9miwudFRP4le+b2dZUkSXVgQYWo3VDsWbn3frJh2EL3\nAQY7SZLqxYIKUbudJ9YD62txL0mSVCULKrpevXeekCRJzWJBRdep284TkiSpySyo6DoGO0mSOtVc\nQUVPjwUVXcKhWEmSOpUFFV3HYCdJUiezoKKrGOwkSeomFlR0NOfYSZLUTSyo6GgGO0mSuokFFR3N\noVhJkrqJBRUdzWAnSVK3saCiYxnsJEnqdhZUdAzn2EmS1O0sqOgYBjtJkrqdBRUdw6FYSZK6nQUV\nHcNgJ0mSLKjoEAY7SZK0kAUVbck5dpIkaSELKtqSwU6SJC1kQUVbcihWkiQtZEFFWzLYSZKk4iyo\naDsGO0mSVB4LKlqec+wkSVJ5LKhoeQY7SZJUHgsqWp5DsZIkqTwWVLQ8g50kSSpfsYIKsKiiRRjs\nJEnS8lhU0TKcYydJkpbHooqWYbCTJEnLY1FFy3AoVpIkLY9FFS3DYCdJkpbPXSpagsFOkiTVngUV\nTeEcO0mSVHsWVDSFwU6SJNWeBRVN4VCsJEmqPQsqmsJgJ0mS6sOCioYz2EmSpMawoKLunGMnSZIa\nw4KKujPYSZKkxrCgou4cipUkSY1hQUXdGewkSVLjWFBRVwY7SZLUPBZU1JRz7CRJUvNYUFFTBjtJ\nktQ8FlTUVEsHu4h4YURsiYhHImJPRGyKiBV55z8aEd+JiMci4nsRsTkiVjXzmSVJUgXmCio2bnQY\ntgZaNthFRA+wBfgWcDIwBJwDvDfvspuBl6aUngG8LPe6ssGPKkmSlmN4GK64Yn6o274dNm3K3lW2\nVi6eOB34OeCslNIB4DsRcROwAbgRIKX09bzrAzgCnNroB5UkSTVkQUXVatJjFxGnRcT1EfFAROyN\niH0R8WBEjEXE8dXetuB97s+DEfGMvO/+o4h4HPh3sh67zVV+nyRJagUWVFStVkOx5wOXA98Gricb\nLt0J3ADcHxGr5y6MiDsjIi3yWpu7dGfufuMRsToiXpj7DoCjwS6l9P6U0glkvXu3AD+o0e8kSZKa\nwYKKqkVKafk3iRgCHkopzRQcvwEYA96dUro5d+xE4NhFbjeTUjqUu/alwE3AK4AfAx8nG4Y9MaX0\nRJHneBPwrpTS2Ys979DQUJqamir315MkSY3mosVHRcRXU0pD5Vxbkzl2KaVSKekusmB3Rt61+4B9\nZd73m8Br5n6OiIuBfyoW6nJWAqeVc29JktTC3KGiKvUunjgl976nmg9HxJnAw8CTwNlkIfHtuXN9\nwH8F7gFmgDXAVcDnl/fIkiSp5VhQUZa6LXeSW67kauAp4JNV3uZNwDRZcHs/8I6U0hdz5xLwVrLg\nt48s4N0LvLvE84xGxFRETO3du7fKx5EkSU1hQUVZ6tljtxkYBq5MKe2s5gYppavJwmGxc48Br6rg\nXhPABGRz7Kp5HkmS1CRzBRVzPXYWVBRVl2AXERuBS4CJlNKmenyHJEnqInM7VDjHblE1D3YRcS3Z\nXLdbgYtqfX9JktSlLKhYUk2DXS7UbQBuBy5MtVhLRZIkqRgLKhaoWfFERFxDFuruAM5PKR2p1b0l\nSZIWsKBigZr02OXWl7sO2A18CTgvIn8nMPbkVbNKkiQtnwUVC9RqKPas3Hs/2TBsofsAg50kSaod\nCyoWqNXOE+uB9bW4lyRJUtksqJin3jtPSJIkNU6XF1TUbecJSZKkhuvyggqDnSRJ6hxzBRU9PV1Z\nUOFQrCRJ6hxdXlBhsJMkSZ2lWEEFdEVRhcFOkiR1vi4pqnCOnSRJ6nxdUlRhsJMkSZ2vS4oqHIqV\nJEmdr0uKKgx2kiSpO3TBLhUGO0mS1J06sKDCOXaSJKk7dWBBhcFOkiR1pw4sqHAoVpIkdacOLKgw\n2EmSpO7VYQUVBjtJkqQ5bV5Q4Rw7SZKkOW1eUGGwkyRJmtPmBRUOxUqSJM1p84IKg50kSVK+Ni6o\nMNhJkiQtpo0KKpxjJ0mStJg2Kqgw2EmSJC2mjQoqHIqVJElaTBsVVBjsJEmSltImBRUGO0mSpEq1\naEGFc+wkSZIq1aIFFQY7SZKkSrVoQYVDsZIkSZVq0YIKg50kSVI1ihVUNJlDsZIkSR3CYCdJktQh\nDHaSJEkdwmAnSZLUIQx2kiRJHcJgJ0mS1CEMdpIkSR3CYCdJktQhDHaSJEkdwmAnSZLUIQx2kiRJ\nHcJgJ0mS1CEMdpIkSR2iqcEuIt4cEV+JiMcjYleR88dExIcj4scR8WhEfDwijss7f1tEHMx9fu71\n6ob+EpIkSS2i2T12PwFuBsZKnL8SOBtYA5wK/BzwgYJrJlJKJ+S9Ple3p5UkSWphTQ12KaUvppTu\nBKZLXHIh8L6U0vdSSnuBa4H1EdHTqGeUJElqF2UFu4g4LSKuj4gHImJvROyLiAcjYiwijq/Hg0XE\nM4EXAA/mHf4acCIwmHdsJDdU+43c8xxTj+eRJElqdeX22J0PXA58G7geeC+wE7gBuD8iVs9dGBF3\nRkRa5LW2zO88Mff+aN6xRwvO/SlwOvBs4G3AemBDmfeXJEnqKOX2bv0VsCmlNJN37JaIeIhsftwF\nZHPlAN4BXLLIvWYWOZdvX+69D/hh7s/PzD+XUvpa3vVTEbEBuA64uszvkCRJ6hhlBbuU0lSJU3eR\nBbsz8q7dx9OhrGoppUcj4jvAy8l6BwF+PnfvXaU+BsRyv1uSJKkdLbd44pTc+55qPhwRPbnlS1Zm\nP8ZxEXFs3iV/DlwRESdHxHPIiiduSykdzn3+LRHRF5kzgWuAu6v9ZSRJktpZpJSq+2BWmfpl4Czg\njJTSziU+Uuwe64FbCw5Pp5QGc+ePAT5ENn9uBdmQ8CUppf258/cBZ5IFwx8Ad5ANGR8q8l2jwGju\nx9N5uhewnp4N/KgB36PK2C6ty7ZpTbZL67JtWlOt22UgpfScci5cTrD7CNlcuitTSpuqukmHi4ip\nlNJQs59D89kurcu2aU22S+uybVpTM9ulqqHYiNhIFuomDHWSJEmtoeJgFxHXAleRDaFeVOsHkiRJ\nUnUqCna5ULcBuB24MFU7jts9Jpr9ACrKdmldtk1rsl1al23TmprWLmXPsYuIa8jWiLsDWJ9SOlLP\nB5MkSVJlygp2EXEx2QLEu8kW/y0MdXtSSl+s/eNJkiSpXOXuPHFW7r2fbBi20H2AwU6SJKmJyppj\nl1Jan1KKRV5r6/ycbSEiVkTE5RHxzYh4MiK+ExF/EhHHN/vZukFEnBYR10fEAxGxNyL2RcSDETFW\nrA0i4vSIuCcifhIRT0TElyPiV5vx7N0mInoj4uHc/tE3Fzlv2zRQRJwUEX8cEf+W+3fX3oj43xHx\nXwqus10aJCJOiIgrI2JH7t9lP4qI+yNifUREwbW2Sx1ExBURcXfev6t2LXF92e1Qz7xQbo+dynMT\n8B7gb4A/AX429/N/jIhXOS+x7s4HLgY+A0wCh4CzgRuAN0fEL+Ytbv1i4H7gKeADZHsYvwP4fES8\nJqX0pSY8fze5Hii62KZt01gRMQBsA04APg58i2yP7jOBn8m7znZpkIhYAXwW+M9ko2QfAXqB3yFb\nkeJngT/MXWu71M/7gB8DX+PpveqLqqId6pcXUkq+avAC/gPZ3MO/Ljj+brI9bM9r9jN2+gsYAvqK\nHL8h1waX5B37S+Aw8PK8YycA02S7kkSzf59OfZHt+fwU8D9y7XJzwXnbprHt8WXgO8Dzl7jOdmlc\nmwzn/m7cVHB8FfAw8Kjt0pB2eFHen/8V2LXItWW3Q73zwnL3itXTfgcIYHPB8Y8Bs8BbG/5EXSal\nNJVSmily6q7c+xkAua7u3wK2pZQezPv842T7E5/G0/NKVUO5rQg/BnwO+FSR87ZNA0XELwO/BHwg\npfSDiFgZEb1FrrNdGusZuffv5x9MKR0k26bqCbBd6i2l9HA511XRDnXNCwa72jmLLIH/Y/7BlNKT\nwIP4l6uZTsm978m9nwkcC2wvcu0DuXfbqz4uB15KtnNNMbZNY52Te98dEVuA/cATEfGtiMj/j4vt\n0lj/CDwK/EFEvCki+iPipRGxCXgFcG3uOtulNVTaDnXNCwa72jkZ+FFK6UCRc98Dnh0Rqxr8TF0v\n10N0NdnQ3ydzh0/OvX+vyEfmjv1MkXNahoh4IdlamNenlHaVuMy2aazTc+8fA04C3k42V/UgcEdE\n/F7uvO3SQCmln5D1AP2YbIhvGvgG2Rzi304pfSx3qe3SGipth7rmBYsnaqcXKNZIAE/mXXOwMY+j\nnM1k81WuTCntzB2bG2oq1l5PFlyj2rmFbH7Qhxa5xrZprBNz7/uAs3NDfUTEPWRt9b6IuB3bpRke\nJ5vX9RmySfknkQW7T0bE61K2dqzt0hoqbYe65gWDXe3MAj9d4txxedeoQSJiI9mQ30RKaVPeqbl2\nOLbIx2yrOsgN6/0a8MsppUOLXGrbNNb+3PtfzIU6yHqMIuIzwO+S9erZLg0UEWvIwtzlKaVb8o7/\nBVnY+1iuCtN2aQ2VtkNd84JDsbXzfbLu02IN+zNk3a721jVIbl/jq8iWBrio4PTchORiQxRzx4p1\nqasKub8THwLuBX4YES+JiJcAA7lL+nLHnolt02jfzb3/sMi5H+TefwrbpdEuJ/sP/N35B1NKs8Df\nkf3dGcR2aRWVtkNd84LBrnb+ieyf5y/kH4yI44CXA1PNeKhulAt1G8jWf7ow5erI8+wg6wYfLvLx\nX8y92161s5pszbrXAg/lvbblzr819/OF2DaNNjd5+5Qi5+aO/Tu2S6PNhYGeIueOyXu3XVpDpe1Q\n17xgsKudu8jWn7ms4Pg7yMbKJxv+RF0oIq4hC3V3AOenIos85krQtwBrI+JleZ89gSxcPERBtZKW\n5QngTUVe78qd/1zu58/YNg13D9n8urfm/hkDEBHPB14PfCul9G+2S8N9Pfe+Pv9grlf7dcBPANul\nRVTRDnXNC7GwM0PVioiPkM3p+huyYae5laT/AfjVYiFDtRMRFwM3A7vJKmEL/3nvyU04JjcU+I9k\nu1PcBDxG9pdqDfDalNLnG/Xc3SoiBoH/B/zPlNIlecdtmwaKiFHgz4D/C3yCbBHc/w48H/jNlNIX\nctfZLg2S2w3ka2TD4JNk/w05ieyf9yBwcUrpo7lrbZc6iYi38fSUkXeT/d34k9zP0ymlO/Kuragd\n6poXmr2ycye9yLrNf59slekDZGPqHwJOaPazdcMLuI3s/4JKvbYVXP+zwKfJ1ouaBb4CvKrZv0e3\nvMj+A7Vg5wnbpilt8Qay9baeIOvB+wLwStulqW3yYrLpJN8lCwuPAX8PvMF2aVgbbCv3vyeVtkM9\n84I9dpIkSR3COXaSJEkdwmAnSZLUIQx2kiRJHcJgJ0mS1CEMdpIkSR3CYCdJktQhDHaSJEkdwmAn\nSZLUIQx2kiRJHcJgJ0mS1CH+P/joemlAiGYsAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(10,10))\n", "plt.semilogy(np.diag(S), 'r.', basey=2, label=\"True Singular Values\")\n", "plt.semilogy(np.diag(RM), 'go', basey=2, label=\"Modified Gram-Shmidt\")\n", "plt.semilogy(np.diag(RC), 'bx', basey=2, label=\"Classic Gram-Shmidt\")\n", "plt.legend()\n", "rcParams.update({'font.size': 18})" ] }, { "cell_type": "code", "execution_count": 94, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(numpy.float64, numpy.float64, numpy.float64)" ] }, "execution_count": 94, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(A[0,0]), type(RC[0,0]), type(S[0,0])" ] }, { "cell_type": "code", "execution_count": 92, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.2204460492503131e-16" ] }, "execution_count": 92, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eps = np.finfo(np.float64).eps; eps" ] }, { "cell_type": "code", "execution_count": 93, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-52.0, -26.0)" ] }, "execution_count": 93, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.log2(eps), np.log2(np.sqrt(eps))" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "### Ex: Numerical loss of orthogonality" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This example is experiment 3 from section 9 of Trefethen." ] }, { "cell_type": "code", "execution_count": 197, "metadata": { "collapsed": true }, "outputs": [], "source": [ "A = np.array([[0.70000, 0.70711], [0.70001, 0.70711]])" ] }, { "cell_type": "code", "execution_count": 198, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0.7 , 0.70711],\n", " [ 0.70001, 0.70711]])" ] }, "execution_count": 198, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Gram-Schmidt:" ] }, { "cell_type": "code", "execution_count": 96, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Q1, R1 = mgs(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Householder:" ] }, { "cell_type": "code", "execution_count": 105, "metadata": { "collapsed": true }, "outputs": [], "source": [ "R2, V, F = householder_lots(A)\n", "Q2T = np.matmul(block_diag(np.eye(1), F[1]), F[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Numpy's Householder:" ] }, { "cell_type": "code", "execution_count": 106, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Q3, R3 = np.linalg.qr(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check that all the QR factorizations work:" ] }, { "cell_type": "code", "execution_count": 107, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0.7 , 0.7071],\n", " [ 0.7 , 0.7071]])" ] }, "execution_count": 107, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.matmul(Q1, R1)" ] }, { "cell_type": "code", "execution_count": 108, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0.7 , 0.7071],\n", " [ 0.7 , 0.7071]])" ] }, "execution_count": 108, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.matmul(Q2T.T, R2)" ] }, { "cell_type": "code", "execution_count": 109, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0.7 , 0.7071],\n", " [ 0.7 , 0.7071]])" ] }, "execution_count": 109, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.matmul(Q3, R3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check how close Q is to being perfectly orthonormal:" ] }, { "cell_type": "code", "execution_count": 110, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.2547268868202263e-11" ] }, "execution_count": 110, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.linalg.norm(np.matmul(Q1.T, Q1) - np.eye(2)) # Modified Gram-Schmidt" ] }, { "cell_type": "code", "execution_count": 111, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.1110522984689321e-16" ] }, "execution_count": 111, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.linalg.norm(np.matmul(Q2T.T, Q2T) - np.eye(2)) # Our implementation of Householder" ] }, { "cell_type": "code", "execution_count": 112, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.5020189909116529e-16" ] }, "execution_count": 112, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.linalg.norm(np.matmul(Q3.T, Q3) - np.eye(2)) # Numpy (which uses Householder)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "GS (Q1) is less stable than Householder (Q2T, Q3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# End" ] } ], "metadata": { "anaconda-cloud": {}, "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.6.1" } }, "nbformat": 4, "nbformat_minor": 1 }