Return to Answer

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) 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, or other usual entry-wise norms.

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

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

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

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) 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, or other usual entry-wise norms.

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

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

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

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