A double series $\frac13 \sum_{j=1}^{\infty}\sum_{i=1}^{\infty}\frac{(i-1)! (j-1)!}{(i+j)!}H_{i+j}$ giving $\zeta(3)$ Here is a symmetric rational double series giving Apery's constant:

$$
\frac13 \sum_{j=1}^{\infty}\sum_{i=1}^{\infty} \displaystyle \frac{(i-1)! (j-1)!}{(i+j)!}  H_{i+j} = \zeta(3)
$$

where $\displaystyle H_{n}:=\sum_{1}^{n} \frac{1}{k}$, $n=1,2,\cdots,$ are the harmonic numbers. 
How would you prove it?
Edit. In 2005, I sent this result to Wolfram MathWorld (see equation 25), I built it as an echo of

$$
 \sum_{j=1}^{\infty}\sum_{i=1}^{\infty} \displaystyle \frac{(i-1)! (j-1)!}{(i+j)!}  = \zeta(2).
$$

See Jack's pretty answer and see my answer below.
 A: This is another path.

Theorem. Let $x$ be a real number such that $-1<x<1$ and set
  $$
Z(x):=\sum_{j=1}^{\infty}\sum_{i=1}^{\infty}\frac{(i-1)! (j-1)!}{(i+j)!}  x^{i+j}.
$$
  Then
  $$
Z(x)= 2 \:\mathrm{Li_2} (x) - \mathrm{Li_2} \left( x(2-x) \right) \tag1
$$

where $\mathrm{Li_2}$ is the dilogarithm function such that $\mathrm{Li_2}(x)=\displaystyle \sum_{n=1}^{\infty}\displaystyle \frac{x^{n}}{n^2}= - \displaystyle \int_{0}^{x}\,\displaystyle \frac{\log(1-t)}{t}\mathrm{d}t$.
Proof. Observe that 
$$
{\frac{(i-1)! (j-1)!}{(i+j)!}=  \frac{\Gamma(i) \Gamma(j)}{(i+j) \Gamma(i+j)} = \frac{1}{i+j}\,B(i,j)=
 \int_{0}^{1}\, \frac{t^{i-1}(1-t)^{j-1}}{i+j}\mathrm{d}t.}
$$
Let  $|x|<1$. We can obtain
$$
\begin{align*}
\sum_{j=1}^{\infty}\sum_{i=1}^{\infty}\frac{(i-1)! (j-1)!}{(i+j)!} x^{i+j} &= \sum_{j=1}^{\infty}\sum_{i=1}^{\infty}\displaystyle \int_{0}^{1} \frac{t^{i-1}(1-t)^{j-1}}{i+j}dt \:x^{i+j} 
\\
 &=  x^2  \int_{0}^{1}  \sum_{j=1}^{\infty}\,(x(1-t))^{j-1}\sum_{i=1}^{\infty}  \frac{(tx)^{i-1}}{i+j}\:\mathrm{d}t \\
 &= x^2  \int_{0}^{1}  \sum_{j=1}^{\infty}\,(x(1-t))^{j-1}(tx)^{-1-j}\sum_{i=j+1}^{\infty}  \frac{(tx)^{i}}{i}\:\mathrm{d}t\\
 &= x^2 \int_{0}^{1}  \sum_{j=1}^{\infty}\,(x(1-t))^{j-1}(tx)^{-1-j}\int_{0}^{tx}  \frac{u^{j}}{1-u}\:\mathrm{d}u\:\mathrm{d}t \\
 &= -x  \int_{0}^{1} \frac{\log(1-x + x u)}{1-x u} \:\mathrm{d}u \\
 &= -x  \int_{1/x - 1}^{1/x} \frac{ \log(1-\frac{x}{2-x }u) + \log(\frac{2-x}{x})}{u}  \:\mathrm{d}u\\
&= 2 \: \mathrm{Li_2} (x) - \mathrm{Li_2} \left( x(2-x) \right)
\end{align*}$$
Example 1. 
$$
\sum_{j=1}^{\infty}\sum_{i=1}^{\infty}\frac{(i-1)! (j-1)!}{(i+j)!}= \frac{\pi^2}{6}
$$ Put $x=1$ in $(1)$.
Example 2. 
$$
\sum_{j=1}^{\infty}\sum_{i=1}^{\infty}\frac{(i-1)! (j-1)!}{(i+j)!} \frac{1}{ \varphi^{2(i+j)}}= \frac{\pi^2}{30}-\ln^2 \varphi
$$ Put $x=1/\varphi^2$ in $(1)$, with $\varphi=\frac{1+\sqrt{5}}{2}$.

Proposition.
  $$
\sum_{j=1}^{\infty}\sum_{i=1}^{\infty} \displaystyle \frac{(i-1)! (j-1)!}{(i+j)!}  H_{i+j} = 3 \: \zeta(3) \tag2
$$

Proof. From the identity
$$
1+x+\cdots+x^{i+j-1} = \frac{1- x^{i+j}}{1 - x\,}, \qquad  i\geq1,\,j\geq1, 0<x<1,
$$
we write
$$
H_{i+j} = \displaystyle \int_{0}^{1} \displaystyle \frac{1- x^{i+j}}{1 - x\,}\: \mathrm{d}x. 
$$
Then, using $(1)$, we obtain
$$
\begin{align} \sum_{j=1}^{\infty}\sum_{i=1}^{\infty} \frac{(i-1)! (j-1)!}{(i+j)!} \: H_{i+j} & =\int_{0}^{1} \frac{ Z(1) - Z(x)}{1 - x}\: \mathrm{d}x 
\\ & =   \int_{0}^{1} \ln(1-x)\left(Z(1) -  2 \: \mathrm{Li_2}  (x) +  \mathrm{Li_2} \left( x(2-x) \right) \right)'\: \mathrm{d}x \\& =  2 \int_{0}^{1} \dfrac{\ln^2 (1-x)}{2-x}\: \mathrm{d}x 
\\ & =  2  \int_{0}^{1} \dfrac{\ln^2 x}{1+x}\: \mathrm{d}x \nonumber 
\\ & = 2 \int_{0}^{1}  \sum_{n=0}^{\infty}(-1)^{n} x^{n} \ln^2 x\: \mathrm{d}x
\\ & = 3 \: \zeta(3). \end{align}
$$
A: Summation by parts seems to work fine.
For $j=1$ we have:
$$\sum_{i=1}^{+\infty}\frac{H_{i+1}}{i(i+1)}=2,$$
since $\sum_{i=1}^{N}\frac{1}{i(i+1)}=1-\frac{1}{1+N}$ and so:
$$\begin{eqnarray*}\sum_{i=1}^{N}\frac{H_{i+1}}{i(i+1)}&=&H_{N+1}\left(1-\frac{1}{N+1}\right)-\sum_{i=1}^{N-1}\left(1-\frac{1}{i+1}\right)\frac{1}{i+2}\\&=&H_2+\sum_{i=1}^{N-1}\frac{1}{(i+1)(i+2)}+O\left(\frac{\log N}{N}\right)\\&=&H_2+\frac{1}{2}+O\left(\frac{\log N}{N}\right),\end{eqnarray*}$$
since:
$$\sum_{i=1}^{+\infty}\frac{1}{(i+1)\cdot\ldots\cdot(i+k)}=\frac{1}{(k-1)k!}.$$
For $j=2$ we have:
$$\begin{eqnarray*}1!\sum_{i=1}^{N}\frac{H_{i+2}}{i(i+1)(i+2)}&=&H_{N+1}\left(\frac{1}{4}-\frac{1}{2(N+1)(N+2)}\right)-\sum_{i=1}^{N-1}\left(\frac{1}{4}-\frac{1}{2(i+1)(i+2)}\right)\frac{1}{i+3}\\&=&\frac{H_3}{4}+\frac{1}{4\cdot 3!}+O\left(\frac{\log N}{N}\right)\end{eqnarray*}$$
and in the general case we have:
$$(j-1)!\sum_{i=1}^{+\infty}\frac{H_{i+j}}{i\cdot\ldots\cdot(i+j)}=(j-1)!\left(\frac{H_{j+1}}{j^2(j-1)!}+\frac{1}{j^2(j+1)!}\right)=\frac{H_{j+1}}{j^2}+\frac{1}{j^3(j+1)},$$
so the original sum equals:
$$\sum_{j=1}^{+\infty}\frac{H_{j+1}}{j^2}+\sum_{j=1}^{+\infty}\frac{1}{j^3(j+1)}=\sum_{j=1}^{+\infty}\frac{H_j}{j^2}+\sum_{j=1}^{+\infty}\frac{1}{j^3}=2\zeta(3)+\zeta(3)=3\zeta(3).$$
