Return 2D array or Matrix in Torsten

Developers
1

I’ve been going back and forth on this one. Any thoughts or advice are welcomed.

Right now, Torsten returns a matrix. This matrix stores quantities of interest at different times, in the same way the returned object of integrate_ode_* does. Except the latter function returns a 2D array.

In most cases, Torsten computes the elements of the matrix one at a time. When it solves linear ODEs, it calls matrix_exp() and computes the matrix one row at a time. Still, no real gain using a matrix (versus an array of vectors; I save a conversion step I would need to do for a 2D array from eigen::vector to std::vector for linear ODEs. That said with the matrix, I need to convert std::vectorstd::vector to eigen::vector for nonlinear ODEs).

My initial thought was that a matrix would be easier to work with in later parts of the model. Usually, I’m interested in one column and I want to compare it to data in my likelihood function. I understand sampling statements are more efficient with vectors. So it pays off to have cHat (the column of interest) be a vector and it easier to extract it from a matrix.

transformed parameters {
  ...
  x = PKModelTwoCpt(theta, time, amt, rate, ii, evid, cmt, addl, ss);
  cHat = col(x, 2) ./ V1;
  ...
  }
}

However, I only want certain elements of cHat: namely predictions at times where a measurement occurs and I have data. So the object that appears in the likelihood function is cHatObs, not cHat:

transformed parameters {
  ...
  cHat = col(x, 2) ./ V1;
  for(i in 1:nObs){
    cHatObs[i] = cHat[iObs[i]]; ## predictions for observed data records
  }
}

model {
  logCObs ~ normal(log(cHatObs), sigma);
}

While cHatObs should be a vector, the type of cHat doesn’t really matter, since I’ll be going through it one element at a time. It might be more difficult to construct from a 2D array or an array of (row-major) vectors.

All that said, if Torsten returns a 2D array, that would be more consistent with what the ODE integrators do. Was there any reason when designing the ODE solvers for returning a 2D array instead of a matrix?

thank you!

2

charlesm93 Developer
November 30
I’ve been going back and forth on this one. Any thoughts or advice are welcomed.

Right now, Torsten returns a matrix. This matrix stores quantities of interest at different times, in the same way the returned object of integrate_ode_* does. Except the latter function returns a 2D array.

In most cases, Torsten computes the elements of the matrix one at a time. When it solves linear ODEs, it calls matrix_exp() and computes the matrix one row at a time. Still, no real gain using a matrix (versus an array of vectors; I save a conversion step I would need to do for a 2D array from eigen::vector to std::vector for linear ODEs. That said with the matrix, I need to convert std::vector to eigen::vector for nonlinear ODEs).

My initial thought was that a matrix would be easier to work with in later parts of the model. Usually, I’m interested in one column and I want to compare it to data in my likelihood function. I understand sampling statements are more efficient with vectors.

It depends. There’s almost no practical difference
between a 1D array (std::vector), row vector and
column vector (Eigen::Matrix types). In more dimensions,
efficiency depends on what you want to do with the object.
Vectors are better for getting partial results as they don’t
involve copying. So if you have

real a[M, N];

then using a[m] doesn’t require any copying and is nicely
localized, whereas doing the same to pull a row out of a matrix
is more costly. Similarly, if you want to pull out column
vectors, it’s best to use

vector[M] a[N];

and if you want to pull out row vectors, use

row_vector[N] a[M];

(note the implicit transpose).

So it pays off to have cHat (the column of interest) be a vector and it easier to extract it from a matrix.

transformed parameters {
…
x = PKModelTwoCpt(theta, time, amt, rate, ii, evid, cmt, addl, ss);
cHat = col(x, 2) ./ V1;
…
}
}

However, I only want certain elements of cHat: namely predictions at times where a measurement occurs and I have data. So the object that appears in the likelihood function is cHatObs, not cHat:

transformed parameters {
…
cHat = col(x, 2) ./ V1;
for(i in 1:nObs){
cHatObs[i] = cHat[iObs[i]]; ## predictions for observed data records
}
}

you can vectorize this as

cHatObs = cHat[iObs];

model {
logCObs ~ normal(log(cHatObs), sigma);
}

While cHatObs should be a vector, the type of cHat doesn’t really matter, since I’ll be going through it one element at a time. It might be more difficult to construct from a 2D array or an array of (row-major) vectors.

All that said, if Torsten returns a 2D array, that would be more consistent with what the ODE integrators do. Was there any reason when designing the ODE solvers for returning a 2D array instead of a matrix?

I think it was just to keep the dependencies minimal, which
obviously doesn’t matter in the Stan language itself.
Eigen::Matrix<T,-1,-1> is more efficient than std::vector<std::vector>
in terms of storage/allocation and better with memory locality for
column-major layouts.

  • Bob
1 Like
3

Hi，I see you reply.Now I have some trouble for the matrix. For a matrix(p*n),I want to extract each column of the matrix as a column vector named as mi(i=1…n）, then calculate mi times transpose(mi), then we will get n matrices.After that ,I want to calculate the sum of n matrices.
Thank you for you help .

4

What exactly is your question?

Are you unsure how to code this in Stan, or do you need some guidance on how to do this efficiently? I can help you best if you share some code.

5

Thank you very much . Now I have solved it. But as a beginner to STAN, Could you please give me some useful advice and information about it? I shall be very grateful to you.

6

Thank you very much. Now I have solved it. But as a beginner to STAN, Could you please give me some useful advice and information about it? I shall be very grateful to you.

Return 2D array or Matrix in Torsten really difficult. You made it easy.