Using matrix routines, how bad would it be to compute the inverse of $X C_{\rm GP} X^T$ ?



In [9]:

    
ncomp = 5
npix = 2300









    Out[9]:





2300



In [10]:

    
pcomps = randn((npix, ncomp));



In [11]:

    
Cgp = eye(ncomp);



In [12]:

    
Cdata = eye(npix);



In [13]:

    
dar = pcomps * Cgp * transpose(pcomps);



In [14]:

    
#Create some fake residuals
ys = 0.01 * randn((npix,));
#S = create_sparse(wl, 1., 100.);



In [15]:

    
function lnprob()
    L = cholfact(dar + Cdata)
    out = L \ ys;
    lnp = -0.5 * (logdet(L) + dot(ys, out[:,1]) + npix/2. * log(2pi))
end









    Out[15]:





lnprob (generic function with 1 method)



In [16]:

    
@time lnprob()









    



elapsed time: 0.566725692 seconds (84729816 bytes allocated, 5.18% gc time)






    Out[16]:





-1076.2286122417768



In [ ]: