Make `ComplexField::norm1` return the sum of the absolute values of all coefficients

[afossa]

For the types that implement ComplexField, the implementation of norm1 delegates to abs.

I would propose to change the implementation such that it computes the sum of the absolute values of all coefficients.

• for Dual this is simply re.abs() + eps.abs().
• for DualVec, this is re.abs() plus the sum of all e.abs(), where e are the entries of eps.
• for Dual2 and Dual2Vec, we need to agree on how to interpret the 1-norm of v2. I would take the sum of the abs values scaled by 1/2, such that it coincides with the sum of the abs values of the second-order coefficients of the Taylor expansion.

This new behavior would then mimic, for example, that of GTPSA.normTPS. It would also allow to use its output as convergence criterion in iterative algorithms that update dual numbers recursively until some convergence is achieved.

[g-bauer]

Hi @afossa. Sorry for the delayed response. We are a bit time constrained at the moment but will look into this as soon as time permits. The application of norm1 to convergence criteria is certainly interesting and relevant.

[prehner]

Hi @afossa, thank you for initiating this discussion! I see where you are coming from, but I am not fully convinced yet.

First of all, we inherit a lot of mess from nalgebra and simba and it doesn't look like dimforge is particularly active or interested in feedback (this issue, which relates to the num-dual crate is three days away from its third anniversary now: dimforge/simba#44). I would honestly love to get rid of all this mess of implementing ComplexField on data structures that do not even represent complex numbers and which somehow has to include SIMD support without actually being able to use it properly (exactly @cormacrelf's point in the issue I linked above). On the other hand it is nice to be able to use the linalg functions in nalgebra directly with dual numbers. Other crates like faer might be a nice alternative, but their analogous trait might be even more inaccessible. So, I guess we're stuck with it for now.

Now, regarding the norm1 function. I'm not really sure, where this is supposed to used in the first place. Is it basically meant as an extension of abs to complex numbers? And then used specifically for checking convergence of iterative linear (!) solvers? In that case your suggestion makes sense, but it could also confuse people, because in the spirit of the trait system of simba all the dual numbers are essentially real numbers, as opposed to complex numbers. That means it could be counterintuitive, if the value of norm1 depends on the derivative you are calculating. In GTPSA the function is calles normTPS which makes it a bit more transparent what it does. So maybe, we should instead have a norm_dual function as part of the DualNum trait?

On a related note: We should really also put all of the nalgebra stuff in macros.

[afossa]

Hi @prehner, thanks for sharing your opinion.

I don't know what norm1 is meant for in simba, but I can expand a bit more on my use case. The idea is to implement algorithms that internally solve one or more nonlinear equations, and have them generic over the input/output types T: ComplexField such that they can be evaluated with e.g. f64 and dual numbers. I know that num-dual provides methods to compute the derivatives of implicit functions, but they are not really applicable here because the "outer" algorithm as a whole is effectively an explicit function (please correct me if I am wrong).

Regarding the implementation of a dedicated method like norm_dual, I would still prefer the other solution and override norm1, but I understand your point. However, there are libraries that implement norm as I am suggesting:

• DACE::DA
• TaylorSeries.jl overrides LinearAlgebra.norm.
• Hipparchus for e.g. DerivativeStructure and Gradient. The default implementation for non-dual-number-like types calls abs just like in simba.

On the side, I fully agree with you about the nalgebra bloat. Most of the linear algebra could be done on a trait that sits in between the too generic Scalar and ComplexField.

[prehner]

That's interesting. Good to see more people going into the automatic differentiation of complex models. Do you need your type to be generic over T: ComplexField or could it also be T: DualNum<f64> + Copy? In either case, what I can suggest is to actually only solve the real part with the non-linear solver. Calculating derivatives is always a linear problem, which implies that if you use a second order solver (at least for the last iterations), the (first) derivatives will converge in a single iteration. So the approach that we have been using very successfully is

• use any non-linear solver that robustly solves your problem for f64 to a reasonable tolerance
• use an undamped Newton solver to obtain very low residuals in a few steps
• do one additional step of the Newton solver, now with dual numbers to recover the (first) derivatives of your variables

That's why, we didn't need a way to check the convergence of the dual parts so far. If you have many iterations, carrying around the dual parts all the time will increase the computational effort significantly without any benefit. The only case where I could see that is when it for some reason is infeasible to use second order methods, maybe that's the case in your problem?

That being said, I am not against implementing a norm metric in num-dual. As you showed nicely, most similar packages have something like that, so we should adapt as well. Code-wise, I would like to generalize the implementation of the ComplexField bloat into a separate macro, so that we have it for every data type. Then, we need to address the definition of the norm for higher-order derivatives, and can change that accordingly in the macro.

[afossa]

Hi @prehner, thanks for the suggestions. What you propose is indeed what I have been doing to compute the high-order Taylor expansion of the solution to an implicit equation using libraries like DACE and more recently GTPSA. In fact, starting from the real solution to the equation, a simple fix point iteration scheme suffices to converge the k-th order derivatives in just k steps. But I sometimes encountered resistance against this method because someone prefers having both f64 and duals following the same code path.

To answer your question, I think that the metric could be implemented for DualNum<F> + Copy if that's easier.