How can I calculate complicated infinite sums with FindIntegerNullVector (or related methods)?

Question

I've recently been very interested in the wonderfully complex world of Euler sums, i.e. (convergent) infinite sums that, roughly speaking, consist of some rational polynomial combination of generalized harmonic numbers $H_p^{(q)} := \sum\limits_{k = 1}^{p}\frac{1}{k^q}$ and the parameters itself. For example $$\sum\limits_{k = 1}^{\infty} \frac{H_k}{k^2}, \sum\limits_{k = 1}^{\infty} \frac{H_k^3H_k^{(2)}}{k^2}, \sum\limits_{k = 1}^{\infty} \frac{H_k H_k^{(2)}H_k^{(4)}}{(k+1)^2}$$

are all convergent sums of this kind that are (with some exceptions) expressible as some polynomial combination of integer $\zeta$-values. Since I can't reasonably sink time into understanding all of the details of doing this for the cases I want to generalize, I turned to Mathematica in the hope that I could somehow squeeze out closed forms only dependent on zeta-values with some educated guessing for $\mathbb{Z}$-basis elements and the function FindIntegerNullVector.

I googled around a bit on how to apply this to my problem; aside from supplying sufficiently high precision of all involved numbers (done that up to 150 digits) and sorting the numbers to be related in certain ways, it seems to be wildly unreliable in outputting simple combinations in more complicated cases that Mathematica can't evaluate directly or are dependent on more than three combinations of $\zeta$-values.

For example:

Take the sum $\sum_{k = 1}^{\infty} \frac{H_k^3}{k^2} = 10\zeta(5) + \zeta(2)\zeta(3)$ - let's pretend we don't know this yet. Having worked with some Euler sums before, you might suspect that $\{\zeta(5),\zeta(2)\zeta(3),1\}$ is indeed a reasonable combination to express this sum as a linear combination of those, and in fact

In[*] := FindIntegerNullVector[N[{NSum[HarmonicNumber[k]^3/k^2,{k,1,Infinity},WorkingPrecision->20],Zeta[5],Zeta[2]Zeta[3],1},20]]
Out[*]= {1,-10,-1,0}

corresponding to the exact value of the sum if we take 20 digits of precision. Let's double the precision:

In[*] := FindIntegerNullVector[N[{NSum[HarmonicNumber[k]^3/k^2,{k,1,Infinity},WorkingPrecision->40],Zeta[5],Zeta[2]Zeta[3],1},40]]
Out[*]= {6293450745760373317489774,-62934507457603733174897742,-6293450745760373317489777,8} 

Ouch. I know that some sort of "overfitting" can occur when using Integer Relation algorithms, but how could one prevent this from happening? To make matters worse, let's look at another example.

Take the sum $\sum\limits_{k = 1}^{\infty} \frac{H_k^3}{k^6} = \frac{521}{24}\zeta(9)+3\zeta(2)\zeta(7)-\frac{97}{8}\zeta(3)\zeta(6) - \frac{51}{4}\zeta(4)\zeta(5) + 2\zeta(3)^3$. Again, with some educated guessing you might suspect the individual irrational terms to be a combination of this sum, and running the same calculations again

In[*] := FindIntegerNullVector[N[{NSum[HarmonicNumber[k]^3/k^6,{k,1,Infinity},WorkingPrecision->20],Zeta[9],Zeta[2]*Zeta[7],Zeta[3]*Zeta[6],Zeta[4]*Zeta[5],Zeta[3]^3},20]]
Out[*]= {55,42,-72,10,60,-35} 

and

In[*] := FindIntegerNullVector[N[{NSum[HarmonicNumber[k]^3/k^6,{k,1,Infinity},WorkingPrecision->40],Zeta[9],Zeta[2]*Zeta[7],Zeta[3]*Zeta[6],Zeta[4]*Zeta[5],Zeta[3]^3},40]]
Out[*]= {15349009,-24000024,-3069275,2446852,-14189761,14812626} 

where coefficients of both of these results are very far from their exact value, i.e. some "simpleness" assumptions would also be fruitless in this case.

