When Domino is not Sufficient

Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 155

Introduction¶

The use of DOMINO to solve linear equations is a particularly pleasant feature of APL to a Statistician as it forms a basis for Regression. However one should be on guard in using any general routine such as DOMINO. Whilst it performs adequately in many contexts, even in comparison to some NAG routines, there are circumstances where it causes concern. On these occasions it would be pleasant to have other general routines.

In this paper the authors outline two main deficiencies they have found using MIP’s APL in Regression and suggest remedies. The first problem is numerical stability due to variables in the regression having low variances. The second occurs when fitting many variables, some of which are collinear.

Briefly, the first problem may be overcome fairly simply by use of transformations. The second requires the introduction of forms of Generalised Matrix Inverse or Psuedo-Inverse. The latter leads to a more natural view of regression which facilitates extensions.

Regression¶

Regression is possibly the most widely used technique in Statistics. It attempts to account for the variance in one variable by another variable or set of variables. In most practical situations the underlying model is taken to be linear of the form;

where \(y\) is the variable one wishes to predict, \(x_1,x_2,\ldots,x_k\) is a set of variables used to predict \(y\), and \(e\) is the residual. The term residual is used to account for any portion of the variation in \(y\) unaccounted for by the variables \(x_1,\ldots,x_k\). Often \(e\) is taken to have a Normal distribution, though this is not a requirement as some practitioners seem to believe, see (1,11,12). (Other forms of generalisation are also mentioned in these papers.) The obvious APL solution for the unknown parameters \(B\) would be;

Y⌹X

where \(X\) is the design matrix with columns composed of \(x_1,\ldots,x_k\) with the first column being a vector of 1’s.

The mathematical representation of the solution is not so attractive being;

Hence the main problems revolve around the invertibility of \(X^TX\). One can divide the problem into two possible classes; Near Singular and Singular.

Near Singular of X^TX due to small variances.¶

A simple example, given by Greenfield and Siday (8) shows how small variances cause many statistical packages to fail. The data are presented in Table 1. Applying DOMINO results in a DOMAIN ERROR for even moderate values of \(n\) and \(m\), see Table 2. Essentially the routine cannot cope with the variance of the predictive variable. A simple solution to the problem is to strip off the mean and then appropriately scale the variable involved. Fitting can then usually be carried out using DOMINO. After fitting one can then rescale if desired.

Table 1 — Greenfield’s Data

(For this data intercept should be \(7+0.3\times10^{-n} - 7.2\times10^{m-n}\), and the gradient should be \(0.9\times10^{m-n}\).)

Obviously there is a wide choice of possible scaling procedures. One might scale by the standard deviation though this would entail squaring of already small values. An alternative is to use average of the absolute deviations from the mean, this is the method employed in the algorithm INVERT.

    ∇ BETA←Y INVERT X;MY;Y1;X1;X2;N;M
[1]    Y1←Y-MY←+/Y÷⍴Y
[2]    X1←X1÷N⍴X2←+⌿|(X1←X-N⍴M←+⌿X÷N[1])÷(N←⍴X)[1]
[3]    N←1↓BETA←Y1⌹1,X1
[4]    BETA←(MY+BETA[1]-+/M×N),N←N÷X2
    ∇

The algorithm has been fairly successful in practice, the result of analysing data in Table 1 is given in Table 2.

Table 2

Number of correct significant figures when applying DOMINO, INVERT and BJORCK to Greenfield’s Data with \(n=0\) and \(m=1,2,\ldots,13\).

(DE means DOMAIN ERROR and ALL means correct answer given.)

Singularity due to Collinearity.¶

Problems due to collinearity are more complex and result in the singularity of \(X^TX\), due to linear dependence between the columns of \(X\). This simply means that one or more of the variables can be expressed in terms of the others for the data in question. This may or may not be due to an intrinsic relationship between the variables. Hence although DOMINO will fail there exists a solution though it will not be unique. It is achieved by removing a dependent subset. The only question to a Statistician in this context would be which variable or set of variables should be removed. The solution depends partially on the context and criteria of efficiency of choice, i.e. in some senses a matter of taste.

The main insight into the Generalized Inverse of a Matrix came from Penrose (13) who reintroduced the concept initially proposed by Moore (10). He was subsequently followed in Statistics by Rao (15) and Wilkinson (18), and summarized both in Chambers (5), and Pringle and Rayner (14). There is a need therefore for an efficient Generalised Matrix Inverse. There are numerous candidates for the role. In this paper we restrict attention primarily to two, but interested readers may consider others, e.g. Given’s (6), Householder’s (4,9) and Rogers and Wilkinson (16) are worthy of study. The two selected are due to Beaton (2) and Bjorck (3). Each of these can be said to have particular merits; Beaton’s is compact; Bjorck is preferred on grounds of accuracy.

Beaton’s Algorithm¶

This routine works directly with the \(X^TX\) matrix, and hence is compact. ( \(X\) is an \(m\times n\) matrix with \(m\) usually much larger than \(n\), so \(X^TX\) is then an \(n\times n\) matrix.) Unfortunately this ‘squaring’ can exacerbate numerical instability when near singular matrices are encountered. Obviously its compactness makes it more appropriate for use when storage or workspace size is limited, or other constraints exist on the size of individual components.

The algorithm applies a sweep-based routine to \(X^TX\), which in the context of regression is either the sum of products matrix or the variance-covariance matrix. Given one has selected a particular pivot, a suitable diagonal element, the pivot is used to remove the effect of the associated variable from the matrix. In the case of using the variance-covariance matrix the pivot is the variance of the variable being fitted, and a matrix of the residual variance-covariances results.

The APL code for a single pivot operation is given below and the matrix is inverted by cycling through all the pivotal elements.

    ∇ M←N BEATON K
[1]    M←N-N[K;]∘.×N[K;]÷N[K;K]
[2]    M[K;]←N[K;]÷N[K;K]
[3]    M[;K]←N[;K]÷N[K;K]
[4]    M[K;K]←-÷N[K;K]
    ∇

In the context of regression it is helpful to consider the following matrix since application of Beaton to it has some pleasant attributes.

At each pivot-sweep the effect of a given variable is removed not only from the other \(X\) variables but also from \(Y\) and hence the side-effects are as follows:

(1) The regression coefficients are produced in the last column.

(2) Dropping off the last row and column produces the \(-(X^TX)\) and this yields the estimators variance covariance matrix.

(3) The entry in the bottom righthand corner is the residual sum of squares in the usual Normal Theory of regression.

What happens when a pivot is 0? Each pivot is examined before it enters the routine. If below a given threshold, set by the user, then that variable is marked. A marked variable is one which is (almost) totally dependent on previous variables, it therefore affords no additional explanation. For the data presented in Table 3 DOMINO would produce DOMAIN ERROR whereas use of Beaton yields the following solution;

Obviously it is not unique — an alternative solution is

In the penultimate section of the paper the authors give some indication of criteria for selecting an appropriate model.

Table 3 — Data with Linear Dependence amongst the X variables

Bjorck’s Algorithm¶

This is a modified version of the Gram-Schmidt Orthogonalisation procedure. It might be thought of as a sweep procedure since each column sweeps through the matrix \(X\) in turn. Hence one removes the effect of a variable from all subsequent variables. It does however result in two matrices usually referred to as \(Q\) and \(R\), hence the generic name for this type of procedure — a \(Q-R\) technique. \(R\) is a square upper triangular matrix whilst \(Q\) is an orthogonalization of \(X\).

The APL code for this algorithm is given below:

    ∇ B←Y BJORCK X;N;J;K;T;S
[1]    J←=/⍴R←(N,N←¯1↑⍴Q←(⍴X)⍴0)⍴0
[2]   L:R[J;J]←(X[;J]+.×X[;J])*0.5
[3]    Q[;J]←X[;J]÷R[J;J]
[4]    →(J=N)/R
[5]    R[J;K]←Q[;J]+.×X[;K←J+⍳N-J]
[6]    X[;K]←X[;K]-Q[;J]∘.×R[J;K]
[7]    →(1<J←J+1)/L
[8]   R:C←Y+.×Q
[9]    B←S\C[T]⌹R[T;T←(S←1E¯9≤ 1 1 ⍉R)/⍳N]
    ∇

Other desired quantities may be calculated easily. For example the residual sum of squares is obtained by

+/(Y-Q+.×C)*2

Chambers indicates in (5) that it is more accurate to use the above approach to calculate the residual than using the estimates of Y themselves.

The properties of Bjorck’s algorithm have been extensively investigated by previous authors (5,16) and found to be very well behaved. As an indication Bjorck’s algorithm was applied to Greenfield’s data and the results presented in Table 2.

