Problem using cov_exp_quad with large matrices solving Log-Gaussian Cox Processes

Thank you for the comment,
then I can try GPflow stuff to get approximate answers and then wait for future release of Stan.

Actually, what I was trying while having computing problems with LGCPs was Sparse Latent Factor Analysis, to reduce the amount of computation.

Rather than ~50k voxels, I set ~500 3D Gaussian Kernels as bases functions and estimate factors as linear combination of bases functions. To reduce the number of factors, I tried to implement sparse-inducing priors for factor-loadings matrix such as as horseshoe priors, but I eventually ended up implementing Dirichlet-Laplace(DL) priors like applied in here.
And I am still having difficulties, may be derived from identifiability I guess…?

Is there any advice on implementing sparse-inducing prior for factor-loadings matrix?
I already opened such questions, but still not answered yet, like this and this.

Thank you for the valuable comments!
You already answered many of my questions that I had struggling with for a long time.
Minho