Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
TB Lecture 7: QR factorization
{
"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
}
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.