If H were dense Nk + p square, symmetric matrix, then, holding k and p constant, the complexity order of the Cholesky decomposition is N 3 (Golub et al. 1996). There are a number of different algorithms that one can use for decomposing sparse Hessians (Davis 2006). The typical strategy is to first permute the rows and columns of H to minimize the number of nonzero elements in L , and then compute the sparsity pattern. This part can be done just once. With the sparsity pattern in hand, the next step is to compute those nonzero elements in L . The time for this step grows with the sum of the squares of the number of nonzero elements in each column of L (Davis 2006). Because each additional household adds k columns to L , with an average of p + 1 $_{2}$( k + 1 ) nonzero elements per column, we can compute the sparse Cholesky decomposition in time that is linear in N . Software for sparse Cholesky decompositions is widely available. Next, we provide some empirical evidence of scalability through a simulation study. For a hypo- thetical dataset with N households, let y$_{i}$ be the number of times household i visits a store during a T week period. The probability of a visit in a single week is p$_{i}$ , where logit p$_{i}$ = b $^{′}$$_{i}$x$_{i}$ , and x$_{i}$ is a vector of k covariates. The distribution of b i across the population is MVN with mean ¯ b and covariance S . In all conditions of the test, we set k = 3 and T = 52, and vary the number of households by setting N to one of 10 discrete values from 500 to 50,000. The “true” values of ¯ b are -10, 0 and 10, and the “true” S is 0.1 I . We place weakly informative priors on both ¯ b and S . In Figure 8, we plot the average time, across 100 replications, to compute the log posterior, the gradient, and the Hessian. As expected, each of these computations grows linearly in N . In Figure 9, we plot average times for the steps involved in sampling from an MVN: adding a vector to columns of a dense matrix, computing a sparse Cholesky decomposition, multiplying a sparse triangular matrix by a dense matrix, sampling standard normal random variates, and solving a sparse triangular linear system. Again, we see that the time for all of these steps is linear in N . Table 2 summarizes the acceptance rates and scale factors when generating 50 samples from the posterior for different values of N . Although there is a weak trend of increasing acceptance rates with N ,we cannot say with any certainty that acceptance rates will always be either larger or smaller for larger datasets. The acceptance rate could be influenced by using different scale 25 4.4 Scalability Test