From: andrew@rentec.com (Andrew Mullhaupt)
Newsgroups: comp.lang.apl
Subject: Re: Dead horse matrix reduction
Date: 23 Jul 92 17:23:24 GMT
Organization: Renaissance Technologies Corp., Setauket, NY.

In article <920722044324_70530.1226_CHC117-1@CompuServe.COM> 70530.1226@CompuServe.COM (Mike Kent) writes:
>OTOH, this is just math, _not_ science, and _not_ engineering (i.e.,
>it's just something you need to do *in the process of modelling something
>REALLY interesting*). 

But "really interesting things" have gone far beyond what you can do with
dense matrices, (or even matrices, in most cases). The point is not that
bandwidth reduction is the sine qua non, but that almost everywhere you
look in the computer pursuit of "really interesting things" you find
important algorithms which are bad candidates for APL.

> And there are interesting things you _can_ do
>effficiently, or efficiently enough, in APL.  

There certainly are. But as time goes by, it is a smaller and smaller
percentage of what gets done in scientific computing that fits nicely into
APL. This is a flaw in the APL language, not in science.

>See, for instance, the "Diamonds in the Sky" paper in APL 89 proceedings
>-- a stellar-interior model for white dwarf stars.  Also Scott
>Kimbrough's kinematics modeling (with graphics) in the APL 89 and APL 90
>Software Exchange packages. 

Let's take an honest survey of the literature. Find me _one percent_ of
articles in Math. Comp., Num. Math., J. Comp. Phys., (or any other similar
selection of the scientific computing literature...) which mention APL,
or were done in APL, or _even done in APL with most of the 'real' work done in
FORTRAN_, and I'll be _extremely surprised_. Find me 5% of a reasonable
survey of the literature which can be implemented effectively in APL alone
and I'll be very surprised. I'll bet you have a much better chance of finding
_Matlab_ than APL in the "really interesting" literature.

Now you will get a little better chance if you stick exclusively to the
statistics literature, which tends to be very linear compared to other
sciences. But you will find that although there are even energetic proponents
of APL in the statistics community, that the academic work is _dominated_ by
the S language (which as an interpreter is a worse language than APL) and
that the industrial stuff has big pieces in SAS, SPSS, etc.

>two largest -- Fischer's Monster, and the _really_ big one).  The grinding
>calculations were done in FORTRAN for speed (several months of full
>weekends of 308x time), but the FORTRAN was just "hand compiled" APL.

Why should APL make anyone do this?

>Or if numerical techniques (of more general interest than bandwidth
>reduction) intrigue you, check out Richard Neidinger's papers on 
>automatic differentiation, zero-finding by arbitrary-order Newton-Raphson
>techniques, expressing DEs as recurrence relations; check out the APL 89
>and APL 92 proceedings [and the citations in those papers].

Automatic differentiation is probably a decent candidate for APL. However,
root finding in APL is not such a good candidate when you take into account
that hybrid (i.e. Brent-like) methods are to be preferred over Newton type
methods in all but exceptional cases. The Brent methods are control bound,
(and I've implemented them personally in APL), and are bad candidates for
APL.

>Anybody who's spent much time doing serious computational work with APL
>knows that there are classes of problems for which APL is not well suited
>-- anything which invloves a lot of data which has to be accessed an item
>at a time (which sparse represenations of matrices, for instance, force
>upon you).  (That's why IBM invented []NA, and IPSA invented AP (RBE, pls
>fill in AP #), and STSC invented []CALL and []XP, etc.)  If you use APL
>for what it's good for, and use one of these interfaces to hand off the
>other 10% of the problem to some FORTRAN or C or assembler code, hopefully
>from a package like LINPACK or ESSL rather tahn purpose-built, you are
>using your resources intelligently:  the human isn't trying to get several
>thousand lines of scalar-oriented code working (and working right), and
>the machine isn't encountering interpretive and space-management overhead
>while doing scalar-oriented computations.

The claim here is that this isn't good enough. Most of these interfaces
force copying - which means that you can't use them for small, control
bound code. You must instead implement the whole of such algorithms in
the alien language. My point is that more and more scientific algorithms
are of this kind, and that there is a reason for this. Maybe I should
explain:

To every instance of an algorithm, let us express the cost of computing the
solution in terms of the size N of the instance; i.e. the cost is C(N).
Since time vs. space can be reasonably traded off in most problems, we
won't worry about distinguishing them - they're both expressed in the
function C(N), having said this, we'll usually mean _time_ rather than
space, to avoid talking about garbage collection while waiting for input,
etc.

We might have (as is very common) that the instances represent a family of
approximations to an exact solution to a non-finitary problem, and that
the accuracy of the approximation increases as N tends to infinity. Another
case of broad interest is when the intances naturally represent some aspect
of the size of the problem specification - e.g. how many galaxies, etc. In
this case, it is more or less experience which tells us that N will go to
infinity, more than anything else.

Now at any given time, you have a fixed amount M of computer resource, and
so you compare feasability of computations by asking if C(N) < M. If more 
than one algorithm is feasable for a problem, say C_1(N) < M and C_2(N) < M,
then you choose on the basis of whether C_1(N) < C_2(N) or not.

Scientific computing doesn't tend to involve many cases where C(N) is
not bounded by some polynomial of N. Let's just assume this. We can argue
it later...

Now it is often the case that two algorithms for the same problem have
costs which satisfy C_1(N) < C_2(N) for N < N* and C_1(N) > C_2(N) for
N >= N*. As time goes on, and your resources increase, it is clear that
C_2(N) will be the chosen algorithm. For example, traditional matrix
multiplication has cost C_1(N) = N^3, and Strassen multiplication has
cost cN^2.81, where 2.81 = log(7)/log(2). This means that no matter what
c is, eventually, we will choose Strassen multiplication over the traditional
approach. Keep in mind that it is only interesting if c > 1, else we would
always choose Strassen multiplication.

Some will argue that the complexity of writing the Strassen algorithm will
prevent us from using it. Yes, but only for a little while. If the problem
is important enough, we will find that since the human time it takes to
write and debug any fixed algorithm is O(1), that this will be overcome
by a scaling advantage, too.

So it is clear. We will choose algorithms for their asymptotic _scaling_
properties, at least eventually. But what does this say about languages?
It turns out that most languages do not get in the way of the _exponent_
in the cost function - if you write an FFT in APL it will still be O(N log N).
The language is principally reflected in the value of c that obtains. If
a FORTRAN FFT and an APL FFT are compared, then the APL value of c is likely
to be much bigger.

FORTRAN, C, and many other languages are often quite close to 'full-bandwidth'
by which I mean that the value of c resulting from an implementation in these
languages is generally within a factor of, say 3, of what the machine can
do. APL is quite a different story. If what you want to do in APL is given
exactly by some convenient string of primitives, and your problem is large
(which we assume here) then the value of c will be very close to the best
possible for the machine. The problem is that some operations in APL are
very far from full-bandwidth.

Now almost all the primitives in APL have 'exact' polynomial complexity -
the exceptions being things like sorting and searching. What this tells us
is that any APL implementation of say, an algorithm with complexity O(N^1.5)
will have a _bad_ constant. (Try the well known O(N^1.5) case of Shell's
sort, and compare it to the upgrade primitive just as an example.)

In fact this is the reason that sorting and searching _should_ have primitives.
but the conclusion of the argument is already obvious: Unless we have good
c constants for operations in which we can implement algorithms of _every
order of complexity_, we will eventually run into problems ill-suited to
the choice of primitives. So sticking in a couple (like diagonals, etc.) is
only a band-aid, because as time goes on, we will want to choose an algorithm
with a complexity for which we will be guaranteed a bad constant.

The only set of primitives which can possibly guarantee fast algorithms of
every order of complexity include fast operations on single data items. I'm
pretty sure that this makes life tough on interpreters in general, but APL
in particular.

I've always been amused at APL programmer's notion that 'One at a time'
programming is always bad. It leads them into what I call 'Too many at a time'
programming. My argument here is that 'one at a time' programming is destined
to be the only efficient programming.



>Oh, and BTW if you want an APL compiler you might contact STSC --
>they sell one which works within their mainframe product -- or
>Interprocess Systems -- they have the Yorktown APL Translator (APL to
>FORTRAN, vector FORTRAN if you like) which IBM Research did originally,
>only now as a nicely packaged product [IS is an IBM Business Partner].
>Or if your code has to run on a Cray, drop a note to Bob Bernecky and ask
>about ACORN -- APL to Cray C.


APl compilation is essential in my way of thinking to saving APL. But will
it _really_ save APL? The crux of the argument is that as time goes by,
the existing primitives will become less and less frequently the foundation
of good algorithms.

I think it is a very good question whether there are any large scale primitives
which do not 'fade away' in this model. It would be very nice to think that
someday we could do better than FORTRAN. I'm sure that we can, but the only
permanent hope is large scale operations which can effectively generate a
very wide range of complexities.

Later,
Andrew Mullhaupt