Refinement

It should be mentioned that for the sake of numerical stability, the measured points $\vec{Y_i}$ should be translated if necessary into a coordinate system whose origin is close to the mean of the points. This stability can be achieved automatically, and the computation simplified somewhat, by reformulating the solution in the following way.

First, solve for $a_m$ from the last of the equations (13):

$$a_m = \overline{Y}_{m} - \sum_{k=1}^{m-1} a_k \overline{Y}_{k}$$ (14)

where

$$\overline{Y}_{k} = { {\textstyle \sum_i W_i Y_{ik}} \over {\textstyle \sum_i W_i}}, \quad k=1, \ldots, m$$ (15)

This shows that the point $\overline{\vec{Y_i}}$ lies on the best-fit hyperplane. Now define $\vec{Y_i}' = \vec{Y_i}- \overline{\vec{Y}}$ and $\vec{z_i}= \vec{y_i}- \overline{\vec{Y}}$. From Equations (8) and (9) we have

$$z_{ij} = Y_{ij}' - W_i a_j \sigma_{ij}^2 f(\vec{Y_i}), \quad i=1, \ldots, n,\ j=1, \ldots m-1$$ (16)

where we can express $f(\vec{Y_i})$ as

$$f(\vec{Y_i}) = \sum_{k=1}^{m-1} a_k Y_{ik}' - Y_{im}' \quad i=1, \ldots, n$$ (17)

Then upon inserting (14) into the remaining equations (13) we find

$$\begin{array}{lclcl} {\displaystyle \sum_{k=1}^{m-1}} (\sum_... ..._i W_i z_{ij} Y_{im}', \quad j = 1, \ldots, m-1 \\ \end{array}$$ (18)

This reformulation has improved the numerical stability and reduced the order of the set of equations that needs to be solved on each iteration by 1.