GHOST : WARNING: Numerical underflow for lh-derivative-mixlen

[roisin-long]

Hi,

When I try to run GHOST (JTT+H4) for some alignments with the -safe flag, massive log files with "WARNING: Numerical underflow for lh-derivative-mixlen" are made. For these alignments, the runs don't finish (at least within the timeframe alignments without this issue do) and I have to stop them because the log files get too big. I have seen this is a reported issue in previous IQ-TREE versions. I get this issue with multiple versions - 2.0.7, 3.0.1 and 3.1.2. Is it possible for this issue to be fixed/ is there a suggested fix on the user end for this?

Thank you!

[StefanFlaumberg]

Hi @roisin-long,

I've recently been working on making GHOST rates compatible with PMSF, which mainly includes tinkering with derivative-computing functions, so it's interesting for me to look into your issue as well.
I'm not sure what's causing the problem, but I'd suggest you append -optlen BFGS to your cmd to use the BFGS algorithm instead of the default expectation maximization (EM) algorithm for heterotachous branch length optimization.
Please write back whether this trick helps. If it does, then the bug is likely in the EM algorithm (the main suspect); your report will help to fix this.

Best,
Stefan

[roisin-long]

Hi @StefanFlaumberg ,

Thank you for your suggestion - I tested appending -optlen BFGS for one of the alignments that previously had this issue and this solves the problem (both when running with and without the -safe flag).

All the best,
Róisín

[bqminh]

That's interesting to know... thanks for figuring that out. We can extend the warning message, telling users to use that option... Or probably better is to switch the optimisation to use BFGS, instead of the EM algorithm, in case that warning occurs...

[StefanFlaumberg]

Hi @roisin-long,
Thanks for the testing! It's good to know that things work well at least with BFGS.
I think for now you can use this trick as a workaround for the problem. But I suggest we keep this issue open until a proper solution to the problem is found.
To this end, could you provide an alignment (ideally, not too large) to reproduce the problem? (I remember having observed the problem myself, but cannot find the relevant data now.)

Hi @bqminh,
Since BFGS worked out fine, I suspect that the problem isn't just due to peculiarities of the data, but a bug in the EM algorithm. Then, I think it should be fixed, not just masked with a workaround suggestion or an automatic switching to BFGS.
Now I'm going through the GHOST code for unrelated reasons (PMSF), so I hope I'll find out soon what causes the problem. Anyway, I'll write back here whether I succeed or fail.

Best regards,
Stefan

[roisin-long]

Hi @StefanFlaumberg ,

I've sent that alignment on now. Thanks for looking into this!

All the best,
Róisín

[StefanFlaumberg]

@roisin-long, great thanks for the alignment! It has helped a lot!

@bqminh, the bug is in the way the categories of the theta_all buffer are scaled in the safe_numeric mode.
In brief:
The text of the warning implies that the problem can arise only in the safe_numeric mode, activated either with the -safe option or automatically due to having too many (>=2000) leaves; otherwise, the warning would also say "Run again with the safe likelihood kernel via -safe option" (see computeLikelihoodDervMixlenSIMD()). In the safe_numeric mode, only the rate-mixture categories with scale <= min_scale+1 are actually used, others being set to zero as very small numbers (see computeLikelihoodBufferSIMD()). The EM algorithm optimizes branch lengths of only one rate (heterotachy) category at a time. Now, if all rate-mixture categories of the current rate category happen to have scale > min_scale+1, hence set to zero, the lh of the current rate category will be zero as well, causing infinite derivative df = log(lh)' = lh'/lh.
The fix is straightforward -- to define min_scale as the minimal scale across the rate-mixture categories belonging only to the current heterotachy category if using the EM algorithm. Also, the theta_all buffer should be recomputed when switching to the next heterotachy category during the EM optimization and right after it, since the buffered values are now rate-category-specific.

I shall add the fix among other changes in an upcoming pull request shortly after my GHOST+PMSF pull request is merged.

Best regards,
Stefan

[bqminh]

I see, makes sense. If you make a PR I'll review later.