Now why do I suspect that it should be theoretically possible to find closed forms of this kind with some integer relation algorithm? I actually got my inspiration for this from Vladimir Reshetnikov over at Math.SE, who has explained his out-of-the-blue thousands-of-digits-correct conjectures of some integrals/sums based on these relation algorithms* - which makes the poor performance of this PSLQ-function in this rather simple case... bizarre.

*See here and here. I do not actually have a link to it, but I vaguely recall him using the term "PSLQ" and presenting a high-caliber approximation in the same post.

My two questions are:

1. How can I avoid "overfitting" (unreasonably large coefficients) when choosing a precision for my problem?
2. What is the trick -other than educated guessing of the irrationals involved- that people like Reshetnikov use to obtain these obscenely good closed forms of really, really complicated expressions? Would this also work for this problem concerning Euler sums? If not, are there alternatives to PSLQ better fit for it?

Answer 1

$Version

(* "13.2.1 for Mac OS X ARM (64-bit) (January 27, 2023)" *)

Clear["Global`*"]

To get high precision with NSum, in addition to increasing the WorkingPrecision, you should use more than the default value for NSumTerms

Options[NSum, NSumTerms]

(* {NSumTerms -> 15} *)

Prepend[(table = Flatten[Outer[{#1, #2,
        NSum[HarmonicNumber[k]^3/k^2, {k, 1, Infinity},
          NSumTerms -> #1, WorkingPrecision -> #2] -
         (10 Zeta[5] + Zeta[2] Zeta[3])} &,
      {15, 30, 60}, {20, 40, 60}], 1]),
  {"#terms", "prec", "error"}] //
 Grid[#, Frame -> All] &

Prepend[(table = Flatten[Outer[{#1, #2,
        NSum[HarmonicNumber[k]^3/k^6, {k, 1, Infinity},
          NSumTerms -> #1, WorkingPrecision -> #2] -
         (521/24 Zeta[9] + 3 Zeta[2] Zeta[7] - 97/8 Zeta[3] Zeta[6] -
           51/4 Zeta[4] Zeta[5] + 2 Zeta[3]^3)} &,
      {15, 30, 60}, {20, 40, 60}], 1]),
  {"#terms", "prec", "error"}] //
 Grid[#, Frame -> All] &

Answer 2

Here is a way to get at least 50 decimal places. One can boost this amount modestly (I did not try to quantify the error term too closely and skip that here).

First sum 10000 terms directly.

ee = HarmonicNumber[k]^3/k^2;
ss = Series[ee, {k, Infinity, 12}];
sum10k = Sum[ee, {k, 10000}];
nsum10k = N[sum10k, 50]

(* Out[901]= 12.217613264921967590943295626274344756881149234350 *)

To estimate the rest we sum the series expansion termwise from 1 to infinity (this can be done symbolically). We'll add this to the value above. Then we brute-force sum the first 10000 terms of the series expansions and subtract that from the sum, to account for the fact that those terms were double-counted.

seriesterms = ss[[3]];
seriestermsums = 
  Table[Sum[seriesterms[[j]]/k^(1 + j), {k, 1, Infinity}], {j, 
    Length[seriesterms]}];
seriestermsumstot = Total[N[seriestermsums, 50]]
(* Out[904]= 12.347098083612797705384266996819903504057847643218 *)

seriestermssumsto10k = 
  Table[Sum[seriesterms[[j]]/k^(1 + j), {k, 1, 10000}], {j, 
    Length[seriesterms]}];
totalsum10k = Total[N[seriestermssumsto10k, 50]]

(* Out[908]= 12.218129446803769914816822317038781008286036215889 *)

approxsum = nsum10k + seriestermsumstot - totalsum10k

(* Out[910]= 12.346581901730995381510740306055467252652960661679 *)

Compare to the proposed value.

N[10*Zeta[5] + Zeta[2]*Zeta[3], 50]

(* Out[867]= 12.346581901730995381510740306055467252652960661679 *)

I would expect that this approach might get one 100 digits without too much trouble. Maybe more, but I think 1000 would be out of reach and require a finer analysis and different (read: better) ideas.