Re: [eigen] Specializing max_coeff_visitor for some number types

[Marc Glisse <marc.glisse@xxxxxxxx>]

On Mon, 7 Apr 2014, Christoph Hertzberg wrote:

> On 07.04.2014 12:30, Marc Glisse wrote:
> > Hello,
> >
> > as far as I understand, Eigen uses max_coeff_visitor internally to
> > determine the best pivot to use. And its notion of "best" is the largest
> > in absolute value, which makes a lot of sense for float/double. For
> > other number types, I may have a different notion of best. For example,
> > for a rational, I may want to pick the "simplest" non-zero coefficient,
> > to avoid an explosion in coefficient size. For interval arithmetic, I
> > may want to compare only the lower bounds of the absolute values, or
> > maybe involve the width somehow, but if we call operator< on [.9,1.1]
> > and [.8,1.2], we have lost.
> >
> > As a quick hack, it looks like I can specialize scalar_abs_op to return
> > something other than the absolute value, but that's ugly and I have to
> > hope that abs is not used anywhere else in Eigen for other purposes.
> >
> > Any advice?
>
> I would object redefining max_coeff_visitor or scalar_abs_op only to make 
> decompositions work.
> E.g., if I have an array of intervals
>  [ (0.15, 0.3), (-0.1, -0.05), (-0.3, 0.2) ]
> I would like .abs() return:
>  [ (0.15, 0.3), (0.05, 0.1), (0, 0.3) ]
> And maxCoeff(Index*) should not be possible, since the type is not 
> well-ordered (either the first or the last could be the biggest).

We agree, that's why I said "quick hack", "ugly" and "hope that".

> maxCoeff() is possible, on the other hand, because
>  max((a,b), (c,d))==(max(a,c), max(b,d))
>
> I would rather suggest to out-factor a method find_best_pivot which then can 
> be specialized for special types.

At what point exactly should it split? The easiest would be that LU calls:

DenseBase<Derived>::bestPivotCoeff(IndexType* index) const

which visits with a best_coeff_visitor, which uses 
internal::scalar_better_coeff_op, which by default compares the absolute 
values. But then we lost the parallelism in computing the absolute values.

We could also always call cwiseAbs before bestPivotCoeff, although that 
might be a waste for some types that don't need to compute the absolute 
value to compare the coefficients.

bestPivotCoeff could just pass the matrix and index to some external 
class, but then users who just want to change the comparison of 2 numbers 
would need to write a visitor and all.

I can imagine some other solutions, but nothing that feels great.

> For interval types, I guess the element with the biggest lower bound is 
> generally to be preferred (however it depends on if you want to determine the 
> best possible decomposition for a given matrix, or the worst-case scenario 
> what could happen when decomposing with standard floating point types).

I want the best in my case (well, not necessarily the best, but one good 
enough). And the biggest lower-bound is just a simple heuristic, given the 
choice between (1,1) and (2,20) I'd pick the first interval as the better 
pivot ;-) But it may not be worth doing anything more subtle.

--
Marc Glisse