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)
2

Welcome! You may get help easier in the Julia Slack channel for Turing or in GitHub issues. It’s best to give a fully runnable example. I strangely hit some other type inference issue.

3

Thank you. I have posted this on the Slack channel for Julia as well. I am also pasting a fully runnable example below.

using Turing, Distributions, LinearAlgebra, Random, Plots, StatsPlots, MatrixLM, DataFrames, StatsModels, PDMats

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

    # Prior on Sigma: Inverse Wishart
    #covY = cov(Y)
    #idx  = diagind(covY)               # indices of the diagonal
    #S    = Psi0 * Diagonal(covY[idx])  # Diagonal matrix with Psi0 * diag(covY)
    #S_dense = Matrix(S) 
    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 ~ MatrixNormal(zeros(p, q), tau2*I(p), tau2*I(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

# Dimensions of matrices 
n = 100
m = 250

# Number of groupings designed in the Z matrix
q = 3

# Generate data with two categorical variables and 3 numerical variables.
dfX = hcat(
    DataFrame(catvar1=string.(rand(0:1, n)), catvar2=rand(["A", "B", "C", "D"], n)), 
    DataFrame(rand(n,3), ["var3", "var4", "var5"])
    );

# Convert dataframe to predicton matrix
X = design_matrix(@mlmformula(catvar1 + catvar2 + var3 + var4 + var5), dfX)

p = size(X)[2];

# Convert dataframe to predicton matrix
my_ctrst = Dict(
             :catvar1 => DummyCoding(base = "0"),
             :catvar2 => DummyCoding(base = "A")
           )
           
X = design_matrix(@mlmformula(catvar1 + catvar2 + var3 + var4 + var5), dfX, my_ctrst);

Z = rand(m,q);
E = randn(n,m).*4;

# (p,q)
B = [
    2.0   3.0   4.0;
    0.01  0.02  0.01;
    -1.0  -0.5  0.02;
    -0.01 0.02  -0.01;
    0.0   0.0   0.0;
    3.0   3.0   3.0;
    0.01  0.0   3.0;
];

Y = X*B*Z' + E;

# Prior hyperparameters
tau2 = 10.0                              # prior std‐dev on vec(B)
Psi0 = Matrix(I, m, m)
#Psi0 = sigma_hat                             # 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)