We gratefully thank the Journal of Statistical Computation and Simulation and copyright owner OPA (Overseas Publishers Association) NV for allowing us to distribute this electronic version of the manuscript.

By

R. Kelley Pace

Ronald P. Barry

This paper develops methods for quickly computing maximum likelihood conditional autoregressions (CARs). By using sparse matrix methods, reorganizing the sum-of-squared errors function to avoid unnecessary calculations, and precomputing a set of determinants, simulations of large CARs become possible. As an illustration of the power of these approaches, a simulation of 250 CARs of 2,905 observations can take fewer than three minutes on a personal computer, despite the necessity of evaluating 100 determinants of 2,905 by 2,905 matrices. The computation of each estimate via examining the profile likelihood sampled at 100 points avoids problems of local optima. Simulating estimates avoids other problems associated with the traditional information matrix approach to inference.

KEY WORDS: Conditionally specified gaussian, Jacobian, GIS, vectorized profile likelihoods, spatial statistics, sparse matrices.

 

The advent of the global positioning system (GPS), which can measure precisely geographic locations, coupled with the continuing development of geographic information systems (GIS) have led to an explosion of spatial data. For example, inexpensive phone books on CD-ROM now provide the latitude and longitude of over 100 million US residences with various demographic variables. The Bureau of the Census provides extensive demographic data for over 250,000 geographic points (census block groups). Remote sensing via satellites and digital imaging yields even vaster quantities of spatial information.

While the availability of spatial information has grown tremendously, the ability to process it has not grown pari passu. The necessary reliance on a n by n matrix, where n represents the number of observations, greatly impedes traditional spatial statistics from handling larger sample sizes. Spatial estimators rely upon either a n by n variance-covariance matrix (e.g., kriging), the inverse of the variance-covariance matrix (e.g., conditional autoregressions (CARs)), or the square root of the inverse variance-covariance matrix (e.g., simultaneous autoregressions (SARs)). As n becomes large, the size of these matrices becomes astronomical. For example, a 50,000 observation spatial regression would require a 50 gigabyte variance-covariance matrix (double precision). Moreover, all of these estimators use computations such as equation solutions, inverses, or determinants which employ order of n cubed operations (O(n³)) when operating on the full matrix.

Various ways of attacking this problem have been proposed. For example, Zimmerman (1989) devised a computational acceleration for the Gaussian covariance structure over a regular parallelogram lattice. Sone and Griffith (1995) discussed trade-offs from approximating the difficult normalizing constant (determinant) in maximum likelihood estimation of spatial statistical models. Li (1995) used a supercomputer to accelerate computations of the determinant term.

Ideally, one would like to quickly compute maximum likelihood spatial statistics for any sample size. To achieve this end, we assume that the direct effect of nearby observations decays to 0 after some distance. This leads to a sparse covariance or inverse covariance matrix, which can potentially aid computations (e.g., Ho and Klotz, 1992). In fact, Pace and Barry (1997) illustrated how to accelerate computations in SAR models via sparse matrices. As an illustration, they computed a simultaneous autoregression (SAR) on 20,640 observations using sparse matrices in under 19 minutes despite evaluating a determinant of 20,640 by 20,640 matrix ten times. Furthermore, Barry and Pace (1997) showed how to apply sparse matrix methods to the kriging problem via a mining data set. In addition, they performed a number of timing experiments showing the potential computational gains for plausible values of the spherical variogram.

In the realm of lattice models, many prefer conditional autoregressions (Besag 1974, 1975), Cressie 1993, p. 407-410). Accordingly, this paper provides a number of computational accelerations applicable to conditional autoregressions (CARs).

The sparse matrix technology greatly accelerates computation of the normalizing constant, the major bottleneck in past algorithms. The acceleration of the difficult normalizing constant highlights other bottlenecks in the computation of CARs. We address these other bottlenecks by modifying the CAR computations. First, since many evaluations of the determinant term occur in the computation of each estimate, we reuse the same determinants by evaluating a set of these prior to estimation. Second, we greatly simplify the sum-of-squared errors (SSE) expression within the likelihood. This substantially reduces the number of computations needed. Third, we vectorize the computations. Many computer architectures favor the use of vectorized code. In addition, many languages such as Fortran 90, Gauss, Matlab, and S-Plus run faster with vectorized code.

