Bilinear Model in Turing

Specific Domains Statistics
1

Hello,

I am new to Julia and I am running Bayesian inference in Turing for a Bilinear Model-

Y = X * B * Z’ + E

where Y is nxm, X is nxp, B is pxq and Z is mxq. E is also nxm and follows multivariate normal (0, Sigma). The object of inference is B and the var-cov matrix Sigma. I am putting a conjugate multivariate normal prior on B and an inverse-wishart prior on Sigma (see code below) but for a problem with high dimension (say m = 250), I end up getting the following error.

PosDefException: matrix is not positive definite; Factorization failed.

Since I can derive the conditional distributions by hand, is there a way I can directly sample from the conditional distributions using Turing?

My code is pasted below:

@model function bilinear_model(X, Z, Y, tau2, Psi0, nu0)
    # dimensions
    n, p = size(X)
    m, q = size(Z)

    # Prior on Sigma: Inverse Wishart
    Sigma ~ InverseWishart(nu0, PDMat(Symmetric(Psi0 + 1e-8*I(m))))

    # Prior on vec(B): Multivariate Normal
    pq = p * q
    beta ~ MvNormal(zeros(pq), tau2 * I(pq))
    B = reshape(beta, p, q)  # Reshape beta to a matrix of size p x q 
   # B is a matrix with independent rows and columns

    # Calculate the mu
    mu = X * B * Z'

    # Likelihood: Y ~ MatrixNormal(mu, In, Sigma)
    Y ~ MatrixNormal(mu, I(n), Sigma)
end

# Prior hyperparameters
tau2 = 10.0                              # prior std‐dev on vec(B)
Psi0 = Matrix(I, m, m)                # Inv‐Wishart scale matrix
nu0  = m + 10                            # Inv‐Wishart degrees of freedom

model = bilinear_model(X, Z, Y, tau2, Psi0, nu0)
chain = sample(model, NUTS(), 2_000)