From 7155c4da3a38516fb3543eaee35a05e7af764cee Mon Sep 17 00:00:00 2001 From: thunder-hammer Date: Thu, 10 Dec 2020 14:03:51 -0500 Subject: [PATCH 1/2] Dallin Skouson Expanded and modified a current write-up to make the explanation more clear and exppand upon the examples in the notebook with python code. --- t19_qr_factorization.ipynb | 352 +++++++++++++++++++++++++++---------- 1 file changed, 256 insertions(+), 96 deletions(-) diff --git a/t19_qr_factorization.ipynb b/t19_qr_factorization.ipynb index d20fc97..1ae5f1c 100644 --- a/t19_qr_factorization.ipynb +++ b/t19_qr_factorization.ipynb @@ -14,7 +14,19 @@ "AUTHOR \n", "\n", "MD SHORIFUL ISLAM; \n", - "EMAIL- ISLAM5@BYU.EDU" + "EMAIL- ISLAM5@BYU.EDU\n", + "\n", + "MODIFIED\n", + "\n", + "DALLIN SKOUSON\n", + "\n", + "changelog \n", + "December 2020:\n", + "Updated code to work with python3 by default\n", + "Added information about other methods of calculating the QR factorization\n", + "Added code example for calcuation using Gram-Schmidt\n", + "Clarified some language\n", + "Added additional exploration to existing code" ] }, { @@ -23,12 +35,19 @@ "source": [ "# Introduction \n", "In linear algebra, a QR factorization of a matrix is a factorization of a matrix A into a product $A = QR$ of an orthogonal matrix Q and an upper triangular matrix R. QR factorization is often used to solve the linear least squares problem and is the basis for a particular eigenvalue algorithm, the QR algorithm.\n", + "\n", + "The Q in the QR factorization will be unitrary (orthogonal when dealing with real numbers). What this means is that $Q^H Q = I$. Additionally it will be shown that $QR$ does yield A.\n", + "\n", "One of the key benefits of using QR factorization over other methods for solving linear least squares is that it is more numerically stable, albeit at the expense of being slower to execute. Hence if we are performing a large quantity of regressions as part of a trading backtest, for instance, we need to consider very extensively whether QR factorization is the best fit.\n", "\n", + "# Methods for Computing the QR factorization\n", + "\n", + "## Calculating the QR matrix using Householder transformations\n", + "\n", "For a square matrix A the QR factorization converts A into the product of an orthogonal matrix Q (i.e. $QTQ=I$) and an upper triangular matrix R. Hence:\n", "\n", "$A=QR$\n", - "The first step is to create the vector x, which is the k-th column of the matrix A, for step k. We define $α=−sgn(xk)(||x||)$. The norm $||⋅||$ used here is the Euclidean norm. Given the first column vector of the identity matrix, I of equal size to A, $e1=(1,0,...,0)T$, we create the vector u:\n", + "The first step is to create the vector x, which is the k-th column of the matrix A, for step k. We define $α=−sgn(xk)(||x||)$ where sgn is the sign operator. The norm $||⋅||$ used here is the Euclidean norm. Given the first column vector of the identity matrix, I of equal size to A, $e1=(1,0,...,0)^T$, we create the vector u:\n", "\n", "$u=x+αe1$\n", "Once we have the vector u, we need to convert it to a unit vector, which we denote as v:\n", @@ -53,12 +72,39 @@ "\n", "$Q=Q^{T}_{1}Q^{T}_{2}...Q^{T}_{t}$\n", "\n", - "This gives $A=QR$, the QR factorization of A." + "This gives $A=QR$, the QR factorization of A.\n", + "\n", + "## Gram-Schmidt and modified Gram-Schmidt algorithms\n", + "\n", + "This method is the most poorly numerically conditioned method of those highlighted here, however it is also the least computationally intensive method of calculating the QR factorization.\n", + "\n", + "The Gram-Schmidt method creates an orthonormal basis to span the columns of A. This orthonormal basis can be used to create an upper triangular R matrix and by extension a unitary Q matrix.\n", + "\n", + "## Givens rotations\n", + "\n", + "Givens rotations are similar to householder transformations, but they zero out individual elements of the matrix one at a time instead of zeroing out portions of columns all at once. Additionally, givens rotations can be implemented using CORDIC rotations which allows the solution to be parallelized and pipleined. Carefully selected rotations can aid in simplifying computational complexity, by ensuring that mulitplications are powers of 2 and that fewer multiplications need to be done.\n" ] }, + { + "source": [ + "# Finding the QR factorization in python\n", + "\n", + "The QR factorization is found using a common library. In order to compare the result for each method the same matrix A is used for each.\n", + "\n", + "A = \n", + "\n", + "\\[\\[12, -51, 4],\n", + "\n", + "\\[6, 167, -68],\n", + "\n", + "\\[-4, 24, -41]]" + ], + "cell_type": "markdown", + "metadata": {} + }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 21, "metadata": {}, "outputs": [], "source": [ @@ -69,7 +115,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 22, "metadata": {}, "outputs": [], "source": [ @@ -78,7 +124,7 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 23, "metadata": {}, "outputs": [], "source": [ @@ -87,7 +133,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 24, "metadata": {}, "outputs": [], "source": [ @@ -98,40 +144,32 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 25, "metadata": {}, "outputs": [ { - "name": "stdout", "output_type": "stream", + "name": "stdout", "text": [ - "A:\n", - "array([[ 12, -51, 4],\n", - " [ 6, 167, -68],\n", - " [ -4, 24, -41]])\n" + "A:\narray([[ 12, -51, 4],\n [ 6, 167, -68],\n [ -4, 24, -41]])\n" ] } ], "source": [ - "#print A matrix\n", - "print (\"A:\")\n", - "#Print the formatted represntation of A on stream, followed by a newline\n", + "print(\"A:\")\n", "pprint.pprint(A)" ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 26, "metadata": {}, "outputs": [ { - "name": "stdout", "output_type": "stream", + "name": "stdout", "text": [ - "Q:\n", - "array([[-0.85714286, 0.39428571, 0.33142857],\n", - " [-0.42857143, -0.90285714, -0.03428571],\n", - " [ 0.28571429, -0.17142857, 0.94285714]])\n" + "Q:\narray([[-0.85714286, 0.39428571, 0.33142857],\n [-0.42857143, -0.90285714, -0.03428571],\n [ 0.28571429, -0.17142857, 0.94285714]])\n" ] } ], @@ -142,17 +180,14 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 27, "metadata": {}, "outputs": [ { - "name": "stdout", "output_type": "stream", + "name": "stdout", "text": [ - "R:\n", - "array([[ -14., -21., 14.],\n", - " [ 0., -175., 70.],\n", - " [ 0., 0., -35.]])\n" + "R:\narray([[ -14., -21., 14.],\n [ 0., -175., 70.],\n [ 0., 0., -35.]])\n" ] } ], @@ -162,47 +197,122 @@ ] }, { + "source": [ + "$QR$ is in fact a factorization of $A$ as demonstrated by multiplying the matricies togehter" + ], "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 28, "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "QR:\narray([[ 12., -51., 4.],\n [ 6., 167., -68.],\n [ -4., 24., -41.]])\n" + ] + } + ], "source": [ - "# Python implementation of QR factorization without any library " + "# demonstrate that QR yields A\n", + "def mult_matrix(X, Y):\n", + " \"\"\"Multiply square matrices X and Y of same dimension\"\"\"\n", + " # Nested list comprehension to calculate matrix multiplication\n", + " return [[sum(a*b for a,b in zip(X_row,Y_col)) for Y_col in zip(*Y)] for X_row in X]\n", + " \n", + "print(\"QR:\")\n", + "pprint.pprint(scipy.array(mult_matrix(Q,R)))" ] }, + { + "source": [ + "Additionally $Q$ is a unitary, and in this case orthogonal matrix as shown by multiplying $Q^T Q$ and $Q Q^T$ both of which should yield $I$" + ], + "cell_type": "markdown", + "metadata": {} + }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 29, "metadata": {}, - "outputs": [], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Q tranpose(Q)\n", + "array([[ 1.00000000e+00, 8.84708973e-17, 0.00000000e+00],\n", + " [ 8.84708973e-17, 1.00000000e+00, -3.46944695e-17],\n", + " [ 0.00000000e+00, -3.46944695e-17, 1.00000000e+00]])\n", + "tranpose(Q) Q\n", + "array([[ 1.00000000e+00, -2.08166817e-17, -5.55111512e-17],\n", + " [-2.08166817e-17, 1.00000000e+00, 2.77555756e-17],\n", + " [-5.55111512e-17, 2.77555756e-17, 1.00000000e+00]])\n" + ] + } + ], "source": [ - "from math import sqrt" + "#show that Q is unitary\n", + "def mult_matrix(X, Y):\n", + " \"\"\"Multiply square matrices X and Y of same dimension\"\"\"\n", + " # Nested list comprehension to calculate matrix multiplication\n", + " return [[sum(a*b for a,b in zip(X_row,Y_col)) for Y_col in zip(*Y)] for X_row in X]\n", + "\n", + "def trans_matrix(M):\n", + " \"\"\"Take the transpose of a matrix.\"\"\"\n", + " n = len(M)\n", + " return [[ M[i][j] for i in range(n)] for j in range(n)]\n", + "\n", + "print(\"Q tranpose(Q)\")\n", + "pprint.pprint(scipy.array(mult_matrix(Q,trans_matrix(Q))))\n", + "print(\"tranpose(Q) Q\")\n", + "pprint.pprint(scipy.array(mult_matrix(trans_matrix(Q),Q)))" + ] + }, + { + "source": [ + "Note that although the values are not exactly zero off the diagonal this can be explained by floating point math error." + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# QR factorization using householder transformations\n", + "\n", + "The QR factorization is computed using householder transformations." ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 30, "metadata": {}, "outputs": [], "source": [ + "from math import sqrt\n", "from pprint import pprint" ] }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 31, "metadata": {}, "outputs": [], "source": [ - "def mult_matrix(M, N):\n", - " \"\"\"Multiply square matrices of same dimension M and N\"\"\"\n", - " # Converts N into a list of tuples of columns \n", - " tuple_N = zip(*N)\n", - " # Nested list comprehension to calculate matrix multiplication \n", - " return [[sum(el_m * el_n for el_m, el_n in zip(row_m, col_n)) for col_n in tuple_N] for row_m in M]" + "def mult_matrix(X, Y):\n", + " \"\"\"Multiply square matrices X and Y of same dimension\"\"\"\n", + " # Nested list comprehension to calculate matrix multiplication\n", + " return [[sum(a*b for a,b in zip(X_row,Y_col)) for Y_col in zip(*Y)] for X_row in X]" ] }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 32, "metadata": {}, "outputs": [], "source": [ @@ -214,7 +324,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 33, "metadata": {}, "outputs": [], "source": [ @@ -225,7 +335,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 34, "metadata": {}, "outputs": [], "source": [ @@ -240,7 +350,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 35, "metadata": {}, "outputs": [], "source": [ @@ -252,29 +362,29 @@ "\n", " # Set R equal to A, and create Q as a zero matrix of the same size\n", " R = A\n", - " Q = [[0.0] * n for i in xrange(n)]\n", + " Q = [[0.0] * n for i in range(n)]\n", "\n", " # The Householder procedure\n", " for k in range(n-1): # We don't perform the procedure on a 1x1 matrix, so we reduce the index by 1\n", " # Create identity matrix of same size as A \n", - " I = [[float(i == j) for i in xrange(n)] for j in xrange(n)]\n", + " I = [[float(i == j) for i in range(n)] for j in range(n)]\n", "\n", " # Create the vectors x, e and the scalar alpha\n", - " # Python does not have a sgn function, so we use cmp instead\n", + " # find the sign of the values. used -cmp() in python2\n", " x = [row[k] for row in R[k:]]\n", " e = [row[k] for row in I[k:]]\n", - " alpha = -cmp(x[0],0) * norm(x)\n", + " alpha = -((x[0] > 0) - (x[0] < 0)) * norm(x)#-cmp(x[0],0) * norm(x)\n", "\n", " # Using anonymous functions, we create u and v\n", - " u = map(lambda p,q: p + alpha * q, x, e)\n", + " u = list(map(lambda p,q: p + alpha * q, x, e))\n", " norm_u = norm(u)\n", - " v = map(lambda p: p/norm_u, u)\n", + " v = list(map(lambda p: p/norm_u, u))\n", "\n", " # Create the Q minor matrix\n", - " Q_min = [ [float(i==j) - 2.0 * v[i] * v[j] for i in xrange(n-k)] for j in xrange(n-k) ]\n", + " Q_min = [ [float(i==j) - 2.0 * v[i] * v[j] for i in range(n-k)] for j in range(n-k) ]\n", "\n", " # \"Pad out\" the Q minor matrix with elements from the identity\n", - " Q_t = [[ Q_i(Q_min,i,j,k) for i in xrange(n)] for j in xrange(n)]\n", + " Q_t = [[ Q_i(Q_min,i,j,k) for i in range(n)] for j in range(n)]\n", "\n", " # If this is the first run through, right multiply by A,\n", " # else right multiply by Q\n", @@ -292,24 +402,11 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 36, "metadata": { "collapsed": true }, - "outputs": [ - { - "ename": "NameError", - "evalue": "name 'xrange' is not defined", - "output_type": "error", - "traceback": [ - "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)", - "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[0mA\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m12\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m-\u001b[0m\u001b[1;36m51\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m4\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;36m6\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m167\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m-\u001b[0m\u001b[1;36m68\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m4\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m24\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m-\u001b[0m\u001b[1;36m41\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mQ\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mR\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mhouseholder\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", - "\u001b[1;32m\u001b[0m in \u001b[0;36mhouseholder\u001b[1;34m(A)\u001b[0m\n\u001b[0;32m 7\u001b[0m \u001b[1;31m# Set R equal to A, and create Q as a zero matrix of the same size\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 8\u001b[0m \u001b[0mR\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mA\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 9\u001b[1;33m \u001b[0mQ\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0.0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mn\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mxrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 10\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 11\u001b[0m \u001b[1;31m# The Householder procedure\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;31mNameError\u001b[0m: name 'xrange' is not defined" - ] - } - ], + "outputs": [], "source": [ "#creating a matrix \n", "A = [[12, -51, 4], [6, 167, -68], [-4, 24, -41]]\n", @@ -319,15 +416,14 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 37, "metadata": {}, "outputs": [ { - "name": "stdout", "output_type": "stream", + "name": "stdout", "text": [ - "A:\n", - "[[12, -51, 4], [6, 167, -68], [-4, 24, -41]]\n" + "A:\n[[12, -51, 4], [6, 167, -68], [-4, 24, -41]]\n" ] } ], @@ -338,17 +434,14 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 38, "metadata": {}, "outputs": [ { - "name": "stdout", "output_type": "stream", + "name": "stdout", "text": [ - "Q:\n", - "array([[-0.85714286, 0.39428571, 0.33142857],\n", - " [-0.42857143, -0.90285714, -0.03428571],\n", - " [ 0.28571429, -0.17142857, 0.94285714]])\n" + "Q:\n[[0.8571428571428571, 0.39428571428571435, -0.33142857142857135],\n [0.4285714285714286, -0.9028571428571429, 0.034285714285714114],\n [-0.28571428571428575, -0.17142857142857126, -0.942857142857143]]\n" ] } ], @@ -359,17 +452,14 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 39, "metadata": {}, "outputs": [ { - "name": "stdout", "output_type": "stream", + "name": "stdout", "text": [ - "R:\n", - "array([[ -14., -21., 14.],\n", - " [ 0., -175., 70.],\n", - " [ 0., 0., -35.]])\n" + "R:\n[[13.999999999999998, 21.00000000000001, -14.000000000000004],\n [-5.506706202140776e-16, -175.00000000000003, 70.0],\n [3.0198066269804245e-16, -3.552713678800501e-14, 35.000000000000014]]\n" ] } ], @@ -379,10 +469,86 @@ ] }, { + "source": [ + "This result is comparable to the result achieved with scipy, being off by some rounding error. Casting the resulting matricies to scipy.array will aid in comparison because it will print with fewer significant digits." + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Q:\n" + ] + }, + { + "output_type": "error", + "ename": "AttributeError", + "evalue": "'function' object has no attribute 'pprint'", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Q:\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mpprint\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mQ\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 4\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"R:\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mpprint\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mR\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mAttributeError\u001b[0m: 'function' object has no attribute 'pprint'" + ] + } + ], + "source": [ + "import scipy\n", + "print(\"Q:\")\n", + "pprint(scipy.array(Q))\n", + "print(\"R:\")\n", + "pprint(scipy.array(R))" + ] + }, + { + "source": [ + "# QR factorization with Gram-Schmidt" + ], "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "# Application of QR Factorization " + "import numpy as np\n", + "def gram_schmidt(A):\n", + " m = A.shape[0]\n", + " n = A.shape[1]\n", + " R=np.zeros([m,n])\n", + " Q=np.zeros([m,m])\n", + " R[0,0]=np.linalg.norm(A[:,0])\n", + " if R[0, 0] == 0:\n", + " print (\"Could not compute the QR factorization\")\n", + " else:\n", + " Q[:,0]=A[:,0]/R[0,0]\n", + " for i in range(1,n):\n", + " Q[:,i]=A[:,i]\n", + " for j in range(i):\n", + " R[j,i]=np.dot(Q[:,j].T,Q[:,i])\n", + " Q[:,i]-=R[j,i]*Q[:,j]\n", + " R[i,i]=np.linalg.norm(Q[:,i])\n", + " if R[i,i]==0:\n", + " print(\"Could not compute the QR factorization\")\n", + " else:\n", + " Q[:, i] = Q[:, i] / R[i, i]\n", + " return Q,R\n", + "\n", + "A = np.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]])\n", + "\n", + "Q, R = gram_schmidt(A)\n", + "\n", + "print(\"Q = \", Q)\n", + "print(\"R = \", R)" ] }, { @@ -391,8 +557,8 @@ "source": [ "# Linear Equations\n", "\n", - "In numerical analysis, different factorization are used to implement efficient matrix algorithms. For instance, when solving a system of linear equations , the matrix A can be factorized via the LU factorization. The LU factorization factorizes a matrix into a lower triangular matrix L and an upper triangular matrix U. The systems and require fewer additions and multiplications to solve, compared with the original system , though one might require significantly more digits in inexact arithmetic such as floating point.\n", - "Similarly, the QR factorization expresses A as QR with Q a unitary matrix and R an upper triangular matrix. The system $Q(Rx) = b$ is solved by $Rx = Q^{T}b = c$, and the system $Rx = c$ is solved by 'back substitution'. The number of additions and multiplications required is about twice that of using the LU solver, but no more digits are required in inexact arithmetic because the QR factorization is numerically stable." + "In numerical analysis, different factorization are used to implement efficient matrix algorithms. For instance, when solving a system of linear equations , the matrix A can be factorized via the LU factorization. The LU factorization factors a matrix into a lower triangular matrix L and an upper triangular matrix U. Solving the system will require fewer additions and multiplications in the LU factorization, though one might require significantly more digits in inexact arithmetic such as floating point.\n", + "Similarly, the QR factorization expresses A as QR with Q a unitary matrix and R an upper triangular matrix. The system $Q(Rx) = b$ is solved by $Rx = Q^{T}b = c$, and the system $Rx = c$ is solved by 'back substitution'. The number of additions and multiplications required is about twice that of using the LU solver, but because the QR factorization is numerically stable, no additional significant digits are required in inexact arithmetic." ] }, { @@ -410,15 +576,9 @@ "metadata": {}, "source": [ "# Conclusion \n", - "The result of the python generated code and the SciPy code should be same for the accurate python code. I got the same result for the python code and SciPy. This code is valid only both for the Python2 and Python3. It may show some error in python3 because python3 prefers range instead of xrange command. xrange need to be change in range for python3 which is a valid command to avoid any kind of error. " + "\n", + "Householder transformations and Gram-Schmidt are an appropriate methods for calculating the QR factorization of a matrix. Gram-Schmidt is much simpler compuatationally but it is less robust to poorly conditioned matricies. All methods yield similar results as well" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { @@ -437,9 +597,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", - "version": "2.7.15" + "version": "3.8.5-final" } }, "nbformat": 4, "nbformat_minor": 2 -} +} \ No newline at end of file From 9e987205748628c73157a06b3b7e79d330468d16 Mon Sep 17 00:00:00 2001 From: thunder-hammer Date: Thu, 10 Dec 2020 14:06:56 -0500 Subject: [PATCH 2/2] forgot a save --- t19_qr_factorization.ipynb | 74 +++++++++++++++++--------------------- 1 file changed, 32 insertions(+), 42 deletions(-) diff --git a/t19_qr_factorization.ipynb b/t19_qr_factorization.ipynb index 1ae5f1c..170a1df 100644 --- a/t19_qr_factorization.ipynb +++ b/t19_qr_factorization.ipynb @@ -104,7 +104,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -115,7 +115,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -124,7 +124,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -133,7 +133,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -144,7 +144,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 5, "metadata": {}, "outputs": [ { @@ -162,7 +162,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 6, "metadata": {}, "outputs": [ { @@ -180,7 +180,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 7, "metadata": {}, "outputs": [ { @@ -205,7 +205,7 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 8, "metadata": {}, "outputs": [ { @@ -236,21 +236,14 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 9, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ - "Q tranpose(Q)\n", - "array([[ 1.00000000e+00, 8.84708973e-17, 0.00000000e+00],\n", - " [ 8.84708973e-17, 1.00000000e+00, -3.46944695e-17],\n", - " [ 0.00000000e+00, -3.46944695e-17, 1.00000000e+00]])\n", - "tranpose(Q) Q\n", - "array([[ 1.00000000e+00, -2.08166817e-17, -5.55111512e-17],\n", - " [-2.08166817e-17, 1.00000000e+00, 2.77555756e-17],\n", - " [-5.55111512e-17, 2.77555756e-17, 1.00000000e+00]])\n" + "Q tranpose(Q)\narray([[ 1.00000000e+00, 8.84708973e-17, 0.00000000e+00],\n [ 8.84708973e-17, 1.00000000e+00, -3.46944695e-17],\n [ 0.00000000e+00, -3.46944695e-17, 1.00000000e+00]])\ntranpose(Q) Q\narray([[ 1.00000000e+00, -2.08166817e-17, -5.55111512e-17],\n [-2.08166817e-17, 1.00000000e+00, 2.77555756e-17],\n [-5.55111512e-17, 2.77555756e-17, 1.00000000e+00]])\n" ] } ], @@ -290,7 +283,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ @@ -300,7 +293,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ @@ -312,7 +305,7 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ @@ -324,7 +317,7 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 13, "metadata": {}, "outputs": [], "source": [ @@ -335,7 +328,7 @@ }, { "cell_type": "code", - "execution_count": 34, + "execution_count": 14, "metadata": {}, "outputs": [], "source": [ @@ -350,7 +343,7 @@ }, { "cell_type": "code", - "execution_count": 35, + "execution_count": 15, "metadata": {}, "outputs": [], "source": [ @@ -402,7 +395,7 @@ }, { "cell_type": "code", - "execution_count": 36, + "execution_count": 16, "metadata": { "collapsed": true }, @@ -416,7 +409,7 @@ }, { "cell_type": "code", - "execution_count": 37, + "execution_count": 17, "metadata": {}, "outputs": [ { @@ -434,7 +427,7 @@ }, { "cell_type": "code", - "execution_count": 38, + "execution_count": 18, "metadata": {}, "outputs": [ { @@ -452,7 +445,7 @@ }, { "cell_type": "code", - "execution_count": 39, + "execution_count": 19, "metadata": {}, "outputs": [ { @@ -477,25 +470,14 @@ }, { "cell_type": "code", - "execution_count": 40, + "execution_count": 20, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ - "Q:\n" - ] - }, - { - "output_type": "error", - "ename": "AttributeError", - "evalue": "'function' object has no attribute 'pprint'", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Q:\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mpprint\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mQ\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 4\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"R:\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mpprint\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mR\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mAttributeError\u001b[0m: 'function' object has no attribute 'pprint'" + "Q:\narray([[ 0.85714286, 0.39428571, -0.33142857],\n [ 0.42857143, -0.90285714, 0.03428571],\n [-0.28571429, -0.17142857, -0.94285714]])\nR:\narray([[ 1.40000000e+01, 2.10000000e+01, -1.40000000e+01],\n [-5.50670620e-16, -1.75000000e+02, 7.00000000e+01],\n [ 3.01980663e-16, -3.55271368e-14, 3.50000000e+01]])\n" ] } ], @@ -516,9 +498,17 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 21, "metadata": {}, - "outputs": [], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Q = [[ 0.85714286 -0.39428571 -0.33142857]\n [ 0.42857143 0.90285714 0.03428571]\n [-0.28571429 0.17142857 -0.94285714]]\nR = [[ 14. 21. -14.]\n [ 0. 175. -70.]\n [ 0. 0. 35.]]\n" + ] + } + ], "source": [ "import numpy as np\n", "def gram_schmidt(A):\n",