{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 1. Fitting with the backslash operator"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"(process:14972): GLib-GIO-WARNING **: 20:56:35.780: Unexpectedly, UWP app `Evernote.Evernote_10.45.18.0_x64__q4d96b2w5wcc2' (AUMId `Evernote.Evernote_q4d96b2w5wcc2!Evernote') supports 1 extensions but has no verbs\n"
]
}
],
"source": [
"using Plots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider as the original data, points on the line\n",
"\n",
"$$\n",
" y = 0.5 x + 1.\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = 0:0.1:1\n",
"y = [0.5*t + 1 for t in x]\n",
"plot(x, y, label=\"original\", legend=:topleft)\n",
"scatter!(x,y, marker=(:blue), label=\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us now add some random noise to it."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ynoisy = [z + 0.2*(rand() - 0.5) for z in y]\n",
"scatter!(x, ynoisy, marker=(:red), label=\"noisy data\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us fit the noisy data with the backslash operator. The slope and the intercept of the line \n",
"\n",
"$$\n",
" y = c_1 + c_2 x\n",
"$$\n",
"\n",
"that fits the data are the solution of an overdetermined linear system:\n",
"\n",
"$$\n",
" \\left[\n",
" \\begin{array}{cc}\n",
" 1 & x_1 \\\\\n",
" 1 & x_2 \\\\\n",
" \\vdots & \\vdots \\\\\n",
" 1 & x_n\n",
" \\end{array}\n",
" \\right]\n",
" \\left[\n",
" \\begin{array}{c}\n",
" c_1 \\\\ c_2\n",
" \\end{array}\n",
" \\right]\n",
" =\n",
" \\left[\n",
" \\begin{array}{c}\n",
" y_1 \\\\ y_2 \\\\ \\vdots \\\\ y_n\n",
" \\end{array}\n",
" \\right].\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Vector{Float64}:\n",
" 0.9921716122792645\n",
" 0.4945273053175294"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = ones(length(x),2)\n",
"A[:,2] = x\n",
"c = A\\ynoisy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us complete the plot with the line with intercept in ``c[1]`` and slope in ``c[2]``."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fity = c[1] .+ c[2]*x\n",
"plot!(x, fity, label=\"fitted data\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The backslash operator solves linear systems in the least squares sense."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# savefig(\"dataplot.png\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Least Squares"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We solve a linear system with more equations than unknows in the least squares sense, minimizing the sum of the squares of the components of the residual vector $b - Ax$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The residual vector is perpendicular to the column space of the matrix."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us illustrate this property in three dimensions with a fit through (0,1), (1,1), and (1,2)."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We set up the linear system to do the fit and solve it with the backslash operator."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Vector{Float64}:\n",
" 0.8333333333333333\n",
" 0.5000000000000001"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = [1 0; 1 1; 1 2]\n",
"b = [1; 1; 2]\n",
"x = A\\b"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us compute the residual vector."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 0.16666666666666674\n",
" -0.3333333333333335\n",
" 0.16666666666666652"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"r = b - A*x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Unlike when we solved linear systems that had as many equations as unknowns, the residual vector is no longer zero, but we can verify that the residual vector is perpendicular to the column space of ``A`` as follows:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 transpose(::Vector{Float64}) with eltype Float64:\n",
" -2.22045e-16 -4.44089e-16"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"transpose(r)*A"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"By ``transpose(r)`` we turn ``r`` into a row vector and the multiplication ``transpose(r)*A`` makes all inner products of ``r`` with the columns of ``A``. We see that the result is a vector with two numbers close to machine precision."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For a good plot, let us normalize the three vectors first, so they all have the same length."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 0.5773502691896258\n",
" 0.5773502691896258\n",
" 0.5773502691896258"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a1 = [A[1,1]; A[2,1]; A[3,1]]\n",
"v1 = a1/norm(a1)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 0.0\n",
" 0.4472135954999579\n",
" 0.8944271909999159"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a2 = [A[1,2]; A[2,2]; A[3,2]]\n",
"v2 = a2/norm(a2)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 0.4082482904638631\n",
" -0.8164965809277261\n",
" 0.4082482904638625"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v3 = r/norm(r)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X = [0, v1[1], v2[1], v3[1]]\n",
"Y = [0, v1[2], v2[2], v3[2]]\n",
"Z = [0, v1[3], v2[3], v3[3]]\n",
"scatter(X,Y,Z, leg=false, color=\"black\")\n",
"V1x = [0, v1[1]]\n",
"V1y = [0, v1[2]]\n",
"V1z = [0, v1[3]]\n",
"plot!(V1x,V1y,V1z, leg=false, color=\"green\")\n",
"V2x = [0, v2[1]]\n",
"V2y = [0, v2[2]]\n",
"V2z = [0, v2[3]]\n",
"plot!(V2x,V2y,V2z, leg=false, color=\"green\")\n",
"V3x = [0, v3[1]]\n",
"V3y = [0, v3[2]]\n",
"V3z = [0, v3[3]]\n",
"plot!(V3x,V3y,V3z, leg=false, color=\"red\", camera=(100,10), aspect=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the right camera angle, we can see that the red vector is perpendicular to the plane spanned by the green vectors."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# savefig(\"leastsquares.png\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.8.0",
"language": "julia",
"name": "julia-1.8"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}