Andreas, Yes, you are correct and your solution works just fine.
Thank you for taking the time to answer my question. I've said it before but the kindness and generosity of this community is remarkable! Thanks, Don On Tuesday, July 8, 2014 5:34:45 AM UTC-4, Andreas Noack wrote: > > I think I-ρW would be even better in this case. I is a generic identity > matrix. A and B are probably positive definite. If so, I think > logdet(cholfact(A)) would be the best solution. > > > 2014-07-08 11:06 GMT+02:00 Simon Byrne <[email protected] <javascript:>>: > >> 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 >>> >>> >>> >>> > > > -- > Med venlig hilsen > > Andreas Noack Jensen >
