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

1  Cholesky Decomposition
1.1  Cholesky decomposition
1.2  Pivoting
1.3  Implementation
1.3.1  Example: positive definite matrix.
1.3.2  Example: positive semi-definite matrix.
1.4  Applications
1.4.1  Multivariate normal density
1.4.2  Linear regression
1.5  Further reading
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Cholesky Decomposition\n", "\n", "\n", "\n", "* A basic tenet in numerical analysis: \n", "\n", "> **The structure should be exploited whenever solving a problem.** \n", "\n", " Common structures include: symmetry, definiteness, sparsity, Kronecker product, low rank, ...\n", "\n", "* LU decomposition (Gaussian Elimination) is **not** used in statistics so often because most of time statisticians deal with positive (semi)definite matrix. (That's why I hate to see `solve()` in R code.)\n", "\n", "* Consider solving the 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", "## 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 $\\ell = \\sqrt{a}$. For $n>1$, the block equation\n", "$$\n", "\\begin{eqnarray*}\n", "\\begin{pmatrix}\n", "a_{11} & \\mathbf{a}^T \\\\ \\mathbf{a} & \\mathbf{A}_{22}\n", "\\end{pmatrix} =\n", "\\begin{pmatrix}\n", "\t\\ell_{11} & \\mathbf{0}_{n-1}^T \\\\ \\mathbf{l} & \\mathbf{L}_{22}\n", "\\end{pmatrix}\n", "\\begin{pmatrix}\n", "\t\\ell_{11} & \\mathbf{l}^T \\\\ \\mathbf{0}_{n-1} & \\mathbf{L}_{22}^T\n", "\\end{pmatrix}\n", "\\end{eqnarray*}\n", "$$\n", "has solution\n", "$$\n", "\\begin{eqnarray*}\n", "\t\\ell_{11} &=& \\sqrt{a_{11}} \\\\\n", "\t\\mathbf{l} &=& \\ell_{11}^{-1} \\mathbf{a}\t\\\\\n", "\t\\mathbf{L}_{22} \\mathbf{L}_{22}^T &=& \\mathbf{A}_{22} - \\mathbf{l} \\mathbf{l}^T = \\mathbf{A}_{22} - a_{11}^{-1} \\mathbf{a} \\mathbf{a}^T.\n", "\\end{eqnarray*}\n", "$$ \n", "Now $a_{11}>0$ (why?), so $\\ell_{11}$ and $\\mathbf{l}$ are uniquely determined. $\\mathbf{A}_{22} - a_{11}^{-1} \\mathbf{a} \\mathbf{a}^T$ is positive definite because $\\mathbf{A}$ is positive definite (why?). By induction hypothesis, $\\mathbf{L}_{22}$ exists and is unique.\n", "\n", "* The constructive proof completely specifies the algorithm: \n", "\n", "\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 by utilizing symmetry.\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", "* 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: [`?potrf`](http://www.netlib.org/lapack/explore-html/d1/d7a/group__double_p_ocomputational_ga2f55f604a6003d03b5cd4a0adcfb74d6.html#ga2f55f604a6003d03b5cd4a0adcfb74d6) (without pivoting), [`?pstrf`](http://www.netlib.org/lapack/explore-html/da/dba/group__double_o_t_h_e_rcomputational_ga31cdc13a7f4ad687f4aefebff870e1cc.html#ga31cdc13a7f4ad687f4aefebff870e1cc) (with pivoting).\n", "\n", "* Julia functions: [`cholfact`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.cholfact), [`cholfact!`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.cholfact!), [`chol`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.chol), or call LAPACK wrapper functions [`potrf!`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.LAPACK.potrf!) and [`pstrf!`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.LAPACK.pstrf!)\n", "\n", "### Example: positive definite matrix." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Int64,2}:\n", " 4 12 -16\n", " 12 37 -43\n", " -16 -43 98" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = [4 12 -16; 12 37 -43; -16 -43 98]" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:\n", "[2.0 6.0 -8.0; 0.0 1.0 5.0; 0.0 0.0 3.0]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Cholesky without pivoting\n", "Achol = cholfact(A)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 LowerTriangular{Float64,Array{Float64,2}}:\n", " 2.0 ⋅ ⋅ \n", " 6.0 1.0 ⋅ \n", " -8.0 5.0 3.0" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Achol[:L]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 UpperTriangular{Float64,Array{Float64,2}}:\n", " 2.0 6.0 -8.0\n", " ⋅ 1.0 5.0\n", " ⋅ ⋅ 3.0" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Achol[:U]" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 28.5833 \n", " -7.66667\n", " 1.33333" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b = [1.0; 2.0; 3.0]\n", "A \\ b # this does LU; wasteful!; 2/3 n^3 + 2n^2" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\\(A::AbstractArray{T,2} where T, B::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) at linalg/generic.jl:820" ], "text/plain": [ "\\(A::AbstractArray{T,2} where T, B::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) in Base.LinAlg at linalg/generic.jl:820" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@which A \\ b" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 28.5833 \n", " -7.66667\n", " 1.33333" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Achol \\ b # two triangular solves; only 2n^2 flops" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\\(F::Factorization, B::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) at linalg/factorization.jl:45" ], "text/plain": [ "\\(F::Factorization, B::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) in Base.LinAlg at linalg/factorization.jl:45" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@which Achol \\ b" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "36.0" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "det(A) # this actually does LU; wasteful!" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "det{T}(A::AbstractArray{T,2}) at linalg/generic.jl:1222" ], "text/plain": [ "det(A::AbstractArray{T,2}) where T in Base.LinAlg at linalg/generic.jl:1222" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@which det(A)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "36.0" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "det(Achol) # cheap" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "det(C::Base.LinAlg.Cholesky) at linalg/cholesky.jl:495" ], "text/plain": [ "det(C::Base.LinAlg.Cholesky) in Base.LinAlg at linalg/cholesky.jl:495" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@which det(Achol)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 49.3611 -13.5556 2.11111 \n", " -13.5556 3.77778 -0.555556\n", " 2.11111 -0.555556 0.111111" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inv(A) # this does LU!" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "inv{T}(A::Union{Base.ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T}, DenseArray{T,2}, SubArray{T,2,A,I,L} where L} where I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex, Int64, Range{Int64}},N} where N} where A<:Union{Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T} where N where T, DenseArray}) at linalg/dense.jl:651" ], "text/plain": [ "inv(A::Union{Base.ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T}, DenseArray{T,2}, SubArray{T,2,A,I,L} where L} where I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex, Int64, Range{Int64}},N} where N} where A<:Union{Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T} where N where T, DenseArray}) where T in Base.LinAlg at linalg/dense.jl:651" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@which inv(A)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 49.3611 -13.5556 2.11111 \n", " -13.5556 3.77778 -0.555556\n", " 2.11111 -0.555556 0.111111" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inv(Achol) " ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "inv(C::Base.LinAlg.Cholesky{#s268,#s267} where #s267<:(Union{Base.ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T}, DenseArray{T,2}, SubArray{T,2,A,I,L} where L} where I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex, Int64, Range{Int64}},N} where N} where A<:Union{Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T} where N where T, DenseArray} where T) where #s268<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64}) at linalg/cholesky.jl:537" ], "text/plain": [ "inv(C::Base.LinAlg.Cholesky{#s268,#s267} where #s267<:(Union{Base.ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T}, DenseArray{T,2}, SubArray{T,2,A,I,L} where L} where I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex, Int64, Range{Int64}},N} where N} where A<:Union{Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T} where N where T, DenseArray} where T) where #s268<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64}) in Base.LinAlg at linalg/cholesky.jl:537" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@which inv(Achol)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: positive semi-definite matrix." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Array{Float64,2}:\n", " 2.06704 -3.10361 2.41766 -0.717719 0.845433 \n", " -3.10361 9.76269 -7.31094 2.14335 0.283818 \n", " 2.41766 -7.31094 6.0473 -0.931321 -0.0610487\n", " -0.717719 2.14335 -0.931321 1.28376 0.115462 \n", " 0.845433 0.283818 -0.0610487 0.115462 0.827396 " ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "srand(280) # seed\n", "A = randn(5, 3)\n", "A = A * A'" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "ename": "LoadError", "evalue": "\u001b[91mBase.LinAlg.PosDefException(5)\u001b[39m", "output_type": "error", "traceback": [ "\u001b[91mBase.LinAlg.PosDefException(5)\u001b[39m", "", "Stacktrace:", " [1] \u001b[1m_chol!\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Array{Float64,2}, ::Type{UpperTriangular}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./linalg/cholesky.jl:55\u001b[22m\u001b[22m", " [2] \u001b[1mcholfact!\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Hermitian{Float64,Array{Float64,2}}, ::Type{Val{false}}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./linalg/cholesky.jl:211\u001b[22m\u001b[22m", " [3] \u001b[1mcholfact\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Array{Float64,2}, ::Symbol\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./linalg/cholesky.jl:350\u001b[22m\u001b[22m", " [4] \u001b[1mcholfact\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Array{Float64,2}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./linalg/cholesky.jl:349\u001b[22m\u001b[22m" ] } ], "source": [ "Achol = cholfact(A)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Base.LinAlg.CholeskyPivoted{Float64,Array{Float64,2}}([3.12453 -3.10361 … -0.717719 0.845433; -0.993306 1.03941 … 2.14335 0.283818; … ; -2.33985 0.0899256 … 5.96046e-8 0.115462; 0.0908354 0.900181 … -1.19908e-8 0.0], 'L', [2, 1, 4, 3, 5], 4, 0.0, 1)" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Achol = cholfact(A, :L, Val{true}) # 3rd argument request partial pivoting" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rank(Achol) # determine rank from Cholesky factor" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rank(A) # determine rank from SVD, which is more numerically stable" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 LowerTriangular{Float64,Array{Float64,2}}:\n", " 3.12453 ⋅ ⋅ ⋅ ⋅ \n", " -0.993306 1.03941 ⋅ ⋅ ⋅ \n", " 0.685974 -0.0349591 0.901097 ⋅ ⋅ \n", " -2.33985 0.0899256 0.751198 5.96046e-8 ⋅ \n", " 0.0908354 0.900181 0.0939082 -1.19908e-8 0.0" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Achol[:L]" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 UpperTriangular{Float64,Array{Float64,2}}:\n", " 3.12453 -0.993306 0.685974 -2.33985 0.0908354 \n", " ⋅ 1.03941 -0.0349591 0.0899256 0.900181 \n", " ⋅ ⋅ 0.901097 0.751198 0.0939082 \n", " ⋅ ⋅ ⋅ 5.96046e-8 -1.19908e-8\n", " ⋅ ⋅ ⋅ ⋅ 0.0 " ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Achol[:U]" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5-element Array{Int64,1}:\n", " 2\n", " 1\n", " 4\n", " 3\n", " 5" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Achol[:p]" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.7763568394002505e-15" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# P A P' = L U\n", "vecnorm(Achol[:P] * A * Achol[:P]' - Achol[:L] * Achol[:U])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Applications\n", "\n", "* **No inversion** mentality: 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 $\\text{MVN}(\\mu, \\Sigma)$, where $\\Sigma$ is p.d., is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "- \\frac{n}{2} \\log (2\\pi) - \\frac{1}{2} \\log \\det \\Sigma - \\frac{1}{2} (\\mathbf{y} - \\mu)^T \\Sigma^{-1} (\\mathbf{y} - \\mu).\n", "$$\n", "\n", "* Method 1: (a) compute explicit inverse $\\Sigma^{-1}$ ($2n^3$ flops), (b) compute quadratic form ($2n^2 + 2n$ flops), (c) compute determinant ($2n^3/3$ flops).\n", " \n", "* Method 2: (a) Cholesky decomposition $\\Sigma = \\mathbf{L} \\mathbf{L}^T$ ($n^3/3$ flops), (b) Solve $\\mathbf{L} \\mathbf{x} = \\mathbf{y} - \\mu$ by forward substitutions ($n^2$ flops), (c) compute quadratic form $\\mathbf{x}^T \\mathbf{x}$ ($2n$ flops), and (d) compute determinant from Cholesky factor ($n$ flops). \n", "\n", "**Which method is better?**" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "logpdf_mvn_3 (generic function with 1 method)" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# this is a person w/o any numerical analsyis training\n", "function logpdf_mvn_1(y::Vector, Σ::Matrix)\n", " n = length(y)\n", " - (n//2) * log(2π) - (1//2) * logdet(Σ) - (1//2) * y' * inv(Σ) * y\n", "end\n", "\n", "# this is an efficiency-savvy person \n", "function logpdf_mvn_2(y::Vector, Σ::Matrix)\n", " n = length(y)\n", " Σchol = cholfact(Symmetric(Σ))\n", " - (n//2) * log(2π) - (1//2) * logdet(Σchol) - (1//2) * sum(abs2, Σchol[:L] \\ y)\n", "end\n", "\n", "# better memory efficiency\n", "function logpdf_mvn_3(y::Vector, Σ::Matrix)\n", " n = length(y)\n", " Σchol = cholfact(Symmetric(Σ))\n", " - (n//2) * log(2π) - (1//2) * logdet(Σchol) - (1//2) * dot(y, Σchol \\ y)\n", "end" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "logpdf_mvn_1(y, Σ) = -4396.173284372774\n", "logpdf_mvn_2(y, Σ) = -4396.173284372773\n", "logpdf_mvn_3(y, Σ) = -4396.173284372773\n" ] } ], "source": [ "using BenchmarkTools, Distributions\n", "\n", "srand(280) # seed\n", "\n", "n = 1000\n", "Σ = randn(n, n); Σ = Σ' * Σ + I # covariance matrix\n", "y = rand(MvNormal(Σ)) # one randdom sample from N(0, Σ)\n", "\n", "# at least they give same answer\n", "@show logpdf_mvn_1(y, Σ)\n", "@show logpdf_mvn_2(y, Σ)\n", "@show logpdf_mvn_3(y, Σ);" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 23.42 MiB\n", " allocs estimate: 26\n", " --------------\n", " minimum time: 68.818 ms (2.48% GC)\n", " median time: 77.068 ms (2.30% GC)\n", " mean time: 81.321 ms (2.25% GC)\n", " maximum time: 131.659 ms (2.40% GC)\n", " --------------\n", " samples: 62\n", " evals/sample: 1" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@benchmark logpdf_mvn_1(y, Σ)" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 15.27 MiB\n", " allocs estimate: 14\n", " --------------\n", " minimum time: 10.761 ms (0.00% GC)\n", " median time: 12.785 ms (13.40% GC)\n", " mean time: 12.987 ms (10.33% GC)\n", " maximum time: 19.108 ms (10.46% GC)\n", " --------------\n", " samples: 385\n", " evals/sample: 1" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@benchmark logpdf_mvn_2(y, Σ)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 7.64 MiB\n", " allocs estimate: 10\n", " --------------\n", " minimum time: 8.802 ms (0.00% GC)\n", " median time: 9.325 ms (0.00% GC)\n", " mean time: 10.130 ms (7.36% GC)\n", " maximum time: 14.764 ms (15.79% GC)\n", " --------------\n", " samples: 493\n", " evals/sample: 1" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@benchmark logpdf_mvn_3(y, Σ)" ] }, { "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\n", "\n", "* Cholesky decomposition is **one** approach to solve linear regression.\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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Further reading\n", "\n", "* Section 7.7 of [Numerical Analysis for Statisticians](http://ucla.worldcat.org/title/numerical-analysis-for-statisticians/oclc/793808354&referer=brief_results) of Kenneth Lange (2010).\n", "\n", "* Section II.5.3 of [Computational Statistics](http://ucla.worldcat.org/title/computational-statistics/oclc/437345409&referer=brief_results) by James Gentle (2010).\n", "\n", "* Section 4.2 of [Matrix Computation](http://catalog.library.ucla.edu/vwebv/holdingsInfo?bibId=7122088) by Gene Golub and Charles Van Loan (2013)." ] }, { "cell_type": "code", "execution_count": 31, "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": "118px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "threshold": 4, "toc_cell": true, "toc_section_display": "block", "toc_window_display": false, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 2 }