I think your direction is very well understood, if the rank is low enough one could apply factorization to speed things up (Cholesky for P being Positive Definite, LDL in case it is only symmetric would be classic choice while KernelMatrices.jl may offer more advanced methods [See KernelMatrices.jl - A Software Package for Working with Hierarchical Matrices]).
What about the case the rank is high to make those method not feasible.
Is there any function which can take advantage of the structure of the problem?