I am creating a model for tensor data, but to my surprise I see there is no native data type for tensors. By this I mean, for example, having a 3-D tensor T of dimension [N,N,N] for which taking a slice T[1, 1:N, 1:N] would return a matrix with all of its properties(especially matrix multiplication and hadamard product).
Currently I create a 3D array, but then have to call a bunch of times the to_matrix and to_vector functions to do the necessary operations. Especially to my surprise these is no element-wise product for 2D arrays. All these changing types I am pretty convinced are not computationally efficient at all and are significantly slowing down the code(is this a reasonable concern?).
I wanted hence to know if there is actually such a data structure and I just did not find it, or if this is not the case if it is planned to be added. Also if anyone would have some suggestions on how to make this more efficient it would be very helpful.
which makes T an array of matrices. Then it’s efficient to access matrix T[n].
The closest we have is elementwise products for matrices. We don’t generally do arithmetic on arrays in Stan, though some of the reductions like sum and softmax and log_sum_exp are overloaded for arrays.
This was a fundamental mistake in the design for which I can take responsibility. We split data types into arrays and linear algebra types. The array types are implemented as C++ vectors, so they don’t even have memory locality for 2D arrays (though it’s efficient to pull out a row). We developed this all before TensorFlow, PyTorch, etc. and before we realized just how important tensor reshaping and batching would be at larger scales. Rewriting now would require more C++ and parser/language devs than we have.
That’s a reasonable concern. Each reshape does an allocation and copy, though we’re working on views so it doesn’t require reallocation.
Great I did not know I could create a list of matrices, that should be definitely good enough! Maybe you could mentioned in the reference manual(I could not find it). I am not sure though you can use a list of matrices as an input to reduce_sum, which is something I could two with 3D arrays to speed up computations; If I remember correctly reduce_sum does not like those composite types. I will try and see if it goes trough.
Unfortunately I remembered correctly, reduce_sum cannot process those composite types, or I am doing something wrong. I created a toy model as such to test
functions {
real loglik_t(matrix At, matrix pit){
return (1.0);
}
real logfull(array[] matrix Ats){
real lpdf_slice = 0.0;
int Tslice = 5;
for (t in 1:Tslice){
/// This is symmetric
lpdf_slice = lpdf_slice + loglik_t(Ats[t]);
}
return lpdf_slice;
}
real partial_sums(array[] matrix slice_Ats, int start, int end){
return (logfull(slice_Ats));
}
}
// The input data is a vector 'y' of length 'N'.
data {
int<lower=1> N; // number of nodes
int<lower=1> T0; // number of time steps
array[T0] matrix[N, N] Ats;
}
parameters {
real mu;
real<lower=0> sigma;
}
model {
int grainsize = 1;
target += reduce_sum(partial_sums, Ats, grainsize,);
}
But when I then try to run I get the following error message
Error in stanc(file = file, model_code = model_code, model_name = model_name, :
0
Syntax error in 'string', line 31, column 52 to column 53, parsing error:
Ill-formed expression. We expect a comma separated list of expressions.