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