Method 1 (requiring starting point)

To compute the Hessian, we take the derivative of (7.7) with respect to $u_j$,

$$H_{db}^{u_i,u_j} = - \sum_{k=1}^N \; \left( \frac{{\rm d}\lambda_k}{{\rm d} \alpha_k}\right)^{-1} B_{k,j} B_{k,i}.$$ (7.9)

The following algorithm is proposed

1. Determine the pseudo-inverse solution ${\bf x}_p$ (5.11).
2. Seek an initial starting point for $\lambda$$_0$, a vector that meets (5.12) and has elements on $[0,1]$. Use an off-the-shelf linear programming solver as explained in Section 5.3.1.
3. Compute free variable ${\bf u}_0$ so that $\lambda$$_0={\bf x}_p + {\bf B}{\bf u}_0$ using ${\bf u}_0 = {\bf B}^\prime$   $\lambda$$_0$.
4. Determine $\alpha$ from $\lambda$ by solving (7.3) for $\alpha_i$, for each $\lambda_i$.
5. Compute entropy (7.6) and first and second derivatives (7.7),(7.9).
6. Take a Newton-Raphson iteration :
   $${\bf u}_{n+1} = {\bf u}_{n} - \left({\bf H}_{db}^{u,u}\right)^{-1} {\bf H}_{db}^{u},$$
   where ${\bf H}_{db}^{u,u}$ and ${\bf H}_{db}^{u}$ are the Hessian and gradient of $H_{db}$ with respect to ${\bf u}$.
7. Re-compute the mean : $\lambda$$_{n+1}={\bf x}_p + {\bf B}{\bf u}_{n+1}$. Check that all elements of $\lambda$$_{n+1}$ are in $[0,1]$. If not, take a smaller step.
8. Go to step 4.

The above algorithm is implemented by software/me_lin_01.m.