Moreover, to avoid problems of multiple local optima reported by Warnes and Ripley (1987), Ripley (1988), and Mardia and Watkins (1989), we evaluate the entire profile likelihood for the spatial differencing parameter. This ensures robustness and aids inference. In fact, we show computing profile likelihoods and simulating the estimator provides a low-cost means of inference superior to the traditional information matrix approach. Finally, we demonstrate an interesting invariance of CARs to a change in the magnitude, but not the direction, of the true errors.

As an illustration of the efficacy of these techniques, we perform a Monte Carlo simulation of the CAR estimator using 2,905 observations. Normally, even one spatial autoregression of this size could prove challenging. For example, Li (1995, p. 130) took 8566.5 seconds on an IBM RS6000 to compute a single 2,500 observation SAR. In contrast, we can compute 250 CARs on 2905 observations in 135.78 seconds (43.28 for 100 determinants and 92.5 for the other computations) using a 133 megahertz Pentium.

Section 2 discusses various facets of CAR computations, section 3 details the Monte Carlo study, and section 4 concludes with the key results.

Section A begins by discussing the CAR model, section B provides details on the construction of the crucial spatial weighting matrix, section C discusses the SSE used in maximum likelihood, section D demonstrates the invariance of the estimate to the magnitude of the errors (but not their direction), section E discusses vectorized ways of computing the sum-of-squares, section F details the advantages of precomputing determinants, section G presents the profile likelihood for an entire simulation, and section H gives ways of conducting inference without computing n by n information matrices.

Suppose the errors arise out of the following spatial conditional autoregressive error process,

As an illustration of how to construct C, compare the distance d_ij between every pair of observations j and i to , the distance from observation i and its mth nearest neighbor. Assign a weight of 1 to all observations in the non-normalized matrix whenever d_ij is less than or equal to as in ,

.

Naturally, this yields a weight of 1 for the observation itself (d_ij= d_ii =0) and 0 for each observation j more than distance from observation i. To prevent the observation from predicting itself, we set . The combination of these two procedures allows for observations j colocated with observation i to have a weight of 1.

It seems reasonable to set to 0 the direct influence of distant observations upon a particular observation. This assumption can lead to rather sparse forms for C which will greatly accelerate the computation of the maximum likelihood estimator, as explored in 3.C.

As it stands, is neither symmetrical nor normalized. Let S represent an n by n diagonal matrix containing the square roots of the reciprocals of the sum of the columns and rows of .

To normalize and make C symmetrical, we pre and post-multiply by S.

This operation will ensure the largest eigenvalue of C equals 1. See Ord (1975) for more details.

The conventional approach to computing CAR estimates relies upon solving the CAR normal equations . However, the solution to the equations depends upon the unobservable parameter .

Due to its favorable statistical properties, maximum likelihood (ML) is the preferred way to determine the optimal value of . Substituting the solution of the normal equations into the likelihood function converts it into the univariate profile likelihood of the unknown parameter, . The optimal value of determines the ML estimate, .

Following the conventional approach becomes somewhat tedious when computing the estimate for different realizations of the dependent variable, Y, or when one desires the entire likelihood as a function of the unknown parameter, , also known as the profile likelihood. In either situation, one may need to solve equation many thousands of times.

To increase the computational performance of the CAR estimator, reparameterize it to orthogonalize the independent variables as in . Hence, let Z represent a set of orthogonal independent variables.

One can write the CAR normal equations in as,

.

Substitute the spectral decomposition of in , where U represents a k by k matrix of eigenvectors and represents the k by k diagonal matrix of eigenvalues. Since the eigenvectors are orthogonal, . Neither the spectral decomposition of or the orthogonal reparameterization of cost much computationally, since they work with a k (as opposed to n) dimensional problem. Substituting both of these relations yields,

.

Solving for yields,

Hence, the raw residuals from the above fit become,

One can form the SSE for ML in as in Cressie (1993, p. 437),

Substituting into yields,

Expanding gives,

Reorganizing highlights the term in brackets within ,

