## 16.3 Data Type and Indexing Efficiency

The underlying matrix and linear algebra operations are implemented in terms of data types from the Eigen C++ library. By having vectors and matrices as basic types, no conversion is necessary when invoking matrix operations or calling linear algebra functions.

Arrays, on the other hand, are implemented as instances of the C++

`std::vector`

class (not to be confused with Eigen’s
`Eigen::Vector`

class or Stan vectors). By implementing arrays
this way, indexing is efficient because values can be returned by
reference rather than copied by value.

### Matrices vs. Two-Dimensional Arrays

In Stan models, there are a few minor efficiency considerations in deciding between a two-dimensional array and a matrix, which may seem interchangeable at first glance.

First, matrices use a bit less memory than two-dimensional arrays. This is because they don’t store a sequence of arrays, but just the data and the two dimensions.

Second, matrices store their data in column-major order. Furthermore, all of the data in a matrix is guaranteed to be contiguous in memory. This is an important consideration for optimized code because bringing in data from memory to cache is much more expensive than performing arithmetic operations with contemporary CPUs. Arrays, on the other hand, only guarantee that the values of primitive types are contiguous in memory; otherwise, they hold copies of their values (which are returned by reference wherever possible).

Third, both data structures are best traversed in the order in which they are stored. This also helps with memory locality. This is column-major for matrices, so the following order is appropriate.

```
matrix[M, N] a;
//...
for (n in 1:N)
for (m in 1:M)
// ... do something with a[m, n] ...
```

Arrays, on the other hand, should be traversed in row-major order (i.e., last index fastest), as in the following example.

```
real a[M, N];
// ...
for (m in 1:M)
for (n in 1:N)
// ... do something with a[m, n] ...
```

The first use of `a[m,~n]`

should bring `a[m]`

into memory.
Overall, traversing matrices is more efficient than traversing arrays.

This is true even for arrays of matrices. For example, the ideal order in which to traverse a two-dimensional array of matrices is

```
matrix[M, N] b[I, J];
// ...
for (i in 1:I)
for (j in 1:J)
for (n in 1:N)
for (m in 1:M)
... do something with b[i, j, m, n] ...
```

If `a`

is a matrix, the notation `a[m]`

picks out row
`m`

of that matrix. This is a rather inefficient operation for
matrices. If indexing of vectors is needed, it is much better to
declare an array of vectors. That is, this

```
row_vector[N] b[M];
// ...
for (m in 1:M)
... do something with row vector b[m] ...
```

is much more efficient than the pure matrix version

```
matrix b[M, N];
// ...
for (m in 1:M)
// ... do something with row vector b[m] ...
```

Similarly, indexing an array of column vectors is more efficient than
using the `col`

function to pick out a column of a matrix.

In contrast, whatever can be done as pure matrix algebra will be the fastest. So if I want to create a row of predictor-coefficient dot-products, it’s more efficient to do this

```
matrix[N, k] x; // predictors (aka covariates)
// ...
vector[K] beta; // coeffs
// ...
vector[N] y_hat; // linear prediction
// ...
y_hat = x * beta;
```

than it is to do this

```
row_vector[K] x[N]; // predictors (aka covariates)
// ...
vector[K] beta; // coeffs
...
vector[N] y_hat; // linear prediction
...
for (n in 1:N)
y_hat[n] = x[n] * beta;
```

### (Row) Vectors vs. One-Dimensional Arrays

For use purely as a container, there is really nothing to decide among
vectors, row vectors and one-dimensional arrays. The
`Eigen::Vector`

template specialization and the
`std::vector`

template class are implemented similarly as
containers of `double`

values (the type `real`

in Stan).
Only arrays in Stan are allowed to store integer values.