# Computing Row-Wise Divergence of Matrix-Valued Functions by Automatic Differentiation

Suppose I have some matrix-valued function $$B : \mathbf{R}^d \to \mathbf{R}^{d \times d}$$, and I want to compute the row-wise divergence of $$B$$, given in vector form as

\begin{align} \alpha &:= \text{div} \,B(x) \in \mathbf{R}^d, \\ \text{where} \quad\alpha_i &= \sum_{j \in [d]} \frac{\partial B_{ij}}{\partial x_j} (x) \quad \text{for} \, i \in [d]. \end{align}

Imagine that $$B$$ is defined by a composition of several smooth functions, and so there is no concern about whether or not these quantities exist. Moreover, this vector can certainly computed, in the sense that I can write out the matrix element-by-element, compute the partial derivatives $$\frac{\partial B_{ij}}{\partial x_j}$$ one-by-one, and then aggregate them.

However, I'd like to avoid this. In my application, while it is easy to compute the action of $$B$$ on a vector, the matrix itself is dense (roughly speaking, it has a diagonal-plus-low-rank form). As such, actually writing down the whole matrix $$B$$ seems like it will be a strain on e.g. memory.

With these concerns in mind, I am (perhaps optimistically) looking for an algorithmic solution which will allow me to exactly compute $$\alpha$$ without ever instantiating the matrix $$B$$, and only ever computing matrix-vector products. Is this possible, and if so, how would such a solution work?