{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Stanford CS205a \n",
"\n",
"## Eigenvalues & SVD\n",
"Supplemental Julia Notebook\n",
"\n",
"Prof. Doug James"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Matrix inverse & eigs"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Array{Float64,2}:\n",
" 0.838168 0.894378 0.860048 \n",
" 0.560628 0.311312 0.287067 \n",
" 0.155911 0.821748 0.0381491"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = rand(3,3)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Complex{Float64},1}:\n",
" 1.61149+0.0im \n",
" -0.211928+0.267408im\n",
" -0.211928-0.267408im"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eigA = eigvals(A)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Complex{Float64},1}:\n",
" 0.620545-0.0im \n",
" -1.82037-2.29692im\n",
" -1.82037+2.29692im"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# 1/eigenvalues match...\n",
"1./eigA"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Complex{Float64},1}:\n",
" 0.620545+0.0im \n",
" -1.82037+2.29692im\n",
" -1.82037-2.29692im"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# ... the eigenvalues of the inverse(A)\n",
"eigvals(inv(A))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Matrix square root"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Array{Float64,2}:\n",
" 0.0968649 0.139353 0.331046\n",
" 0.139353 0.507817 0.728221\n",
" 0.331046 0.728221 1.53727 "
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"## Symmetric positive definite matrix:\n",
"B = rand(3,3);\n",
"A = B'*B + 0.0001*eye(3)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Array{Float64,2}:\n",
" 0.204495 0.0575041 0.227464\n",
" 0.0575041 0.576255 0.415259\n",
" 0.227464 0.415259 1.1459 "
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Has a matrix square root\n",
"sqrtA = sqrtm(A)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Float64,1}:\n",
" 0.147315\n",
" 0.371388\n",
" 1.40795 "
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Notice that the eigenvalues of sqrt(A)...\n",
"eigvals(sqrtA)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Float64,1}:\n",
" 0.147315\n",
" 0.371388\n",
" 1.40795 "
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Are equal to the sqrt of the eigs of A:\n",
"sqrt.(eigvals(A))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Polar Decomposition"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Array{Float64,2}:\n",
" 0.184173 0.0912081 0.749521\n",
" 0.383199 0.643191 0.870778\n",
" 0.667197 0.195219 0.323629"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Given A, find R s.t. A = R U:\n",
"A = rand(3,3)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Array{Float64,2}:\n",
" 0.678215 0.218509 0.343789\n",
" 0.218509 0.51118 0.388681\n",
" 0.343789 0.388681 1.07495 "
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Find U first, then R:\n",
"U2 = A'*A;\n",
"U = sqrtm(U2) # sqrt root of U*U"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Array{Float64,2}:\n",
" -0.0180764 -0.480526 0.876794 \n",
" 0.0390546 0.875929 0.480856 \n",
" 0.999074 -0.042935 -0.00293309"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Now A=RU --> R=A invU\n",
"R = A*inv(U) # not efficient algorithm (for now), but true"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Array{Float64,2}:\n",
" 1.0 1.63064e-15 3.33934e-17\n",
" 1.63064e-15 1.0 -1.00546e-15\n",
" 3.33934e-17 -1.00546e-15 1.0 "
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Verify R is orthogonal:\n",
"R'*R"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1.3207184582579895e-16"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Factorization A=RU is accurate:\n",
"norm(A-R*U)/norm(A)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Thin SVD using `svd()`"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"search: \u001b[1ms\u001b[22m\u001b[1mv\u001b[22m\u001b[1md\u001b[22m \u001b[1ms\u001b[22m\u001b[1mv\u001b[22m\u001b[1md\u001b[22ms \u001b[1ms\u001b[22m\u001b[1mv\u001b[22m\u001b[1md\u001b[22mvals \u001b[1ms\u001b[22m\u001b[1mv\u001b[22m\u001b[1md\u001b[22mfact \u001b[1ms\u001b[22m\u001b[1mv\u001b[22m\u001b[1md\u001b[22mvals! \u001b[1ms\u001b[22m\u001b[1mv\u001b[22m\u001b[1md\u001b[22mfact! i\u001b[1ms\u001b[22m\u001b[1mv\u001b[22mali\u001b[1md\u001b[22m\n",
"\n"
]
}
],
"source": [
"?(svd)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"4×2 Array{Float64,2}:\n",
" 0.537532 0.575281\n",
" 0.792717 0.185602\n",
" 0.699253 0.996808\n",
" 0.791883 0.72888 "
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = rand(4,2)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(\n",
"[-0.409043 0.0895816; -0.364472 -0.865563; -0.621606 0.490516; -0.559866 -0.0465768],\n",
"\n",
"[1.92195,0.479357],\n",
"[-0.721563 -0.692349; -0.692349 0.721563])"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Thin SVD:\n",
"(U,S,V) = svd(A)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4×2 Array{Float64,2}:\n",
" -0.409043 0.0895816\n",
" -0.364472 -0.865563 \n",
" -0.621606 0.490516 \n",
" -0.559866 -0.0465768"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"U"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 1.92195 \n",
" 0.479357"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"S"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Array{Float64,2}:\n",
" 1.92195 0.0 \n",
" 0.0 0.479357"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diagm(S)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Array{Float64,2}:\n",
" -0.721563 -0.692349\n",
" -0.692349 0.721563"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"V'"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Verifying `A = V*diagm(S)*V'`"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4×2 Array{Float64,2}:\n",
" 0.537532 0.575281\n",
" 0.792717 0.185602\n",
" 0.699253 0.996808\n",
" 0.791883 0.72888 "
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4×2 Array{Float64,2}:\n",
" 0.537532 0.575281\n",
" 0.792717 0.185602\n",
" 0.699253 0.996808\n",
" 0.791883 0.72888 "
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"U*diagm(S)*V'"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4.32189512827081e-16"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"## Relative Error:\n",
"norm(U*diagm(S)*V'-A)/norm(A)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Vandermonde Matrix & SVD"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"vander (generic function with 1 method)"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function vander(N)\n",
" V = [(i/N)^(j-1) for i=1:N, j=1:N];\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"10×10 Array{Float64,2}:\n",
" 1.0 0.1 0.01 0.001 0.0001 … 1.0e-7 1.0e-8 1.0e-9 \n",
" 1.0 0.2 0.04 0.008 0.0016 1.28e-5 2.56e-6 5.12e-7 \n",
" 1.0 0.3 0.09 0.027 0.0081 0.0002187 6.561e-5 1.9683e-5 \n",
" 1.0 0.4 0.16 0.064 0.0256 0.0016384 0.00065536 0.000262144\n",
" 1.0 0.5 0.25 0.125 0.0625 0.0078125 0.00390625 0.00195313 \n",
" 1.0 0.6 0.36 0.216 0.1296 … 0.0279936 0.0167962 0.0100777 \n",
" 1.0 0.7 0.49 0.343 0.2401 0.0823543 0.057648 0.0403536 \n",
" 1.0 0.8 0.64 0.512 0.4096 0.209715 0.167772 0.134218 \n",
" 1.0 0.9 0.81 0.729 0.6561 0.478297 0.430467 0.38742 \n",
" 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 "
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N = 10;\n",
"V = vander(N)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"10-element Array{Float64,1}:\n",
" 4.65476 \n",
" 2.10904 \n",
" 0.650574 \n",
" 0.16112 \n",
" 0.0325356 \n",
" 0.00526855 \n",
" 0.000656848\n",
" 5.8534e-5 \n",
" 3.25623e-6 \n",
" 8.30914e-8 "
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sv = svdvals(V)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4.654763958770751"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"## Matrix two-norm is largest singval:\n",
"norm(V)"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4.654763958770751"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sv[1]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Condition number in singular values:"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"5.6019784622350566e7"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cond(V)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"5.6019784622350566e7"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"## Definition of condition number in singular values:\n",
"sv[1]/sv[length(sv)]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Plotting singular values of Vandermonde matrix"
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "-"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 61,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots\n",
"N = 1000;\n",
"V = vander(N) / norm(vander(N)); ## normalize V for clarity\n",
"sv = svdvals(V);\n",
"plot(log10.(sv))"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Tikhonov Regularization & SVD"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A0 = V'*V; \n",
"alpha = 1.e-10; \n",
"A = A0 + alpha*eye(N);\n",
"sv0 = svdvals(A0);\n",
"sv = svdvals(A);\n",
"plot(log10.([sv0, sv]))"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Augmented matrix approach to regularization"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"A = [V; sqrt(alpha)*eye(N)];"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 64,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(log10.([svdvals(V), svdvals(A)]))"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 65,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(log10.([sv0, svdvals(A'*A)]))"
]
},
{
"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
}