Calculating the gradient of Log-Euclidean distance between SPD matrices on Riemannian manifold

Question

In the paper Log-Euclidean metrics for fast and simple calculus on diffusion tensors, the geodesic distance between SPD matrices $A,B$ is defined as $$d(A,B)=||\log A- \log B||_F,$$ where $F$ is the Frobenius norm. I did not find more information on the other aspects of the Riemannian structure.

I am interested in the gradient (on Riemannian manifold) of the distance function $$ f(X) = d(X,A_0). $$ I am not familiar with geometry other than some basics, $\langle grad f,X\rangle = Xf $ for all smooth vector field $X$.

How can I derive $grad f$?

Answer 1

$ \def\op#1{\operatorname{#1}} \def\bbR#1{{\mathbb R}^{#1}} \def\e{\varepsilon} \def\o{{\tt1}} \def\f{\frac{\tt1}{f}} \def\p{\partial} \def\A{{\cal A}}\def\L{{\cal L}} \def\LR#1{\left(#1\right)} \def\BR#1{\Big(#1\Big)} \def\trace#1{\op{Tr}\LR{#1}} \def\Diag#1{\op{Diag}\LR{#1}} \def\diag#1{\op{diag}\LR{#1}} \def\qiq{\quad\implies\quad} \def\grad#1#2{\frac{\p #1}{\p #2}} \def\m#1{\left[\begin{array}{r}#1\end{array}\right]} \def\mb#1{\left[\begin{array}{c|c}#1\end{array}\right]} \def\c#1{\color{red}{#1}} \def\CLR#1{\c{\LR{#1}}} \def\fracLR#1#2{\LR{\frac{#1}{#2}}} \def\gradLR#1#2{\LR{\grad{#1}{#2}}} \def\z#1{\op{\zeta}\!\LR{#1}} \def\smA{{\small A}} \def\smB{{\small B}} $Introduce the matrix variables $$\eqalign{ A &= X,\quad &L = \log(A), \quad &L_0 = \log(A_0) \\ }$$ Since $A$ is positive definite it can be diagonalized $$A=QBQ^T,\quad B=\Diag{b},\quad Q^TQ=I$$

First calculate the differential of the distance function $f$ $$\eqalign{ f^2 &= \|L-L_0\|^2_F \\ &= \LR{L-L_0}:\LR{L-L_0} \\ 2f\;df &= 2\LR{L-L_0}:dL \\ df &= \f\LR{L-L_0}:dL \\ df &= G:\c{dL} \qiq G \doteq \fracLR{L-L_0}{f} \\ }$$ Now invoke the Daleckii-Krein theorem $\:\big(\odot$ denotes the Hadamard product$\big)$ $$\eqalign{ \c{dL} &= Q\BR{R\odot\LR{Q^TdA\,Q}}Q^T \\ }$$ Substituting this into the previous differential leads to the desired gradient $$\eqalign{ df &= G:\c{Q\BR{R\odot\LR{Q^TdA\,Q}}Q^T} \\ &= Q\BR{R\odot\LR{Q^TG\,Q}}Q^T:dA \\ \grad{f}{A} &= Q\BR{R\odot\LR{Q^TG\,Q}}Q^T \\ }$$ The final task is to evaluate the $R$ matrix which lies at the heart of the theorem. This can be done using the log function and its derivative $\left\{\log(x),\frac{1}{x}\right\}$ evaluated at $B$, an all-ones matrix $J$, and $\z{X}$ which is an elementwise $\c{\sf zero\:indicator}$ function $$\eqalign{ \z{X}_{ij} &= \begin{cases} \o\qquad {\rm if}\;X_{ij}=0 \\ 0\qquad {\rm otherwise} \\ \end{cases} \\ Z &= \z{BJ-JB},\qquad L_\smB = \log(B),\qquad L_\smB' = B^{-1} \\ R &= {\frac{L_\smB J-JL_\smB+ZL_\smB'}{BJ-JB+Z}} \qquad \big({\rm Hadamard\;division}\big) \\\\ }$$ or in terms of the components of the $b$ vector $$\eqalign{ R_{jk} &= \begin{cases} {\Large\frac{\log(b_j)\,-\,\log(b_k)}{b_j\,-\,b_k}} \quad{\rm if}\; b_j\ne b_k \\ \\ \qquad\quad {\Large\frac{\o}{b_k}} \qquad\qquad {\rm otherwise} \\ \end{cases}\\ \\ }$$

