Lipschitz constant of a function of matrix

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$.

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$).