{
"cells": [
{
"cell_type": "markdown",
"id": "4159be0e",
"metadata": {
"hideCode": false,
"hidePrompt": false,
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 2023-03-01 Singular Value Decomposition\n",
"\n",
"## Last time\n",
"\n",
"* Householder QR\n",
"* Composition of reflectors\n",
"\n",
"## Today\n",
"\n",
"* Comparison of interfaces\n",
"* Profiling\n",
"* Cholesky QR\n",
"* Matrix norms and conditioning\n",
"* Geometry of the SVD"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "eb781a7a",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"data": {
"text/plain": [
"qr_householder (generic function with 1 method)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using LinearAlgebra\n",
"using Plots\n",
"default(linewidth=4, legendfontsize=12)\n",
"\n",
"function vander(x, k=nothing)\n",
" if isnothing(k)\n",
" k = length(x)\n",
" end\n",
" m = length(x)\n",
" V = ones(m, k)\n",
" for j in 2:k\n",
" V[:, j] = V[:, j-1] .* x\n",
" end\n",
" V\n",
"end\n",
"\n",
"function gram_schmidt_classical(A)\n",
" m, n = size(A)\n",
" Q = zeros(m, n)\n",
" R = zeros(n, n)\n",
" for j in 1:n\n",
" v = A[:,j]\n",
" R[1:j-1,j] = Q[:,1:j-1]' * v\n",
" v -= Q[:,1:j-1] * R[1:j-1,j]\n",
" R[j,j] = norm(v)\n",
" Q[:,j] = v / R[j,j]\n",
" end\n",
" Q, R\n",
"end\n",
"\n",
"function qr_householder(A)\n",
" m, n = size(A)\n",
" R = copy(A)\n",
" V = [] # list of reflectors\n",
" for j in 1:n\n",
" v = copy(R[j:end, j])\n",
" v[1] += sign(v[1]) * norm(v) # <---\n",
" v = normalize(v)\n",
" R[j:end,j:end] -= 2 * v * v' * R[j:end,j:end]\n",
" push!(V, v)\n",
" end\n",
" V, R\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "0df4de0d",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Composition of reflectors\n",
"\n",
"\\begin{align}\n",
"(I - 2 v v^T) (I - 2 w w^T) &= I - 2 v v^T - 2 w w^T + 4 v (v^T w) w^T \\\\\n",
"&= I - \\Bigg[v \\Bigg| w \\Bigg] \\begin{bmatrix} 2 & -4 v^T w \\\\ 0 & 2 \\end{bmatrix} \\begin{bmatrix} v^T \\\\ w^T \\end{bmatrix}\n",
"\\end{align}\n",
"\n",
"This turns applying reflectors from a sequence of vector operations to a sequence of (smallish) matrix operations. It's the key to high performance and the native format (`QRCompactWY`) returned by Julia `qr()`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "c21e1e5e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}}\n",
"Q factor:\n",
"20×20 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}}:\n",
" -0.223607 -0.368394 -0.430192 0.437609 … -3.23545e-5 -5.31905e-6\n",
" -0.223607 -0.329616 -0.294342 0.161225 0.000550027 0.000101062\n",
" -0.223607 -0.290838 -0.173586 -0.0383868 -0.00436786 -0.000909558\n",
" -0.223607 -0.252059 -0.067925 -0.170257 0.0214511 0.00515416\n",
" -0.223607 -0.213281 0.0226417 -0.243417 -0.0726036 -0.0206167\n",
" -0.223607 -0.174503 0.0981139 -0.266901 … 0.178209 0.06185\n",
" -0.223607 -0.135724 0.158492 -0.24974 -0.323416 -0.144317\n",
" -0.223607 -0.0969458 0.203775 -0.200966 0.429021 0.268016\n",
" -0.223607 -0.0581675 0.233964 -0.129612 -0.386119 -0.402025\n",
" -0.223607 -0.0193892 0.249058 -0.0447093 0.157308 0.491364\n",
" -0.223607 0.0193892 0.249058 0.0447093 … 0.157308 -0.491364\n",
" -0.223607 0.0581675 0.233964 0.129612 -0.386119 0.402025\n",
" -0.223607 0.0969458 0.203775 0.200966 0.429021 -0.268016\n",
" -0.223607 0.135724 0.158492 0.24974 -0.323416 0.144317\n",
" -0.223607 0.174503 0.0981139 0.266901 0.178209 -0.06185\n",
" -0.223607 0.213281 0.0226417 0.243417 … -0.0726036 0.0206167\n",
" -0.223607 0.252059 -0.067925 0.170257 0.0214511 -0.00515416\n",
" -0.223607 0.290838 -0.173586 0.0383868 -0.00436786 0.000909558\n",
" -0.223607 0.329616 -0.294342 -0.161225 0.000550027 -0.000101062\n",
" -0.223607 0.368394 -0.430192 -0.437609 -3.23545e-5 5.31905e-6\n",
"R factor:\n",
"20×20 Matrix{Float64}:\n",
" -4.47214 0.0 -1.64763 0.0 … -0.514468 2.22045e-16\n",
" 0.0 2.71448 1.11022e-16 1.79412 -2.498e-16 0.823354\n",
" 0.0 0.0 -1.46813 5.55112e-17 -0.944961 -2.23779e-16\n",
" 0.0 0.0 0.0 -0.774796 3.83808e-17 -0.913056\n",
" 0.0 0.0 0.0 0.0 0.797217 -4.06264e-16\n",
" 0.0 0.0 0.0 0.0 … -3.59496e-16 0.637796\n",
" 0.0 0.0 0.0 0.0 -0.455484 -1.3936e-15\n",
" 0.0 0.0 0.0 0.0 4.40958e-16 -0.313652\n",
" 0.0 0.0 0.0 0.0 -0.183132 1.64685e-15\n",
" 0.0 0.0 0.0 0.0 4.82253e-16 0.109523\n",
" 0.0 0.0 0.0 0.0 … 0.0510878 5.9848e-16\n",
" 0.0 0.0 0.0 0.0 -2.68709e-15 0.0264553\n",
" 0.0 0.0 0.0 0.0 -0.0094344 -2.94383e-15\n",
" 0.0 0.0 0.0 0.0 2.08514e-15 0.00417208\n",
" 0.0 0.0 0.0 0.0 0.0010525 -2.24994e-15\n",
" 0.0 0.0 0.0 0.0 … -1.64363e-15 -0.000385264\n",
" 0.0 0.0 0.0 0.0 -5.9057e-5 7.69025e-16\n",
" 0.0 0.0 0.0 0.0 1.76642e-16 -1.66202e-5\n",
" 0.0 0.0 0.0 0.0 -1.04299e-6 -1.68771e-16\n",
" 0.0 0.0 0.0 0.0 0.0 1.71467e-7"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = LinRange(-1, 1, 20)\n",
"A = vander(x)\n",
"Q, R = qr(A)\n",
"#norm(Q[:,1])"
]
},
{
"cell_type": "markdown",
"id": "57039d94",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# This works even for very nonsquare matrices"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "790bb615",
"metadata": {
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"size(Q) = (1000000, 1000000)\n",
"norm(Q * R - A) = 1.5627435002979355e-12\n"
]
},
{
"data": {
"text/plain": [
"1000000-element Vector{Float64}:\n",
" -0.000710716962490326\n",
" -0.0006221486003295463\n",
" -0.001481395515703071\n",
" 0.0006986145984357538\n",
" 0.00030719398438148966\n",
" -3.7070273376922167e-6\n",
" -1.846993064721867e-6\n",
" 5.892273227768894e-7\n",
" -3.4792383435730715e-6\n",
" 2.2249702781554393e-6\n",
" -1.8148213625435382e-6\n",
" -1.3842243836198498e-8\n",
" -2.943731575863994e-6\n",
" ⋮\n",
" 5.4302791596705405e-8\n",
" -1.2863812388719472e-6\n",
" -3.1178852037441024e-6\n",
" -1.454678366309665e-6\n",
" 1.129578508676748e-6\n",
" -1.0958261567418687e-6\n",
" -1.4438579316234545e-7\n",
" 4.451050785148942e-7\n",
" -1.2256371019533118e-6\n",
" 0.9999963354732083\n",
" 4.843751131958308e-7\n",
" 2.761092851686931e-7"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = rand(1000000, 5)\n",
"Q, R = qr(A)\n",
"@show size(Q)\n",
"@show norm(Q*R - A)\n",
"R\n",
"Q * vcat(zeros(1000000 - 3), [1,0,0])"
]
},
{
"cell_type": "markdown",
"id": "cd14dfaf",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"This is known as a \"full\" (or \"complete\") QR factorization, in contrast to a reduced QR factorization in which $Q$ has the same shape as $A$.\n",
"* How much memory does $Q$ use?"
]
},
{
"cell_type": "markdown",
"id": "8aed5e91",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Compare to [`numpy.linalg.qr`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.qr.html)\n",
"\n",
"* Need to decide up-front whether you want full or reduced QR.\n",
"* Full QR is expensive to represent."
]
},
{
"cell_type": "markdown",
"id": "7ec2f869",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Cholesky QR\n",
"\n",
"$$ R^T R = (QR)^T QR = A^T A$$\n",
"\n",
"so we should be able to use $L L^T = A^T A$ and then $Q = A L^{-T}$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "d8bba8a1",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"norm(Q' * Q - I) = 1.3547676815713508e-13\n",
"norm(Q * R - A) = 1.289504394296693e-15\n"
]
},
{
"data": {
"text/plain": [
"1.289504394296693e-15"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function qr_chol(A)\n",
" R = cholesky(A' * A).U\n",
" Q = A / R\n",
" Q, R\n",
"end\n",
"\n",
"A = rand(20,20)\n",
"Q, R = qr_chol(A)\n",
"@show norm(Q' * Q - I)\n",
"@show norm(Q * R - A)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "5b4167c8",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"ename": "LoadError",
"evalue": "PosDefException: matrix is not positive definite; Cholesky factorization failed.",
"output_type": "error",
"traceback": [
"PosDefException: matrix is not positive definite; Cholesky factorization failed.",
"",
"Stacktrace:",
" [1] checkpositivedefinite",
" @ /usr/share/julia/stdlib/v1.7/LinearAlgebra/src/factorization.jl:18 [inlined]",
" [2] cholesky!(A::Hermitian{Float64, Matrix{Float64}}, ::Val{false}; check::Bool)",
" @ LinearAlgebra /usr/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:266",
" [3] cholesky!(A::Matrix{Float64}, ::Val{false}; check::Bool)",
" @ LinearAlgebra /usr/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:298",
" [4] #cholesky#143",
" @ /usr/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:394 [inlined]",
" [5] cholesky (repeats 2 times)",
" @ /usr/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:394 [inlined]",
" [6] qr_chol(A::Matrix{Float64})",
" @ Main ./In[9]:2",
" [7] top-level scope",
" @ In[14]:3",
" [8] eval",
" @ ./boot.jl:373 [inlined]",
" [9] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
" @ Base ./loading.jl:1196"
]
}
],
"source": [
"x = LinRange(-1, 1, 20)\n",
"A = vander(x)\n",
"Q, R = qr_chol(A)\n",
"@show norm(Q' * Q - I)\n",
"@show norm(Q * R - A);"
]
},
{
"cell_type": "markdown",
"id": "bc52e84f",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Can we fix this?\n",
"\n",
"Note that the product of two triangular matrices is triangular."
]
},
{
"cell_type": "code",
"execution_count": 43,
"id": "640aa2ec",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"5×5 Matrix{Float64}:\n",
" 0.937018 0.166589 1.28394 0.856513 1.91001\n",
" 0.0 0.147171 0.422197 0.410025 1.03557\n",
" 0.0 0.0 0.466456 0.381229 0.947622\n",
" 0.0 0.0 0.0 0.690276 0.506928\n",
" 0.0 0.0 0.0 0.0 0.0989786"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"R = triu(rand(5,5))\n",
"R * R"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "b81eb674",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"norm(Q' * Q - I) = 1.1333031064724956e-15\n",
"norm(Q * R - A) = 1.1896922215066695e-15\n"
]
}
],
"source": [
"function qr_chol2(A)\n",
" Q, R = qr_chol(A)\n",
" Q, R1 = qr_chol(Q)\n",
" Q, R1 * R\n",
"end\n",
"\n",
"x = LinRange(-1, 1, 18)\n",
"A = vander(x)\n",
"Q, R = qr_chol2(A)\n",
"@show norm(Q' * Q - I)\n",
"@show norm(Q * R - A);"
]
},
{
"cell_type": "markdown",
"id": "d845fd38",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# How fast are these methods?"
]
},
{
"cell_type": "code",
"execution_count": 34,
"id": "ce4c6d9a",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.507001 seconds (7 allocations: 77.393 MiB, 0.34% gc time)\n"
]
}
],
"source": [
"m, n = 5000, 2000\n",
"A = randn(m, n)\n",
"\n",
"@time qr(A);"
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "2bbae91e",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.067684 seconds (21 allocations: 305.176 MiB, 0.90% gc time)\n"
]
}
],
"source": [
"A = randn(m, n)\n",
"@time qr_chol2(A);"
]
},
{
"cell_type": "markdown",
"id": "4c375979",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Profiling"
]
},
{
"cell_type": "code",
"execution_count": 37,
"id": "8c9fcd68",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
"ProfileSVG.FGConfig(Node(FlameGraphs.NodeData(ip:0x0, 0x00, 1:2139)), Dict{Symbol, Any}(), FlameGraphs.FlameColors(RGB{FixedPointNumbers.N0f8}[RGB{N0f8}(0.882,0.698,1.0), RGB{N0f8}(0.435,0.863,0.569), RGB{N0f8}(0.0,0.71,0.545), RGB{N0f8}(0.173,0.639,1.0)], RGB{N0f8}(1.0,1.0,1.0), RGB{N0f8}(0.0,0.0,0.0), RGB{FixedPointNumbers.N0f8}[RGB{N0f8}(0.953,0.0,0.302), RGB{N0f8}(0.894,0.0,0.255), RGB{N0f8}(0.831,0.129,0.216), RGB{N0f8}(0.773,0.192,0.184)], RGB{FixedPointNumbers.N0f8}[RGB{N0f8}(1.0,0.627,0.0), RGB{N0f8}(1.0,0.643,0.0), RGB{N0f8}(0.965,0.651,0.039), RGB{N0f8}(0.894,0.655,0.11)]), :fcolor, :fcolor, 1.0, false, 50, 2000, 960.0, 0.0, 2.0, \"inherit\", 12.0, false, :none, 0.001)"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using ProfileSVG\n",
"\n",
"@profview qr(A)"
]
},
{
"cell_type": "markdown",
"id": "1a6f0f24",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Condition number of a matrix\n",
"\n",
"We may have informally referred to a matrix as \"ill-conditioned\" when the columns are nearly linearly dependent, but let's make this concept for precise. Recall the definition of (relative) condition number:\n",
"\n",
"$$ \\kappa = \\max_{\\delta x} \\frac{|\\delta f|/|f|}{|\\delta x|/|x|} . $$\n",
"\n",
"We understood this definition for scalar problems, but it also makes sense when the inputs and/or outputs are vectors (or matrices, etc.) and absolute value is replaced by vector (or matrix) norms. Consider matrix-vector multiplication, for which $f(x) = A x$.\n",
"\n",
"$$ \\kappa(A) = \\max_{\\delta x} \\frac{\\lVert A (x+\\delta x) - A x \\rVert/\\lVert A x \\rVert}{\\lVert \\delta x\\rVert/\\lVert x \\rVert}\n",
"= \\max_{\\delta x} \\frac{\\lVert A \\delta x \\rVert}{\\lVert \\delta x \\rVert} \\, \\frac{\\lVert x \\rVert}{\\lVert A x \\rVert} = \\lVert A \\rVert \\frac{\\lVert x \\rVert}{\\lVert A x \\rVert} . $$\n",
"\n",
"There are two problems here:\n",
"\n",
"* I wrote $\\kappa(A)$ but my formula depends on $x$.\n",
"* What is that $\\lVert A \\rVert$ beastie?\n"
]
},
{
"cell_type": "markdown",
"id": "4028a14e",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Matrix norms induced by vector norms\n",
"\n",
"Vector norms are built into the linear space (and defined in term of the inner product). Matrix norms are *induced* by vector norms, according to\n",
"\n",
"$$ \\lVert A \\rVert = \\max_{x \\ne 0} \\frac{\\lVert A x \\rVert}{\\lVert x \\rVert} . $$\n",
"\n",
"* This equation makes sense for non-square matrices -- the vector norms of the input and output spaces may differ.\n",
"* Due to linearity, all that matters is direction of $x$, so it could equivalently be written\n",
"\n",
"$$ \\lVert A \\rVert = \\max_{\\lVert x \\rVert = 1} \\lVert A x \\rVert . $$\n"
]
},
{
"cell_type": "markdown",
"id": "5c953dc0",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# The formula makes sense, but still depends on $x$.\n",
"\n",
"$$\\kappa(A) = \\lVert A \\rVert \\frac{\\lVert x \\rVert}{\\lVert Ax \\rVert}$$\n",
"\n",
"Consider the matrix\n",
"\n",
"$$ A = \\begin{bmatrix} 1 & 0 \\\\ 0 & 0 \\end{bmatrix} . $$\n",
"\n",
"* What is the norm of this matrix?\n",
"* What is the condition number when $x = [1,0]^T$?\n",
"* What is the condition number when $x = [0,1]^T$?"
]
},
{
"cell_type": "markdown",
"id": "7d148680",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Condition number of the matrix\n",
"\n",
"The condition number of matrix-vector multiplication depends on the vector. The condition number of the matrix is the worst case (maximum) of the condition number for any vector, i.e.,\n",
"\n",
"$$ \\kappa(A) = \\max_{x \\ne 0} \\lVert A \\rVert \\frac{\\lVert x \\rVert}{\\lVert A x \\rVert} .$$\n",
"\n",
"If $A$ is invertible, then we can rephrase as\n",
"\n",
"$$ \\kappa(A) = \\max_{x \\ne 0} \\lVert A \\rVert \\frac{\\lVert A^{-1} (A x) \\rVert}{\\lVert A x \\rVert} =\n",
"\\max_{A x \\ne 0} \\lVert A \\rVert \\frac{\\lVert A^{-1} (A x) \\rVert}{\\lVert A x \\rVert} = \\lVert A \\rVert \\lVert A^{-1} \\rVert . $$\n",
"\n",
"Evidently multiplying by a matrix is just as ill-conditioned of an operation as solving a linear system using that matrix."
]
},
{
"cell_type": "markdown",
"id": "80852e6a",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Least squares and the normal equations\n",
"\n",
"A **least squares problem** takes the form: given an $m\\times n$ matrix $A$ ($m \\ge n$), find $x$ such that\n",
"$$ \\lVert Ax - b \\rVert $$\n",
"is minimized. If $A$ is square and full rank, then this minimizer will satisfy $A x - b = 0$, but that is not the case in general because $b$ is not in the range of $A$.\n",
"The residual $A x - b$ must be orthogonal to the range of $A$.\n",
"\n",
"* Is this the same as saying $A^T (A x - b) = 0$?\n",
"* If $QR = A$, is it the same as $Q^T (A x - b) = 0$?\n",
"\n",
"In the quiz, we showed that $QQ^T$ is an orthogonal projector onto the range of $Q$. If $QR = A$,\n",
"$$ QQ^T (A x - b) = QQ^T(Q R x - b) = Q (Q^T Q) R x - QQ^T b = QR x - QQ^T b = A x - QQ^T b . $$\n",
"So if $b$ is in the range of $A$, we can solve $A x = b$. If not, we need only *orthogonally* project $b$ into the range of $A$."
]
},
{
"cell_type": "markdown",
"id": "d959c7a6",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Solution by QR (Householder)\n",
"\n",
"Solve $R x = Q^T b$.\n",
"\n",
"* QR factorization costs $2 m n^2 - \\frac 2 3 n^3$ operations and is done once per matrix $A$.\n",
"* Computing $Q^T b$ costs $4 (m-n)n + 2 n^2 = 4 mn - 2n^2$ (using the elementary reflectors, which are stable and lower storage than naive storage of $Q$).\n",
"* Solving with $R$ costs $n^2$ operations. Total cost per right hand side is thus $4 m n - n^2$.\n",
"\n",
"This method is stable and accurate."
]
},
{
"cell_type": "markdown",
"id": "dbdda7ac",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Solution by Cholesky\n",
"\n",
"The mathematically equivalent form $(A^T A) x = A^T b$ are called the **normal equations**. The solution process involves factoring the symmetric and positive definite $n\\times n$ matrix $A^T A$.\n",
"\n",
"* Computing $A^T A$ costs $m n^2$ flops, exploiting symmetry.\n",
"* Factoring $A^T A = R^T R$ costs $\\frac 1 3 n^3$ flops. The total factorization cost is thus $m n^2 + \\frac 1 3 n^3$.\n",
"* Computing $A^T b$ costs $2 m n$.\n",
"* Solving with $R^T$ costs $n^2$.\n",
"* Solving with $R$ costs $n^2$. Total cost per right hand side is thus $2 m n + 2 n^2$.\n",
"\n",
"The product $A^T A$ is ill-conditioned: $\\kappa(A^T A) = \\kappa(A)^2$ and can reduce the accuracy of a least squares solution."
]
},
{
"cell_type": "markdown",
"id": "a92bf79f",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Solution by Singular Value Decomposition\n",
"\n",
"Next, we will discuss a factorization\n",
"$$ U \\Sigma V^T = A $$\n",
"where $U$ and $V$ have orthonormal columns and $\\Sigma$ is diagonal with nonnegative entries.\n",
"The entries of $\\Sigma$ are called **singular values** and this decomposition is the **singular value decomposition** (SVD).\n",
"It may remind you of an eigenvalue decomposition $X \\Lambda X^{-1} = A$, but\n",
"* the SVD exists for all matrices (including non-square and deficient matrices)\n",
"* $U,V$ have orthogonal columns (while $X$ can be arbitrarily ill-conditioned).\n",
"Indeed, if a matrix is symmetric and positive definite (all positive eigenvalues), then $U=V$ and $\\Sigma = \\Lambda$.\n",
"Computing an SVD requires a somewhat complicated iterative algorithm, but a crude estimate of the cost is $2 m n^2 + 11 n^3$. Note that this is similar to the cost of $QR$ when $m \\gg n$, but much more expensive for square matrices.\n",
"Solving with the SVD involves\n",
"* Compute $U^T b$ at a cost of $2 m n$.\n",
"* Solve with the diagonal $n\\times n$ matrix $\\Sigma$ at a cost of $n$.\n",
"* Apply $V$ at a cost of $2 n^2$. The total cost per right hand side is thus $2 m n + 2n^2$."
]
},
{
"cell_type": "markdown",
"id": "7716016b",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Activity: Geometry of the Singular Value Decomposition"
]
},
{
"cell_type": "code",
"execution_count": 38,
"id": "b5209a1c",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"text/plain": [
"Aplot (generic function with 1 method)"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"default(aspect_ratio=:equal)\n",
"\n",
"function peanut()\n",
" theta = LinRange(0, 2*pi, 50)\n",
" r = 1 .+ .4*sin.(3*theta) + .6*sin.(2*theta)\n",
" r' .* [cos.(theta) sin.(theta)]'\n",
"end\n",
"\n",
"function circle()\n",
" theta = LinRange(0, 2*pi, 50)\n",
" [cos.(theta) sin.(theta)]'\n",
"end\n",
"\n",
"function Aplot(A)\n",
" \"Plot a transformation from X to Y\"\n",
" X = peanut()\n",
" Y = A * X\n",
" p = scatter(X[1,:], X[2,:], label=\"in\")\n",
" scatter!(p, Y[1,:], Y[2,:], label=\"out\")\n",
" X = circle()\n",
" Y = A * X\n",
" q = scatter(X[1,:], X[2,:], label=\"in\")\n",
" scatter!(q, Y[1,:], Y[2,:], label=\"out\")\n",
" plot(p, q, layout=2)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 39,
"id": "03165240",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Aplot(1.1*I)"
]
},
{
"cell_type": "markdown",
"id": "d8e6d03e",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Diagonal matrices\n",
"\n",
"Perhaps the simplest transformation is a scalar multiple of the identity.\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "445cfa3f",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Aplot([2 0; 0 2])"
]
},
{
"cell_type": "markdown",
"id": "3838847e",
"metadata": {
"cell_style": "split"
},
"source": [
"The diagonal entries can be different sizes.\n",
"\n",
"$$ A = \\begin{bmatrix} 2 & 0 \\\\ 0 & .5 \\end{bmatrix}$$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "dbcd9900",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Aplot([2 0; 0 .5])"
]
},
{
"cell_type": "markdown",
"id": "901eaf4f",
"metadata": {},
"source": [
"The circles becomes an **ellipse** that is **aligned with the coordinate axes**"
]
},
{
"cell_type": "markdown",
"id": "cdd57284",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Givens Rotation (as example of orthogonal matrix)\n",
"\n",
"We can rotate the input using a $2\\times 2$ matrix, parametrized by $\\theta$. Its transpose rotates in the opposite direction.\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "6502493f",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function givens(theta)\n",
" s = sin(theta)\n",
" c = cos(theta)\n",
" [c -s; s c]\n",
"end\n",
"\n",
"G = givens(0.3)\n",
"Aplot(G)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "de42c619",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Aplot(G')"
]
},
{
"cell_type": "markdown",
"id": "7e031441",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Reflection\n",
"\n",
"We've previously seen that reflectors have the form $F = I - 2 v v^T$ where $v$ is a normalized vector. Reflectors satisfy $F^T F = I$ and $F = F^T$, thus $F^2 = I$."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "e347c10a",
"metadata": {
"cell_style": "center"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function reflect(theta)\n",
" v = [cos(theta), sin(theta)]\n",
" I - 2 * v * v'\n",
"end\n",
"\n",
"Aplot(reflect(0.3))"
]
},
{
"cell_type": "markdown",
"id": "e196bafd",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Singular Value Decomposition\n",
"\n",
"The SVD is $A = U \\Sigma V^T$ where $U$ and $V$ have orthonormal columns and $\\Sigma$ is diagonal with nonnegative entries. It exists for any matrix (non-square, singular, etc.). If we think of orthogonal matrices as reflections/rotations, this says any matrix can be represented as reflect/rotate, diagonally scale, and reflect/rotate again.\n",
"\n",
"Let's try a random symmetric matrix."
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "b6415fe7",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"det(A) = -1.6820338013412035\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = randn(2, 2)\n",
"A += A' # make symmetric\n",
"@show det(A) # Positive means orientation is preserved\n",
"Aplot(A)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "bec2a378",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"norm(U * diagm(S) * V' - A) = 4.002966042486721e-16\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"U, S, V = svd(A)\n",
"@show norm(U * diagm(S) * V' - A) # Should be zero\n",
"Aplot(V') # Rotate/reflect in preparation for scaling"
]
},
{
"cell_type": "markdown",
"id": "53ea375d",
"metadata": {},
"source": [
"* What visual features indicate that this is a symmetric matrix?\n",
"* Is the orthogonal matrix a reflection or rotation?\n",
" * Does this change when the determinant is positive versus negative (rerun the cell above as needed).\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "13e054eb",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Parts of the SVD"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "a0a9cdc4",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Aplot(diagm(S)) # scale along axes"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "39ced49d",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Aplot(U) # rotate/reflect back"
]
},
{
"cell_type": "markdown",
"id": "3086892c",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Putting it together"
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "237dfc09",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Aplot(U * diagm(S) * V')"
]
},
{
"cell_type": "markdown",
"id": "254d918d",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Observations\n",
"\n",
"* The circle always maps to an ellipse\n",
"* The $U$ and $V$ factors may reflect even when $\\det A > 0$"
]
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"celltoolbar": "Slideshow",
"hide_code_all_hidden": false,
"kernelspec": {
"display_name": "Julia 1.7.2",
"language": "julia",
"name": "julia-1.7"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.4"
},
"rise": {
"enable_chalkboard": true
}
},
"nbformat": 4,
"nbformat_minor": 5
}