# Bilinear Model in Turing
**URL:** https://discourse.julialang.org/t/bilinear-model-in-turing/131614
**Category:** Statistics
**Tags:** question
**Created:** [August 14, 2025, 8:37pm UTC](https://discourse.julialang.org/t/bilinear-model-in-turing/131614 "2025-08-14T20:37:55Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@hdubey](https://discourse.julialang.org/u/hdubey)
#### Post date: [August 14, 2025, 8:37pm UTC](https://discourse.julialang.org/t/bilinear-model-in-turing/131614/1 "2025-08-14T20:37:55Z")
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:
```julia-auto
@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)
```
---
### Author:  [@eteppo](https://discourse.julialang.org/u/eteppo)
#### Post date: [August 16, 2025, 3:15pm UTC](https://discourse.julialang.org/t/bilinear-model-in-turing/131614/2 "2025-08-16T15:15:49Z")
Welcome! You may get help easier in the [Julia Slack](https://julialang.org/slack/) channel for [Turing](https://julialang.slack.com/archives/CCYDC34A0) or in [GitHub](https://github.com/TuringLang) issues. It’s best to give a fully runnable example. I strangely hit some other type inference issue.
---
### Author:  [@hdubey](https://discourse.julialang.org/u/hdubey)
#### Post date: [August 18, 2025, 3:23pm UTC](https://discourse.julialang.org/t/bilinear-model-in-turing/131614/3 "2025-08-18T15:23:02Z")
Thank you. I have posted this on the Slack channel for Julia as well. I am also pasting a fully runnable example below.
```julia-auto
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)
```