11 - Factorizations and other fun.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Factorizations and other fun\n",
    "Based on work by Andreas Noack\n",
    "\n",
    "## Outline\n",
    " - Factorizations\n",
    " - Special matrix structures\n",
    " - Generic linear algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we get started, let's set up a linear system and use `LinearAlgebra` to bring in the factorizations and special matrix structures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 2.205296596065241\n",
       " 1.8551229792583508\n",
       " 0.8533886927959882"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "A = rand(3, 3)\n",
    "x = fill(1, (3,))\n",
    "b = A * x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Factorizations\n",
    "\n",
    "#### LU factorizations\n",
    "In Julia we can perform an LU factorization\n",
    "```julia\n",
    "PA = LU\n",
    "``` \n",
    "where `P` is a permutation matrix, `L` is lower triangular unit diagonal and `U` is upper triangular, using `lufact`.\n",
    "\n",
    "Julia allows computing the LU factorization and defines a composite factorization type for storing it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LU{Float64, Matrix{Float64}}\n",
       "L factor:\n",
       "3×3 Matrix{Float64}:\n",
       " 1.0       0.0       0.0\n",
       " 0.513615  1.0       0.0\n",
       " 0.295738  0.271404  1.0\n",
       "U factor:\n",
       "3×3 Matrix{Float64}:\n",
       " 0.83195  0.576718   0.446455\n",
       " 0.0      0.687695   0.564783\n",
       " 0.0      0.0       -0.0351695"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Alu = lu(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LU{Float64, Matrix{Float64}}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(Alu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The different parts of the factorization can be extracted by accessing their special properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 0.0  1.0  0.0\n",
       " 1.0  0.0  0.0\n",
       " 0.0  0.0  1.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Alu.P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 1.0       0.0       0.0\n",
       " 0.513615  1.0       0.0\n",
       " 0.295738  0.271404  1.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Alu.L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 0.83195  0.576718   0.446455\n",
       " 0.0      0.687695   0.564783\n",
       " 0.0      0.0       -0.0351695"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Alu.U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia can dispatch methods on factorization objects.\n",
    "\n",
    "For example, we can solve the linear system using either the original matrix or the factorization object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 0.9999999999999999\n",
       " 1.0000000000000002\n",
       " 1.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A\\b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 0.9999999999999999\n",
       " 1.0000000000000002\n",
       " 1.0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Alu\\b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, we can calculate the determinant of `A` using either `A` or the factorization object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "det(A) ≈ det(Alu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### QR factorizations\n",
    "\n",
    "In Julia we can perform a QR factorization\n",
    "```\n",
    "A=QR\n",
    "``` \n",
    "\n",
    "where `Q` is unitary/orthogonal and `R` is upper triangular, using `qrfact`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}}\n",
       "Q factor:\n",
       "3×3 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}}:\n",
       " -0.441843   0.858892  -0.258997\n",
       " -0.860261  -0.487536  -0.149194\n",
       " -0.254412   0.156884   0.954286\n",
       "R factor:\n",
       "3×3 Matrix{Float64}:\n",
       " -0.96709  -1.02174   -0.798571\n",
       "  0.0       0.619937   0.503618\n",
       "  0.0       0.0       -0.0335617"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Aqr = qr(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly to the LU factorization, the matrices `Q` and `R` can be extracted from the QR factorization object via"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}}:\n",
       " -0.441843   0.858892  -0.258997\n",
       " -0.860261  -0.487536  -0.149194\n",
       " -0.254412   0.156884   0.954286"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Aqr.Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " -0.96709  -1.02174   -0.798571\n",
       "  0.0       0.619937   0.503618\n",
       "  0.0       0.0       -0.0335617"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Aqr.R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Eigendecompositions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results from eigendecompositions, singular value decompositions, Hessenberg factorizations, and Schur decompositions are all stored in `Factorization` types.\n",
    "\n",
    "The eigendecomposition can be computed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}\n",
       "values:\n",
       "3-element Vector{Float64}:\n",
       " -0.871021215213241\n",
       " -0.02678184303514708\n",
       "  3.406139917482376\n",
       "vectors:\n",
       "3×3 Matrix{Float64}:\n",
       "  0.768088  -0.0906084  -0.633901\n",
       " -0.596429  -0.461515   -0.656716\n",
       " -0.233051   0.882493   -0.408526"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Asym = A + A'\n",
    "AsymEig = eigen(Asym)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The values and the vectors can be extracted from the Eigen type by special indexing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " -0.871021215213241\n",
       " -0.02678184303514708\n",
       "  3.406139917482376"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AsymEig.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       "  0.768088  -0.0906084  -0.633901\n",
       " -0.596429  -0.461515   -0.656716\n",
       " -0.233051   0.882493   -0.408526"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AsymEig.vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again, when the factorization is stored in a type, we can dispatch on it and write specialized methods that exploit the properties of the factorization, e.g. that $A^{-1}=(V\\Lambda V^{-1})^{-1}=V\\Lambda^{-1}V^{-1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       "  1.0           5.32907e-15  4.88498e-15\n",
       "  2.84217e-14   1.0          2.22045e-14\n",
       " -5.68434e-14  -6.03961e-14  1.0"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(AsymEig)*Asym"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Special matrix structures\n",
    "Matrix structure is very important in linear algebra. To see *how* important it is, let's work with a larger linear system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 1000\n",
    "A = randn(n,n);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia can often infer special matrix structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Asym = A + A'\n",
    "issymmetric(Asym)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "but sometimes floating point error might get in the way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.9631367153816814"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Asym_noisy = copy(Asym)\n",
    "Asym_noisy[1,2] += 5eps()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "false"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "issymmetric(Asym_noisy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Luckily we can declare structure explicitly with, for example, `Diagonal`, `Triangular`, `Symmetric`, `Hermitian`, `Tridiagonal` and `SymTridiagonal`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "Asym_explicit = Symmetric(Asym_noisy);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare how long it takes Julia to compute the eigenvalues of `Asym`, `Asym_noisy`, and `Asym_explicit`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.258496 seconds (134.03 k allocations: 15.456 MiB, 6.10% gc time, 38.03% compilation time)\n"
     ]
    }
   ],
   "source": [
    "@time eigvals(Asym);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1.050047 seconds (13 allocations: 7.920 MiB)\n"
     ]
    }
   ],
   "source": [
    "@time eigvals(Asym_noisy);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.149938 seconds (5.88 k allocations: 8.358 MiB, 2.11% compilation time)\n"
     ]
    }
   ],
   "source": [
    "@time eigvals(Asym_explicit);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, using `Symmetric()` on `Asym_noisy` made our calculations about `5x` more efficient :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A big problem\n",
    "Using the `Tridiagonal` and `SymTridiagonal` types to store tridiagonal matrices makes it possible to work with potentially very large tridiagonal problems. The following problem would not be possible to solve on a laptop if the matrix had to be stored as a (dense) `Matrix` type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.669902 seconds (823.18 k allocations: 229.291 MiB, 12.15% gc time, 26.00% compilation time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "6.7249639294089665"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 1_000_000;\n",
    "A = SymTridiagonal(randn(n), randn(n-1));\n",
    "@time eigmax(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generic linear algebra\n",
    "The usual way of adding support for numerical linear algebra is by wrapping BLAS and LAPACK subroutines. For matrices with elements of `Float32`, `Float64`, `Complex{Float32}` or `Complex{Float64}` this is also what Julia does.\n",
    "\n",
    "However, Julia also supports generic linear algebra, allowing you to, for example, work with matrices and vectors of rational numbers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Rational numbers\n",
    "Julia has rational numbers built in. To construct a rational number, use double forward slashes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1//2"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1//2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example: Rational linear system of equations\n",
    "The following example shows how linear system of equations with rational elements can be solved without promoting to floating point element types. Overflow can easily become a problem when working with rational numbers so we use `BigInt`s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Rational{BigInt}}:\n",
       " 1//10  3//5   2//5\n",
       " 1//5   1//10  1//5\n",
       " 1//5   1//5   3//10"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Arational = Matrix{Rational{BigInt}}(rand(1:10, 3, 3))/10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Rational{BigInt}}:\n",
       " 11//10\n",
       "  1//2\n",
       "  7//10"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = fill(1, 3)\n",
    "b = Arational*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Rational{BigInt}}:\n",
       " 1//1\n",
       " 1//1\n",
       " 1//1"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Arational\\b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LU{Rational{BigInt}, Matrix{Rational{BigInt}}}\n",
       "L factor:\n",
       "3×3 Matrix{Rational{BigInt}}:\n",
       " 1//1  0//1   0//1\n",
       " 1//2  1//1   0//1\n",
       " 1//1  2//11  1//1\n",
       "U factor:\n",
       "3×3 Matrix{Rational{BigInt}}:\n",
       " 1//5   1//10  1//5\n",
       " 0//1  11//20  3//10\n",
       " 0//1   0//1   1//22"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lu(Arational)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "#### 11.1\n",
    "What are the eigenvalues of matrix A?\n",
    "\n",
    "```\n",
    "A =\n",
    "[\n",
    " 140   97   74  168  131\n",
    "  97  106   89  131   36\n",
    "  74   89  152  144   71\n",
    " 168  131  144   54  142\n",
    " 131   36   71  142   36\n",
    "]\n",
    "```\n",
    "and assign it a variable `A_eigv`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "hide_input": true,
    "nbgrader": {
     "checksum": "f9f16fdef201ed372323a291f1dd1346",
     "grade": true,
     "grade_id": "cell-4d5f60c8a814c789",
     "locked": true,
     "points": 0,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "@assert A_eigv ==  [-128.49322764802145, -55.887784553056875, 42.7521672793189, 87.16111477514521, 542.4677301466143]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 11.2 \n",
    "Create a `Diagonal` matrix from the eigenvalues of `A`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "hide_input": true,
    "nbgrader": {
     "checksum": "3ca676f6282c1a7c214ab2cb9f9b322d",
     "grade": true,
     "grade_id": "cell-3b000a3710c9c263",
     "locked": true,
     "points": 1,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "@assert A_diag ==  [-128.493    0.0      0.0      0.0       0.0;\n",
    "    0.0    -55.8878   0.0      0.0       0.0;\n",
    "    0.0      0.0     42.7522   0.0       0.0;\n",
    "    0.0      0.0      0.0     87.1611    0.0;\n",
    "    0.0 0.0      0.0      0.0     542.468]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 11.3 \n",
    "Create a `LowerTriangular` matrix from `A` and store it in `A_lowertri`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "hide_input": true,
    "nbgrader": {
     "checksum": "b3b1a272343a05082f378a5e1aa3426d",
     "grade": true,
     "grade_id": "cell-b76cee2b4a8777da",
     "locked": true,
     "points": 0,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "@assert A_lowertri ==  [140    0    0    0   0;\n",
    "  97  106    0    0   0;\n",
    "  74   89  152    0   0;\n",
    " 168  131  144   54   0;\n",
    " 131   36   71  142  36]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.0",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language": "Julia",
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}