5.4 Vector and Matrix Data Types
Stan provides three types of container objects: arrays, vectors, and matrices. Vectors and matrices are more limited kinds of data structures than arrays. Vectors are intrinsically one-dimensional collections of reals, whereas matrices are intrinsically two dimensional. Vectors, matrices, and arrays are not assignable to one another, even if their dimensions are identical. A \(3 \times 4\) matrix is a different kind of object in Stan than a \(3 \times 4\) array.
The intention of using matrix types is to call out their usage in the code. There are three situations in Stan where only vectors and matrices may be used,
- matrix arithmetic operations (e.g., matrix multiplication)
- linear algebra functions (e.g., eigenvalues and determinants), and
- multivariate function parameters and outcomes (e.g., multivariate normal distribution arguments).
Vectors and matrices cannot be typed to return integer values. They
are restricted to real values.2
Indexing from 1
Vectors and matrices, as well as arrays, are indexed starting from one in Stan. This follows the convention in statistics and linear algebra as well as their implementations in the statistical software packages R, MATLAB, BUGS, and JAGS. General computer programming languages, on the other hand, such as C++ and Python, index arrays starting from zero.
Vectors
Vectors in Stan are column vectors; see the next subsection for
information on row vectors. Vectors are declared with a size (i.e., a
dimensionality). For example, a 3-dimensional vector is declared with
the keyword vector, as follows.
vector[3] u;Vectors may also be declared with constraints, as in the following declaration of a 3-vector of non-negative values.
vector<lower=0>[3] u;Similarly, they may be declared with a offset and/or multiplier, as in the following example
vector<offset=42,multiplier=3>[3] u;Unit Simplexes
A unit simplex is a vector with non-negative values whose entries sum
to 1. For instance, \([0.2,0.3,0.4,0.1]^{\top}\) is a unit 4-simplex.
Unit simplexes are most often used as parameters in categorical
or multinomial distributions, and they are also the sampled variate in
a Dirichlet distribution. Simplexes are declared with their full
dimensionality. For instance, theta is declared to
be a unit \(5\)-simplex by
simplex[5] theta;Unit simplexes are implemented as vectors and may be assigned to other vectors and vice-versa. Simplex variables, like other constrained variables, are validated to ensure they contain simplex values; for simplexes, this is only done up to a statically specified accuracy threshold \(\epsilon\) to account for errors arising from floating-point imprecision.
In high dimensional problems, simplexes may require smaller step sizes in the inference algorithms in order to remain stable; this can be achieved through higher target acceptance rates for samplers and longer warmup periods, tighter tolerances for optimization with more iterations, and in either case, with less dispersed parameter initialization or custom initialization if there are informative priors for some parameters.
Unit Vectors
A unit vector is a vector with a norm of one. For instance, \([0.5, 0.5, 0.5, 0.5]^{\top}\) is a unit 4-vector. Unit vectors are sometimes
used in directional statistics. Unit vectors are declared with their
full dimensionality. For instance, theta is declared to be a unit
\(5\)-vector by
unit_vector[5] theta;Unit vectors are implemented as vectors and may be assigned to other vectors and vice-versa. Unit vector variables, like other constrained variables, are validated to ensure that they are indeed unit length; for unit vectors, this is only done up to a statically specified accuracy threshold \(\epsilon\) to account for errors arising from floating-point imprecision.
Ordered Vectors
An ordered vector type in Stan represents a vector whose entries are sorted in ascending order. For instance, \((-1.3,2.7,2.71)^{\top}\) is an ordered 3-vector. Ordered vectors are most often employed as cut points in ordered logistic regression models (see section).
The variable c is declared as an ordered 5-vector by
ordered[5] c;After their declaration, ordered vectors, like unit simplexes, may be assigned to other vectors and other vectors may be assigned to them. Constraints will be checked after executing the block in which the variables were declared.
Positive, Ordered Vectors
There is also a positive, ordered vector type which operates similarly to ordered vectors, but all entries are constrained to be positive. For instance, \((2,3.7,4,12.9)\) is a positive, ordered 4-vector.
The variable d is declared as a positive, ordered 5-vector by
positive_ordered[5] d;Like ordered vectors, after their declaration, positive ordered vectors may be assigned to other vectors and other vectors may be assigned to them. Constraints will be checked after executing the block in which the variables were declared.
Row Vectors
Row vectors are declared with the keyword row_vector.
Like (column) vectors, they are declared with a size. For example,
a 1093-dimensional row vector u would be declared as
row_vector[1093] u;Constraints are declared as for vectors, as in the following example of a 10-vector with values between -1 and 1.
row_vector<lower=-1,upper=1>[10] u;Offset and multiplier are also similar as for the following 3-row-vector with offset -42 and multiplier 3.
row_vector<offset=-42,multiplier=3>[3] u;Row vectors may not be assigned to column vectors, nor may column vectors be assigned to row vectors. If assignments are required, they may be accommodated through the transposition operator.
Matrices
Matrices are declared with the keyword matrix along with a
number of rows and number of columns. For example,
matrix[3, 3] A;
matrix[M, N] B;declares A to be a \(3 \times 3\) matrix and B to be a \(M \times N\) matrix. For the second declaration to be well formed, the
variables M and N must be declared as integers in either
the data or transformed data block and before the matrix declaration.
Matrices may also be declared with constraints, as in this (\(3 \times 4\)) matrix of non-positive values.
matrix<upper=0>[3, 4] B;Similarly, matrices can be declared to have a set offset and/or multiplier, as in this matrix with multiplier 5.
matrix<multiplier=5>[3, 4] B;Assigning to Rows of a Matrix
Rows of a matrix can be assigned by indexing the left-hand side of an assignment statement. For example, this is possible.
matrix[M, N] a;
row_vector[N] b;
// ...
a[1] = b;This copies the values from row vector b to a[1], which
is the first row of the matrix a. If the number of columns in
a is not the same as the size of b, a run-time error is
raised; the number of columns of a is N, which is also the
number of columns of b.
Assignment works by copying values in Stan. That means any subsequent
assignment to a[1] does not affect b, nor does an
assignment to b affect a.
Covariance Matrices
Matrix variables may be constrained to represent covariance matrices. A matrix is a covariance matrix if it is symmetric and positive definite. Like correlation matrices, covariance matrices only need a single dimension in their declaration. For instance,
cov_matrix[K] Omega;declares Omega to be a \(K \times K\) covariance matrix, where
\(K\) is the value of the data variable K.
Correlation Matrices
Matrix variables may be constrained to represent correlation matrices. A matrix is a correlation matrix if it is symmetric and positive definite, has entries between \(-1\) and \(1\), and has a unit diagonal. Because correlation matrices are square, only one dimension needs to be declared. For example,
corr_matrix[3] Sigma;declares Sigma to be a \(3 \times 3\) correlation matrix.
Correlation matrices may be assigned to other matrices, including unconstrained matrices, if their dimensions match, and vice-versa.
Cholesky Factors of Covariance Matrices
Matrix variables may be constrained to represent the Cholesky factors of a covariance matrix. This is often more convenient or more efficient than representing covariance matrices directly.
A Cholesky factor \(L\) is an \(M \times N\) lower-triangular matrix (if \(m < n\) then \(L[m, n] =0\)) with a strictly positive diagonal (\(L[k, k] > 0\)) and \(M \geq N\). If \(L\) is a Cholesky factor, then \(\Sigma = L \, L^{\top}\) is a covariance matrix (i.e., it is positive definite). The mapping between positive definite matrixes and their Cholesky factors is bijective—every covariance matrix has a unique Cholesky factorization.
The typical case of a square Cholesky factor may be declared with a single dimension,
cholesky_factor_cov[4] L;Cholesky factors of positive semi-definite matrices
In general, two dimensions may be declared, with the above being equal to
cholesky_factor_cov[4, 4]. The
type cholesky_factor_cov[M, N] may be used for the general
\(M \times N\) case to produce positive semi-definite matrices of rank \(M\).
Cholesky Factors of Correlation Matrices
Matrix variables may be constrained to represent the Cholesky factors of a correlation matrix.
A Cholesky factor for a correlation matrix \(L\) is a \(K \times K\) lower-triangular matrix with positive diagonal entries and rows that are of length 1 (i.e., \(\sum_{n=1}^K L_{m,n}^2 = 1\)). If \(L\) is a Cholesky factor for a correlation matrix, then \(L\,L^{\top}\) is a correlation matrix (i.e., symmetric positive definite with a unit diagonal).
To declare the variable L to be a K by K Cholesky factor of a
correlation matrix, the following code may be used.
cholesky_factor_corr[K] L;Assigning Constrained Variables
Constrained variables of all types may be assigned to other variables
of the same unconstrained type and vice-versa. Matching is interpreted strictly
as having the same basic type and number of array dimensions.
Constraints are not considered, but basic data types are. For instance, a
variable declared to be real<lower=0,upper=1> could be assigned
to a variable declared as real and vice-versa. Similarly, a
variable declared as matrix[3, 3] may be assigned to a variable
declared as cov_matrix[3] or
cholesky_factor_cov[3], and vice-versa.
Checks are carried out at the end of each relevant block of statements to ensure constraints are enforced. This includes run-time size checks. The Stan compiler isn’t able to catch the fact that an attempt may be made to assign a matrix of one dimensionality to a matrix of mismatching dimensionality.
Expressions as Size Declarations
Variables may be declared with sizes given by expressions. Such expressions are constrained to only contain data or transformed data variables. This ensures that all sizes are determined once the data is read in and transformed data variables defined by their statements. For example, the following is legal.
data {
int<lower=0> N_observed; int<lower=0> N_missing;
// ...
transformed parameters {
vector[N_observed + N_missing] y;
// ...Accessing Vector and Matrix Elements
If v is a column vector or row vector, then v[2] is the
second element in the vector. If m is a matrix, then
m[2, 3] is the value in the second row and third column.
Providing a matrix with a single index returns the specified row. For
instance, if m is a matrix, then m[2] is the second row.
This allows Stan blocks such as
matrix[M, N] m;
row_vector[N] v;
real x;
// ...
v = m[2];
x = v[3]; // x == m[2][3] == m[2, 3]The type of m[2] is row_vector because it is the second
row of m. Thus it is possible to write m[2][3] instead
of m[2, 3] to access the third element in the second row. When
given a choice, the form m[2, 3] is preferred.
Array index style
The form m[2, 3] is more efficient because it does not require the creation and use of an intermediate expression template for m[2]. In later versions, explicit calls to m[2][3] may be optimized to be as efficient as m[2, 3] by the Stan compiler.
Size Declaration Restrictions
An integer expression is used to pick out the sizes of vectors,
matrices, and arrays. For instance, we can declare a vector of size
M + N using
vector[M + N] y;Any integer-denoting expression may be used for the size declaration, providing all variables involved are either data, transformed data, or local variables. That is, expressions used for size declarations may not include parameters or transformed parameters or generated quantities.
This may change if Stan is called upon to do complicated integer matrix operations or boolean matrix operations. Integers are not appropriate inputs for linear algebra functions.↩