This really hammers you for high dimensional settings (e.g. brain imaging n>100k).
Solutions:
* Pick a very special kernel with an known / approximable solution.
* Reduce the dimension of the problem (sparse gps).
Chapter 8 of the Rasmussen reference does talk about approximation methods for large datasets, but I couldn't find any reference implementations to go off of.
For brain imaging, 100k was 'the dimensionality of the input' (but also the number of datapoints). This could easily be O(10^6) though with higher res imaging.
I'd check out, for discussion:
https://www.prowler.io/blog/sparse-gps-approximate-the-poste...
For code:
https://github.com/GPflow/GPflow
Based on the people involved knowing what they are talking about, rather than on experience using that particular work.