go to Trig home page   Guide to GAUSS Programming - a basic introduction


 
Matrix algebra and manipulation
1 Matrix algebra

Algebra involving matrices translates almost directly from the page into GAUSS. At bottom, most mathematical statements can be directly transcribed, with some small changes.

1.1 The basic operators

GAUSS has eight mathematical operators and six relational ones. The mathematical ones are

+
Addition
-
Subtraction
*
Multiplication
/
Division
'
Transposition
%
Modulo division
!
Factorial
^
Exponentiation

and the six relational operators are:

==
EQ
equals
/=
NE
does not equal
>
GT
greater than
<
LT
less than
>=
GE
greater than/equals
<=
LE
less than/equals

Either the symbols or the two-letter acronyms may be used.

Warning

Note the double-equals sign for equivalence. This must not be confused with the single-equals sign implying assignment. The two return very different results:
mat = 5; mat is assigned the value 5; the "result" of this operation is 5
mat ==5; mat is compared to the value 5; the "result" of this operation is "true" if mat is equal to 5, "false" otherwise

With respect to logical results, GAUSS standard procedures use the convention

"false" = 0    "true" /= 0

and there are four logical operators for these which all return true or false

NOT var1

true if var1 false, and vice-versa
var1 AND var2

true if var1 true and var2 true, else false
var1 OR var2

true if var1 true or var2 true, else false
var1 XOR var2

true if var1 true or var2 true but not both, else false
var1 EQV var2

true if var1 is equivalent to var2 i.e. both true or both false


Warning

The GAUSS manuals state that procedures set variables to 1 to signify true and 0 to false, but this is not strictly necessary - nor is it adhered to, despite several functions depending upon it. Do not rely on true==1 (eg if x==1 then...). Instead, use true/=0 (eg, if x /= 0 then...) Better still, do not rely on a particular mathematical value for true or false.

GAUSS is a "strict" language: if a logical expression has several elements, all the elements of the expression will be checked even if the program has enough information to return true or false. Thus using these logical statements may be less efficient then, for example, using nested IF statements. This is also different from the way some other programs operate.

Operators work in the usual way. Thus these operations on matrices a to e are, subject to conformability requirements, all valid operations:

a = b+c-d;
a = b'*c';
a = (b+c)*(d-e);
a = ((b+c)*(d+e))/((b-c)*(d-e));
a = (b*c)';

Notice from this that matrix algebra translates almost directly into GAUSS commands. This is one of GAUSS's strong points. GAUSS will check the conformability of the above operations and reject those it finds impossible to carry out; however, see section 1.2 below.

The order of operation is complex; see the section on operators in the manual for details. But essentially the order is left to right with the following rough precedence:

brackets
transposition
factorial
exponentiation
negation
multiplication and division
addition and subtraction
dot relational operators
dot logical operators
relational operators
logical operators
row and column indices

See the next section for an explanation of dot operators.

The division operator can be used like any other. When one or other variable is a scalar, then the division operation will be carried on an element-by-element basis (see below). However, when the variables are both matrices then GAUSS will compute a generalised inverse; that is, a = b/c is deemed to be the solution to ca = b which leads to the equations

