{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"toc": "true"
},
"source": [
"# Table of Contents\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Julia Version 1.1.0\n",
"Commit 80516ca202 (2019-01-21 21:24 UTC)\n",
"Platform Info:\n",
" OS: macOS (x86_64-apple-darwin14.5.0)\n",
" CPU: Intel(R) Core(TM) i7-6920HQ CPU @ 2.90GHz\n",
" WORD_SIZE: 64\n",
" LIBM: libopenlibm\n",
" LLVM: libLLVM-6.0.1 (ORCJIT, skylake)\n",
"Environment:\n",
" JULIA_EDITOR = code\n"
]
}
],
"source": [
"versioninfo()"
]
},
{
"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/v1/stdlib/LinearAlgebra/#LinearAlgebra.cond) computes $\\kappa_p$ for $p=1$, 2 (default), or $\\infty$.\n",
"\n",
"* Julia function [`condskeel`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.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": 2,
"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": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using DelimitedFiles\n",
"\n",
"# Base.download(\"http://hua-zhou.github.io/teaching/biostatm280-2019spring/slides/14-cond/longley.txt\", \"./longley.txt\")\n",
"longley = readdlm(\"longley.txt\")"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"