Add low-rank multivariate normal distribution

Description

If a low-rank structure exists (or is assumed to exist) in the covariance matrix of MVN, i.e., $$\Sigma=W W^\top + D$$, where $$W$$ is $$n\times m$$ ($$m<<n$$), $$D$$ is a positive definite diagonal matrix, then the determinant and the inverse can be reduced from $$O(n^3)$$ to $$O(nm^2+m^3)$$. A continuation of #3038.

For implementation, can refer to
https://en.wikipedia.org/wiki/Woodbury_matrix_identity
https://en.wikipedia.org/wiki/Matrix_determinant_lemma

numpyro
Low-rank Multivariate Normal distribution · Issue #387 · pyro-ppl/numpyro
https://github.com/pyro-ppl/numpyro/blob/master/numpyro/distributions/continuous.py

pytorch
Probability distributions - torch.distributions — PyTorch 2.10 documentation
pytorch/torch/distributions/lowrank_multivariate_normal.py at 14e348b7ad1b3472812f2b077020d80deaf6a787 · pytorch/pytorch

Applications:
Intermediate between mean-field and full-rank ADVI
AutoLowRankMultivariateNormal-numpyro

Expected Output

Maybe follow the multi_normal_XXX series?

real multi_normal_lowrank_lpdf(vector y | vector mu, matrix cov_factor, vector cov_diag)
also
real multi_normal_lowrank_lpdf(vectors y | vectors mu, matrix cov_factor, vector cov_diag)

Start by examining the existing multi_normal_XXX distribution series in Stan Math and compare the referenced NumPyro continuous.py and PyTorch lowrank_multivariate_normal.py implementations. Use the Woodbury matrix identity and matrix determinant lemma references to define the low-rank covariance behavior, with completion demonstrated by the two scalar and vectorized multi_normal_lowrank_lpdf signatures described in the issue.