Differentiate expressions involving symmetric matrix $\mathbf{D}=\mathrm{diag}(\tau)\Omega\mathrm{diag}(\tau)$ with respect to element of $\tau$

Question

I have a square $q$ by $q$ symmetric matrix $\mathbf{D}=\operatorname{diag}(\tau)\Omega\operatorname{diag}(\tau)$ where $\Omega$ is a square $q$ by $q$ matrix, and $\tau$ is a vector of length $q$.

Basically $\mathbf{D}$ is a covariance matrix that I am decomposing into a correlation matrix $\Omega$ and a scale vector $\tau$.

I need to differentiate various functions that contain $\mathbf{D}$ with respect to $\tau_g$, for $g$ in $1, \ldots, q$. Specifically I need to find:

$$\frac{d(\log(|\mathbf{D}|)}{d\tau_g} $$

and

$$\frac{d(\mathbf{b}^T\mathbf{D}^{-1}\mathbf{b})}{d\tau_g} $$

I've been following the matrix cook book (https://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf) but it mentions under the section for derivatives of matrices, that many of the derivaties don't apply in general for matrices with structure (e.g symmetric matrices such as $\mathbf{D}$).

Based on that, how would I go about finding the above two derivatives with respect to $\tau_g$ of an expression involving $\mathbf{D}$? Any pointers / tips / useful identities would be appreciated.

Answer 1

For the first formula, we seem to have: \begin{gather*} |\mathbf{D}| = |\Omega|\tau_1^2\tau_2^2\cdots\tau_q^2, \\ \therefore\ \log|\mathbf{D}| = \log|\Omega| + 2\log\tau_1 + 2\log\tau_2 + \cdots + 2\log\tau_q, \\ \therefore\ \frac{\partial\log|\mathbf{D}|}{\partial\tau_g} = \frac2{\tau_g} \quad (g = 1, \ldots, q). \end{gather*} For the second formula, let $\mathbf{b}^{\mathrm{T}} = (b_1, \ldots, b_q),$ and write $\Omega^{-1} = (c_{hg})_{1\leqslant h \leqslant q, 1 \leqslant g \leqslant q}.$ That is, for $g = 1, \ldots, q,$ let the $g^\text{th}$ column of $\Omega^{-1}$ be $(c_{hg})_{1\leqslant h \leqslant q}^{\mathrm{T}}.$ After a page of messy calculation, I arrive at the formula: $$ \frac{\partial(\mathbf{b}^{\mathrm{T}}\mathbf{D}^{-1}\mathbf{b})}{\partial\tau_g} = -\frac{2b_g}{\tau_g^2}\sum_{h=1}^q\frac{b_hc_{hg}}{\tau_h}. $$

In its dependence on $\tau_g,$ the expression $\mathbf{b}^{\mathrm{T}}\mathbf{D}^{-1}\mathbf{b}$ determines a function $f \colon \mathbb{R}_{>0} \to \mathbb{R},$ where \begin{align*} f(\tau_g) & = \mathbf{b}^{\mathrm{T}}\mathbf{D}^{-1}\mathbf{b} \\ & = \mathbf{b}^{\mathrm{T}}\operatorname{diag}(\tau)^{-1} \Omega^{-1}\operatorname{diag}(\tau)^{-1}\mathbf{b} \\ & = \left(\frac{\mathbf{b}}{\tau}\right)^{\mathrm{T}} \!\Omega^{-1}\left(\frac{\mathbf{b}}{\tau}\right). \end{align*} Here, for the sake of brevity, I have written: $$ \frac{\mathbf{b}}{\tau} = \left( \frac{b_1}{\tau_1}, \ldots, \frac{b_q}{\tau_q}\right)^{\operatorname{T}}. $$

Probably the most sensible way to prove the result is by applying the Chain Rule to the decomposition $f(t) = \beta(\alpha(t)),$ where $\alpha \colon \mathbb{R}_{>0} \to \mathbb{R},$ $t \mapsto b_g/t,$ and $\beta \colon \mathbb{R} \to \mathbb{R},$ $x \mapsto (\mathbf{a} + x\mathbf{d})^{\mathrm{T}}\Omega^{-1}(\mathbf{a} + x\mathbf{d}),$ where $\mathbf{a} = (b_1/\tau_1, \ldots, 0, \ldots, b_q/\tau_q)^{\mathrm{T}}$ and $\mathbf{d} = (0, \ldots, 1, \ldots, 0)^{\mathrm{T}}.$ Such a proof only requires differentiating a quadratic function of $x$, and it has enough structure to inspire a feeling of confidence in the result. Personally, however, I still prefer a more advanced proof, using only the familiar formula $\frac{d}{dt}\frac1t = -\frac1{t^2},$ in conjunction with general results about the Fréchet derivatives of linear and bilinear maps. This gives the final formula more directly, but it requires careful handling, in order to avoid creating thickets of LISP-like nested parentheses (as in my original "messy" handwritten proof). Although it probably can't be recommended objectively, I can't resist giving it here. The function $f$ is expressed in quite a simple way as a composite of four functions, \begin{gather*} f \colon \mathbb{R}_{>0} \xrightarrow{\alpha} \mathbb{R} \xrightarrow{\gamma} \mathbb{R}^q \xrightarrow{\delta} \mathbb{R}^q \times \mathbb{R}^q \xrightarrow{\epsilon} \mathbb{R}, \\ \alpha(t) = \frac{b_g}t \quad (t > 0), \\ \gamma(u) = \mathbf{a} + u\mathbf{d} = \left(\frac{b_1}{\tau_1}, \ldots, u, \ldots, \frac{b_q}{\tau_q}\right)^{\mathrm{T}} \quad (u \in \mathbb{R}), \\ \delta(x) = (x, x) \quad (x \in \mathbb{R}^q), \\ \epsilon(x, y) = x^{\mathrm{T}}\Omega^{-1}y, \quad (x, y \in \mathbb{R}^q). \end{gather*} Clearly, $$ \gamma(\alpha(\tau_g)) = \frac{\mathbf{b}}{\tau}. $$ By the usual rules of differentiation for functions $\mathbb{R}_{>0} \to \mathbb{R},$ $$ \alpha'(t)(k) = -\frac{b_gk}{t^2} \quad (t > 0, \ k \in \mathbb{R}). $$ Because $\gamma$ is the sum of constant and linear mappings, $$ \gamma'(u)(s) = s\mathbf{d} \quad (u, s \in \mathbb{R}). $$ By the Chain Rule, $$(\gamma \circ \alpha)'(t) = \gamma'(\alpha(t)) \circ \alpha'(t) \quad (t > 0). $$ That is, \begin{align*} (\gamma \circ \alpha)'(t)(k) & = \gamma'(\alpha(t))(\alpha'(t)(k))\\ & = \alpha'(t)(k)\mathbf{d} \\ & = -\frac{b_gk}{t^2}\mathbf{d} \quad (t > 0, \ k \in \mathbb{R}). \end{align*} Because $\delta$ is linear, $$ \delta'(x)(v) = \delta(v) \quad (x, v \in \mathbb{R}^q). $$ Because $\epsilon$ is bilinear and symmetric, \begin{align*} \epsilon'(x, y)(v, w) & = \epsilon(x, w) + \epsilon(v, y) \\ & = \epsilon(x, w) + \epsilon (y, v) \quad (x, y, u, v \in \mathbb{R}^q). \end{align*} By the Chain Rule, $$ (\epsilon \circ \delta)'(x) = \epsilon'(\delta(x)) \circ \delta'(x) \quad (x \in \mathbb{R}^q). $$ That is, \begin{align*} (\epsilon \circ \delta)'(x)(v) & = \epsilon'(\delta(x))(\delta'(x)(v)) \\ & = \epsilon'(\delta(x))(\delta(v)) \\ & = \epsilon'(x, x)(v, v) \\ & = 2\epsilon(x, v) \quad (x, v \in \mathbb{R}^q). \end{align*} By the Chain Rule again, $$ f'(t) = (\epsilon \circ \delta)'(\gamma(\alpha(t))) \circ (\gamma \circ \alpha)'(t) \quad (t > 0). $$ That is, \begin{align*} f'(t)(k) & = (\epsilon \circ \delta)'(\gamma(\alpha(t))) ((\gamma \circ \alpha)'(t)(k)) \quad (t > 0, \ k \in \mathbb{R}). \end{align*} Therefore, \begin{align*} f'(\tau_g)(k) & = (\epsilon \circ \delta)' \left(\frac{\mathbf{b}}{\tau}\right) \left(-\frac{b_gk}{t^2}\mathbf{d}\right) \\ & = 2\epsilon\left(\frac{\mathbf{b}}{\tau}, -\frac{b_gk}{\tau_g^2}\mathbf{d}\right) \\ & = -\frac{2b_gk}{\tau_g^2}\epsilon\left( \frac{\mathbf{b}}{\tau}, \mathbf{d}\right) \\ & = -\frac{2b_gk}{\tau_g^2} \left( \frac{b_1}{\tau_1}, \ldots, \frac{b_q}{\tau_q}\right) \Omega^{-1}\mathbf{d} \\ & = -\frac{2b_gk}{\tau_g^2} \left(\frac{b_1}{\tau_1}, \ldots, \frac{b_q}{\tau_q}\right) (c_{1g}, \ldots, c_{qg})^{\mathrm{T}} \\ & = -\frac{2b_gk}{\tau_g^2}\sum_{h=1}^q\frac{b_hc_{hg}}{\tau_h} \quad (k \in \mathbb{R}). \end{align*}

To conclude, here are the remaining details of the more "sensible" (although in my opinion less intuitive) proof that was outlined earlier. By the Chain Rule (it is no longer expressed in terms of Fréchet derivatives, of course): $$ f'(t) = \beta'(\alpha(t)) \cdot \alpha'(t) \quad (t > 0). $$ We have already differentiated $\alpha,$ thus: $$ \alpha'(t) = -\frac{b_g}{t^2} \quad (t > 0). $$ We can also easily differentiate $\beta,$ because it is a quadratic function. Using the symmetry of $\Omega,$ we have: \begin{align*} \beta(x) & = (\mathbf{a} + x\mathbf{d})^{\mathrm{T}}\Omega^{-1} (\mathbf{a} + x\mathbf{d}) \\ & = (\mathbf{d}^{\mathrm{T}}\Omega^{-1}\mathbf{d})x^2 + (\mathbf{d}^{\mathrm{T}}\Omega^{-1}\mathbf{a})x + (\mathbf{a}^{\mathrm{T}}\Omega^{-1}\mathbf{d})x + \mathbf{a}^{\mathrm{T}}\Omega^{-1}\mathbf{a} \\ & = (\mathbf{d}^{\mathrm{T}}\Omega^{-1}\mathbf{d})x^2 + 2(\mathbf{a}^{\mathrm{T}}\Omega^{-1}\mathbf{d})x + \mathbf{a}^{\mathrm{T}}\Omega^{-1}\mathbf{a} \quad (x \in \mathbb{R}), \end{align*} therefore \begin{align*} \beta'(x) & = 2(\mathbf{d}^{\mathrm{T}}\Omega^{-1}\mathbf{d})x + 2(\mathbf{a}^{\mathrm{T}}\Omega^{-1}\mathbf{d}) \\ & = 2(\mathbf{a} + x\mathbf{d})^{\mathrm{T}}\Omega^{-1}\mathbf{d} \quad (x \in \mathbb{R}). \end{align*} But $$ \mathbf{a} + \alpha(\tau_g)\mathbf{d} = \frac{\mathbf{b}}{\tau}, $$ therefore $$ f'(\tau_g) = -\frac{2b_g}{\tau_g^2} \left(\frac{\mathbf{b}}{\tau}\right)^{\mathrm{T}} \!\!\Omega^{-1}\mathbf{d}, $$ the same expression that emerged more transparently from the other proof.

Answer 2

Let's use a colon as a convenient product notation for the trace, i.e. $$A:B = {\rm Tr}(A^TB)$$ and note that the terms in such a product can be rearranged in a variety of ways in accordance with the cyclic property of the trace, e.g. $$\eqalign{ A:B &= A^T:B^T &= B:A \\ A:BC &= AC^T:B &= B^TA:C \\ }$$ The $D$ matrix has an alternative form in terms of the elementwise/Hadamard product. $$\eqalign{ D &= {\rm Diag}(\tau)\,\Omega\,{\rm Diag}(\tau) = \Omega\odot \tau\tau^T \\ dD &= \Omega\odot d(\tau\tau^T) = \Omega\odot(d\tau\,\tau^T+\tau\,d\tau^T) \\ }$$ The latter expression is the differential of $D$ and it will be substituted in several places below.

Calculate the differential and gradient of the first function function. $$\eqalign{ \phi &= \log\det D \\ d\phi &= D^{-1}:dD \\ &= D^{-1}:(\Omega\odot(d\tau\,\tau^T+\tau\,d\tau^T)) \\ &= \Omega\odot D^{-1}:(d\tau\,\tau^T+\tau\,d\tau^T) \\ &= 2(\Omega\odot D^{-1}):d\tau\,\tau^T \\ &= 2(\Omega\odot D^{-1})\tau:d\tau \\ \frac{\partial \phi}{\partial \tau} &= 2(\Omega\odot D^{-1})\tau \\ }$$ Multiply by the cartesian base vector $e_k$ to extract the derivative with respect to $\tau_k$ $$\eqalign{ \frac{d\phi}{d\tau_k} = e_k^T\left(\frac{\partial \phi}{\partial \tau}\right) &= 2e_k^T(\Omega\odot D^{-1})\tau \\ &= 2\tau^T(\Omega\odot D^{-1})e_k \quad&({\rm symmetric\, matrices\, yay!}) \\ &= 2(\Omega\odot D^{-1}):\tau e_k^T \quad&({\rm a\, better\, way\, to\, write\, it?}) \\ }$$ Calculate the gradient of the second function. $$\eqalign{ \psi &= bb^T:D^{-1} \\ d\psi &= bb^T:dD^{-1} \\ &= -bb^T:D^{-1}\,dD\,D^{-1} \\ &= -D^{-1}bb^TD^{-1}:dD \\ &= -D^{-1}bb^TD^{-1}:(\Omega\odot(d\tau\,\tau^T+\tau\,d\tau^T)) \\ &= -\Omega\odot(D^{-1}bb^TD^{-1}):(d\tau\,\tau^T+\tau\,d\tau^T) \\ &= -2\Omega\odot(D^{-1}bb^TD^{-1}):d\tau\,\tau^T \\ &= -2\Big(\Omega\odot(D^{-1}bb^TD^{-1})\Big)\tau:d\tau \\ \frac{\partial \psi}{\partial \tau} &= -2\Big(\Omega\odot(D^{-1}bb^TD^{-1})\Big)\tau \\ }$$ Once again, the derivative with respect to $\tau_k$ is obtained by multiplying the vector-valued gradient with the basis vector $e_k$

NB: The Hadamard and colon products commute with themselves and each other. $$\eqalign{ A:B &= B:A \\ A\odot B &= B\odot A \\ C:A\odot B &= C\odot A:B \\ }$$ Also a symmetric matrix like $D$ permits the following useful reduction $$\eqalign{ D:(AB^T+BA^T) &= D:AB^T + D:BA^T \\ &= D:AB^T + D^T:AB^T \\ &= (D+D^T):AB^T \\ &= 2D:AB^T \\ }$$ which was employed in several steps above.