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