While answering a question involving Harmonic numbers $H_n$, I wanted to simplify the terms $$f_n = \sum_{i=0}^{n-1}{H_i}^2 - \frac{1}{n}\sum_{0\le i,j\le n-1}H_i H_j. $$ To do so, I used SageMath to compute $f_n$ for $1\le n \le 100$ and then found the numerators of the resulting sequence to be OEIS sequence A187487; i.e., the $n$th numerator is the numerator of $n-H_n$. Indeed, all cases computed are consistent with the following being an identity for $n\ge 1$:
$$\sum_{i=0}^{n-1}{H_i}^2 - \frac{1}{n}\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}H_i H_j = n - H_n $$ which is the same as $$\sum_{i=0}^{n-1}{H_i}^2 - \frac{1}{n}\left(\sum_{i=0}^{n-1}H_i\right)^2 = n - H_n. $$ Does anyone have an idea of how to prove this? Or a reference?
After looking at this some more, I think a better way of approaching it is via summation by parts (or Abel's lemma), the discrete analog of integration by parts. This isn't Abel's formula in full generality, but it's all we need, and it's easy to verify:
$$\sum_{k=0}^{n-1} a_k = na_n-\sum_{k=1}^{n} k(a_{k}-a_{k-1}).$$
From this we can compute
\begin{align}\sum_{k=0}^{n-1} H_k &= nH_n-\sum_{k=1}^{n} k(H_{k}-H_{k-1}) \\&=nH_n-\sum_{k=1}^{n} k \frac1{k} \\&=nH_n-n \\&=n(H_n-1) \end{align}
and
\begin{align}\sum_{k=0}^{n-1} H_k^{\,2} &= nH_n^{\,2}-\sum_{k=1}^{n} k(H_{k}^{\,2}-H_{k-1}^{\,2}) \\&=nH_n^{\,2}-\sum_{k=1}^{n}k\left((H_{k-1}+\frac1{k})^2 -H_{k-1}^{\,2}\right) \\&=nH_n^{\,2}-\sum_{k=1}^{n}k\left(\frac1{k^2}+\frac{2H_{k-1}}{k} \right) \\&=nH_n^{\,2}-\sum_{k=1}^{n}\big(\frac1{k}+2H_{k-1}\big) \\&=nH_n^{\,2}-\sum_{k=1}^{n}\frac1{k}-2\sum_{k=1}^{n}H_{k-1} \\&=nH_n^{\,2}-H_n-2\sum_{k=0}^{n-1}H_{k} \\&=nH_n^{\,2}-H_n-2n(H_n-1) \\&=\frac{\big(n(H_n-1)\big)^2}{n}+n-H_n \\&=\frac1{n}\left(\sum_{k=0}^{n-1}H_k\right)^2+n-H_n, \end{align}
so
$$ \sum_{k=0}^{n-1} H_k^{\,2} - \frac1{n}\left(\sum_{k=0}^{n-1}H_k\right)^2=n-H_n,$$
as desired.
By the way, essentially the same argument, but using integration by parts instead of Abel's lemma, yields the formulas
$$\int_e^x (\ln t)^2 \,dt-\frac1{x}\left(\int_e^x \ln t\; dt\right)^2=x-e$$
and
$$\int_1^x (\ln t)^2 \,dt-\frac1{x}\left(\int_1^x \ln t\; dt-1\right)^2=x-2.$$