RE your stacking Q, what I meant was the definition of the vec operator here at the top where it says vec(A): Vectorization (mathematics) - Wikipedia
In this regard I have two follow up questions. So in the final application there will be about 5000 observations with a test structure of 4 different condition types (t, S1, S2, S3) yielding a y shape of (5000 x 1) and X with (5000x4). This would result in a covariance matrix in the dimensionality region of somewhat (1e15 x 1e15). Would it be a suitable approach (however that my work) to first model an individual GP for each f_j(t) at a given set of (S1,S2,S3) and finally model a GP of those individual GPs. This would drastically reduce the dimensionality of the covariance matrix and thus spead up the process.
As annoying as it is to answer a question with a question, I think it depends on whether there are any relationships between the 4 condition types. Is there covariance between the S_i, or can you treat each time series as an independent draw from the same GP, like: S_i ~ GP(0, K(t, t’)).
Depends on your computer, but with that dataset size, you may want to look into approximations too, or other time series models. MAP estimation may work fine, but NUTS will start to be a bit slow…