Internal helper to compute \(\mathrm{tr}(A B)\) without forming A %*% B:
$$\mathrm{tr}(A B) = \sum (A \circ B^\top).$$
This identity is used repeatedly in the fully exponential (FE) trace
corrections (Karl, Yang, and Lohr, 2014, Appendix B), where a naive
diag(A %*% B) would allocate the full product.
glmmfe_trAB(A, B)A single numeric scalar equal to tr(A %*% B).
Numeric matrices with identical dimensions.