Parsing Matrix Expressions

A parser for evaluating matrix expressions

By Bill Murray, Annandale, Virginia

About the author

Bill Murray is retired from NASA, having worked at the Goddard Space Flight Center in Greenbelt, MD for 22 years as an applied mathematician. In the early days of the Apollo program he was involved with orbital calculations and statistical studies. Later he worked in the area of mathematical modelling of satellite imagery data from passive microwave radiometers. He has a BA in Mathematics from Duke University (1954) and a MS in Applied Mathematics from Catholic University (1962). His main “outside” activity over the past 18 years has been long distance running, (presently running shorter distances), swimming, and hiking. He also enjoys playing some of the old favorites on the piano for the elderly.

Introduction

This article describes a parser for evaluating matrix expressions. A parser for calculating real numbers has been described in [1]. The present discussion extends the ideas introduced in [1], using a matrix as a token. Since a real number is defined as a matrix with one row and column, real number calculations can be made with the present algorithm. Similar to [1], a Pascal function, eval, is calculated which takes as its only input a str255 variable containing the matrix expression, and returns the value of a real number, a pointer to a matrix, or an error message. In addition, the Pascal code for eval can be incorporated in software code to make an existing program more flexible.

Equally as important as developing the tools for implementing a matrix parser is finding a way to efficiently allocate space and store a matrix array. A good bit of the following discussion, therefore, will be devoted to arrays, pointers, handles, and arrays of handles.

In order to conserve space on the heap, a handle is used to point to a single (as opposed to double) subscripted array of reals dynamically allocated with the NewHandle function. Using NewHandle and long integers for array indices, the maximum number of elements in an array in the type declaration can be made quite large, and space need only be allocated for the actual number of elements in a particular array. If the matrix is large, there is the option to store the elements in an external disk file. In addition, in order to work with large numbers of matrices, arrays of handles are used, with each handle pointing to either a single subscripted array or an external disk file. Space for a matrix array, therefore, is allocated only as the need arises.

The value of a matrix parser can be seen in many scientific and engineering areas, especially where linear algebra is involved - solving systems of linear equations, obtaining the solution to a least squares problem, making vector computations, and a myriad of other applications. Linear algebra is an exciting field in its own right and working with matrices can be fun.

The Pascal program listed in the appendix was written using THINK Pascal 2.0. It is a stand-alone interactive program allowing a user to enter and evaluate matrix expressions at the keyboard. The listing can be obtained by writing the author or Xplain Corporation (the publishers of MacTutor).

The typical matrix operations of transpose, multiplication, addition, inverse, etc., are discussed, as well as array operations and using functions for elements of a matrix. Examples are given and it will be shown how easily the solution to a least squares problem can be obtained.

Enjoy! Read on!

Matrix

Visually we represent a matrix as a rectangular array of reals [arranged in rows and columns] which possesses certain properties. An m x n matrix A [m rows and n columns] will look like the following

where the ijth element (ith row, jth column) is (aij).