A More Natural View of Regression¶

The traditional view of regression is through the solution of the Normal Equations, see (7,17). Whilst mathematically sound we have found it has little intuitive appeal even to mathematics students. The idea of ‘sweeping’ in both Beaton’s and Bjorck’s algorithms suggests the geometric viewpoint of projecting a new variable onto the present variable space. Also such procedures imply that the order of introducing covariates may be important.

Choosing the Solution in case of Singularity¶

As indicated before, where no unique solution exists the selection of variables to include is a matter for the individual. First of all, the variables included in the model should be understandable within the physical background of the problem. Secondly, models should be as simple as possible, because they are easier to understand and tend to be more robust. Hence one should include as few variables as possible. Also the variables should themselves be uncomplicated in form.

The selection of variables in these two techniques depends on the order of their consideration, (selection of order of sweeping). This means either the order has to be specified by the user, or some criterion must be enforced by the routine. In Statistics it would be controversial to lay down a definitive solution to this problem. However, the authors feel that the selection based on partial correlation, while not always optimum, might be a suitable way around the problem.

Extensions of the Problem¶

The paper has discussed some numerical problems associated with standard Regression Analysis, where traditionally one assumes (amongst other things) that the error is Normal and has constant variance. The latter may be relaxed by using weighted least squares, replacing \(X\) by \(WX\) and \(Y\) by \(WY\) where \(W\) is a diagnonal matrix of weights. The amended algorithm forms the core of the iterative solution of such Generalized Linear Models which allows different error distributions to be contemplated, see (11,12). The authors are at present working on an APL version of the Generalized Linear Models (1).

References¶

(1) Ansell, J.I., and Sykes, A.M., (1984), ‘APLIM — Applied Linear Interactive Models’, Tech. Rep. 84/2, Department of Man. Sc. and Stats, University College of Swansea, Wales.

(2) Beaton, A.E., (1964), ‘The use of special matrix operators in statistical calculus’, Research Bulletin RB-64-57 Educational Testing Service, Princeton, New Jersey.

(3) Bjorck, A., (1967), ‘Solving linear least squares problems using Gram-Schmidtt orthogonalization’, BIT, 7, 257-78.

(4) Bussinger, P.A., and Golub, G.H., (1965), ‘Linear least squares solution by Householder transformation’, Numer. Math., 7, 269-276.

(5) Chambers, J.M., (1977), Computational methods for data analysis, Wiley, New York.

(6) Clarke, M.R.B., (1981), ‘Algorithm AS163: A Givens algorithm for moving from one linear model to another without going back to the data’, Appl. Stat., 30, 166-8.

(7) Edwards, A.L., (1979), Multiple regression and the analysis of variance and covariance, W.H. Freeman & Co., San Francisco.

(8) Greenfield, T., and Siday, S., (1980), ‘Statistical computing for business and industry’, Statistician, 29, 1, 33-55.

(9) Lawson, C.L., and Hanson, R.J., (1974), Linear least squares solutions, Prentice-Hall, Englewood Cliffs, NJ.

(10) Moore, E.H., (1935), ‘General analysis’, Mem. Amer. Phil. Soc., 1, 197.

(11) Nelder, J.A., and McCullagh, P., (1983), Generalized Linear Models, Methuen, London.

(12) Nelder, J.A., and Wedderburn, R.W.M., (1972), ‘Generalized Linear Models’, JRSS, A, 135, 370-84.

(13) Penrose, R., (1955), ‘A generalized inverse of matrices’, Proc. Camb. Phil. Soc., 51, 406-11.

(14) Pringle, R.M., and Rayner, A.A., (1971), Generalized inverse matrices with applications to statistics, Griffin, London.

(15) Rao, C.R., (1962), ‘A note on a generalized inverse of a matrix with applications to problems in mathematical statistics’, JRSS, B, 24, 152-8.

(16) Rogers, C.E., and Wilkinson, G.N., (1974), ‘Regression, curve fitting and smoothing numerical problems in recursive analysis of variance algorithms’, IMA Bull., June, 141-143.

(17) Stoodley, K.D.C., Lewis, T., and Stainton, C.L.S., (1980), Applied Statistical Techniques, Ellis Horwood Ltd, Chichester.

(18) Wilkinson, G.N., (1958), ‘Estimation of missing values for the analysis of incomplete data’, Biometrics, 14, 257-86.