Derivative of Network Loglikelihood for brlen optimization

[lutteropp]

(from #43)

Use Newton-Raphson instead of Brent for optimizing branches (this will give us some niiiiiice additional speedup). This requires some Maths now, as I need to understand how to compute the derivative of network loglikelihood. Given the derivatives of the displayed tree loglikelihoods, this should be easy. I then need to understand how pll_modules computes likelihood derivatives for trees (something something sumtable), and how this can be generalized to networks.

Luckily, given the derivatives of displayed tree loglikelihoods, the derivative of the network loglikelihood is easy to compute (little mistake in the image, but not relevant here - the displayed tree loglikelihoods of course depend on the partition, also on brlens, model parameters, etc):

(EDIT: THIS IS WRONG, SEE LATER COMMENT)

Thus the main challenge here is to understand these sumtables used for tree loglh derivatives in pll_modules, being especially careful about how this works with virtual rerooting. And then I need to do this for the displayed trees displayed by the network and voilà, then we can switch from Brent brlen-opt to Newton-Raphson brlen-opt! :-)

[lutteropp]

One tricky thing: There can be displayed trees for which the current branch we want to optimize is inactive.
What happens to their sumtables? Since here we care about derivatives, we can safely skip those trees, not compute sumtables for them, and the derivative and second derivative of these tree loglikelihoods is zero.

There is even some more CLV saving potential in this (to be explored if later enough runtime gets spent in CLV updates again): If we can safely skip a tree, then we don't need to update its CLVs after virtual rerooting.

[lutteropp]

I made a mistake in the partition loglikelihood formulas above. New image coming soon.

[lutteropp]

Do we really need the correct derivative of phylogenetic network likelihood? It is apparently pretty ugly: https://en.wikipedia.org/wiki/Product_rule#Product_of_more_than_two_factors

[lutteropp]

Now I did the Math correctly, and it does not look good for Newton-Raphson on phylogenetic networks :-( We can still do it as it is known to converge quicker, but it does not help in getting rid of pll_compute_edge_loglikelihood.

[lutteropp]

Or maybe we can still do something here? I need to dig into how the Math is done for derivative of tree likelihood. After all, there is a product involved there, too (the product over per-site likelihoods ... plus, also a product over partitions).

[lutteropp]

Whichever trick is used for trees there should also work for the networks. Because with the approach I used here, I get into the same problem if there is only a single displayed tree. There has to be some further trick involved!

[lutteropp]

Okay... now I understand the text in Alexeys thesis (https://cme.h-its.org/exelixis/pubs/dissAlexey.pdf, page 26).
What this sumtable gives us is fast computation of:

• L(b)
• L'(b)
• L''(b)

My problem was that I thought it only gives L'(b) and L''(b)

[lutteropp]

My confusion originated from this: For trees, we end up getting and using these:

• (loglL(T))' -> first derivative of loglikelihood of tree T
• (logL(T))'' -> second derivative of loglikelihood of tree T

But what the trick with the sumtable actually gives us are fast computations of (with respect to branch b):

• L(b) -> persite likelihood of tree T
• L'(b) -> first derivative of persite likelihood of tree T
• L''(b) -> second derivative of persite likelihood of tree T

Instead of directly reusing the final result ((loglL(T))' and (loglL(T))'') from the trees case, I need to plug these faster-computed tree persite likelihoods and their derivatives into the network loglikelihood and loglikelihood derivative formulas.

[lutteropp]

I got it all together now.
The main difference (compared to the tree case) is that in order to compute network loglikelihood derivatives, we also need to use the sumtables to compute the displayed tree loglikelihoods. Only the tree loglikelihood derivatives are not enough for us.

[lutteropp]

Wait... do I really need the fully computed displayed tree loglikelihoods here?

For LikelihoodModel.BEST, we can do without explicitly computing them. Because there we essentially have one tree per partition...
For LikelihoodModel.AVERAGE, I currently see no way around actually computing the displayed tree loglikelihoods... But I am still not convinced that they are really needed.

--> Looks like there is no way around actually inserting the displayed tree persite loglikelihoods (and their derivatives) into the formulas for the networks case. Thus, another drawing will follow soon here. :-)

[lutteropp]

Turns out that for networks, we can't go down to a per-site computation basis if we have more than one displayed tree.

• I am convinced about this in the LikelihoodType.AVERAGE case.
• For the LikelihoodType.BEST case it depends on whether max{f', g'} / max{f, g} = max{f'/f, g'/g} always holds or not.

[lutteropp]

Okay, found a counterexample (for x > 0):

• f(x) = 2x + 3
• g(x) = 2x + 5
• f'(x) = 2
• g'(x) = 2
• max{f'(x), g'(x)} = 2
• max {f(x), g(x)} = 2x+5
• max{f'(x), g'(x)} / max{f(x), g(x)} = 2/(2x+5)
• max {f'/f, g'/g} = max{2/(2x+3), 2/(2x+5)} = 2/(2x+3)

----> max{f', g'} / max{f, g} != max{f'/f, g'/g}

---> Thus, no matter which likelihood model (out of AVERAGE and BEST) we are using, as soon as we have more than one displayed tree, we cannot profit from going down to the persite-level.

[lutteropp]

This makes it clear that for networks, we need to compute all three of logL(T|P), (logL(T|P))', (logL(T|P))'' with help of the sumtable.

The only question remains whether it makes sense to go down to a per-site level in case the network has only one displayed tree. Special treatment for the single-tree case has led to numerical quirks when computing network loglikelihood before (there, while mathematically the same, the different operations (like, adding a log(exp()) around it or not) led to a slight numerical difference that confused BIC).

[lutteropp]

Computation of second network loglikelihood derivative caused numerical problems (underflow/overflow). Thus, I had to play a bit with the formula:

[lutteropp]

Damn the Newton Raphson in pll-modules is awful, also the second network loglikelihood derivative is such a high number that it does not fit into a double. Luckily, what Newton-Raphson actually needs is just the quotient of first and second derivative, as well as the sign of the second derivative...

I'll just switch to implementing vanilla Newton-Raphson myself, using https://www.cup.uni-muenchen.de/ch/compchem/geom/nr.html and Alexis lecture slides ...