Symbolic Computation and Recursion
Vector 6:3, page 80
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 80
Background¶
I recently had occasion to indulge myself in some interesting APL, exploiting APL’s power in recursion and in manipulating character information. For the purposes of this article, I have greatly simplified the problem, and hence the solution, but even in this reduced form I think it raises some interesting issues.
The Problem¶
I was interested in the cost, C, of a certain, rather large project. C was not known precisely and was therefore regarded as a random quantity. In order to estimate it, I required its expectation E(C) and variance V(C). The following (fictitious but not entirely unrealistic) example may help to explain how such a problem may arise.
Suppose that a new section is planned for an existing motorway, connecting it to another. C is the total construction cost. In thinking about C, we will naturally decompose it into identifiable sources of cost. There will be a new junction at the connection with the other motorway, at a cost J. We will also need to replace the existing terminal of this motorway with a proper junction, but the existing work makes the cost of this junction less than J, say J × P, where P is the proportion of new work required. There will also be a bridge over the new section at a cost B. Then the cost of actually laying the road is K × L, where K is the number of kilometres involved and L is the laying cost per kilometre. (We structure this term in this way because we have good information from previous motorway projects about L.)
Finally, we could suppose that a rebate of proportion Q of the last cost is available (perhaps from the European community), so we subtract K × L × Q. We have constructed a model for C:
C = J + (J × P) + B + (K × L) - (K × L × Q)
Typically, none of the quantities J, P, B, K, L and Q will be known exactly in advance, but we will have information in the form of means and variances like E(J) and V(J). Some quantities, like Q or K, may be more or less predictable, and so have small variances, but others will be very difficult to predict and will have relatively much larger variances. We can also assume them to be essentially independent of each other. (We can define C as a cost in ‘constant pounds’ to eliminate inflation, which would otherwise cause quantities like J and B to be correlated.)
In general, then, we have a model relating C to various other random variables, which are statistically independent and whose expectations and variances are known. Then finding E(C) and V(C) reduces to applying some standard formulae for sums, differences and products of random variables. First, for any X and Y:
where C(X,Y) is the covariance between X and Y, and is zero if X and Y are independent. Also, if X and Y are independent:
(The last is not very well known but quite easy to prove.) In applying these formulae to our example model, we find only one complication. Although individual variables in the model are independent, we want to apply the formula for V(X+Y) when X and Y may be products. Now X=J and Y=J × P are not independent because they have J in common. We need another general result …
Finally, in the APL solution later, we will actually use V(X+Y) when Y is something like J+(J×P), and we will need two more formulae:
Armed with all these expressions, it is possible to find E(C) and V(C) in any formula like the example. It could be very complicated, provided it only uses addition, subtraction and multiplication. By carefully breaking it down into individual additions, subtractions and multiplications we can eventually arrive at the result. Actually, E(C) is always found simply by substituting the expectation of every variable in the model. Thus, in the example,
However, V(C) quickly becomes a very messy formula indeed.
I needed to do this not just once but often. I had different models for costs in various parts of the task I was handling, and these were all subjected to modification in the course of the work. Furthermore, the expectations and variances of variables within each formula were changing, requiring the resulting formulae for E(C) and V(C) to be applied to new values each time. There were many other complications:
- the models were typically much longer than the example above;
- all the random variables were actually random vectors, which in particular meant that variances were actually variance-covariance matrices, with perhaps over a hundred rows and columns;
- means and variances of variables in the formula were specified in various different ways.
We ignore all of these complications hereafter, but they must be mentioned to explain the form of solution that I chose.
I decided that I needed an APL system that would accept just a symbolic definition of the model as a character vector, and do all the work for me. First, it would analyse the model to produce new character vectors, one for E(C) and one for V(C), in the form of APL expressions which could simply be executed to compute values for E(C) and V(C). It was clear that analysing the model to apply the above formulae would be a highly recursive process.
The Basic Solution¶
The example model becomes the APL character vector
The rules are that the MODEL is a series of sections each beginning with a plus or minus sign, then one or more variable names. A further complication in the real problem (which changed the details of the coding quite substantially) was that variable names could be more than one character long. Notice that the model is understood in the usual algebraic sense, not as an APL expression. A minus sign, for instance, applies only to the product of variables it precedes. The rules were framed in this way because I did not want to allow parentheses. To have done so would have introduced numerous extra complications.
A function EGET was written to take as argument a model and to return a character vector with the APL expression for E(C). This is a simple recursive function, the only complication being the need to add parentheses so that when the result is executed by APL the operations are done in the right order:
∇ E←EGET M;S;Y;Z
[1] ⍝ Returns expectation of model M
[2] SYZ
[3] E←'+(',((S='-')/S),(EPROD Y),')'
[4] →(0=⍴Z)/0
[5] E←E,EGET Z
∇
∇ E←EPROD P
[1] ⍝ Returns expectation of product P
[2] →(1<⍴P)/L
[3] E←'E',P
[4] →0
[5] L:E←'E',(1↑P),'×',EPROD 1↓P
∇
For instance, the result of EGET MODEL, with MODEL defined as above, is the character vector:
There are surplus parentheses here, but it is better to have some that are not needed than to leave out any that are! The variables EA, EB, etc. will be in the workspace, containing the expectations of the variables in the model. Similarly, their variances will be in variables VA, VB, etc.
A More Complex and Interesting Recursion¶
Computing V(C) requires more complex and interesting recursion, so I have given a listing of the key functions as an appendix. I will describe their operation in a bottom-up way, since this should be clearer than the mainly top-down way in which they were originally written. The basic idea of the recursion is to split the model into three parts:
(a) the initial sign, ‘+’ or ‘-’;
(b) the sequence of one or more letters following the sign, called a product, and
(c) the rest of the model, which itself is a model.
There is a further recursion on splitting a product into its first letter and the remainder, which is another product.
Function VPROD takes as its argument a product and returns an expression for its variance, using the formula for V(X×Y). It first checks whether its argument is a single character. If so, it returns that character preceded by 'V'. Otherwise it applies the formula with X set to the first letter and Y to the rest. Since Y is in general a product, it needs to call itself for V(Y). A simple recursion results. Notice that a local variable YV is defined as the result of VPROD Y and used elsewhere in the expression generated by the last line. The variable YV is created when we execute the character result of VPROD. Because of the recursion, executing the result of applying VPROD to a product of three or more variables means YV will be redefined during execution of the result. For instance VPROD'ABD' gives the result:
The assignment of YV first to VD then to a much more complicated expression causes no problems in execution, but does not make the expression easy to read. In general, the APL expressions resulting from even quite simple models in my system are quite unreadable - but they work!
The only other thing to notice about VPROD is the call to EGET. Since the argument of EGET is a model, which must begin with a sign, a plus sign is added, then removed from the result.
Next, the function CPP applies the formula for C((X×Y),(X×Z)). Its arguments are two products. It first calls a function XYZ whose role is to identify X, Y and Z. It sets X to the vector of letters common to P1 and P2. Then Y and Z are the results of removing those letters from P1 and P2 respectively. VPP is not recursive. The only complication is the need to handle the cases when any of X, Y or Z is empty.
The next function up the system is CPM, which calculates C(X,M) when X is a product and M is a model. It first calls a function SYZ to decompose M into sign (S), product (Y) and remaining model (Z), and if the last is empty just calls CPP. Otherwise it applies the formulae for C(X,Y+Z) and C(X,Y-Z). It calls itself recursively when it needs C(X,Z). Finally, the same recursion is applied in the top level function, VGET, which calls VPROD, CPM and itself in computing V(C) by the formulae for V(X+Y) and V(X-Y).
The introduction of these few simple functions leads to high levels of recursion, even in quite short models. The expression for V(C) generated by VGET for our example model has 197 characters.
Comments¶
Symbolic computation, or ‘computer algebra’, is currently attracting considerable interest. Ordinary computation puts numerical values into a formula and produces the numerical result. Symbolic computation manipulates the formulae themselves, using specified algebraic rules, to produce new formulae. There exist an increasing number of symbolic computation packages, and this is another growing field of computing activity with which APL needs to compete. This article has given just a simple application of symbolic computation in a particular specialised context. It does, however, show that APL is quite capable of handling such problems elegantly, and with little programming effort. APL’s treatment of character data and its recursion capability are particularly useful.
For me, the exercise had many benefits. First, it saved me a lot of time. For the first few models I had to deal with, I worked out expressions for E(C) and V(C) by hand and programmed each resulting formula separately. Automating this process soon repaid the development effort. It had a second, unexpected, benefit when I was checking its results against my earlier hand-programmed models - I found two errors in my earlier work. I hate to think how many other errors I would have had to track down if I had not changed to an automatic system. The third benefit was in enjoyment. It is a lot more fun to work on something new than to do old things repetitively.
The one drawback was in speed. My hand-programmed solutions ran noticeably faster than executing the expressions generated by EGET and VGET. I made two modifications to improve execution speed, the first of which was to write a simple function to remove most of the redundant parentheses and plus signs from the expressions. The second was a compromise over allowing parentheses in the model itself. I allowed ‘local’ variables within the model to stand for submodels, provided that all ‘local’ variables remained independent of each other and of the remaining original variables. These two tricks enabled me to reduce the length of each expression considerably and so greatly increase execution speed. Symbolic computation will never produce code that runs faster than a good hand-programmed solution, but in the end I was very satisfied with my system.
Warwick University is participating in international research in symbolic computation, as well as having an ongoing interest in the use of APL in Statistics. I would be very interested to hear of anyone else who has used APL for symbolic computing, in however simple a form.
Appendix (remaining function listings)¶
∇ V←VPROD P;X;Y
[1] ⍝ Returns variance of product P
[2] →(1<⍴P)/L
[3] V←'V',P
[4] →0
[5] L:X←1↑P
[6] Y←1↓P
[7] V←'(V',X,'×',(1↓EGET '+',Y),'*2)+(YV×E',X,'*2)+V',X,'×YV←',VPROD Y
∇
∇ C←P1 CPP P2;X;Y;Z
[1] ⍝ Returns covariance of products P1 and P2
[2] XYZ
[3] C←''
[4] →(0=⍴X)/0
[5] C←VPROD X
[6] →(0=⍴Y)/L
[7] C←(1↓EGET '+',Y),'×',C
[8] L:→(0=⍴Z)/0
[9] C←(1↓EGET '+',Z),'×',C
∇
∇ C←P CPM M;S;Y;Z;L1;C1;L2;C2
[1] ⍝ Returns covariance of product P and model M
[2] SYZ
[3] →(0=⍴Z)/L
[4] L1←0<⍴C1←P CPP Y
[5] L2←0<⍴C2←P CPM Z
[6] C←((L1∧L2)/'('),(L1/S,C1),((L1∧L2)/')+'),L2/C2
[7] →0
[8] L:L1←0<⍴C1←P CPP Y
[9] C←L1/S,C1
∇
∇ V←VGET M;S;Y;Z;L1;C1
[1] ⍝ Returns variance of model M
[2] SYZ
[3] →(0=⍴Z)/L
[4] L1←0<⍴C1←Y CPM Z
[5] V←'(',(VPROD Y),')+(',(VGET Z),')',L1/'+2×',S,C1
[6] →0
[7] L:V←VPROD Y
∇
∇ SYZ;N
[1] ⍝ Splits model M into sign S, product Y and model Z
[2] S←M[1]
[3] Y←1↓M
[4] N←¯1+(Y⍳'+')⌊Y⍳'-'
[5] Z←N↓Y
[6] Y←N↑Y
∇
∇ XYZ
[1] ⍝ Identifies common terms X in P1 and P2,
[2] ⍝ setting Y and Z to the residues
[3] I←P1⍳P2
[4] X←(I≤⍴P1)/P2
[5] Y←P1[(~P1∊X)/⍳⍴P1]
[6] Z←(I>⍴P1)/P2
∇