On 0.3.x there is a very expensive error bounds calculation in the triangular solve which the reason for the surprisingly slow calculation. This is not acceptable and we have therefore removed the error bounds calculation in 0.4. On my machine I get
julia> @time L\B; elapsed time: 2.535437796 seconds (8096672 bytes allocated) on 0.3.7 and julia> @time L\B; elapsed time: 0.022492837 seconds (15 MB allocated, 11.99% gc time in 1 pauses with 0 full sweep) on 0.4. 2015-03-16 18:22 GMT-04:00 matt <katzf...@gmail.com>: > Hi, > > Assume that L is the lower Cholesky factor of a symmetric p.d. matrix A, > such that A = LL'. I am trying to calculate L^(-1)B for some matrix B, but > my current way of doing it is extremely slow. > > Here is an example: > using Distances > n=1000 > A=exp(-pairwise(Euclidean(),[1:n]',[1:n]')/n*50) > B=rand(n,n) > L=chol(A,:L) # lower cholesky factor of A > @time L\B # this is what I want, but it is slow (1.2 seconds on my > computer) > > Compare this to the time it takes to solve the following (harder) problem > of computing A^(-1)B=L'^(-1)(L^(-1)B): > Achol=cholfact(A,:L) # cholesky representation of A = LL' > @time Achol\B # this is much faster (.04 seconds), even though it solves > the harder problem L'^(-1)(L^(-1)B) > @time Achol[:L]\B # again, this is what I want, but computation is slow > (1.3 seconds) > > I suspect this discrepancy in timing is because Achol\B uses > back-/forward-solves that exploit the triangular structure of Cholesky > factors, whereas L\B does not. Is there any way I can obtain L^(-1)B much > faster? > > Thank you! >
