{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "

1  Gaussian Elimination and LU Decomposition
1.1  Elementary operator matrix
1.2  Gauss transformations
1.3  LU decomposition
1.4  Matrix inversion
1.5  Pivoting
1.6  Implementation
1.7  Further reading
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Gaussian Elimination and LU Decomposition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Goal: solve linear equation\n", "$$\n", "\\mathbf{A} \\mathbf{x} = \\mathbf{b}.\n", "$$\n", "For simplicity we consider a square matrix $\\mathbf{A} \\in \\mathbb{R}^{n \\times n}$.\n", "\n", "* History: Chinese mathematical text [The Nine Chapters on the Mathematical Art](https://en.wikipedia.org/wiki/The_Nine_Chapters_on_the_Mathematical_Art), Issac Newton and Carl Friedrich Gauss.\n", "\n", "* A [toy example](https://en.wikipedia.org/wiki/Gaussian_elimination#Example_of_the_algorithm)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 2.0 1.0 -1.0\n", " -3.0 -1.0 2.0\n", " -2.0 1.0 2.0" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = [2.0 1.0 -1.0; -3.0 -1.0 2.0; -2.0 1.0 2.0]" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 8.0\n", " -11.0\n", " -3.0" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b = [8.0; -11.0; -3.0]" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 2.0\n", " 3.0\n", " -1.0" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# solve linear equation\n", "A \\ b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What happens when we call `A \\ b` to solve a linear equation?" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Elementary operator matrix\n", "\n", "* **Elementary operator matrix** is the identity matrix with the 0 in position $(j,k)$ replaced by $c$\n", "$$\n", " \\mathbf{E}_{jk}(c) = \\begin{pmatrix}\n", " 1 & & \\\\\n", " & \\ddots & \\\\\n", " & & 1 & \\\\\n", " & & & \\ddots & \\\\\n", " & & c & & 1 & \\\\\n", " & & & & & \\ddots \\\\\n", " & & & & & & 1\n", " \\end{pmatrix} = \\mathbf{I} + c \\mathbf{e}_j \\mathbf{e}_k^T.\n", "$$\n", "\n", "* $\\mathbf{E}_{jk}(c)$ is unit triangular, full rank, and its inverse is\n", "$$\n", " \\mathbf{E}_{jk}^{-1}(c) = \\mathbf{E}_{jk}(-c).\n", "$$\n", "\n", "* $\\mathbf{E}_{jk}(c)$ left-multiplies an $n \\times m$ matrix $\\mathbf{X}$ effectively replace its $j$-th row $\\mathbf{x}_{j\\cdot}$ by $c \\mathbf{x}_{k \\cdot} + \\mathbf{x}_{j \\cdot}$\n", "$$\n", " \\mathbf{E}_{jk}(c) \\times \\mathbf{X} = \\mathbf{E}_{jk}(c) \\times \\begin{pmatrix}\n", " & & \\\\\n", " \\cdots & \\mathbf{x}_{k\\cdot} & \\cdots \\\\\n", " & & \\\\\n", " \\cdots & \\mathbf{x}_{j\\cdot} & \\cdots \\\\\n", " & & \n", " \\end{pmatrix} = \\begin{pmatrix}\n", " & & \\\\\n", " \\cdots & \\mathbf{x}_{k\\cdot} & \\cdots \\\\\n", " & & \\\\\n", " \\cdots & c \\mathbf{x}_{k\\cdot} + \\mathbf{x}_{j\\cdot} & \\cdots \\\\\n", " & & \n", " \\end{pmatrix}.\n", "$$\n", "$2m$ flops.\n", "\n", "* Gaussian elimination applies a sequence of elementary operator matrices to transform the linear system $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$ to an upper triangular system\n", "$$\n", "\\begin{eqnarray*}\n", " \\mathbf{E}_{n,n-1}(c_{n,n-1}) \\cdots \\mathbf{E}_{21}(c_{21}) \\mathbf{A} \\mathbf{x} &=& \\mathbf{E}_{n,n-1}(c_{n,n-1}) \\cdots \\mathbf{E}_{21}(c_{21}) \\mathbf{b} \\\\\n", " \\mathbf{U} \\mathbf{x} &=& \\mathbf{b}_{\\text{new}}.\n", "\\end{eqnarray*}\n", "$$\n", " \n", " Column 1:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 0.0 0.0\n", " 1.5 1.0 0.0\n", " 0.0 0.0 1.0" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "E21 = [1.0 0.0 0.0; 1.5 1.0 0.0; 0.0 0.0 1.0]" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 2.0 1.0 -1.0\n", " 0.0 0.5 0.5\n", " -2.0 1.0 2.0" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# zero (2, 1) entry\n", "E21 * A" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 0.0 0.0\n", " 0.0 1.0 0.0\n", " 1.0 0.0 1.0" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "E31 = [1.0 0.0 0.0; 0.0 1.0 0.0; 1.0 0.0 1.0]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 2.0 1.0 -1.0\n", " 0.0 0.5 0.5\n", " 0.0 2.0 1.0" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# zero (3, 1) entry\n", "E31 * E21 * A" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Column 2:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 0.0 0.0\n", " 0.0 1.0 0.0\n", " 0.0 -4.0 1.0" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "E32 = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 -4.0 1.0]" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 2.0 1.0 -1.0\n", " 0.0 0.5 0.5\n", " 0.0 0.0 -1.0" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# zero (3, 2) entry\n", "E32 * E31 * E21 * A" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gauss transformations\n", "\n", "* For the first column,\n", "$$\n", " \\mathbf{M}_1 = \\mathbf{E}_{n1}(c_{n1}) \\cdots \\mathbf{E}_{31}(c_{31}) \\mathbf{E}_{21}(c_{21}) = \\begin{pmatrix}\n", " 1 & \\\\\n", " c_{21} & \\\\\n", " & \\ddots & \\\\\n", " c_{n1} & & 1\n", " \\end{pmatrix}\n", "$$ \n", "For the $k$-th column,\n", "$$\n", "\t\\mathbf{M}_k = \\mathbf{E}_{nk}(c_{nk}) \\cdots \\mathbf{E}_{k+1,k}(c_{k+1,k}) = \\begin{pmatrix}\n", "\t1 & \\\\\n", "\t& \\ddots \\\\\n", "\t& & 1 & \\\\\n", "\t& & c_{k+1,k} & 1\\\\\n", "\t& & \\vdots & & \\ddots \\\\\n", "\t& & c_{n,k} & & & 1\n", "\t\\end{pmatrix}.\n", "$$\n", "\n", "* $\\mathbf{M}_1, \\ldots, \\mathbf{M}_{n-1}$ are called the **Gauss transformations**." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 0.0 0.0\n", " 1.5 1.0 0.0\n", " 1.0 0.0 1.0" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "M1 = E31 * E21" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 0.0 0.0\n", " 0.0 1.0 0.0\n", " 0.0 -4.0 1.0" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "M2 = E32" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Gauss transformations $\\mathbf{M}_k$ are unit triangular, full rank, with inverse\n", "$$\n", "\t\\mathbf{M}_k^{-1} = \\mathbf{E}_{k+1,k}^{-1}(c_{k+1,k}) \\cdots \\mathbf{E}_{nk}^{-1}(c_{nk}) = \\begin{pmatrix}\n", "\t1 & \\\\\n", "\t& \\ddots \\\\\n", "\t& & 1 & \\\\\n", "\t& & - c_{k+1,k}\\\\\n", "\t& & \\vdots & & \\ddots \\\\\n", "\t& & - c_{n,k} & & & 1\n", "\t\\end{pmatrix}.\n", "$$" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 0.0 0.0\n", " -1.5 1.0 0.0\n", " -1.0 0.0 1.0" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inv(M1)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 0.0 0.0\n", " 0.0 1.0 0.0\n", " 0.0 4.0 1.0" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inv(M2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## LU decomposition\n", "\n", "Gaussian elimination does\n", "$$\n", "\\mathbf{M}_{n-1} \\cdots \\mathbf{M}_1 \\mathbf{A} = \\mathbf{U}.\n", "$$ \n", "Let" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{equation*}\n", "\\mathbf{L} = \\mathbf{M}_1^{-1} \\cdots \\mathbf{M}_{n-1}^{-1} = \\begin{pmatrix} \n", "\t1 & & & & \\\\\n", "\t- c_{21} & \\ddots & & & \\\\\n", "\t& & 1 & & \\\\\n", "\t- c_{k+1,1} & & - c_{k+1,k} & & \\\\\n", "\t& & \\vdots & & \\ddots \\\\\n", "\t- c_{n1} & & - c_{nk} & & & 1\n", "\t\\end{pmatrix}.\n", "\\end{equation*}\n", "Thus we have the **LU decomposition**\n", "$$\n", "\\mathbf{A} = \\mathbf{L} \\mathbf{U},\n", "$$\n", "where $\\mathbf{L}$ is unit lower triangular and $\\mathbf{U}$ is upper triangular.\n", "\n", "* The whole LU algorithm is done in place, i.e., $\\mathbf{A}$ is overwritten by $\\mathbf{L}$ and $\\mathbf{U}$.\n", "\n", "* LU decomposition exists if the principal sub-matrix $\\mathbf{A}[1:k, 1:k]$ is non-singular for $k=1,\\ldots,n-1$.\n", "\n", "* If the LU decomposition exists and $\\mathbf{A}$ is non-singular, then the LU decmpositon is unique and\n", "$$\n", " \\det(\\mathbf{A}) = \\det(\\mathbf{L}) \\det(\\mathbf{U}) = \\prod_{k=1}^n u_{kk}.\n", "$$\n", "\n", "* The LU decomposition costs\n", "$$\n", " 2(n-1)^2 + 2(n-2)^2 + \\cdots + 2 \\cdot 1^2 \\approx \\frac 23 n^3 \\quad \\text{flops}.\n", "$$\n", "\n", "\n", "\n", "* Actual implementations can differ: outer product LU ($kij$ loop), block outer product LU (higher level-3 fraction), Crout's algorithm ($jki$ loop).\n", "\n", "* Given LU decomposition $\\mathbf{A} = \\mathbf{L} \\mathbf{U}$, solving $(\\mathbf{L} \\mathbf{U}) \\mathbf{x} = \\mathbf{b}$ for one right hand side costs $2n^2$ flops:\n", " - One forward substitution ($n^2$ flops) to solve\n", " $$\n", " \\mathbf{L} \\mathbf{y} = \\mathbf{b}\n", " $$\n", " - One back substitution ($n^2$ flops) to solve\n", " $$\n", " \\mathbf{U} \\mathbf{x} = \\mathbf{y}\n", " $$\n", " \n", "* Therefore to solve $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$ via LU, we need a total of\n", "$$\n", " \\frac 23 n^3 + 2n^2 \\quad \\text{flops}.\n", "$$\n", "\n", "* If there are multiple right hand sides, LU only needs to be done once.\n", "\n", "\n", "## Matrix inversion\n", "\n", "* For matrix inversion, there are $n$ right hand sides $\\mathbf{e}_1, \\ldots, \\mathbf{e}_n$. Taking advantage of 0s reduces $2n^3$ flops to $\\frac 43 n^3$ flops. So **matrix inversion via LU** costs\n", "$$\n", "\\frac 23 n^3 + \\frac 43 n^3 = 2n^3 \\quad \\text{flops}.\n", "$$\n", "\n", "* **No inversion mentality**: \n", " > Whenever we see matrix inverse, we should think in terms of solving linear equations. \n", "\n", " We do not compute matrix inverse unless \n", " 0. it is necessary to compute standard errors \n", " 0. number of right hand sides is much larger than $n$ \n", " 0. $n$ is small\n", " \n", "## Pivoting \n", "\n", "* Let\n", "$$\n", " \\mathbf{A} = \\begin{pmatrix}\n", " 0 & 1 \\\\\n", " 1 & 0 \\\\\n", " \\end{pmatrix}.\n", "$$\n", "Is there a solution to $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$ for an arbitrary $\\mathbf{b}$? Does GE/LU work for $\\mathbf{A}$?\n", "\n", "* What if, during LU procedure, the **pivot** $a_{kk}$ is 0 or nearly 0 due to underflow? \n", " Solution: pivoting.\n", "\n", "* **Partial pivoting**: before zeroing the $k$-th column, the row with $\\max_{i=k}^n |a_{ik}|$ is moved into the $k$-th row. \n", "\n", "* LU with partial pivoting yields\n", "$$\n", "\t\\mathbf{P} \\mathbf{A} = \\mathbf{L} \\mathbf{U},\n", "$$\n", "where $\\mathbf{P}$ is a permutation matrix, $\\mathbf{L}$ is unit lower triangular with $|\\ell_{ij}| \\le 1$, and $\\mathbf{U}$ is upper triangular.\n", "\n", "* Complete pivoting: Do both row and column interchanges so that the largest entry in the sub matrix `A[k:n, k:n]` is permuted to the $(k,k)$-th entry. This yields the decomposition \n", "$$\n", "\\mathbf{P} \\mathbf{A} \\mathbf{Q} = \\mathbf{L} \\mathbf{U},\n", "$$\n", "where $|\\ell_{ij}| \\le 1$.\n", "\n", "* LU decomposition with partial pivoting is the most commonly used methods for solving **general** linear systems. Complete pivoting is the most stable but costs more computation. Partial pivoting is stable most of times.\n", "\n", "## Implementation\n", "\n", "* LAPACK: [`?GETRF`](http://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga0019443faea08275ca60a734d0593e60.html#ga0019443faea08275ca60a734d0593e60) does $\\mathbf{P} \\mathbf{A} = \\mathbf{L} \\mathbf{U}$ (LU decomposition with partial pivoting) in place.\n", "\n", "* R: `solve()` implicitly performs LU decomposition (wrapper of LAPACK routine `DGESV`). `solve()` allows specifying a single or multiple right hand sides. If none, it computes the matrix inverse. The `matrix` package contains `lu()` function that outputs `L`, `U`, and `pivot`.\n", "\n", "* Julia: \n", " - [`lufact`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.lufact), [`lufact!`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.lufact!), [`lu`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.lu)\n", " - Or call LAPACK wrapper function [`getrf!`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.LAPACK.getrf!) directly\n", " - Other LU-related LAPACK wrapper functions: [`gesv`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.LAPACK.gesv!), [`gesvx`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.LAPACK.gesvx!), [`trtri`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.LAPACK.trtri!) (inverse of triangular matrix), [`trtrs`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.LAPACK.trtrs!)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 2.0 1.0 -1.0\n", " -3.0 -1.0 2.0\n", " -2.0 1.0 2.0" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Base.LinAlg.LU{Float64,Array{Float64,2}}" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# second argument indicates partial pivoting (default) or not\n", "alu = lufact(A, Val{true})\n", "typeof(alu)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 1.0 0.0 0.0\n", " 0.666667 1.0 0.0\n", " -0.666667 0.2 1.0" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "alu[:L]" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " -3.0 -1.0 2.0 \n", " 0.0 1.66667 0.666667\n", " 0.0 0.0 0.2 " ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "alu[:U]" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Int64,1}:\n", " 2\n", " 3\n", " 1" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "alu[:p]" ] }, { "cell_type": "code", "execution_count": 19, "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": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "alu[:P]" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " -3.0 -1.0 2.0\n", " -2.0 1.0 2.0\n", " 2.0 1.0 -1.0" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "alu[:L] * alu[:U]" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " -3.0 -1.0 2.0\n", " -2.0 1.0 2.0\n", " 2.0 1.0 -1.0" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A[alu[:p], :]" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 2.0\n", " 3.0\n", " -1.0" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# this is doing two triangular solves\n", "alu \\ b" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.9999999999999998" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "det(A) # this does LU!" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.9999999999999998" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "det(alu) # this is cheap" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 4.0 3.0 -1.0\n", " -2.0 -2.0 1.0\n", " 5.0 4.0 -1.0" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inv(A) # this does LU!" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 4.0 3.0 -1.0\n", " -2.0 -2.0 1.0\n", " 5.0 4.0 -1.0" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inv(alu) # this is cheap" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Further reading\n", "\n", "* Sections II.5.2 and II.5.3 of [Computational Statistics](http://ucla.worldcat.org/title/computational-statistics/oclc/437345409&referer=brief_results) by James Gentle (2010).\n", "\n", "* A definite reference is Chapter 3 of the book [Matrix Computation](http://catalog.library.ucla.edu/vwebv/holdingsInfo?bibId=7122088) by Gene Golub and Charles Van Loan.\n", "\n", "\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Julia Version 0.6.2\n", "Commit d386e40c17 (2017-12-13 18:08 UTC)\n", "Platform Info:\n", " OS: macOS (x86_64-apple-darwin14.5.0)\n", " CPU: Intel(R) Core(TM) i7-6920HQ CPU @ 2.90GHz\n", " WORD_SIZE: 64\n", " BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)\n", " LAPACK: libopenblas64_\n", " LIBM: libopenlibm\n", " LLVM: libLLVM-3.9.1 (ORCJIT, skylake)\n" ] } ], "source": [ "versioninfo()" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.6.2", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.2" }, "toc": { "colors": { "hover_highlight": "#DAA520", "running_highlight": "#FF0000", "selected_highlight": "#FFD700" }, "moveMenuLeft": true, "nav_menu": { "height": "30px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "threshold": 4, "toc_cell": true, "toc_section_display": "block", "toc_window_display": true, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 2 }