It is well known (arXiv) in continuous variable (CV) quantum information processing that an $N$-mode CV system whose dynamics are entirely Gaussian can be simulated efficiently with a classical computer as the dynamics can be simulated by computing how the desired Gaussian unitary operations transform the $2N \times 2N$ covariance matrix and $2N$-dimensional vector of means (first moments) that contain the covariances and expectations of the field quadratures of each mode in the system, respectively. These Gaussian unitaries are those that are generated by Hamiltonians which consist of terms that are at most quadratic in the quadrature field operators $\hat{q}$ and $\hat{p}$. For example
$$ \hat{H} = \sum_{i=1}^N{\hat{q}_i^2+\chi_i \hat{p}_{i-1}\hat{p}_i + w_i\hat{p}_i}. $$ However I have not been able to find information on how to efficiently compute the covariance matrix and vector of means of the ground state (or any eigenstate for that matter) of an arbitrary Gaussian Hamiltonian. I would imagine this should also be able to be computed efficiently with a classical computer. Is anyone aware of techniques on how to do this?
This question was also asked here.