From 829a5aba92c6a5802f5c0b3a6610e95cd5dc3055 Mon Sep 17 00:00:00 2001 From: Vincent Tembo Date: Fri, 10 Sep 2021 15:56:06 -0400 Subject: [PATCH] a_diag only approximately equal to matrix because of roundoff --- 11 - Factorizations and other fun.ipynb | 1927 +++++++++++------------ 1 file changed, 899 insertions(+), 1028 deletions(-) diff --git a/11 - Factorizations and other fun.ipynb b/11 - Factorizations and other fun.ipynb index 443f17d..02c0dc6 100644 --- a/11 - Factorizations and other fun.ipynb +++ b/11 - Factorizations and other fun.ipynb @@ -1,1032 +1,903 @@ { - "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 + "cells": [ + { + "cell_type": "markdown", + "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" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "Before we get started, let's set up a linear system and use `LinearAlgebra` to bring in the factorizations and special matrix structures." + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "using LinearAlgebra\n", + "A = rand(3, 3)\n", + "x = fill(1, (3,))\n", + "b = A * x" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 1, + "data": { + "text/plain": "3-element Vector{Float64}:\n 2.205296596065241\n 1.8551229792583508\n 0.8533886927959882" + }, + "metadata": {} + } + ], + "execution_count": 1, + "metadata": {} + }, + { + "cell_type": "markdown", + "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." + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "Alu = lu(A)" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 2, + "data": { + "text/plain": "LU{Float64, Matrix{Float64}}\nL factor:\n3×3 Matrix{Float64}:\n 1.0 0.0 0.0\n 0.513615 1.0 0.0\n 0.295738 0.271404 1.0\nU factor:\n3×3 Matrix{Float64}:\n 0.83195 0.576718 0.446455\n 0.0 0.687695 0.564783\n 0.0 0.0 -0.0351695" + }, + "metadata": {} + } + ], + "execution_count": 2, + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "typeof(Alu)" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 3, + "data": { + "text/plain": "LU{Float64, Matrix{Float64}}" + }, + "metadata": {} + } + ], + "execution_count": 3, + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "The different parts of the factorization can be extracted by accessing their special properties" + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "Alu.P" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 4, + "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" + }, + "metadata": {} + } + ], + "execution_count": 4, + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "Alu.L" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 5, + "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" + }, + "metadata": {} + } + ], + "execution_count": 5, + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "Alu.U" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 6, + "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" + }, + "metadata": {} + } + ], + "execution_count": 6, + "metadata": {} + }, + { + "cell_type": "markdown", + "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." + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "A\\b" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 7, + "data": { + "text/plain": "3-element Vector{Float64}:\n 0.9999999999999999\n 1.0000000000000002\n 1.0" + }, + "metadata": {} + } + ], + "execution_count": 7, + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "Alu\\b" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 8, + "data": { + "text/plain": "3-element Vector{Float64}:\n 0.9999999999999999\n 1.0000000000000002\n 1.0" + }, + "metadata": {} + } + ], + "execution_count": 8, + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "Similarly, we can calculate the determinant of `A` using either `A` or the factorization object" + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "det(A) ≈ det(Alu)" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 9, + "data": { + "text/plain": "true" + }, + "metadata": {} + } + ], + "execution_count": 9, + "metadata": {} + }, + { + "cell_type": "markdown", + "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`. " + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "Aqr = qr(A)" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 10, + "data": { + "text/plain": "LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}}\nQ factor:\n3×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\nR factor:\n3×3 Matrix{Float64}:\n -0.96709 -1.02174 -0.798571\n 0.0 0.619937 0.503618\n 0.0 0.0 -0.0335617" + }, + "metadata": {} + } + ], + "execution_count": 10, + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "Similarly to the LU factorization, the matrices `Q` and `R` can be extracted from the QR factorization object via" + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "Aqr.Q" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 11, + "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" + }, + "metadata": {} + } + ], + "execution_count": 11, + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "Aqr.R" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 12, + "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" + }, + "metadata": {} + } + ], + "execution_count": 12, + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "#### Eigendecompositions" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "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" + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "Asym = A + A'\n", + "AsymEig = eigen(Asym)" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 13, + "data": { + "text/plain": "Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}\nvalues:\n3-element Vector{Float64}:\n -0.871021215213241\n -0.02678184303514708\n 3.406139917482376\nvectors:\n3×3 Matrix{Float64}:\n 0.768088 -0.0906084 -0.633901\n -0.596429 -0.461515 -0.656716\n -0.233051 0.882493 -0.408526" + }, + "metadata": {} + } + ], + "execution_count": 13, + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "The values and the vectors can be extracted from the Eigen type by special indexing" + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "AsymEig.values" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 14, + "data": { + "text/plain": "3-element Vector{Float64}:\n -0.871021215213241\n -0.02678184303514708\n 3.406139917482376" + }, + "metadata": {} + } + ], + "execution_count": 14, + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "AsymEig.vectors" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 15, + "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" + }, + "metadata": {} + } + ], + "execution_count": 15, + "metadata": {} + }, + { + "cell_type": "markdown", + "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}$." + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "inv(AsymEig)*Asym" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 16, + "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" + }, + "metadata": {} + } + ], + "execution_count": 16, + "metadata": {} + }, + { + "cell_type": "markdown", + "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" + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "n = 1000\n", + "A = randn(n,n);" + ], + "outputs": [], + "execution_count": 17, + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "Julia can often infer special matrix structure" + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "Asym = A + A'\n", + "issymmetric(Asym)" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 18, + "data": { + "text/plain": "true" + }, + "metadata": {} + } + ], + "execution_count": 18, + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "but sometimes floating point error might get in the way." + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "Asym_noisy = copy(Asym)\n", + "Asym_noisy[1,2] += 5eps()" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 19, + "data": { + "text/plain": "2.9631367153816814" + }, + "metadata": {} + } + ], + "execution_count": 19, + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "issymmetric(Asym_noisy)" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 20, + "data": { + "text/plain": "false" + }, + "metadata": {} + } + ], + "execution_count": 20, + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "Luckily we can declare structure explicitly with, for example, `Diagonal`, `Triangular`, `Symmetric`, `Hermitian`, `Tridiagonal` and `SymTridiagonal`." + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "Asym_explicit = Symmetric(Asym_noisy);" + ], + "outputs": [], + "execution_count": 21, + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "Let's compare how long it takes Julia to compute the eigenvalues of `Asym`, `Asym_noisy`, and `Asym_explicit`" + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "@time eigvals(Asym);" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + " 0.258496 seconds (134.03 k allocations: 15.456 MiB, 6.10% gc time, 38.03% compilation time)\n" + ] + } + ], + "execution_count": 22, + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "@time eigvals(Asym_noisy);" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + " 1.050047 seconds (13 allocations: 7.920 MiB)\n" + ] + } + ], + "execution_count": 23, + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "@time eigvals(Asym_explicit);" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + " 0.149938 seconds (5.88 k allocations: 8.358 MiB, 2.11% compilation time)\n" + ] + } + ], + "execution_count": 24, + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "In this example, using `Symmetric()` on `Asym_noisy` made our calculations about `5x` more efficient :)" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "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." + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "n = 1_000_000;\n", + "A = SymTridiagonal(randn(n), randn(n-1));\n", + "@time eigmax(A)" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + " 0.669902 seconds (823.18 k allocations: 229.291 MiB, 12.15% gc time, 26.00% compilation time)\n" + ] + }, + { + "output_type": "execute_result", + "execution_count": 25, + "data": { + "text/plain": "6.7249639294089665" + }, + "metadata": {} + } + ], + "execution_count": 25, + "metadata": {} + }, + { + "cell_type": "markdown", + "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." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "#### Rational numbers\n", + "Julia has rational numbers built in. To construct a rational number, use double forward slashes:" + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "1//2" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 26, + "data": { + "text/plain": "1//2" + }, + "metadata": {} + } + ], + "execution_count": 26, + "metadata": {} + }, + { + "cell_type": "markdown", + "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." + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "Arational = Matrix{Rational{BigInt}}(rand(1:10, 3, 3))/10" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 27, + "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" + }, + "metadata": {} + } + ], + "execution_count": 27, + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "x = fill(1, 3)\n", + "b = Arational*x" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 28, + "data": { + "text/plain": "3-element Vector{Rational{BigInt}}:\n 11//10\n 1//2\n 7//10" + }, + "metadata": {} + } + ], + "execution_count": 28, + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "Arational\\b" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 29, + "data": { + "text/plain": "3-element Vector{Rational{BigInt}}:\n 1//1\n 1//1\n 1//1" + }, + "metadata": {} + } + ], + "execution_count": 29, + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "lu(Arational)" + ], + "outputs": [ + { + "output_type": "execute_result", + "execution_count": 30, + "data": { + "text/plain": "LU{Rational{BigInt}, Matrix{Rational{BigInt}}}\nL factor:\n3×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\nU factor:\n3×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" + }, + "metadata": {} + } + ], + "execution_count": 30, + "metadata": {} + }, + { + "cell_type": "markdown", + "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`" + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "using LinearAlgebra" + ], + "outputs": [], + "execution_count": null, + "metadata": {} + }, + { + "cell_type": "code", + "source": [], + "outputs": [], + "execution_count": null, + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "@assert A_eigv == [-128.49322764802145, -55.887784553056875, 42.7521672793189, 87.16111477514521, 542.4677301466143]" + ], + "outputs": [], + "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 + } + } + }, + { + "cell_type": "markdown", + "source": [ + "#### 11.2 \n", + "Create a `Diagonal` matrix from the eigenvalues of `A`." + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [], + "outputs": [], + "execution_count": null, + "metadata": {} + }, + { + "cell_type": "code", + "source": [ + "@assert isapprox(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], atol=1e-3)" + ], + "outputs": [], + "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 + } + } + }, + { + "cell_type": "markdown", + "source": [ + "#### 11.3 \n", + "Create a `LowerTriangular` matrix from `A` and store it in `A_lowertri`" + ], + "metadata": {} + }, + { + "cell_type": "code", + "source": [], + "outputs": [], + "execution_count": null, + "metadata": {} + }, + { + "cell_type": "code", + "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]" + ], + "outputs": [], + "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 + } + } + }, + { + "cell_type": "markdown", + "source": [ + "### Please let us know how we're doing!\n", + "https://tinyurl.com/introJuliaFeedback" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "Please click on `Validate` on the top, once you are done with the exercises." + ], + "metadata": {} } - }, - "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 + ], + "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.0" + }, + "nteract": { + "version": "0.28.0" } - }, - "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]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Please let us know how we're doing!\n", - "https://tinyurl.com/introJuliaFeedback" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Please click on `Validate` on the top, once you are done with the exercises." - ] - } - ], - "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.0" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file