{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"toc": "true"
},
"source": [
"# Table of Contents\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Condition Number\n",
"\n",
"* Assume $\\mathbf{A} \\in \\mathbb{R}^{n \\times n}$ is nonsingular and consider the system of linear equation \n",
"$$\n",
"\\mathbf{A} \\mathbf{x} = \\mathbf{b}.\n",
"$$\n",
"The solution is \n",
"$$\n",
"\\mathbf{x} = \\mathbf{A}^{-1} \\mathbf{b}.\n",
"$$\n",
"We want to know how the solution changes with a small perturbation of the input $\\mathbf{b}$ (or $\\mathbf{A}$).\n",
"\n",
"* Let \n",
"$$\n",
"\\mathbf{b}_{\\text{new}} = \\mathbf{b} + \\Delta \\mathbf{b}.\n",
"$$\n",
"Then \n",
"$$\n",
"\\mathbf{x}_{\\text{new}} = \\mathbf{A}^{-1} (\\mathbf{b} + \\Delta \\mathbf{b}) = \\mathbf{x} + \\Delta \\mathbf{x}.\n",
"$$\n",
"Thus \n",
"$$\n",
"\\|\\Delta \\mathbf{x}\\| = \\|\\mathbf{A}^{-1} \\Delta \\mathbf{b}\\| \\le \\|\\mathbf{A}^{-1}\\| \\|\\Delta \\mathbf{b}\\|.\n",
"$$\n",
"Because $\\mathbf{b} = \\mathbf{A} \\mathbf{x}$, \n",
"$$\n",
"\\frac{1}{\\|\\mathbf{x}\\|} \\le \\|\\mathbf{A}\\| \\frac{1}{\\|\\mathbf{b}\\|}.\n",
"$$\n",
"This results \n",
"$$\n",
"\\frac{ \\|\\Delta \\mathbf{x}\\|}{\\|\\mathbf{x}\\|} \\le \\|\\mathbf{A}\\|\\|\\mathbf{A}^{-1}\\| \\frac{\\|\\Delta \\mathbf{b}\\|}{\\|\\mathbf{b}\\|}.\n",
"$$\n",
"\n",
"* $\\kappa(\\mathbf{A}) = \\|\\mathbf{A}\\| \\|\\mathbf{A}^{-1}\\|$ is called the **condition number for linear equation**. It depends on the matrix norm being used. \n",
" * $\\kappa_p$ means condition number defined from matrix-$p$ norm.\n",
"\n",
"* Large condition number means \"bad\".\n",
"\n",
"* Some facts: \n",
"$$\n",
"\\begin{eqnarray*}\n",
"\\kappa(\\mathbf{A}) &=& \\kappa(\\mathbf{A}^{-1})\t\\\\\n",
"\\kappa(c\\mathbf{A}) &=& \\kappa(\\mathbf{A})\t\\\\\n",
"\\kappa(\\mathbf{A}) &\\ge& 1\t\\\\\n",
"%\\kappa_1(\\mathbf{A}) &=& \\kappa_\\infty (\\mathbf{A}^T)\t\\\\\n",
"\\kappa_2 (\\mathbf{A}) &=& \\kappa_2 (\\mathbf{A}^T) = \\frac{\\sigma_1(\\mathbf{A})}{\\sigma_n(\\mathbf{A})}\t\\\\\n",
"\\kappa_2(\\mathbf{A}^T \\mathbf{A}) &=& \\frac{\\lambda_1(\\mathbf{A}^T \\mathbf{A})}{\\lambda_n(\\mathbf{A}^T \\mathbf{A})} = \\kappa_2^2(\\mathbf{A}) \\ge \\kappa_2(\\mathbf{A}).\n",
"\\end{eqnarray*}\n",
"$$\n",
"The last fact says that the condition number of $\\mathbf{A}^T \\mathbf{A}$ can be much larger than that of $\\mathbf{A}$.\n",
"\n",
"* The smallest singular value $\\sigma_n$ indicates the _distance to the trouble_.\n",
"\n",
"* **Condition number for the least squares problem** is more complicated. Roughly speaking, \n",
" - the method based on normal equation (Cholesky, sweep) has a condition depending on $[\\kappa_2(\\mathbf{X})]^2$ \n",
" - QR for a _small residuals_ least squares problem has a condition depending on $\\kappa_2(\\mathbf{X})$\n",
" \n",
"* Consider the simple case\n",
"$$\n",
"\\begin{eqnarray*}\n",
"\t\\mathbf{X} = \\begin{pmatrix}\n",
"\t1 & 1 \\\\\n",
"\t10^{-3} & 0 \\\\\n",
"\t0 & 10^{-3}\n",
"\t\\end{pmatrix}.\n",
"\\end{eqnarray*}\n",
"$$\n",
"Forming normal equation yields a singular Gramian matrix\n",
"$$\n",
"\\begin{eqnarray*}\n",
"\t\\mathbf{X}^T \\mathbf{X} = \\begin{pmatrix}\n",
"\t1 & 1 \\\\\n",
"\t1 & 1\n",
"\t\\end{pmatrix}\n",
"\\end{eqnarray*}\n",
"$$\n",
"if executed with a precision of 6 decimal digits."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Implementation\n",
"\n",
"* Julia function [`cond`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.cond) computes $\\kappa_p$ for $p=1$, 2 (default), or $\\infty$.\n",
"\n",
"* Julia function [`condskeel`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.condskeel) computes the Skeel condition number."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Longley example\n",
"\n",
"The [Longley (1967)](https://www.jstor.org/stable/2283673?seq=1#page_scan_tab_contents) macroeconomic data set is a famous test example for numerical software in early dates."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"16×7 Array{Float64,2}:\n",
" 60323.0 83.0 234289.0 2356.0 1590.0 107608.0 1947.0\n",
" 61122.0 88.5 259426.0 2325.0 1456.0 108632.0 1948.0\n",
" 60171.0 88.2 258054.0 3682.0 1616.0 109773.0 1949.0\n",
" 61187.0 89.5 284599.0 3351.0 1650.0 110929.0 1950.0\n",
" 63221.0 96.2 328975.0 2099.0 3099.0 112075.0 1951.0\n",
" 63639.0 98.1 346999.0 1932.0 3594.0 113270.0 1952.0\n",
" 64989.0 99.0 365385.0 1870.0 3547.0 115094.0 1953.0\n",
" 63761.0 100.0 363112.0 3578.0 3350.0 116219.0 1954.0\n",
" 66019.0 101.2 397469.0 2904.0 3048.0 117388.0 1955.0\n",
" 67857.0 104.6 419180.0 2822.0 2857.0 118734.0 1956.0\n",
" 68169.0 108.4 442769.0 2936.0 2798.0 120445.0 1957.0\n",
" 66513.0 110.8 444546.0 4681.0 2637.0 121950.0 1958.0\n",
" 68655.0 112.6 482704.0 3813.0 2552.0 123366.0 1959.0\n",
" 69564.0 114.2 502601.0 3931.0 2514.0 125368.0 1960.0\n",
" 69331.0 115.7 518173.0 4806.0 2572.0 127852.0 1961.0\n",
" 70551.0 116.9 554894.0 4007.0 2827.0 130081.0 1962.0"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Requests\n",
"\n",
"longleyurl = get(\"http://hua-zhou.github.io/teaching/biostatm280-2018spring/slides/13-cond/longley.txt\")\n",
"longley = readdlm(IOBuffer(longleyurl.data))"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mKeyword argument match_dimensions not supported with Plots.GRBackend(). Choose from: Set(Symbol[:top_margin, :group, :background_color, :yforeground_color_text, :yguidefontcolor, :seriesalpha, :legendfontcolor, :seriescolor, :ztick_direction, :zlims, :overwrite_figure, :xguidefonthalign, :normalize, :linestyle, :xflip, :fillcolor, :ygrid, :background_color_inside, :zguidefonthalign, :bins, :yscale, :xtickfontcolor, :xguide, :fillalpha, :tick_direction, :yguidefontsize, :legendfontfamily, :foreground_color, :xtickfonthalign, :x, :ytickfontrotation, :legend, :discrete_values, :ytick_direction, :xguidefontrotation, :ribbon, :tickfontrotation, :xdiscrete_values, :legendtitle, :xgridstyle, :orientation, :gridstyle, :markersize, :camera, :xforeground_color_grid, :quiver, :zticks, :markerstrokecolor, :ztickfontrotation, :ztickfonthalign, :legendfonthalign, :xtickfontsize, :levels, :zgridstyle, :foreground_color_border, :zguidefontvalign, :marker_z, :markerstrokealpha, :markeralpha, :tickfontvalign, :zguidefontcolor, :ygridlinewidth, :zlink, :zscale, :smooth, :xticks, :zguidefontsize, :y, :margin, :ytickfontcolor, :yforeground_color_border, :zguidefontfamily, :zgridalpha, :yguidefontvalign, :yguidefonthalign, :ztickfontcolor, :html_output_format, :tickfontcolor, :titlefontrotation, :legendfontvalign, :tickfontsize, :z, :yforeground_color_axis, :xtickfontrotation, :xerror, :contour_labels, :xguidefontcolor, :primary, :guidefonthalign, :aspect_ratio, :link, :yguide, :guidefontvalign, :yguidefontfamily, :layout, :polar, :right_margin, :xlink, :series_annotations, :inset_subplots, :ytickfontsize, :tickfontfamily, :xgrid, :ygridalpha, :xtick_direction, :colorbar, :zflip, :ticks, :legendfontrotation, :linealpha, :arrow, :xtickfontvalign, :zgrid, :bar_width, :zguide, :zforeground_color_text, :weights, :xgridalpha, :ygridstyle, :fill_z, :ztickfontfamily, :markershape, :background_color_subplot, :xguidefontvalign, :markerstrokewidth, :xguidefontfamily, :gridlinewidth, :foreground_color_subplot, :xgridlinewidth, :foreground_color_text, :titlefonthalign, :yerror, :zgridlinewidth, :grid, :xguidefontsize, :xforeground_color_axis, :background_color_outside, :titlefontcolor, :line_z, :size, :projection, :zguidefontrotation, :ydiscrete_values, :seriestype, :yflip, :fillrange, :ztickfontvalign, :xlims, :xforeground_color_border, :markercolor, :ylink, :yforeground_color_grid, :color_palette, :lims, :xscale, :left_margin, :annotations, :window_title, :foreground_color_axis, :yguidefontrotation, :guidefontsize, :zdiscrete_values, :tickfonthalign, :bottom_margin, :framestyle, :scale, :zforeground_color_border, :background_color_legend, :linecolor, :foreground_color_legend, :title, :subplot_index, :flip, :titlefontvalign, :foreground_color_grid, :linewidth, :ztickfontsize, :gridalpha, :guidefontfamily, :ylims, :xtickfontfamily, :ytickfontvalign, :ytickfontfamily, :xforeground_color_text, :show, :guidefontrotation, :legendfontsize, :subplot, :label, :ytickfonthalign, :guide, :guidefontcolor, :titlefontsize, :titlefontfamily, :zforeground_color_axis, :zforeground_color_grid, :yticks])\u001b[39m\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"