RFC: stability of DenseVector::two_norm()
This is the current implementation of DenseVector::two_norm()
which you will find here
//! two norm sqrt(sum over squared values of entries)
typename FieldTraits<value_type>::real_type two_norm () const
{
typename FieldTraits<value_type>::real_type result( 0 );
for (size_type i=0; i<size(); i++)
result += fvmeta::abs2((*this)[i]);
return fvmeta::sqrt(result);
}
There's the potential for an avoidable under- or overflow in this code as explained e.g. here in detail or proven by the following snippet:
#include <cmath>
#include <cstdio>
#include <tuple>
int
main() {
double x, y;
std::tie(x, y) = std::make_pair(3e-200, 4e-200);
printf("%5g vs. %5g\n",
std::sqrt(x*x + y*y),
std::hypot(x, y));
std::tie(x, y) = std::make_pair(3e+200, 4e+200);
printf("%5g vs. %5g\n",
std::sqrt(x*x + y*y),
std::hypot(x, y));
}
which prints
0 vs. 5e-200
inf vs. 5e+200
for me: here, when std::sqrt(x^2 + y^2)
under/overflows, std::hypot(x, y)
does not.
Avoiding the under/overflow comes at a cost: The hypot
-based version appears to be slower than the sqrt
version by roughly a factor of 4 as illustrated by this R/C++ microbenchmark:
> Rcpp::cppFunction('double cheap(DoubleVector x) { double ret = 0; for (int i=0; i < x.size(); ++i) ret += x[i] * x[i]; return std::sqrt(ret); }')
> Rcpp::cppFunction('double expensive(DoubleVector x) { double ret = 0; for (int i=0; i < x.size(); ++i) ret = std::hypot(ret,x[i]); return ret; }')
> x <- runif(1000000); microbenchmark::microbenchmark(cheap(x),expensive(x))
Unit: milliseconds
expr min lq mean median uq max neval
cheap(x) 4.229765 4.244638 4.433328 4.361079 4.42705 6.30066 100
expensive(x) 16.178626 16.219625 16.542381 16.291846 16.55490 20.25674 100
>
So what I'd like to ask is:
- Does this stability issue seem relevant to anyone?
- If so, is the decrease in performance by a factor of 4 for this method acceptable?
An implementation could look like this:
//! two norm
typename FieldTraits<value_type>::real_type two_norm () const
{
typename FieldTraits<value_type>::real_type result( 0 );
for (size_type i=0; i<size(); i++)
result = fvmeta::hypot(result, (*this)[i]);
return result;
}
except that this does not answer the question of what fvmeta::hypot
would forward to for a real_type
other than float
, double
, or long double
(since e.g. std::hypot
could be used for those three but not e.g. a GMP field).