amull@Morgan.COM (Andrew P. Mullhaupt) (11/27/89)
I would like to solicit suggestions for primitives to be added to the APL lexicon which would update the representation of linear algebraic concepts. In particular I have in mind a way to deal with the concept of matrix factorizations; Cholesky, QR, and SVD. The point is that if you use quad-divide, you have to factor the matrix every time you come up with a right hand side (left hand side in APL terms) that needs evaluating. Having access to the factored form of the matrix, especially in a storage efficient form is quite desirable. In particular, since I am using SVD to implement quad-divide in a version of APL, I think it's a stupid shame to have to junk the matrix factors. It is also desirable to have the choice of more than one decomposition available to the programmer, and in the case of the LU; why should double the storage be required? At present, the implementation is going to bail out of APL for almost all of this stuff. This will have the same result that it had for us on APL2 on the mainframe: there, quad divide is not competitive with the ESSL SVD, so we never use it for much anymore. (In fact I checked around our group of programmers, and I think it is never used in new code.) Now HP calculators keep track of whether or not a given matrix has been factored in place or not, but this doesn't give much insight into what to do with the SVD. Thanks in advance, Andrew Mullhaupt P.S. Sparse and band matrix tricks, and updates (Cholesky and QR) would seem very desirable if they can be embedded comfortably in the APL syntax.
jeff@zem.uchicago.edu (Jeff Adler) (11/29/89)
In article <539@s5.Morgan.COM> amull@Morgan.COM (Andrew P. Mullhaupt) writes: >I would like to solicit suggestions for primitives to be added to >the APL lexicon which would update the representation of linear >algebraic concepts... > >At present, the implementation is going to bail out of APL for >almost all of this stuff. This will have the same result that >it had for us on APL2 on the mainframe: there, quad divide is >not competitive with the ESSL SVD... > >P.S. Sparse and band matrix tricks, and updates (Cholesky and >QR) would seem very desirable if they can be embedded comfortably >in the APL syntax. Rather than trying to get APL to do everything, wouldn't it be better to have APL 1. do a small number of things really well, 2. have hooks into other programs with other capabilities? Isn't this the point of quadNA? If you need some fancy mathematical functions, what's wrong with calling standard routines from a library? In general, the issue of primitive vs defined vs external routines applies not only to mathematical functions, but to editors, graphics, communication, and who-knows-what-else. My view: If another product does something well, and if APL can access that product through a reasonable hook, then the hook is preferable to adding a feature to APL. Does anyone else have any opinions on this? --------------------------------------------------------------------------- Jeff Adler jeff@zaphod.uchicago.edu Math Department jeff%zaphod@gargoyle.uchicago.edu University of Chicago jda3@tank.uchicago.edu
amull@Morgan.COM (Andrew P. Mullhaupt) (12/01/89)
In article <6458@tank.uchicago.edu>, jeff@zem.uchicago.edu (Jeff Adler) writes: > In article <539@s5.Morgan.COM> amull@Morgan.COM (Andrew P. Mullhaupt) writes: > >I would like to solicit suggestions for primitives to be added to > >the APL lexicon which would update the representation of linear > >algebraic concepts... > > > >in the APL syntax. > > Rather than trying to get APL to do everything, wouldn't it be better > to have APL > 1. do a small number of things really well, > 2. have hooks into other programs with other capabilities? > > Isn't this the point of quadNA? If you need some fancy mathematical > functions, what's wrong with calling standard routines from a library? Well for one thing, most of those wonderful FORTRAN oriented libraries require spurious transposes which are another thing that APL doesn't exactly do instantaneously. Efficient use of stride-oriented routines from APL requires some intriguing juggling and often extra data movement and copying. And why must this stuff get re-examined when you port your APL to a machine where the addressing and calling conventions are not the same? You try writing 13 versions of quadNA. Not me. Then there is the issue of quad divide itself. The days of dense matrices are long over, and there is a great opportunity for APL. The concept of the matrix divide primitive is simple and clean, but it is at odds with the state of the art in numerical linear algebra. You'd really like to split the call across the factorization and the solve because that is what you really want mathematically as well as for the sake of speed. It seems a shame not to add to APL primitives which could radically improve on it's scope and give it an edge on all other languages for scientific computing, right where it has always excelled: matrix operations. > > In general, the issue of primitive vs defined vs external routines > applies not only to mathematical functions, but to editors, graphics, > communication, and who-knows-what-else. > > My view: If another product does something well, and if APL can > access that product through a reasonable hook, then the hook is > preferable to adding a feature to APL. > Can't always be done. Frankly, APL is terrible at control-bound code. Now one very well-understood case of this is minimization of smooth functions of many variables, (for example the BFGS Quasi-Newton method with finite difference approximate derivatives). Now you really can't write a nice version in APL: I'm at a place where I wrote one myself as a challenge to several of our fairly expert programmers. (I figured they'd show me how to do it better). Well, they use my code because it's too disgusting to write it in APL. Many good versions of this algorithm exist in FORTRAN, so why not use them? One good reason: The algorithm needs to call the user-supplied APL function, and as everybody knows, calling APL from FORTRAN is uh - as IBM puts it "impossible". So I leave this hugely important example as a lasting testimony to why you either need to tie down variables in the workspace (an APL philosophical "no-no") or else let me put stuff I need (structured constructs for low overhead control flow) in APL. (I think this is an APL philosophy "lose-lose" situation.) Your move. Also: One really cool primitive I'd like to see would be something like this; Let M be a matrix. The diagonal primitive returns the vector of diagonal elements of M (easy in APL2, a small modification works in ISO). For example: if M is 3 3 rho iota 9 then the diagonal primitive returns 7 4 8 1 5 9 2 6 3 . You can write this function in APL but it'll be too slow, and you won't be able to do the best part: Make the dyadic version so the K diagonal M is the Kth diagonal. Then in APL2 you'd expect to be able to do the selective specification (K diagonal M) gets V and have the diagonal assigned IN PLACE and RAPIDLY. The interpreter can keep a matrix constructed this way as a band matrix until it comes to something which causes it to densify the matrix. That is, if you start with a zero matrix and put diagonals in it, why take up more space than you need? Then you can make very long band matrices and allow an equipped APL interpreter to speed up and save space, and allow a less complex interpreter to implement it in dense form. It would be pretty useful to let plus dot times and quad divide in on the joke, too, else you don't gain much. Later, Andrew Mullhaupt
jeff@trillian.uchicago.edu (Jeff Adler) (12/02/89)
In article <558@s5.Morgan.COM> amull@Morgan.COM (Andrew P. Mullhaupt) writes: >In article <6458@tank.uchicago.edu>, jeff@zem.uchicago.edu (Jeff Adler) writes: >> In article <539@s5.Morgan.COM> amull@Morgan.COM (Andrew P. Mullhaupt) writes: >> >I would like to solicit suggestions for primitives to be added to >> >the APL lexicon which would update the representation of linear >> >algebraic concepts... >> > >> >in the APL syntax. >> >> Rather than trying to get APL to do everything, wouldn't it be better >> to have APL >> 1. do a small number of things really well, >> 2. have hooks into other programs with other capabilities? >> >> Isn't this the point of quadNA? > >Well for one thing, most of those wonderful FORTRAN oriented libraries >require spurious transposes [...] >And why must this stuff get re-examined when you >port your APL to a machine where the addressing and calling conventions >are not the same? You try writing 13 versions of quadNA. Not me. Are you porting APL code or the APL interpreter? If the former, then the vendor should handle this. (I realize that currently only one implementation of APL has quadNA. But if you're going to ask other implementors to add more matrix functions, why not ask them to add quadNA?) If the latter, see my proposal further below. > [...] It seems a shame not to add to APL >primitives which could radically improve on it's scope and give >it an edge on all other languages for scientific computing, right >where it has always excelled: matrix operations. > >> My view: If another product does something well, and if APL can >> access that product through a reasonable hook, then the hook is >> preferable to adding a feature to APL. > >Can't always be done. Frankly, APL is terrible at control-bound code. >Now one very well-understood case of this is [...] >Many good versions of this algorithm exist in FORTRAN, so why not >use them? One good reason: The algorithm needs to call the user-supplied >APL function, and as everybody knows, calling APL from FORTRAN is >uh - as IBM puts it "impossible". So I leave this hugely important >example as a lasting testimony to why you either need to tie down >variables in the workspace (an APL philosophical "no-no") or else > [...] > >Your move. Ok, I see your general point. I don't think your example contradicts my point, though. Rather, you've shown that reasonable hooks don't always exist. But maybe better ones can be found. More on that below. In the meantime, you propose adding these functions as primitives. But if we added every useful matrix function that anyone considers basic, we'd have a proliferation of new primitive functions, all of which would be useless to those not doing matrix algebra. (BTW: doesn't it seem a little bit strange that every ISO-conforming interpreter, no matter how dinky the machine it runs on, must provide the gamma and matrix inversion functions?) Instead of primitives, we could add "magic" functions (or operators, in your case), like IBM's POLYZ. Objection #1: "We" can't do that. The vendor has to do that. Answer: True. Just like a primitive. Objection #2: Unlike primitives, magic functions aren't available unless you )COPY them into your workspace. Answer: C has standard libraries of useful things which are not part of the language. To use them, you must explicitly "include" them. Maybe APL should have a standard library. Objection #3: If you are porting to Brand X APL, which doesn't support the standard library, then you're sunk. Answer: If Brand X doesn't support all the matrix primitives, then you're in the same situation. Objection #4: Magic functions are a form of cheating. Answer: I agree. Finally, here's an idea for a better 'hook'. Actually it's not really a hook, since the functions you call have to be specially written to be called from APL. Also, it depends on the interpreter being constructed in a certain way. But the facility described below allows one to write one's own magic functions, and allows the vendor to provide source code for the standard library, making it portable. Suppose you had a facility which was a cross between APL2's quadNA and Dyalog APL's auxilliary processors. The latter are programs (typically written in C) which set up associated functions in the workspace, and which communicate with the interpreter through those functions, not through shared variables. An AP always receives APL data in the form of a C structure (of a predefined type, called ARRAY), which contains information about shape, rank, type, etc. I believe that the Dyalog interpreter itself represents arrays in this way. Thus, no data conversion takes place between the interpreter and the AP. Also, the AP is just as portable as the interpreter itself. If the vendor makes slight changes to the internal representation of data, then it is only necessary to change the definition of ARRAY, and APs will still compile correctly. If the vendor provides a new APL data type (sparse matrix, maybe?) a properly written AP should be able to handle it correctly (though not necessarily efficiently) after being recompiled. quadNA currently allows one to call functions which have access to the internal representation of their arguments, but the functions have to be written in assembler. What I've described above is similar, except that one can use a higher-level language. You also want to be able to pass functions as operands. This sounds like a feasible extension of the Dyalog scheme. If we are to call functions from APL2 in this way, then IBM needs to supply another 'associated processor.' Then we can use quadNA. Objections, arguments, and criticism welcome. >Also: One really cool primitive I'd like to see would be something >like this: [description of "diagonal" primitive]. I agree. --------------------------------------------------------------------------- Jeff Adler jeff@zaphod.uchicago.edu Math Department jeff%zaphod@gargoyle.uchicago.edu University of Chicago jda3@tank.uchicago.edu
amull@Morgan.COM (Andrew P. Mullhaupt) (12/03/89)
In response to your response (which for some reason won't include in this response (?!)). We are writing a protable APL interpreter. That's why I want to know if there are good matrix primitives for modern numerical linear algebra. I think the only hope to get over the sluggishness of 'one at a time' operations in APL is to devise elegant primitives for the new language. The AP approach seems better at first, but then you get the 'quality of implementation' and 'differential portability' problems. Since we are the originators of our own APL (with most of ISO plus some Dictionary and even a few other things) we have to ask the question of what should be in it? Iverson's first book on APL included data types and primitives not yet included in APL, so it seems clear that rethinking the primitives has been done before. My proposed diagonal primitive turns out to be close to something proposed at SIGAPL 1981. We must keep in mind that APL became a language before the advent of widespread acceptance of orthogonal matrix factorizations for solving least squares problems, so in some sense quad divide is not far off the state of the art circa 1965. Why should APL remain interested in a primitive which can be almost always dominated in usefulness and efficiency by a different set? Any new primitives would completely embrace the current functionality of quad divide and machine updating existing code would be nearly trivial, so what's the big deal? Later, Andrew Mullhaupt
hafer@tubsibr.uucp (Udo Hafermann) (12/05/89)
How about discussing a "standard" way of having user-defined functions written in a "foreign" language, say, C? (I myself found C to be a fairly convenient environment when implementing auxiliary processors for APL.68000, as a good deal of interface details can be "hidden" in macros or functions, and as C is rather good at accessing "foreign" data structures.) Or do the internal representations of APL arrays vary too much across implementations to allow even thinking of this kind of thing?
raulmill@usc.edu (Raul Deluth Rockwell) (12/06/89)
In article <539@s5.Morgan.COM> amull@Morgan.COM (Andrew P. Mullhaupt) writes:
I would like to solicit suggestions for primitives to be added to
the APL lexicon which would update the representation of linear
algebraic concepts. In particular I have in mind a way to deal
with the concept of matrix factorizations; Cholesky, QR, and
SVD. The point is that if you use quad-divide, you have to factor
the matrix every time you come up with a right hand side (left
hand side in APL terms) that needs evaluating.
Seems to me that factorization tends to yield multiple results. This
could be dealt with by 'inventing' a new type, but how much would be
saved by this approach?
Comparing factorization with Guass-Jordan reduction, how many
operations are saved over the course of an algorithm? (Assume some
sort of ideal efficiency in the implementation of factorization). I
suppose I ought to do some of this myself, but I don't know what kind
of algorithms you are trying to optimize.
--amull@Morgan.COM (Andrew P. Mullhaupt) (12/09/89)
In article <RAULMILL.89Dec5214208@aludra.usc.edu>, raulmill@usc.edu (Raul Deluth Rockwell) writes: > > Seems to me that factorization tends to yield multiple results. This > could be dealt with by 'inventing' a new type, but how much would be >saved by this approach? If you want to escape from the incredibly painful limitation of APL to dense matrices, quite a bit. Also; there are different factorizations to fit different needs. You may find occasion to use Cholesky's method and at other times need for the SVD. That there is a pretty big factor. > > Comparing factorization with Guass-Jordan reduction, how many > operations are saved over the course of an algorithm? (Assume some > sort of ideal efficiency in the implementation of factorization). I > suppose I ought to do some of this myself, but I don't know what kind > of algorithms you are trying to optimize. Gauss-Jordan is not normally acceptable as such for a wide-spectrum quad-divide primitive. APL2, for example specifies that Householder QR (Hanson-Lawson method) is used. You actually lose operations to Gauss-Jordan on each call, and that's one of the motivations to save the factors. Matrix factorizations address a wide spectrum of numerical linear algebra problems, among them we are interested in the QR, SVD and Shur decompositions, and symmetric, polydiagonal, and sparse implementations. Good algorithms for these are well known, but not easily translated into time and space efficient APL. What might be the biggest saving is not having to worry about the same stuff each time you have to get some of the information from a column pivoted QR, for example. Later, Andrew Mullhaupt