Computation: Cleaner, Faster, Seamless, and Higher Quality Alan Edelman Massachusetts Institute of Technology

1
Computation Cleaner, Faster, Seamless, and
Higher QualityAlan EdelmanMassachusetts
Institute of Technology
2
Unique Opportunity

• Raise the bar for all functions
• Linear Algebra!
• Elementary Functions
• Special Functions
• Hard Functions Optimization/Diff Eqs
• Combinatorial and Set Functions

3
Unique Opportunity

• Raise the bar for all functions
• Linear Algebra!
• Elementary Functions
• Special Functions
• Hard Functions Optimization/Diff Eqs
• Combinatorial and Set Functions
• 0-d, 1-d, n-d arrays, empty arrays
• Inf, Not Available (NaN)
• Types (integer, double, complex) and embeddings
• (real inside complex, matrices insidse n-d
  arrays)
• Methodology
• Work in Progress Create a web experience that
  guides and educates users in the amount of time
  that it takes to say Wikipedia
• Futuristic? The Semantic Web and Other
  Automations?

4
What the world understands about floating point

• Typical user In grade school, math was either
  right or wrong. Getting it wrong meant bad
  boy/bad girl. (Very shameful!)
• By extension, if a computer gets it wrong Shame,
  Bad computer! Bad programmer!

5
What the world understands about floating point

• Typical user In grade school, math was either
  right or wrong. Getting it wrong meant bad
  boy/bad girl. (Very shameful!)
• By extension, if a computer gets it wrong Shame,
  Bad computer! Bad programmer!
• People can understand a bad digit or two in the
  insignficant part at the end. We might say
  small forward error, but they think something
  more like fuzzy measurement error.

6
The Linear Algebra World

• Some of the worlds best practices.
• Still Inconsistent, still complicated.
• Never go wrong reading a book authored by someone
  named Nick, a paper or key example by Velvel, all
  the good stuff in terms of error bounds of LAPACK
  and papers from our community

7
Landscape of Environments

• Traditional Languages plus libraries
• (Fortran, C, C, Java?)
• Proprietary Environments
• MATLAB, Mathematica, Maple
• Freely Available Open Source
• Scilab, Python(SciPy)(Enthought!), R, Octave
• Other Major Players
• Excel, Labview, PVWave

8
Computing Environments whats appropriate?

• Good numerics?
• Who defines good?
• How do we evaluate?
• Seamless portability?
• Rosetta stone?
• Consumer Reports?
• Easy to learn?
• Good documentation?
• Proprietary or Free?
• Download, next, next, next, enter, and compute?

9
Numerical Errors Whats appropriate?

• The exact answer?

10
Numerical Errors Whats appropriate?

• The exact answer?
• What MATLAB gives!

11
Numerical Errors Whats appropriate?

• The exact answer?
• What MATLAB gives!
• float(f(x))??? (e.g. some trig libraries)
• Small Forward Error
• Small Backward Error
• Nice Mathematical Properties
• What people expect, right or wrong
• Consistent?
• Nice mathematics?

12
Simple Ill Conditioning

• Matrix gtgt amagic(4)inv(a)
• Warning Matrix is close to singular or
  badly scaled.
• Results may be inaccurate. RCOND
  1.306145e-17.
• Statistics Function erfinv
• erfinv(1-eps)
• Exact 5.8050186831934533001812
• MATLAB 5.805365688510648
• erfinv(1-eps-eps/2) is exactly 5.77049185355333287
  054782 anderfinv(1-epseps/2) is exactly
  5.86358474875516792720766
• ?condition number is O(1e14)

13
Some Conclusions?

• MATLAB should have had a small forward error and
  its been bad?
• Nobody would ever ask for erfinv at such a
  strange argument?
• MATLAB should have been consistent and warned
  us, like it has so nicely done with inv which
  raises the bar high, but has let us down on all
  the other functions?
• We have a good backward error so we need not
  worry?
• Users either are in extremely bad shape or
  extremely good shape anyway
• Have such huge problems due to ill-conditioning
  they should have formulated the problem
  completely differently (as in they should have
  asked for the pseudospectra not the spectra!)
• The ill-conditioning is not a problem because
  later they apply the inverse function or project
  onto a subspace where it doesnt matter

14
Discontinuity

• octavegt atan(1i264(10 eps))ans
•    1.5708  -1.5708

15
Discontinuity

• octavegt atan(1i264(10 eps))ans
•    1.5708  -1.5708
• Exact? tan(xpi) is tan(x) after all so is it
  okay to return atan(x) pi any random integer?
• There is an official atan with definition on the
  branch cuts i(1,8)
• Continuous everywhere else!
• MATLAB atan(i2509) is pi/2 and atan(i2510)
  is pi/2
• atan(1i2509) is pi/2 as is
  atan(1i2510)

16
QR?

• No different from atan
• Good enough to take atanpi(random integer?)
• Is it okay to take Qdiag(random(signs?))?
• There is an argument to avoid gratuitous
  discontinuity! Just like atan
• Many people take Q,Rqr(randn(n)) to get
  uniformly distributed orthogonal matrices! But
  broken because Q is only piecewise continuous.
• Backward error analysis has a subtlety. We
  guarantee that Q,R is a QR decomposition of a
  nearby matrix, not the QR decomposition of a
  nearby matrix. Even without roundoff, there may
  be no algorithm that gives THIS QR for the
  rounded output.
• Is it okay to have any QR for rank deficiency?

17
Snap to grid?

• Sin( float(pi) n)
• 0??? Vs. float(sin(float(pi)n))
• Determinant of a matrix of integers?
• Should it be an integer?
• It is in MATLAB! Cleve said so on line

18
But is it right?

• asign(randn(27)),det(a)
• a
• 1 1 1 1 -1 -1 -1 1
  1 -1 1 -1 -1 1 -1 1
  -1 1 -1 -1 -1 -1 -1 1
  -1 -1 -1
• 1 -1 -1 1 -1 1 -1 -1
  -1 -1 -1 1 -1 -1 -1 -1
  -1 1 -1 -1 -1 -1 1 -1
  1 1 1
• -1 1 -1 -1 1 1 1 1
  1 -1 1 -1 1 1 1 1
  -1 -1 1 -1 -1 1 1 -1
  1 1 -1
• -1 -1 1 -1 -1 -1 -1 1
  -1 -1 -1 -1 1 -1 1 -1
  -1 1 1 -1 -1 1 -1 -1
  -1 -1 -1
• -1 1 1 -1 -1 -1 -1 -1
  1 -1 -1 1 1 1 -1 -1
  -1 -1 -1 1 1 1 -1 1
  1 1 -1
• 1 1 1 1 1 1 -1 1
  -1 -1 1 1 -1 1 -1 -1
  1 1 -1 -1 1 1 1 1
  -1 1 -1
• -1 -1 -1 -1 -1 -1 1 -1
  1 -1 1 -1 1 -1 1 1
  1 1 1 -1 1 -1 1 -1
  -1 1 -1
• 1 1 1 -1 1 1 -1 -1
  1 1 -1 -1 1 1 -1 1
  -1 1 1 1 1 1 -1 -1
  1 1 -1
• -1 1 -1 1 1 1 -1 -1
  1 -1 1 1 1 1 -1 1
  1 1 -1 1 1 -1 1 1
  1 -1 -1
• -1 1 -1 1 1 -1 -1 -1
  -1 1 1 1 1 1 -1 1
  -1 1 1 -1 -1 1 1 1
  1 -1 -1
• -1 -1 -1 -1 1 1 -1 -1
  -1 1 1 -1 -1 1 -1 -1
  1 1 1 1 -1 1 1 1
  -1 -1 1
• 1 1 -1 -1 1 -1 -1 -1
  1 -1 -1 -1 -1 -1 -1 1
  -1 1 -1 1 1 -1 1 1
  1 -1 -1
• -1 1 -1 -1 -1 -1 1 -1
  -1 -1 -1 1 -1 -1 1 1
  -1 1 -1 -1 1 1 -1 -1
  1 -1 -1
• 1 1 -1 1 1 -1 1 1
  -1 -1 -1 1 1 1 1 -1
  -1 1 -1 1 1 -1 1 -1
  -1 1 -1
• -1 1 -1 1 -1 1 1 1
  -1 -1 1 1 -1 -1 -1 -1
  -1 1 1 1 1 1 -1 1
  1 -1 -1
• -1 -1 1 -1 1 -1 1 -1
  -1 1 -1 -1 -1 -1 -1 -1
  1 -1 1 -1 1 1 -1 1
  -1 -1 1
• 1 1 -1 1 1 -1 1 -1
  -1 1 1 -1 -1 1 -1 -1
  1 1 -1 1 -1 1 1 -1
  1 1 -1

19
MPI Based Libraries

• Typical sentence we enjoy using parallel
  computing libraries such as Scalapack
• What else? you know, such as scalapack
• And ? Well, there is scalapack
• (petsc, superlu, mumps, trilinos, )
• Very few users, still many bugs, immature
• Highly Optimized Libraries? Yes and No

20
Star-P
21
Some of the hardest parallel computing
problems(not what you think)

• A row distributed array (or worse)
• B column distributed array(or worse)
• CAB
• Indexing
• CA(I,J)

22
LU Example
23
Scientific Software

• Proprietary
• MATLAB, Mathematica, Maple
• Open Source
• Octave and Scilab
• Python and R
• Also
• IDL, PVWave, Labview
• What should I use and whats the difference
  anyway?
• Seamless how to avoid the tendency to get stuck
  with one package and believe its the best
• Seamless portability,
• Higher Quality standards, accuracy

24
Documentation I would like

• Tell me on page 1 how to type in a 2x2 matrix
  such as

