# HarmonicNumber problem

I am looking for all the roots of HarmonicNumber, in the domain [-30, 1] and [0, 6000] for the real and imaginary parts, respectively, and where the parameter n of HarmonicNumber is as high as 200.

As I increase the size of the domains where I am looking for roots, the following warnings show up:

FindRoot::jsing: "Encountered a singular Jacobian at the point {x,y} = {56.9935,171.987}. Try perturbing the initial point(s)."

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances "

Could you help me fix this? Or suggest another way to solve this equation.

I got the same warnings when I tried to apply another technique as described in this answer

Here is the code:

Options[FindRoots2D] = {PlotPoints -> Automatic, MaxRecursion -> Automatic};

FindRoots2D[funcs_, {x_, a_, b_}, {y_, c_, d_}, opts___] := Module[
{fZero, seeds, signs, fy},fy = Compile[{x, y}, Evaluate[funcs[[2]]]];
fZero = Cases[Normal[
ContourPlot[
funcs[[1]] == 0,
{x, a-(b-a)/97, b+(b-a)/103}, {y, c-(d-c)/98, d+(d-c)/102},
Evaluate[FilterRules[{opts}, Options[ContourPlot]]]]],
Line[z_] :> z, Infinity];

seeds = Flatten[(
(signs = Sign[Apply[fy, #1, {1}]];
#1[[1 + Flatten[Position[Rest[signs*RotateRight[signs]], -1]]]]) &
) /@ fZero, 1];
If[seeds == {}, {},
Select[
Union[({x, y} /.
FindRoot[{funcs[[1]], funcs[[2]]}, {x, #1[[1]]}, {y, #1[[2]]},
Evaluate[FilterRules[{opts}, Options[FindRoot]]]] & ) /@ seeds,
SameTest -> (Norm[#1 - #2] < 10^(-6) & )],
a <= #1[[1]] <= b && c <= #1[[2]] <= d & ]]]
n = 20;
zeros = FindRoots2D[{Re[HarmonicNumber[n, x + I y]], Im[HarmonicNumber[n, x + I y]]}, {x, -30, 1}, {y, 0, 6000}]
mediumSeaGreen = RGBColor[0.235298, 0.702002, 0.443098];
violet = RGBColor[0.559999, 0.370006, 0.599994];
orangeRed = RGBColor[1., 0.270608, 0.];
ContourPlot[{Re[HarmonicNumber[n, x + I y]] == 0, Im[HarmonicNumber[n, x + I y]] == 0},
{x, -30, 1}, {y, 0, 6000},
MaxRecursion -> 3,ContourStyle -> {Directive[Thick, mediumSeaGreen],
Directive[Thick, violet]},
Epilog -> {PointSize[0.02], orangeRed, Point[zeros]}]


The following was my trial of n=20 and range x={-32,1} and y={0,4000}. n = 20; seeds = findSeeds[myHN[n, x + I y], {x, -32, 1}, {y, 0, 4000}, PlotPoints -> 25, MeshFunctions -> {myHNReImCx[n][#1, #2] &, myHNReImCy[n][#1, #2] &}]; // AbsoluteTiming allfoundzeros = Quiet@myFRC[myHN[n, z], {z, -32, 1 + 4000 I}, seeds.{1, I}]; // AbsoluteTiming uniqueszeros = #[[First@Ordering[Abs@myHN[n, #]]]] & /@ findclusters[allfoundzeros, 10.^-6]; // AbsoluteTiming zeros = Pick[uniqueszeros, With[{zz = DeveloperToPackedArray@N@uniqueszeros}, UnitStep[(#1 - (-32)) (1 - #1) (#2) (4000 - #2)] &[Re[zz], Im[zz], 1]], 1]; // AbsoluteTiming Length[zeros] zeroplot = Graphics[{PointSize[0.0035], Red, Point@Transpose@{Re[zeros], Im[zeros]}}, PlotRange -> {{-32, 1}, {0, 4000}}, AspectRatio -> 1, Frame -> True]

-
I ask you in all earnestness: do you think anybody can help you when don't show the code you have written. –  m_goldberg Jun 23 '14 at 8:03
I have updated the code. –  user16023 Jun 23 '14 at 16:36
Oh. No problem. Thanks for your reply. Enjoy your vacation. –  user16023 Jun 26 '14 at 23:11
–  Artes Jul 14 '14 at 23:06

The problems with the method seems to arise from two sources, the objective function (and its components Re, Im, HarmonicNumber) and the working precision (MachinePrecision).

The objective function

The Jacobian: The functions Re and Im do not have derivatives. Inspecting the Jacobian at {56.9935, 171.987} shows unevaluated derivatives:

D[{Re[HarmonicNumber[n, x + I y]], Im[HarmonicNumber[n, x + I y]]}, {{x, y}}] /.
((* {{(-4.7632*10^-18 - 8.08894*10^-19 I) Derivative[1][Re][1. + 1.16699*10^-18 I],... *)


The derivatives of HarmonicNumber evaluate to a complex number but the derivatives of Re and Im do not. In such a case, I believe FindRoot will use a finite difference method for computing the Jacobian. Insufficient working precision can lead to issues in computing the Jacobian, which is not the hint provided by FindRoot::jsing message. In particular, function evaluations that mathematically should produce distinct numbers might be equal at MachinePrecision and cancel out.

Precision: Mathematically, HarmonicNumber[n, z] and Sum[1/k^z, {k, n}] are equivalent. One would like to think that HarmonicNumber would be at least as precise and fast as the sum. This may be true over some domains. For small values of n (certainly up to 200 and more), the sum is more precise. In fact, over the domain of this problem (n <= 200, z in the complex rectangle from -30 + 0 I to 1 + 6000 I), the sum is much more precise. For n = 20 and 1000 random, arbitrary-precision complex numbers with precision \$MachinePrecision (15.95), the mean loss of precision for the sum was 4.3 digits and for HarmonicSum, it was 14.1 digits. For n = 20 and 1000 random, arbitrary-precision complex numbers with precision 100, the mean loss of precision for the sum was 4.3 digits and for HarmonicSum, it was 26.6 digits.

Speed: The sum is also considerable faster than HarmonicNumber for small values of n. At about n == 1800, the speed of the sum and HarmonicNumber are about equal; for larger values of n, HarmonicNumber is faster.

Compilability: HarmonicNumber is not compilable. This is not a serious problem, but the FindRoots2D code compiles the objective functions.

Working precision

MachinePrecision is good enough for many cases, but numerical errors can accumulate and (1) cause a Jacobian computed with finite differences to vanish erroneously or (2) confuse the root-finding algorithm with erratic function evaluations. These seem to be happening in this case. The first may be the reason for the FindRoot::jsing error; the second may be the reason for the FindRoot::lstol warning. For a function with isolated critical points, a FindRoot::jsing ought to be rare. A FindRoot::lstol warning can happen when FindRoot gets stuck near a local extremum, but I was not able to find such an example in the OP's problem after a (very) brief search.

### Some improvements

Increase WorkingPrecision

Setting WorkingPrecision -> 100 took care of the FindRoot::jsing errors. It did not take care of all the FindRoot::lstol warnings. However, it is much slower to use arbitrary-precision arithmetic on all solutions.

Use the explicit sum instead HarmonicNumber

As described above, using an explicit sum (or Plus) expression results in greater accuracy and speed. Below myHN is my version of HarmonicNumber. If real and imaginary parts are desired, then myHNReIm can be used. It can be compiled, which will result in much faster computations at MachinePrecision. (For a specific integer n, each definition generates the sum to be computed and memoizes it. In this way the Sum is expanded only once.)

myHN[n_Integer, z0_] := Block[{z},                  (* complex-valued *)
myHN[n, z_] := Evaluate@Sum[1/k^z, {k, n}];
myHN[n, z0]
];
myHNReIm[n_Integer, x0_, y0_] :=                    (* real and imaginary parts *)
Block[{x, y, k},
With[{reHN = Simplify[Re[1/k^(x + I y)] // ComplexExpand,
x ∈ Reals && y ∈ Reals && k > 0],
imHN = Simplify[Im[1/k^(x + I y)] // ComplexExpand,
x ∈ Reals && y ∈ Reals && k > 0]},
myHNReIm[n, x_, y_] := Evaluate@{Sum[reHN, {k, n}], Sum[imHN, {k, n}]};
myHNReIm[n, x0, y0]
]];
myHNReImCx[n_Integer] :=                            (* real part, compiled *)
myHNReImCx[n] =
Compile @@ Hold[{x, y}, Evaluate@First@myHNReIm[n, x, y]];
myHNReImCy[n_Integer] :=                            (* imaginary part, compiled *)
myHNReImCy[n] =
Compile @@ Hold[{x, y}, Evaluate@Last@myHNReIm[n, x, y]];
myHNReImC[n_Integer, x_?NumericQ,  y_?NumericQ] :=  (* interface to compiled parts *)
{myHNReImCx[20][x, y], myHNReImCy[20][x, y]};


Use the complex-valued function

Using myHN, which returns complex values is faster than using the individual components (such as in myHNReIm). This is true for HarmonicNumber, too. First of all, FindRoot can deal with such functions. For either myHN or HarmonicNumber, the Jacobian can be computed symbolically, which helps FindRoot to be efficient. And the way the real and imaginary parts are used in the call to FindRoots2D, the harmonic number has to be computed twice for each input.

Use Check to catch troublesome seeds

One can use Check to catch error/warning messages and re-execute FindRoot with adjusted settings. The following catches singular Jacobians, line-step tolerance warnings, and nonconvergence after the maximum number of interations and restarts FindRoot with a higher working precision. If messages are generated on the second try, an empty solution {} is returned on the assumption that a zero is not going to be found near the seed (starting point) x0. One could instead store them in some other variable and inspect them later. They might correspond to an actual root.

Check[{sol = FindRoot[f, {x, x0}]
Check[FindRoot[f, {x, x /. sol}, WorkingPrecision -> 100, AccuracyGoal -> 8],
{}],
{FindRoot::jsing, FindRoot::lstol, FindRoot::cvmit}]


Whether to restart FindRoot where it stopped x /. sol or from the initial starting point x0 is open to debate. If the solution sol is closer to the true zero, it will save some time. If round-off error has moved the solution far away from the true zero, the zero may be missed. In the OP's case, it seems best to start at the original seed (see below).

Use MeshFunctions to look for seed points

The following an arguable improvement. One can use ContourPlot in the form

ContourPlot[f1, {x, 0, 1}, {y, 0, 1}, MeshFunctions -> {f2}, Mesh -> {{0}}]


to find the seed points, which I've used before (see also 1, 2). Here f2 needs to be a pure function; it also ought to be continuous for good behavior. ContourPlot takes time and adding mesh functions increases the time it takes. It actually will search for the point on the contour line where f2 == 0.

Use the real and imaginary parts both ways in ContourPlot

Using one ContourPlot of the real part and searching every contour line for a change in sign of the imaginary part, as in FindRoots2D, is generally a good way to find zeros. One could also switch the roles of the real and imaginary parts. It turns out, more zeros are discovered if you do both. Since ContourPlot takes time, both plots could be computed in parallel, which won't take any more real time, if you have at least two cores.

Avoid Union

Unless you need the zeros sorted, Union tends to be slow. Here is a way to delete repeated roots that has two advantages, it is a bit faster and you can pick the root that gives the smallest value of HarmonicNumber. While HarmonicNumber is rather oscillatory, this probably gives the root closest to the actual zero. The function findClusters gathers, for each data point, the data points within a certain radius and then deletes the duplicate clusters:

findclusters[data_, r_: 10^-12] := With[{nf = Nearest[data]},
DeleteDuplicates[Sort[nf[#, {All, (* use Length[data] instead of All for pre-V10 *)
r}]] & /@ data]
];


One can then select the member of each cluster that produces the smallest function value like this:

reducedzeros = #[[First@Ordering[Abs @ myHN[20, #]]]] & /@ findclusters[zeros, 10^-6]


Overall this was about 6 times faster than Union on 27000 zeros.

Stopping runaway roots

If we use FindRoot like this

FindRoot[f, {z, seed, -30 + 0. I, 1 + 6000 I}]


we can catch the message FindRoot::reged with Check. This message indicates that the search for the root has left the region. If the region is slightly larger than {-30 + 0. I, 1 + 6000 I}, then perhaps such a seed can be safely discarded. It saves a lot of time, since FindRoot will pursue a root well beyond the region of interest.

Parallelize

There are a few things that can be parallelized, and since multi-core machines are now common, it makes sense to do so. This is particularly so for the FindRoot calls on the thousands of seeds.

### Some code

I'll separate FindRoots2D into two parts, finding the initial seeds and finding the zeros.

Clear[findSeeds];
findSeeds[f_, {x_, a_, b_}, {y_, c_, d_}, opts : OptionsPattern[ContourPlot]] :=
Module[{plotf, mf1, mf2},
plotf = Function @@ Hold[{x, y}, Evaluate@f];
If[OptionValue[MeshFunctions] === {},
{mf1, mf2} = {
With[{mf =
Compile[{x0, y0}, Evaluate[Re[f] /. {x -> x0, y -> y0}]]},
mf[##] &],
With[{mf =
Compile[{x0, y0}, Evaluate[Im[f] /. {x -> x0, y -> y0}]]},
mf[##] &]
},
{mf1, mf2} = OptionValue[MeshFunctions]
];
Join @@ Parallelize[
{Cases[Normal@#, Point[p_] :> p, Infinity] &[
ContourPlot[
Re @ f == 0,
{x, a - (b - a)/97, b + (b - a)/103}, {y, c - (d - c)/98, d + (d - c)/102},
ContourStyle -> None, MeshFunctions -> {mf2}, Mesh -> {{0}}, opts]
],
Cases[Normal@#, Point[p_] :> p, Infinity] &[
ContourPlot[
Im @ f == 0,
{x, a - (b - a)/97, b + (b - a)/103}, {y, c - (d - c)/98, d + (d - c)/102},
ContourStyle -> None, MeshFunctions -> {mf1}, Mesh -> {{0}}, opts]
]}
]
]


The tolerance to reject a seed as a nonroot is set at 10^-6. This could be changed or made an option. It appears in two places.

Clear[myFRC];
myFRC[f_, {z_, a_, b_}, seeds_, opts___] := Module[{z1, z2},
z1 = Re[a - (b - a)/97] + Im[a - (b - a)/98] I;
z2 = Re[b + (b - a)/103] + Im[b + (b - a)/102] I;
z /. Flatten[
ParallelMap[
Module[{sol},
Quiet[       (* did not suppress FindRoot::jsing *)
Check[      (* for FindRoot::reged : no solution *)
Check[{    (* for FindRoot::jsing, FindRoot::lstol, FindRoot::cvmit : restart *)
sol = FindRoot[f, {z, #1, z1, z2},
Evaluate[FilterRules[{opts}, Options[FindRoot]]]]},
Check[    (* for any message : no solution *)
If[Abs[f /. sol] < 10^-6,   (* tolerance 10^-6 for nonzero *)
{sol},
If[Abs[f /. #] < 10^-6, {#}, {}] &@
FindRoot[f, {z, SetPrecision[#1, 100], z1, z2},
WorkingPrecision -> 100, AccuracyGoal -> 8,
Evaluate[FilterRules[{opts}, Options[FindRoot]]]]],
{}],
{FindRoot::jsing, FindRoot::lstol, FindRoot::cvmit}],
{}, {FindRoot::reged}],
{FindRoot::reged, FindRoot::lstol, FindRoot::cvmit}
]
] &,
seeds
],
1]
]


Example

The steps in this example could replace the body of the code for FindRoots2D, if desired.

I used Pick with UnitStep instead of Select. In this case it is about ten times faster. Packing the array (DeveloperToPackedArray@N@... makes only a little difference, on my machine); there is also a small time savings from reducing the arbitrary precision zeros to MachinePrecision, whether the list is packed or not.

n = 20;
seeds = findSeeds[myHN[n, x + I y], {x, -30, 1}, {y, 0, 6000},
PlotPoints -> 25,
MeshFunctions -> {myHNReImCx[n][#1, #2] &,
myHNReImCy[n][#1, #2] &}]; // AbsoluteTiming
allfoundzeros =
Quiet @ myFRC[myHN[n, z], {z, -30, 1 + 6000 I},
seeds.{1, I}]; // AbsoluteTiming
uniqueszeros = #[[First @ Ordering[Abs @ myHN[20, #]]]] & /@
findclusters[allfoundzeros, 10^-6]; // AbsoluteTiming
zeros = Pick[
uniqueszeros,
With[{zz = DeveloperToPackedArray @ N @ uniqueszeros},
UnitStep[(#1 - (-30)) (1 - #1) #2 (6000 - #2)] &[Re[zz], Im[zz], 1]],
1]; // AbsoluteTiming
Length[zeros]
zeroplot = Graphics[{PointSize[Tiny], Red, Point @ Transpose@{Re[zeros], Im[zeros]}},
PlotRange -> {{-30, 1}, {0, 6000}}, AspectRatio -> 1, Frame -> True]

(*
{8.560900, Null}
{32.787730, Null}
{1.368914, Null}
{0.011737, Null}
{0.000582, Null}
2429
*)


Verification:

myHN[n, zeros] // Abs // Max
(* 5.11995*10^-7 *)


### Caveat

Being sure you have found all the zeros seems problematic. Changing the PlotPoints option gave different numbers of seeds and zeros. There seem to be about ten times as many seeds generated as zeros, which strikes me as a very poor ratio. Increasing the WorkingPrecision on ContourPlot did not change the result.

-
Thank you so much for your time on my problem. I just get back to school, so I did not know your update. –  user16023 Jun 30 '14 at 20:02
@user16023 You're quite welcome. No worries about the delay. –  Michael E2 Jun 30 '14 at 23:39
I could not have a chance to run you code until today. Unfortunately, I can not get the plot as you did. I check the code carefully and confuse that where the "foundzeros" in the example code came from. I guess it would be "allfoundzeros". But the code still does not excute. Could you give me some advices or post the entire code that produced your plot? Thank you very much. Have a nice day. –  user16023 Jul 10 '14 at 1:18
@user16023 See if it works now. Aside from changing the variable names, I must have had a findSeeds defined for complex z in Mma as well as the one posted for x + I y. Hopefully the update makes everything consistent. Worked when I copy-pasted it just now. –  Michael E2 Jul 10 '14 at 2:36
Yes. It works now. My slow computer takes roughly 30mins to run the code. But I am happy to see the plot showing up. And it has only one warning, "Nearest::near3: {All,1/1000000} is neither a non-negative integer nor a list of the form {n, r} with n a positive integer and r a non-negative real number. >> , and General::stop: Further output of Nearest::near3 will be suppressed during this calculation. >> ". I really appreciate if you could give me a hint to fix this warning as well. Thanks again and have a good one. –  user16023 Jul 10 '14 at 21:28