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