Last active
November 25, 2024 11:41
-
-
Save tobydriscoll/8d87997704e9fd400e96ea860d9f6a34 to your computer and use it in GitHub Desktop.
Trefethen & Bau lecture notes in MATLAB
This file contains hidden or 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": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "A =\n", | |
| "\n", | |
| " 35 1 6 26\n", | |
| " 3 32 7 21\n", | |
| " 31 9 2 22\n", | |
| " 8 28 33 17\n", | |
| " 30 5 34 12\n", | |
| " 4 36 29 13\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "A = magic(6); A = A(:,1: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": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "R =\n", | |
| "\n", | |
| " 56.3471\n", | |
| "\n", | |
| "\n", | |
| "Q =\n", | |
| "\n", | |
| " 0.6211\n", | |
| " 0.0532\n", | |
| " 0.5502\n", | |
| " 0.1420\n", | |
| " 0.5324\n", | |
| " 0.0710\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "R = norm(A(:,1)), Q = A(:,1)/R" | |
| ] | |
| }, | |
| { | |
| "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": [ | |
| "innerprod =\n", | |
| "\n", | |
| " -2.2204e-15\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "P = Q*Q'; \n", | |
| "v = A(:,2) - P*A(:,2);\n", | |
| "innerprod = v'*Q" | |
| ] | |
| }, | |
| { | |
| "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 =\n", | |
| "\n", | |
| " 56.3471 16.4693\n", | |
| " 0 54.2196\n", | |
| "\n", | |
| "\n", | |
| "Q =\n", | |
| "\n", | |
| " 0.6211 -0.1702\n", | |
| " 0.0532 0.5740\n", | |
| " 0.5502 -0.0011\n", | |
| " 0.1420 0.4733\n", | |
| " 0.5324 -0.0695\n", | |
| " 0.0710 0.6424\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "R(1,2) = Q(:,1)'*A(:,2);\n", | |
| "v = A(:,2) - Q(:,1)*R(1,2);\n", | |
| "R(2,2) = norm(v)\n", | |
| "Q(:,2) = v/R(2,2)" | |
| ] | |
| }, | |
| { | |
| "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": [ | |
| "innerprod =\n", | |
| "\n", | |
| " 1.0e-14 *\n", | |
| "\n", | |
| " 0.2228 -0.2821\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "P = Q(:,1:2)*Q(:,1:2)';\n", | |
| "v = A(:,2) - P*A(:,2);\n", | |
| "innerprod = v'*Q" | |
| ] | |
| }, | |
| { | |
| "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 =\n", | |
| "\n", | |
| " 56.3471 16.4693 30.0459\n", | |
| " 0 54.2196 34.8797\n", | |
| " 0 0 32.4907\n", | |
| "\n", | |
| "\n", | |
| "Q =\n", | |
| "\n", | |
| " 0.6211 -0.1702 -0.2070\n", | |
| " 0.0532 0.5740 -0.4500\n", | |
| " 0.5502 -0.0011 -0.4460\n", | |
| " 0.1420 0.4733 0.3763\n", | |
| " 0.5324 -0.0695 0.6287\n", | |
| " 0.0710 0.6424 0.1373\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "R(1:2,3) = Q(:,1:2)'*A(:,3);\n", | |
| "v = A(:,3) - Q(:,1:2)*R(1:2,3);\n", | |
| "R(3,3) = norm(v)\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": [ | |
| "ans =\n", | |
| "\n", | |
| " 3.5527e-15\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "norm( Q*R - 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": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "R = norm(A(:,1));\n", | |
| "Q = A(:,1)/R(1,1);\n", | |
| "for j = 2:size(A,2)\n", | |
| " R(1:j-1,j) = Q'*A(:,j);\n", | |
| " v = A(:,j) - Q*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 =\n", | |
| "\n", | |
| " 0.6211 -0.1702 -0.2070 0.4998\n", | |
| " 0.0532 0.5740 -0.4500 0.2106\n", | |
| " 0.5502 -0.0011 -0.4460 -0.4537\n", | |
| " 0.1420 0.4733 0.3763 0.5034\n", | |
| " 0.5324 -0.0695 0.6287 -0.2096\n", | |
| " 0.0710 0.6424 0.1373 -0.4501\n", | |
| "\n", | |
| "\n", | |
| "R =\n", | |
| "\n", | |
| " 56.3471 16.4693 30.0459 39.0969\n", | |
| " 0 54.2196 34.8797 23.1669\n", | |
| " 0 0 32.4907 -8.9182\n", | |
| " 0 0 0 7.6283\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "Q,R" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 10, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "ans =\n", | |
| "\n", | |
| " 1.2510e-15\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "norm(Q'*Q-eye(4))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 11, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "ans =\n", | |
| "\n", | |
| " 7.9441e-15\n" | |
| ] | |
| } | |
| ], | |
| "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 =\n", | |
| "\n", | |
| " 0.0024\n", | |
| " 0.0024\n", | |
| " 0.0024\n", | |
| " 0.0024\n", | |
| " 0.1024\n", | |
| " 0.0024\n", | |
| " 0.0024\n", | |
| " 0.0024\n", | |
| " 0.0024\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "A = magic(9); b = (1:9)';\n", | |
| "[Q,R] = qr(A);\n", | |
| "z = Q'*b;\n", | |
| "x(9,1) = z(9)/R(9,9);\n", | |
| "for i = 8:-1:1\n", | |
| " x(i) = (z(i) - R(i,i+1:9)*x(i+1:9)) / R(i,i);\n", | |
| "end\n", | |
| "x" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 13, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "ans =\n", | |
| "\n", | |
| " 3.1792e-15\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "norm(A*x-b)" | |
| ] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Matlab", | |
| "language": "matlab", | |
| "name": "matlab" | |
| }, | |
| "language_info": { | |
| "codemirror_mode": "octave", | |
| "file_extension": ".m", | |
| "help_links": [ | |
| { | |
| "text": "MetaKernel Magics", | |
| "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" | |
| } | |
| ], | |
| "mimetype": "text/x-octave", | |
| "name": "matlab", | |
| "version": "0.11.0" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 1 | |
| } |
This file contains hidden or 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": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "A =\n", | |
| "\n", | |
| " 35 1 6 26\n", | |
| " 3 32 7 21\n", | |
| " 31 9 2 22\n", | |
| " 8 28 33 17\n", | |
| " 30 5 34 12\n", | |
| " 4 36 29 13\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "A = magic(6); A = A(:,1: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": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "R =\n", | |
| "\n", | |
| " 56.3471\n", | |
| "\n", | |
| "\n", | |
| "Q =\n", | |
| "\n", | |
| " 0.6211\n", | |
| " 0.0532\n", | |
| " 0.5502\n", | |
| " 0.1420\n", | |
| " 0.5324\n", | |
| " 0.0710\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "R = norm(A(:,1)), Q = A(:,1)/R" | |
| ] | |
| }, | |
| { | |
| "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": [ | |
| "innerprod =\n", | |
| "\n", | |
| " -2.2204e-15\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "P = Q*Q'; \n", | |
| "v = A(:,2) - P*A(:,2);\n", | |
| "innerprod = v'*Q" | |
| ] | |
| }, | |
| { | |
| "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 =\n", | |
| "\n", | |
| " 56.3471 16.4693\n", | |
| " 0 54.2196\n", | |
| "\n", | |
| "\n", | |
| "Q =\n", | |
| "\n", | |
| " 0.6211 -0.1702\n", | |
| " 0.0532 0.5740\n", | |
| " 0.5502 -0.0011\n", | |
| " 0.1420 0.4733\n", | |
| " 0.5324 -0.0695\n", | |
| " 0.0710 0.6424\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "R(1,2) = Q(:,1)'*A(:,2);\n", | |
| "v = A(:,2) - Q(:,1)*R(1,2);\n", | |
| "R(2,2) = norm(v)\n", | |
| "Q(:,2) = v/R(2,2)" | |
| ] | |
| }, | |
| { | |
| "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": [ | |
| "innerprod =\n", | |
| "\n", | |
| " 1.0e-14 *\n", | |
| "\n", | |
| " 0.2228 -0.2821\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "P = Q(:,1:2)*Q(:,1:2)';\n", | |
| "v = A(:,2) - P*A(:,2);\n", | |
| "innerprod = v'*Q" | |
| ] | |
| }, | |
| { | |
| "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 =\n", | |
| "\n", | |
| " 56.3471 16.4693 30.0459\n", | |
| " 0 54.2196 34.8797\n", | |
| " 0 0 32.4907\n", | |
| "\n", | |
| "\n", | |
| "Q =\n", | |
| "\n", | |
| " 0.6211 -0.1702 -0.2070\n", | |
| " 0.0532 0.5740 -0.4500\n", | |
| " 0.5502 -0.0011 -0.4460\n", | |
| " 0.1420 0.4733 0.3763\n", | |
| " 0.5324 -0.0695 0.6287\n", | |
| " 0.0710 0.6424 0.1373\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "R(1:2,3) = Q(:,1:2)'*A(:,3);\n", | |
| "v = A(:,3) - Q(:,1:2)*R(1:2,3);\n", | |
| "R(3,3) = norm(v)\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": [ | |
| "ans =\n", | |
| "\n", | |
| " 3.5527e-15\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "norm( Q*R - 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": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "R = norm(A(:,1));\n", | |
| "Q = A(:,1)/R(1,1);\n", | |
| "for j = 2:size(A,2)\n", | |
| " R(1:j-1,j) = Q'*A(:,j);\n", | |
| " v = A(:,j) - Q*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 =\n", | |
| "\n", | |
| " 0.6211 -0.1702 -0.2070 0.4998\n", | |
| " 0.0532 0.5740 -0.4500 0.2106\n", | |
| " 0.5502 -0.0011 -0.4460 -0.4537\n", | |
| " 0.1420 0.4733 0.3763 0.5034\n", | |
| " 0.5324 -0.0695 0.6287 -0.2096\n", | |
| " 0.0710 0.6424 0.1373 -0.4501\n", | |
| "\n", | |
| "\n", | |
| "R =\n", | |
| "\n", | |
| " 56.3471 16.4693 30.0459 39.0969\n", | |
| " 0 54.2196 34.8797 23.1669\n", | |
| " 0 0 32.4907 -8.9182\n", | |
| " 0 0 0 7.6283\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "Q,R" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 10, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "ans =\n", | |
| "\n", | |
| " 1.2510e-15\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "norm(Q'*Q-eye(4))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 11, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "ans =\n", | |
| "\n", | |
| " 7.9441e-15\n" | |
| ] | |
| } | |
| ], | |
| "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 =\n", | |
| "\n", | |
| " 0.0024\n", | |
| " 0.0024\n", | |
| " 0.0024\n", | |
| " 0.0024\n", | |
| " 0.1024\n", | |
| " 0.0024\n", | |
| " 0.0024\n", | |
| " 0.0024\n", | |
| " 0.0024\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "A = magic(9); b = (1:9)';\n", | |
| "[Q,R] = qr(A);\n", | |
| "z = Q'*b;\n", | |
| "x(9,1) = z(9)/R(9,9);\n", | |
| "for i = 8:-1:1\n", | |
| " x(i) = (z(i) - R(i,i+1:9)*x(i+1:9)) / R(i,i);\n", | |
| "end\n", | |
| "x" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 13, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "ans =\n", | |
| "\n", | |
| " 3.1792e-15\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "norm(A*x-b)" | |
| ] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Matlab", | |
| "language": "matlab", | |
| "name": "matlab" | |
| }, | |
| "language_info": { | |
| "codemirror_mode": "octave", | |
| "file_extension": ".m", | |
| "help_links": [ | |
| { | |
| "text": "MetaKernel Magics", | |
| "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" | |
| } | |
| ], | |
| "mimetype": "text/x-octave", | |
| "name": "matlab", | |
| "version": "0.11.0" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 1 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment