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.



Lipschitz constant of a function of matrix
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$).


