TL;DR: Is there a way to make Schur's (elegant) proof of Hilbert's inequality feel like less of a trick/miracle?
Longer version: Let me go quickly over Schur's proof to show what I mean. Actually, let me do it for Montgomery and Vaughan's generalized Hilbert inequality, since the idea is the same, and let me also change it a little, so as to have fewer sums floating around. (Perhaps I've already taken the first step towards answering my question by asking it.)
Let $\{r_i\}_{i\in I}$, $r_i$ real, $I$ an index set, be given. Our task is to bound the norm of the operator $A = (a_{i,j})_{i,j\in I}$ given by $$a_{i,j} = \frac{1}{r_i-r_j}$$ for $i\ne j$, $a_{i,i}=0$, under the assumption that $|r_i-r_j|\geq \delta$ for all distinct $i,j\in I$ and some $\delta>0$.
Write $A^2 = (b_{i,j})_{i,j\in I}$. Of course $b_{i,j} = \sum_{k\ne i,j} a_{i,k} a_{k,j}$, and so, for $i\ne j$, $$\begin{aligned}b_{i,j} &= \sum_{k\ne i,j} \frac{1}{(r_i-r_k) (r_k - r_j)} = \frac{1}{r_i-r_j} \sum_{k\ne i,j} \left(\frac{1}{r_i-r_k} + \frac{1}{r_k - r_j}\right)\\ &=\frac{1}{r_i-r_j} \left(\sum_{k\ne i} \frac{1}{r_i-r_k} + \sum_{k\ne j}\frac{1}{r_k - r_j}\right) - \frac{2}{(r_i-r_j)^2} \end{aligned}$$ (the second equality here is the tiny miracle) whereas $b_{i,i} = - \sum_k 1/(r_i-r_k)^2$. In other words, $$A^2 = J A - A J - D - S,$$ where $J$ and $D$ are diagonal matrices with $J_{i,i}=\sum_{k\ne i} 1/(r_i-r_k)$ and $D_{i,i} = \sum_{k\ne i} 1/(r_i-r_k)^2$, and $S$ is a matrix with $S_{i,j} = 1/(r_i-r_j)^2$ for $i\ne j$ and $S_{i,i}=0$.
Because $A$ is real antisymmetric, the spectral theorem holds. Hence, to bound the norm of $A$, it is enough to bound $|A v|_2/|v|_2$ for all eigenvectors $v$ of $A$. Again because $A$ is real antisymmetric, all of its eigenvalues are pure imaginary. Hence, for $v$ an eigenvector with eigenvalue $\lambda$, $$\langle v, J A v\rangle - \langle v, A J v\rangle = \lambda \langle v, J v\rangle + \langle A v, J v\rangle = (\lambda + \overline{\lambda}) \langle v, J v\rangle = 0,$$ and so $$|A v|_2^2 = -\langle v, A^2 v\rangle = \langle v, D v\rangle + \langle v, S v\rangle.$$
The norm of $D$ is clearly $\max_i |D_{i,i}| = \max_i \sum_{k\ne i} 1/(r_i-r_k)^2 \leq 2 \zeta(2)/\delta^2 = \pi^2/3 \delta^2$, whereas the norm of $S$ is $\leq \max_i \sum_j |S_{i,j}| = \max_i \sum_{j\ne i} 1/(r_i-r_j)^2\leq \pi^2/3 \delta^2$. Hence, $$|A v|_2^2 \leq \frac{\pi^2}{\delta^2} |v|_2^2$$ for any eigenvector $v$ of $A$. We conclude that the norm of $A$ is $$\leq \frac{\pi}{\delta}.$$
So, questions:
(a) Can this variant of the proof be shortened further, or made more natural?
(b) The above was closely inspired by Montgomery-Vaughan and the exposition in Iwaniec-Kowalski. I take the above may literally appear elsewhere?
(c) Are there other interesting matrices for which we have expressions such as $A^2 = J A - A J - D - S,$ with similar consequences?