Hi Don,
I think the reason they're not sparse is that in the line
A = (In - rho*W);
the matrix In is not sparse: if you replace the line
In = eye(n);
by
In = speye(n);
the result should then be sparse. However at the moment there doesn't seem
to be a sparse det method (I just filed issue #7543
<https://github.com/JuliaLang/julia/issues/7543>): you can get around this
by using log(det(lufact(A))) instead of log(det(A)).
You might also want to have a quick look at the performance tips
<http://docs.julialang.org/en/latest/manual/performance-tips/>, in
particular:
* You should use randn() instead of randn(1): randn() samples a single
Float64 scalar, randn(1) samples a vector of length 1. On my machine (using
0.3 pre-release) your code doesn't work, as this promotes some vectors to
matrices, which doesn't work with dot. You also then won't need to call
variables like rho_proposed[1].
* Try wrapping everything in a function: this makes it easier for the JIT
compiler to work its magic, as well as making it easier to use profiling
tools.
* If you're using dense matrices, you can use logdet(A) instead of
log(det(A)): this can avoid some underflow and overflow problems
(unfortunately this doesn't work for sparse lu factorizations yet, see the
issue mentioned above).
* Purely from a style point of view, there's no need to keep functions in
separate files, or end lines with semi-colons.
Simon
On Tuesday, 8 July 2014 01:48:56 UTC+1, Donald Lacombe wrote:
>
> Dear Julia Users,
>
> I am currently developing/converting several Bayesian spatial econometric
> models from MATLAB into Julia. I have successfully coded the spatial
> autoregressive model (SAR) with diffuse priors in Julia but have a question
> regarding the use of sparse matrix algorithms. The models use a spatial
> weight matrix which is usually sparse in nature and this matrix appears in
> many of the expressions, especially in the random walk Metropolis algorithm
> used to obtain the marginal distribution for the spatial autocorrelation
> parameter rho. Here is a code snippet and the complete code is attached:
>
> # Draw for rho
>
> A = (In - rho*W);
>
> denominator = log(det(A))-.5*sigma^-2*dot(A*y-x*beta,A*y-x*beta);
>
> accept = 0;
>
> rho_proposed = rho + zta*randn(1);
>
> while accept == 0
>
> if ((rho_proposed[1] > -1) & (rho_proposed[1] < 1));
>
> accept = 1;
>
> else
>
> rho_proposed = rho + zta*randn(1);
>
> end
>
> end
>
> B = (In - rho_proposed[1]*W);
>
> numerator = log(det(B))-.5*sigma^-2*dot(B*y-x*beta,B*y-x*beta);
>
> u = rand(1);
>
> if ((numerator[1] - denominator[1]) > exp(1))
>
> pp = 1;
>
> else
>
> ratio = exp(numerator[1] - denominator[1]);
>
> pp = min(1,ratio);
>
> end
>
> if (u[1] < pp[1])
>
> rho = rho_proposed[1];
>
> acc = acc + 1;
>
> end
>
> ar = acc/gibbs;
>
>
> if ar < 0.4
>
> zta = zta/1.1;
>
> end
>
> if ar > 0.6
>
> zta = zta*1.1;
>
> end
>
>
> A = (In - rho*W);
>
>
> Basically, I'm wondering if there is any way to make the "A" and "B" matrices
> sparse and possibly make it run faster, especially in the log(det(A)) terms.
> Currently, 110 draws (with 10 burn) takes approximately 8 seconds on my 64
> bit, core i7 laptop. The computational speed decreases with the sample size n
> because the weight matrices are treated as full.
>
>
> Any help would be greatly appreciated and if anyone is interested in running
> the code, the Distributions and Distance packages must be included and
> initiated first.
>
>
> Regards,
>
> Don
>
>
>
>
