I want to estimate parameters using the simulated likelihood approach in a particular context that I describe below.
Let $G_1,\dots,G_{S}$, where $S=20000$, be matrices of size $N\times N$ with $N=10000$, $\mathbf{Y} \in \mathbb{R}^{N}$, $ \mathbf{X} \in \mathcal{M}_{N,K}$ a matrices of size $N\times K$. The simulated likelihood is defined by: $$\hat{f}(\alpha,\beta,\sigma^2)=\frac{1}{S}\sum_{s=1}^S\dfrac{|A_s|}{\sqrt{\sigma^{2N}}}\exp{\big\{-\dfrac{1}{2\sigma^2}(A_s Y - X\beta)'(A_s Y - X\beta)\big\}}$$ $$\text{where}\quad A_s=I - \alpha G_s$$
Theoretically, it seems simple. But the problem is that, I do not have the 20000 matrices $G_s$ simultaneously. They are computed by another process and I can not save them simultaneously on my computer. I will need 15 TB to do that. Is there a way to estimate $\alpha$, $\beta$ and $\sigma^2$ without using simultaneously the matrices?
Perhaps an iterative approach that only require one $G_s$ at each iteration. I can compute even 8 or 10 matrix and save in memory but not 20000 matrices.
Thank!
Finally I found a solution. I post it hoping that It will be useful later. Instead of saving the matrices $G_s$, I only saved what is needed to compute the simulated likelihood. In fact, we normally need $A_s Y$ and $|A_s|$ for $s=1,2,\dots,S$. However, these objects depends on $\alpha$ which is unknown.
We knows that $A_sY = Y - \alpha G_s Y$. Thus, to compute this expression, I only need to save $G_s Y$. One $G_s Y$ contains $N$ entries and costs less than 1 MB and the 20000 cost 1.5 GB. As a consequence, It is possible to save them. The most complicated is $|A_s|$. But I do that:
\begin{align*}|A_s| &= |I-\alpha G_s| = (-\alpha)^N\big|G_s-\frac{1}{\alpha}I\big|\\ &=(-\alpha)^NP\Big(\frac{1}{\alpha}\Big) \quad \text{where $P$ is the characteristic polynomial of $G_s$}\\ &=(-\alpha)^N \prod_{k=1}^N\Big(\frac{1}{\alpha}-\lambda_k\Big) \quad \text{where $\lambda_1,\dots,\lambda_N$ are $N$ eigenvalues of $G_s$}\\ &=\prod_{k=1}^N\big(\alpha\lambda_k-1\big)\end{align*}
Therefore, to compute $|A_s|$, I need the $N$ eigenvalues of $G_s$. The problem then becomes tractable. When I compute one $G_s$, I only save in memory $G_s Y$ and the $N$ eigenvalues of $G_s$. This latter is not difficult numerically, even for a large matrix. With this approach, I need 3 GB (instead of 15000 GB) to do the work.