Computing the diagonal of a matrix product

Naively, one would compute the matrix product and then takes it diagonal, but this is inefficient.

Is there a way, in C++ Stan or in Eigen, to only compute the diagonal elements of a matrix product?

Just break up the calculation analytically – any element of the output matrix is the dot product of two vectors – a row from one of the input matrices and a column from the other,

A_{ii} = \sum_{j = 1}^{J} B_{ij} C_{ji},

or in Eigen,

B.col(j) * C.row(j)

Eigen’s lazy evaluation takes care of this automatically, where (B * C).diagonal() will compute the diagonal coefficients.

See this SO comment from one the developers:

1 Like

Thanks for the insight. Note that using the matrix operations in Stan doesn’t do the lazy evaluation. Indeed:

diagonal(multiply(A, B))

is much slower than:

(A * B).diagonal()

where the latter uses only Eigen’s functions.

That’s right. Stan’s matrix functions evaluate the entire matrix as they’re called, whereas Eigen’s functions use expression templates to identify which parts of the matrix to evaluate