Hi there,

I’ve built a multilevel, multi-variable linear regression model inspired by the example in section 1.13 of the Stan manual. We have N samples organised into J groups. K is the number of predictors. The model looks as follows:

```
# X_train = (centered) training data
# y_train = observed regression targets
# QR reparametrization
Q_train_, R_train_ = np.linalg.qr(X_train, mode='reduced')
Q_train = Q_train_ * np.sqrt(n_samples - 1)
R_train = R_train_ / np.sqrt(n_samples - 1)
R_train_inv = np.linalg.inv(R_train) # i know...
# (...)
with pm.Model() as linear_model:
# std for likelihood
s = pm.HalfCauchy('sigma', 2.5)
# covariance matrix
packed_L_Omega = pm.LKJCholeskyCov('packed_L_Omega', n=n_predictors, eta=3, sd_dist=pm.HalfCauchy.dist(2.5))
L_Omega = pm.expand_packed_triangular(n_predictors, packed_L_Omega)
# group mean
mu = pm.Normal('mu', 0., 5., shape=(n_predictors, n_groups))
# group slopes
beta_mat = [pm.MvNormal(f'_beta_{i_grp:04d}', mu[:, i_grp], chol=L_Omega, shape=n_predictors) for i_grp in range(n_groups)]
beta = pm.Deterministic('beta', tt.stack(beta_mat, axis=1))
# group intercepts
alpha = pm.Normal('alpha', 0., 5., shape=(n_samples, n_groups))
# expected value
y_hat = alpha[:, group_idx] + tt.dot(Q_train, beta[:, group_idx])
# data likelihood
y = pm.Normal('y', mu=y_hat, sd=s, observed=y_train)
# sample
trace = pm.sample(draws=1000, chains=2, cores=1, tune=1000)
# PPC
posterior_predictive = pm.sample_posterior_predictive(trace)
```

Two questions:

**1/** the sampler starts up, but it is *very* slow. N \approx 3500, `pymc3`

estimates one chain will take approx 2 days to run

Is there an obvious way to speed this up? I was hoping to scale this model up to 350K samples and 200 predictors, is there any real hope of doing that? CPU power is not an issue, as I am planning to run this on a supercomputer.

**2/** I get some weird Theano warnings

```
python3.7/site-packages/theano/tensor/basic.py:6611: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
result[diagonal_slice] = x
python3.7/site-packages/theano/tensor/subtensor.py:2339: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
out[0][inputs[2:]] = inputs[1]
python3.7/site-packages/theano/tensor/subtensor.py:2197: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
rval = inputs[0].__getitem__(inputs[1:])
python3.7/site-packages/theano/tensor/subtensor.py:2339: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
out[0][inputs[2:]] = inputs[1]
python3.7/site-packages/theano/tensor/basic.py:6611: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
result[diagonal_slice] = x
python3.7/site-packages/theano/tensor/subtensor.py:2197: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
rval = inputs[0].__getitem__(inputs[1:])
```

Any idea what part of the code is generating this?

Thank you!