Method 1 (requiring starting point)

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

$\displaystyle 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}.$ (6.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 (6.3) for $ \alpha_i$, for each $ \lambda_i$.
  5. Compute entropy (6.6) and first and second derivatives (6.7),(6.9).
  6. Take a Newton-Raphson iteration :

    $\displaystyle {\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.

Baggenstoss 2017-05-19