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