Re: nchoosek warning, checks and tests

[Jaroslav Hajek]

On Wed, Dec 10, 2008 at 4:56 PM, Francesco Potortì <address@hidden> wrote:
>>> Discussion:
> ...
>>> - log of sum is used instead of product so that overflow is further away
> ...
>>>   if (n == 1)
>>> -    if (k > v/2)
>>> -      k = v - k;
>>> +    k = min (k, v-k);          # improve precision at next step
>>> +    A = round (exp (sum (log (v-k+1:v)) - sum (log (2:k))));
>>> +    if (A > 1/(exp (2*k*eps) - 1))
>>> +      warning ("nchoosek", "nchoosek: possible loss of precision");
>>>     endif
>>> -    A = round (prod ((v-k+1:v)./(1:k)));
>
>>Why did you replace
>>round (prod ((v-k+1:v)./(1:k)))
>>by
>>round (exp (sum (log (v-k+1:v)) - sum (log (2:k))))
>>?
>
> I had included a "Discussion" section at the beginning of my mail
> describing a short rationale for every change.  Please have a look at it
> also for discussion of other changes before installing this changeset,
> as some changes may be objectionable, and feel free to change what you
> think is wrong or suboptimal.
>
> For this specific change, here is a more in-depth discussion.
>
> When comparing nchoosek to bincoeff, it turns out that the first one is
> thought for integer and not very big numbers, while the second is
> tailored for fractional and big numbers.  In fact, bincoeff uses the
> gamma function, which allows real arguments, and is faster for big
> numbers, because nchoosek creates an array to be summed or multiplied.
>
> So I think it is important for nchoosek to emit a warning when the
> result is not exact (Matlab does the same).  However, as far as I can
> tell, the error bounds for the two expressions above (the prod and the
> sum of logs) is the same, so using one or the other is the same, as long
> as we want to have an exact result.  On the other hand, the sum of logs
> avoids overflow in some cases, such as v=300, k=140.

No, the existing expression doesn't overflow for these numbers. In
fact it simply cannot overflow if the result doesn't overflow.

>
> Important: this assumes that the log and exp functions produce a
> relative error near eps, which I am not sure is true in Octave's
> implementation.  If this assumption is wrong, then the prod
> implementation is preferable.
>

Whether true or not, the exp-log implementation calculates 2*k
logarithms and one exponential, which is certainly slower, not
speaking about more room for errors to accumulate. On the other hand,
my expression does k divisions and multiplications.


> To summarise, both implementations give the same results as long as we
> are interested in an exact result.  If we can live with an approximate
> result, the prod implementation is on the average slightly more accurate
> while the sum of logs can handle some more cases without overflow.
>



> Some more considerations:
>
> rethinking about it, the above expression
>  if (A > 1/(exp (2*k*eps) - 1))
> is equivalent with the much simpler
>  if (A*2*k*eps > 1)
> moreover, this is not entirely accurate and should be instead:
>  if (A*2*k*eps > 0.5)
>
> Note that Matlab is not so picky about rounding errors, and as long as
> its help page is accurate, gives a warning only when (A > 1e15).
>
> Please let me know your comments.
>

Actually, I don't exactly get the reason for the warning, because just
about every computation in Octave (and Matlab) is inexact. That's just
a property of floating-point numbers. Surely we don't warn at "sin(1)"
that the result is only accurate to 15-16 decimal places...
but if it's in Matlab, then why not. Maybe it's for naive users who
could expect the result to be computed exactly.


> --
> Francesco Potortì (ricercatore)        Voice: +39 050 315 3058 (op.2111)
> ISTI - Area della ricerca CNR          Fax:   +39 050 315 2040
> via G. Moruzzi 1, I-56124 Pisa         Email: address@hidden
> (entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/
>



-- 
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz