Dear Doug, My apologies for misstating something. The variance-covariance matrix for the errors in these spatial models is positive definite, not the (In - rho*W) term.
The W matrix is not symmetric in general. My goal is to create a "lookup table" of values of log(det(In-rho*W)) for values of rho from lmin to lmax and then simply look these values up during the Metropolis step. This is why I asked if the code could be used in the case of a non-symmetric W matrix. Thank you again for taking the time to look at this. Don On Friday, August 1, 2014 3:20:11 PM UTC-4, Douglas Bates wrote: > > On Thursday, July 31, 2014 8:33:49 PM UTC-5, Donald Lacombe wrote: >> >> Dear Doug, >> >> Thank you for taking the time to look at this. In actuality, the W matrix >> is a spatial weight matrix which defines who is a neighbor to whom. I think >> the "symmetric" verbiage in the code is not technically correct. The >> "standardized" means that the rows of the W matrix sum to 1. This >> ensures that the (In -rho*W) matrix is positive definite as you state. >> > > I'm confused about the terminology then. To me it doesn't make sense to > say that (I-rho*W) is positive definite unless W is symmetric (or > Hermitian, for complex arithmetic) as stated in > http://en.wikipedia.org/wiki/Positive-definite_matrix. > > >> >> The following matrix would not technically be a proper spatial weight >> matrix because the 4th row has 1 neighbor and the 5th row has three >> neighbors when in fact they should all have 2 neighbors (actually there can >> be a different number of neighbors but I'm choosing 2 just for an >> illustration). >> >> julia> full(W) >> 5x5 Array{Float64,2}: >> 0.0 0.5 0.5 0.0 0.0 >> 0.5 0.0 0.0 0.0 0.5 >> 0.5 0.0 0.0 0.0 0.5 >> 0.0 0.0 0.0 0.0 0.5 >> 0.0 0.5 0.5 0.5 0.0 >> >> The original matrix I posted was a proper nearest neighbor spatial weight >> matrix: >> >> >> full(W) >> >> ans = >> >> 0 0.5000 0.5000 0 0 >> 0.5000 0 0 0 0.5000 >> 0.5000 0 0 0 0.5000 >> 0 0.5000 0 0 0.5000 >> 0 0 0.5000 0.5000 0 >> >> Each row represent a geographic entity and the columns represent their >> neighbors. The diagonal elements are zero because no geographic entity is a >> neighbor to itself. The 0.5's represent the row-normalization and in this >> case the value is .5 because there are 2 neighbors. >> >> I guess the bottom line is that the code you provided works but I would >> need it to work for a non-symmetric W matrix. >> >> Can the above code be altered to deal with the non-symmetric W case? >> >> Thank you again for taking the time to look at this. I truly appreciate >> the help and education. >> >> Regards, >> Don >> >> >> On Thursday, July 31, 2014 4:49:21 PM UTC-4, Douglas Bates wrote: >>> >>> Rats, I had a reply with code and examples then Google groups croaked >>> and when it reloaded my draft reply was gone. >>> >>> On Thursday, July 31, 2014 2:13:36 PM UTC-5, Donald Lacombe wrote: >>>> >>>> Dear Doug, >>>> >>>> I've created a small example to show what the MATLAB code is actually >>>> doing. Here is the output with some comments: >>>> >>>> % Create random coordinates >>>> >> xc = randn(5,1); >>>> >> yc = randn(5,1); >>>> >>>> % Create 2 nearest neighbor weight matrix >>>> >> W = make_neighborsw(xc,yc,2); >>>> >> W >>>> >>>> W = >>>> >>>> (2,1) 0.5000 >>>> (3,1) 0.5000 >>>> (1,2) 0.5000 >>>> (4,2) 0.5000 >>>> (1,3) 0.5000 >>>> (5,3) 0.5000 >>>> (5,4) 0.5000 >>>> (2,5) 0.5000 >>>> (3,5) 0.5000 >>>> (4,5) 0.5000 >>>> >>>> % Full version >>>> >> full(W) >>>> >>>> ans = >>>> >>>> 0 0.5000 0.5000 0 0 >>>> 0.5000 0 0 0 0.5000 >>>> 0.5000 0 0 0 0.5000 >>>> 0 0.5000 0 0 0.5000 >>>> 0 0 0.5000 0.5000 0 >>>> >>> >>> Notice that this matrix is not symmetric but the description of the >>> arguments for indetfull below says that W should be symmetric. I am going >>> to assume that W should be symmetric and the allowable values of rho are >>> such that I - rho*W is positive-definite. I will use the matrix >>> >>> julia> W = sparse([2,3,1,5,1,5,5,2,3,4],[1,1,2,2,3,3,4,5,5,5],0.5) >>> 5x5 sparse matrix with 10 Float64 entries: >>> [2, 1] = 0.5 >>> [3, 1] = 0.5 >>> [1, 2] = 0.5 >>> [5, 2] = 0.5 >>> [1, 3] = 0.5 >>> [5, 3] = 0.5 >>> [5, 4] = 0.5 >>> [2, 5] = 0.5 >>> [3, 5] = 0.5 >>> [4, 5] = 0.5 >>> >>> julia> issym(W) >>> true >>> >>> julia> full(W) >>> 5x5 Array{Float64,2}: >>> 0.0 0.5 0.5 0.0 0.0 >>> 0.5 0.0 0.0 0.0 0.5 >>> 0.5 0.0 0.0 0.0 0.5 >>> 0.0 0.0 0.0 0.0 0.5 >>> 0.0 0.5 0.5 0.5 0.0 >>> >>> >>> for illustration. >>> >>> >>>> % Find colamd >>>> >> z = speye(n) - 0.1*sparse(W) >>>> >>>> z = >>>> >>>> (1,1) 1.0000 >>>> (2,1) -0.0500 >>>> (3,1) -0.0500 >>>> (1,2) -0.0500 >>>> (2,2) 1.0000 >>>> (4,2) -0.0500 >>>> (1,3) -0.0500 >>>> (3,3) 1.0000 >>>> (5,3) -0.0500 >>>> (4,4) 1.0000 >>>> (5,4) -0.0500 >>>> (2,5) -0.0500 >>>> (3,5) -0.0500 >>>> (4,5) -0.0500 >>>> (5,5) 1.0000 >>>> >>>> % Full z >>>> >> full(z) >>>> >>>> ans = >>>> >>>> 1.0000 -0.0500 -0.0500 0 0 >>>> -0.0500 1.0000 0 0 -0.0500 >>>> -0.0500 0 1.0000 0 -0.0500 >>>> 0 -0.0500 0 1.0000 -0.0500 >>>> 0 0 -0.0500 -0.0500 1.0000 >>>> >>>> >> p = colamd(z) >>>> >>>> p = >>>> >>>> 4 3 1 2 5 >>>> >>>> What I'm looking for is the "p" vector above. The full MATLAB code for >>>> the routine is as follows: >>>> >>> >>> You don't need colamd for the calculations in indetfull. Julia has a >>> logdet function with a method for the sparse Cholesky factorization so the >>> core calculation becomes >>> >>> julia> rho = 0.1 >>> 0.1 >>> >>> julia> cf = cholfact(speye(W) - rho*W); >>> >>> julia> logdet(cf) >>> -0.012566124059126464 >>> >>> You don't need it but for interest, the 0-based permutation is stored as >>> the .Perm member of the cf object. To make it 1-based we add 1 to it >>> >>> julia> show(cf.Perm .+ 1) >>> [4,5,1,3,2] >>> >>> As with most sparse matrix operations the Cholesky factorization is done >>> as a symbolic phase followed by a numeric phase. If you are just modfiying >>> rho you don't need to do the symbolic phase over and over. You can save >>> the factor from the first one and do the subsequent factorizations in place. >>> >>> function indetfull(W::SparseMatrixCSC,lmin,lmax) >>> issym(W) || error("Sparse matrix W must be symmetric") >>> rhos = lmin:0.01:lmax >>> I = speye(W) >>> cf = cholfact(I - lmin*W) >>> hcat(rhos, [logdet(cholfact!(cf,I-rho*W))::Float64 for rho in rhos]) >>> end >>> >>> (By the way, be careful not to use lmin = 0. or else change the code so >>> that cf is defined using lmax instead of lmin. The pattern of I - 0.*W is >>> the pattern of the identity, after you have dropped the zero values, and is >>> different from the pattern of I - rho*W for non-zero rho. It is important >>> that the pattern be the same for all values of rho.) >>> >>> Here is an example of its use >>> >>> julia> indetfull(W,0.01,0.1) >>> 10x2 Array{Float64,2}: >>> 0.01 -0.000125007 >>> 0.02 -0.000500105 >>> 0.03 -0.00112553 >>> 0.04 -0.00200168 >>> 0.05 -0.00312911 >>> 0.06 -0.00450853 >>> 0.07 -0.00614082 >>> 0.08 -0.00802701 >>> 0.09 -0.0101683 >>> 0.1 -0.0125661 >>> >>> >>> >>> >>> >>>> function out=lndetfull(W,lmin,lmax) >>>> % PURPOSE: computes Pace and Barry's grid for log det(I-rho*W) using >>>> sparse matrices >>>> % >>>> ----------------------------------------------------------------------- >>>> % USAGE: out = lndetfull(W,lmin,lmax) >>>> % where: >>>> % W = symmetric spatial weight matrix (standardized) >>>> % lmin = lower bound on rho >>>> % lmax = upper bound on rho >>>> % >>>> ----------------------------------------------------------------------- >>>> % RETURNS: out = a structure variable >>>> % out.lndet = a vector of log determinants for 0 < rho < 1 >>>> % out.rho = a vector of rho values associated with lndet >>>> values >>>> % >>>> ----------------------------------------------------------------------- >>>> % NOTES: should use 1/lambda(max) to 1/lambda(min) for all possible rho >>>> values >>>> % >>>> ----------------------------------------------------------------------- >>>> % References: % R. Kelley Pace and Ronald Barry. 1997. ``Quick >>>> % Computation of Spatial Autoregressive Estimators'', Geographical >>>> Analysis >>>> % >>>> ----------------------------------------------------------------------- >>>> >>>> % written by: >>>> % James P. LeSage, Dept of Economics >>>> % University of Toledo >>>> % 2801 W. Bancroft St, >>>> % Toledo, OH 43606 >>>> >>>> >>>> rvec = lmin:.01:lmax; >>>> spparms('tight'); >>>> [n junk] = size(W); >>>> z = speye(n) - 0.1*sparse(W); >>>> p = colamd(z); >>>> niter = length(rvec); >>>> dettmp = zeros(niter,2); >>>> for i=1:niter; >>>> rho = rvec(i); >>>> z = speye(n) - rho*sparse(W); >>>> [l,u] = lu(z(:,p)); >>>> dettmp(i,1) = rho; >>>> dettmp(i,2) = sum(log(abs(diag(u)))); >>>> end; >>>> >>>> out.lndet = dettmp(:,2); >>>> out.rho = dettmp(:,1); >>>> >>>> The code will be very easy to convert once I figure out this one >>>> command, >>>> >>>> Thank you in advance for any help that you are able or willing to >>>> provide. I'm not an expert as sparse matrices but these routines come in >>>> very handy for these types of calculations. >>>> >>>> Thanks, >>>> Don >>>> >>>> >>>> On Thursday, July 31, 2014 10:42:55 AM UTC-4, Douglas Bates wrote: >>>>> >>>>> On Wednesday, July 30, 2014 8:10:37 PM UTC-5, Donald Lacombe wrote: >>>>>> >>>>>> Dear Julia Users: >>>>>> >>>>>> I am attempting to speed up some code regarding some Bayesian spatial >>>>>> econometric models and would like to translate some MATLAB code for >>>>>> computing the log determinant term in the log-likelihood function. >>>>>> Specifically, I'm interested in the Pace and Barry approximation to the >>>>>> log-determinant term as in R. Kelley Pace and Ronald Barry. (1997) >>>>>> "Quick Computation >>>>>> of Spatial Autoregressive Estimators''. The term in question is det(In - >>>>>> rho*W) which has to be calculated through each iteration of the >>>>>> Metropolis-Hastings algorithm I'm using. >>>>>> >>>>>> The MATLAB code for this has the command "colamd" which MATLAB's help >>>>>> defines as "P = colamd(S) returns the column approximate minimum >>>>>> degree permutation vector for the sparse matrix S". >>>>>> >>>>> >>>>> colamd is indirectly available in the Cholmod module but I don't think >>>>> there is a direct way of calling it. It is a method for determining a >>>>> fill-reducing permutation for a sparse, symmetric matrix of the form X'X >>>>> directly from X. There are several other approaches implemented in amd, >>>>> Metis, Scotch, MUMPS, etc. for dealing with the symmetric matrix itself. >>>>> All of these just use the pattern of nonzeros in the matrix as they are >>>>> part of the symbolic stage of the sparse Cholesky or sparse QR >>>>> factorization. >>>>> >>>>> Can you be more specific about how it is to be used? In particular, >>>>> is it important to work with the pattern of X rather than the pattern of >>>>> X'X? >>>>> >>>>> >>>>>> Is there an equivalent Julia command? I've tried looking at "symperm" >>>>>> but I do not think I am using it correctly and more importantly I have >>>>>> no >>>>>> real clue. >>>>>> >>>>>> Any help would be appreciated and if I successfully code this I >>>>>> would be happy to share with anyone who is interested. >>>>>> >>>>>> Thanks, >>>>>> Don >>>>>> >>>>>
