GaussNewton params and cost out of sync in state

[gmilleramilar]

    fn next_iter(
        &mut self,
        problem: &mut Problem<O>,
        state: IterState<P, (), J, (), F>,
    ) -> Result<(IterState<P, (), J, (), F>, Option<KV>), Error> {
        let param = state.get_param().ok_or_else(argmin_error_closure!(
            NotInitialized,
            concat!(
                "`GaussNewton` requires an initial parameter vector. ",
                "Please provide an initial guess via `Executor`s `configure` method."
            )
        ))?;
        let residuals = problem.apply(param)?;
        let jacobian = problem.jacobian(param)?;

        let p = jacobian
            .clone()
            .t()
            .dot(&jacobian)
            .inv()?
            .dot(&jacobian.t().dot(&residuals));

        let new_param = param.sub(&p.mul(&self.gamma));

        Ok((state.param(new_param).cost(residuals.l2_norm()), None))
    }

The code above from gaussnewton_method.rs. Here, residuals are calculated calling apply(params). Those residuals are used to create a new parameter vector. However there is no new set of residuals calculated. So at the end of the function, the new parameter vector and the old residuals (reduced to a cost) are stored in the state. This causes problems for early termination schemes dependent on a correct cost value (like tolerance).

[stefan-k]

Hi,

thanks for reporting this! That's a very good point. I'd like to avoid computing the residuals again at the end of the function just for the cost. I think IterState needs to be able to carry the residuals as well. I'll have to think about this.

[gmilleramilar]

Yeah, I think that's probably right. In our application re-computing the residuals is quite expensive. Would this require you to make IterState generic over the type of the residuals?

[stefan-k]

Sorry that I wasn't able to work on this, I'm swamped with other responsibilities lately.

Would this require you to make IterState generic over the type of the residuals?

Yes, probably. Well, at least I think that this would be the approach that would make the most sense. Alternatively one could just say that residuals must be of the same type as the parameter vectors. This will probably work well in most cases but will be a difficult to justify limitation in other cases.

[stefan-k]

Thanks again for the bug report and for the initial implementation of residual handling in IterState (#343). I've extended upon your work and fixed the bug in GaussNewton in #392. The stopping criteria should now work properly. Let me know if you run into any issues :)