I have a model of the form x_{t+1}=A x_{t} (a dynamical system) and want to regress the coefficients of A \in \mathbb{R}^{n\times n} onto x_{t+1} \in \mathbb{R}^n given m data pairs (x_{t+1},x_{t}). For this task, I have yet to find an approach to implement Bayesian regression. Help getting started with this would be greatly appreciated.
On top of that task, I would like to add a regularizing term which penalizes the Frobenius norm of A (equivalently an L2 penalization on both the rows and columns of A).
Here’s a model that should implement what you’re looking for. I would normally suggest using one of the timeseries distributions since you are describing a vector autoregressive order-1 system, but there are some rough points in dealing with shapes there.
import numpy as np
import pymc3 as pm
import theano.tensor as tt
def frobenius_norm(X):
return tt.sum(tt.nlinalg.trace(A@A.T))**0.5
# Create simulated data via forward evolution of system
K = 3
T = 20
error_sd_true = 0.5
A_true = np.random.randn(K,K)
x = np.zeros([T, K])
x[0] = np.random.randn(K)
for t in range(1, T):
jump = np.random.randn(K)
x[t] = A_true@x[t-1] + jump*error_sd_true
penalty_weight = 0.1
with pm.Model():
A_scale = pm.HalfCauchy('A_scale', beta=1)
A = pm.Normal('A', sd=A_scale, shape=[K,K])
error_sd = pm.HalfCauchy('error_sd', beta=1)
# We regress the shifted observations on the earlier observations
# with the regression weights provided by A. Note that the right-multiplication
# with A is only because of the shape we set up earlier for x.
pm.Normal('x_end', mu=x[0:-1]@A, observed=x[1:], sd=error_sd)
pm.Potential('A_penalty', frobenius_norm(A) * penalty_weight)
trace = pm.sample()
Now, you may want to play around with the choice of prior for A. Right now, there are two sources of regularization - one from the pm.Normal prior placed on its entries (essentially an element-wise L2 loss) and the Frobenius norm calculated using the trace of A. I also tried out a version omitting the pm.Normal prior and instead replacing it with pm.Flat which applies no prior. It works, but it just takes longer to sample because the problem is more poorly conditioned.
Would you be able to share your PyMC and Theano versions with pm.version and theano.version? Mine was running on PyMC 3.11.4 and Theano 1.1.2. Also, I edited the example because for some random seeds, the data were not sufficiently informative to make the problem tractable.
I can try downgrading PyMC and see if that solves the problem.
Another question I had if I may ask:
I read through the paper on Probabilistic Matrix Factorization (Probabilistic Matrix Factorization) and noticed that simply by taking normal priors over the columns of the factor matrices (U and V), the Frobenius norms were directly penalized when considering the sum-of-squared errors objective function (eq. 4). Does this imply that I shouldn’t be explicitly penalizing the Frobenius norm of the matrix A if there is already a normal prior? Sorry for the unrelated question, but thought I’d ask and continue to pick your brain if you have the time.
I haven’t upgraded to PyMC v4, so that would likely be the source of the error.
With regard to the penalization, I read the part you mentioned in that paper where the authors refer to an IID normal prior being equivalent to a Frobenius norm on A. I didn’t know that ahead of time but it’s not surprising to me that it’s the case - it looks like you can omit the Frobenius penalty from the code I wrote.