The transpose of a matrix is a matrix with the rows and columns interchanged. The (i,j)th element of the transpose of matrix A is the (j,i)th element of A. We can write the transpose of A [ or A'] as

A vector is a matrix having one column (or row). A column vector with m rows can be written

And, V' will be a row vector with m columns. In the following, when we talk about a vector we will consider it to be a column vector.

We have become so accustomed to this two-dimensional representation that we tend to impose the structure in our code by double dimensioning matrix arrays. The computer, however, works with and stores the elements of a matrix sequentially. In the following section, we will show how valuable RAM can be saved by using single subscripted arrays.

Representation of a matrix

We represent a matrix with a handle as this allows us to dynamically allocate and store the elements in a relocatable block of memory on the heap (and not the stack) using NewHandle. The stack is a good place for allocating a few objects at a time, using them, and then removing them. The heap, however, is the best place for a data structure such as a matrix array. Since we don’t know in advance how many matrices will be required in our computations, it seems logical to use the NewHandle function to allocate space for our arrays, creating and disposing of relocatable blocks of memory as needed.

Let us compare the difference between allocating space for a singly subscripted array and allocating space for a doubly subscripted array using the NewHandle function.

In allocating space for a double subscripted matrix array, the allocation must be done in blocks equal to the number of columns. For example, let us define the following variable types,

matrixdouble = array[1..maxrows, 1..maxcols] of extended;
ptrmatrixdouble = ^matrixdouble;
hdlmatrixdouble = ^ptrmatrixdouble;

where maxrows and maxcols are global constants. Then, if ‘matrix’ is an (m x n) matrix variable of type ‘hdlmatrixdouble’, space [in number of bytes] is dynamically allocated by the following

blocksize:= m * maxcols * 10;
matrix:=hdlmatdouble(NewHandle(blocksize));

where blocksize is the Pascal type, Size [longint], and is equal to the number of required bytes, and m and n are long integers. The elements, matrix^^[i,j], are extended real numbers [10 bytes for an extended real data type].

Now, let us allocate space for a singly dimensioned matrix variable, and define the following variable types

matrixsingle = array[1..maxelements] of extended;
ptrmatrixsingle = ^ matrixsingle;
hdlmatrixsingle = ^ptrmatrixsingle;

where ‘maxelements’ is a global constant. Then, if ‘matrix’ is an (m x n) matrix variable of type ‘hdlmatrixsingle’, space is allocated by the following

blocksize:= m * n * 10;
matrix:=hdlmatrixsingle(NewHandle(blocksize));

where m*n <= maxelements.

It can be seen that no matter how large the dimension of the argument in the type declaration (maxelements), space need only be allocated for the actual number of elements in a particular array. (Actually, we allocate an additional 20 bytes, since the number of rows and columns are embedded in the first two elements of the array). The percent savings using single instead of double subscripted arrays is 100(1- (n/maxcols)). For maxcols large and n small, the savings can be considerable. If n = maxcols, there is no savings.

We note that in declaring a variable, we pre-allocate on the stack that amount of storage that that type of variable demands. If it is an array of extended reals, for example, we have to allocate 10 times the number of elements in the array (10 bytes for an extended real). If our variable, however, is a handle to an array of extended reals, we only allocate 4 bytes of storage on the stack (4 bytes for a handle). It is the NewHandle function which then allows us to allocate just the right amount of space needed for a particular array. We must, however, make the transformation from double subscripting to single subscripting.

Thus, we are making provision to work with large matrices, but at the same time are not being penalized for working with smaller ones. Also, we are not held to a prescribed number of rows and columns. The limiting factor is the total number of elements in an array. Using NewHandle and long integers to index an array, the global constant ‘maxelements’ can be made quite large. [In Globals, this is set equal to 2,000,000].

Allocating space for a matrix using the NewHandle function also helps us when we wish to work with a large number of matrices. In this case we define an array of handles each of which points to either a singly subscripted array of reals or to an external disk file of reals. We then define a handle to this array of handles. (This may be confusing, but there’s a method to our madness!). For N handles in an array, only 4*N bytes of storage need be allocated on the heap [a handle is worth 4 bytes]. If, instead of a handle, we had used an array of extended reals as an element of the array, the memory requirements might easily exceed that needed to run the program. Thus, we can store our matrices in an indexed array, allocating space for each new one as the need arises. And, we don’t pay a penalty in the process.

Summing up, we conserve on RAM by: (1) representing a matrix with a handle which points to a single [as opposed to double] subscripted array of reals dynamically allocated on the heap using the NewHandle function; (2) storing the matrix [if large] in an external disk file; (3) using an array of handles [to single dimensioned arrays of reals or to external disk files of reals] in order to work with [large numbers] of matrices.

Housekeeping

Since the elements of a matrix are stored either in RAM or in an external file, there are some housekeeping chores associated with each. The following may appear a bit cumbersome, but it does the job and keeps track of things pretty well.

Associated with each matrix are a name and a pointer [plus other ancillary information]. The variable, strvar, is a handle to an array of handles, the ith one, strvar^^[i], pointing to the string [20] variable, strvar^^[i]^^, the name of the ith matrix. The long integer, i, is a pointer to either the ith handle to an external disk file of extended reals, matfile^^[i], or to the ith handle to a singly subscripted array of extended reals, storematrix^^[i], depending upon the value of the boolean variable, matrixstoredinfile^^[i].

If the number of elements in, strvar^^[i]^^, is greater than or equal to, ‘bignumber’, (input long integer with a default of 100), matrixstoredinfile^^[i] is set to true, and the elements are stored in, matfile^^[i]^^. If the number is less than ‘bignumber’, matrixstoredinfile^^[i] is set to false, and the elements are stored in, storematrix^^[i]^^[j], a relocatable block of memory on the heap [j runs from 1 to the number of elements plus two, with the number of rows and number of columns being stored in the first two elements of a matrix array].

If the matrix elements are stored in an external disk file, there are two additional booleans to contend with- mfilenew^^[i], and mfileopen^^[i]. Mfilenew^^[i] is true if NewHandle has been called to allocate space for matfile^^[i], and, mfileopen^^[i] is true if matfile^^[i]^^ is open.

If the matrix elements are stored in RAM [on the heap], there is only one boolean to contend with - matrixnew^^[i]. Matrixnew^^[i] is true if NewHandle has been called to allocate space for, storematrix^^[i].

It will be seen later [in the discussion of EvaluateNodes] that node matrices are stored on the heap, as well as matrices [amat, bmat, and cmat] which are used in binary and function calculations. Once we are through with these matrices, their handles are disposed of, and RAM space is restored. It is only when we wish to save a matrix that we have to be concerned with “where it goes” - in RAM or out of RAM.

And that’s all there is to the housekeeping!

Matrix multiplication and code.

Since a number of published algorithms for matrix operations use doubly subscripted variables for matrix arrays, we will show how the Pascal code for matrix multiplication can be modified to use singly subscripted arrays, and will give the linear transformation which takes us from doubly to singly dimensioned arrays.

Matrix multiplication (the symbol * is used to distinguish this type of multiplication from array multiplication #) for two matrices A and B is defined only if the number of columns of A equals the number of rows of B. The (i,j)th element of the product A*B is given by the sum of the products of the individual elements of the ith row of A with the elements of the jth column of B [In vector language, the (i,j)th element in the product is the dot product of the ith row of A with jth column of B]. Matrix multiplication is not commutative. If A and B are both square with the same number of rows and columns, A*B <> B*A.

The following illustrates the difference in code between handling matrix multiplication with singly dimensioned arrays and matrix multiplication with doubly dimensioned arrays.

Using doubly subscripted arrays, the [usual] Pascal code for multiplying A and B to obtain C, where A is an (m x l) and B an (l x n) matrix, is given by

{1}

for i:=1 to m do
 for j:=1 to n do
 begin
 sum:=0;
 for k:=1 to l do
 sum:=sum + a^^[i,k]*b^^[k,j];
 c^^[i,j]:=sum;
 end;

The above loop can be modified using singly dimensioned arrays.

{2}

for i:=1 to m do 
 for j:=1 to n do 
 begin  
 sum:=0;
 for k:=1 to l do
 sum:=sum + a^^[(i-1)*n + k]*b^^[(k-1)*n + j];
 c^^[(i-1)*n + j]:=sum;
 end;

It can be seen that the linear transformation to get the double subscripted (i,j)th element as a single subscripted variable, say, k, is given by, k = (i-1)*n + j, where n is the number of columns in the matrix.

As an example of matrix multiplication, we have the following,

In general, matrix multiplication is not commutative.

Array Operations

Array operations are element-by-element arithmetic operations performed on matrices having the same numbers of rows and columns. The operations which can be performed on these matrices are: add (+), subract (-), multiply (#), divide-by (/), divide-into (\), raise to a power (^). As a multiplication example we have the following.

Structure of the parser algorithm

The present algorithm identifies the basic elements in a matrix expression, transforming them into an ordered set of tokens from which a node table is constructed. The node table is then used to calculate the result of the expression - a real matrix. Besides including matrices as tokens, this algorithm differs from that described in [1] in that the values for operand variables are pointers (long integers) to matrices, instead of real numbers. However, since a real number is defined as a matrix with one row and one column, all real number calculations can be performed with the present code.

In the main driver, ParserDriver, the user inputs a str255 variable, line, containing either a command or a matrix expression. The commands are: ‘quit’ (quit the program), ‘changebig’ (change ‘bignumber’), ‘creatematrix’ (create either a random matrix or input one from the keyboard), ‘readmatrix’ (display a stored matrix on the screen), ‘dec’ (change the number of decimal places), ‘cls’ (clear the screen), ‘clm’ (clear memory of all variable names and values), ‘delete’ (delete a variable), ‘listv’ (list the variable names in storage and their values - real numbers or pointers to matrices). If there are no commands, eval is called, and a real number, a pointer to a matrix, or an error message is returned in ‘result’, a str255 variable.

Eval calls the following procedures: LexicalAnalysis, Parser, SetValues, and EvaluateNodes.

LexicalAnalysis inputs, line, transforming the matrix expression into an ordered set of tokens, their types, and precedence values. The tokens, sy^^[i]^^, [symbols, words, numbers, etc.] are string[20] variables, and their types, tokentype^^[i]^^, are ‘binary’, ‘unary’, ‘constant’, ‘variable’, and ‘function’. The precedence values [long integers] associated with each token, pr^^[i]^^, [ i = 1 to ntot], impose a certain hierarchy or ordering of functions and operations, determining which are performed first, etc., e.g., multiplication before addition or subtraction.

Eval next checks the ‘variable’ tokens against a stored list of matrix names. If there is no match for a particular ‘variable’ token, an error message is displayed on the screen. If all the ‘variable’ names in the expression, sy^^[i]^^, match stored names, strvar^^[j]^^, for some i and j, tokentype^^[i]^^, is changed from ‘variable’ to ‘matrix’, and the revised set of tokens, types, and precedence values, are then input into Parser.

With the input set of ordered tokens, types, and precedence values, Parser constructs a node table - an indexed array of node records (junction points or nodes) in the evaluation of the expression. Each indexed entry, a node description, contains sufficient information for a ‘node’ matrix to be calculated and stored. (Node matrices are temporarily stored in RAM until the final expression has been evaluated. The handles to these matrix arrays are then disposed of, freeing up RAM.). The ith indexed record contains the following six fields:

nodetable^^[i]^^.optype,  type of operation [function, binary, etc.]
nodetable^^[i]^^.roptype, token type for right operand               
                    nodetable^^[i]^^.loptype,                  token 
type for left operand                                                
                    nodetable^^[i]^^.op.index,                 operator 
token                                 nodetable^^[i]^^.rop.index,
 right operand token                                 nodetable^^[i]^^.lop.index,
 left operand token
                               

Eval next inputs the node table into SetValues where pointers to indexed arrays of matrices [long integers] are substituted for operand matrix names. Constants appearing in operand fields are later embedded in matrices in EvaluateNodes.

The final call in Eval is to the procedure, EvaluateNodes. If there are no errors in the construction of the matrix expression, t^^[numnodes], is set equal to either a real number or a pointer to a matrix, and is returned to Eval, where it is embedded in ‘result’.

EvaluateNodes

This procedure illustrates how the parsing algorithm uses constants and pointers [embedded in operand fields of the node table] to calculate matrices. The pointer to a matrix is a long integer and is equal to (or “points to”) the index of a particular matrix stored in an array of matrices, in RAM or in an external disk file.

EvaluateNodes runs through the indexed records of the node table, from i = 1, numnodes, calculating at the ith node a matrix, nodematrix^^[i], which is stored temporarily on the heap, and given the pointer, t^^[i] = i. If the input expression (line) contains an equals sign (assignment statement), the final node matrix is stored in RAM or an external disk file, depending upon the size of the matrix (the number of elements is compared with the global variable, ‘bignumber’).

There are four procedures which are part of EvaluateNodes and which are called from within the main procedure: OpenMatrixFile, GetMatrix, GetNodeMatrix, and GetConstantMatrix (the procedure, GetMatrix, calls OpenMatrixFile). In each procedure, the NewHandle function is called to allocate space for, dummymatrix, a handle to a singly subsripted array of extended reals.

OpenMatrixFile opens an external disk file of extended reals, matfile^^[matpointer]^^, given the pointer, matpointer. As mentioned above, mfileopen^^[matpointer], keeps track of whether the file is open or closed. If closed, the procedure opens it, and sets mfileopen^^[matpointer] to true. The file is reset, and the values of the matrix are read into, dummymatrix. The file is then closed, and mfileopen^^[matfile] set to false.

Depending upon the value of, matrixstoredinfile^^[matpointer], GetMatrix reads the contents of, matfile^^[matpointer]^^ or storematrix^^[matpointer]^^, into dummymatrix.

GetNodeMatrix stores the elements of, nodematrix^^[matpointer], a previously calculated node matrix, into dummymatrix.

GetConstantMatrix stores the constant, real number, into dummymatrix, a matrix with one row and one column.

Let’s look at the structure of EvaluateNodes.

After initializing ‘error’ to the null string, the matrix handles, amat, bmat, and cmat, are set to nil, and their corresponding boolean variables, anew, bnew, and cnew, are set to false. These boolean variables are true if the NewHandle function has been called to allocate space for the corresponding matrices. The matrices, amat, bmat, and cmat, are used later in binary and function calculations.

As we step through the ith record in the node table, we first allocate space for, t^^[i], the pointer to, nodematrix^^[i]. We read ‘rop.index’ (using ‘readstring’) to obtain, b2, and round this to, mn, as we will need a long integer if we have to get a stored matrix with GetMatrix or a node matrix with GetNodeMatrix. [The variable, t^^[i], is an extended real, since, if the result of our calculations is real, we set t^^[i] equal to it].

We check the ‘roptype’ field next. If this is ‘node’, we obtain the pointer to the mnth matrix, t^^[i]^^, setting, b2 = t^^[mn], mn = round(b2), and call GetNodeMatrix. This puts the contents of, nodematrix^^[mn], into ‘dummymatrix’. If ‘roptype’ is ‘constant’, we call GetConstantMatrix with b2, and embed it in ‘dummymatrix’. If ‘roptype’ is ‘matrix’, we call GetMatrix which puts the contents of the mnth stored matrix into ‘dummymatrix’.

As we need a right operand matrix to operate with, we allocate space for bmat, set bnew to true, and put the contents of dummymatrix into bmat. We set m2 equal to the number of rows, and n2 equal to the number of columns of bmat.

We next check for an equals sign in the ‘op.index’ field (an assignment). If there is one, the contents of the matrix pointed to in the ‘rop.index’ field (which are now in dummymatrix), are read into, nodematrix^^[i], and the index, i, is upped one.

The variable, matrixoper, is set equal to ‘op.index’, and the ‘optype’ field is checked for a ‘unary’ or ‘function’ value.

If the ‘optype’ field is ‘unary’ or ‘function’, there are a number of cases to examine. If, matrixoper, is ‘minus’, dummymatrix is set equal to the negative of, dummymatrix. If a ‘quote’, the transpose of dummymatrix, is calculated and returned as, dummymatrix. If ‘matrixoper’ is ‘inv’, the pseudo inverse of dummymatrix is obtained, only if the number of rows is greater than or equal to the number of columns. If ‘matrixoper’ is none of these, dummymatrix is input into ‘matrixfunctions’, and becomes on output, a matrix with each element operated on by ‘matrixoper’. This latter matrix then becomes the node matrix for the ith node, and i is upped one.

If there is neither an equality, nor a ‘function’, nor a ‘unary’ operation, we read the left operand field, lop.index. Similar to the procedure discussed above for the, rop.index field, the extended real, b1, is read from, lop.index, then rounded to, lm, and we obtain a left operand matrix, amat.

The left operand matrix, amat, operates on the right operand matrix, bmat, and we obtain on output, cmat, which then becomes our ith node matrix. The blocksize for, cmat, depends upon matrixoper. If the latter is an asterisk, we set blocksize equal to the number of elements in the matrix resulting from the matrix multiplication of, amat, times, bmat. If not an asterisk, the blocksize is set equal to the product of the maximum number of rows (in, amat, or, bmat), and the maximum number of columns (in, amat, or bmat). The procedure, matrixoperations, is called, the resulting matrix returned in, cmat, and the contents of, cmat, are read into, nodematrix^^[i].

At the end of our ith loop, we dispose of the handles, amat, bmat, and cmat, if they are not nil, and reset anew, bnew, and cnew, to false, accordingly.

After the last node calculation, matrices and pointers for all previous nodes are disposed of and we get back some heap space.

If there is an assignment statment in, line, (save[2] is an equals token), the calculated matrix, nodematrix^^[numnodes], depending upon its size, is stored in RAM, or in an external disk file. If, nodematrix^^[numnodes], has one row and one column, t^^[numnodes] is set equal to the value of the real number in the third element of the array, and is printed to the screen (in ParserDriver).

Finally, space taken up on the heap by, dummymatrix, and, nodematrix^^[numnodes], is restored with DisposHandle, and these handles are set to nil. All variables are “cleaned up” (duplicate names eliminated) by the call to, cleanupvariables.

Examples

We will now give some examples to illustrate how the program may be used interactively to calculate matrices. It should be noted that although matrix examples are given, the program can be used to make calculations with real numbers. The commands typed in by the user as well as the responses from the program will be displayed. It is rather easy to follow.

In the first example we will create a random (4 x 3) matrix, ‘a’, where the elements of the matrix are random integers between -9 and 9 inclusive. We print it to the screen by typing, a, then set, b, equal to the “transpose of a,” i.e., b = a'. The newly calculated matrix, b, is automatically displayed on the screen.

We next multiply the “transpose of a”, times, a, and list our variables.

In the following, we type in the result of the last calculation, ans, which will be equal to the matrix, c, just calculated.

Next, we take the natural logarithm of the absolute value of matrix, c, which is displayed immediately.

By raising, ans, to the base of the natural logarithms, we get back the absolute value of, c.

In the next example we will recover the solution to a linear system of equations, ax = b. We first create a (3 x 1) vector, x, [where x' = (1 2 3)], multiply the matrix, a, times x to obtain the righthand side, b. We will then ‘recover’ the solution, lsqsol, using the pseudo inverse, inv(a), and calculate the difference vector, (lsqsol - x), which will contains zeros.

From the above it can be seen that the solution has been recovered.

Our last example shows how we can use functions as elements within a matrix.

Finally, we quit our program, saving variables and files.

Conclusion

This article has described and illustrated an efficient parser for evaluating matrix expressions. One of the main features has been the development of a method for working with large matrix arrays. Using long integers to index arrays and handles to point to single subscripted (as opposed to double) arrays, we are able to dynamically allocate large relocatable blocks of memory on the heap with the NewHandle function. If the number of matrix elements is large, there is the option to store the elements in an external disk file. Using arrays of matrix handles to point to either single subscripted arrays or to external disk files, we are able to work with large numbers of matrices.

Since a real number is treated as a matrix with one row and one column, calculations with real numbers can be made using the parsing algorithm.

The Pascal code can be incorporated within a larger program to make it more flexible. The key function is eval which returns the results of a matrix or real calculation, given an input matrix expression (str255 Pascal variable type).

The code listing for the stand-alone program described in this article can be obtained by writing Xplain or the author. A more extensive program which incorporates building one’s own text files (using a simple programming language), matrix partitioning, eigenvalue/eigenvector computations, singular value decomposition, and many other linear algebra routines, is obtainable from the author. The program which is well-documented is available for $75. (shipping $3). The complete code listing for the matrix parser is $155 (shipping $3). These can be obtained from: Greer Software Products, Box 268, Annandale, Virginia 22003, (703) 978-3327.

References

[1] “Pascal Procedures, A Practical Parser”, MacTutor, May 1991, by Bill Murray.

Acknowledgement

The author wishes to thank Richard F. Thompson, Engineering Systems, Inc., Vienna, Virginia, for many helpful suggestions and discussions, and, especially, for much of the motivation behind this article.