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))
}