10 - Basic linear algebra.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic linear algebra in Julia\n",
    "Author: Andreas Noack Jensen (MIT & JuliaComputing) (https://twitter.com/anoackjensen?lang=en)\n",
    "(with edits from Jane Herriman)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First let's define a random matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       " 3  1  3\n",
       " 1  1  2\n",
       " 3  2  2"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = rand(1:4,3,3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a vector of ones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 1.0\n",
       " 1.0\n",
       " 1.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = fill(1.0, (3,)) # = fill(1.0, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that $A$ has type Array{Int64,2} but $x$ has type Array{Float64,1}. Julia defines the aliases Vector{Type}=Array{Type,1} and Matrix{Type}=Array{Type,2}. \n",
    "\n",
    "Many of the basic operations are the same as in other languages\n",
    "#### Multiplication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 7.0\n",
       " 4.0\n",
       " 7.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = A*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transposition\n",
    "As in other languages `A'` is the conjugate transpose, or adjoint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 adjoint(::Matrix{Int64}) with eltype Int64:\n",
       " 3  1  3\n",
       " 1  1  2\n",
       " 3  2  2"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and we can get the transpose with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 transpose(::Matrix{Int64}) with eltype Int64:\n",
       " 3  1  3\n",
       " 1  1  2\n",
       " 3  2  2"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "transpose(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transposed multiplication\n",
    "Julia allows us to write this without *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       " 19  10  17\n",
       " 10   6   9\n",
       " 17   9  17"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A'A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solving linear systems \n",
    "The problem $Ax=b$ for ***square*** $A$ is solved by the \\ function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 1.0\n",
       " 1.0\n",
       " 1.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A\\b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`A\\b` gives us the *least squares solution* if we have an overdetermined linear system (a \"tall\" matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×2 Matrix{Float64}:\n",
       " 0.01259    0.562346\n",
       " 0.571931   0.208287\n",
       " 0.0761351  0.507003"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Atall = rand(3, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       "  2.3722998636787898\n",
       " 12.85556853172185"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Atall\\b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the *minimum norm least squares solution* if we have a rank-deficient least squares problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×2 Matrix{Float64}:\n",
       " 0.682093  0.682093\n",
       " 0.745287  0.745287\n",
       " 0.698837  0.698837"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = rand(3)\n",
    "rankdef = hcat(v, v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 4.190527048390038\n",
       " 4.190527048390039"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rankdef\\b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia also gives us the minimum norm solution when we have an underdetermined solution (a \"short\" matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×3 Matrix{Float64}:\n",
       " 0.692218  0.193236   0.153774\n",
       " 0.127852  0.0359394  0.063034"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bshort = rand(2)\n",
    "Ashort = rand(2, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " -4.139674407706176\n",
       " -0.9797338415813472\n",
       " 23.563074056363842"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ashort\\bshort"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The LinearAlgebra library\n",
    "\n",
    "While much of linear algebra is available in Julia by default (as shown above), there's a standard library named `LinearAlgebra` that brings in many more relevant names and functions. In particular, it provides factorizations and some structured matrix types.  As with all packages, you can bring these additional features into your session with a `using LinearAlgebra`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "#### 10.1 \n",
    "Take the inner product (or \"dot\" product) of a vector `v` with itself and assign it to variable `dot_v`.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Int64}:\n",
       " 1\n",
       " 2\n",
       " 3"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = [1,2,3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "deletable": false,
    "editable": false,
    "hide_input": true,
    "nbgrader": {
     "checksum": "b93dad361f66498eb2460d708f674220",
     "grade": true,
     "grade_id": "cell-913fef9b0d19cd52",
     "locked": true,
     "points": 1,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "@assert dot_v == 14"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 10.2 \n",
    "Take the outer product of a vector v with itself and assign it to variable `outer_v`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "@assert outer_v == [1 2 3\n",
    "                    2 4 6\n",
    "                    3 6 9]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 10.3 \n",
    "Use [LinearAlgebra.cross](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.cross) to compute the cross product of a vector v with itself and assign it to variable `cross_v`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "deletable": false,
    "editable": false,
    "hide_input": true,
    "nbgrader": {
     "checksum": "01642581e27c0ac19752cd90d11ac2ae",
     "grade": true,
     "grade_id": "cell-e6b6970ffe104df5",
     "locked": true,
     "points": 1,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "@assert cross_v == [0, 0, 0]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.7.1",
   "language": "julia",
   "name": "julia-1.7"
  },
  "language": "Julia",
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}