Julia equivalent of MATLAB "colamd"

classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|

Julia equivalent of MATLAB "colamd"

Donald Lacombe
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".

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 
Reply | Threaded
Open this post in threaded view
|

Re: Julia equivalent of MATLAB "colamd"

Douglas Bates
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 
Reply | Threaded
Open this post in threaded view
|

Re: Julia equivalent of MATLAB "colamd"

Donald Lacombe
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 
Reply | Threaded
Open this post in threaded view
|

Re: Julia equivalent of MATLAB "colamd"

Douglas Bates
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 
Reply | Threaded
Open this post in threaded view
|

Re: Julia equivalent of MATLAB "colamd"

Donald Lacombe
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 
Reply | Threaded
Open this post in threaded view
|

Re: Julia equivalent of MATLAB "colamd"

Viral Shah
In reply to this post by Douglas Bates
I think some of this is already in the documentation in bits and pieces, but we should probably add Doug's explanation into the manual. 

Regarding colamd, like Doug said, it is already included, and it should be easy enough to expose it with a few ccalls.

-viral

On Friday, August 1, 2014 2:19:21 AM UTC+5:30, 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 
Reply | Threaded
Open this post in threaded view
|

Re: Julia equivalent of MATLAB "colamd"

Douglas Bates
In reply to this post by Donald Lacombe
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 
Reply | Threaded
Open this post in threaded view
|

Re: Julia equivalent of MATLAB "colamd"

Donald Lacombe
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 <a href="http://en.wikipedia.org/wiki/Positive-definite_matrix" target="_blank" onmousedown="this.href='http://www.google.com/url?q\75http%3A%2F%2Fen.wikipedia.org%2Fwiki%2FPositive-definite_matrix\46sa\75D\46sntz\0751\46usg\75AFQjCNFpBxIroI1141Y1bviIgSWH-nYIXA';return true;" onclick="this.href='http://www.google.com/url?q\75http%3A%2F%2Fen.wikipedia.org%2Fwiki%2FPositive-definite_matrix\46sa\75D\46sntz\0751\46usg\75AFQjCNFpBxIroI1141Y1bviIgSWH-nYIXA';return true;">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