Q: How to prove this Harmonic numbers identity? 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?

A: Your formula is correct.   The algebra is kind of messy, but here it is:
\begin{align}
\sum_{i=0}^{n-1}H_i&=\sum_{i=0}^{n-1} \sum_{j=1}^i \frac1{j}
\\&=\sum_{j=1}^{n-1}\frac{n-j}{j} \quad\quad\quad\scriptsize{\text{(since }\frac1{j}\text{ occurs in the above sum }n-j\text{ times)}}
\\&=n\sum_{j=1}^{n-1}\frac1{j}-\sum_{j=1}^{n}1
\\&=nH_{n-1}-n+1.
\\
\\
\sum_{i=0}^{n-1}H_i^{\,2}&=\sum_{i=0}^{n-1}\left( \sum_{j=1}^i \frac1{j} \cdot\sum_{k=1}^i \frac1{k}\right)
\\&=\sum_{\substack{0\le i \le n-1\\1\le j \le i\\1 \le k \le i}}\frac1{jk}
\\&=\sum_{j=1}^{n-1}\sum_{k=1}^{n-1} \frac{n-\max(j,k)}{jk} \quad\quad\quad\scriptsize{\text{(since }\frac1{jk}\text{ occurs in the above sum }n-\max(j,k)\text{ times)}}
\\&=\sum_{j=1}^{n-1}\left(\frac1{j}\left(\sum_{k=1}^{j-1}\frac{n-j}{k}+\sum_{k=j}^{n-1}\frac{n-k}{k} \right) \right)
\\&=\sum_{j=1}^{n-1}\left(\frac1{j}\left((n-j)\,H_{j-1}+\sum_{k=j}^{n-1}\frac{n}{k} - \sum_{k=j}^{n-1}1\right) \right)
\\&=\sum_{j=1}^{n-1}\frac{(n-j)\,H_{j-1}+n(H_{n-1}-H_{j-1}) - (n-j) }{j}
\\&=\sum_{j=1}^{n-1}\frac{-jH_{j-1}+nH_{n-1} - (n-j)}{j}
\\&=\sum_{j=1}^{n-1}\frac{nH_{n-1} - n}{j} -
\sum_{j=1}^{n-1}\frac{jH_{j-1} - j}{j}
\\&=n(H_{n-1}-1)H_{n-1}-\sum_{j=1}^{n-1}H_{j-1}+n-1
\\&=n(H_{n-1}-1)H_{n-1}-\left((n-1)H_{n-2}-(n-1)+1\right)+n-1\quad\quad\quad\scriptsize\text{(by the first formula)}
\\&=n(H_{n-1}-1)H_{n-1}-(n-1)H_{n-2}+2n-3
\\&=n H_{n-1}^{\,2}-nH_{n-1}-(n-1)(H_{n-1}-\frac1{n-1})+2n-3
\\&=n H_{n-1}^{\,2}-nH_{n-1}-(n-1)H_{n-1}+(n-1)\frac1{n-1}+2n-3
\\&=n H_{n-1}^{\,2}+(1-2n)H_{n-1}+2n-2
\\&=H_{n-1}(nH_{n-1}+1-2n)+2n-2.
\end{align}
It follows that
\begin{align}
\sum_{i=0}^{n-1}{H_i}^2 - \frac{1}{n}\left(\sum_{i=0}^{n-1}H_i\right)^2 &=H_{n-1}(nH_{n-1}+1-2n)+2n-2-\frac1{n}(nH_{n-1}-n+1)^2
\\&=nH_{n-1}^{\,2} +H_{n-1}(1-2n)+2n-2-\frac{n^2H_{n-1}^{\,2}+n^2+1-2n^2H_{n-1}+2nH_{n-1}-2n}{n}
\\&=nH_{n-1}^{\,2} +H_{n-1}(1-2n)+2n-2-nH_{n-1}^{\,2}-n-\frac1{n}+2nH_{n-1}-2H_{n-1}+2
\\&=H_{n-1}(1-2n+2n-2) +2n-2 -n-\frac1{n}+2
\\&=n-H_{n-1}-\frac1{n}
\\&=n-H_n.
\end{align}

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