a=b/c   =>a=c-1b (c square)   ora=(c'c)-1c'b (c non-square)

Therefore, if two matrices are divided, then it may be preferable to do the inverse explicitly rather than leave the calculation to GAUSS. Division is a common source of unnoticed errors, because GAUSS will try as hard as possible to find an appropriate inverse.

There are two concatenation operators:

~    horizontal concatenation
|    vertical concatenation

These add one matrix to the right or bottom of another. Obviously, the relevant rows and columns must match. Consider the following operations on two matrices, a and b, with ra and rb rows and ca and cb columns, and the result placed in the matrix c:

dimensions of a dimensions of b operation dimensions of c condition
ra x ca rb x cb c = a ~ b ra x (ca + cb) ra = rb
ra x ca rb x cb c = a | b (ra + rb) x ca ca = cb

Parts of matrices may be used, and results may be assigned to matrices or to parts:

a = b*c;
a = b[r1:r2,c1]*c[r3, c2:c3];
a[r1, c1:c2] = b[r1,.]*c;

subject to, in the last case, the recipient area being of the correct size.

These operations are available on all variables, but obviously "a=b*c" is nonsensical when b and c are strings or character matrices. However, the relational operators may be used; and there is one useful numerical operator - addition:

a = b $+ c;

This appends c to b. Note that the operator needs the string signifier "$" to inform GAUSS to do a string concatenation rather than a numerical addition. If you omit the $ GAUSS will carry out a normal addition.

For example,
b = "hello";
c = "mum";
a = b $+ " " $+ c;
PRINT $a;

will lead to "hello mum" being printed.

With character matrices, the rules for the conformability of matrices and the application of the operator are the same as for mathematical operators (see the next section). Note that, in contrast to the matrix concatenation operators, the overall matrix remains the same size (strings grow) but each of the elements in the matrix will be changed. Thus if a is an r by c matrix of file names,

a = a $+ ".RES";

will add the extension ".RES" to all the names in the matrix (subject to the eight-character limit) but a will still be an r by c matrix. If any of the cells then have more than eight characters, the extra ones are cut off.

String concatenation applied to strings and string arrays will cause these to grow.

Strings and character matrices may be compared using the relational operators. The string signifier $ is not always necessary, but it makes the program more readable and may avoid unexpected results.

In the eight bytes of data used for each matrix cell, characters and numbers are stored in different ways. GAUSS uses the $ symbol to signify the byte order, but otherwise makes no distinction between characters and numbers. So if you mix data types, omit a $ sign or put one in where it shouldn't be, GAUSS will not complain but the result will be gibberish.

1.2 Conformability and the "dot" operators

GAUSS generally operates in an expected way. If a scalar operand is applied to a matrix, then the operation will be applied to every element of the matrix. If two matrices are involved, the usual conformability rules apply:

Operation Dimensions of b Dimensions of c Dimensions of a
a = b * c; scalar 4x2 4x2
a = b * c; 3x2 4x2 illegal
a = b * c'; 3x2 4x2 3x4
a = b + c; scalar 4x2 4x2
a = b - c; 3x2 4x2 illegal
a = b - c; 3x2 3x2 3x2

and so on. However, GAUSS allows most of the mathematical and logical operators to be prefixed by a dot:

a = b.>c; a = (b+c).*d'; a = b.==c;

This tells the machine that operations are to be carried out on an "element by element" basis (or ExE, as the oracular manual so succinctly puts it). This means that the operands are essentially broken down into the smallest conformable elements and then the scalar operators are applied. How this works in practice depends on the matrices.

To give an example, suppose that mat1 is a 5x4 matrix. Then the following results occur for multiplication:

Operation mat2 r x c Result
mat1 * mat2 scalar 5x4; mat2 times each element of mat1
mat1 .* mat2 5x4 5x4; mat1[i,j] * mat2[i,j] for all i, j (Hadamard product)
mat1 .* mat2 5x1 5x4; the ith element in mat2 is multiplied by each element in the ith row of mat1
mat1 .* mat2 1x4 5x4; the jth element in mat2 is multiplied by each element in the jth column of mat1
mat1 .* mat2 anything else illegal

Similarly for the other numerical operators:

Operation mat2 r x c Result
mat1 ./ mat2 5x4 5x4; mat1[i,j] / mat2[i,j] for all i, j
mat1 .% mat2 1x4 5x4; modulus mat1[i,j] / mat2[j] for all i, j
mat1 .*. mat2 5x4 25x16; mat1[i, j] * mat2 for all i,j (Kronecker product)

Warning

The dot operators do not work consistently across all operands. In particular, for addition and subtraction no dot is needed.

1.3 Relational operators and dot operators

For the relational operators, the results are slightly different. These operators return a scalar 0 or 1 in normal circumstances; for example, compare two conformable matrices:

mat1 /= mat2    mat1 GT mat2

The first returns "true" if every element of mat1 is not equal to every corresponding element of mat2; the second returns "true" if every element of mat1 is greater than every corresponding element of mat2. If either variable is a scalar than the result will reflect whether every element of the matrix variable is not equal to, or greater than, the scalar. These are all scalar results.

Prefixing the operator by a dot means that the element-by-element result is returned. If mat1 and mat2 are both r by c matrices, then the results of

mat1 ./= mat2    mat1 .GT mat2

will be a r by c matrix reflecting the element-by-element result of the comparison: each cell in the result will be set to "true" or "false". If either variable is a scalar than the result will still be a r by c matrix, except that each cell will reflect whether the corresponding element of the matrix variable is not equal to, or greater than, the scalar.

1.4 Fuzzy operators

In complex calculations, there will always be some element of rounding. This can lead to erroneous results from the relational operators. To avoid this, fuzzy operators are available. These are procedures which carry out comparisons within tolerance limits, rather than the exact results used by the non-fuzzy operators. The commands are

FEQFNEFGTFLTFGEFLE

with corresponding dot operators

DOTFEQDOTFNEDOTFGTDOTFLTDOTFGEDOTFGE

and are used, for example FEQ, by

result = FEQ (mat1, mat2);

This will compare mat1 and mat2 to see whether they are equal within the tolerance limit, returning "true" or "false". Apart from this, the fuzzy operators (and their dot equivalents) operate as the exact relational operators.

The tolerance limit is held in a variable called _fcmptol which can be changed at any time. The default tolerance limit is 1.0x10-15. To change the limit simply involves giving this variable a new value:

_fcmptol = newValue;


2 Set operations

Column vectors can be treated like sets for some purposes. GAUSS provides three standard procedures for set operation:

unVec = UNION (vec1, vec2, flag);
intVec = INTRSECT (vec1, vec2, flag);
difVec = SETDIF (vec1, vec2, flag);

where unVec, intVec, and difVec are the results of union, intersection, and difference operations on the two column vectors vec1 and vec2. The scalar flag is used to indicate whether the data is character or numeric: 1 for numeric data, 0 for character. The difference operator returns the elements of vec1 not in vec2, but not the elements of vec2 not in vec1.

These commands will only work on column vectors (and obviously scalars). The two vectors can be of different sizes. A related command to the set operators is

unVec = UNIQUE (vec, flag);

which returns the column vector vec with all its duplicate elements removed and the remaining elements sorted into ascending order.


3 Special matrix operations

GAUSS provides methods to create and manipulate a number of useful matrix forms. The commonest are covered in this section. A fuller description is to be found in the GAUSS Command Reference.

3.1 Some useful matrix types

Firstly, three useful matrix creating operations:

identMat = EYE (iSize);
onesMat = ONES (onesRows, onesCols);
zerosMat = ZEROS (zeroRows, zeroCols);

These create, respectively: an identity matrix of size iSize; a matrix of ones of size onesRows by onesCols; and a matrix of zeroes of size zeroRows by zeroCols. Note the US spelling.

3.2 Special operations

A number of common mathematical operations have been coded in GAUSS. These are simple to use to use and more efficient then building them up from scratch. They are

invMat = INV (mat);
invPDMat = INVPD (mat);
momMat = MOMENT (mat, missFlag);
determ = DET (mat);
determ = DETL;
matRank = RANK (mat);

The first two of these invert matrices. The matrices must be square and non-singular. INVPD and INV are almost identical except that the input matrix for INVPD must be symmetric and positive definite, such as a moment matrix. INV will work on any square invertible matrix; however, if the matrix is symmetric, then INVPD will work almost twice as fast because it uses the symmetry to avoid calculation. Of course, if a non-symmetric matrix is given to INVPD, then it will produce the wrong result because it will not check for symmetry.

GAUSS determines whether a matrix is non-singular or not using another tolerance variable. However, even if it decides that a matrix is invertible, the INV procedure may fail due to near-singularity. This is most likely to be a problem on large matrices with a high degree of multicollinearity. The GAUSS manual suggests a simple way to test for singularity to machine precision, although I have found it necessary to augment their solution with fuzzy comparisons to ensure a workable result (for an example, see the file SingColl.GL on the code page).

The MOMENT function calculates the cross-product matrix from mat; that is, mat'*mat. For anything other than small matrices, MOMENT(x, flag) is much quicker than using x'x explicitly as GAUSS uses the symmetry of the result to avoid unnecessary operations. The missFlag instructs GAUSS what to do about missing values (see below) - whether to ignore them (missFlag=0) or excise them (missFlag=1 or 2).

DET and DETL compute the determinants of matrices. DET will return the determinant of mat. DETL, however, uses the last determinant created by one of the standard functions; for example, INV, DET itself, decomposition functions all create determinants along the way. DETL simply reads this value. Thus DETL can avoid repeating calculations. The obvious drawback is that it is easy to lose track of the last matrix passed to the decomposition routines, and so determinants should be read as soon as possible after the relevant decomposition function has been called. See the Command Reference for details of which procedures create the DETL variable.

RANK calculates the rank of mat.

3.3 Manipulating matrices

There are a number of functions which perform useful little operations on matrices. Commonly-used ones are:

vec = DIAG (mat);
mat = DIAGRV (vec);
newMat = DELIF (oldMat, flagVec);
newMat = SELIF (oldMat, flagVec);
newMat = RESHAPE (oldMat, newRows, newCols);
nRows = ROWS (mat);
nCols = COLS (mat);
maxVec = MAXC (mat);
minVec = MINC (mat);
sumVec = SUMC (mat);

DIAG and DIAGRV abstract and insert, respectively, a column vector from or into the diagonal of a matrix.

DELIF and SELIF allow certain rows and columns to be deleted from the matrix oldMat. The column vector flagVec has the same number of rows as oldMat and contains a series of ones and zeros. DELIF will delete all the rows from the matrix for which there is a corresponding one in flagVec, while SELIF will select all those rows and throw away the rest. Therefore DELIF and SELIF will, between themselves, cover the whole matrix.

DELIF and SELIF must have only ones and zeros in flagVec for the function to work properly. This is something to consider as the vector flagVec is often created as a result of some logical operation. For example, to delete all the rows from matrix mat1 whose first two columns are negative would involve

flags = (mat1[1,.] .< 0) .AND (mat1[2,.] .< 0);
mat2 = DELIF (mat1, flags);

This particular example should work on most systems, as the logical operator AND only returns 1 or 0. But because true is really non-zero (not 1) some operations could lead to unexpected results.

DELIF and SELIF also use a lot of memory to run. A program calling these procedures often would be improved by rewriting them (versions can be downloaded from the Web; see the appendix).

ROWS and COLS return the number of rows and columns in the matrix of interest.

MAXC, MINC, and SUMC produce information on the columns in a matrix. MAXC creates a vector with the number of elements equal to the number of columns in the matrix. The elements in the vector are the maximum numbers in the corresponding columns of the matrix. MINC does the same for minimum values, while SUMC sums all the elements in the column. However, note that all these functions return column vectors. So, to concatenate onto the bottom of a matrix the sum of elements in each column would require an additional transposition:

sums = SUMC(mat1);
mat1 = mat1 | sums';

On the other hand, because these functions work on columns, then calling the functions again on the column vectors produced by the first call allows for matrix-wide numbers to be calculated:

maxMat=MAXC(MAXC(mat1));
minMat=MINC(MINC(mat1));
sumMat=SUMC(SUMC(mat1));

will return the largest value in mat1, the smallest value, and the total sum of the elements.

4 Missing values

GAUSS has a number of "non-numbers" which can be used to signify missing values, faulty operations, maths overflow, and so on. These NANs (in GAUSS's terms) are not values or numbers in the usual sense; although all the usual operations could be carried out with them, the results make no sense. These are just identifiers which GAUSS recognises and acts upon.

Generally GAUSS will not accept these values in numerical calculations, and will stop the program. However, the string operators can be used on these values to test for equalities. To see if the variable var is one of these odd values or not, the code

var $== TestValue    orvar $/= TestValue

would work. The other relational operators would work as well, but the result is meaningless. The TestValues are scattered around the GAUSS manual in excitingly unpredictable places.

With empirical datasets, the largest problem is likely to be with missing values. These missing values will invalidate any calculation involving them. If one number in a sequence is a missing value, then the sum of the whole sequence will be a missing value; similarly for the other operators. Thus checking for missing values is an important part of most programs.

Missing values can have their uses. They can indicate that a program must stop rather than go any further; they can also be used as flags to identify cells. To this end we have three functions

newMat = MISS (oldMat, badValue);
newMat = MISSRV (oldMat, newValue);
newMat = MISSEX (oldMat, mask);

The first of these converts all the cells in oldMat with badValue into the missing value code. MISSRV does the opposite, replacing missing values in oldMat with newValue. The second can be used to remove missing values from a matrix; however, in conjunction with the first, it can be used to convert one value into another. For example, to convert all the ones in mat1 into twos could be done by:

tempMat = MISS (mat1, 1);
mat1 = MISSRV (tempMat, 2);

This of course assumes that mat1 had no prior missing values to be erroneously convered into twos. MISSEX is similar to MISS, except that instead of checking to see which elements of the matrix mat1 match badValue, GAUSS takes instructions from mask, a matrix of ones and zeros of the same size as mat1. Any ones in mask will lead to the corresponding values in mat1 being changed into missing values. MISS and MISSEX are thus very similar in that

MISS (mat1, 2); is virtually equivalent to MISSEX (mat1, mat1.==2);

To test for missing values, use

missing = ISMISS (mat);
missing = SCALMISS (mat);

The first of these tests to see whether mat contains any missing values, returning one if it finds any and zero otherwise; the second returns one only if mat is a scalar and a missing value.

4.1 Non-fatal use of missing values - DOS versions of GAUSS

This section relates to DOS versions of GAUSS. Unix and NT-based Windows software isolate system exceptions, and so GAUSS no longer stops on maths processor overflows or underflows. Thus in newer versions of GAUSS DISABLE (see below) is effectively always on. You can access the system interrupts if you desperately want to but there is little need. ENABLE, DISABLE, NDPCNTRL and other system settings are now deprecated (that is, don't use them any more because they are being phased out).

Generally, whenever GAUSS it comes across missing values, the program fails. This is so that missing values will not cascade through the program and cause erroneous results. However, in that case, none of the above code will work.

The way to get round this is to use

ENABLE;
DISABLE;

These two commands enable and disable checking for missing values. If GAUSS is ENABLEd, then any missing values will cause the program to crash. When GAUSS is DISABLEd, the checking is switched off and all the above operations with GAUSS can be carried out - along with the inclusion of missing values in calculations and the havoc that could wreak.

Whether to switch off missing value checking depends on the situation. If a missing value is not expected but would have a devastating effect on the program, then clearly GAUSS should be ENABLEd. Alternatively, if the program encounters lots of missing data which play no significant part in the results, then GAUSS should probably be DISABLEd. Intermediate cases require more thought. However, ENABLE and DISABLE can be used at any point, and so a program could DISABLE GAUSS while it checks for missing values and then ENABLE GAUSS again when it has dealt with them. There are no firm rules.

5 Other functions

GAUSS has a large repertoire of functions to perform operations on matrices. For most mathematical operations on or manipulations of a matrix (as opposed to altering the data) there will be a GAUSS function. Generally, these functions will be much faster than the equivalent user-written code.

To find a function, the GAUSS manuals have commands and operations organised into groups, as does the GAUSS Help system. In addition, each GAUSS function in the Command Reference will indicate what related functions are available.

[ previous page ] [ next page ]