---

Note that some of the steps above use the Frobenius product $(A:B)$, which is an extremely useful product notation for the trace function $$\eqalign{ A:B &= \trace{A^TB} \;=\; \sum_{i=1}^m\sum_{j=1}^n A_{ij}B_{ij} \\ A:A &= \|A\|^2_F \qquad \big({\rm hence\;the\;name}\big) \\ }$$ Note that the Frobenius and Hadamard products commute $$\eqalign{ A:(B\odot C) \;=\; (A\odot B):C \;=\; \sum_{i=1}^m\sum_{j=1}^n A_{ij}B_{ij}C_{ij} \\ }$$ There are also easily derived rules for rearranging product terms, e.g. $$\eqalign{ A\odot B &= B\odot A \\ A:B &= B:A \\ (AC^T):B &= A:(BC) = (B^TA):C \\ }$$

Answer 2

This could be derived step by step, using chain rules

1. Frobenious norm. Let's deal with the squared norm first! taking sqrt is easy. $$ d(A,B)^2=\|A-B\|^2_F = tr((A-B)^T(A-B)) = \sum_{ij} (A_{ij}-B_{ij})^2 $$ So if one of the matrices $A$ is a function of a scalar $t$ $$ \frac{\partial d(A(t),B)^2}{\partial t} = 2\sum_{ij} (A_{ij}-B_{ij})\frac{\partial A_{ij}}{\partial t}\\ =2tr((A-B)^T\frac{\partial A}{\partial t}) $$ One specific case is $$ \frac{\partial d(A,B)^2}{\partial A_{ij}} = 2(A_{ij}-B_{ij})\\ $$ So we can write $$\frac{\partial d(A,B)^2}{\partial A}=2(A-B)$$

Getting rid of the square $$ \frac{\partial d(A(t),B)}{\partial t} = \frac{1}{2} d^{-1/2}\cdot 2tr((A-B)^T\frac{\partial A}{\partial t}) = d^{-1/2}tr((A-B)^T\frac{\partial A}{\partial t}) $$

2. matrix logrithm (See Derivative of matrix logarithm) $$ \frac{\partial \log(X(t))}{\partial t} = X(t)^{-1} \frac{\partial X(t)}{\partial t} $$ We can think of the entry $X_{ij}$ as a parameter of $X$, then we get $$\frac{\partial X(t)}{\partial X_{ij}}=E_{ij}$$ .$E_{ij}$ is a matrix which only $i,j$ element is 1, others are all $0$.

3. Using chain rule

$$ \frac{\partial d(X,A_0)^2}{\partial X_{ij}} =\partial_{X_{ij}} \|\log X -\log A_0\|^2_F \\ = 2tr\left((\log X -\log A_0)^T \frac{\partial \log X}{\partial X}\frac{\partial X} {\partial X_{ij}}\right)\\ =2 tr\left((\log X -\log A_0)^T (X^{-1}E_{ij})\right) $$

Note that the effect $E_{ij}$ in the trace operator is quite unique Since $$ tr(A^TB) = \sum{ij} A_{ij}B_{ij}\\ tr(M^TE_{kl}) = \sum{ij} M_{ij}E_{ij} = M_{kl} $$ $E_{ij}$ works as a delta function in getting the entry.

Thus taking the element wise derivative together and form a matrix $$ \left(\frac{\partial d(X,A_0)^2}{\partial X}\right)_{ij} = 2 tr((\log X -\log A_0)^T (X^{-1}E_{ij}))\\ \frac{\partial d(X,A_0)^2}{\partial X} = 2 X^{-1}(\log X -\log A_0)\\ \frac{\partial d(X,A_0)}{\partial X} = d^{-1/2} X^{-1}(\log X -\log A_0) $$