Lipschitz constant of a function of matrix The function is given by
$f(X) = (AX^{-1}A^\top + B)^{-1}$ where $X$, $A$, and $B$ are $n \times n$ positive definite matrices. 
I'm trying to find the Lipschitz constant such that $\| f(X)-f(Y) \| \leq L \|X-Y\|$ where $X \geq 0$ and $Y \geq 0$. Motivated by Lemma 3.1 in Nonlinear Systems(H. Khalil, 3rd Ed.), I tried to find the derivative of $f(X)$ (i.e. $\| \frac{ \partial f(X)}{\partial X} \|$) but it's not easy to find the derivative of a function of a matrix over a matrix. 
How can I find the Lipschitz constant? or Please let me know if there exists the way to calculate the derivative of a function of a matrix.
 A: Let $h=Y-X$. Using the first order expansion of the matrix inverse,
$$
f(X+h)=(A(X^{-1}-X^{-1}hX^{-1})A^\top+B)^{-1}+O(\|h\|^2)
$$
Now let $Z=AX^{-1}A^\top+B$ and let $g = AX^{-1}hX^{-1}A^\top$. Then
$$
f(X+h)=(Z-g)^{-1}+O(\|h\|^2)=Z^{-1}+Z^{-1}gZ^{-1}+O(\|h\|^2).
$$
Since $Z^{-1}=f(X)$, it follows that 
$$
f(X+h)-f(X)=Z^{-1}AX^{-1}hX^{-1}A^\top Z^{-1}+O(\|h\|^2).
$$
Using $Z^{-1}AX^{-1}=(XA^{-1}Z)^{-1}=(A^T+XA^{-1}B)^{-1}$, we obtain that
$$
\|Z^{-1}AX^{-1}\|\leq \|A^{-1}\|.
$$
Similarly, $\|X^{-1}A^\top Z^{-1}\|\leq \|A^{-1}\|$.
Consequently,
$$
\|f(X+h)-f(X)\|\leq \|A^{-1}\|^2\|h\|+O(\|h\|^2),
$$
yielding a Lipschitz constant of $L=\|A^{-1}\|^2$.
A: I assume $X\ge0$ means $u^\top X u\ge0$, and that $B$ is definite positive $$\inf_{\|u\|=1} u^\top B\, u:=\beta>0.$$ I also assume matrix norms are the Euclidean operator norms.  
Compute the differential by the chain rule, as suggested in comments by F.Poloni:
$$Df(X)H=(AX^{-1}A^{\top}+B)^{-1}AX^{-1}\cdot H\cdot X^{-1}A^{\top}(AX^{-1}A^{\top}+B)^{-1}$$
$$=(A^{\top}+XA^{-1}B)^{-1}\cdot H \cdot(A+BA^{-\top}X)^{-1}=$$
$$=\big(B^\top+Y\big)^{-1}B^\top A^{-\top}\cdot H\cdot A^{-1}B^\top\big(B^\top+Z)^{-1}, $$
where  $Y:=B^{\top}A^{-\top} X A^{-1}B\ge0 $ and $Z:=BA^{-\top}XA^{-1}B^\top\ge0$, conjugated to $X\ge0$. Thus for any unit norm vector $u\in\mathbb{R}^n $
$$\big\|\big(B^\top+Y\big)u\big\|\ge u^\top\big(B^\top+Y\big)u\ge u^\top B  u \ge\beta$$ and
$$\big\|\big(B^\top+Z)u\big\|\ge u^\top\big(B^\top+Z)u \ge u^\top B  u \ge\beta.$$
Hence 
$$\big\|\big(B^\top+Y\big)^{-1}\big\|\le \beta^{-1}$$ and
$$\big\|\big(B^\top+Z)^{-1}\big\|\le \beta^{-1}.$$
Therefore 
$$\|Df\|_\infty\le\|B\|^2\|A^{-1}\|^2\beta^{-2}$$
which is also a Lipschitz constant for $f$, since its domain is convex, $\{X\ge0\}$.
$$*$$
Rmk. The above bounds on the $L_2$ operator norms (or others matrix norms) could be improved, but not up to  $$\big\|\big(A^{\top}+XA^{-1}B\big)^{-1}\big\|\le\big\|A^{-1}\big\|,$$ 
even for symmetric definite positive matrices. 
Take  e.g. $n=2$ and 
$$A=:I=\begin{bmatrix} 
1  & 0 \\
0 & 1
\end{bmatrix}\quad X:=\begin{bmatrix} 
5/2 & 1 \\
1 & 1/2 
\end{bmatrix}\quad B:=\begin{bmatrix} 
1 & -1/2 \\
-1/2 & 1/2 
\end{bmatrix}$$
Then 
$$\big(A^{\top}+XA^{-1}B\big)^{-1}=\big(I+XB\big)^{-1}= \begin{bmatrix} 
4/15 & 4/15\\
-4/15 & 16/15 
\end{bmatrix}$$
whose maximum singular value is $\displaystyle{2\over 15}\sqrt{29}+{2\over 5}>1=\|A^{-1}\|$. The same holds for the Frobenius, and other common entry-wise norms (due to the coefficient $16/15>1$).
