{ "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 Array{Float64,1}:\n", " 0.9343211185345781\n", " 1.5152135778152167\n", " 2.018287744747687" ] }, "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,Array{Float64,2}}\n", "L factor:\n", "3×3 Array{Float64,2}:\n", " 1.0 0.0 0.0\n", " 0.878194 1.0 0.0\n", " 0.0938958 0.70389 1.0\n", "U factor:\n", "3×3 Array{Float64,2}:\n", " 0.963615 0.168704 0.382895\n", " 0.0 0.644811 0.0428252\n", " 0.0 0.0 0.308029" ] }, "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,Array{Float64,2}}" ] }, "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 Array{Float64,2}:\n", " 0.0 1.0 0.0\n", " 0.0 0.0 1.0\n", " 1.0 0.0 0.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 Array{Float64,2}:\n", " 1.0 0.0 0.0\n", " 0.878194 1.0 0.0\n", " 0.0938958 0.70389 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 Array{Float64,2}:\n", " 0.963615 0.168704 0.382895\n", " 0.0 0.644811 0.0428252\n", " 0.0 0.0 0.308029" ] }, "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 Array{Float64,1}:\n", " 0.9999999999999999\n", " 1.0000000000000004\n", " 0.9999999999999997" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A\\b" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 0.9999999999999999\n", " 1.0000000000000004\n", " 0.9999999999999997" ] }, "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,Array{Float64,2}}\n", "Q factor:\n", "3×3 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:\n", " -0.0703771 0.655875 -0.751581\n", " -0.749523 -0.531943 -0.394021\n", " -0.658227 0.535597 0.52903\n", "R factor:\n", "3×3 Array{Float64,2}:\n", " -1.28564 -0.681455 -0.56284\n", " 0.0 0.643045 0.244737\n", " 0.0 0.0 -0.231509" ] }, "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,Array{Float64,2}}:\n", " -0.0703771 0.655875 -0.751581\n", " -0.749523 -0.531943 -0.394021\n", " -0.658227 0.535597 0.52903" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Aqr.Q" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " -1.28564 -0.681455 -0.56284\n", " 0.0 0.643045 0.244737\n", " 0.0 0.0 -0.231509" ] }, "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,Array{Float64,2},Array{Float64,1}}\n", "values:\n", "3-element Array{Float64,1}:\n", " -1.18266124303562\n", " -0.5242093942277856\n", " 2.98340018375136\n", "vectors:\n", "3×3 Array{Float64,2}:\n", " 0.755935 -0.346775 0.555257\n", " -0.649515 -0.503283 0.569944\n", " -0.0818095 0.791488 0.605685" ] }, "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 Array{Float64,1}:\n", " -1.18266124303562\n", " -0.5242093942277856\n", " 2.98340018375136" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "AsymEig.values" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 0.755935 -0.346775 0.555257\n", " -0.649515 -0.503283 0.569944\n", " -0.0818095 0.791488 0.605685" ] }, "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 Array{Float64,2}:\n", " 1.0 9.99201e-16 1.44329e-15\n", " 8.88178e-16 1.0 1.44329e-15\n", " -8.88178e-16 -6.66134e-16 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": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1000×1000 Array{Float64,2}:\n", " -0.54116 0.723425 -0.454392 … 0.153621 0.49022 2.57411\n", " 1.06068 -2.12598 0.578722 -1.75703 -0.694434 0.0692548\n", " -0.975224 -0.897821 -1.38972 -1.66046 1.5895 -2.62675\n", " -0.795238 -0.637253 -0.414014 -0.620459 -2.76297 -1.68696\n", " 1.19269 -0.409572 0.664132 0.397247 1.43729 1.20937\n", " 0.0898626 -0.573358 -0.64753 … 0.5808 0.44763 -0.463222\n", " 2.60072 -1.40748 -1.50094 1.21125 -1.40829 1.07843\n", " -0.0354091 0.21075 -0.344279 2.0649 -0.0658153 0.437195\n", " 2.13631 0.837857 1.33064 1.02112 -1.99292 -1.03814\n", " 0.327448 -0.690153 0.194605 -2.01814 0.152014 -0.640831\n", " -0.229009 -1.23092 1.62548 … 0.00552164 -0.283409 0.49161\n", " 1.1002 -0.722856 -2.43117 0.639876 -0.775289 0.348139\n", " -0.277119 -0.87718 2.16643 -0.84017 -0.177649 0.356069\n", " ⋮ ⋱ \n", " -1.14988 0.208948 0.347219 1.03618 -3.71452 -0.831069\n", " -1.41506 0.219713 -0.292957 1.38114 -1.97914 -0.257589\n", " 0.135981 0.653265 -1.67877 … 0.640539 -2.38161 0.250336\n", " 0.487141 0.353631 0.209629 -1.12091 -0.123502 1.14744\n", " 0.299761 0.353541 0.368271 0.595472 -0.12664 0.74589\n", " -0.152775 0.408049 0.134503 -0.951019 1.78102 0.105328\n", " 0.544877 -0.199899 -0.115924 -1.51412 0.280126 -2.05524\n", " 2.22375 -1.07693 0.543478 … 0.920955 1.13643 0.576951\n", " 1.17054 1.56297 -0.208841 -0.562418 0.327011 0.299262\n", " 0.619724 -0.595299 -0.458258 0.768808 -2.68816 0.173424\n", " 0.566441 -0.847296 -0.697741 0.656236 1.23601 0.366564\n", " 0.343866 -0.296927 -1.5821 0.586209 1.86195 0.574866" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "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": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "true" ] }, "execution_count": 20, "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": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.7841028548025857" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Asym_noisy = copy(Asym)\n", "Asym_noisy[1,2] += 5eps()" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "false" ] }, "execution_count": 22, "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": 23, "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": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.206384 seconds (168.04 k allocations: 16.062 MiB)\n" ] } ], "source": [ "@time eigvals(Asym);" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.507267 seconds (14 allocations: 7.921 MiB)\n" ] } ], "source": [ "@time eigvals(Asym_noisy);" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.141559 seconds (12 allocations: 7.989 MiB, 22.12% gc 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": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.603548 seconds (473.51 k allocations: 206.850 MiB, 22.09% gc time)\n" ] }, { "data": { "text/plain": [ "6.430820033538295" ] }, "execution_count": 28, "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": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1//2" ] }, "execution_count": 29, "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": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Rational{BigInt},2}:\n", " 9//10 1//5 1//2\n", " 1//1 1//10 1//2\n", " 1//2 9//10 4//5" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Arational = Matrix{Rational{BigInt}}(rand(1:10, 3, 3))/10" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Rational{BigInt},1}:\n", " 8//5\n", " 8//5\n", " 11//5" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = fill(1, 3)\n", "b = Arational*x" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Rational{BigInt},1}:\n", " 1//1\n", " 1//1\n", " 1//1" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Arational\\b" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "LU{Rational{BigInt},Array{Rational{BigInt},2}}\n", "L factor:\n", "3×3 Array{Rational{BigInt},2}:\n", " 1//1 0//1 0//1\n", " 1//2 1//1 0//1\n", " 9//10 11//85 1//1\n", "U factor:\n", "3×3 Array{Rational{BigInt},2}:\n", " 1//1 1//10 1//2\n", " 0//1 17//20 11//20\n", " 0//1 0//1 -9//425" ] }, "execution_count": 33, "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": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Array{Int64,2}:\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" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using LinearAlgebra\n", "A = [\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", "]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5-element Array{Float64,1}:\n", " -128.49322764802145\n", " -55.887784553056875\n", " 42.7521672793189\n", " 87.16111477514521\n", " 542.4677301466143" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A_eigv = eigvals(A)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "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` and store this in the variable `A_diag`." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Diagonal{Float64,Array{Float64,1}}:\n", " -128.493 ⋅ ⋅ ⋅ ⋅ \n", " ⋅ -55.8878 ⋅ ⋅ ⋅ \n", " ⋅ ⋅ 42.7522 ⋅ ⋅ \n", " ⋅ ⋅ ⋅ 87.1611 ⋅ \n", " ⋅ ⋅ ⋅ ⋅ 542.468" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A_diag = Diagonal(A_eigv)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "ename": "AssertionError", "evalue": "AssertionError: A_diag == [-128.493 0 0 0 0; 0 -55.8878 0 0 0; 0 0 42.7522 0 0; 0 0 0 87.1611 0; 0 0 0 0 542.468]", "output_type": "error", "traceback": [ "AssertionError: A_diag == [-128.493 0 0 0 0; 0 -55.8878 0 0 0; 0 0 42.7522 0 0; 0 0 0 87.1611 0; 0 0 0 0 542.468]", "", "Stacktrace:", " [1] top-level scope at In[22]:1" ] } ], "source": [ "@assert A_diag == [-128.493 0 0 0 0 \n", " 0 -55.8878 0 0 0 \n", " 0 0 42.7522 0 0 \n", " 0 0 0 87.1611 0 \n", " 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": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 LowerTriangular{Int64,Array{Int64,2}}:\n", " 140 ⋅ ⋅ ⋅ ⋅\n", " 97 106 ⋅ ⋅ ⋅\n", " 74 89 152 ⋅ ⋅\n", " 168 131 144 54 ⋅\n", " 131 36 71 142 36" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A_lowertri = LowerTriangular(A)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "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": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.4", "language": "julia", "name": "julia-1.3" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.4.0" } }, "nbformat": 4, "nbformat_minor": 4 }