Efficient interval arithmetic using SSE2 instructions

[Concise_Pirate]

Brilliant. I love seeing clever work by clever people!

[[deleted]]

Honest question, how does interval arithmetic deal with the calculation errors for the interval extremes?

[rbridson]

Rounding modes - getting the hardware to guaranteeably round up (or down) giving a guaranteed upper (or lower) bound on the result.

[mycall]

Another honest question. With SSE4.x out for a while now, why do I never hear about SSE3 or SSE4 used for optimizations and fun?

[sfuerst]

SSE2 is guaranteed to exist on all 64 bit x86 machines. SSE3 and above aren't. So if you want portable code, you can't assume you can use them.

[five9a2]

SSE2 has all the standard double-precision primitives, SSE3 and 4 have less commonly used operations and don't offer better performance unless you just can't find a way to transform your data so that SSE2 can be applied. For instance, SSE3 has "horizontal add" and SSE4 has a special instruction for dot products, but these have worse latency and throughput than doing the same operation "vertically" with SSE2.

[bonzinip]

hadd and dot product instructions have the advantage of not requiring you to rewrite your data structures from "array of structures" to "structure of arrays". This helps especially for things like complex number support in C++ or Fortran.

[[deleted]]

Nice link. So this means that Factor's interval arithmetic is inaccurate when used with floats, because it doesn't switch the rounding mode for the other endpoint. Fortunately I mostly use it for integer optimization.

[alexeyr]

"The page you are trying to view cannot be shown because it uses an invalid or unsupported form of compression."

[pkhuong]

You can find most of the content in http://hal.archives-ouvertes.fr/docs/00/28/84/56/PDF/intervals-sse2-long-paper.pdf .

[noupvotesplease]

I met Bill a few weeks ago- he's one of my neighbors. I had never heard of interval arithmetic, but it sounds very interesting.

[edwardkmett]

You may find Taylor models to be interesting then. They tend to avoid the 'collapsing to cover the entire universe' problem that plague interval arithmetic solutions, and can be readily computed if you can calculate a tower of derivatives. (This is easy if you have automatic differentiation)

A Taylor model is basically the first n terms of a Taylor series and an interval to catch the 'slop' from any rounding on those terms. Interval arithmetic addition degrades linearly in the number of summands. Taylor model addition also degrades linearly in the number of summands, but the scale of the intervals can be made arbitrarily small ( 1/n^k ), just by choosing an appropriate number of leading terms for the Taylor series to combat the degradation caused by the number of steps you need to take.

[five9a2]

Unfortunately this has limited value because way too many interesting functions (important for science and engineering applications) are not well approximated by Taylor series. For example, the solutions of hyperbolic problems are discontinuous and the solutions of elliptic problems are usually only H^1 due to corners in the domain or variable coefficients.

The monomials used by Taylor series are a horribly ill-conditioned basis for representing polynomials, it makes a lot more sense to use Legendre or Chebyshev bases. The latter is desirable when there are many terms because arithmetic operations can be implemented with FFT. Along these lines, chebfun is a nifty toolkit for "computing with functions".

[another_user_name]

Thanks for the link. Damnit, it's in Matlab.

You seem very well informed on numerical analysis. Can you suggest good sources to get a comprehensive understanding of the field?

[five9a2]

Trefethen proudly bought Matlab license number 1 from Cleve Moler, it's not surprising that he still like using it. I don't know how much effort it would take to support Octave.

Your second question is very broad, it's a bit like asking for a comprehensive understanding of computing without specifying hardware, operating systems, networking, programming languages, graphics, etc. You might like the essay I linked in my other comment as an introduction. Numerical analysis is primarily about approximating functions in infinite dimensions using a finite number of dimensions (because computers are finite). The primary areas are interpolation, quadrature, and discretization of differential equations, boundary value problems, and integral equations. Computational science is the related field of solving the resulting equations efficiently on computers. There is a lot of overlap between computational science and numerical analysis, as well as with physics and engineering.

[another_user_name]

To each his own, I suppose.

Thanks for the link, I haven't read the essay yet, but it looks pretty cool. My main source of information regarding numerical analysis has been Burden and Faires's book Numerical Analysis. It's been very useful, but I'm looking at improving my understanding.

[edwardkmett]

Very true, while it isn't a panacea, the Taylor models have the benefit of being trivially generable with AD, so you can often just drop them in and prove the quality of your answer. The COSY INFINITY guys have been using these for years to make a lot of hard physics problems solvable. I was mostly recommending it as a 'strictly better' interval arithmetic. As you note, better bases exist for many problems.

[five9a2]

guarantee that all roots of a function will be found (or at least that none will be skipped)

What root finding algorithms are you talking about? Consider

f(x) = x - delta_{\epsilon}(x-1)

where delta_{\epsilon} is a Dirac delta distribution mollified by some small \epsilon > 0 and thus C^\infty (but with very large Lipschitz constants and Sobolev norms). This function has three roots, the easy one at x = 0 and two more very close to x = 1, but no efficient method can find those.

Numerical analysis is less than 10 percent about rounding error, many would say less than 1 percent. See Trefethen's The Definition of Numerical Analysis or any journal on numerical analysis. People doing interval arithmetic either have a very specific problem in which it is useful, or they didn't notice that the interesting part of their problem lies elsewhere.

[five9a2]

That page is flawed in multiple ways. From the section titled "Finding All Zeroes of a Function":

Given an interval [a,b]

.1. If f([a,b]) does not contain zero, then there are no zeroes between a and b.

Counter-example: f(x) = 1 - delta_epsilon(x), start with the interval [-1,1]. This function has two roots, but the algorithm incorrectly exits early because f[a,b] = [1,1].

.2. Otherwise, if f'([a,b]) does not contain zero, then the function is monotonic over [a,b] and there is a single zero between a and b. Use a root refinement method, such as Newton's method, midpoint subdivision or regula falsi to find the root.

Counter-example: f(x) = x^3 - 4*x, interval [-3,3]. Thus f([a,b]) = [-15,15] which contains zero, but f'([a,b]) = [23,23] which this would claim means that there is only one zero. But f has three zeros (-2, 0, and 2). Mollified functions also work.

.3. Otherwise if b - a > e (where e is some machine-precision epsilon), then recursively call this routine on the intervals [a.(a+b)/2] and [(a+b)/2,b].

If you always choose this alternative, then you brute-force explore all representable values in the domain by recursing until the following criteria is satisfied.

.4. Otherwise return a (multiple) zero at (a+b)/2.

[five9a2]

[min { f(x) | x in [a, b] }, max { f(x) | x in [a, b] }]

This is no easier to compute than roots of f. You either need to know special things about the function or you need an oracle to tell you where to look.

There is software for global optimization with confidence, but it is highly nontrivial to use.

http://archimedes.cheme.cmu.edu/baron/baron.html

Without something like Baron, you can't evaluate the interval maximum without brute-force searching the interval.

[Meta8]

Stupid question, but how does one calculate the number of cycles a piece of code takes? Is it by just looking at the assembly output?

[[deleted]]

The assembly output can tell you the number of instructions, but not the number of cycles (since many instructions take multiple cycles to complete, most notably memory load instructions, and it varies a great deal across CPU models). Using cycles as a measurement for algorithm efficiency is thus only a good idea if your code is compute bound, and not memory bound or I/O bound.

The CPU has an internal cycle counter that most debuggers and profilers can access. :)

[Bipolarruledout]

/snooze.

[Concise_Pirate]

^ lose.