Parallel backfits

Variance Structure

For now, I’ve restricted parallel backfits to the case where residual variance is assumed to be constant across all entries (var.type = 0). Whereas sequential backfitting can take advantage of the fact that the update to expected residuals is rank-one, parallel backfits must re-estimate the residual variance from scratch at each iteration. Recall that the most useful variance structures (row-wise, column-wise, and constant) can be estimated as simple functions of expected squared residuals (row-wise means, column-wise means, and the overall mean). Recall also that flashier doesn’t usually store a full matrix of residuals \(R\), so that the expected squared residual \(R_{ij}^2\) must be calculated as:

\[ \mathbb{E} R_{ij}^2 = \mathbb{E} (Y_{ij} - \sum_k L_{ik} F_{jk})^2 = Y_{ij}^2 - 2 Y_{ij} \sum_k \mathbb{E} L_{ik} \mathbb{E} F_{jk} + \sum_{k \ne \ell} \mathbb{E} L_{ik} \mathbb{E} F_{jk} \mathbb{E} L_{i\ell} \mathbb{E} F_{j\ell} + \sum_k \mathbb{E} L_{ik}^2 \mathbb{E} F_{jk}^2 \]

When residual variance is constant across all entries, we only need \(\sum_{i, j} R_{ij}^2\), and each of the above terms can be efficiently summed over \(i\) and \(j\). The trick, of course, is to move the summation over \(i\) and \(j\) to the inside (and to pre-compute \(\sum_{i, j} Y_{ij}^2\)). For example,

\[ \sum_{i, j} \sum_{k \ne \ell} \mathbb{E} L_{ik} \mathbb{E} F_{jk} \mathbb{E} L_{i\ell} \mathbb{E} F_{j\ell} = \sum_{k, \ell} \sum_i \mathbb{E} L_{ik} \mathbb{E} L_{i\ell} \sum_j \mathbb{E} F_{jk} \mathbb{E} F_{j\ell} - \sum_k \sum_i (\mathbb{E} L_{ik})^2 \sum_j (\mathbb{E} F_{jk})^2 \]

The first term on the RHS can be computed as sum(crossprod(EL) * crossprod(EF)); the second can be computed as sum(colSums(EL^2) * colSums(EF^2)).

For row-wise or column-wise variance structures, however, the first term is much more difficult to compute. Instead of simply taking crossproducts, one must form a \(n \times k^2\) (or \(p \times k^2\)) matrix, so that unless \(k^2 \ll n\) (or \(k^2 \ll p\)), one would not be much worse off by simply storing the matrix of expected residuals. But we only stand to benefit from parallelization when we are doing large backfits on large data matrices; that is, when \(k\) is not small and when storing a matrix of residuals is expensive.

Fixed factors

Session information

sessionInfo()

#> R version 3.5.3 (2019-03-11)
#> Platform: x86_64-apple-darwin15.6.0 (64-bit)
#> Running under: macOS Mojave 10.14.5
#> 
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_3.2.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.1        knitr_1.22        whisker_0.3-2    
#>  [4] magrittr_1.5      workflowr_1.2.0   tidyselect_0.2.5 
#>  [7] munsell_0.5.0     colorspace_1.4-1  R6_2.4.0         
#> [10] rlang_0.3.1       dplyr_0.8.0.1     stringr_1.4.0    
#> [13] tools_3.5.3       grid_3.5.3        data.table_1.12.2
#> [16] gtable_0.3.0      xfun_0.6          withr_2.1.2      
#> [19] git2r_0.25.2      htmltools_0.3.6   assertthat_0.2.1 
#> [22] yaml_2.2.0        lazyeval_0.2.2    rprojroot_1.3-2  
#> [25] digest_0.6.18     tibble_2.1.1      crayon_1.3.4     
#> [28] purrr_0.3.2       fs_1.2.7          glue_1.3.1       
#> [31] evaluate_0.13     rmarkdown_1.12    labeling_0.3     
#> [34] stringi_1.4.3     pillar_1.3.1      compiler_3.5.3   
#> [37] scales_1.0.0      backports_1.1.3   pkgconfig_2.0.2