1.2.1. Computing the direction \(d_k\)

A desirable property for \(d_k\) is that \(d_k\) ensures a descent \(f(x_{k+1}) < f(x_{k})\). Newton methods are such that \(d_k\) minimizes a local quadratic approximation of \(f\) based on a Taylor expansion, that is \(q_f(d) = f(x_k) + g(x_k)^Td +\frac{1}{2} d^T H(x_k) d\) where \(g\) denotes the gradient and \(H\) denotes the Hessian.

The consists in using the exact solution of local minimization problem \(d_k = - H(x_k)^{-1} g(x_k)\).
In practice, other methods are preferred (at least to ensure positive definiteness). The method approximates the Hessian by a matrix \(H_k\) as a function of \(H_{k-1}\), \(x_k\), \(f(x_k)\) and then \(d_k\) solves the system \(H_k d = - g(x_k)\). Some implementation may also directly approximate the inverse of the Hessian \(W_k\) in order to compute \(d_k = -W_k g(x_k)\). Using the Sherman-Morrison-Woodbury formula, we can switch between \(W_k\) and \(H_k\).

To determine \(W_k\), first it must verify the secant equation \(H_k y_k =s_k\) or \(y_k=W_k s_k\) where \(y_k = g_{k+1}-g_k\) and \(s_k=x_{k+1}-x_k\). To define the \(n(n-1)\) terms, we generally impose a symmetry and a minimum distance conditions. We say we have a rank 2 update if \(H_k = H_{k-1} + a u u^T + b v v^T\) and a rank 1 update if $H_k = H_{k-1} + a u u^T $. Rank \(n\) update is justified by the spectral decomposition theorem.

There are two rank-2 updates which are symmetric and preserve positive definiteness

• DFP minimizes \(\min || H - H_k ||_F\) such that \(H=H^T\): \[ H_{k+1} = \left (I-\frac {y_k s_k^T} {y_k^T s_k} \right ) H_k \left (I-\frac {s_k y_k^T} {y_k^T s_k} \right )+\frac{y_k y_k^T} {y_k^T s_k} \Leftrightarrow W_{k+1} = W_k + \frac{s_k s_k^T}{y_k^{T} s_k} - \frac {W_k y_k y_k^T W_k^T} {y_k^T W_k y_k} . \]
  
• BFGS minimizes \(\min || W - W_k ||_F\) such that \(W=W^T\): \[ H_{k+1} = H_k - \frac{ H_k y_k y_k^T H_k }{ y_k^T H_k y_k } + \frac{ s_k s_k^T }{ y_k^T s_k } \Leftrightarrow W_{k+1} = \left (I-\frac {y_k s_k^T} {y_k^T s_k} \right )^T W_k \left (I-\frac { y_k s_k^T} {y_k^T s_k} \right )+\frac{s_k s_k^T} {y_k^T s_k} . \]

In R, the so-called BFGS scheme is implemented in optim.

Another possible method (which is initially arised from quadratic problems) is the nonlinear conjugate gradients. This consists in computing directions \((d_0, \dots, d_k)\) that are conjugate with respect to a matrix close to the true Hessian \(H(x_k)\). Directions are computed iteratively by \(d_k = -g(x_k) + \beta_k d_{k-1}\) for \(k>1\), once initiated by \(d_1 = -g(x_1)\). \(\beta_k\) are updated according a scheme:

• \(\beta_k = \frac{ g_k^T g_k}{g_{k-1}^T g_{k-1} }\): Fletcher-Reeves update,
• \(\beta_k = \frac{ g_k^T (g_k-g_{k-1} )}{g_{k-1}^T g_{k-1}}\): Polak-Ribiere update.

There exists also three-term formula for computing direction \(d_k = -g(x_k) + \beta_k d_{k-1}+\gamma_{k} d_t\) for \(t<k\). A possible scheme is the Beale-Sorenson update defined as \(\beta_k = \frac{ g_k^T (g_k-g_{k-1} )}{d^T_{k-1}(g_{k}- g_{k-1})}\) and \(\gamma_k = \frac{ g_k^T (g_{t+1}-g_{t} )}{d^T_{t}(g_{t+1}- g_{t})}\) if \(k>t+1\) otherwise \(\gamma_k=0\) if \(k=t\). See Yuan (2006) for other well-known schemes such as Hestenses-Stiefel, Dixon or Conjugate-Descent. The three updates (Fletcher-Reeves, Polak-Ribiere, Beale-Sorenson) of the (non-linear) conjugate gradient are available in optim.

