Last active
September 22, 2016 15:28
-
-
Save tobydriscoll/4d1c8856da0c718e1d99d8171e1babec to your computer and use it in GitHub Desktop.
TB Lecture 7: QR factorization
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"# Lecture 7: The QR factorization" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## Gram-Schmidt orthogonalization\n", | |
"\n", | |
"The classical way to prove the existence of a QR factorization is to construct it via the Gram-Schmidt process. Let's look at it for a small example. " | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 1, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"6x4 Array{Float64,2}:\n", | |
" 0.0 2.0 9.0 1.0\n", | |
" 7.0 8.0 1.0 4.0\n", | |
" 8.0 9.0 9.0 1.0\n", | |
" 1.0 5.0 5.0 1.0\n", | |
" 7.0 10.0 4.0 2.0\n", | |
" 0.0 5.0 4.0 4.0" | |
] | |
}, | |
"execution_count": 1, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"A = round(10*rand(6,4))" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"First we want an orthonormal basis for the span of the first column. That's pretty easy." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 2, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [], | |
"source": [ | |
"(Q,R) = (zeros(6,4),zeros(4,4));\n", | |
"R[1,1] = norm(A[:,1]);\n", | |
"Q[:,1] = A[:,1]/R[1,1];" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Next we add another orthonormal vector to become a basis for the span of the first _two_ columns. We can do that via orthogonal projection..." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 3, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"v' * Q[:,1] = [-4.6629367034256575e-15]\n" | |
] | |
} | |
], | |
"source": [ | |
"P = Q[:,1]*Q[:,1]'; \n", | |
"v = A[:,2] - P*A[:,2];\n", | |
"@show v'*Q[:,1];" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"However, we'll organize the process in a particular way. \n", | |
"\n", | |
"$$\n", | |
"r_{22}q_2 = v = a_2 - (QQ^*)a_2 = a_2 - Q(Q^*a_2) = a_2 - q_1(r_{12})\n", | |
"$$\n", | |
"\n", | |
"Doing so allows us to rearrange to\n", | |
"\n", | |
"$$\n", | |
"a_2 = q_1r_{12} + q_2 r_{22}, \n", | |
"$$\n", | |
"\n", | |
"which is the second column of $A=QR$. " | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 4, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"R = [12.767145334803704 15.900187134755534 0.0 0.0\n", | |
" 0.0 6.7958847164850145 0.0 0.0\n", | |
" 0.0 0.0 0.0 0.0\n", | |
" 0.0 0.0 0.0 0.0]\n", | |
"Q = [0.0 0.29429575153747534 0.0 0.0\n", | |
" 0.5482823149915702 -0.10562148138001427 0.0 0.0\n", | |
" 0.6266083599903659 -0.14173138954412173 0.0 0.0\n", | |
" 0.07832604499879574 0.5524815949108433 0.0 0.0\n", | |
" 0.5482823149915702 0.18867427015746108 0.0 0.0\n", | |
" 0.0 0.7357393788436885 0.0 0.0]\n" | |
] | |
} | |
], | |
"source": [ | |
"R[1,2] = dot(Q[:,1],A[:,2]);\n", | |
"v = A[:,2] - Q[:,1]*R[1,2];\n", | |
"\n", | |
"R[2,2] = norm(v);\n", | |
"@show R;\n", | |
"Q[:,2] = v/R[2,2];\n", | |
"@show Q;" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Now the relevant projector is" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 5, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"v' * Q[:,1:2] = [-2.852267617604043e-15 3.2753572043548736e-15]\n" | |
] | |
} | |
], | |
"source": [ | |
"P = Q[:,1:2]*Q[:,1:2]';\n", | |
"v = A[:,2] - P*A[:,2];\n", | |
"@show v'*Q[:,1:2];" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"But we break it into steps as before." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 6, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"R = [12.767145334803704 15.900187134755534 8.772517039865122 0.0\n", | |
" 0.0 6.7958847164850145 7.727520347118983 0.0\n", | |
" 0.0 0.0 9.128437657679335 0.0\n", | |
" 0.0 0.0 0.0 0.0]\n" | |
] | |
} | |
], | |
"source": [ | |
"(R[1,3],R[2,3]) = dot(Q[:,1],A[:,3]),dot(Q[:,2],A[:,3]);\n", | |
"v = A[:,3] - Q[:,1:2]*R[1:2,3];\n", | |
"R[3,3] = norm(v);\n", | |
"@show R;\n", | |
"Q[:,3] = v/R[3,3];" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"$Q$ is unitary and $R$ is upper triangular by construction, and we have" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 7, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"norm(Q * R[:,1:3] - A[:,1:3]) = 1.275549143317629e-15\n" | |
] | |
} | |
], | |
"source": [ | |
"@show norm( Q*R[:,1:3] - A[:,1:3] );" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Rather than finish things manually, let's start over and automate everything in a loop." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 8, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [], | |
"source": [ | |
"(m,n) = size(A);\n", | |
"(Q,R) = (zeros(m,n),zeros(n,n));\n", | |
"R[1,1] = norm(A[:,1]);\n", | |
"Q[:,1] = A[:,1]/R[1,1];\n", | |
"for j = 2:n\n", | |
" R[1:j-1,j] = Q[:,1:j-1]'*A[:,j];\n", | |
" v = A[:,j] - Q[:,1:j-1]*R[1:j-1,j];\n", | |
" R[j,j] = norm(v);\n", | |
" Q[:,j] = v/R[j,j];\n", | |
"end" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 9, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"Q = [0.0 0.29429575153747534 0.7367989840260739 0.3107599849234479\n", | |
" 0.5482823149915702 -0.10562148138001454 -0.32794481561090993 0.6096183641478815\n", | |
" 0.6266083599903659 -0.14173138954412198 0.5037334814158353 -0.11637297026426416\n", | |
" 0.07832604499879574 0.5524815949108433 0.004773069733776502 -0.44844297969083946\n", | |
" 0.5482823149915702 0.18867427015746083 -0.24843245882629786 -0.41255740103288735\n", | |
" 0.0 0.7357393788436885 -0.184636307262656 0.3833357516263034]\n" | |
] | |
} | |
], | |
"source": [ | |
" \n", | |
"@show Q;" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 10, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"2.5207720726002076e-15" | |
] | |
}, | |
"execution_count": 10, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"norm(Q'*Q-eye(n))" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 11, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"1.5036864803660879e-15" | |
] | |
}, | |
"execution_count": 11, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"norm(A-Q*R)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"This only produces the \"thin\" version of the QR factorization. Conceptually at least, the full version just tacks on an orthonormal basis for whatever remains of $\\mathbb{C}^m$." | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## Linear system solution\n", | |
"\n", | |
"We can apply QR to solve a linear system $Ax=b$ for $m\\times m$ matrix $A$. We have $QRx=b$, so $Rx=Q^*b$ is an equivalent triangular system, solvable by back substitution. " | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 12, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"x = [0.40145622419481397,0.6999528823342192,1.4627482274682255,-0.4412557347859269,1.155606318206312,1.1186980122017582,-2.074662732712543,-0.17133966158239913,-0.9745719879262823]\n" | |
] | |
} | |
], | |
"source": [ | |
"A = round(10*rand(9,9)); b = collect(1:9);\n", | |
"m = 9;\n", | |
"(Q,R) = qr(A);\n", | |
"z = Q'*b;\n", | |
"x = zeros(m);\n", | |
"x[m] = z[m]/R[m,m];\n", | |
"for i = m-1:-1:1\n", | |
" y = z[i] - R[i,i+1:m]*x[i+1:m];\n", | |
" x[i] = y[1] / R[i,i];\n", | |
"end\n", | |
"@show x;" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 13, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"1.0880185641326534e-14" | |
] | |
}, | |
"execution_count": 13, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"norm(A*x-b)" | |
] | |
} | |
], | |
"metadata": { | |
"kernelspec": { | |
"display_name": "Julia 0.4.3", | |
"language": "julia", | |
"name": "julia-0.4" | |
}, | |
"language_info": { | |
"file_extension": ".jl", | |
"mimetype": "application/julia", | |
"name": "julia", | |
"version": "0.4.3" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 1 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment