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. 

## 4.4 Scalability Test 

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 