1.2.2. Computing the stepsize \(t_k\)

Let \(\phi_k(t) = f(x_k + t d_k)\) for a given direction/iterate \((d_k, x_k)\). We need to find conditions to find a satisfactory stepsize \(t_k\). In literature, we consider the descent condition: \(\phi_k'(0) < 0\) and the Armijo condition: \(\phi_k(t) \leq \phi_k(0) + t c_1 \phi_k'(0)\) ensures a decrease of \(f\). Nocedal & Wright (2006) presents a backtracking (or geometric) approach satisfying the Armijo condition and minimal condition, i.e. Goldstein and Price condition.

• set \(t_{k,0}\) e.g. 1, \(0 < \alpha < 1\),
• Repeat until Armijo satisfied,
  • \(t_{k,i+1} = \alpha \times t_{k,i}\).
• end Repeat

This backtracking linesearch is available in optim.

2.1.1. Theoretical value

The density of the beta distribution is given by \[ f(x; \delta_1,\delta_2) = \frac{x^{\delta_1-1}(1-x)^{\delta_2-1}}{\beta(\delta_1,\delta_2)}, \] where \(\beta\) denotes the beta function, see the NIST Handbook of mathematical functions https://dlmf.nist.gov/. We recall that \(\beta(a,b)=\Gamma(a)\Gamma(b)/\Gamma(a+b)\). There the log-likelihood for a set of observations \((x_1,\dots,x_n)\) is \[ \log L(\delta_1,\delta_2) = (\delta_1-1)\sum_{i=1}^n\log(x_i)+ (\delta_2-1)\sum_{i=1}^n\log(1-x_i)+ n \log(\beta(\delta_1,\delta_2)) \] The gradient with respect to \(a\) and \(b\) is \[ \nabla \log L(\delta_1,\delta_2) = \left(\begin{matrix} \sum\limits_{i=1}^n\ln(x_i) - n\psi(\delta_1)+n\psi( \delta_1+\delta_2) \\ \sum\limits_{i=1}^n\ln(1-x_i)- n\psi(\delta_2)+n\psi( \delta_1+\delta_2) \end{matrix}\right), \] where \(\psi(x)=\Gamma'(x)/\Gamma(x)\) is the digamma function, see the NIST Handbook of mathematical functions https://dlmf.nist.gov/.

2.1.2. R implementation

As in the fitdistrplus package, we minimize the opposite of the log-likelihood: we implement the opposite of the gradient in grlnL. Both the log-likelihood and its gradient are not exported.

lnL <- function(par, fix.arg, obs, ddistnam) 
  fitdistrplus:::loglikelihood(par, fix.arg, obs, ddistnam) 
grlnlbeta <- fitdistrplus:::grlnlbeta

3.1.1. Theoretical value

The p.m.f. of the Negative binomial distribution is given by \[ f(x; m,p) = \frac{\Gamma(x+m)}{\Gamma(m)x!} p^m (1-p)^x, \] where \(\Gamma\) denotes the beta function, see the NIST Handbook of mathematical functions https://dlmf.nist.gov/. There exists an alternative representation where \(\mu=m (1-p)/p\) or equivalently \(p=m/(m+\mu)\). Thus, the log-likelihood for a set of observations \((x_1,\dots,x_n)\) is \[ \log L(m,p) = \sum_{i=1}^{n} \log\Gamma(x_i+m) -n\log\Gamma(m) -\sum_{i=1}^{n} \log(x_i!) + mn\log(p) +\sum_{i=1}^{n} {x_i}\log(1-p) \] The gradient with respect to \(m\) and \(p\) is \[ \nabla \log L(m,p) = \left(\begin{matrix} \sum_{i=1}^{n} \psi(x_i+m) -n \psi(m) + n\log(p) \\ mn/p -\sum_{i=1}^{n} {x_i}/(1-p) \end{matrix}\right), \] where \(\psi(x)=\Gamma'(x)/\Gamma(x)\) is the digamma function, see the NIST Handbook of mathematical functions https://dlmf.nist.gov/.

3.1.2. R implementation

As in the fitdistrplus package, we minimize the opposite of the log-likelihood: we implement the opposite of the gradient in grlnL.

grlnlNB <- function(x, obs, ...)
{
  m <- x[1]
  p <- x[2]
  n <- length(obs)
  c(sum(psigamma(obs+m)) - n*psigamma(m) + n*log(p),
    m*n/p - sum(obs)/(1-p))
}