# Whittaker-Shannon interpolation: Accuracy dies with speedup; can it be fixed?

With a truncated Whitaker-Shannon series (cardinal series) $$f(t) = \sum_{j = 0}^{n-1} y_{j} \frac{\sin\left(\pi( \frac{t-t_0}{h} -j)\right)}{\pi\left(\frac{t-t_0}{h}-j\right)}$$ we can naively evaluate the sum by repeated calls to sinc routines, such as the following code does:

   Real naive_sum(Real t) const {
using boost::math::constants::pi;
Real f = 0;
for (size_t i = 0; i < m_y.size(); ++i) {
Real arg = pi<Real>()*( (t-m_t0)/m_h - i);
f += m_y[i]*boost::math::sinc_pi(arg);
}
return y;
}


However, repeated calls to sinc are expensive, so we can use the identity $$\sin(\theta - j\pi) = (-1)^{j}\sin(\theta)$$ to write this as $$f(t) = \frac{\sin(\pi(t-t_0)/h)}{\pi} \sum_{j=0}^{n-1} (-1)^{j}\frac{y_{j}}{(t-t_0)/h -j}$$

which can be implemented as follows:

    Real operator()(Real t) const {
using boost::math::constants::pi;
using std::sin;
Real y = 0;
Real x = (t - m_t0)/m_h;

for (size_t i = 0; i < m_y.size(); ++i)
{
Real denom = (x - i);
if (denom == 0) {
return m_y[i];
}
if (i & 1) {
y -= m_y[i]/denom;
}
else {
y += m_y[i]/denom;
}
}
return y*sin(pi<Real>()*x)/pi<Real>();
}


However, I have observed a vast decrease in accuracy using the fast method over the slow method. Can the speed of the fast method be preserved without a massive decrease in accuracy?

Working code, for those that care:

#ifndef BOOST_MATH_INTERPOLATORS_WHITAKKER_SHANNON_HPP
#define BOOST_MATH_INTERPOLATORS_WHITAKKER_SHANNON_HPP
#include <boost/math/special_functions/sinc.hpp>
#include <boost/math/constants/constants.hpp>

namespace boost::math::interpolators {

template<class RandomAccessContainer>
class whittaker_shannon {
public:

using Real = typename RandomAccessContainer::value_type;
whittaker_shannon(RandomAccessContainer&& y, Real t0, Real h) : m_y{std::move(y)}, m_t0{t0}, m_h{h}
{
}

Real operator()(Real t) const {
using boost::math::constants::pi;
using std::sin;
Real y = 0;
Real x = (t - m_t0)/m_h;

for (size_t i = 0; i < m_y.size(); ++i)
{
Real denom = (x - i);
if (denom == 0) {
return m_y[i];
}
if (i & 1) {
y -= m_y[i]/denom;
}
else {
y += m_y[i]/denom;
}
}
return y*sin(pi<Real>()*x)/pi<Real>();
}

Real naive_sum(Real t) const {
using boost::math::constants::pi;
Real y = 0;
Real s = pi<Real>()*(t-m_t0)/m_h;
for (size_t i = 0; i < m_y.size(); ++i) {
Real arg = pi<Real>()*( (t-m_t0)/m_h - i);
y += m_y[i]*sinc_pi(arg);
}
return y;
}

Real operator[](size_t i) const {
return m_y[i];
}

private:
RandomAccessContainer m_y;
Real m_t0;
Real m_h;
};
}
#endif


Here's a test that reproduces the phenomenon:

template<class Real>
void test_bump()
{
using std::exp;
using std::abs;
auto bump = [](Real x) { if (abs(x) >= 1) { return Real(0); } return exp(-Real(1)/(Real(1)-x*x)); };

Real t0 = -1;
size_t n = 2049;
Real h = Real(2)/Real(n-1);

std::vector<Real> v(n);
for(size_t i = 0; i < n; ++i) {
Real t = t0 + i*h;
v[i] = bump(t);
}

auto ws = whittaker_shannon(std::move(v), t0, h);

std::mt19937 gen(323723);
std::uniform_real_distribution<long double> dis(-0.95, 0.95);

size_t i = 0;
while (i++ < 1000)
{
Real t = static_cast<Real>(dis(gen));
Real expected = bump(t);
if(!CHECK_MOLLIFIED_CLOSE(expected, ws(t), 50*std::numeric_limits<Real>::epsilon())) {
std::cerr << "  Problem occured at abscissa " << t << "\n";
}
}
}


• Have you tried (1) Accumulating positive and negative terms separately and subtracting at the end (2) Compensated addition (Kahan summation)? Apr 2 '19 at 18:20
• I believe (1) is discussed by Higham where he demonstrates it doesn't work; see doi.org/10.1137/0914050, I just tried (2), and it doesn't help. The summation condition number is not large at the abscissas where the unit tests fail (I computed it around to be ~5 at the problem abscissas.) There are abscissas where the condition number is ~3000, at those points, the unit tests pass. Apr 2 '19 at 19:48
• (1) Make sure you are compiling with strictest adherence to IEEE-754 (for my Intel compiler that is /fp:strict for example). (2) Try explicitly using of fma() as extensively as possible (convert divisions to multiplication if need be) to guard against subtractive cancellation when combining products. (3) Be skeptical of your reference function, as it may have numerical issues, too. I would suggest using triple the precision of the actual computation you are trying to implement. Apr 2 '19 at 20:38
• The compilation is with g++-8 -O3 -march=native; none of these flags break IEEE compliance. I feel like my reference function has no need for skepticism, esp. since I draw the test values away from it's problematic range as |x| -> 1. Apr 2 '19 at 23:16
• I have repro and I think I have identified the culprit. sin of large arguments polluted by rounding error. Suggest switching sin to sin_pi: return y*sin(pi<Real>()*x)/pi<Real>(); --> return y*boost::math::sin_pi(x)/pi<Real>(); Note: I used my own sinpi from this answer instead of Boost. Apr 3 '19 at 3:55

I was able to reproduce the behavior reported in the question, and traced the observed inaccuracies to the following line:

return y*sin(pi<Real>()*x)/pi<Real>();


The explicit multiplication with a floating-point approximation of π introduces a small error into the argument to sin, which comprises the representational error in the constant and the rounding error added by the multiply and is on the order of 1 ulp. This small error in the argument represents a phase error in sin which grows with the magnitude of x. In this case |x| is on the order of 1000, resulting in quite a bit of error magnification.

Because this is a relatively frequent scenario, various platforms offer a sinpi function that computes sin (π x) with the multiplication by π happening inside sinpi after argument reduction, minimizing the phase error. In the specific case of Boost the desired function is boost::math::sin_pi. Often the sinpi function is also more efficient than the regular sin function, since its internal argument reduction is simpler.

Replacing the original source line with the line below should fix the observed accuracy issue completely:

return y*boost::math::sin_pi(x)/pi<Real>();

• Not only did this solve the problem, the resulting accuracy increases by a factor of 5 over the naive sum. I'll need to keep sin_pi in mind. Apr 3 '19 at 13:05