{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Linear Algebra in Julia" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The linear algebra facilities in [Julia](http://julialang.org) are a superset of the base facilities in [R](http://www.R-project.org) and Matlab/Octave. As in `R` a vector is considered as a column vector for purposes of linear algebra. Its transpose is a 1$\\times$n matrix. The equivalent of `dim` in `R` is `size`." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(6,)" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "v1 = collect(1:6);\n", "size(v1) # size always returns a tuple" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "6" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "length(v1)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3x2 Array{Int64,2}:\n", " 1 4\n", " 2 5\n", " 3 6" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m1 = reshape(collect(1:6),(3,2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Because `size` returns a tuple the number of rows and columns can be assigned in one statement. A common idiom is" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(3,2)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m,n = size(m1)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING: readbytes is deprecated, use read instead.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Number of rows = 3 and number of columns = 2\n" ] } ], "source": [ "println(\"Number of rows = $m and number of columns = $n\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As in `R`, matrices and higher-dimensional arrays are stored in _column-major_ ordering. Transpose is indicated by `'`. The `*` operator applied to arrays denotes matrix multiplication. Element-wise multiplication is indicated by `.*`. In general a prefix `.` to an operator makes it behave element-wise. Constructors of arrays with specific content include `ones()`, `zeros()` and `eye()` (the identity matrix)." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 5.0\n", " 7.0\n", " 9.0" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ " in depwarn" ] } ], "source": [ "m1 * ones(2)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "(" ] }, { "data": { "text/plain": [ "2x2 Array{Int64,2}:\n", " 14 32\n", " 32 77" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "::" ] } ], "source": [ "m1'*m1" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3x3 Array{Float64,2}:\n", " 1.0 0.0 0.0\n", " 0.0 1.0 0.0\n", " 0.0 0.0 1.0" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "ASCIIString, " ] } ], "source": [ "eye(3)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "::" ] }, { "data": { "text/plain": [ "3x4 Array{Float64,2}:\n", " 1.0 0.0 0.0 0.0\n", " 0.0 1.0 0.0 0.0\n", " 0.0 0.0 1.0 0.0" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "Symbol" ] } ], "source": [ "eye(3,4)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2x3 Array{Int64,2}:\n", " 1 1 1\n", " 1 1 1" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ones(Int,(2,3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The equivalent of `R`'s `cbind` and `rbind` functions are `hcat` and `vcat`." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3x2 Array{Float64,2}:\n", " 1.0 -1.0\n", " 1.0 0.0\n", " 1.0 1.0" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ ") at ./deprecated.jl:" ] } ], "source": [ "m2 = hcat(ones(3),collect(-1.:1.))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`vcat` is used for concatenating vectors." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "64" ] }, { "data": { "text/plain": [ "5-element Array{Float64,1}:\n", " 1.0\n", " 2.0\n", " 3.0\n", " 1.0\n", " 1.0" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "vcat([1:3;],ones(2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In some circumstances you can omit the `*` in a multiplication." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2x2 Array{Float64,2}:\n", " 3.0 0.0\n", " 0.0 2.0" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ " in readbytes" ] } ], "source": [ "m22 = m2'm2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Unlike in `R`, the `Julia` expression `X'X` is evaluated using only one copy of `X` and will always produce a symmetric matrix. In `R` the function `crossprod` must be called to do this. (Be careful, in `Julia` there is an operation called a _cross product_ but it is the cross product from Physics.) " ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "(::" ] }, { "data": { "text/latex": [ "\\begin{verbatim}\n", "cross(x, y)\n", "×(x,y)\n", "\\end{verbatim}\n", "Compute the cross product of two 3-vectors.\n" ], "text/markdown": [ "```\n", "cross(x, y)\n", "×(x,y)\n", "```\n", "\n", "Compute the cross product of two 3-vectors.\n" ], "text/plain": [ "```\n", "cross(x, y)\n", "×(x,y)\n", "```\n", "\n", "Compute the cross product of two 3-vectors.\n" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "Base.PipeEndpoint, ::Vararg" ] } ], "source": [ "?×" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solving linear systems of equations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The backslash operator, `\\`, denotes the solution of a linear system of equations of the form $Ax=b$." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "{" ] }, { "data": { "text/plain": [ "2-element Array{Float64,1}:\n", " 1.0\n", " 0.0" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "Any" ] } ], "source": [ "m22\\(m2'ones(3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With a rectangular left-hand side `\\` denotes a least-squares solution" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2-element Array{Float64,1}:\n", " 1.0\n", " -0.0" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m2\\ones(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that this solution is a bit different from the previous solution because of numerical imprecision in the computation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Factorizations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Those who haven't studied numerical linear algebra assume that a linear system of equations is solved as $A^{-1}b$. In practice it is almost never necessary to evaluate the inverse of a matrix. Instead the matrix is factored into a product of other matrices that have desirable properties such as being triangular.\n", "\n", "Such factorizations include the singular value decomposition (SVD), the QR decomposition (with or without pivoting), the Cholesky factorization (with or without pivoting) for a positive-semidefinite symmetric matrix, the LU factorization and the Eigen decomposition. The names of the constructors of the decomposition usually end in `fact`." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3x2 Array{Float64,2}:\n", " 1.0 1.0\n", " 1.0 2.0\n", " 1.0 3.0" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "}) at ./deprecated.jl" ] } ], "source": [ "m3 = hcat(ones(3),[1:3;])" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":" ] }, { "data": { "text/plain": [ "2x2 Array{Float64,2}:\n", " 3.0 6.0\n", " 6.0 14.0" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "30" ] } ], "source": [ "m33 = m3'm3" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:\n", "2x2 UpperTriangular{Float64,Array{Float64,2}}:\n", " 1.73205 3.4641 \n", " ⋅ 1.41421" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cf33 = cholfact(m33)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\n", " in send_stream" ] }, { "data": { "text/plain": [ "2x2 LowerTriangular{Float64,Array{Float64,2}}:\n", " 1.73205 ⋅ \n", " 3.4641 1.41421" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "(" ] } ], "source": [ "cf33[:L] # extract the lower Cholesky factor" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "::" ] }, { "data": { "text/plain": [ "2x2 UpperTriangular{Float64,Array{Float64,2}}:\n", " 1.73205 3.4641 \n", " ⋅ 1.41421" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cf33[:U] # extract the upper triangular Cholesky factor, cf33[:L]'" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Base." ] }, { "data": { "text/plain": [ "2x2 Array{Float64,2}:\n", " 3.0 6.0\n", " 6.0 14.0" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "PipeEndpoint" ] } ], "source": [ "full(cf33) # the matrix (up to round-off) represented by the factorization" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ", " ] }, { "data": { "text/plain": [ "5.999999999999983" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "det(cf33) # determinant of full(cf33), calculated from det(cf33[:U])" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1.7917594692280523" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "::ASCIIString)" ] } ], "source": [ "logdet(cf33) # more practical for positive semidefinite" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Although one rarely needs the inverse of a general matrix, it is sometimes useful to get the inverse of the matrix represented by the Cholesky factorization. This is feasible because the factor is triangular and determining the inverse of a triangular matrix is easier than doing so for a general matrix. Suppose we want the estimated covariance matrix of the coefficient estimates." ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ " at " ] }, { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 0.268839\n", " 0.129757\n", " 0.348011" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "/home/bates/.julia/v0.5/IJulia/src/stdio.jl" ] } ], "source": [ "y = rand(3)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2-element Array{Float64,1}:\n", " 0.169697 \n", " 0.0395858" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "β = cf33\\(m3'y)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 0.209283\n", " 0.248869\n", " 0.288455" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ŷ = m3*β # note you have to construct ŷ as y\\hat not \\haty" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 0.0595559\n", " -0.119112 \n", " 0.0595559" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "resid = y - ŷ" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.021281438491229213" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s² = sumabs2(resid)/(3-2)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2x2 Array{Float64,2}:\n", " 0.0496567 -0.0212814\n", " -0.0212814 0.0106407" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ ":25\n", " in watch_stream(::Base.PipeEndpoint" ] } ], "source": [ "vcv = s²*inv(cf33)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The pivoted Cholesky decomposition has better numerical properties than does the unpivoted version and is a _rank-revealing_ decomposition." ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ", " ] }, { "data": { "text/plain": [ "Base.LinAlg.CholeskyPivoted{Float64,Union{DenseArray{Float64,2},SubArray{Float64,2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}}}(2x2 Array{Float64,2}:\n", " 3.74166 1.60357 \n", " 6.0 0.654654,'U',Int32[2,1],2,0.0,0)" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "::" ] } ], "source": [ "cf33 = cholfact(m33, :U, Val{true})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At present this factorization doesn't have an attractive `show` method." ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2x2 UpperTriangular{Float64,Array{Float64,2}}:\n", " 3.74166 1.60357 \n", " ⋅ 0.654654" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "ASCIIString)" ] } ], "source": [ "cf33[:U]" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ " at " ] }, { "data": { "text/plain": [ "2-element Array{Int32,1}:\n", " 2\n", " 1" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "/home/bates/.julia/v0.5/IJulia/src/stdio.jl" ] } ], "source": [ "cf33[:p] # the pivot vector" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2x2 Array{Float64,2}:\n", " 0.0 1.0\n", " 1.0 0.0" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cf33[:P] # the pivot matrix" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rank(cf33)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The QR decomposition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The QR decomposition is a direct (i.e. non-iterative) decomposition used, among other things, to solve least squares problems. It decomposes a rectangular `X` into an orthogonal matrix `Q` and an upper triangular matrix `R`. The matrix `Q`, which could be very large if the number of rows of `X` is large, is stored implicitly." ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":41\n", " in (::" ] }, { "data": { "text/plain": [ "Base.LinAlg.QRCompactWY{Float64,Array{Float64,2}}(3x2 Array{Float64,2}:\n", " -1.73205 -3.4641 \n", " 0.366025 -1.41421 \n", " 0.366025 0.767327,2x2 Array{Float64,2}:\n", " 1.57735 -1.28446\n", " 6.93333e-310 1.25882)" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "IJulia" ] } ], "source": [ "qr3 = qrfact(m3)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "." ] }, { "data": { "text/plain": [ "2x2 Array{Float64,2}:\n", " -1.73205 -3.4641 \n", " 0.0 -1.41421" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "qr3[:R]" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3x3 Base.LinAlg.QRCompactWYQ{Float64,Array{Float64,2}}:\n", " -0.57735 0.707107 0.408248\n", " -0.57735 1.11022e-16 -0.816497\n", " -0.57735 -0.707107 0.408248" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "qr3[:Q] # the implicit representation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The QR factorization can be used directly for least squares solutions." ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "##4#8)()" ] }, { "data": { "text/plain": [ "2-element Array{Float64,1}:\n", " 0.169697 \n", " 0.0395858" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ " at " ] } ], "source": [ "β = qr3\\y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What are called the `effects` in an `R` `lm` object are `Q'y`" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "./task.jl" ] }, { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " -0.431054 \n", " -0.0559827\n", " 0.145882 " ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "qr3[:Q]'y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is also a pivoted QR decomposition which is rank-revealing." ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Base.LinAlg.QRPivoted{Float64,Array{Float64,2}}(3x2 Array{Float64,2}:\n", " -3.74166 -1.60357 \n", " 0.421793 0.654654\n", " 0.63269 0.859768,[1.26726,1.14995],Int32[2,1])" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "qr3 = qrfact(m3, Val{true})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The SVD" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The singular value decomposition represents `X` as `U*S*V'` where `U` and `V` are orthogonal and `S` is diagonal. It is related to the eigen decomposition in that the eigen decomposition of `X'X` is `V*S*S*V'`. The eigenvalues of `X'X` are the squares of the singular values of `X`." ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":431\n", "while loading In[5], in expression starting on line 1\n" ] }, { "data": { "text/plain": [ "Base.LinAlg.SVD{Float64,Float64,Array{Float64,2}}(3x2 Array{Float64,2}:\n", " -0.32311 0.853776\n", " -0.547507 0.18322 \n", " -0.771904 -0.487337,[4.07914,0.600491],2x2 Array{Float64,2}:\n", " -0.402663 -0.915348\n", " 0.915348 -0.402663)" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" }, { "name": "stdout", "output_type": "stream", "text": [ "search: ×\n", "\n" ] } ], "source": [ "sv3 = svdfact(m3)" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3x2 Array{Float64,2}:\n", " -0.32311 0.853776\n", " -0.547507 0.18322 \n", " -0.771904 -0.487337" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sv3[:U]" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2x2 Array{Float64,2}:\n", " -0.402663 0.915348\n", " -0.915348 -0.402663" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sv3[:V]" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2-element Array{Float64,1}:\n", " 4.07914 \n", " 0.600491" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sv3[:S] # vector of singular values" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2x2 Array{Float64,2}:\n", " -0.402663 -0.915348\n", " 0.915348 -0.402663" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sv3[:Vt]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is the _thin_ SVD. If you want the full matrix `U` use" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Base.LinAlg.SVD{Float64,Float64,Array{Float64,2}}(3x3 Array{Float64,2}:\n", " -0.32311 0.853776 0.408248\n", " -0.547507 0.18322 -0.816497\n", " -0.771904 -0.487337 0.408248,[4.07914,0.600491],2x2 Array{Float64,2}:\n", " -0.402663 -0.915348\n", " 0.915348 -0.402663)" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sv3 = svdfact(m3; thin=false)" ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3x3 Array{Float64,2}:\n", " -0.32311 0.853776 0.408248\n", " -0.547507 0.18322 -0.816497\n", " -0.771904 -0.487337 0.408248" ] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sv3[:U]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The eigen decomposition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `eigfact` function has a special method for real symmetric matrices" ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Base.LinAlg.Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}([0.36059,16.6394],2x2 Array{Float64,2}:\n", " -0.915348 0.402663\n", " 0.402663 0.915348)" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ef3 = eigfact(m33)" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2-element Array{Float64,1}:\n", " 0.36059\n", " 16.6394 " ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ef3[:values]" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2-element Array{Float64,1}:\n", " 16.6394 \n", " 0.36059" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "abs2(sv3[:S]) # same as the eigenvalues but in the opposite order" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2x2 Array{Float64,2}:\n", " -0.915348 0.402663\n", " 0.402663 0.915348" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ef3[:vectors]" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2x2 Array{Float64,2}:\n", " -0.402663 0.915348\n", " -0.915348 -0.402663" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sv3[:V] # same as the eigenvectors of m33 with permuted rows and columns" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For general matrices we must allow for complex eigenvalues and eigenvectors. If all the eigenvalues and eigenvectors are real they are returned as such" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Base.LinAlg.Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}([16.1168,-1.11684,-5.70069e-16],3x3 Array{Float64,2}:\n", " -0.464547 -0.882906 0.408248\n", " -0.570796 -0.23952 -0.816497\n", " -0.677044 0.403865 0.408248)" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eigfact(reshape([1:9;],(3,3)))" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}}(Complex{Float64}[0.0+1.0im,0.0-1.0im],2x2 Array{Complex{Float64},2}:\n", " 0.707107+0.0im 0.707107-0.0im \n", " 0.0+0.707107im 0.0-0.707107im)" ] }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eigfact([0 1;-1 0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The LU factorization" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The LU factorization is the most general direct factorization. It is related to Gaussian elimination." ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Base.LinAlg.LU{Float64,Array{Float64,2}}(3x3 Array{Float64,2}:\n", " 3.0 6.0 9.0\n", " 0.333333 2.0 4.0\n", " 0.666667 0.5 0.0,Int32[3,3,3],3)" ] }, "execution_count": 56, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lu4 = lufact(reshape([1:9;],(3,3)))" ] }, { "cell_type": "code", "execution_count": 57, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3x3 Array{Float64,2}:\n", " 1.0 0.0 0.0\n", " 0.333333 1.0 0.0\n", " 0.666667 0.5 1.0" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lu4[:L] # unit lower triangular factor" ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3x3 Array{Float64,2}:\n", " 3.0 6.0 9.0\n", " 0.0 2.0 4.0\n", " 0.0 0.0 0.0" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lu4[:U] # upper triangular factor. Note that it is singular." ] }, { "cell_type": "code", "execution_count": 59, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Int32,1}:\n", " 3\n", " 1\n", " 2" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lu4[:p] # permutation vector" ] }, { "cell_type": "code", "execution_count": 60, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3x3 Array{Float64,2}:\n", " 0.0 0.0 1.0\n", " 1.0 0.0 0.0\n", " 0.0 1.0 0.0" ] }, "execution_count": 60, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lu4[:P] # permutation matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Extracting submatrices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As in `R` elements or subarrays are extracted with using indices. To indicate all values of a given index use `:`" ] }, { "cell_type": "code", "execution_count": 61, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " -1.0\n", " 0.0\n", " 1.0" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m2[:,2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are two alternatives to indexing, the use of `sub` or `slice`." ] }, { "cell_type": "code", "execution_count": 62, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element SubArray{Float64,1,Array{Float64,2},Tuple{Colon,Int64},2}:\n", " -1.0\n", " 0.0\n", " 1.0" ] }, "execution_count": 62, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sub(m2,:,2)" ] }, { "cell_type": "code", "execution_count": 63, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element SubArray{Float64,1,Array{Float64,2},Tuple{Colon,Int64},2}:\n", " -1.0\n", " 0.0\n", " 1.0" ] }, "execution_count": 63, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slice(m2, :, 2)" ] }, { "cell_type": "code", "execution_count": 64, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1x2 SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},Colon},2}:\n", " 1.0 0.0" ] }, "execution_count": 64, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sub(m2, 2, :)" ] }, { "cell_type": "code", "execution_count": 65, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2-element SubArray{Float64,1,Array{Float64,2},Tuple{Int64,Colon},2}:\n", " 1.0\n", " 0.0" ] }, "execution_count": 65, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slice(m2, 2, :)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The difference between `sub` and `slice` is that `slice` always drops singleton dimensions whereas `sub` only drops trailing singleton dimensions.\n", "\n", "Both `sub` and `slice` return views into the original array. That is, they do not copy. At present using indices produces a copy but that will change." ] }, { "cell_type": "code", "execution_count": 66, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Ptr{Float64} @0x00007fa1a0221000" ] }, "execution_count": 66, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pointer(m2)" ] }, { "cell_type": "code", "execution_count": 67, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Ptr{Float64} @0x00007fa1a26e9360" ] }, "execution_count": 67, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pointer(m2[:,1])" ] }, { "cell_type": "code", "execution_count": 68, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Ptr{Float64} @0x00007fa1a0221000" ] }, "execution_count": 68, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pointer(sub(m2, :, 1))" ] }, { "cell_type": "code", "execution_count": 69, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Ptr{Float64} @0x00007fa1a0221000" ] }, "execution_count": 69, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pointer(slice(m2, 1, :))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To convert a matrix or higher-order array to a vector use `vec`" ] }, { "cell_type": "code", "execution_count": 70, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "6-element Array{Float64,1}:\n", " 1.0\n", " 1.0\n", " 1.0\n", " -1.0\n", " 0.0\n", " 1.0" ] }, "execution_count": 70, "metadata": {}, "output_type": "execute_result" } ], "source": [ "vec(m2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result of an indexing operation is a copy of the elements from the array." ] }, { "cell_type": "code", "execution_count": 71, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "999.0" ] }, "execution_count": 71, "metadata": {}, "output_type": "execute_result" } ], "source": [ "v2 = m2[:,2]\n", "v2[2] = 999." ] }, { "cell_type": "code", "execution_count": 72, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3x2 Array{Float64,2}:\n", " 1.0 -1.0\n", " 1.0 0.0\n", " 1.0 1.0" ] }, "execution_count": 72, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m2" ] }, { "cell_type": "code", "execution_count": 73, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " -1.0\n", " 999.0\n", " 1.0" ] }, "execution_count": 73, "metadata": {}, "output_type": "execute_result" } ], "source": [ "v2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `sub` function produces an array view which accesses the actual elements in the original array." ] }, { "cell_type": "code", "execution_count": 74, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element SubArray{Float64,1,Array{Float64,2},Tuple{Colon,Int64},2}:\n", " -1.0\n", " 0.0\n", " 1.0" ] }, "execution_count": 74, "metadata": {}, "output_type": "execute_result" } ], "source": [ "v2 = sub(m2,:,2)" ] }, { "cell_type": "code", "execution_count": 75, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1.0" ] }, "execution_count": 75, "metadata": {}, "output_type": "execute_result" } ], "source": [ "v2[2] = 1." ] }, { "cell_type": "code", "execution_count": 76, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3x2 Array{Float64,2}:\n", " 1.0 -1.0\n", " 1.0 1.0\n", " 1.0 1.0" ] }, "execution_count": 76, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `sub` function is preferred over explicit indexing if you need a subarray that you are not going to modify (`sub` does not create unnecessary copies) and for when you really do want to change the elements of a subarray." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Underlying numerical codes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Most linear algebra in `Julia` is performed on `Matrix{Float64}` objects using the LAPACK and BLAS software. By default `Julia` is compiled against OpenBLAS, an accelerated multi-threaded BLAS that also provides accelerated versions of some of the LAPACK functions. For example the LU decomposition and solutions of linear systems using LU are accelerated because they constitute the LAPACK benchmark. `Julia` can also be compiled with Intel's MKL (Math Kernel Library) BLAS, when they are available, to give accelerated performance. Understandably these are tuned for best performance on Intel processors.\n", "\n", "Julia provides thin wrappers around many LAPACK and BLAS routines if you want maximal control. These mimic the BLAS and LAPACK routines without the annoyingly redundant arguments. Note that many internal linear algebra functions provide mutating versions." ] }, { "cell_type": "code", "execution_count": 77, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " BLAS 2515 KB Module\n", " asum 0 bytes Base.LinAlg.BLAS.#asum\n", " blascopy! 0 bytes Base.LinAlg.BLAS.#blascopy!\n", " dotc 0 bytes Base.LinAlg.BLAS.#dotc\n", " dotu 0 bytes Base.LinAlg.BLAS.#dotu\n", " gbmv 0 bytes Base.LinAlg.BLAS.#gbmv\n", " gbmv! 0 bytes Base.LinAlg.BLAS.#gbmv!\n", " gemm 0 bytes Base.LinAlg.BLAS.#gemm\n", " gemm! 0 bytes Base.LinAlg.BLAS.#gemm!\n", " gemv 0 bytes Base.LinAlg.BLAS.#gemv\n", " gemv! 0 bytes Base.LinAlg.BLAS.#gemv!\n", " ger! 0 bytes Base.LinAlg.BLAS.#ger!\n", " hemm 0 bytes Base.LinAlg.BLAS.#hemm\n", " hemm! 0 bytes Base.LinAlg.BLAS.#hemm!\n", " hemv 0 bytes Base.LinAlg.BLAS.#hemv\n", " hemv! 0 bytes Base.LinAlg.BLAS.#hemv!\n", " her! 0 bytes Base.LinAlg.BLAS.#her!\n", " her2k 0 bytes Base.LinAlg.BLAS.#her2k\n", " her2k! 0 bytes Base.LinAlg.BLAS.#her2k!\n", " herk 0 bytes Base.LinAlg.BLAS.#herk\n", " herk! 0 bytes Base.LinAlg.BLAS.#herk!\n", " iamax 0 bytes Base.LinAlg.BLAS.#iamax\n", " nrm2 0 bytes Base.LinAlg.BLAS.#nrm2\n", " sbmv 0 bytes Base.LinAlg.BLAS.#sbmv\n", " sbmv! 0 bytes Base.LinAlg.BLAS.#sbmv!\n", " scal 0 bytes Base.LinAlg.BLAS.#scal\n", " scal! 0 bytes Base.LinAlg.BLAS.#scal!\n", " symm 0 bytes Base.LinAlg.BLAS.#symm\n", " symm! 0 bytes Base.LinAlg.BLAS.#symm!\n", " symv 0 bytes Base.LinAlg.BLAS.#symv\n", " symv! 0 bytes Base.LinAlg.BLAS.#symv!\n", " syr! 0 bytes Base.LinAlg.BLAS.#syr!\n", " syr2k 0 bytes Base.LinAlg.BLAS.#syr2k\n", " syr2k! 0 bytes Base.LinAlg.BLAS.#syr2k!\n", " syrk 0 bytes Base.LinAlg.BLAS.#syrk\n", " syrk! 0 bytes Base.LinAlg.BLAS.#syrk!\n", " trmm 0 bytes Base.LinAlg.BLAS.#trmm\n", " trmm! 0 bytes Base.LinAlg.BLAS.#trmm!\n", " trmv 0 bytes Base.LinAlg.BLAS.#trmv\n", " trmv! 0 bytes Base.LinAlg.BLAS.#trmv!\n", " trsm 0 bytes Base.LinAlg.BLAS.#trsm\n", " trsm! 0 bytes Base.LinAlg.BLAS.#trsm!\n", " trsv 0 bytes Base.LinAlg.BLAS.#trsv\n", " trsv! 0 bytes Base.LinAlg.BLAS.#trsv!\n" ] } ], "source": [ "whos(BLAS)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Operations like `A'B` and `A\\b` call functions with names like `Ac_mul_B` and `A_ldiv_B`. There are mutating versions of these functions. The `Ac` denotes `A'` which, in general, is the conjugate transpose of `A`. For complex `A` there is a distinction between `A'`, the conjugate transpose, and `A.'`, the transpose. Usually you want `A'`. For real matrices `A'` and `A.'` are the same." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Special Matrix types" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that the type of `cf33[:U]` from the first Cholesky example is" ] }, { "cell_type": "code", "execution_count": 78, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "UpperTriangular{Float64,Array{Float64,2}}" ] }, "execution_count": 78, "metadata": {}, "output_type": "execute_result" } ], "source": [ "typeof(cholfact(m33)[:U])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Many operations can be steamlined if it is known that one of the operands is a triangular matrix or, even more specific, a unit triangular matrices. There are also special types for `Symmetric` and `Diagonal` matrices. Ongoing work will incorporate these matrix types more deeply into the `Julia` system." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A more complex example" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When fitting a linear mixed model, the core computation can be expressed as solving a penalized least squares problem. In the `MixedModels` package there are several different `PLSSolver` types defined to facilitate this computation. When there is only a single random-effects term the system to be solved is block diagonal.\n", "\n", "The `PLSOne` type is defined as \n", "```julia\n", "type PLSOne <: PLSSolver # Solver for models with a single random-effects term\n", " Ad::Array{Float64,3} # diagonal blocks\n", " Ab::Array{Float64,3} # base blocks\n", " At::Symmetric{Float64} # lower right block\n", " Ld::Array{Float64,3} # diagonal blocks\n", " Lb::Array{Float64,3} # base blocks\n", " Lt::Base.LinAlg.Cholesky{Float64} # lower right triangle\n", " Zt::SparseMatrixCSC\n", "end\n", "```\n", "and is updated from a triangular matrix λ as\n", "```julia\n", "## update!(s,lambda)->s : update Ld, Lb and Lt\n", "function update!(s::PLSOne,λ::Triangular)\n", " m,n,l = size(s)\n", " n == size(λ,1) || error(\"Dimension mismatch\")\n", " Lt = copy!(s.Lt.UL,s.At.S)\n", " Lb = copy!(s.Lb,s.Ab)\n", " if n == 1 # shortcut for 1×1 λ\n", " lam = λ[1,1]\n", " Ld = map!(x -> sqrt(x*lam*lam + 1.), s.Ld, s.Ad)\n", " Lb = scale!(reshape(Lb,(m,n*l)),lam ./ vec(Ld))\n", " BLAS.syrk!('L','N',-1.0,Lb,1.0,Lt)\n", " else\n", " Ac_mul_B!(λ,reshape(copy!(s.Ld,s.Ad),(n,n*l)))\n", " for k in 1:l\n", " wL = A_mul_B!(sub(s.Ld,:,:,k),λ)\n", " for j in 1:n # Inflate the diagonal\n", " wL[j,j] += 1.\n", " end\n", " _, info = LAPACK.potrf!('L',wL) # i'th diagonal block of L\n", " info == 0 || error(\"Cholesky failure at L diagonal block $k\")\n", " Base.LinAlg.A_rdiv_Bc!(A_mul_B!(sub(s.Lb,:,:,k),λ),Triangular(wL,:L,false))\n", " end\n", " BLAS.syrk!('L','N',-1.0,reshape(Lb,(m,n*l)),1.,Lt)\n", " end\n", " _, info = LAPACK.potrf!('L',Lt)\n", " info == 0 || error(\"downdated X'X is not positive definite\")\n", " s\n", "end\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sparse Matrices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sparse matrices and sparse vectors are available in the base `Julia` system but their description will need to wait for another time. " ] }, { "cell_type": "code", "execution_count": 79, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "search: " ] }, { "data": { "text/latex": [ "\\begin{verbatim}\n", "sparse(A)\n", "\\end{verbatim}\n", "Convert an AbstractMatrix \\texttt{A} into a sparse matrix.\n", "\\begin{verbatim}\n", "sparse(I,J,V,[m,n,combine])\n", "\\end{verbatim}\n", "Create a sparse matrix \\texttt{S} of dimensions \\texttt{m x n} such that \\texttt{S[I[k], J[k]] = V[k]}. The \\texttt{combine} function is used to combine duplicates. If \\texttt{m} and \\texttt{n} are not specified, they are set to \\texttt{maximum(I)} and \\texttt{maximum(J)} respectively. If the \\texttt{combine} function is not supplied, \\texttt{combine} defaults to \\texttt{+} unless the elements of \\texttt{V} are Booleans in which case \\texttt{combine} defaults to \\texttt{|}. All elements of \\texttt{I} must satisfy \\texttt{1 <= I[k] <= m}, and all elements of \\texttt{J} must satisfy \\texttt{1 <= J[k] <= n}.\n" ], "text/markdown": [ "```\n", "sparse(A)\n", "```\n", "\n", "Convert an AbstractMatrix `A` into a sparse matrix.\n", "\n", "```\n", "sparse(I,J,V,[m,n,combine])\n", "```\n", "\n", "Create a sparse matrix `S` of dimensions `m x n` such that `S[I[k], J[k]] = V[k]`. The `combine` function is used to combine duplicates. If `m` and `n` are not specified, they are set to `maximum(I)` and `maximum(J)` respectively. If the `combine` function is not supplied, `combine` defaults to `+` unless the elements of `V` are Booleans in which case `combine` defaults to `|`. All elements of `I` must satisfy `1 <= I[k] <= m`, and all elements of `J` must satisfy `1 <= J[k] <= n`.\n" ], "text/plain": [ "```\n", "sparse(A)\n", "```\n", "\n", "Convert an AbstractMatrix `A` into a sparse matrix.\n", "\n", "```\n", "sparse(I,J,V,[m,n,combine])\n", "```\n", "\n", "Create a sparse matrix `S` of dimensions `m x n` such that `S[I[k], J[k]] = V[k]`. The `combine` function is used to combine duplicates. If `m` and `n` are not specified, they are set to `maximum(I)` and `maximum(J)` respectively. If the `combine` function is not supplied, `combine` defaults to `+` unless the elements of `V` are Booleans in which case `combine` defaults to `|`. All elements of `I` must satisfy `1 <= I[k] <= m`, and all elements of `J` must satisfy `1 <= J[k] <= n`.\n" ] }, "execution_count": 79, "metadata": {}, "output_type": "execute_result" }, { "name": "stdout", "output_type": "stream", "text": [ "sparse sparsevec SparseVector SparseArrays SparseMatrixCSC issparse\n", "\n" ] } ], "source": [ "?sparse" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.5.0-dev", "language": "julia", "name": "julia-0.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.5.0" } }, "nbformat": 4, "nbformat_minor": 0 }