Strategies For Summing Harmonic Numbers Lately, I have found several interesting problems involving Harmonic numbers such as \begin{equation*}\sum_{n=1}^{\infty}\frac{H_n^{(2)}}{n^2}=\frac{7\pi^4}{360}\end{equation*}
I am not familiar with computing sums involving Harmonic numbers. Is there a general strategy for tackling such problems? 
How can this series be evaluated by operating on the generating function for $H_n^{(2)}$?
 A: Here is a general version I discovered some time ago.

Proposition : $$\sum_{r=1}^{n} \dfrac{H_{r} ^{(m)}}{r^m} = \dfrac{1}{2} \left( [H_{n}^{(m)}]^2 + H_{n}^{(2m)} \right) \quad ; \quad m \geq 1$$   

Proof : Expand the summation, using the definition of $\displaystyle H_{r}^{(m)}$ , to see that $$\displaystyle \sum_{r=1}^{n} \dfrac{H_{r}^{(m)}}{r^m} = \sum_{r=1}^{n} \dfrac{1}{r^{2m}} + \sum_{r < j} \sum_{j=2}^{n} \dfrac{1}{rj} = H_{n}^{(2m)} + \sum_{r < j} \sum_{j=2}^{n} \dfrac{1}{rj}  $$    
And,  
$\displaystyle [H_{n}^{(m)}]^2 = H_{n}^{(2m)} + 2\sum_{r < j} \sum_{j=2}^{n} \dfrac{1}{rj} $   
Eliminating $\displaystyle \sum_{r < j} \sum_{j=2}^{n} \dfrac{1}{rj}$ from the above equations, we have,    
$$\displaystyle \sum_{r=1}^{n} \dfrac{H_{r}^{(m)}}{r^m} =\dfrac{1}{2} \left( [H_{n}^{(m)}]^2 + H_{n}^{(2m)} \right) \quad \square \tag{*} $$   
Now put $m=2$ take limit to infinity to get the desired result.   
We can also derive the result using Summation By Parts.
A: You may evaluate it without using generating functions. 
$$S=\sum\limits_{n=1}^{\infty} \frac{H_n^{(2)}}{n^2} = \sum\limits_{n=1}^{\infty}\sum\limits_{k=1}^{n}\frac{1}{n^2k^2}$$
By changing the order of summation, you may write it as:
$$\begin{align}S &= \sum\limits_{k=1}^{\infty}\sum\limits_{n=k}^{\infty}\frac{1}{n^2k^2}\\&= \sum\limits_{k=1}^{\infty}\frac{1}{k^2}\left(\sum\limits_{n=k}^{\infty}\frac{1}{n^2}\right)\\&= \sum\limits_{k=1}^{\infty}\frac{1}{k^2}\left(\sum\limits_{n=1}^{\infty}\frac{1}{n^2} - \sum\limits_{n=1}^{k-1}\frac{1}{n^2}\right) \\&=  \sum\limits_{k=1}^{\infty}\frac{1}{k^2}\left(\zeta(2) - \sum\limits_{n=1}^{k-1}\frac{1}{n^2}\right)\\&= \sum\limits_{k=1}^{\infty}\frac{1}{k^2}\left(\zeta(2) + \frac{1}{k^2} - H_k^{(2)}\right)\\&= \zeta(2)\sum\limits_{k=1}^{\infty}\frac{1}{k^2} + \sum\limits_{k=1}^{\infty}\frac{1}{k^4} - \sum\limits_{k=1}^{\infty}\frac{H_k^{(2)}}{k^2} \\&= \zeta^2(2)+\zeta(4) - S\end{align}$$
Hence, $\displaystyle S = \frac{\zeta^2(2)+\zeta(4)}{2}$.
In general for any sequence $(a_n)_{n \ge 1}$ such that $\displaystyle \sum_{n=1}^{\infty} a_n$ converges absolutely,
We have $$\sum_{n=1}^{\infty}\left(a_n\sum_{k=1}^{n}a_k\right) = \frac{1}{2}\left(\sum_{n=1}^{\infty}a_n^2 + \left(\sum_{n=1}^{\infty}a_n\right)^2\right)$$ by using the same method as above.
A: By using the formula in the answer to question Calculating alternating Euler sums of odd powers we have:
\begin{equation}
{\bf H}_2^{(2)}(t) = Li_4(t) + \frac{1}{2} [Li_2(t)]^2  +\underbrace{\int\limits_0^1 \frac{\log(\xi)}{\xi} [\log(1-t \xi)]^2 d\xi}_{I(t)}
\end{equation}
Now from Compute an integral containing a product of powers of logarithms. we have:
\begin{equation}
I(1) = \frac{1}{2}\left( -3! \zeta(4) + 2[\zeta(2)]^2\right)
\end{equation}
where $\Psi$ is the polygamma function. If we want to find the value at $t=-1$ it is more complicated. Here we use the identity $\log(1+\xi) = \log(1-\xi^2) - \log(1-\xi)$ and change variable in one of the integrals accordingly and use the approach developed in Calculating alternating Euler sums of odd powers. Then we have:
\begin{eqnarray}
I(-1) &=& \left(\frac{1}{4}-1\right) \cdot \int\limits_0^1 \frac{\log(\xi)}{\xi} \log(1-\xi)^2 d\xi - 2 \int\limits_0^1 \frac{\log(\xi)}{\xi}  \log(1-\xi) \log(1+\xi) d\xi\\
I(-1)&=&\left(\frac{1}{4}-1\right) \cdot I(1) + 2 \left({\bf H}^{(1)}_3(-1) + {\bf H}_2^{(2)}(-1) +\frac{1}{2} \zeta(2)^2\right)
\end{eqnarray}
Bringing all together we have:
\begin{eqnarray}
{\bf H}^{(2)}_2(1) &=& \frac{7 \pi^4}{360}\\
{\bf H}^{(2)}_2(-1) &=& - \frac{37 \pi^4}{1440} - 2 {\bf H}^{(1)}_3(-1)
\end{eqnarray}
