{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Stanford CS205a \n", "\n", "## Column Spaces and QR\n", "Supplemental Julia Notebook\n", "\n", "Prof. Doug James" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Condition number of \n", "\\begin{align}\n", " \\textrm{cond} \\; A^T A & = \\|A^T A\\| \\|(A^T A)^{-1}\\| \\\\\n", " & \\approx \\|A^T\\| \\|A\\| \\|A^{-1}\\| \\|(A^T)^{-1}\\| \\qquad \\mbox{(esp. in two-norm)}\\\\\n", " & = \\|A\\|^2 \\|A^{-1}\\|^2 \\\\\n", " & = (\\textrm{cond} \\; A)^2 \n", "\\end{align}\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Experiment: $\\textrm{cond} \\; A^T A \\approx (\\textrm{cond} \\; A)^2$" ] }, { "cell_type": "code", "execution_count": 285, "metadata": { "collapsed": true }, "outputs": [], "source": [ "N = 10;\n", "A = rand(N,N);\n", "AtA = A'*A;" ] }, { "cell_type": "code", "execution_count": 286, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3222.353911452386" ] }, "execution_count": 286, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cond(A)" ] }, { "cell_type": "code", "execution_count": 287, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1.038356473065249e7" ] }, "execution_count": 287, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(cond(A))^2" ] }, { "cell_type": "code", "execution_count": 288, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1.0383564733178914e7" ] }, "execution_count": 288, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cond(AtA)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Gram-Schmidt Algorithm (Example 5.1)" ] }, { "cell_type": "code", "execution_count": 289, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 1.0 1.0\n", " 0.0 1.0 1.0\n", " 0.0 1.0 0.0" ] }, "execution_count": 289, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Input independent vectors\n", "v1 = [1.; 0.; 0.];\n", "v2 = [1.; 1.; 1.];\n", "v3 = [1.; 1.; 0.];\n", "A = [v1 v2 v3]" ] }, { "cell_type": "code", "execution_count": 290, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 1.0\n", " 0.0\n", " 0.0" ] }, "execution_count": 290, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# VECTOR v1.\n", "# Normalize v1:\n", "a1 = v1/norm(v1) ## book calls these a_i but q_i is better" ] }, { "cell_type": "code", "execution_count": 291, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 0.0\n", " 1.0\n", " 1.0" ] }, "execution_count": 291, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# VECTOR v2.\n", "# Remove the span of a1 from v2:\n", "a2 = v2 - a1*a1'*v2" ] }, { "cell_type": "code", "execution_count": 292, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 0.0 \n", " 0.707107\n", " 0.707107" ] }, "execution_count": 292, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# normalize a2:\n", "a2 /= norm(a2)" ] }, { "cell_type": "code", "execution_count": 293, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 0.0\n", " 0.5\n", " -0.5" ] }, "execution_count": 293, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# VECTOR v3.\n", "# Remove the span of a1 & a2 from v3: \n", "a3 = v3 - a1*a1'*v3 - a2*a2'*v3" ] }, { "cell_type": "code", "execution_count": 294, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 0.0 \n", " 0.707107\n", " -0.707107" ] }, "execution_count": 294, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# normalize a3:\n", "a3 /= norm(a3)" ] }, { "cell_type": "code", "execution_count": 295, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 0.0 0.0 \n", " 0.0 0.707107 0.707107\n", " 0.0 0.707107 -0.707107" ] }, "execution_count": 295, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Orthogonal basis matrix formed by Gram-Schmidt:\n", "Q = [a1 a2 a3]" ] }, { "cell_type": "code", "execution_count": 296, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 0.0 0.0 \n", " 0.0 1.0 2.77556e-16\n", " 0.0 2.77556e-16 1.0 " ] }, "execution_count": 296, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Q'*Q" ] }, { "cell_type": "code", "execution_count": 297, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 1.0 1.0 \n", " 0.0 1.41421 0.707107\n", " 0.0 3.33067e-16 0.707107" ] }, "execution_count": 297, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Upper triangle matrix (encoded in our operations above)\n", "R = Q'*A" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## QR Implementations in Julia" ] }, { "cell_type": "code", "execution_count": 298, "metadata": { "collapsed": false, "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "search: \u001b[1mq\u001b[22m\u001b[1mr\u001b[22m \u001b[1mq\u001b[22m\u001b[1mr\u001b[22mfact \u001b[1mq\u001b[22m\u001b[1mr\u001b[22mfact! s\u001b[1mq\u001b[22m\u001b[1mr\u001b[22mt s\u001b[1mq\u001b[22m\u001b[1mr\u001b[22mtm is\u001b[1mq\u001b[22m\u001b[1mr\u001b[22mt \u001b[1mQ\u001b[22muickSo\u001b[1mr\u001b[22mt Partial\u001b[1mQ\u001b[22muickSo\u001b[1mr\u001b[22mt\n", "\n" ] }, { "data": { "text/markdown": [ "```\n", "qr(A [,pivot=Val{false}][;thin=true]) -> Q, R, [p]\n", "```\n", "\n", "Compute the (pivoted) QR factorization of `A` such that either `A = Q*R` or `A[:,p] = Q*R`. Also see `qrfact`. The default is to compute a thin factorization. Note that `R` is not extended with zeros when the full `Q` is requested.\n", "\n", "```\n", "qr(v::AbstractVector)\n", "```\n", "\n", "Computes the polar decomposition of a vector.\n", "\n", "Input:\n", "\n", " * `v::AbstractVector` - vector to normalize\n", "\n", "Outputs:\n", "\n", " * `w` - A unit vector in the direction of `v`\n", " * `r` - The norm of `v`\n", "\n", "See also:\n", "\n", "`normalize`, `normalize!`, `LinAlg.qr!`\n" ], "text/plain": [ "```\n", "qr(A [,pivot=Val{false}][;thin=true]) -> Q, R, [p]\n", "```\n", "\n", "Compute the (pivoted) QR factorization of `A` such that either `A = Q*R` or `A[:,p] = Q*R`. Also see `qrfact`. The default is to compute a thin factorization. Note that `R` is not extended with zeros when the full `Q` is requested.\n", "\n", "```\n", "qr(v::AbstractVector)\n", "```\n", "\n", "Computes the polar decomposition of a vector.\n", "\n", "Input:\n", "\n", " * `v::AbstractVector` - vector to normalize\n", "\n", "Outputs:\n", "\n", " * `w` - A unit vector in the direction of `v`\n", " * `r` - The norm of `v`\n", "\n", "See also:\n", "\n", "`normalize`, `normalize!`, `LinAlg.qr!`\n" ] }, "execution_count": 298, "metadata": {}, "output_type": "execute_result" } ], "source": [ "?(qr)" ] }, { "cell_type": "code", "execution_count": 299, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "search: \u001b[1mq\u001b[22m\u001b[1mr\u001b[22m\u001b[1mf\u001b[22m\u001b[1ma\u001b[22m\u001b[1mc\u001b[22m\u001b[1mt\u001b[22m \u001b[1mq\u001b[22m\u001b[1mr\u001b[22m\u001b[1mf\u001b[22m\u001b[1ma\u001b[22m\u001b[1mc\u001b[22m\u001b[1mt\u001b[22m!\n", "\n" ] }, { "data": { "text/markdown": [ "```\n", "qrfact(A) -> SPQR.Factorization\n", "```\n", "\n", "Compute the QR factorization of a sparse matrix `A`. A fill-reducing permutation is used. The main application of this type is to solve least squares problems with `\\`. The function calls the C library SPQR and a few additional functions from the library are wrapped but not exported.\n", "\n", "```\n", "qrfact(A [,pivot=Val{false}]) -> F\n", "```\n", "\n", "Computes the QR factorization of `A`. The return type of `F` depends on the element type of `A` and whether pivoting is specified (with `pivot==Val{true}`).\n", "\n", "| Return type | `eltype(A)` | `pivot` | Relationship between `F` and `A` |\n", "|:------------- |:--------------- |:------------ |:-------------------------------- |\n", "| `QR` | not `BlasFloat` | either | `A==F[:Q]*F[:R]` |\n", "| `QRCompactWY` | `BlasFloat` | `Val{false}` | `A==F[:Q]*F[:R]` |\n", "| `QRPivoted` | `BlasFloat` | `Val{true}` | `A[:,F[:p]]==F[:Q]*F[:R]` |\n", "\n", "`BlasFloat` refers to any of: `Float32`, `Float64`, `Complex64` or `Complex128`.\n", "\n", "The individual components of the factorization `F` can be accessed by indexing:\n", "\n", "| Component | Description | `QR` | `QRCompactWY` | `QRPivoted` |\n", "|:--------- |:----------------------------------------- |:--------------- |:------------------ |:--------------- |\n", "| `F[:Q]` | `Q` (orthogonal/unitary) part of `QR` | ✓ (`QRPackedQ`) | ✓ (`QRCompactWYQ`) | ✓ (`QRPackedQ`) |\n", "| `F[:R]` | `R` (upper right triangular) part of `QR` | ✓ | ✓ | ✓ |\n", "| `F[:p]` | pivot `Vector` | | | ✓ |\n", "| `F[:P]` | (pivot) permutation `Matrix` | | | ✓ |\n", "\n", "The following functions are available for the `QR` objects: `size`, `\\`. When `A` is rectangular, `\\` will return a least squares solution and if the solution is not unique, the one with smallest norm is returned.\n", "\n", "Multiplication with respect to either thin or full `Q` is allowed, i.e. both `F[:Q]*F[:R]` and `F[:Q]*A` are supported. A `Q` matrix can be converted into a regular matrix with [`full`](:func:`full`) which has a named argument `thin`.\n", "\n", "!!! note\n", " `qrfact` returns multiple types because LAPACK uses several representations that minimize the memory storage requirements of products of Householder elementary reflectors, so that the `Q` and `R` matrices can be stored compactly rather as two separate dense matrices.\n", "\n", " The data contained in `QR` or `QRPivoted` can be used to construct the `QRPackedQ` type, which is a compact representation of the rotation matrix:\n", "\n", " $$\n", " Q = \\prod_{i=1}^{\\min(m,n)} (I - \\tau_i v_i v_i^T)\n", " $$\n", "\n", " where $\\tau_i$ is the scale factor and $v_i$ is the projection vector associated with the $i^{th}$ Householder elementary reflector.\n", "\n", " The data contained in `QRCompactWY` can be used to construct the `QRCompactWYQ` type, which is a compact representation of the rotation matrix\n", "\n", " $$\n", " Q = I + Y T Y^T\n", " $$\n", "\n", " where `Y` is $m \\times r$ lower trapezoidal and `T` is $r \\times r$ upper triangular. The *compact WY* representation [^Schreiber1989] is not to be confused with the older, *WY* representation [^Bischof1987]. (The LAPACK documentation uses `V` in lieu of `Y`.)\n", "\n", " [^Bischof1987]: C Bischof and C Van Loan, \"The WY representation for products of Householder matrices\", SIAM J Sci Stat Comput 8 (1987), s2-s13. [doi:10.1137/0908009](http://dx.doi.org/10.1137/0908009)\n", "\n", " [^Schreiber1989]: R Schreiber and C Van Loan, \"A storage-efficient WY representation for products of Householder transformations\", SIAM J Sci Stat Comput 10 (1989), 53-57. [doi:10.1137/0910005](http://dx.doi.org/10.1137/0910005)\n", "\n" ], "text/plain": [ "```\n", "qrfact(A) -> SPQR.Factorization\n", "```\n", "\n", "Compute the QR factorization of a sparse matrix `A`. A fill-reducing permutation is used. The main application of this type is to solve least squares problems with `\\`. The function calls the C library SPQR and a few additional functions from the library are wrapped but not exported.\n", "\n", "```\n", "qrfact(A [,pivot=Val{false}]) -> F\n", "```\n", "\n", "Computes the QR factorization of `A`. The return type of `F` depends on the element type of `A` and whether pivoting is specified (with `pivot==Val{true}`).\n", "\n", "| Return type | `eltype(A)` | `pivot` | Relationship between `F` and `A` |\n", "|:------------- |:--------------- |:------------ |:-------------------------------- |\n", "| `QR` | not `BlasFloat` | either | `A==F[:Q]*F[:R]` |\n", "| `QRCompactWY` | `BlasFloat` | `Val{false}` | `A==F[:Q]*F[:R]` |\n", "| `QRPivoted` | `BlasFloat` | `Val{true}` | `A[:,F[:p]]==F[:Q]*F[:R]` |\n", "\n", "`BlasFloat` refers to any of: `Float32`, `Float64`, `Complex64` or `Complex128`.\n", "\n", "The individual components of the factorization `F` can be accessed by indexing:\n", "\n", "| Component | Description | `QR` | `QRCompactWY` | `QRPivoted` |\n", "|:--------- |:----------------------------------------- |:--------------- |:------------------ |:--------------- |\n", "| `F[:Q]` | `Q` (orthogonal/unitary) part of `QR` | ✓ (`QRPackedQ`) | ✓ (`QRCompactWYQ`) | ✓ (`QRPackedQ`) |\n", "| `F[:R]` | `R` (upper right triangular) part of `QR` | ✓ | ✓ | ✓ |\n", "| `F[:p]` | pivot `Vector` | | | ✓ |\n", "| `F[:P]` | (pivot) permutation `Matrix` | | | ✓ |\n", "\n", "The following functions are available for the `QR` objects: `size`, `\\`. When `A` is rectangular, `\\` will return a least squares solution and if the solution is not unique, the one with smallest norm is returned.\n", "\n", "Multiplication with respect to either thin or full `Q` is allowed, i.e. both `F[:Q]*F[:R]` and `F[:Q]*A` are supported. A `Q` matrix can be converted into a regular matrix with [`full`](:func:`full`) which has a named argument `thin`.\n", "\n", "!!! note\n", " `qrfact` returns multiple types because LAPACK uses several representations that minimize the memory storage requirements of products of Householder elementary reflectors, so that the `Q` and `R` matrices can be stored compactly rather as two separate dense matrices.\n", "\n", " The data contained in `QR` or `QRPivoted` can be used to construct the `QRPackedQ` type, which is a compact representation of the rotation matrix:\n", "\n", " $$\n", " Q = \\prod_{i=1}^{\\min(m,n)} (I - \\tau_i v_i v_i^T)\n", " $$\n", "\n", " where $\\tau_i$ is the scale factor and $v_i$ is the projection vector associated with the $i^{th}$ Householder elementary reflector.\n", "\n", " The data contained in `QRCompactWY` can be used to construct the `QRCompactWYQ` type, which is a compact representation of the rotation matrix\n", "\n", " $$\n", " Q = I + Y T Y^T\n", " $$\n", "\n", " where `Y` is $m \\times r$ lower trapezoidal and `T` is $r \\times r$ upper triangular. The *compact WY* representation [^Schreiber1989] is not to be confused with the older, *WY* representation [^Bischof1987]. (The LAPACK documentation uses `V` in lieu of `Y`.)\n", "\n", " [^Bischof1987]: C Bischof and C Van Loan, \"The WY representation for products of Householder matrices\", SIAM J Sci Stat Comput 8 (1987), s2-s13. [doi:10.1137/0908009](http://dx.doi.org/10.1137/0908009)\n", "\n", " [^Schreiber1989]: R Schreiber and C Van Loan, \"A storage-efficient WY representation for products of Householder transformations\", SIAM J Sci Stat Comput 10 (1989), 53-57. [doi:10.1137/0910005](http://dx.doi.org/10.1137/0910005)\n", "\n" ] }, "execution_count": 299, "metadata": {}, "output_type": "execute_result" } ], "source": [ "?(qrfact)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## QR Example" ] }, { "cell_type": "code", "execution_count": 300, "metadata": { "collapsed": false, "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "5×2 Array{Float64,2}:\n", " 0.260917 0.555757\n", " 0.68724 0.477956\n", " 0.862184 0.323348\n", " 0.230656 0.841944\n", " 0.405023 0.239465" ] }, "execution_count": 300, "metadata": {}, "output_type": "execute_result" } ], "source": [ "M = 5;\n", "N = 2;\n", "A = rand(M,N)" ] }, { "cell_type": "code", "execution_count": 301, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Base.LinAlg.QRCompactWY{Float64,Array{Float64,2}}([-1.22515 -0.851695; 0.462457 0.826251; … ; 0.155213 -0.624; 0.272548 0.144251],[1.21297 -1.01404; 6.94162e-310 1.20929])" ] }, "execution_count": 301, "metadata": {}, "output_type": "execute_result" } ], "source": [ "QR1 = qrfact(A, Val{false}) ## QR with pivot=false (for demo purposes)" ] }, { "cell_type": "code", "execution_count": 302, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "5×5 Base.LinAlg.QRCompactWYQ{Float64,Array{Float64,2}}:\n", " -0.212968 0.453099 -0.203175 -0.821032 -0.184315 \n", " -0.560945 0.00024438 -0.690906 0.374908 -0.259679 \n", " -0.70374 -0.334067 0.58745 -0.103852 -0.193046 \n", " -0.188268 0.824927 0.340959 0.401694 0.0802437\n", " -0.330592 -0.0509507 -0.141485 -0.114919 0.924602 " ] }, "execution_count": 302, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Q1 = QR1[:Q] ## \"Fat\" QR with null-space basis, too" ] }, { "cell_type": "code", "execution_count": 303, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2×2 Array{Float64,2}:\n", " -1.22515 -0.851695\n", " 0.0 0.826251" ] }, "execution_count": 303, "metadata": {}, "output_type": "execute_result" } ], "source": [ "R1 = QR1[:R]" ] }, { "cell_type": "code", "execution_count": 304, "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "text/plain": [ "5×2 Array{Float64,2}:\n", " 0.260917 0.555757\n", " 0.68724 0.477956\n", " 0.862184 0.323348\n", " 0.230656 0.841944\n", " 0.405023 0.239465" ] }, "execution_count": 304, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Note: Julia permits dimension mismatch here with QRCompactWYQ\n", "Q1 * R1 " ] }, { "cell_type": "code", "execution_count": 305, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "5×2 Array{Float64,2}:\n", " 0.260917 0.555757\n", " 0.68724 0.477956\n", " 0.862184 0.323348\n", " 0.230656 0.841944\n", " 0.405023 0.239465" ] }, "execution_count": 305, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 306, "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "text/plain": [ "(\n", "[-0.212968 0.453099; -0.560945 0.00024438; … ; -0.188268 0.824927; -0.330592 -0.0509507],\n", "\n", "[-1.22515 -0.851695; 0.0 0.826251])" ] }, "execution_count": 306, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(Q2, R2) = qr(A) ## Default: unpivoted \"thin QR\"" ] }, { "cell_type": "code", "execution_count": 307, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "5×2 Array{Float64,2}:\n", " -0.212968 0.453099 \n", " -0.560945 0.00024438\n", " -0.70374 -0.334067 \n", " -0.188268 0.824927 \n", " -0.330592 -0.0509507 " ] }, "execution_count": 307, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Q2 ## thin basis" ] }, { "cell_type": "code", "execution_count": 308, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2×2 Array{Float64,2}:\n", " -1.22515 -0.851695\n", " 0.0 0.826251" ] }, "execution_count": 308, "metadata": {}, "output_type": "execute_result" } ], "source": [ "R2" ] }, { "cell_type": "code", "execution_count": 309, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "5×2 Array{Float64,2}:\n", " 0.260917 0.555757\n", " 0.68724 0.477956\n", " 0.862184 0.323348\n", " 0.230656 0.841944\n", " 0.405023 0.239465" ] }, "execution_count": 309, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Q2*R2" ] }, { "cell_type": "code", "execution_count": 310, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "5×2 Array{Float64,2}:\n", " 0.260917 0.555757\n", " 0.68724 0.477956\n", " 0.862184 0.323348\n", " 0.230656 0.841944\n", " 0.405023 0.239465" ] }, "execution_count": 310, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Solve Ax=b with $\\vec{b}=\\hat{1}$" ] }, { "cell_type": "code", "execution_count": 311, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2-element Array{Float64,1}:\n", " 0.392681\n", " 0.483478" ] }, "execution_count": 311, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b = ones(M)/sqrt(M); ## unit vector\n", "x1 = QR1 \\ b ## QR1 from qrfact (fat Q)" ] }, { "cell_type": "code", "execution_count": 312, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2-element Array{Float64,1}:\n", " 0.392681\n", " 0.483478" ] }, "execution_count": 312, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x2 = R2 \\ (Q2'*b) ## QR2 from qr (thin Q, no pivoting)" ] }, { "cell_type": "code", "execution_count": 313, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(0.20786446598748637,0.2078644659874864)" ] }, "execution_count": 313, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Residual:\n", "norm(A*x1-b), norm(A*x2-b)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Julia 0.5.1", "language": "julia", "name": "julia-0.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.5.1" } }, "nbformat": 4, "nbformat_minor": 1 }