Suppose a user has a Cholesky decomposition of a matrix R and a vector of positive entries σ and wishes to create Σ = PDMat(cholesky(σ .* Matrix(R) .* σ')). This is common in Bayesian inference. e.g. a Cholesky factor of a correlation matrix is drawn from a LKJCholesky distribution, while a vector of standard deviations is drawn from some other distribution, and one wants to construct the Cholesky factorization of the covariance matrix for use in MvNormal.
This can be done much more efficiently with Σ = PDMat(Cholesky(R.uplo == 'U' ? R.U * Diagonal(σ) : Diagonal(σ) * R.L, R.uplo, 0)). But this is still quite long. It would be ideal to have either a constructor for PDMat or some other utility function like the following:
function PDMat(fac::Cholesky, scale::AbstractVector)
factors_scaled = fac.uplo == 'U' ? fac.factors .* scale' : scale .* fac.factors
return PDMat(Cholesky(factors_scaled, fac.uplo, 0))
end
function PDMat(fac::Cholesky, scale::UniformScaling)
return PDMat(Cholesky(fac.factors * scale, fac.uplo, 0))
end
PDMat(fac::Cholesky, scale::Number) = PDMat(fac, scale * I)
Suppose a user has a
Choleskydecomposition of a matrixRand a vector of positive entriesσand wishes to createΣ = PDMat(cholesky(σ .* Matrix(R) .* σ')). This is common in Bayesian inference. e.g. a Cholesky factor of a correlation matrix is drawn from aLKJCholeskydistribution, while a vector of standard deviations is drawn from some other distribution, and one wants to construct the Cholesky factorization of the covariance matrix for use inMvNormal.This can be done much more efficiently with
Σ = PDMat(Cholesky(R.uplo == 'U' ? R.U * Diagonal(σ) : Diagonal(σ) * R.L, R.uplo, 0)). But this is still quite long. It would be ideal to have either a constructor forPDMator some other utility function like the following:See also discussions in https://github.com/JuliaStats/Distributions.jl/issues/1336 and https://github.com/TuringLang/Turing.jl/issues/1629#issuecomment-851959832.