Learning

This chapter describes how the global parameters space is explored to maximize aphid likelihood.

Heuristic starting point

If not every initial parameter has been specified by user, aphid uses the following procedure to automatically choose initial parameters.

(this reuses notations from the likelihood formula)

First, calculate the mean branches counts among the selected trees: \[ % Factorize common mathjax commands. \newcommand{gt}{gt} \newcommand{pab}{{p_{ab}}} \newcommand{pac}{{p_{ac}}} \newcommand{pbc}{{p_{bc}}} \newcommand{po}{{p_{o}}} \newcommand{gt}[1]{{gt_{#1}}} \newcommand{sl}{{sl}} % \newcommand{So}{{S_\varnothing}} \newcommand{SI}[1]{{S_{ILS(#1)}}} \newcommand{SG}[2]{{S_{GF(#1,\ #2)}}} % \newcommand{pgf}{{p_{GF}}} \newcommand{pils}{{p_{ILS}}} % \newcommand{ABC}{{\scriptsize \text{ABC}}} \newcommand{ACB}{{\scriptsize \text{ACB}}} \newcommand{BCA}{{\scriptsize \text{BCA}}} \newcommand{CAB}{{\scriptsize \text{CAB}}} \newcommand{CBA}{{\scriptsize \text{CBA}}} % % Expected branches lengths. \newcommand{exb}[3]{{\bar{#1}({#2},\ {#3})}} \newcommand{exu}[1]{{\bar{u}({#1})}} \newcommand{exv}[1]{{\bar{v}({#1})}} \begin{align} m_a &= \frac{1}{G}\sum_g{n_g(a_g)}\qquad\qquad m_b = \frac{1}{G}\sum_g{n_g(b_g)}\\ m_c &= \frac{1}{G}\sum_g{n_g(c_g)}\qquad\qquad m_d = \frac{1}{G}\sum_g{n_g(d_g)} \end{align} \]

Then, calculate the frequency of discordant topologies among the selected trees: \[ \begin{align} f_{ac} &= \frac{1}{G}\sum_g{𝟙(T_g = \ACB)}\\ f_{bc} &= \frac{1}{G}\sum_g{𝟙(T_g = \BCA)} \end{align} \]

Then use these as default starting points:

\[ \begin{align} \tau_1 &= \frac{m_a + m_b - \theta}{2} \\ \tau_2 &= \frac{m_c + m_d + \tau_1 - \theta}{2} \\ \pils &= e^{2\frac{\tau_1 - \tau_2}{\theta}} \\ \pac &= f_{ac} - \pils / 3 \\ \pbc &= f_{bc} - \pils / 3 \\ \pab &= \frac{\pac - \pbc}{2} \\ \po &= \frac{1}{|gt|} \end{align} \]

Note that the above formula may yield negative values for \(\tau_1\) and \(\tau_2\), in which case their starting value is set to 0 instead.

(see aphid's default values for theta = \(𝜃\))

Reparametrization

Exploring a constrained space is more difficult than exploring a space without constraints. For this reason, aphid's constrained parameters are rewritten as functions of unconstrained scores, ranging freely over \(ℝ\).

The "sigmoid" function and its inverse are useful in this respect: \[ \newcommand{rev}{𝜎^{-1}} % Factorize common mathjax commands. \newcommand{gt}{gt} \newcommand{pab}{{p_{ab}}} \newcommand{pac}{{p_{ac}}} \newcommand{pbc}{{p_{bc}}} \newcommand{po}{{p_{o}}} \newcommand{gt}[1]{{gt_{#1}}} \newcommand{sl}{{sl}} % \newcommand{So}{{S_\varnothing}} \newcommand{SI}[1]{{S_{ILS(#1)}}} \newcommand{SG}[2]{{S_{GF(#1,\ #2)}}} % \newcommand{pgf}{{p_{GF}}} \newcommand{pils}{{p_{ILS}}} % \newcommand{ABC}{{\scriptsize \text{ABC}}} \newcommand{ACB}{{\scriptsize \text{ACB}}} \newcommand{BCA}{{\scriptsize \text{BCA}}} \newcommand{CAB}{{\scriptsize \text{CAB}}} \newcommand{CBA}{{\scriptsize \text{CBA}}} % % Expected branches lengths. \newcommand{exb}[3]{{\bar{#1}({#2},\ {#3})}} \newcommand{exu}[1]{{\bar{u}({#1})}} \newcommand{exv}[1]{{\bar{v}({#1})}} \begin{align} 𝜎 &↦ \begin{cases} ℝ &→ (0, 1) \\ x &↦ \frac{1}{1 + e^{-x}} \end{cases} \\ \rev &↦ \begin{cases} (0, 1) &→ ℝ \\ y &↦ \ln{\left(\frac{y}{1 - y}\right)} \end{cases} \end{align} \]

Considering that the following scores live in ℝ with no constraint:

\[ \begin{align} s_𝜃 &= \ln(𝜃) \\ s_{𝜏_1} &= \ln(𝜏_1) \\ s_{𝛥𝜏} &= \ln(𝜏_2 - 𝜏_1) \\ s_{gf} &= \rev(\pab + \pac + \pbc) \\ s_{ac} &= \ln\left(\pac × \frac{e^1}{\pab}\right) \\ s_{bc} &= \ln\left(\pbc × \frac{e^1}{\pab}\right) \\ s_\po &= \rev(\po) \\ s_{t_i} &= \begin{cases} i = 1 &\implies \rev(\gt{1}) \\ i > 1 &\implies \rev\left(\frac{\gt{i}}{\gt{i-1}}\right) \end{cases} \end{align} \]

then here is the reparametrization formula:

\[ \begin{align} 𝜃 &= e^{s_ 𝜃} \\ 𝜏_1 &= e^{s_{𝜏_1}} \\ 𝜏_2 &= 𝜏_1 + e^{s_{𝛥𝜏}} \\ \pab &= \frac{e^1}{e^1 + e^{s_{ac}} + e^{s_{bc}}} × 𝜎(s_{gf}) \\ \pac &= \frac{e^{s_{ac}}}{e^1 + e^{s_{ac}} + e^{s_{bc}}} × 𝜎(s_{gf}) \\ \pbc &= \frac{e^{s_{bc}}}{e^1 + e^{s_{ac}} + e^{s_{bc}}} × 𝜎(s_{gf}) \\ \po &= 𝜎(s_\po) \\ \gt{i} &= \begin{cases} i = 1 &\implies 𝜎(s_{t_1}) \\ i > 1 &\implies \gt{i-1} × 𝜎(s_{t_i}) \end{cases} \end{align} \]

Instead of the constrained \(𝜃\), \(𝜏_1\), \(𝜏_2\) etc., the concrete optimisation targets for aphid are the unconstrained \(s_𝜃\), \(s_{𝜏_1}\), \(s_{𝛥𝜏}\) etc.

Maximizing likelihood

Once reparametrized, the likelihood is viewed as a function \(f\) of the \(s\) scores \(X\):

\[ f:\left\{\begin{align} ℝ^s & → ℝ \\ X & ↦ f(X) = f(s_𝜃, s_{𝜏_1}, …) = -\ln(ℒ) \end{align}\right. \]

The logarithm is used to improve numerical stability of calculations. Maximizing the likelihood boils down to finding a value of \(X\) in \(ℝ^s\) that minimizes \(f\).

Aphid uses a heuristic approach to minimization commonly referred to as Broyden–Fletcher–Goldfarb–Shanno algorithm (BFGS) as described in Nocedal & Wright (2006). The method uses successive first-order derivatives of \(f\) to construct a running estimate of its second-order derivative and exploit the local curvature of the likelihood surface to greedily converge towards local minima.

This section specifies aphid's BFGS implementation detail and configuration.

Derivatives

Given a sample \(X_k\), aphid uses the likelihood formula to calculate \(f(X_k)\). If required, it then uses automatic differentiation to calculate the exact local gradient \(∇f(X_k)\). Both steps are performed with tensors from the external C++ library pytorch, bound to aphid via the Rust crate tch-rs.

This dependency to automatic differentiation may be transfered to candle in the future.

Parameters forcing

Any parameter forced by user erases and replaces the value calculated from the scores during likelihood calculation, in a way that makes the likelihood gradient null with respect to these parameters, and prevents automatic differentiation to flow upstream to the optimized scores.

The consequence is that both likelihood formula and the "scores → parameters" formulae are unchanged when forcing, and that all scores variables will always appear in the output even when some parameters have been force-fixed.

BFGS Update

Given an initial starting point \(X_0\), the procedure starts by initializing a running estimate of the function's inverse Hessian matrix to the identity matrix:

\[ H_0 = 𝕀_{s, s} \]

On every step \(k\), a search direction \(P\) is chosen according to:

\[ P_k = -H_k\,. ∇f(X_k) \]

Samples of the following form are taken along this search direction, for various "step size" values \(𝛼 \in ℝ^+\):

\[ X_\bullet = X_k + 𝛼_\bullet\,P_k \]

until a sample verifying the strict Wolfe conditions is found. That is, a step size \(𝛼_k\) such that:

\[ X_{k+1} = X_k + 𝛼_k\,P_k \]

\[ \left\{\begin{align} f(X_{k+1}) \leqslant f(X_k) + c_1 \, 𝛼 \, ∇f(x_k)^{\intercal} P_k \qquad &\text{(Armijo's rule: sufficient decrease)} \\ ∇f(X_{k+1})^{T} \, P_k \geqslant c_2 \, ∇f(x_k)^{\intercal} P_k \qquad &\text{(curvature condition)} \end{align}\right. \]

for some fixed, positive values \(c_1\) and \(c_2\). These values can be chosen from the configuration table [search.bfgs.wolfe] with parameters c1 and c2, and they default to Nocedal & Wright (2006) recommendation:

\[ \left\{\begin{align} c_1 &= 10^{-4} \\ c_2 &= 0.1 \end{align}\right. \]

The exact "linear search" procedure to find \(𝛼_k\) is specified in the next section.

Once \(X_{k+1}\) has been chosen, the inverse Hessian approximation is updated with:

\[ \begin{alignat}{2} S_k &= X_{k+1} - X_k = 𝛼_k\,P_k &&(\in ℝ^s)\\ 𝛥_k &= ∇f(X_{k+1}) - ∇f(X_k) && (\in ℝ^s) \\ 𝜌_k &= (𝛥_k^{\intercal}\,.S_k)^{-1} &&(\in ℝ) \\ H_{k+1} &= (𝕀 - 𝜌_k\,S_k\,𝛥_k^{\intercal}) \, H_k \, (𝕀 - 𝜌_k\,&&𝛥_k\,S_k^{\intercal}) - 𝜌_k\,S_k\,S_k^{\intercal} \end{alignat} \]

Exploration continues until a termination criterion, as described in the upcoming termination section, is reached.

Linear Search

Given a sample \(X_k\) and a search direction \(P_k\), "linear search" is the process of exploring various successive candidate step sizes \(𝛼\) until we find one satisfying the Wolfe conditions.

Again, aphid's linear search strategy is taken from Nocedal & Wright (2006), section 3.5. Its implementation, specified here, has been refined to handle possible degenerated non-finite floating-point values i.e. inf or NaN.

The linear search problem considers the one-dimensional function \(𝜑\) mapping every candidate step size to the value of its corresponding candidate sample:

\[ 𝜑:\left\{\begin{align} ℝ^+ & → ℝ \\ 𝛼 & ↦ 𝜑(𝛼) = f(X_k + 𝛼\,P_k) \end{align}\right. \]

Aphid assumes that all floating-point values of \(𝜑(𝛼)\) are finite between \(𝛼 = 0\) and some "horizon" value \(𝛼 = h > 0\). It also assumes that the range \(𝛼 ∈ [0, h]\) does contain acceptable step sizes meeting the strict Wolfe conditions. The goal of the algorithm is to find such a step size while estimating \(h\) if necessary and bisecting within \([0, h]\) if found.

The algorithm is specified as a flowchart below, (configuration options \((d, i, c)\) highlighted) and described in the next few sections.

1. Reduction phase

Starting from \(𝛼 = 1\), and if necessary, decrease the step size by a constant factor \(d\) until 𝜑(𝛼) becomes finite.

This factor can be chosen from the configuration table [search.bfgs.wolfe] with parameter step_decrease.

If decreasing was necessary, record the lowest value of \(𝛼\) yielding non-finite \(𝜑(𝛼)\) as an upper-bound estimate \(\hat{h}\) of the horizon \(h\).

2. Bracketing phase

Increase \(𝛼\) by a constant factor \(i\) until either a satisfying candidate is found, or an interval within \([0, h]\) is proven to contain such a candidate.

This factor can be chosen from the configuration table [search.bfgs.wolfe] with parameter step_increase.

If increasing yields non-finite values for \(𝜑(a)\), then bisect within \([0, \hat{h}]\) while updating \(\hat{h}\) instead.

3. Zoom phase

Once an interval \([u, v] \subset [0, h]\) has been proven to contain a satisfying candidate, use cubic interpolation to quickly shrink \([u, v]\) and converge towards \(𝛼_k\).

In the case of aphid at least, I observed that cubic interpolation hindered convergence when cubic candidates values for \(𝛼\) where too close from either \(u\) or \(v\). For this reason, aphid falls back to naive bisection steps within \([u, v]\) whenever the cubic interpolation step would yield such a value.

The threshold for bisection is expressed as a fraction \(c\) of the interval length \(v - u\). It can be chosen from the configuration table [search.bfgs.wolfe] with parameter bisection_threshold.

Termination

The linear search terminates as a success when either of the following events occurs:

• A candidate satisfying Wolfe conditions is found.
• During the bracketting or the zooming phase, the binary search interval reduces below floating-point precision.
• During the zooming phase, a candidate \((𝜑(𝛼), 𝜑'(𝛼))\) cannot be distinguished from \((𝜑(u), 𝜑'(u))\) or \((𝜑(v), 𝜑'(v))\) within floating-point precision, and \(𝜑'(𝛼)\) is lower than the threshold defined in the configuration table [search.bfgs.wolfe] by parameter flat_gradient.

It terminates as a failure when either of the following events occurs:

• No finite floating-point \(𝜑(𝛼)\) value is found during the reduction phase.
• A non-finite value for \(𝜑'(𝛼)\) is obtained whereas \(𝜑(𝛼)\) was finite.
• A non-finite value for \(𝜑(𝛼)\) is found within \([u, v]\).

The whole BFGS optimisation procedure terminates as a success when either of the following events occurs:

• The best step resulting from linear search is numerically null: \(𝛼_k = 0\).
• The corresponding step norm is numerically null: \(||S_k|| = 0\).
• There is no difference between \(∇f(X_{k+1})\) and \(∇f(X_k)\) within floating-point precision, and the step norm is smaller than the threshold defined in the configuration table [search.bfgs] by parameter step_size_threshold.
• The mean slope over the last few \((f(X_{k - i + 1}))_{i \in \{ 1, .., n \}}\) values is lower than the threshold defined in the configuration table [search.bfgs.slope_tracking] by parameter threshold. The number of samples \(n\) is determined by parameter sample_size in the same table, and the slope is estimated every grain steps.
• The maximum number of iterations is reached, defined in table [search.bfgs.max_iter].

And as a failure when:

• Either initial \(f(X_0)\) or \(∇f(X_0)\) yields non-finite floating-point results.

Exploration

BFGS is greedily attracted towards the closest local optimum. To ensure exploration of the parameters space, aphid issues one parallel BFGS search per initial condition defined within the [init] table. The result yielding the best likelihood among all searches is the one kept for output, even if it does not correspond to a terminal search step.