Re: [eigen] matrix::linspace

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2010/1/27 Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>:
> On Tue, 26 Jan 2010, Benoit Jacob wrote:
>
>> 2010/1/26 Hauke Heibel <hauke.heibel@xxxxxxxxxxxxxx>:
>>>
>>> This depends on your application. If you want something like this to work
>>>
>>> VectorXd x = VectorXd::LinSpaced(-1.0,1.0,10000);
>>> std:acos(x(0));
>>> std::acos(x(end));
>>>
>>> it better be in range or you'll end up with NaN when x(0)<-1 or
>>> x(end)>-1. But maybe that's too application specific.
>>
>> I'd call that a problem specific to acos(), not to the problem at hand!
>
> That depends how you define the problem at hand :)

OK let me explain. Normally, when you have a floating point number
that should be roughly 1, it shouldn't matter whether it is exactly 1
or is slightly larger. Or if it is <=1 or >1. Exact comparisons
shouldn't matter. There's one exception, that's exact comparisions
with 0, x==0, due to the special role of 0 in floating point
arithmetic, but all other floating point numbers fit in a continuum
inside of which small errors are considered the rule and floating
point arithmetic, being continuous (or better, you could say,
Lipschitz) means that small errors in the input give small errors in
the output. But acos() is special in that it is sensitive to arbitrary
small errors. acos(1)=0, acos(1+epsilon)=NaN. This is exceptional in
floating point arithmetic: almost every usual other operation/functor
is continuous away from 0, that is, f(1+epsilon) ~= f(1) up to (a
small multiple of) machine precision. This special behavior of acos()
forces people to keep special rules in mind just for that function,
which is a bit crazy! Fortunately its use is not too frequent.

If it were up to me, I'd have implemented acos() just returning the
pade approximant no matter if x>1 or x<-1. That would even be a bit
more efficient as it would mean two if's less. My acos() would be just
spec'd as "any continuous function on R that agrees with acos on
[-1,1]". The current behavior of returning NaN outside of [-1,1] is a
bit useless in practice, if theoretically satisfying. People who want
to be careful could just do the test on x, not on the result.

>
> I'd define it as, LinSpaced(start, end, size) returns a vector with 'size'
> elements whose first element is 'start' and last element is 'end' such that
> the distance between consecutive elements is constant up to round-off error.

There is no problem with static LinSpaced, as that uses the
random-access version of the functor!
The problem is only with setLinSpaced(), the optimized sequential version.

> I wonder, in what practical situation is the performance of linspace a
> bottleneck? If there are none - and I cannot think of any - then perhaps
> accuracy is more important than performance.

Good observation. LinSpaced() indeed has good accuracy and not optimal
performance. For setLinSpaced(), since Hauke spent so much effort
vectorizing it, I'd assume he has goods reasons to. Fortunately, see
my other mail, I still think that we can have good performance and
good precision.

Benoit

>
>
> Jitse
>
>
>