Recall . Because of the orthogonality of U, a number of terms cancel, thus simplifying greatly.

Take the expression for in and factor out of the main body of . Hence,

Assume, as CAR does, that,

Substituting into the yields,

The terms involving vanish, thus leaving,

Hence, scaling the errors by a positive constant will not affect where attains its minimum or where the likelihood attains its maximum.

Expanding the basic equation in yields,

One has the choice of using a vector of values of but using only one vector of Y each iteration or using a matrix of values of Y and looping over . In an estimation context with a particular Y, the first approach seems more natural. In a simulation context with multiple realizations of Y, the second approach seems more natural and thus we follow this approach below.

For use with multiple realizations of Y as in simulations redefine Y_r as [Y₁, Y₂,... Y_iter], an n by iter matrix. Define,

where the symbol denotes the use of Hadamard element by element multiplication.

Let take on values . Loop over l iterations of to define .

The log-determinant, , plays an important part in the likelihood function. If takes on only the values , we can precompute an l by 1 vector of log-determinants. Several advantages accrue to computing these as a group. First, the ordering algorithms for sparse matrices apply to all matrices with a given pattern of sparsity. This pattern of zeros and non-zeros for would apply for all positive . Hence, computing the determinants as a group amortizes over l iterations the cost of computing the optimal ordering. Second, as this could prove a time consuming step for some problems, computing them as a group allows the computations to occur during non-peak usage. Third, most fitting exercises involve exploratory or diagnostic variations in the model which would affect the specification of X and hence Z but not affect . Thus, changes in the specification of X further amortize the cost of precomputing .

Given the various computations, we can define twice the profile log-likelihood over the entire simulated Y_r.

For each column of twice the profile-likelihood, find the value of which maximizes each column, . Hence, represents the maximum likelihood estimate, , for the dependent variable associated with that column. Each row provides information for inference on across different realizations Y. Hence, the matrix of profile likelihoods across and Y provides a wealth of information.

Computation of the observed or expected information matrix becomes expensive in a spatial context. Inspection of the information matrix in Cressie (1993, p. 484) shows it requires, among other computations, an n by n inverse (O(n³)) and multiplication of n by n matrices (O(n³)). The inversion of the n by n matrix defeats CAR’s computational advantage of modeling the inverse of the variance-covariance matrix as opposed to kriging’s approach of modeling the variance-covariance matrix directly.

This section examines a simulation of 2,50 regressions each using 2,905 observations. Section A discusses the data, section B defines the role of sparsity, section C contains the timings of the simulation computations, section D discusses the statistical results from the simulation, while section E illustrates the use of simulation for inference.

To provide verisimilitude to the simulation, we chose an actual set of locations for use in forming C, the spatial weighting matrix. Specifically, we used the geographic centroids from all the census block groups in Connecticut from the 1990 Census. This yielded a matrix C with 2,905 rows and 2,905 columns.

In the simulation, we:

We subsequently generated the autocorrelated dependent variable, Y, according to ,

where u represents an n by iter matrix of N(0,1) random variates. In actuality, we solved the corresponding equation system in for the coefficients A_S rather than computing the inverse as this goes much faster.

Compare this to the usual Cholesky decomposition and inverse formulation,

Relative to , requires solving a larger system (n by n instead of n by iter and subsequently multiplying an n by n matrix by a n by iter matrix (O(n²iter)). Since iter usually is much smaller than n, the inverse method takes substantially longer for the same results.

Since efficient computation of the determinant required by maximum likelihood estimation uses the Cholesky decomposition of , this reduces the additional cost of simulating the random numbers. However, solving the 2,905 equations by 2,905 unknowns for iter right hand sides would prove quite difficult without resorting to sparse matrix techniques as described in 3.C below.

If differencing an observation with its nearby neighbors removes most of the effects of autocorrelation, the spatial weighting matrix C can be quite sparse. For example, if an observation displays error dependency only with its nearest m neighbors, only m non-zero entries exist per row of C. Thus, C will contain nm non-zero elements out of n² possible ones. This produces a m/n proportion of non-zero elements, a popular measure of sparsity. For example, with this problem we used four neighbors for each observation in computing . This resulted in C having on average 5.08 neighbors and a sparsity of 5.08/2905 (0.21%). This represents a very high level of sparsity which grows as n increases.

Sparsity results in a number of computational gains. First, it dramatically decreases the storage needed for C and . Using traditional dense techniques, C requires 67.5 MB of storage (double precision). Using sparse matrix techniques, C requires less than 200 KB of storage. Naturally, this divergence grows with n.

Second, sparsity greatly accelerates computations. For example, multiplying the n by n matrix C with the n by k matrix X requires O(kn²) operations using dense matrices. Barring computational bookkeeping, the equivalent sparse operation requires O(knm) operations, a much smaller number. The real benefits come when computing determinants, inverses, or solving systems of equations. All of these operations can build upon the Cholesky decomposition of a matrix and all require O(n³) operations when using dense matrices. Sparsity, however, can totally change the order of the number of operations required in these computations. For example, if had a band structure with half-bandwidth p, the Cholesky decomposition of , would require O() operations (Golub and Van Loan (1989, p. 154)). Hence, for fixed bandwidths the computations grow linearly with n, the number of observations.

Unfortunately, the existence of a pure band structure does not arise very often. Figure 1a shows the actual plot of the non-zero elements in . The existence of such dispersed off-diagonals could make it difficult to achieve computational gains.

However, one can permute the rows or columns of to reduce bandwidth or to achieve other optimizations. A variety of such permutations exist (see George and Liu (1981) for more on the various orderings). For example, the reverse Cuthill-McKee algorithm attempts to reorder the rows and columns of the matrix to create a variable band matrix as shown in Figure 1b. Figure 1b makes the gains of exploiting sparsity obvious. Less obviously, Figure 1c shows the plot of permuted using the column minimum degree algorithm while Figure 1d shows the plot of permuted using the symmetric minimum degree algorithm.

Table I shows the timings associated with computing the log-determinants using the original, random, reverse Cuthill-McKee, column minimum degree orderings, and symmetric column minimum degree orderings for . As Table I makes clear, the ordering of the rows and columns matters, with the symmetric minimum degree ordering reducing execution times by 86% over the original ordering. Intriguingly, the more intuitive reverse Cuthill-McKee ordering actually performed worse than the original ordering but much better than the random ordering. As a worst case scenario, the random ordering produced computational times worse than the optimal ordering by a factor of 572. All computations used the Matlab language running on a 133Mhz Pentium computer.

To place these results in perspective, Li (1995) took the eigenvalue route to computing determinants. Li used an IBM RS6000 Model 550 and a CM5 parallel processing supercomputer. The CM5 had 32 processors each with 32MB of local memory and four vector units. For a 2500 by 2500 spatial weight matrix the RS6000 required 8515.07 seconds while the CM5 required 45.78 seconds. In contrast, it took only 43.28 seconds to compute 100 different determinants of a 2,905 by 2,905 matrix on a 133Mhz Pentium computer using Matlab. Hence, the sparse technology employed here allowed a personal computer on a larger problem to exceed the performance of a supercomputer on a smaller problem!

Finally, even supercomputers have memory limitations. The use of sparse matrix technology has allowed us to handle problems with 20,640 observations (Pace and Barry (forthcoming)). A dense spatial weight matrix would have required 3.58 gigabytes (double precision), which would have further exacerbated the overall computational difficulties.

In a spatial context, maximum likelihood techniques have a reputation of being slow to compute, prone to error through local optima, and offering only asymptotic inference. The approach adopted here to computing conditional autoregressions attacks each of these problems. First, the use of sparse matrices and restructuring the sum-of-squared errors function greatly accelerates computations. Second, the evaluation of the likelihood over a grid of values of the important spatial differencing parameter greatly enhances robustness to local optima. Third, the quick simulation of the conditional autoregression coupled with the profile likelihood information provides a wealth of information for conducting inference in all sample sizes.

As an illustration of the efficacy of these techniques, we computed (a) 250 simulated conditionally autoregressive vectors of Y of 2905 observations; (b) 100 determinants of a 2,905 by 2,905 matrix; and (c) computed 250 estimates with the profile likelihood values. This took 135.8 seconds (43.3 seconds for the determinants and 92.5 for the estimates). In contrast, Li (1995) on a smaller problem took longer on a supercomputer to calculate the determinants than we accomplished on a personal computer.

Hopefully, additions like these to the spatial statistics toolbox will enable users to help cope with the increasing flow of spatial information.

We both would like to thank the University of Alaska for its generous research support. In addition, Pace would like to acknowledge support from the Center for Real Estate and Urban Economic Studies, University of Connecticut. Finally, we would like the thank Bill Hughes for having discussed the paper at the American Real Estate Society meetings.

Barry, Ronald P., and R. Kelley Pace. (1997). "Kriging with Large Data Sets Using Sparse Matrix Techniques," Communications in Statistics: Computation and Simulation, 26.

Barry, Ronald P., and Jay M. Ver Hoef. (1997). "Blackbox Kriging — Kriging without Specifying Variogram Models," Journal of Agricultural, Biological, and Environmental Statistics.

Besag, Julian E. (1974). "Spatial Interaction and the Statistical Analysis of Lattice Systems," Journal of the Royal Statistical Society, B 36, 192-236.

Besag, Julian E. (1975). "Statistical Analysis of Non-lattice Data," The Statistician, 24, 179-95.

Cressie, Noel A.C. (1993). Statistics for Spatial Data, Revised ed. New York: John Wiley.

George, Alan, and Joseph Liu. (1981). Computer Solution of Large Sparse Positive Definite Systems, Englewood Cliffs: Prentice-Hall.

Gilley, O.W., and R. Kelley Pace. (1995). "Using Inequality Restrictions to Tame Multicollinearity in Hedonic Pricing Models," Review of Economics and Statistics, 77, 609-621.

Golub, G. H., and C. F. Van Loan. (1989). Matrix Computations, second edition, John Hopkins.

Haining, Robert. (1990). Spatial Data Analysis in the Social and Environmental Sciences, Cambridge University Press.

Ho, Shu-Yen, and Jerome Klotz. (1992). "Sparse Matrix Methods for Unbalanced Multifactor Analysis of Variance and Covariance," Journal of Statistical Computation and Simulation, 41, 55-72.

Li, Bin. (1995). "Implementing Spatial Statistics on Parallel Computers," in: Arlinghaus, S., ed. Practical Handbook of Spatial Statistics, (CRC Press, Boca Raton), 107-148.

Mardia, K.V., and A.J. Watkins. (1989). "On Multimodality of the Likelihood in the Spatial Linear Model," Biometrika, 76, 289-295.

Meeker, W.Q., and L.A. Escobar. (1995). "Teaching about Approximate Confidence Regions Based on Maximum Likelihood Estimates," The American Statistician, 49, 48-53.

Ord, J.K. (1975). Estimation Methods for Models of Spatial Interaction, Journal of the American Statistical Association, 70, 120-126.

Pace, R. Kelley, and Ronald Barry. (1997). "Sparse Spatial Autoregressions," Statistics and Probability Letters, 33, 291-297.

Pace, R. Kelley, and O.W. Gilley. (1993). "Improving Prediction and Assessing Specification Quality in Non-Linear Statistical Valuation Models," Journal of Business and Economics Statistics, 11, 301-310.

Ripley, Brian D. (1981). Spatial Statistics, New York: John Wiley.

Ripley, Brian D., (1988). Statistical Inference for Spatial Processes, Cambridge: Cambridge University Press.

Rubinstein, R.V. (1981). Simulation and the Monte Carlo Method, New York: Wiley.

Sone, Akio, and Daniel Griffith. (1995). "Trade-Offs Associated with Normalizing Constant Computational Simplifications for Estimating Spatial Statistical Models," Journal of Statistical Computation and Simulation, 51, 165-184.

Warnes, J.J., and Brian Ripley. (1987). "Problems with Likelihood Estimation of Covariance Functions of Spatial Gaussian Processes," Biometrika, 74, 640-642.

Zimmerman, Dale. (1989). "Computationally Exploitable Structure of Covariance Matrices and Generalized Covariance Matrices in Spatial Models," Journal of Statistical Computation and Simulation, 32, 1-15.