{
"cells": [
{
"cell_type": "markdown",
"id": "4159be0e",
"metadata": {
"hideCode": false,
"hidePrompt": false,
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 2022-02-07 Linear Algebra\n",
"\n",
"\n",
"## Last time\n",
"\n",
"* Forward and backward stability\n",
"* Beyond IEEE double precision\n",
"* Discuss [portfolios](https://classroom.github.com/a/strHwkaV)\n",
"\n",
"## Today\n",
"\n",
"* Algebra of linear transformations\n",
"* Polynomial evaluation and fitting\n",
"* Orthogonality"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "eb781a7a",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"using Plots\n",
"default(linewidth=4, legendfontsize=12)"
]
},
{
"cell_type": "markdown",
"id": "1f18c04a",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Matrices as linear transformations\n",
"\n",
"Linear algebra is the study of linear transformations on vectors, which represent points in a finite dimensional space. The matrix-vector product $y = A x$ is a linear combination of the columns of $A$. The familiar definition,\n",
"\n",
"$$ y_i = \\sum_j A_{i,j} x_j $$\n",
"\n",
"can also be viewed as\n",
"\n",
"$$ y = \\Bigg[ A_{:,1} \\Bigg| A_{:,2} \\Bigg| \\dotsm \\Bigg] \\begin{bmatrix} x_1 \\\\ x_2 \\\\ \\vdots \\end{bmatrix}\n",
"= \\Bigg[ A_{:,1} \\Bigg] x_1 + \\Bigg[ A_{:,2} \\Bigg] x_2 + \\dotsb . $$"
]
},
{
"cell_type": "markdown",
"id": "47d22601",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Math and Julia Notation\n",
"\n",
"The notation $A_{i,j}$ corresponds to the Julia syntax `A[i,j]` and the colon `:` means the entire range (row or column). So $A_{:,j}$ is the $j$th column and $A_{i,:}$ is the $i$th row. The corresponding Julia syntax is `A[:,j]` and `A[i,:]`.\n",
"\n",
"Julia has syntax for row vectors, column vectors, and arrays."
]
},
{
"cell_type": "code",
"execution_count": 44,
"id": "ff3361f2",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"1×3 Matrix{Float64}:\n",
" 1.0 2.0 3.0"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[1. 2 3]"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "d35b1571",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Int64}:\n",
" 1\n",
" 2\n",
" 3"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[1, 2, 3]"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "a32da848",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"3×2 Matrix{Int64}:\n",
" 1 0\n",
" 0 2\n",
" 10 3"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[1 0; 0 2; 10 3]"
]
},
{
"cell_type": "code",
"execution_count": 46,
"id": "bc208ae0",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"1×3 adjoint(::Vector{Int64}) with eltype Int64:\n",
" 1 2 3"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[1; 2; 3]' # transpose"
]
},
{
"cell_type": "markdown",
"id": "5ac73ba3",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Implementing multiplication by row"
]
},
{
"cell_type": "code",
"execution_count": 51,
"id": "32665740",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"3×4 reshape(::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, 3, 4) with eltype Float64:\n",
" 1.0 4.0 7.0 10.0\n",
" 2.0 5.0 8.0 11.0\n",
" 3.0 6.0 9.0 12.0"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function matmult1(A, x)\n",
" m, n = size(A)\n",
" y = zeros(m)\n",
" for i in 1:m\n",
" for j in 1:n\n",
" y[i] += A[i,j] * x[j]\n",
" end\n",
" end\n",
" y\n",
"end\n",
"\n",
"A = reshape(1.:12, 3, 4) # 3x4 matrix\n",
"#x = [10., 0, 0, 0]\n",
"#matmult1(A, x)"
]
},
{
"cell_type": "code",
"execution_count": 57,
"id": "1e791ea6",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"20.0"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Dot product\n",
"A[2,:]' * x"
]
},
{
"cell_type": "code",
"execution_count": 58,
"id": "a0cf9482",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 10.0\n",
" 20.0\n",
" 30.0"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function matmult2(A, x)\n",
" m, n = size(A)\n",
" y = zeros(m)\n",
" for i in 1:m\n",
" y[i] = A[i,:]' * x\n",
" end\n",
" y\n",
"end\n",
"\n",
"matmult2(A, x)"
]
},
{
"cell_type": "markdown",
"id": "4f3af7b3",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Implementing multiplication by column\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "12c1dea5",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 10.0\n",
" 20.0\n",
" 30.0"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function matmult3(A, x)\n",
" m, n = size(A)\n",
" y = zeros(m)\n",
" for j in 1:n\n",
" y += A[:, j] * x[j]\n",
" end\n",
" y\n",
"end\n",
"\n",
"matmult3(A, x)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "b573daf7",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 10.0\n",
" 20.0\n",
" 30.0"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A * x # We'll use this version "
]
},
{
"cell_type": "markdown",
"id": "ffa46c2c",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Polynomial evaluation is (continuous) linear algebra\n",
"We can evaluate polynomials using matrix-vector multiplication.\n",
"For example,\n",
"$$ 5x^3 - 3x = \\Bigg[ 1 \\Bigg|\\, x \\Bigg|\\, x^2 \\,\\Bigg|\\, x^3 \\Bigg] \\begin{bmatrix}0 \\\\ -3 \\\\ 0 \\\\ 5 \\end{bmatrix} . $$"
]
},
{
"cell_type": "code",
"execution_count": 59,
"id": "a7e1f3dc",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"f = Polynomial(1 - x + 3*x^2 + 9*x^3)\n",
"P(p + q) = Polynomial(1 - x + 3*x^2 + 9*x^3)\n"
]
},
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 1.0\n",
" 12.0\n",
" 83.0"
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Polynomials\n",
"P(x) = Polynomial(x)\n",
"\n",
"p = [0, -3, 0, 5]\n",
"q = [1, 2, 3, 4]\n",
"f = P(p) + P(q)\n",
"@show f\n",
"@show P(p+q)\n",
"x = [0., 1, 2]\n",
"f.(x)"
]
},
{
"cell_type": "code",
"execution_count": 65,
"id": "3c90d1af",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 65,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(f, legend=:bottomright)"
]
},
{
"cell_type": "markdown",
"id": "68418fd7",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Polynomial evaluation is (discrete) linear algebra"
]
},
{
"cell_type": "code",
"execution_count": 66,
"id": "6fdc3793",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"3×4 Matrix{Float64}:\n",
" 1.0 0.0 0.0 0.0\n",
" 1.0 1.0 1.0 1.0\n",
" 1.0 2.0 4.0 8.0"
]
},
"execution_count": 66,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"V = [one.(x) x x.^2 x.^3]"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "a7bd6f28",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 1.0\n",
" 12.0\n",
" 83.0"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"V * p + V * q"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "c18070be",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 1.0\n",
" 12.0\n",
" 83.0"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"V * (p + q)"
]
},
{
"cell_type": "markdown",
"id": "bcac861a",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Vandermonde matrices\n",
"\n",
"A Vandermonde matrix is one whose columns are functions evaluated at discrete points.\n",
"\n",
"$$V(x) = \\begin{bmatrix} 1 \\Bigg| x \\Bigg| x^2 \\Bigg| x^3 \\Bigg| \\dotsb \\end{bmatrix}$$"
]
},
{
"cell_type": "code",
"execution_count": 67,
"id": "55291bbf",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"vander (generic function with 2 methods)"
]
},
"execution_count": 67,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"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"
]
},
{
"cell_type": "code",
"execution_count": 68,
"id": "6151d57b",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = LinRange(-1, 1, 50)\n",
"V = vander(x, 4)\n",
"plot(x, V)"
]
},
{
"cell_type": "markdown",
"id": "ce88943a",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Fitting is linear algebra\n",
"\n",
"$$ \\underbrace{\\begin{bmatrix} 1 \\Bigg| x \\Bigg| x^2 \\Bigg| x^3 \\Bigg| \\dotsb \\end{bmatrix}}_{V(x)} \\Big[ p \\Big] = \\Bigg[ y \\Bigg]$$"
]
},
{
"cell_type": "code",
"execution_count": 69,
"id": "451865ac",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 69,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x1 = [-.9, 0.1, .5, .8]\n",
"y1 = [1, 2.4, -.2, 1.3]\n",
"scatter(x1, y1, markersize=8)"
]
},
{
"cell_type": "code",
"execution_count": 82,
"id": "62a84cfe",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"size(V) = (4, 4)\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 82,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"V = vander(x1, 4)\n",
"@show size(V)\n",
"p = V \\ y1 # write y1 in the polynomial basis\n",
"scatter(x1, y1, markersize=8, xlims=(-1, 1))\n",
"#plot!(P(p), label=\"P(p)\")\n",
"plot!(x, vander(x, 4) * p, label=\"\\$ V(x) p\\$\", linestyle=:dash)"
]
},
{
"cell_type": "markdown",
"id": "17d2ea90",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Some common terminology\n",
"\n",
"* The **range** of $A$ is the space spanned by its columns. This definition coincides with the range of a function $f(x)$ when $f(x) = A x$.\n",
"* The (right) **nullspace** of $A$ is the space of vectors $x$ such that $A x = 0$.\n",
"* The **rank** of $A$ is the dimension of its range.\n",
"* A matrix has **full rank** if the nullspace of either $A$ or $A^T$ is empty (only the 0 vector). Equivalently, if all the columns of $A$ (or $A^T$) are linearly independent.\n",
"* A **nonsingular** (or **invertible**) matrix is a square matrix of full rank. We call the inverse $A^{-1}$ and it satisfies $A^{-1} A = A A^{-1} = I$."
]
},
{
"cell_type": "markdown",
"id": "e5269948",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"$\\DeclareMathOperator{\\rank}{rank} \\DeclareMathOperator{\\null}{null} $\n",
"If $A \\in \\mathbb{R}^{m\\times m}$, which of these doesn't belong?\n",
"1. $A$ has an inverse $A^{-1}$\n",
"2. $\\rank (A) = m$\n",
"3. $\\null(A) = \\{0\\}$\n",
"4. $A A^T = A^T A$\n",
"5. $\\det(A) \\ne 0$\n",
"6. $A x = 0$ implies that $x = 0$"
]
},
{
"cell_type": "code",
"execution_count": 86,
"id": "ecf60702",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"-0.08969850225286818"
]
},
"execution_count": 86,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = rand(4,4)\n",
"#A' * A - A * A'\n",
"det(A)"
]
},
{
"cell_type": "markdown",
"id": "efa97478",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# What is an inverse?\n",
"\n",
"When we write $x = A^{-1} y$, we mean that $x$ is the unique vector such that $A x = y$.\n",
"(It is rare that we explicitly compute a matrix $A^{-1}$, though [it's not as \"bad\"](https://arxiv.org/abs/1201.6035) as people may have told you.)\n",
"A vector $y$ is equivalent to $\\sum_i e_i y_i$ where $e_i$ are columns of the identity.\n",
"Meanwhile, $x = A^{-1} y$ means that we are expressing that same vector $y$ in the basis of the columns of $A$, i.e., $\\sum_i A_{:,i} x_i$.\n"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "f44d4a40",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Matrix{Float64}:\n",
" 0.529731 0.00344266 0.0249711 0.446597\n",
" 0.65379 0.874346 0.629848 0.296424\n",
" 0.308001 0.687619 0.654854 0.627373\n",
" 0.252958 0.033781 0.86887 0.405012"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using LinearAlgebra\n",
"A = rand(4, 4)"
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "6d5610a3",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Matrix{Float64}:\n",
" 1.0 1.69813e-16 6.45103e-17 -1.4185e-18\n",
" -0.0 1.0 7.87401e-17 -2.28895e-17\n",
" 0.0 0.0 1.0 3.32473e-17\n",
" 0.0 0.0 5.25774e-17 1.0"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A \\ A"
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "33f6511a",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Matrix{Float64}:\n",
" 1.0 2.80573e-16 4.18345e-17 1.09056e-17\n",
" 6.36277e-17 1.0 9.18807e-17 -1.04138e-16\n",
" 1.07459e-16 6.39634e-17 1.0 4.11885e-17\n",
" 4.48719e-17 7.42799e-17 7.86144e-17 1.0"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"inv(A) * A"
]
},
{
"cell_type": "markdown",
"id": "7a8eb495",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Inner products and orthogonality\n",
"\n",
"The **inner product**\n",
"$$ x^T y = \\sum_i x_i y_i $$\n",
"of vectors (or columns of a matrix) tell us about their magnitude and about the angle.\n",
"The **norm** is induced by the inner product,\n",
"$$ \\lVert x \\rVert = \\sqrt{x^T x} $$\n",
"and the angle $\\theta$ is defined by\n",
"$$ \\cos \\theta = \\frac{x^T y}{\\lVert x \\rVert \\, \\lVert y \\rVert} . $$\n",
"Inner products are **bilinear**, which means that they satisfy some convenient algebraic properties\n",
"$$ \\begin{split}\n",
"(x + y)^T z &= x^T z + y^T z \\\\\n",
"x^T (y + z) &= x^T y + x^T z \\\\\n",
"(\\alpha x)^T (\\beta y) &= \\alpha \\beta x^T y \\\\\n",
"\\end{split} . $$"
]
},
{
"cell_type": "markdown",
"id": "343e1be5",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Examples with inner products"
]
},
{
"cell_type": "code",
"execution_count": 90,
"id": "8d795390",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x' * y = 1\n",
"y' * x = 1\n"
]
}
],
"source": [
"x = [0, 1]\n",
"y = [1, 1]\n",
"@show x' * y\n",
"@show y' * x;"
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "58505af3",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"cos_θ = 0.49999999999999994\n",
"cos(ϕ - pi / 2) = 0.4999999999999999\n"
]
}
],
"source": [
"ϕ = pi/6\n",
"y = [cos(ϕ), sin(ϕ)]\n",
"cos_θ = x'*y / (norm(x) * norm(y))\n",
"@show cos_θ\n",
"@show cos(ϕ-pi/2);"
]
},
{
"cell_type": "markdown",
"id": "239ef511",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Polynomials can be orthogonal too!\n"
]
},
{
"cell_type": "code",
"execution_count": 39,
"id": "95341495",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = LinRange(-1, 1, 50)\n",
"A = vander(x, 4)\n",
"M = A * [.5 0 0 0; # 0.5\n",
" 0 1 0 0; # x\n",
" 0 0 1 0]' # x^2\n",
"scatter(x, M)\n",
"plot!(x, 0*x, label=:none, color=:black)"
]
},
{
"cell_type": "markdown",
"id": "fe22f024",
"metadata": {
"cell_style": "split"
},
"source": [
"* Which inner product will be zero?\n",
"\n",
" * Which functions are even and odd?"
]
},
{
"cell_type": "markdown",
"id": "ffae59ed",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Polynomial inner products"
]
},
{
"cell_type": "code",
"execution_count": 40,
"id": "3c8c7db8",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"-2.220446049250313e-16"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M[:,1]' * M[:,2]"
]
},
{
"cell_type": "code",
"execution_count": 42,
"id": "9ee5e1f7",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"8.673469387755102"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M[:,1]' * M[:,3]"
]
},
{
"cell_type": "code",
"execution_count": 41,
"id": "68ad1630",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"-4.440892098500626e-16"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M[:,2]' * M[:,3]"
]
}
],
"metadata": {
"@webio": {
"lastCommId": "cd1892c1e2064686902fbdf399e6e3e3",
"lastKernelId": "6eae1736-6531-4224-a6c2-38e7ff163d19"
},
"celltoolbar": "Slideshow",
"hide_code_all_hidden": false,
"kernelspec": {
"display_name": "Julia 1.7.1",
"language": "julia",
"name": "julia-1.7"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.7.1"
},
"rise": {
"enable_chalkboard": true
}
},
"nbformat": 4,
"nbformat_minor": 5
}