I have implemented 2 main additions to the standard SRIF. An inherent downside to the SRIF is the number of operations necessary to generate the Householder matrix \(H\) grows like \(m^2(n-m)\). My implementation of the SRIF runs the matrix triangularization simultaneously on multiple cores, which results in a significant speed-up. This can be done in a thread-safe manner because each Householder operation functions on a column-by-column basis.
The second addition I have made to the standard SRIF is the inclusion of a secondary \(\chi^2\) minimization for the scaling of \(\vec{\delta{x}}\). So \[Q'' = Q(\vec{x} + \alpha \vec{\delta x})\] is minimized. This minimization is done with a grid search over 6 decades. The additional minimization adds a trivial additional computation cost to the overall minimization of \(\chi^2\), but allows for faster convergence, and the possibility of skipping over local minima in the \(\chi^2\)-space.
The SPF routine can currently fit models to data while taking into account a suite of “penalty functions”, to ensure that the models that are generated are physical. In a way, these penalty functions serve to make the fits operate in a more global context – there may be a local minimum in \(\chi^2\)-space which the fitting algorithm would want to tend towards, but that minimum can be ruled out a priori thanks to physical limitations. These penalty functions include limits on ellipsoid axis ratios, constraints to surface concavities, limits on the model center of mass distance from the image center of mass, as well as others. I have implemented these penalty functions in the SRIF framework. Specifically, I have done so by redefining the residual vector as \[\vec{R}'' = \left( \begin{array}{c} \vec{z} - \vec{m} \\ \vec{p_w} \end{array} \right)\] Where \[\vec{p_w} = \left( \begin{array}{c} p_1\times w_1 \\ \vdots \\ p_N\times w_N \end{array} \right)\] for \(\{p_i\}\) , \(\{w_i\}\) are the set of penalty functions and penalty weights, respectively, and \[A'' = -\frac{\partial \vec{R''}}{\partial \vec{x}}\] The algorithm then progresses as described in section ?, with \(\vec{R}\) replaced with \(\vec{R}''\) and \(A''\) replacing \(A\).