Bounding the norm of the product of random PSD matrices

Question

Consider the following setup,

1. $(X, \hat{X}, Y, \hat{Y})$ are four $n \times n$ real, symmetric, full-rank, positive-definite matrices with entries between zero and one and operator norm $O(n)$.
2. The matrix entries are $X_{ij} = k(x_i,x_j)$, $\hat{X}_{ij} = k(\hat{x}_i, \hat{x}_j)$, $Y_{ij} = k(y_i,y_j)$, and $\hat{Y}_{ij} = k(\hat{y}_i, \hat{y}_j)$, where $k : \mathbb{R}^d \times \mathbb{R}^d \mapsto \mathbb{R}$ is a smooth, positive-definite kernel function such that $|k(a,b)| \leq 1$ for all $a$, $b$.
3. $(\hat{X}, \hat{Y})$ are random and such that $\mathbb{E}[\hat{X}] = X$, $\mathbb{E}[\hat{Y}] = Y$.
4. I have $O(n\log n)$ upper bounds on the expected operator norms $\mathbb{E}[\| \hat{X} - X\|]$ and $\mathbb{E}[\| \hat{Y}-Y\|]$.

My question is: Given $\gamma > 0$, is

$$\mathbb{E}[\|(\hat{X}+\gamma I)^{-1}\hat{Y}-(X + \gamma I)^{-1}Y\|] \leq C_x\, \mathbb{E}[\| \hat{X} - X\|] + C_{y}\, \mathbb{E}[\| \hat{Y}-Y\|] = O(n\log n)$$

for some pair of constants $C_x$, $C_y$, independent from $n$? Numerical experiments reveal that this might be the case.