Method 2 (not requiring starting point)

To solve (5.15), define $\alpha$ as the element-wise inverse of $\lambda$. Then, (5.15) can be written ${\bf B}$$\alpha$$={\bf0}$, which means that $\alpha$ must be fully in the column space of ${\bf A}$, or that

$$\mbox{\boldmath$\alpha$} = {\bf A} {\bf v},$$ (5.17)

for some $D\times 1$ vector ${\bf v}$. Therefore, to find the mean of the surrogate density, we solve for the free variable ${\bf v}$ such that

$${\bf A}^\prime \mbox{\boldmath$\lambda$} ({\bf A} {\bf v}) = {\bf z}^*,$$ (5.18)

where

$$\mbox{\boldmath$\lambda$} ( \mbox{\boldmath$\alpha$} ) = [1/\alpha_1, 1/\alpha_2, \ldots 1/\alpha_N]^\prime.$$ (5.19)

We can solve this by driving the square error to zero:

$$\rho({\bf v})=({\bf A}^\prime \mbox{\boldmath$\lambda$} ({\bf A} {\bf v}) - {\bf z})^\prime ({\bf A}^\prime \mbox{\boldmath$\lambda$} ({\bf A} {\bf v}) - {\bf z}).$$ (5.20)

We can easily find derivative of $\rho({\bf v})$ with respect to ${\bf v}$ :

$$\left[ \frac{\partial \rho}{\partial v_k}\right] = 2 ({\bf A}^\prime \mbox{\boldmath$\lambda$} ({\bf A} {\bf v}) - {\bf z})^\prime {\bf A}^\prime {\mbox{\boldmath$\Lambda$}} {\bf A},$$ (5.21)

where $\Lambda$$($$\alpha$$)$ is the diagonal matrix with diagonal elements

$$\Lambda_i = -1/\alpha_i^2.$$ (5.22)

We have found that if we use the negative-definite Hessian approximation

$$\left[ \frac{\partial^2 \rho}{\partial v_k \partial v_l}\right] \... ...A}^\prime {\mbox{\boldmath$\Lambda$}(\mbox{\boldmath$\alpha$})} {\bf A}\right),$$ (5.23)

the resulting Newton-Raphson algorithm has excellent convergence properties when starting with $\alpha$$= [1, 1, \ldots 1]^\prime$. Let $\hat{\mbox{\boldmath $\lambda$}}({\bf z}^*)$ be the value of $\lambda$ at the solution to (5.18). The algorithm to find $\hat{\mbox{\boldmath $\lambda$}}({\bf z}^*)$ can be summarized as follows.

1. Set iteration counter $n=0$.
2. To initialize, let ${\bf v}_0=\left({\bf A}^\prime {\bf A}\right)^{-1} {\bf A}^\prime {\bf 1},$ where ${\bf 1}$ is the vector of ones.
3. $\alpha$$_n = {\bf A} {\bf v}_n$. Initially, $\alpha$ will be the vector of ones.
4. $\lambda$$_n($$\alpha$$_n) = [1/\alpha_1, 1/\alpha_2, \ldots 1/\alpha_N]^\prime.$
5. Compute derivative and Hessian according to (5.21),(5.22) and (5.23), then update ${\bf v}$:
   $${\bf v}_{n+1} = {\bf v}_n + \left[ \frac{\partial^2 \rho}{\partial u_k \partial u_l}\right]^{-1} \left[ \frac{\partial \rho}{\partial u_k}\right].$$

6. Increment $n$ and go to step 3.

The algorithm is implemnted by software/lam_solve.m with dbound=0. Although this method corresponds to classical methods, it is based on a novel sampling argument and can be extend to other manifolds as we will see.