{ "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)\n", "library(mvtnorm)\n", "library(microbenchmark)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Cholesky Decomposition\n", "\n", "* Named after [André-Louis Cholesky, 1875-1918](https://en.wikipedia.org/wiki/André-Louis_Cholesky).\n", "\n", "* Our mantra 2:\n", "\n", "> **The structure should be exploited whenever possible in solving a problem.** \n", "\n", " Common structures include: symmetry, positive (semi)definiteness, sparsity, Kronecker product, low rank, ...\n", "\n", "* Recall that a matrix $M$ is positive (semi)definite if\n", "$$\n", " \\mathbf{x}^T\\mathbf{M}\\mathbf{x} > 0 ~(\\ge 0), \\quad \\forall \\mathbf{x}\\neq\\mathbf{0}.\n", "$$\n", "\n", "* We can always assume that a postive (semi)definite matrix is symmetric (why?)\n", "\n", "* Example: normal equation \n", "$$\n", " \\mathbf{X}^T \\mathbf{X} \\beta = \\mathbf{X}^T \\mathbf{y}\n", "$$\n", "for linear regression, the coefficient matrix $\\mathbf{X}^T \\mathbf{X}$ is symmetric and positive semidefinite. How to exploit this structure?\n", "\n", "* Most of time statisticians deal with positive (semi)definite matrices." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cholesky decomposition\n", "\n", "* **Theorem**: Let $\\mathbf{A} \\in \\mathbb{R}^{n \\times n}$ be symmetric and positive definite. Then $\\mathbf{A} = \\mathbf{L} \\mathbf{L}^T$, where $\\mathbf{L}$ is lower triangular with positive diagonal entries and is unique. \n", "**Proof** (by induction): \n", "If $n=1$, then $0 < a = \\sqrt{a}\\sqrt{a}$. For $n>1$, the block equation\n", "$$\n", "\\begin{eqnarray*}\n", "\\begin{bmatrix}\n", "a_{11} & \\mathbf{a}^T \\\\ \\mathbf{a} & \\mathbf{A}_{22}\n", "\\end{bmatrix} =\n", "\\begin{bmatrix}\n", "\t\\ell_{11} & \\mathbf{0}_{n-1}^T \\\\ \\mathbf{b} & \\mathbf{L}_{22}\n", "\\end{bmatrix}\n", "\\begin{bmatrix}\n", "\t\\ell_{11} & \\mathbf{b}^T \\\\ \\mathbf{0}_{n-1} & \\mathbf{L}_{22}^T\n", "\\end{bmatrix}\n", "\\end{eqnarray*}\n", "$$\n", "entails\n", "$$\n", "\\begin{eqnarray*}\n", "\ta_{11} &=& \\ell_{11}^2 \\\\\n", "\t\\mathbf{a} &=& \\ell_{11} \\mathbf{b}\t\\\\\n", "\t\\mathbf{A}_{22} &=& \\mathbf{b} \\mathbf{b}^T + \\mathbf{L}_{22} \\mathbf{L}_{22}^T\n", " .\n", "\\end{eqnarray*}\n", "$$ \n", "Since $a_{11}>0$ (why?), $\\ell_{11}=\\sqrt{a_{11}}$ and $\\mathbf{b}=a_{11}^{-1/2}\\mathbf{a}$ are uniquely determined. \n", "$\\mathbf{A}_{22} - \\mathbf{b} \\mathbf{b}^T = \\mathbf{A}_{22} - a_{11}^{-1} \\mathbf{a} \\mathbf{a}^T$ is positive definite of size $(n-1)\\times(n-1)$ because $\\mathbf{A}_{22}$ is positive definite (why? homework). By the induction hypothesis, lower triangular $\\mathbf{L}_{22}$ such that $\\mathbf{A}_{22} - \\mathbf{b} \\mathbf{b}^T = \\mathbf{L}_{22}^T\\mathbf{L}_{22}$ exists and is unique." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* This proof is constructive and completely specifies the algorithm: \n", "```r\n", "for (k in seq_len(n)) {\n", " for (j in k + seq_len(n - k)) {\n", " a_jk_div_a_kk <- A[j, k] / A[k, k] \n", " for (i in j - 1 + seq_len(n - j + 1)) {\n", " A[i, j] <- A[i, j] - A[i, k] * a_jk_div_a_kk # L_22\n", " }\n", " }\n", " sqrt_akk <- sqrt(A[k, k])\n", " for (i in k - 1 + seq_len(n - k + 1)) {\n", " A[i, k] = A[i, k] / sqrt_akk # l_11 and b\n", " }\n", "}\n", "```\n", "\n", "\n", "\n", "Source: \n", "\n", "* Computational cost: \n", "$$\n", "\\frac{1}{2} [2(n-1)^2 + 2(n-2)^2 + \\cdots + 2 \\cdot 1^2] \\approx \\frac{1}{3} n^3 \\quad \\text{flops}\n", "$$ \n", "plus $n$ square roots. **Half** the cost of LU decomposition.\n", "\n", "* In general Cholesky decomposition is very stable. Failure of the decomposition simply means $\\mathbf{A}$ is not positive definite. **It is an efficient way to test positive definiteness.**\n", "\n", "\n", "## Pivoting\n", "\n", "* Pivoting is only needed if you want the Cholesky factor of a positive *semidefinite* matrix.\n", "\n", "* When $\\mathbf{A}$ does not have full rank, e.g., $\\mathbf{X}^T \\mathbf{X}$ with a non-full column rank $\\mathbf{X}$, we encounter $a_{kk} = 0$ during the procedure.\n", "\n", "* **Symmetric pivoting**. At each stage $k$, we permute both row and column such that $\\max_{k \\le i \\le n} a_{ii}$ becomes the pivot. If we encounter $\\max_{k \\le i \\le n} a_{ii} = 0$, then $\\mathbf{A}[k:n,k:n] = \\mathbf{0}$ (why?) and the algorithm terminates.\n", "\n", "* With symmetric pivoting: \n", "$$\n", "\\mathbf{P} \\mathbf{A} \\mathbf{P}^T = \\mathbf{L} \\mathbf{L}^T,\n", "$$\n", "where $\\mathbf{P}$ is a permutation matrix and $\\mathbf{L} \\in \\mathbb{R}^{n \\times r}$, $r = \\text{rank}(\\mathbf{A})$.\n", "\n", "## Implementation\n", "\n", "* LAPACK functions: [`dpotrf`](http://www.netlib.org/lapack/explore-html/d1/d7a/group__double_p_ocomputational_ga2f55f604a6003d03b5cd4a0adcfb74d6.html#ga2f55f604a6003d03b5cd4a0adcfb74d6) (without pivoting), [`dpstrf`](http://www.netlib.org/lapack/explore-html/da/dba/group__double_o_t_h_e_rcomputational_ga31cdc13a7f4ad687f4aefebff870e1cc.html#ga31cdc13a7f4ad687f4aefebff870e1cc) (with pivoting).\n", "\n", "* R functions: `base::chol()` or `Matrix::chol()` (gives an **upper** triangular factor: $A = R^T R$)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example: positive definite matrix" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dsyMatrix\"\n", " [,1] [,2] [,3]\n", "[1,] 4 12 -16\n", "[2,] 12 37 -43\n", "[3,] -16 -43 98" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "A <- t(Matrix(c(4, 12, -16, 12, 37, -43, -16, -43, 98), nrow=3)) # check out this is pd!\n", "A" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"Cholesky\"\n", " [,1] [,2] [,3]\n", "[1,] 2 6 -8\n", "[2,] . 1 5\n", "[3,] . . 3" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Cholesky without pivoting\n", "Achol <- Matrix::chol(Matrix::symmpart(A))\n", "Achol" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/html": [ "'Cholesky'" ], "text/latex": [ "'Cholesky'" ], "text/markdown": [ "'Cholesky'" ], "text/plain": [ "[1] \"Cholesky\"\n", "attr(,\"package\")\n", "[1] \"Matrix\"" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "class(Achol)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/html": [ "
Dim
'integer'
Dimnames
'list'
x
'numeric'
uplo
'character'
diag
'character'
\n" ], "text/latex": [ "\\begin{description*}\n", "\\item[Dim] 'integer'\n", "\\item[Dimnames] 'list'\n", "\\item[x] 'numeric'\n", "\\item[uplo] 'character'\n", "\\item[diag] 'character'\n", "\\end{description*}\n" ], "text/markdown": [ "Dim\n", ": 'integer'Dimnames\n", ": 'list'x\n", ": 'numeric'uplo\n", ": 'character'diag\n", ": 'character'\n", "\n" ], "text/plain": [ " Dim Dimnames x uplo diag \n", " \"integer\" \"list\" \"numeric\" \"character\" \"character\" " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "getSlots(\"Cholesky\")" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Formal class 'Cholesky' [package \"Matrix\"] with 5 slots\n", " ..@ Dim : int [1:2] 3 3\n", " ..@ Dimnames:List of 2\n", " .. ..$ : NULL\n", " .. ..$ : NULL\n", " ..@ x : num [1:9] 2 0 0 6 1 0 -8 5 3\n", " ..@ uplo : chr \"U\"\n", " ..@ diag : chr \"N\"\n" ] } ], "source": [ "str(Achol)" ] }, { "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,] 28.583333\n", "[2,] -7.666667\n", "[3,] 1.333333" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "b <- c(1.0, 2.0, 3.0)\n", "solve(A, b) # this does LU; wasteful!; 2/3 n^3 + 2n^2" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 1 Matrix of class \"dgeMatrix\"\n", " [,1]\n", "[1,] 28.583333\n", "[2,] -7.666667\n", "[3,] 1.333333" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "y <- solve(t(Achol), b) # two triangular solves; only 2n^2 flops\n", "solve(Achol, y)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/html": [ "36" ], "text/latex": [ "36" ], "text/markdown": [ "36" ], "text/plain": [ "[1] 36" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "det(A) # this actually does LU; wasteful!" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/html": [ "36" ], "text/latex": [ "36" ], "text/markdown": [ "36" ], "text/plain": [ "[1] 36" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "det(Achol)^2 # cheap" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dsyMatrix\"\n", " [,1] [,2] [,3]\n", "[1,] 49.361111 -13.5555556 2.1111111\n", "[2,] -13.555556 3.7777778 -0.5555556\n", "[3,] 2.111111 -0.5555556 0.1111111" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "solve(A) # this does LU!" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dgeMatrix\"\n", " [,1] [,2] [,3]\n", "[1,] 49.361111 -13.5555556 2.1111111\n", "[2,] -13.555556 3.7777778 -0.5555556\n", "[3,] 2.111111 -0.5555556 0.1111111" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "solve(Achol, solve(t(Achol))) # 2n triangular solves" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example: positive *semi*definite matrix" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "3 x 3 Matrix of class \"dsyMatrix\"\n", " [,1] [,2] [,3]\n", "[1,] 1.0 1.0 0.5\n", "[2,] 1.0 1.0 0.5\n", "[3,] 0.5 0.5 1.0" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
  1. 2.36602540378444
  2. 0.633974596215564
  3. 0
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 2.36602540378444\n", "\\item 0.633974596215564\n", "\\item 0\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 2.36602540378444\n", "2. 0.633974596215564\n", "3. 0\n", "\n", "\n" ], "text/plain": [ "[1] 2.3660254 0.6339746 0.0000000" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "(A <- Matrix(c(1,1,.5,1,1,.5,.5,.5,1), nrow=3))\n", "eigen(A, symmetric = TRUE)$values" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "ename": "ERROR", "evalue": "Error in value[[3L]](cond): chol(x) is undefined: 'x' is not positive definite\n", "output_type": "error", "traceback": [ "Error in value[[3L]](cond): chol(x) is undefined: 'x' is not positive definite\nTraceback:\n", "1. Matrix::chol(A, pivot = FALSE)", "2. Matrix::chol(A, pivot = FALSE)", "3. tryCatch(.Call(dpoMatrix_trf, if (x@uplo == \"U\") x else t(x), \n . 2L), error = function(e) stop(\"chol(x) is undefined: 'x' is not positive definite\"))", "4. tryCatchList(expr, classes, parentenv, handlers)", "5. tryCatchOne(expr, names, parentenv, handlers[[1L]])", "6. value[[3L]](cond)", "7. stop(\"chol(x) is undefined: 'x' is not positive definite\")" ] } ], "source": [ "Matrix::chol(A, pivot=FALSE) # 2nd argument requests pivoting" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Warning message in chol.default(A, pivot = TRUE):\n", "“the matrix is either rank-deficient or indefinite”\n" ] }, { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\n", "
A matrix: 3 × 3 of type dbl
10.50000001
00.86602540
00.00000000
\n" ], "text/latex": [ "A matrix: 3 × 3 of type dbl\n", "\\begin{tabular}{lll}\n", "\t 1 & 0.5000000 & 1\\\\\n", "\t 0 & 0.8660254 & 0\\\\\n", "\t 0 & 0.0000000 & 0\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "A matrix: 3 × 3 of type dbl\n", "\n", "| 1 | 0.5000000 | 1 |\n", "| 0 | 0.8660254 | 0 |\n", "| 0 | 0.0000000 | 0 |\n", "\n" ], "text/plain": [ " [,1] [,2] [,3]\n", "[1,] 1 0.5000000 1 \n", "[2,] 0 0.8660254 0 \n", "[3,] 0 0.0000000 0 " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
  1. 1
  2. 3
  3. 2
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 1\n", "\\item 3\n", "\\item 2\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 1\n", "2. 3\n", "3. 2\n", "\n", "\n" ], "text/plain": [ "[1] 1 3 2" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "2" ], "text/latex": [ "2" ], "text/markdown": [ "2" ], "text/plain": [ "[1] 2" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "(Achol <- base::chol(A, pivot=TRUE)) # turn off checking pd\n", "attr(Achol, \"pivot\")\n", "attr(Achol, \"rank\")" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/html": [ "2" ], "text/latex": [ "2" ], "text/markdown": [ "2" ], "text/plain": [ "[1] 2\n", "attr(,\"method\")\n", "[1] \"tolNorm2\"\n", "attr(,\"useGrad\")\n", "[1] FALSE\n", "attr(,\"tol\")\n", "[1] 6.661338e-16" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Matrix::rankMatrix(A) # determine rank from SVD, which is more numerically stable" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\n", "
A matrix: 3 × 3 of type dbl
110.5000000
000.8660254
000.0000000
\n" ], "text/latex": [ "A matrix: 3 × 3 of type dbl\n", "\\begin{tabular}{lll}\n", "\t 1 & 1 & 0.5000000\\\\\n", "\t 0 & 0 & 0.8660254\\\\\n", "\t 0 & 0 & 0.0000000\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "A matrix: 3 × 3 of type dbl\n", "\n", "| 1 | 1 | 0.5000000 |\n", "| 0 | 0 | 0.8660254 |\n", "| 0 | 0 | 0.0000000 |\n", "\n" ], "text/plain": [ " [,1] [,2] [,3] \n", "[1,] 1 1 0.5000000\n", "[2,] 0 0 0.8660254\n", "[3,] 0 0 0.0000000" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Achol[, attr(Achol, \"pivot\")]" ] }, { "cell_type": "code", "execution_count": 19, "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": [ "(P <- as(attr(Achol, \"pivot\"), \"pMatrix\")) # permutation matrix, actually P' in our notation" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/html": [ "1.11022302462516e-16" ], "text/latex": [ "1.11022302462516e-16" ], "text/markdown": [ "1.11022302462516e-16" ], "text/plain": [ "[1] 1.110223e-16" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# P A P' = L L'\n", "Matrix::norm(t(P) %*% A %*% P - t(Achol) %*% Achol, type=\"F\") # Frobenius norm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Applications\n", "\n", "* **No inversion** principle: whenever we see matrix inverse, we should think in terms of solving linear equations. If the matrix is positive (semi)definite, use Cholesky decomposition, which is twice cheaper than LU decomposition.\n", "\n", "### Multivariate normal density \n", "\n", "Multivariate normal density $\\mathcal{N}(0, \\Sigma)$, where $\\Sigma$ is $n\\times n$ positive definite, is\n", "$$\n", " \\frac{1}{(2\\pi)^{n/2}\\det(\\Sigma)^{1/2}}\\exp\\left(-\\frac{1}{2}\\mathbf{y}^T\\Sigma^{-1}\\mathbf{y}\\right).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hence the log likelihood is\n", "$$\n", "- \\frac{n}{2} \\log (2\\pi) - \\frac{1}{2} \\log \\det \\Sigma - \\frac{1}{2} \\mathbf{y}^T \\Sigma^{-1} \\mathbf{y}.\n", "$$\n", "\n", "* Method 1: \n", " 1. compute explicit inverse $\\Sigma^{-1}$ ($2n^3$ flops), \n", " 2. compute quadratic form ($2n^2 + 2n$ flops), \n", " 3. compute determinant ($2n^3/3$ flops).\n", " \n", "* Method 2: \n", " 1. Cholesky decomposition $\\Sigma = \\mathbf{L} \\mathbf{L}^T$ ($n^3/3$ flops), \n", " 2. Solve $\\mathbf{L} \\mathbf{x} = \\mathbf{y}$ by forward substitutions ($n^2$ flops), \n", " 3. compute quadratic form $\\mathbf{x}^T \\mathbf{x}$ ($2n$ flops), and \n", " 4. compute determinant from Cholesky factor ($n$ flops). \n", "\n", "**Which method is better?**" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [], "source": [ "# word-by-word transcription of mathematical expression\n", "# Matrix::determinant computes the logarithm of the determinant\n", "logpdf_mvn_1 <- function(y, Sigma) {\n", " n <- length(y)\n", " - (n / 2) * log(2 * pi) - (1 / 2) * Matrix::determinant(Sigma)$modulus - (1 / 2) * t(y) %*% solve(Sigma) %*% y\n", "}\n", "\n", "# you learnt why you should avoid inversion\n", "logpdf_mvn_2 <- function(y, Sigma) {\n", " n <- length(y)\n", " Sigmachol <- Matrix::chol(Matrix::symmpart(Sigma))\n", " - (n / 2) * log(2 * pi) - Matrix::determinant(Sigmachol)$modulus - (1 / 2) * sum(abs(solve(t(Sigmachol), y))^2)\n", "}\n", "\n", "# this is for the efficiency-concerned \n", "logpdf_mvn_3 <- function(y, Sigma) {\n", " n <- length(y)\n", " Sigmachol <- Matrix::chol(Matrix::symmpart(Sigma))\n", " - 0.5 * n * log(2 * pi) - sum(log(diag(Sigmachol))) - 0.5 * sum(abs(solve(t(Sigmachol), y))^2)\n", "}" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [ { "data": { "text/plain": [ "1 x 1 Matrix of class \"dgeMatrix\"\n", " [,1]\n", "[1,] -5272.565" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "-5272.56471529637" ], "text/latex": [ "-5272.56471529637" ], "text/markdown": [ "-5272.56471529637" ], "text/plain": [ "[1] -5272.565\n", "attr(,\"logarithm\")\n", "[1] TRUE" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "-5272.56471529637" ], "text/latex": [ "-5272.56471529637" ], "text/markdown": [ "-5272.56471529637" ], "text/plain": [ "[1] -5272.565" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "set.seed(123) # seed\n", "\n", "n <- 1000\n", "# a pd matrix\n", "Sigma_half <- Matrix(rnorm(n * (2 * n)), nrow=n)\n", "Sigma <- Sigma_half %*% t(Sigma_half)\n", "Sigma <- Matrix::symmpart(Sigma)\n", "Sigmachol <- Matrix::chol(Sigma)\n", "\n", "y <- Sigmachol %*% rnorm(n) # one random sample from N(0, Sigma). Homework\n", "\n", "# at least they give same answer\n", "(logpdf_mvn_1(y, Sigma))\n", "(logpdf_mvn_2(y, Sigma))\n", "(logpdf_mvn_3(y, Sigma))" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Unit: milliseconds\n", " expr min lq mean median uq\n", " logpdf_mvn_1(y, Sigma) 194.263957 197.941488 202.391593 202.178377 205.562591\n", " logpdf_mvn_2(y, Sigma) 2.657588 3.105828 5.582544 3.476632 4.333807\n", " logpdf_mvn_3(y, Sigma) 2.536665 3.143556 5.134440 3.367640 3.741760\n", " max neval cld\n", " 223.3571 100 b\n", " 111.4274 100 a \n", " 111.8951 100 a \n" ] } ], "source": [ "res <- microbenchmark::microbenchmark(logpdf_mvn_1(y, Sigma), logpdf_mvn_2(y, Sigma), logpdf_mvn_3(y, Sigma))\n", "print(res)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* To evaluate same multivariate normal density at many observations $y_1, y_2, \\ldots$, we pre-compute the Cholesky decomposition ($n^3/3$ flops), then each evaluation costs $n^2$ flops." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Linear regression by Cholesky\n", "\n", "* Cholesky decomposition is **one** approach to solve linear regression. Assume $\\mathbf{X} \\in \\mathbb{R}^{n \\times p}$ and $\\mathbf{y} \\in \\mathbb{R}^n$. \n", " - Compute $\\mathbf{X}^T \\mathbf{X}$: $np^2$ flops \n", " - Compute $\\mathbf{X}^T \\mathbf{y}$: $2np$ flops \n", " - Cholesky decomposition of $\\mathbf{X}^T \\mathbf{X}$: $\\frac{1}{3} p^3$ flops \n", " - Solve normal equation $\\mathbf{X}^T \\mathbf{X} \\beta = \\mathbf{X}^T \\mathbf{y}$: $2p^2$ flops \n", " - If need standard errors, another $(4/3)p^3$ flops\n", "\n", "Total computational cost is $np^2 + (1/3) p^3$ (without s.e.) or $np^2 + (5/3) p^3$ (with s.e.) flops." ] } ], "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": "118px", "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 }