2021-07-05

To serve as a demonstration and test of modelling using the *largeScaleR* system, the canonical statistical model, the linear model, was chosen for implementation. While simple to derive the mathematics detailing a linear model, the naive method doesn’t map straightforwardly to a distributed system. This is due to matrix inversion and multiplication operations requiring more than isolated chunkwise information. Conveniently, the common method of decomposing a matrix, typically for computationally simplified inversion, is given in a chunked form as part of Algorithm AS274 [1]. This algorithm is in turn implemented and wrapped through the *biglm* package in R with several additional features such as such as sandwiching [2]. The *biglm* package creates a linear model object from some initial row-wise chunk of data with a space complexity of *m**a**t**h**c**a**l**O*(*p*^{2}) for *p* variables. The linear model object is then updated sequentially with further chunks of the data, until all of the data has been read, yielding a final linear model object with methods available for standard linear model object inspection, including summaries and predictions. This sequential rolling update over chunked data is able to be captured succinctly in a `Reduce()`

pattern typical of the functional programming paradigm. Using the dataset as an example, a linear model of the form, *P**e**t**a**l*.*L**e**n**g**t**h*_{i} = *S**e**p**a**l*.*W**i**d**t**h*_{i} + *S**e**p**a**l*.*L**e**n**g**t**h*_{i} + *S**e**p**a**l*.*W**i**d**t**h* × *S**e**p**a**l*.*L**e**n**g**t**h* + *ϵ*_{i} is fitted using *biglm* as part of a demonstrative non-distributed `Reduce()`

in lst. 1.

Listing 1: Splitting the iris dataframe into 15 chunks stored as elements of a list and reducing over the list with the biglm update function.

```
library(biglm)
<- 15
nc <- split(iris, cumsum((seq(nrow(iris)) - 1) %% nc == 0)) # split to chunks
cs <- Reduce(f = update, x = cs[-1],
model init = biglm(Petal.Length ~ Sepal.Width * Sepal.Length, cs[[1]]))
```

The *largeScaleR* package provides the distributed Reduce, `dreduce`

, as discussed in the dreduce document. The `dreduce`

is therefore to be used as an equivalent structural backbone for the implementation of a distributed linear model.

The implementation of the distributed linear model is succinct enough to be given in it’s entireity in lst. 2.

Listing 2: Full listing of distributed linear model implementation.

```
<- function(formula, data, weights=NULL, sandwich=FALSE) {
dlm stopifnot(largeScaleR::is.distObjRef(data))
<- largeScaleR::chunkRef(data)
chunks stopifnot(length(chunks) > 0L)
<- dbiglm(formula, chunks[[1]], weights, sandwich)
dblm if (length(chunks) != 1L)
::dreduce(f="biglm::update.biglm",
largeScaleRx=largeScaleR::distObjRef(chunks[-1]),
init=dblm)
else init
}
<- function(formula, data, weights=NULL, sandwich=FALSE) {
dbiglm stopifnot(largeScaleR::is.chunkRef(data))
<- largeScaleR::currCallFun(-1)
sys.call ::do.ccall(what="biglm::biglm",
largeScaleRargs=list(formula=largeScaleR::envBase(formula),
data=data,
weights=largeScaleR::envBase(weights),
sandwich=sandwich),
target=data,
insert=list(sys.call=largeScaleR::envBase(sys.call)))
}
```

This implementation includes several important aspects. Mapping to the list of chunks, an existing distributed object must be transformed to a list of chunks for the reduce function. The reduce function itself is accessed through a light wrapper `dreduce()`

, with all of the internal code operating transparently on chunks without concern for the type of object. A major divergence is given in the generation of the initial reduction object. This uses a `do.ccall`

function to create the initial biglm linear model object, which is then passed to the `dreduce()`

function.

A more significant divergence than different initialisation, with severe performance implications, is shown through the intercept and insertion of the function call. This is made necessary by the fact that *biglm* captures the call it was called with, and stores it as part of the model object – this is not unique to *biglm*, and this behaviour is common to most modelling functions in R. Due to actually enacting the call to construct the inital *biglm* linear model object on the worker process, rather than the master process, with worker processes evaluating requests through construction with the `do.call()`

function, rather than exact replication of the initial call, the call as seen by the function is not necessarily the same as that issued on the master process. This presents two problems: inaccuracy, and unbounded call sizes. Inaccuracy is not an enormous problem, as the call isn’t typically used for anything other than rendering a portion of the string representation of the model object. The greater issue is that as calls are constructed on the worker, all arguments are evaluated, and the captured call will include fully expanded objects in a `dump()`

-like form. This object may very well have a larger memory footprint than all of the arguments to the call combined, and will lead to memory limitations and slowdowns, particularly when transferring the model object from process to process.

The solution to this is twofold: allow new insertion environments to be inserted to the requested function on the worker for the purpose of non-destructive masking; and capture the call on the master end, wrapped into a function returning the call, and insert into the new insertion environment to take the place of the previous call capture function.

An insertion environment follows the simple concept of being placed between some function and the function’s original enclosing environment for the purpose of having it’s objects first on the search path, for masking or perhaps making previously global variables available. Its form is given in the diagram at Figure .

The call is captured on the master end, with a function constructed to return this call, given by the *largeScaleR*-provided `currCallFun()`

, and this is then sent to the worker to be inserted as `sys.call()`

in the insertion environment, thus effectively masking the call capture function at the top level of the requested *biglm* function. Depending on perspective, the fact that this only works at the top level can be a feature, as it doesn’t messily mask further along the call stack, however it has the associated limitation.

Worth noting is that call capture is notoriously messy, with *biglm* itself featuring source code directly manipulating the call as part of some functions.

[1]

A. J. Miller, “Algorithm AS 274: Least squares routines to supplement those of gentleman,” *Journal of the Royal Statistical Society. Series C (Applied Statistics)*, vol. 41, no. 2, pp. 458–478, 1992,Available: http://www.jstor.org/stable/2347583

[2]

T. Lumley, *Biglm: Bounded memory linear and generalized linear models*. 2013.