algorithms of this kind are described in Nocedal et al. (2006), and are available for many technical computing platforms. Once the algorithm finds the posterior mode, there remains the need to compute the Hessian exactly. 

The second approach is to run a quasi-Newton algorithm, and compute the Hessian at each it- eration explicitly, but store the Hessian in a sparse format. The trustOptim package for R (Braun 2014) implements a trust region algorithm that exploits the sparsity of the Hessian. The user can sup- ply a Hessian that is derived analytically, computed using AD, or estimated numerically using sparseHessianFD . Since memory requirements and matrix computation costs will grow only lin- early in N , finding the posterior mode becomes feasible for large problems, compared to similar algorithms that ignore sparsity. 

We should note that we cannot predict the time to convergence for general problems. Log posteriors with ridges or plateaus, or that require extensive computation themselves, may still take a long time to find the local optimum. Whether any mode-finding algorithm is “fast enough” depends on the specific application. However, if one optimization algorithm has difficulty finding a mode, another algorithm may do better. 

## 4.3 Sampling from an MVN distribution 

Once we find the posterior mode, and the Hessian at the mode, generating proposal samples from an MVN( q $^{∗}$, sH − $^{1}$) distribution is straightforward. Let 1 $_{s}$H = LL ′ represent the Cholesky decomposition of the precision of the proposal, and let z be a vector of Nk + p samples from a standard normal distribution. To sample q from an MVN, solve the triangular linear system L ′ q = x , and then add q $^{∗}$. Since E ( z ) = 0, E ($_{q}$) = q $^{∗}$, and since E ( zz $^{′}$) = I , cov ($_{q}$) = L ′− $^{1}$L − 1 = ( LL $^{′}$)− 1 = sH − $^{1}$. 

Because L is sparse, the costs of both solving the triangular system, and the premultiplication, grow linearly with the number of nonzero elements, which itself grows linearly in N (Davis 2006). If L were dense, then the cost of solving the triangular system would grow quadratically in N . Furthermore, computing the MVN density would involve premultiplying z by a triangular matrix, whose cost is cubic in N (Golub et al. 1996). 

Computation of the Cholesky decomposition can also benefit from the sparsity of the Hessian. 