1 2
1i nan
25
function functobj_func(v,T) sigalphav(1) T0v
(2) betav(3) fv(4) R0.008314 aNsigalphaf/
2 aPsigalpha(2-f)./2 H-100021000' Hr
Hcsize(H) G0zeros(Hr,1) for i1Hr if
H(i)lt0 G0(i)-2beta(H(i)./aN).2abs(bet
a)(H(i)./aN).4 else
G0(i)-2beta(H(i)./aP).2abs(beta)(H(i)./aP).
4 end end sthexp(-G0./(RT0)) PT0sth./(tr
apz(sth(H(2)-H(1)))) Tr Tcsize(T) for
i1Tr sth2PT0.exp(-(1/R)((1/T(i))-(1/T0)).
H) PT(,i)sth2./(trapz(sth2(H(2)-H(1))))
Hi(i)(1./(trapz(sth2(H(2)-H(1)))))trapz(H.
sth2(H(2)-H(1))) H2(i)(1./(trapz(sth2(H(2)
-H(1)))))trapz((H.2).sth2(H(2)-H(1)))
Cpex(i)(H2(i)-Hi(i).2)./(R.(T(i)2)) end funct
Cpex'
function functobj_func(v,T) betav(3)fv(4)R
0.008314H-100021000' HH 2 
-10002-2'/f021000'/(2-f)/v(1) G0-2bet
aHH.2 abs(beta)HH.4 sthexp(-G0/(Rv(2)))
PT0sth/sum(sth)  EEexp((1/v(2)-1./T)H'/R)tt
EEPT0Hi(EE(H.PT0))./ttH2(EE(H.2.PT0)).
/ttfunct(H2-Hi.2)./(RT.2)
26
A Factor of 10 or so in MATLAB

• Life is too short to write for loops ?
• Trapz is just a sum if the boundaries are 0
• H(2)-H(1) never needed (cancels out)
• Sums are sometimes dot products
• Dot products are sometimes matvecs

27
Higher Quality
28
New Standards for Quality of Computation

• Associative Law
• (ab)ca(bc)
• Not true in roundoff
• Mostly didnt matter in serial
• Parallel computation reorganizes computation
• Lawyers get very upset!

29
Accuracy Sine Function Extreme Argument
sin(264)0.02359850990443955863436
59 One IEEE double ulp higher
sin(2644096)0.6134493700154282634382759 over
651 wavelengths higher. 264 is a good test case
Maple 10 ?evalhf(sin(264))
0.0235985099044395581 note
evalf(sin(264)) does not accurately use 264 but
rather computes evalf(sin(evalf(264,10)))
Mathematica
6.0 ?? PrintNSin264 PrintNSin264,2
0 0.312821 0.0235985099044395
58634 MATLAB ? sin(264)
0.023598509904440 sin(single(264))
0.0235985 Octave 2.1.72 ? sin(264)
0.312821314503348 python 2.5 numpy ?
sin(264) 0.24726064630941769 R 2.4.0
? format(sin(264),digits17)
0.24726064630941769 Scilab
? format(25)sin(264) 0.023598509904439558121
4
Mathematica and python have been observed to
give the correct answer on certain 64 bit linux
machines. The python and R above numbers
listed here were taken with Windows XP. The .3128
number was seen with python on a 32 bit linux
machine. The correct answer was seen with R on a
sun4 machine. Its very likely that default N in
MMA, Octave, Python, and R pass thru to various
libraries, such as LIBM, and MSVCRT.DLL.
Notes Maple and Mathematica use both hardware
floating point and custom vpa MATLAB discloses
use of FDLIBM. MATLAB and extra precision
MAPLE and Mathematica can claim small forward
error All others can claim small backward error
Excel doesnt compute 264 accurately and sin
is an error. Google gives 0.312821315 Ruby on
http//tryruby.hobix.com/ gives -0.35464 for
Math.sin(264). On a sun4 the correct was given
with Ruby. Relative Backward Error lt (2pi/264)
3e-19 Relative Condition Number
abs(cos(264)dx/sin(264))/(dx/264) 8e20 One
can see accuracy drainage going back to around
221 or so, by powers of 2
30
Coverage Sine Function Complex Support

• All packages seem good
• Maple evalf(sin(I))
• Mathematica NSinI
• MATLAB sin(i)
• Octave
• Numpy sin(1j)
• R sin(1i)
• Ruby require 'complex' ComplexI
  Mathsin(ComplexI)
• Scilab sin(i)
• EXCEL IMSIN(COMPLEX(0,1))
• Note special command in excel

31
Coverage Erf Function Complex Support

• python/scipy is the only numerical package good
• Maple evalf(erf(I)) 1.650425759 I
• Mathematica NErfI 0. 1.65043 I
• MATLAB erf(i) ??? Input must be real.
• Octave error erf unable to handle
  complex arguments
• Python/Scipy erf(1j) 1.6504257588j
• R pnorm(1i) Error in pnorm Non-numeric
  argument
• Scilab erf(i) !--error 52 argument must be a
  real

32
Acceptable or not?

• tan(zkp)tan(z)
• atan officially discontinuous on branch cut
• atan discontinuous off branch cut?
• MATLAB does it
• MATLAB highly inaccurate in some zones
• Python, octave, R worse in some ways
• Scilab very Good
• QR does it? (unacceptable in my opinion)

33
Clones

• Just because Octave is a MATLAB wannabe (want
  to be) doesnt mean they are the same
• Better some ways, worse in other ways

34
Type Casting

• Mathematica and Maple, type cast as the symbolic
  languages have excellent numerics
• MATLAB the matrix laboratory uses the same
  lapack and math libraries as everyone else
• MATLAB has 1500 functions some terrific some
  trivial

35
Optimization

• Would love to see the Consumer Reports Style
  Article