C functions for precomputing normalization constant do not work for zero entries

[gerhardkurz]

I noticed something odd about the C functions for precomputing the normalization constant. When I call the function

bingham_F_3d(double z1, double z2, double z3)

with parameters z1 <= z2 <= z3 < 0, I get the correct normalization constant. However, if z3==0 or z2==0 and z3==0, I get an incorrect result.

For example, z1=-3, z2=-2, z3=-1 yields 5.4011, which is correct. However, z1=-3, z2=-2, z3=0 yields 3.5491, but the correct result would be 7.3201. Also, z1=-3, z2=0, z3=0 yields 2.9292, but the correct result would be 11.5765.

I think this is somehow related to the fact that compute_1F1_3d reduces the computation to the two-dimensional case if z3==0:

if (fabs(z3) < EPSILON)  // z3 = 0
    return sqrt(M_PI)*compute_1F1_2d(dim, z1, z2, iter);

Although this might be possible in principle, the result does not seem to be correct for some reason. In the case of z2==0, the calculation is further reduced to the one-dimensional case.

[SebastianRiedel]

Hi Gerhard, I can’t reproduce the wrong constants using the higher-level functions of Jared (not using e.g. bingham_F_3d directly), see outputs below. Jared uses pre-computed constants from a look-up table, so when initializing a bingham object, I'm pretty sure the functions you tried are not called, but rather the normalization constant is looked up. I don’t know exactly how the look-up table is generated or if there is any magic happening which prevents the wrong results you saw propagating into the look-up table.

Anyway, here is the results I got using the bingham_sample executable (I just added a print_bingham(&B); statement to it. They are pretty close to your correct values:

~/software/bingham/c [master]$ ./bingham_sample -3 -2 -1 5
Sampling from bingham:
B->F = 5.404276
B->Z = [ -3.000000 -2.000000 -1.000000 ]
B->V[0] = [ 0.000000 1.000000 0.000000 0.000000 ]
B->V[1] = [ 0.000000 0.000000 1.000000 0.000000 ]
B->V[2] = [ 0.000000 0.000000 0.000000 1.000000 ]
********* seed = 1437923666
0.882952 0.376986 0.122291 -0.251640
0.882952 0.376986 0.122291 -0.251640
0.882952 0.376986 0.122291 -0.251640
0.801891 0.252364 0.327529 0.431287
-0.007494 -0.484809 -0.332580 -0.808885

~/software/bingham/c [master]$ ./bingham_sample -3 -2 0 5
Sampling from bingham:
B->F = 7.323641
B->Z = [ -3.000000 -2.000000 0.000000 ]
B->V[0] = [ 0.000000 1.000000 0.000000 0.000000 ]
B->V[1] = [ 0.000000 0.000000 1.000000 0.000000 ]
B->V[2] = [ 0.000000 0.000000 0.000000 1.000000 ]
********* seed = 1437923873
0.320445 0.667523 -0.435943 -0.511549
-0.275925 0.459170 -0.036546 0.843619
0.248626 -0.529573 -0.671189 0.455238
0.201899 0.529786 -0.369366 -0.736297
0.260022 0.340124 -0.779548 0.457174

~/software/bingham/c [master]$ ./bingham_sample -3 0 0 5
Sampling from bingham:
B->F = 11.579555
B->Z = [ -3.000000 0.000000 0.000000 ]
B->V[0] = [ 0.000000 1.000000 0.000000 0.000000 ]
B->V[1] = [ 0.000000 0.000000 1.000000 0.000000 ]
B->V[2] = [ 0.000000 0.000000 0.000000 1.000000 ]
********* seed = 1437923912
0.885661 -0.219993 0.258030 -0.317220
-0.829918 0.183327 -0.463881 -0.249882
0.257433 0.895597 -0.030125 0.361563
0.144841 -0.368585 0.913107 0.096966
-0.189069 0.636853 0.308998 0.680581

It would be helpful to know how the look-up table is generated. I wanted to do that at some point anyway, but at the moment I need to focus on other projects at work. Feel free to investigate on your own.

[gerhardkurz]

Thank you very much. You are correct, the error does not occur when using the functions that do a table look-up.

As far as I know, the tables are generated from MATLAB using the scripts: bingham_compute_F_cache.m export_bingham_constants.m

In principle it should be possible to generate the tables from C as well, but I suspect that this might fail because of the bug described above.