{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "R version 4.2.2 (2022-10-31)\n", "Platform: x86_64-apple-darwin17.0 (64-bit)\n", "Running under: macOS Catalina 10.15.7\n", "\n", "Matrix products: default\n", "BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib\n", "LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib\n", "\n", "locale:\n", "[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8\n", "\n", "attached base packages:\n", "[1] stats graphics grDevices utils datasets methods base \n", "\n", "loaded via a namespace (and not attached):\n", " [1] fansi_1.0.4 crayon_1.5.2 digest_0.6.31 utf8_1.2.2 \n", " [5] IRdisplay_1.1 repr_1.1.5 lifecycle_1.0.3 jsonlite_1.8.4 \n", " [9] evaluate_0.20 pillar_1.8.1 rlang_1.0.6 cli_3.6.0 \n", "[13] uuid_1.1-0 vctrs_0.5.2 IRkernel_1.3.2 tools_4.2.2 \n", "[17] glue_1.6.2 fastmap_1.1.0 compiler_4.2.2 base64enc_0.1-3\n", "[21] pbdZMQ_0.3-9 htmltools_0.5.4" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sessionInfo()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Required R packages" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [], "source": [ "library(Matrix)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Computer methods for solving linear systems of equations\n", "\n", "We consider computer algorithms for solving linear equations $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$, a ubiquitous task in statistics. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Our mantra\n", "\n", "> **The form of a mathematical expression and the way the expression should be evaluated in actual practice may be quite different.**\n", "\\\n", " -- James Gentle, *Matrix Algebra*, Springer, New York (2007)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Algorithms\n", "\n", "* Algorithm is loosely defined as a set of instructions for doing something, which terminates in finite time. An algorithm requires input and output.\n", "\n", "\n", "### Measure of algorithm efficiency \n", "\n", "* A basic unit for measuring algorithmic efficiency is **flop**. \n", "> A flop (**floating point operation**) consists of a floating point addition, subtraction, multiplication, division, or comparison, and the usually accompanying fetch and store. \n", "\n", "* How to measure efficiency of an algorithm? Big O notation. If $n$ is the size of a problem, an algorithm has order $O(f(n))$, where the leading term in the number of flops is $c \\cdot f(n)$. For example,\n", " - matrix-vector multiplication `A * b`, where `A` is $m \\times n$ and `b` is $n \\times 1$, takes $2mn$ or $O(mn)$ flops \n", " - matrix-matrix multiplication `A * B`, where `A` is $m \\times n$ and `B` is $n \\times p$, takes $2mnp$ or $O(mnp)$ flops\n", "\n", "* Typical orders of computational complexity: \n", " - Exponential order: $O(b^n)$ (\"horrible\") \n", " - Polynomial order: $O(n^q)$ (doable) \n", " - $O(n \\log n )$ (fast) \n", " - Linear order $O(n)$ (fast) \n", " - Logarithmic order $O(\\log n)$ (super fast) \n", " \n", "* Difference of $O(n^2)$ and $O(n\\log n)$ on massive data. Suppose we have a teraflops supercomputer capable of doing $10^{12}$ flops per second. For a problem of size $n=10^{12}$, $O(n \\log n)$ algorithm takes about \n", "$$10^{12} \\log (10^{12}) / 10^{12} \\approx 27 \\text{ seconds}.$$ \n", "$O(n^2)$ algorithm takes about $10^{12}$ seconds, which is approximately 31710 years!\n", "\n", "* The Quicksort and Fast Fourier Transform (invented by John Tukey--recall \"bit\"--) are celebrated algorithms that turn $O(n^2)$ operations into $O(n \\log n)$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Triangular systems\n", "\n", "Idea: turning original problem into an **easy** one, e.g., triangular system." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Lower triangular system\n", "\n", "To solve $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$, where $\\mathbf{A} \\in \\mathbb{R}^{n \\times n}$ is **lower triangular**\n", "\n", "$$\n", "\\begin{bmatrix}\n", " a_{11} & 0 & \\cdots & 0 \\\\\n", " a_{21} & a_{22} & \\cdots & 0 \\\\\n", " \\vdots & \\vdots & \\ddots & \\vdots \\\\\n", " a_{n1} & a_{n2} & \\cdots & a_{nn}\n", "\\end{bmatrix}\n", "\\begin{bmatrix}\n", "x_1 \\\\ x_2 \\\\ \\vdots \\\\ x_n\n", "\\end{bmatrix} = \\begin{bmatrix}\n", "b_1 \\\\ b_2 \\\\ \\vdots \\\\ b_n\n", "\\end{bmatrix}.\n", "$$\n", "\n", "* **Forward substitution**: \n", "$$\n", "\\begin{eqnarray*}\n", " x_1 &=& b_1 / a_{11} \\\\\n", " x_2 &=& (b_2 - a_{21} x_1) / a_{22} \\\\\n", " x_3 &=& (b_3 - a_{31} x_1 - a_{32} x_2) / a_{33} \\\\\n", " &\\vdots& \\\\\n", " x_n &=& (b_n - a_{n1} x_1 - a_{n2} x_2 - \\cdots - a_{n,n-1} x_{n-1}) / a_{nn}.\n", "\\end{eqnarray*}\n", "$$\n", "\n", "* $1 + 3 + 5 + \\cdots + (2n-1) = n^2$ flops. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Upper triangular system\n", "\n", "To solve $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$, where $\\mathbf{A} \\in \\mathbb{R}^{n \\times n}$ is upper triangular \n", "$$\n", "\\begin{bmatrix}\n", " a_{11} & \\cdots & a_{1,n-1} & a_{1n} \\\\\n", " \\vdots & \\ddots & \\vdots & \\vdots \\\\\n", " 0 & \\cdots & a_{n-1,n-1} & a_{n-1,n} \\\\\n", " 0 & 0 & 0 & a_{nn}\n", "\\end{bmatrix}\n", "\\begin{bmatrix}\n", "x_1 \\\\ \\vdots \\\\ x_{n-1} \\\\ x_n\n", "\\end{bmatrix} = \\begin{bmatrix}\n", "b_1 \\\\ \\vdots \\\\ b_{n-1} \\\\ b_n\n", "\\end{bmatrix}.\n", "$$\n", "\n", "* **Back substitution** (backsolve): \n", "$$\n", "\\begin{eqnarray*}\n", " x_n &=& b_n / a_{nn} \\\\\n", " x_{n-1} &=& (b_{n-1} - a_{n-1,n} x_n) / a_{n-1,n-1} \\\\\n", " x_{n-2} &=& (b_{n-2} - a_{n-2,n-1} x_{n-1} - a_{n-2,n} x_n) / a_{n-2,n-2} \\\\\n", " &\\vdots& \\\\\n", " x_1 &=& (b_1 - a_{12} x_2 - a_{13} x_3 - \\cdots - a_{1,n} x_{n}) / a_{11}.\n", "\\end{eqnarray*}\n", "$$\n", "\n", "* $n^2$ flops." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementation\n", "\n", "* Most numerical computing software packages (including R) are built on top of the BLAS (_basic linear algebra subprograms_). \n", " - [Netlib](http://www.netlib.org/blas/) provides a reference implementation. \n", " - Vanila R (the binary that you can download from ) uses its own implementation of BLAS (see `sessionInfo()` above)\n", "\n", "* The BLAS triangular system solve is done *in place*, i.e., $\\mathbf{b}$ is **overwritten** by $\\mathbf{x}$.\n", "```r\n", " # forward substitution\n", " for (i in seq_len(n)) {\n", " for (j in seq_len(i-1)) {\n", " b[i] <- b[i] - A[i, j] * b[j]\n", " }\n", " b[i] <- b[i] / A[i, i]\n", " }\n", " # backsolve\n", " for (i in rev(seq_len(n))) {\n", " for (j in i + seq_len(max(n - i, 0))) {\n", " b[i] <- b[i] - A[i, j] * b[j]\n", " }\n", " b[i] <- b[i] / A[i, i]\n", " }\n", "```\n", "\n", "* R\n", " - `?backsolve`, `?forwardsolve`\n", " - This is a wrapper for the level-3 BLAS routine [`dtrsm`](http://www.netlib.org/lapack/explore-html/d1/d54/group__double__blas__level3_ga6a0a7704f4a747562c1bd9487e89795c.html#ga6a0a7704f4a747562c1bd9487e89795c)\n", " - The `d` stands for *double precision*." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "5 x 5 Matrix of class \"dgeMatrix\"\n", " [,1] [,2] [,3] [,4] [,5]\n", "[1,] -0.56047565 1.7150650 1.2240818 1.7869131 -1.0678237\n", "[2,] -0.23017749 0.4609162 0.3598138 0.4978505 -0.2179749\n", "[3,] 1.55870831 -1.2650612 0.4007715 -1.9666172 -1.0260044\n", "[4,] 0.07050839 -0.6868529 0.1106827 0.7013559 -0.7288912\n", "[5,] 0.12928774 -0.4456620 -0.5558411 -0.4727914 -0.6250393" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "set.seed(123) # seed\n", "n <- 5\n", "A <- Matrix(rnorm(n * n), nrow = n)\n", "b <- rnorm(n)\n", "# a random matrix\n", "A" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "5 x 5 Matrix of class \"dgeMatrix\"\n", " [,1] [,2] [,3] [,4] [,5]\n", "[1,] -0.5604756 1.7150650 1.2240818 1.7869131 -1.0678237\n", "[2,] 0.0000000 0.4609162 0.3598138 0.4978505 -0.2179749\n", "[3,] 0.0000000 0.0000000 0.4007715 -1.9666172 -1.0260044\n", "[4,] 0.0000000 0.0000000 0.0000000 0.7013559 -0.7288912\n", "[5,] 0.0000000 0.0000000 0.0000000 0.0000000 -0.6250393" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Al <- lower.tri(A) # Returns a matrix of logicals the same size of a given matrix with entries TRUE in the lower triangle\n", "Aupper <- A\n", "Aupper[Al] <- 0\n", "Aupper" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/html": [ "\n", "
  1. 14.6233350267841
  2. 22.7861639334402
  3. -22.9457487682144
  4. -3.70749941033203
  5. -2.00597784101513
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 14.6233350267841\n", "\\item 22.7861639334402\n", "\\item -22.9457487682144\n", "\\item -3.70749941033203\n", "\\item -2.00597784101513\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 14.6233350267841\n", "2. 22.7861639334402\n", "3. -22.9457487682144\n", "4. -3.70749941033203\n", "5. -2.00597784101513\n", "\n", "\n" ], "text/plain": [ "[1] 14.623335 22.786164 -22.945749 -3.707499 -2.005978" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "backsolve(Aupper, b) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Some algebraic facts of triangular system\n", "\n", "* Eigenvalues of a triangular matrix $\\mathbf{A}$ are diagonal entries $\\lambda_i = a_{ii}$. \n", "\n", "* Determinant $\\det(\\mathbf{A}) = \\prod_i a_{ii}$.\n", "\n", "* The product of two upper (lower) triangular matrices is upper (lower) triangular.\n", "\n", "* The inverse of an upper (lower) triangular matrix is upper (lower) triangular.\n", "\n", "* The product of two unit upper (lower) triangular matrices is unit upper (lower) triangular.\n", "\n", "* The inverse of a unit upper (lower) triangular matrix is unit upper (lower) triangular." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Gaussian Elimination and LU Decomposition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* A [toy example](./gelu.pdf).\n", "\n", "* [History](https://en.wikipedia.org/wiki/Gaussian_elimination#History): the method is named after Carl Friedrich Gauss (1777–1855), although it stems from the notes of Isaac Newton. If fact, GE was known to [Chinese mathematicians](https://en.wikipedia.org/wiki/The_Nine_Chapters_on_the_Mathematical_Art) as early as 179 AD." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dgeMatrix\"\n", " [,1] [,2] [,3]\n", "[1,] 2 -4 2\n", "[2,] 4 -9 7\n", "[3,] 2 1 3" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "A <- t(Matrix(c(2.0, -4.0, 2.0, 4.0, -9.0, 7.0, 2.0, 1.0, 3.0), ncol=3)) # R uses column major ordering\n", "A" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [], "source": [ "b <- c(6.0, 20.0, 14.0)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 1 Matrix of class \"dgeMatrix\"\n", " [,1]\n", "[1,] 2\n", "[2,] 1\n", "[3,] 3" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# R's way tp solve linear equation\n", "\n", "solve(A, b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What happens when we call `solve(A, b)`?" ] }, { "cell_type": "markdown", "metadata": {}, "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{bmatrix}\n", " 1 & & \\\\\n", " & \\ddots & \\\\\n", " & & 1 & \\\\\n", " & & & \\ddots & \\\\\n", " & & c & & 1 & \\\\\n", " & & & & & \\ddots \\\\\n", " & & & & & & 1\n", " \\end{bmatrix} = \\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{bmatrix}\n", " & & \\\\\n", " \\cdots & \\mathbf{x}_{k\\cdot} & \\cdots \\\\\n", " & & \\\\\n", " \\cdots & \\mathbf{x}_{j\\cdot} & \\cdots \\\\\n", " & & \n", " \\end{bmatrix} = \\begin{bmatrix}\n", " & & \\\\\n", " \\cdots & \\mathbf{x}_{k\\cdot} & \\cdots \\\\\n", " & & \\\\\n", " \\cdots & c \\mathbf{x}_{k\\cdot} + \\mathbf{x}_{j\\cdot} & \\cdots \\\\\n", " & & \n", " \\end{bmatrix}.\n", "$$\n", "$2m$ flops (why?).\n", "\n", "* Gaussian elimination can be understood as applying 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", "Let's get back to our toy example.\n", " \n", "### Column 1" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "3 x 3 sparse Matrix of class \"dtCMatrix\"\n", " \n", "[1,] 1 . .\n", "[2,] -2 1 .\n", "[3,] . . 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "E21 <- t(Matrix(c(1.0, 0.0, 0.0, -2.0, 1.0, 0.0, 0.0, 0.0, 1.0), nrow=3)) # Step 1, inv(L1)\n", "E21" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dgeMatrix\"\n", " [,1] [,2] [,3]\n", "[1,] 2 -4 2\n", "[2,] 0 -1 3\n", "[3,] 2 1 3" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# zero (2, 1) entry \n", "E21 %*% A # Step 1, A'" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true }, "tags": [] }, "outputs": [ { "data": { "text/plain": [ "3 x 3 sparse Matrix of class \"dtCMatrix\"\n", " \n", "[1,] 1 . .\n", "[2,] . 1 .\n", "[3,] -1 . 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "E31 <- t(Matrix(c(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, 1.0), nrow=3)) # Step 2, find the corresponding matrix!\n", "E31" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dgeMatrix\"\n", " [,1] [,2] [,3]\n", "[1,] 2 -4 2\n", "[2,] 0 -1 3\n", "[3,] 0 5 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# zero (3, 1) entry \n", "E31 %*% E21 %*% A # Step2, A''" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Column 2" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 sparse Matrix of class \"dtCMatrix\"\n", " \n", "[1,] 1 . .\n", "[2,] . 1 .\n", "[3,] . 5 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "E32 <- t(Matrix(c(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 5.0, 1.0), nrow=3)) # Step 3, find the corresponding matrix!\n", "E32" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dgeMatrix\"\n", " [,1] [,2] [,3]\n", "[1,] 2 -4 2\n", "[2,] 0 -1 3\n", "[3,] 0 0 16" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# zero (3, 2) entry\n", "E32 %*% E31 %*% E21 %*% A # Step 3, A'''" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Thus we have arrived at an upper triangular matrix $\\mathbf{U}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gauss transformations\n", "\n", "* We can combine elementary operations by column.\n", "\n", "* For column 1,\n", "$$\n", " \\mathbf{M}_1 = \\mathbf{E}_{n1}(c_{n1}) \\cdots \\mathbf{E}_{31}(c_{31}) \\mathbf{E}_{21}(c_{21}) = \\begin{bmatrix}\n", " 1 & \\\\\n", " c_{21} & \\\\\n", " & \\ddots & \\\\\n", " c_{n1} & & 1\n", " \\end{bmatrix}\n", "$$ \n", "\n", "* For column $k$,\n", "$$\n", "\t\\mathbf{M}_k = \\mathbf{E}_{nk}(c_{nk}) \\cdots \\mathbf{E}_{k+1,k}(c_{k+1,k}) = \\begin{bmatrix}\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{bmatrix}.\n", "$$\n", "\n", "* Matrices $\\mathbf{M}_1, \\ldots, \\mathbf{M}_{n-1}$ are called the **Gauss transformations**." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 sparse Matrix of class \"dtCMatrix\"\n", " \n", "[1,] 1 . .\n", "[2,] -2 1 .\n", "[3,] -1 . 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "M1 <- E31 %*% E21 # inv(L1 * L2)\n", "M1" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 sparse Matrix of class \"dtCMatrix\"\n", " \n", "[1,] 1 . .\n", "[2,] . 1 .\n", "[3,] . 5 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "M2 <- E32 # inv(L3)\n", "M2" ] }, { "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{bmatrix}\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{bmatrix}.\n", "$$" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 sparse Matrix of class \"dtCMatrix\"\n", " \n", "[1,] 1 . .\n", "[2,] 2 1 .\n", "[3,] 1 . 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "solve(M1) # L2; solve(A) gives the inverse of A, but use it with caution (see below)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 sparse Matrix of class \"dtCMatrix\"\n", " \n", "[1,] 1 . .\n", "[2,] . 1 .\n", "[3,] . -5 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "solve(M2) # L3" ] }, { "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{bmatrix} \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{bmatrix}.\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", "* 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", " - This is basically how R (and MATLAB and Julia) compute determinants." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* The whole LU algorithm is done *in place*, i.e., $\\mathbf{A}$ is **overwritten** by $\\mathbf{L}$ and $\\mathbf{U}$.\n", "```r\n", "for (k in seq_len(n-1)) {\n", " for (i in k +seq_len(max(n - k, 0))) {\n", " A[i, k] <- A[i, k] / A[k, k] # computes L\n", " for (j in k +seq_len(max(n - k, 0))) {\n", " A[i, j] <- A[i, j] - A[i, k] * A[k, j] # computes U\n", " }\n", " }\n", "}\n", "```\n", "\n", "\n", "\n", "Source: \n", "\n", "* The LU decomposition costs (ignoring the cost of division)\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", "* 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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 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", " > In short, no ```solve(A) %*% b```, but ```solve(A, b)```.\n", "\n", " We do not compute matrix inverse unless \n", " 1. it is necessary to compute standard errors \n", " 2. number of right hand sides is much larger than $n$ \n", " 3. $n$ is small\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Pivoting \n", "\n", "* Let\n", "$$\n", " \\mathbf{A} = \\begin{bmatrix}\n", " 0 & 1 \\\\\n", " 1 & 0 \\\\\n", " \\end{bmatrix}.\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", " **Example**:\n", "$$\n", "\\begin{split}\n", "0.001x_1 + x_2 &= 1 \\\\\n", " x_1 + x_2 &= 2 \n", "\\end{split}\n", "$$\n", "with solution\n", "\n", "$$\n", "\\begin{bmatrix} x_1 \\\\ x_ 2 \\end{bmatrix}\n", "=\n", "\\begin{bmatrix} 1.001001 \\\\ 0.998999 \\end{bmatrix}\n", "\\approx\n", "\\begin{bmatrix} 1.0 \\\\ 1.0 \\end{bmatrix}\n", ".\n", "$$\n", "\n", "* With two significant digits, GE via LU yields \n", "$$\n", "\\begin{split}\n", "0.001x_1 + x_2 &= 1 \\\\\n", " (1-1000)x_2 &= 2 - 1000 \n", "\\end{split}\n", "$$\n", "or \n", "$$\n", "\\begin{split}\n", "0.001x_1 + x_2 &= 1 \\\\\n", " -1000 x_2 &= -1000 \n", "\\end{split}\n", ",\n", "$$\n", "yielding\n", "$$\n", "\\begin{bmatrix} x_1 \\\\ x_ 2 \\end{bmatrix}\n", "=\n", "\\begin{bmatrix} 0.0 \\\\ 1.0 \\end{bmatrix}\n", "!\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Condition number and stability\n", "\n", "- The problem of solving $\\mathbf{A}\\mathbf{x}=\\mathbf{b}$ for fixed $b$ can be viewed as a map $f: \\mathbb{R}^{n\\times n} \\ni \\mathbf{A} \\mapsto \\mathbf{A}^{-1}\\mathbf{b} \\in \\mathbb{R}^n$.\n", "\n", "- A *well-conditioned* problem is one such that *all* small perturbations of a given input $\\mathbf{A}$ lead to only small changes in $f(\\mathbf{A})$.\n", " + An *ill-conditioned* problem is one such that *some* small perturbation of $\\mathbf{A}$ leads to a large change in $f(\\mathbf{A})$.\n", "\n", "- The conditioning of the problem of computing $\\mathbf{A}^{-1}\\mathbf{b}$ given $\\mathbf{A}$ is usually measured by the **condition number** \n", "$$\n", "\\kappa(\\mathbf{A}) = \\sigma_{\\max}(\\mathbf{A})/\\sigma_{\\min}(\\mathbf{A})\n", "$$\n", " + $\\sigma_{\\max}(\\mathbf{A}) = \\max_{\\mathbf{x}\\neq 0} \\|\\mathbf{A}\\mathbf{x}\\|/\\|\\mathbf{x}\\|$ is the maximum singular value of $\\mathbf{A}$;\n", " + $\\sigma_{\\min}(\\mathbf{A}) = \\min_{\\mathbf{x}\\neq 0} \\|\\mathbf{A}\\mathbf{x}\\|/\\|\\mathbf{x}\\|$ is the minimum singular value of $\\mathbf{A}$.\n", " \n", "- If $\\mathbf{A}$ is singular, then $\\sigma_{\\min}(\\mathbf{A}) = 0$, hence $\\kappa(\\mathbf{A}) = \\infty$. So a large condition number indicates that $\\mathbf{A}$ is close to singular and the problem is inherently hard." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Stability\n", "\n", "- An *algorithm* for computing $\\mathbf{A}^{-1}\\mathbf{b}$ using floating-point arithmetric can be viewed as another map $\\tilde{f}: \\mathbb{R}^{n\\times n} \\to \\mathbb{R}^{n}$.\n", "\n", "- For given input $\\mathbf{A}$, the solution computed by algorithm $\\tilde{f}$ is $\\hat{\\mathbf{x}} = \\tilde{f}(\\mathbf{A})$. \n", " + Example 1: solve $\\mathbf{A}\\mathbf{x}=\\mathbf{b}$ by GE/LU followed by forward and backward substitutions.\n", " + Example 2: solve $\\mathbf{A}\\mathbf{x}=\\mathbf{b}$ by Gauss-Seidel (an iterative method to come).\n", " \n", "- Algorithms will be affected by at least rounding errors.\n", "\n", "- Algorithm $\\tilde{f}$ is called *stable* if for any perturbed input $\\tilde{\\mathbf{A}}$ from $\\mathbf{A}$ with relative error bounded by $\\epsilon$, the relative error between $\\hat{\\mathbf{x}}$ and $\\tilde{\\mathbf{A}}^{-1}\\mathbf{b}$ is bounded by $O(\\epsilon)$.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Analysis of the example\n", "\n", "- The condition number of matrix $\\mathbf{A} = \\begin{bmatrix} 0.001 & 1 \\\\ 1 & 1 \\end{bmatrix}$ is 2.262, so the problem itself is well-conditioned: a small change in $\\mathbf{A}$ won't yield a large change in the solution $\\mathbf{x} = \\mathbf{A}^{-1}\\mathbf{b}$. \n", "\n", "- It is the instability of the GE/LU algorithm that yields the large error between the true solution $\\mathbf{x}\\approx(1, 1)$ and the computed solution (or the output of the algorithm) $\\hat{\\mathbf{x}}=(0, 1)$. This example suggests that with GE/LU, the worst-case relative error is $\\|\\mathbf{x}-\\hat{\\mathbf{x}}\\| / \\|\\mathbf{x}\\| \\ge 1/\\sqrt{2}$ (zero relative error in input $\\tilde{\\mathbf{A}}$), not small.\n", "\n", "- The root cause is the vastly different scale of $a_{11}$ with other coefficients and the order of ellimination. In particular, any $a_{22}$ such that $[a_{22} - 1000] = -1000$ and $b_2$ such that $[b_2 - 1000] = 1000$ will yield the same solution $\\hat{\\mathbf{x}}=(0, 1)$. \n", "That is, the $1000 = a_{21}/a_{11}$ causes the loss of information in $a_{22}$ and $b_2$.\n", "(Recall Scenario 1 of Catastrophic Cancellation in Lecture 1.)\n", "\n", "- Indeed, the LU decomposition of $\\mathbf{A}$ within the precision is\n", "$$\n", " \\mathbf{L} = \\begin{bmatrix} 1 & 0 \\\\ 1000 & 1 \\end{bmatrix},\n", " \\quad\n", " \\mathbf{U} = \\begin{bmatrix} 0.001 & 1 \\\\ 0 & -1000 \\end{bmatrix},\n", "$$\n", "resulting in $\\mathbf{L}\\mathbf{U} = \\begin{bmatrix} 0.001 & 1 \\\\ 1 & 0 \\end{bmatrix}$, so $\\mathbf{L}\\mathbf{U}$ is not even closed to $\\mathbf{A} = \\begin{bmatrix} 0.001 & 1 \\\\ 1 & 1 \\end{bmatrix}$!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Pivoting (cont'd)\n", "\n", "* Solution: fix the algorithm by **pivoting** (i.e., change the order of elimination)\n", "$$\n", "\\begin{split}\n", " x_1 + x_2 &= 2 \\\\\n", "0.001x_1 + x_2 &= 1 \n", "\\end{split}\n", "$$\n", "\n", "* With two significant digits, GE yields \n", "$$\n", "\\begin{split}\n", " x_1 + x_2 &= 2 \\\\\n", " (1-0.001)x_2 &= 1 - 0.002 \n", "\\end{split}\n", "$$\n", "or\n", "$$\n", "\\begin{split}\n", " x_1 + x_2 &= 2 \\\\\n", " 1.0 x_2 &= 1.0\n", "\\end{split}\n", ",\n", "$$\n", "yielding\n", "$$\n", "\\begin{bmatrix} x_1 \\\\ x_ 2 \\end{bmatrix}\n", "=\n", "\\begin{bmatrix} 1.0 \\\\ 1.0 \\end{bmatrix}\n", ".\n", "$$\n", "\n", "- Note that $0.001=a_{11}/a_{22}$ does not cause a loss of information in $a_{12}$ and $b_1$. The LU decomposition within the precision is\n", "$$\n", " \\mathbf{L} = \\begin{bmatrix} 1 & 0 \\\\ 0.001 & 1 \\end{bmatrix},\n", " \\quad\n", " \\mathbf{U} = \\begin{bmatrix} 1 & 1 \\\\ 0 & 1 \\end{bmatrix},\n", " \\quad\n", " \\mathbf{L}\\mathbf{U} = \\begin{bmatrix} 1 & 1 \\\\ 0.001 & 1.0 \\end{bmatrix}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Partial pivoting**: before zeroing the $k$-th column, move the row with $\\max_{i=k}^n |a_{ik}|$ is into the $k$-th row (called GEPP).\n", " - Cost: $(n-1) + (n-2) + \\dotsb + 1 = O(n^2)$.\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 (GECP): 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", "* GEPP is the most commonly used methods for solving **general** linear systems. GECP is more stable but costs more computation. Partial pivoting is stable most of times (which is partially because GEPP guarantees $|\\ell_{ij}| \\le 1$)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementation\n", "\n", "* Except for the iterative methods, most of these numerical linear algebra tasks are implemented in the BLAS and LAPACK libraries. LAPACK is built on top of BLAS. They form the **building blocks** of most statistical computing tasks.\n", "\n", "> LAPACK (Linear Algebra PACKage) is a standard software library for numerical linear algebra. It provides routines for solving systems of linear equations and linear least squares, eigenvalue problems, and singular value decomposition. -[Wikipedia](https://en.wikipedia.org/wiki/LAPACK)\n", "\n", "* LAPACK routine [`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`." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dgeMatrix\"\n", " [,1] [,2] [,3]\n", "[1,] 2 -4 2\n", "[2,] 4 -9 7\n", "[3,] 2 1 3" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "A # toy example, just to recall" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/html": [ "'denseLU'" ], "text/latex": [ "'denseLU'" ], "text/markdown": [ "'denseLU'" ], "text/plain": [ "[1] \"denseLU\"\n", "attr(,\"package\")\n", "[1] \"Matrix\"" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "alu <- Matrix::lu(A)\n", "class(alu)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [], "source": [ "ealu <- expand(alu) # expand the decomposition into Factors" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dtrMatrix\" (unitriangular)\n", " [,1] [,2] [,3] \n", "[1,] 1.00000000 . .\n", "[2,] 0.50000000 1.00000000 .\n", "[3,] 0.50000000 0.09090909 1.00000000" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ealu$L" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dtrMatrix\"\n", " [,1] [,2] [,3] \n", "[1,] 4.000000 -9.000000 7.000000\n", "[2,] . 5.500000 -0.500000\n", "[3,] . . -1.454545" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ealu$U" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 sparse Matrix of class \"pMatrix\"\n", " \n", "[1,] . . |\n", "[2,] | . .\n", "[3,] . | ." ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ealu$P" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dgeMatrix\"\n", " [,1] [,2] [,3]\n", "[1,] 4 -9 7\n", "[2,] 2 1 3\n", "[3,] 2 -4 2" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "with(ealu, L %*% U)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dgeMatrix\"\n", " [,1] [,2] [,3]\n", "[1,] 2 -4 2\n", "[2,] 4 -9 7\n", "[3,] 2 1 3" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "with(ealu, P %*% L %*% U)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/html": [ "-32" ], "text/latex": [ "-32" ], "text/markdown": [ "-32" ], "text/plain": [ "[1] -32" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "det(A) # this does LU!" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/html": [ "-32" ], "text/latex": [ "-32" ], "text/markdown": [ "-32" ], "text/plain": [ "[1] -32" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "det(ealu$U) # this is cheap" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dgeMatrix\"\n", " [,1] [,2] [,3]\n", "[1,] 1.0625 -0.4375 0.3125\n", "[2,] -0.0625 -0.0625 0.1875\n", "[3,] -0.6875 0.3125 0.0625" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "solve(A) # this does LU!" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dgeMatrix\"\n", " [,1] [,2] [,3]\n", "[1,] 1.0625 -0.4375 0.3125\n", "[2,] -0.0625 -0.0625 0.1875\n", "[3,] -0.6875 0.3125 0.0625" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "with(ealu, solve(U) %*% solve(L) %*% t(P) ) # this is cheap???" ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "R", "language": "R", "name": "ir" }, "language_info": { "codemirror_mode": "r", "file_extension": ".r", "mimetype": "text/x-r-source", "name": "R", "pygments_lexer": "r", "version": "4.2.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, "skip_h1_title": true, "threshold": 4, "toc_cell": true, "toc_section_display": "block", "toc_window_display": true, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 4 }