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 >>>>> >>>>
