Suppose $X_i , X_j\in \mathbb{R}^d$ are gaussian vectors and $A$ is an $n\times n$ symmetric PSD matrix where $A_{ij} = f(\|X_i-X_j\|_2), \quad i,j\in 1,\ldots,n\;$ for some non-negative Lipschitz function $f$. Suppose we know that $A_{ij}$ concentrates around its mean, i.e. $\mathbb{P}(|A_{ij}-\mathbb{E}[A_{ij}]| > \delta)\leq \epsilon$, I'm trying to prove that $A$ concentrates around $E[A]$ in operator norm (looking for a tight bound). Here is my attempt:

\begin{align*} \mathbb{P}(\|A - \mathbb{E}[A]\|_{\text{op}} \geq t) & \leq \mathbb{P}(\|A-\mathbb{E}[A]\|_{\infty} \geq t)\\ & = \mathbb{P}(\max_i \sum_{j\neq i} |A_{ij}-\mathbb{E}[A_{ij}]| \geq t) \end{align*}

I'm not sure how to proceed from here. Obviously the entries in each row are (weakly) dependent by definition of $A_{ij}$. But can we assume $B_i := \sum_{j\neq i} |A_{ij}-\mathbb{E}[A_{ij}]| \quad \forall i=1,\ldots n$ are independent? If not, how to bound the last term without independence assumption? Any other approach to prove the matrix concentration is also appreciated.