Overflow of entropy (again)

[jeanluct]

Miriam Ritchie (Dundee) is having some trouble computing entropy for braids with a large entropy. The issue is with overflow of intermediate steps. This has been a tricky issue to get right in an elegant way (see #39, #42, #61, #69).

As a first step here's a code that shows when overflow occurs, and two simple but inelegant methods to get around it:

% Test program to check when the internal algorithm for entropy overflows.

import braidlab.*

% A simple pA braid and its entropy.
b0 = braid([1 -2]);
entr0 = entropy(b0);

% Estimate max number of repetitions before we overflow.
Nrepmax = ceil(log(realmax)/entr0);

for divNrepmax = [10 5 3 2]
  % Number of repetitions of basic braid: inch our way towards Nrepmax.
  rep = ceil(Nrepmax/divNrepmax);
  fprintf('b0^%d\thas entropy %.3e',rep,entropy(b0^rep));
  fprintf(' (exact=%.3e)\n',rep*entr0)
end

% Direct method: compute generator-by-generator, 
% Do the last case above.
% (This is much slower!  Use larger "chunks" to make it faster.)
b = b0^rep;
l = loop(b.n); l = loop(l.coords/minlength(l));  % normalized initial loop
loglen = 0;
for i = 1:length(b)
  l1 = braid(b.word(i),b.n)*l;
  len = minlength(l1);
  loglen = loglen + log(len); % keep track of log-growth
  l = loop(l1.coords/len);    % keep loop coordinates normalized to avoid
                              % overflow.
end
fprintf('\nentropy from per-generator computation=%.3e\n',loglen)

% Yet another way: use complexity, but with arbitrary precision.
l = loop(b.n,@vpi);
l1 = b*l;
% The minlength function doesn't seem to work well with VPI, so use L^1
% norm of entries instead.
fprintf('\nentropy using Variable Precision Integers (VPI)=%.3e\n', ...
        log(sum(abs(double(l1.coords)))) - log(sum(abs(double(l.coords)))))

To be clear, here the overflow comes from the unrealistic repetition of the same braid, but the same could happen with a very long random braid.

This code was added to the develop branch as devel/tests/test_entropyOverflow.m. Maybe eventually create a feature branch.

[jeanluct]

John Lynch has hit the same issue, I think. He wrote a code called bigTEPO.m which generates large braids:

>> type bigTEPO.m

function [TEPO,b] = bigTEPO(n)
% Computes the TEPO (with respect to the new generating set) of a family of high-TEPO braids.

import braidlab.*

% First operation: swap punctures 1 and n, 2 and n-1, etc.
b = braid([],n);
for i = 1:floor(n/2)
    b = b*bigGen(i,n-i+1,n,'new');
end

% Second operation: two more inverse twists to make it pseudo-Anosov.
twist1 = bigGen(1,floor(n/2),n,'new');
twist2 = bigGen(floor(n/2)+1,n,n,'new');
b = b*inv(twist1)*inv(twist2);

TEPO = entropy(b)/2;

Then

>> bigTEPO(1894)

ans =

    4.4657

and

>> bigTEPO(1895)

ans =

  3.4664e-310

Oops. Overflow it seems.

Turn on debugging:

>> global BRAIDLAB_debuglvl
>> BRAIDLAB_debuglvl=2;

>> bigTEPO(1895)

TOL = 1.0e-06    MAXIT = 2148944991      NCONV = 3   LENGTH = l2norm

Converged 3 time(s) in a row after 1 iterations

ans =

  3.4664e-310

Now weird things happen when we disable the MEX file:

>> global BRAIDLAB_braid_nomex
>> BRAIDLAB_braid_nomex=true;

>> bigTEPO(1895)

TOL = 1.0e-06    MAXIT = 2148944991      NCONV = 3   LENGTH = l2norm

Number of threads auto-set to 4 using java.lang.Runtime.
Using MEX loopsigma with Matlab data structures.
loopsigma_helper: multiplication running UNTHREADED.
  iteration 1  entr=1.4462026612e+01  diff=1.5462e+01
Previously set ComputedThreads=4. Run 'clear getAvailableThreadNumber'  to recompute
Using MEX loopsigma with Matlab data structures.
loopsigma_helper: multiplication running UNTHREADED.
  iteration 2  entr=8.9095155465e+00  diff=-5.5525e+00
Previously set ComputedThreads=4. Run 'clear getAvailableThreadNumber'  to recompute
Using MEX loopsigma with Matlab data structures.
loopsigma_helper: multiplication running UNTHREADED.
  iteration 3  entr=8.9317451772e+00  diff=2.2230e-02
Previously set ComputedThreads=4. Run 'clear getAvailableThreadNumber'  to recompute
Using MEX loopsigma with Matlab data structures.
loopsigma_helper: multiplication running UNTHREADED.
  iteration 4  entr=8.9317501970e+00  diff=5.0198e-06
Previously set ComputedThreads=4. Run 'clear getAvailableThreadNumber'  to recompute
Using MEX loopsigma with Matlab data structures.
loopsigma_helper: multiplication running UNTHREADED.
  iteration 5  entr=8.9317501980e+00  diff=9.9344e-10
Previously set ComputedThreads=4. Run 'clear getAvailableThreadNumber'  to recompute
Using MEX loopsigma with Matlab data structures.
loopsigma_helper: multiplication running UNTHREADED.
  iteration 6  entr=8.9317501980e+00  diff=6.1462e-13
Previously set ComputedThreads=4. Run 'clear getAvailableThreadNumber'  to recompute
Using MEX loopsigma with Matlab data structures.
loopsigma_helper: multiplication running UNTHREADED.
  iteration 7  entr=8.9317501980e+00  diff=-8.3489e-14
Converged 3 time(s) in a row after 7 iterations

ans =

    4.4659

It seems the C++ code only does one iteration, whereas the pure Matlab code converges fine after 7 iterations. Not sure what is going on, but now upgrading (downgrading?) this to an actual bug.

[jeanluct]

Ok, some progress. I think this is actually two bugs! The one that John uncovered is easy to fix: maxit overflowed a 32-bit int, so I set its max at a 32-bit int. It prints a warning, but really it should realistically never be able to do that many iterations. Therefore it is perfectly fine to issue

>> warning('off','BRAIDLAB:braid:entropy:maxitint32')

to suppress this warning.

BUT, that still leaves the other original problem found by Miriam. Her problem is that the braid actually does have larger entropy, so the internal calculation is overflowing. What would be needed is a renormalization of the braid as we go, but this is not always a good idea. Maybe have a flag for renormalizing?

[jeanluct]

Ok! I think this is fixed. There's a more detailed explanation in the code, but the idea is to compute growth 300 generators at a time. This means that even for the max TEPG braid the growth will not overflow.

There's a new unit test in entropyTest.m that tests a massive braid, with entropy 400 or so. (That's right, entropy 400, not dilatation! So the dilatation is e^400.)

Closing! This issue has been a real pain. I suspect we've not seen the last of it.