Beaton’s Recipe for Least Squares
Vector 8:1, page 36
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 36
Introduction¶
An exhaustive search through academic papers on the subject of least-squares algorithms would no doubt reveal enough titles to fill the space allotted for this article! So why another one? The aim here is to look at the least-squares algorithm that replaces ‘domino’ in the generalized linear regression models volume (called ASLGREG) in ASL. As such, it updates the article by Sykes and Ansell (published in Vector Volume 1 Number 4, - see [3]) entitled ‘When domino is not sufficient’. As this title suggests, it looks at a particular implementation of Beaton’s least-squares operator, which in turn involves the notion of sweeping on successive predictor variables. For the sake of completeness, some aspects covered in that paper are also mentioned here.
Beaton’s Method¶
Householder’s least-squares operator (alias ‘domino’ to APLers) operates on the design matrix of a linear statistical model in which each column consists of the values of one predictor variable. In contrast, Beaton’s method manipulates a symmetric matrix of sums of squares and products formed from the predictor variables followed by the y variable, by performing successive sweeps (see [1]). The sweep on the ‘k’th variable is given by:
Where S denotes the sum of squares and products matrix. The values of S after the sweep on the kth variable all have appropriate statistical interpretations. In particular, having ‘swept’ variables 1,2,…k, the top-left k by k matrix is minus the covariance matrix of the regression coefficients; the first k elements in the last column (corresponding to the y variable) are the regression coefficients, and the bottom right-hand corner element is the residual sum of squares of the y variable after regression on variables 1 to k. Also note that repeating the same sweep twice is equivalent to not having swept on that variable at all, a feature useful for stepwise procedures.
As an example, consider the regression data used in our previous paper.
⎕←D←1,X1,X2,X3,[1.5]Y
1 2.5 1.05 0.625 2.6
1 2.375 1.96 3 1.3
1 3.125 1.7 1.75 2.5
1 1 0.87 1.375 1.5
1 1.875 1.3 1.75 3
1 1.5 1.73 3.125 2.2
1 0.875 1.68 3.5 2.2
1 0.125 1.04 2.5 2.4
1 3.5 1.47 0.875 1.7
1 1.625 1.57 2.625 2.2
Now form the sums of squares and products matrix:
⎕←S←(⍉D)+.×D
10 18.5 14.37 21.125 21.6
18.5 44.094 27.791 34.203 39.575
14.37 27.791 21.804 32.276 30.684
21.125 34.203 32.276 53.328 45.05
21.6 39.575 30.684 45.05 49.12
Now sweep on the first three columns, i.e. regress y on the constant vector, x1 and x2
3 SWEEP 2 SWEEP 1 SWEEP S
¯1.9034 0.040355 1.203 ¯4.885E¯15 2.6029
0.040355 ¯0.11619 0.1215 ¯0.8 ¯1.5763E¯3
1.203 0.1215 ¯0.99359 2.5 ¯0.30615
¯4.885E¯15 ¯0.8 2.5 1.5099E¯14 5.9952E¯15
2.6029 ¯1.5763E¯3 ¯0.30615 5.9952E¯15 2.3546
Hence the regression coefficients are given by 2.6029, ¯1.5763E¯3 and ¯0.30615, as verified by using ‘domino’:
Also, the residual sum of squares of the y data after regression is 2.3546, which may also be verified by the raw APL code:
the kth diagonal element (the pivot) plays a crucial role in the sweep operation since it is used to divide other quantities. A pivot close to zero may lead to considerable numerical inaccuracies. however, the sweep makes no distinction between the y+variable and other variables not so far swept on. Hence the pivot represents the residual sum of squares for that variable when regressed on the variables swept so far. For example, the 4 4 element of the result of the sweeps above, equal to 1.5099E¯14 shows that the variable x3 can be represented as a linear combination of the previous variables: in fact the figures in column 4 show that X3=(¯.8×X1)+2.5×X2.
Collinearity¶
When predictor variables are (near-) linearly dependent, as in the example above (a situation labelled collinearity), the use of domino either produces a DOMAIN ERROR or large regression coefficients which have huge variances and which are highly correlated. In simple terms, duplication of information confuses the algorithm, and results can sometimes owe more to the small details of the algorithm than the original problem.
A simple way out of this dilemma is to miss out the variable(s) which duplicate(s) information. This is where the pivot plays its role. Before pivoting on the kth variable, s[k;k] can be related to the correlation coefficient of the kth variable on those that have been pivoted. Hence it is possible to decide to not include the kth variable if 1 minus its r-squared value is, say, less than .001 - the criticality parameter.
Frane in his paper on checking the tolerance in regression calculations (see [2]) showed that prior to using the kth pivot, it is also easy to see whether, if the kth variable had been included previously, variables currently included would have failed the inclusion test. In either case, the kth variable is omitted.
Frane’s illustrative example uses the correlation matrix C of three predictor variables given below.
If the sweep operator is applied to variables 1 and 2, then the pivot for the third sweep is .01, which is above the criticality parameter .001.
However, if variables 2 and 3 were swept first, then the pivot for sweeping on the variable 1 is then less than .001.
The Basic Design¶
It is clear from the above discussion, that Beaton’s method can thus be applied to produce the covariance matrix of the estimates, the estimates, the residual sum of squares, and information on what variables are nearly linearly dependent on previously included variables. In addition, as the consecutive sweeps are performed, the residual sum of squares of the intermediate models can be ‘taken down and used in evidence’.
Hence the result of the weighted least squares routine ∆WLS takes the form of an n by [n+3] matrix m such that:
(a) M[⍳N;⍳N] is the covariance matrix of the estimates
(b) M[;N+1] is a column of the regression coefficients
(c) M[;N+2] is a column of the consecutive residual sum of squares
(d) M[;N+3] is a boolean vector with a zero value indicating that the corresponding variable has been omitted, for one of two possible reasons.
The function BEATON (also supplied in ASLGREG) illustrates this form of output. It operates, using a right argument equal to the sum of squares and products matrix, and a two-element vector left argument specifying the criticality parameter and the number of variables to be swept.
.001 3 BEATON S
¯1.9034 0.040355 1.203 ¯4.885E¯15 2.6029 1 2.464
0.040355 ¯0.11619 0.1215 ¯0.8 ¯1.5763E¯3 1 2.449
1.203 0.1215 ¯0.99359 2.5 ¯0.30615 1 2.3546
¯4.885E¯15 ¯0.8 2.5 1.5099E¯14 5.9952E¯15 0 0
.001 4 BEATON S
¯1.9034 0.040355 1.203 ¯4.885E¯15 2.6029 1 2.464
0.040355 ¯0.11619 0.1215 ¯0.8 ¯1.5763E¯3 1 2.449
1.203 0.1215 ¯0.99359 2.5 ¯0.30615 1 2.3546
¯4.885E¯15 ¯0.8 2.5 1.5099E¯14 5.9952E¯15 0 2.3546
Because of the collinearity already noted, these two statements provide nearly identical results. Note the succession of residual sum of squares, namely 2.464 (after fitting the constant term) 2.449 after additionally fitting x1, and so on. Applying BEATON to Frane’s correlation matrix of predictor variables gives the following result.
Thus the function BEATON incorporates a check for both possible occurrences of collinearity. As such it follows the methodology suggested by Frane, and also, according to the paper by Clarke (see [1]), that employed in the Royal Statistical Society’s generalized linear modelling package GLIM.
The Constant Vector¶
Most (but not necessarily all) regression models include an intercept term (the constant vector). Sweeping on it essentially converts the sums of squares and products matrix S to the equivalent one where each variable is replaced by the same variable centred about its mean. However, this conversion is subject to potentially large rounding errors, and should be avoided. (This is the essential reason why, in our previous paper, Beaton’s method was out-performed by, for example Bjorck’s method.) Hence to apply Beaton’s method successfully this factor must be taken into account.
The Final Recipe - ∆WLS¶
The result of applying ∆WLS has already been described above. But what are its parameters? First the left argument is a two-element numeric vector consisting of the criticality parameter and a boolean value, with 1 indicating that the constant vector is included, 0 that it is excluded. The effect of this flag, if set, is that the sums of squares and products matrix calculated is that centred around the mean. All regression results are then transformed back to standard form before delivery.
Secondly, the right argument is a segmented string indicating in order the names of the design matrix, the y-variable, (optionally) a prior-weight vector, and (optionally) a boolean vector indicating predictor variables not to be included. As an example illustrating its use consider the following quadratic regression problem.
DES←(⍳10)∘.*0 1 2
YV←(DES+.×1 ¯1 2)+¯50+?10⍴100
YV
21 16 60 64 49 27 108 113 175 233
YV⌹DES
49.483 ¯17.46 3.4583
.001 1 ∆WLS 'DES/YV/'
1.3833 ¯0.525 0.041667 49.483 1 45874
¯0.525 0.24129 ¯0.020833 ¯17.46 1 10926
0.041667 ¯0.020833 1.8939E¯3 3.4583 1 4611.6
Notice how the last column gives sufficient information to test that the quadratic term is necessary - this is exploited in ASLGREG by the provision of a TESTHYPOTHESIS function.
Now repeat the fit missing the last point out.
PW←(⍳10)<10
.001 1 ∆WLS 'DES/YV/PW/'
1.619 ¯0.67857 0.059524 39.833 1 22060
¯0.67857 0.34134 ¯0.032468 ¯11.173 1 6507.4
0.059524 ¯0.032468 3.2468E¯3 2.7273 1 4216.5
As a final example, consider the data used already, contained in the matrix D (with the y-variable in the last column).
.001 1 ∆WLS 'D[;1 2 3 4]/D[;5]/'
1.9034 ¯0.040355 ¯1.203 0 2.6029 1 2.464
¯0.040355 0.11619 ¯0.1215 0 ¯1.5763E¯3 1 2.449
¯1.203 ¯0.1215 0.99359 0 ¯0.30615 1 2.3546
0 0 0 0 0 0 2.3546
The penultimate column shows that variable x3 has not been included, because it is linearly dependent on x1 and x2.
Accuracy Considerations¶
How well does the proposed algorithm ∆WLS perform according to the benchmark introduced in our previous paper? The answer is: very well indeed. In that paper ([3]), Greenfield’s data (depending on two parameters n, and m) was analysed by different methods and the number of significant figures correct quoted. For n equal to 0 and m bigger than 4, domino gave domain error, and other methods gave varying numbers of significant figures.
The function GREENFIELD constructs x,y data, applies ∆WLS, and compares the answer with the true answer, outputting a 2 by 3 matrix such that each row consists of the estimate, true estimate and percentage relative error.
∇ R←GREENFIELD N_M
[1] ⍝ Creates Greenfield's data given N_M as two non -ve intgrs
[2] ⍝ Compares least-squares output from ∆WLS with true values
[3] ⍝ Outputs both and the relative % error to the nearest intgr
[4] ⍝ See Sykes and Ansell Vector Vol.1 No.4
[5] Y←7+ 1 2 3 4 5 ×10*-N_M[1]
[6] X←1,[1.5]8+ 1 3 2 4 5 ×10*-N_M[2]
[7] R← 1E¯3 1 ∆WLS 'X/Y/'
[8] R←R[;3],(7+(0.3×10*-N_M[1])-7.2×10*--/N_M),0.9×10*--/N_M
[9] R←⍉ 3 2 ⍴R,100×|(R[1 2]-R[3 4])÷R[1 2]
∇
GREENFIELD 0 5
¯719992.699978941 ¯719992.7 2.92493227896128E¯9
89999.9999973676 90000 2.92490262339789E¯9
GREENFIELD 0 13
¯72148849499738.5 ¯71999999999992.7 0.206308902744647
9018606187468.23 9000000000000 0.206308902744626
As m increases, the size of error increases. By m equal to 13, the previous methods used gave no correct significant figures. As examples such as the ones above show, ∆WLS has a relative error of at most .2% for the range of m values from 0 to 13.
Speed¶
How quick is ∆WLS, compared with domino? To make such a comparison is slightly unfair on ∆WLS as the latter gives so much more information than domino. It is annoying that Householder’s method has the potential to give all the information that ∆WLS does. It was such a consideration that prompted one of us to write a statistician’s form of domino - see [4]; nevertheless it is the obvious comparison to take.
The following table gives timings of domino and ∆WLS in milliseconds for a variety of design matrix sizes: the timings were made using Dyalog APL on an IBM Model 70.
+------------------------------------------------------------------+
| | Number of Columns |
| | 5 10 20 |
|----------------|-------------------------------------------------|
| | ⌹/∆WLS ⌹/∆WLS ⌹/∆WLS |
| | ------ ------ ------ |
| 50 | 60/270 160/660 550/1760 |
| Number of 100 | 110/440 330/880 1150/2470 |
| Rows 1000 | 1210/2750 3410/6480 11200/17800 |
+------------------------------------------------------------------+
As expected, the larger the problem, the more ∆WLS gains on domino, as the former only operates on columns to produce the sum of squares and products matrix; thereafter, the sweep operator is much faster than the equivalent Householder operation on a column. For small problems the ratio seems to be about 4:1, but for only medium sized problems (e.g. 100 by 20) this has dropped almost to 2:1.
Conclusion¶
Beaton’s recipe can be programmed efficiently in APL to provide a reliable, accurate, and flexible algorithm that is reasonably fast. It is superior to domino in coping with problems of collinearity. It is particularly gratifying to be able to analyse a particular model and to be able to extract residual sums of squares of intermediate models without further calculation.
References¶
[1] Clarke M.R.B. The Gauss-Jordan Sweep Operator with Detection of Collinearity. AS178, Applied Statistics Volume 31 No. 2, 1982, pp166-168.
[2] Frane J.W. A Note on Checking Tolerance in Matrix Inversion and Regression. Technometrics Volume 19, No. 4, November 1977, pp 513-514.
[3] Sykes A.M. and Ansell J.I. When Domino Is Not Sufficient. Vector Volume 1, no. 4, pp 155-160.
[4] Sykes A.M. A Second Generation Domino for Statisticians. Proceedings of APL87, Dallas, pp 369-380.