{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Stanford CS205a \n",
"\n",
"## Norms, Sensitivity & Conditioning\n",
"Supplemental Julia Notebook\n",
"\n",
"Prof. Doug James"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Vector p-norms in Julia"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"search: \u001b[1mn\u001b[22m\u001b[1mo\u001b[22m\u001b[1mr\u001b[22m\u001b[1mm\u001b[22m \u001b[1mn\u001b[22m\u001b[1mo\u001b[22m\u001b[1mr\u001b[22m\u001b[1mm\u001b[22mpath \u001b[1mn\u001b[22m\u001b[1mo\u001b[22m\u001b[1mr\u001b[22m\u001b[1mm\u001b[22malize \u001b[1mn\u001b[22m\u001b[1mo\u001b[22m\u001b[1mr\u001b[22m\u001b[1mm\u001b[22malize! \u001b[1mn\u001b[22m\u001b[1mo\u001b[22m\u001b[1mr\u001b[22m\u001b[1mm\u001b[22malize_string vec\u001b[1mn\u001b[22m\u001b[1mo\u001b[22m\u001b[1mr\u001b[22m\u001b[1mm\u001b[22m issub\u001b[1mn\u001b[22m\u001b[1mo\u001b[22m\u001b[1mr\u001b[22m\u001b[1mm\u001b[22mal\n",
"\n"
]
},
{
"data": {
"text/markdown": [
"```\n",
"norm(A, [p])\n",
"```\n",
"\n",
"Compute the `p`-norm of a vector or the operator norm of a matrix `A`, defaulting to the `p=2`-norm.\n",
"\n",
"For vectors, `p` can assume any numeric value (even though not all values produce a mathematically valid vector norm). In particular, `norm(A, Inf)` returns the largest value in `abs(A)`, whereas `norm(A, -Inf)` returns the smallest.\n",
"\n",
"For matrices, the matrix norm induced by the vector `p`-norm is used, where valid values of `p` are `1`, `2`, or `Inf`. (Note that for sparse matrices, `p=2` is currently not implemented.) Use [`vecnorm`](:func:`vecnorm`) to compute the Frobenius norm.\n"
],
"text/plain": [
"```\n",
"norm(A, [p])\n",
"```\n",
"\n",
"Compute the `p`-norm of a vector or the operator norm of a matrix `A`, defaulting to the `p=2`-norm.\n",
"\n",
"For vectors, `p` can assume any numeric value (even though not all values produce a mathematically valid vector norm). In particular, `norm(A, Inf)` returns the largest value in `abs(A)`, whereas `norm(A, -Inf)` returns the smallest.\n",
"\n",
"For matrices, the matrix norm induced by the vector `p`-norm is used, where valid values of `p` are `1`, `2`, or `Inf`. (Note that for sparse matrices, `p=2` is currently not implemented.) Use [`vecnorm`](:func:`vecnorm`) to compute the Frobenius norm.\n"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"?(norm)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Float64,1}:\n",
" 1.0\n",
" 1.0\n",
" 1.0"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a=ones(3)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1.7320508075688772"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"norm(a)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1.7320508075688772"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"norm(a,2)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1.0"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"norm(a,Inf)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3.0"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"norm(a,1)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Visualizing Vector p-norms\n",
"Let us look at contour plots of \n",
"\\begin{equation}\n",
" \\| (x,y) \\|_p\n",
"\\end{equation}\n",
"for different values of $p$."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using Plots"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"f (generic function with 1 method)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(x,y) = vecnorm([x y], p) # parameter-dependent p-norm function"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"x = y = linspace(-1, 1, 100);"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p=1; \n",
"contourf(x, y, f, xlabel=\"x\", ylabel=\"y\", title=string(\"Norm, p=\",p))"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p=2; \n",
"contourf(x, y, f, xlabel=\"x\", ylabel=\"y\", title=string(\"Norm, p=\",p))"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p=5; \n",
"contourf(x, y, f, xlabel=\"x\", ylabel=\"y\", title=string(\"Norm, p=\",p))"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p=Inf; \n",
"contourf(x, y, f, xlabel=\"x\", ylabel=\"y\", title=string(\"Norm, p=\",p))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Matrix-transforming vectors on a circle"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: Method definition unitVec(Any) in module Main at In[14]:1 overwritten at In[15]:1.\n"
]
},
{
"data": {
"text/plain": [
"unitVec (generic function with 1 method)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"unitVec(θ)=[cos(θ), sin(θ)] # Unit vector (in 2-norm) in direction θ"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 1.0\n",
" 0.0"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"unitVec(0)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"2×360 Array{Float64,2}:\n",
" 0.949766 0.80411 0.577666 0.293185 … 0.654252 0.3847 0.076498\n",
" 0.312962 0.594481 0.816273 0.956056 0.756277 0.923042 0.99707 "
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"UV=zeros(2,360) # Matrix of unit 2-vectors as columns\n",
"for i=1:360\n",
" UV[:,i] = unitVec(i/π)\n",
"end\n",
"UV"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"R = randn(2,2) # randn matrix\n",
"B = R * UV # Multiply UV by random matrix\n",
"plot(B[1,:], B[2,:], title=string(\"RandnMtx * UnitCircle, ||R||_2=\",norm(R,2)))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## 2D Visualization of Induced Matrix Norms\n",
"\n",
"For a 2D unit vector $\\vec{x}(\\theta) = (\\cos\\theta, \\sin\\theta)$, and a random matrix $\\mathbf{R}$, let's visualize the induced matrix $p$-norm ratios\n",
"\\begin{equation}\n",
" ratio(\\theta)=\\frac{\\| \\mathbf{R} \\; \\vec{x}(\\theta) \\|_p}{\\|\\vec{x}(\\theta)\\|_p}\n",
"\\end{equation}\n",
"for different values of $p$. "
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"nRUVRatio (generic function with 1 method)"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nRUVRatio(θ) = norm(R*unitVec(θ), p)/norm(unitVec(θ),p) # p-norm of R*unitVec / p-norm of unitVec"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"R = [0.504342 0.139108; -0.0371534 -0.0523837]\n"
]
}
],
"source": [
"R=randn(2,2) # a randn matrix to measure\n",
"println(\"R = \",R)"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"||R||_Inf = 0.6434504854188444\n"
]
},
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p=Inf;\n",
"println(\"||R||_\",p,\" = \",norm(R,p))\n",
"plot(nRUVRatio, linspace(0,π,360), xlabel=\"θ\", title=string(\"ratio(θ) for p=\",p))"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"||R||_2 = 0.5255487353898564\n"
]
},
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p=2;\n",
"println(\"||R||_\",p,\" = \",norm(R,p))\n",
"plot(nRUVRatio, linspace(0,π,360), xlabel=\"θ\", title=string(\"ratio(θ) for p=\",p))"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"||R||_1 = 0.5414954937506968\n"
]
},
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p=1;\n",
"println(\"||R||_\",p,\" = \",norm(R,p))\n",
"plot(nRUVRatio, linspace(0,π,360), xlabel=\"θ\", title=string(\"ratio(θ) for p=\",p))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Experiments with an ill-conditioned Vandermonde matrix"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"vander (generic function with 1 method)"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# A famous ill-conditioned matrix: Vandermonde matrix\n",
"vander(N) = [x^i for x=linspace(1./N,1,N), i=1:N] # Vandermonde matrix"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Array{Float64,2}:\n",
" 0.25 0.0625 0.015625 0.00390625\n",
" 0.5 0.25 0.125 0.0625 \n",
" 0.75 0.5625 0.421875 0.316406 \n",
" 1.0 1.0 1.0 1.0 "
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N = 4\n",
"V = vander(N)"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"size = (4,4)\n"
]
},
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"println(\"size = \",size(V))\n",
"plot(V, title=string(\"Vandermonde matrix. size=\",size(V)))"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"524.7924487965224"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cond(V) # Condition number (expensive computation but...)"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4-element Array{Float64,1}:\n",
" 2.3134 \n",
" 0.447495 \n",
" 0.0601855 \n",
" 0.00440821"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"svdvals(V)"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4-element Array{Float64,1}:\n",
" 1.0\n",
" 1.0\n",
" 1.0\n",
" 1.0"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b=ones(N)"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4-element Array{Float64,1}:\n",
" 8.33333\n",
" -23.3333 \n",
" 26.6667 \n",
" -10.6667 "
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"LU = lufact(V)\n",
"x=LU\\b"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"8.899114524108741e-16"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vecnorm(b-V*x)/vecnorm(b) # relative error of matrix solve"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## How does solve rel-error depend on size of Vandermonde matrix?"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"vanderSolveRelErr (generic function with 1 method)"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function vanderSolveRelErr(N)\n",
" V = [x^i for x=linspace(1.0/N,1.0,N), i=1:N]\n",
" b = ones(N) # a predictable RHS, \\vec{1}\n",
" LU = lufact(V)\n",
" x = LU \\ b\n",
" return vecnorm(b-V*x)/vecnorm(b)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"100-element Array{Float64,1}:\n",
" 0.0 \n",
" 0.0 \n",
" 2.56395e-16\n",
" 8.89911e-16\n",
" 4.19071e-15\n",
" 1.07055e-14\n",
" 3.29288e-14\n",
" 1.09224e-13\n",
" 3.2904e-13 \n",
" 2.90792e-12\n",
" 1.02942e-11\n",
" 4.54925e-11\n",
" 1.48135e-10\n",
" ⋮ \n",
" 0.0906789 \n",
" 1.94849 \n",
" 0.198944 \n",
" 0.363464 \n",
" 3.10338 \n",
" 0.321879 \n",
" 0.299673 \n",
" 0.614996 \n",
" 0.0526596 \n",
" 0.0704173 \n",
" 0.0933069 \n",
" 0.306217 "
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vErr = [vanderSolveRelErr(n) for n=1:100] # solve all problems from size n=1 to 100"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(log10.(vErr+eps()), title=\"RelError solving V x = 1\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Testing $V^{-1} V = I$:"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n=2: relErr(invV*V-I): 0.0\n",
"n=4: relErr(invV*V-I): 5.704378753739517e-15\n",
"n=6: relErr(invV*V-I): 7.206346992435156e-13\n",
"n=8: relErr(invV*V-I): 4.232127554641132e-11\n",
"n=10: relErr(invV*V-I): 3.494164228312279e-9\n",
"n=12: relErr(invV*V-I): 2.67328213816742e-7\n",
"n=14: relErr(invV*V-I): 1.13656866037886e-5\n",
"n=16: relErr(invV*V-I): 0.0005713609123286228\n",
"n=18: relErr(invV*V-I): 0.03460692390723406\n",
"n=20: relErr(invV*V-I): 2.7953035818633856\n",
"n=22: relErr(invV*V-I): 32.06809436167812\n",
"n=24: relErr(invV*V-I): 2981.979751536043\n",
"n=26: relErr(invV*V-I): 74.36798538478568\n",
"n=28: relErr(invV*V-I): 146.86182945261046\n",
"n=30: relErr(invV*V-I): 61.98487651414396\n",
"n=32: relErr(invV*V-I): 632.6816419885412\n",
"n=34: relErr(invV*V-I): 980.6993965209957\n",
"n=36: relErr(invV*V-I): 86.15079649282958\n",
"n=38: relErr(invV*V-I): 205.17510463084815\n",
"n=40: relErr(invV*V-I): 1207.7221166111324\n",
"n=42: relErr(invV*V-I): 215.92511095558865\n",
"n=44: relErr(invV*V-I): 1380.6409470134085\n",
"n=46: relErr(invV*V-I): 170.36879116352844\n",
"n=48: relErr(invV*V-I): 867.0948304744451\n",
"n=50: relErr(invV*V-I): 286.17763388790155\n"
]
}
],
"source": [
"for n=2:2:50\n",
" println(\"n=\",n,\": relErr(invV*V-I): \",norm(inv(vander(n))*vander(n)-eye(n))/norm(eye(n)))\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"PROBLEM: Vandermonde matrix has many near-zero singular values for large $N$."
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"V = vander(100)\n",
"plot(log10.(svdvals(V)))"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"-0.0"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"det(V) # numerically singular"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"5.205505037741025e19"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cond(V) # awful condition number"
]
},
{
"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
}