Solutions to Multilevel Sparse Matrix Problems

TL;DR

This paper addresses multilevel sparse matrix problems common in hierarchical data analysis, providing solutions that enable efficient computation of standard errors and Bayesian inference for complex models.

Contribution

It offers novel matrix inverse sub-block results and comprehensive solutions for two- and three-level problems, serving as blueprints for higher-level cases.

Findings

Streamlined computation of standard errors in multilevel models

Novel matrix inverse sub-block results for efficient inference

Blueprints for solving higher-level multilevel matrix problems

Abstract

We define and solve classes of sparse matrix problems that arise in multilevel modeling and data analysis. The classes are indexed by the number of nested units, with two-level problems corresponding to the common situation in which data on level 1 units are grouped within a two-level structure. We provide full solutions for two-level and three-level problems and their derivations provide blueprints for the challenging, albeit rarer in applications, higher level versions of the problem. Whilst our linear system solutions are a concise recasting of existing results, our matrix inverse sub-block results are novel and facilitate streamlined computation of standard errors in frequentist inference as well as allowing streamlined mean field variational Bayesian inference for models containing higher level random effects.

Tables

Table 1. Table 1: Averages (standard deviations) of elapsed computing times in seconds for solving the two-level sparse matrix problem naïvely versus with the streamlined approach provided by Theorem 2 . The fourth column lists the ratios of the median computing times.

$m$	naïve	streamlined	naïve/streamlined
100	0.103 (0.0173)	0.00253 (0.000541)	40.8
200	0.768 (0.0214)	0.00429 (0.000574)	179.0
400	6.050 (0.1110)	0.00798 (0.000586)	758.0
800	48.300 (0.7830)	0.01530 (0.001650)	3150.0
1600	386.000 (6.0200)	0.03000 (0.001900)	12900.0

Equations

$$\bm{A}=\left[\begin{array}[]{c | c | c | c | c}\bm{A}_{11}&\bm{A}_{12,1}&\bm{A}_{12,2}&\ \ \dots&\bm{A}_{12,m}\\ \hline\cr\bm{A}_{12,1}^{T}&\bm{A}_{22,1}&\bm{O}&\dots&\bm{O}\\ \hline\cr\bm{A}_{12,2}^{T}&\bm{O}&\bm{A}_{22,2}&\dots&\bm{O}\\ \hline\cr\vdots&\vdots&\vdots&\ddots&\vdots\\ \hline\cr\bm{A}_{12,m}^{T}&\bm{O}&\bm{O}&\dots&\bm{A}_{22,m}\end{array}\right]$$

$${\mathbf{A}}_{11}\text{is}p\times p\text{and, for each 1≤i≤m,}{\mathbf{A}}_{12,i}\text{is}p\times q\text{and}{\mathbf{A}}_{22,i}\text{is}q\times q.$$

$$\mathbf{A}\mathbf{x}=\mathbf{a},$$

$$\bm{A}^{-1}=\left[\begin{array}[]{c | c | c | c | c}\bm{A}^{11}&\bm{A}^{12,1}&\bm{A}^{12,2}&\ \ \dots&\bm{A}^{12,m}\\ \hline\cr\bm{A}^{12,1T}&\bm{A}^{22,1}&{\LARGE\mbox{$\times$}}&\dots&{\LARGE\mbox{$\times$}}\\ \hline\cr\bm{A}^{12,2T}&{\LARGE\mbox{$\times$}}&\bm{A}^{22,2}&\dots&{\LARGE\mbox{$\times$}}\\ \hline\cr\vdots&\vdots&\vdots&\ddots&\vdots\\ \hline\cr\bm{A}^{12,mT}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&\dots&\bm{A}^{22,m}\end{array}\right],\qquad\bm{a}\equiv\left[\begin{array}[]{c}\bm{a}_{1}\\ \hline\cr\bm{a}_{2,1}\\ \hline\cr\bm{a}_{2,2}\\ \hline\cr\vdots\\ \hline\cr\bm{a}_{2,m}\end{array}\right]\qquad\text{and}\qquad\bm{x}\equiv\left[\begin{array}[]{c}\bm{x}_{1}\\ \hline\cr\bm{x}_{2,1}\\ \hline\cr\bm{x}_{2,2}\\ \hline\cr\vdots\\ \hline\cr\bm{x}_{2,m}\end{array}\right]$$

$$\mbox{both}\ \bm{a}_{1}\ \mbox{and}\ \bm{x}_{1}\ \mbox{are}\ p\times 1\ \mbox{and, for $1\leq i\leq m$,}\ \mbox{both}\ \bm{a}_{2,i}\ \mbox{and}\ \bm{x}_{2,i}\ \mbox{are}\ q\times 1.$$

$$\mbox{obtain the matrix}\ \bm{A}^{11}\ \mbox{and, for each $1\leq i\leq m$, obtain the matrices}\ \bm{A}^{12,i}\ \mbox{and}\ \bm{A}^{22,i}.$$

$${\mathbf{A}}^{11}$$

$$\text{and}{\mathbf{A}}^{12,i}$$

$$|\bm{A}|=\left|\left(\bm{A}^{11}\right)^{-1}\right|\,\prod_{i=1}^{m}\big{|}\bm{A}_{22,i}\big{|}.$$

$${\mathbf{x}}_1={\mathbf{A}}^{11}\left({\mathbf{a}}_1-\sum _{i=1}^m{\mathbf{A}}_{12,i}{\mathbf{A}}_{22,i}^{-1}{\mathbf{a}}_{2,i}\right)\text{and}{\mathbf{x}}_{2,i}={\mathbf{A}}_{22,i}^{-1}({\mathbf{a}}_{2,i}-{\mathbf{A}}_{12,i}^T{\mathbf{x}}_1),1\le i\le m.$$

$$\mathbf{X}=\mathbf{Q}\left[\begin{array}{c}\mathbf{R}\\ \mathbf{O}\end{array}\right],$$

$$\Vert \mathbf{b}-\mathbf{B}\mathbf{x}{\Vert }^2\equiv (\mathbf{b}-\mathbf{B}\mathbf{x}{)}^T(\mathbf{b}-\mathbf{B}\mathbf{x})$$

$$\bm{B}\equiv\left[\begin{array}[]{c | c | c | c | c}\bm{B}_{1}&\accentset{\mbox{\mbox{{$\bullet$}}}}{\bm{B}}_{1}&\bm{O}&\dots&\bm{O}\\ \hline\cr\bm{B}_{2}&\bm{O}&\accentset{\mbox{\mbox{{$\bullet$}}}}{\bm{B}}_{2}&\dots&\bm{O}\\ \hline\cr\vdots&\vdots&\vdots&\ddots&\vdots\\ \hline\cr\bm{B}_{m}&\bm{O}&\bm{O}&\dots&\accentset{\mbox{\mbox{{$\bullet$}}}}{\bm{B}}_{m}\end{array}\right]\quad\text{and}\quad\bm{b}\equiv\left[\begin{array}[]{c}\bm{b}_{1}\\ \hline\cr\bm{b}_{2}\\ \hline\cr\vdots\\ \hline\cr\bm{b}_{m}\end{array}\right]$$

$$\bm{B}_{i}\ \mbox{is}\ n_{i}\times p,\quad\accentset{\mbox{\mbox{{$\bullet$}}}}{\bm{B}}_{i}\ \mbox{is}\ n_{i}\times q\quad\mbox{and}\quad\bm{b}_{i}\ \mbox{is}\ n_{i}\times 1\quad\mbox{for $1\leq i\leq m$}$$

$$\mathbf{A}\mathbf{x}=\mathbf{a}\text{mbox}where\mathbf{A}={\mathbf{B}}^T\mathbf{B}\text{mbox}and\mathbf{a}={\mathbf{B}}^T\mathbf{b}.$$

$${\mathbf{A}}_{11}=\sum _{i=1}^m{\mathbf{B}}_i^T{\mathbf{B}}_i,{\mathbf{a}}_1=\sum _{i=1}^m{\mathbf{B}}_i^T{\mathbf{b}}_i$$

$$\bm{A}_{12,i}=\bm{B}_{i}^{T}\accentset{\mbox{\mbox{{$\bullet$}}}}{\bm{B}}_{i},\quad\bm{A}_{22,i}=\accentset{\mbox{\mbox{{$\bullet$}}}}{\bm{B}}_{i}^{T}\accentset{\mbox{\mbox{{$\bullet$}}}}{\bm{B}}_{i},\quad\bm{a}_{2,i}=\accentset{\mbox{\mbox{{$\bullet$}}}}{\bm{B}}_{i}^{T}\bm{b}_{i},\quad 1\leq i\leq m.$$

$$\operatorname*{stack}\limits_{1\leq i\leq d}(\bm{M}_{i})\equiv\left[\begin{array}[]{c}\bm{M}_{1}\\ \vdots\\ \bm{M}_{d}\end{array}\right].$$

$${\mathbf{c}}_{0i}\equiv {\mathbf{Q}}_i^T{\mathbf{b}}_i,{\mathbf{c}}_{1i}\equiv \text{first q rows of}{\mathbf{c}}_{0i},{\mathbf{c}}_{2i}\equiv \text{remaining rows of}{\mathbf{c}}_{0i},$$

$${\mathbf{C}}_{0i}\equiv {\mathbf{Q}}_i^T{\mathbf{B}}_i,{\mathbf{C}}_{1i}\equiv \text{first q rows of}{\mathbf{C}}_{0i}\text{and}{\mathbf{C}}_{2i}\equiv \text{remaining rows of}{\mathbf{C}}_{0i}.$$

$$\mathbf{c}\equiv \text{first p rows of}{\mathbf{Q}}^T\left\{\underset{1\le i\le m}{\mathrm{stack}}({\mathbf{c}}_{2i})\right\}.$$

$${\mathbf{x}}_1={\mathbf{R}}^{-1}\mathbf{c},{\mathbf{A}}^{11}={\mathbf{R}}^{-1}{\mathbf{R}}^{-T}$$

$${\mathbf{x}}_{2,i}={\mathbf{R}}_i^{-1}({\mathbf{c}}_{1i}-{\mathbf{C}}_{1i}{\mathbf{x}}_1),{\mathbf{A}}^{12,i}=-{\mathbf{A}}^{11}({\mathbf{R}}_i^{-1}{\mathbf{C}}_{1i}{)}^T,{\mathbf{A}}^{22,i}={\mathbf{R}}_i^{-1}({\mathbf{R}}_i^{-T}-{\mathbf{C}}_{1i}{\mathbf{A}}^{12,i}).$$

$$|\bm{A}|=\left\{\big{(}\mbox{product of the diagonal entries of $\bm{R}$}\big{)}\,\prod_{i=1}^{m}\,\big{(}\mbox{product of the diagonal entries of $\bm{R}_{i}$}\big{)}\right\}^{2}.$$

$$\bm{A}=\left[\begin{array}[]{c | c | c | c | c | c | c | c}\bm{A}_{11}&\bm{A}_{12,1}&\bm{A}_{12,11}&\bm{A}_{12,12}&\bm{A}_{12,2}&\bm{A}_{12,21}&\bm{A}_{12,22}&\bm{A}_{12,23}\\ \hline\cr\bm{A}_{12,1}^{T}&\bm{A}_{22,1}&\bm{A}_{12,1,1}&\bm{A}_{12,1,2}&\bm{O}&\bm{O}&\bm{O}&\bm{O}\\ \hline\cr\bm{A}_{12,11}^{T}&\bm{A}_{12,1,1}^{T}&\bm{A}_{22,11}&\bm{O}&\bm{O}&\bm{O}&\bm{O}&\bm{O}\\ \hline\cr\bm{A}_{12,12}^{T}&\bm{A}_{12,1,2}^{T}&\bm{O}&\bm{A}_{22,12}&\bm{O}&\bm{O}&\bm{O}&\bm{O}\\ \hline\cr\bm{A}_{12,2}^{T}&\bm{O}&\bm{O}&\bm{O}&\bm{A}_{22,2}&\bm{A}_{12,2,1}&\bm{A}_{12,2,2}&\bm{A}_{12,2,3}\\ \hline\cr\bm{A}_{12,21}^{T}&\bm{O}&\bm{O}&\bm{O}&\bm{A}_{12,2,1}^{T}&\bm{A}_{22,21}&\bm{O}&\bm{O}\\ \hline\cr\bm{A}_{12,22}^{T}&\bm{O}&\bm{O}&\bm{O}&\bm{A}_{12,2,2}^{T}&\bm{O}&\bm{A}_{22,22}&\bm{O}\\ \hline\cr\bm{A}_{12,23}^{T}&\bm{O}&\bm{O}&\bm{O}&\bm{A}_{12,2,3}^{T}&\bm{O}&\bm{O}&\bm{A}_{22,23}\end{array}\right].$$

$$\bm{A}=\left[\begin{array}[]{cc}\bm{A}_{11}&\bm{A}_{12}\\[4.30554pt] \bm{A}_{12}^{T}&\bm{A}_{22}\end{array}\right]$$

$$\bm{A}_{12}=\left\{\left[\begin{array}[]{c| c|c|c}\bm{A}_{12,i}&\bm{A}_{12,ij}&\dots&\bm{A}_{12,in_{i}}\end{array}\right]\right\}_{1\leq i\leq m}$$

$$\bm{A}_{22}=\operatorname*{blockdiag}_{1\leq i\leq m}\left(\left[\begin{array}[]{c | c | c | c | c}\bm{A}_{22,i}&\bm{A}_{12,\,i,\,1}&\bm{A}_{12,\,i,\,2}&\ \ \dots&\bm{A}_{12,\,i,\,n_{i}}\\ \hline\cr\bm{A}_{12,\,i,\,1}^{T}&\bm{A}_{22,i1}&\bm{O}&\dots&\bm{O}\\ \hline\cr\bm{A}_{12,\,i,\,2}^{T}&\bm{O}&\bm{A}_{22,i2}&\dots&\bm{O}\\ \hline\cr\vdots&\vdots&\vdots&\ddots&\vdots\\ \hline\cr\bm{A}_{12,\,i,\,n_{i}}^{T}&\bm{O}&\bm{O}&\dots&\bm{A}_{22,in_{i}}\end{array}\right]\right)$$

$$\mathbf{A}\mathbf{x}=\mathbf{a},$$

$$\bm{A}^{-1}=\left[\begin{array}[]{c | c | c | c | c | c | c | c}\bm{A}^{11}&\bm{A}^{12,1}&\bm{A}^{12,11}&\bm{A}^{12,12}&\bm{A}^{12,2}&\bm{A}^{12,21}&\bm{A}^{12,22}&\bm{A}^{12,23}\\ \hline\cr\bm{A}^{12,1T}&\bm{A}^{22,1}&\bm{A}^{12,1,1}&\bm{A}^{12,1,2}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}\\ \hline\cr\bm{A}^{12,11T}&\bm{A}^{12,1,1T}&\bm{A}^{22,11}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}\\ \hline\cr\bm{A}^{12,12T}&\bm{A}^{12,1,2T}&{\LARGE\mbox{$\times$}}&\bm{A}^{22,12}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}\\ \hline\cr\bm{A}^{12,2T}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&\bm{A}^{22,2}&\bm{A}^{12,2,1}&\bm{A}^{12,2,2}&\bm{A}^{12,2,3}\\ \hline\cr\bm{A}^{12,21T}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&\bm{A}^{12,2,1T}&\bm{A}^{22,21}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}\\ \hline\cr\bm{A}^{12,22T}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&\bm{A}^{12,2,2T}&{\LARGE\mbox{$\times$}}&\bm{A}^{22,22}&{\LARGE\mbox{$\times$}}\\ \hline\cr\bm{A}^{12,23T}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&\bm{A}^{12,2,3T}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&\bm{A}^{22,23}\end{array}\right],$$

$$\bm{a}\equiv\left[\begin{array}[]{c}\bm{a}_{1}\\ \hline\cr\bm{a}_{2,1}\\ \hline\cr\bm{a}_{2,11}\\ \hline\cr\bm{a}_{2,12}\\ \hline\cr\bm{a}_{2,2}\\ \hline\cr\bm{a}_{2,21}\\ \hline\cr\bm{a}_{2,22}\\ \hline\cr\bm{a}_{2,23}\end{array}\right]\quad\mbox{and}\quad\bm{x}\equiv\left[\begin{array}[]{c}\bm{x}_{1}\\ \hline\cr\bm{x}_{2,1}\\ \hline\cr\bm{x}_{2,11}\\ \hline\cr\bm{x}_{2,12}\\ \hline\cr\bm{x}_{2,2}\\ \hline\cr\bm{x}_{2,21}\\ \hline\cr\bm{x}_{2,22}\\ \hline\cr\bm{x}_{2,23}\end{array}\right].$$

Full text

Solutions to Multilevel Sparse Matrix Problems

By Tui H. Nolan and Matt P. Wand

University of Technology Sydney

11th March, 2020

Abstract

We define and solve classes of sparse matrix problems that arise in multilevel modeling and data analysis. The classes are indexed by the number of nested units, with two-level problems corresponding to the common situation in which data on level 1 units are grouped within a two-level structure. We provide full solutions for two-level and three-level problems and their derivations provide blueprints for the challenging, albeit rarer in applications, higher level versions of the problem. Whilst our linear system solutions are a concise recasting of existing results, our matrix inverse sub-block results are novel and facilitate streamlined computation of standard errors in frequentist inference as well as allowing streamlined mean field variational Bayesian inference for models containing higher level random effects.

Keywords: Best linear unbiased prediction; Linear mixed models; Longitudinal data analysis; Panel data; Small area estimation; Variational inference.

1 Introduction

Higher level sparse matrices arise in statistical models for multilevel data, such as units grouped according to geographical sub-regions or repeated measures on medical study patients (e.g. Goldstein, 2010). Other areas of statistics and econometrics that use essentially the same types of models are longitudinal data analysis (e.g. Fitzmaurice et al., 2008), panel data analysis (e.g. Baltagi, 2013) and small area estimation (e.g. Rao & Molina, 2015). Linear mixed models (e.g. McCulloch, Searle & Neuhaus, 2008) are the main vehicle for modeling, fitting and inference. While they can be extended to generalised linear mixed models to cater for skewed, categorical and count response data, ordinary linear mixed models for Gaussian responses have the most relevance for the sparse matrix results presented here.

Both frequentist and Bayesian estimates of the fixed and random effects can be expressed succinctly in terms of ridge regression-type expressions involving design matrices (e.g. Henderson, 1975). However, typically the design matrices are sparse and naïve computation of the fixed and random effects estimates for large numbers of groups is inefficient and storage-greedy. For example, a random intercept linear mixed model for data with $1,000$ groups and $100$ observations per group involves a random effects design matrix containing 100 million entries of which 99.9% are zeroes. Streamlined computation of the best linear unbiased predictors of the fixed and random effects are well-documented with Section 2.2 of Pinheiro & Bates (2000) being a prime example. The matrix algebraic notion of QR decomposition plays a central role in numerically stable least squares-based fitting of linear models (e.g. Gentle, 2007) and also arises in the current context. It is important to note that such computations are performed after estimates of the covariance matrix parameters have been obtained via approaches such as minimum norm quadratic unbiased estimation or restricted maximum likelihood. Streamlined computation of covariance matrix estimates is tackled in, for example, Longford (1987). Given the covariance matrix estimates, implementation-ready matrix algebraic results for streamlined standard error calculations are not, to the best of our knowledge, present in the existing literature. These rely on efficient extraction of sub-blocks of the inverses of potentially very large sparse symmetric matrices. Presentation of these results, in the form of four theorems, is our main novel contribution. In the interests of conciseness and digestibility, we do not delve into the linear mixed model ramifications here – which are long-winded due to the various cases that require separate treatment. This article is purely concerned with generic matrix algebraic facts and, whilst motivated by statistical analysis, is totally free of statistical concepts in its main results and derivations. Ramifications for statistical inference are described in Nolan, Menictas & Wand (2020).

There is a related, but essentially non-overlapping, literature concerning inversion of so-called arrowhead matrices, which are invertible matrices that have all entries equal to zero except for those on the main diagonal and in one row and one column. Hołubowski et al. (2015) explain that such matrices “often appear in areas of applied science and engineering such as head-positioning systems of hard disk drives or kinematic chains of industrial robots” and then develop a fast method for inversion of a generalisation of arrowhead matrices known as block arrowhead matrices. Other recent contributions of this type are Sebari Nejafi et al. (2015) and Stanimirović et al. (2019). It is important to note that, whilst block arrowhead matrices coincide with two-level sparse matrices, the class of problems treated here and their motivating statistical applications are different from the central goal of the arrowhead matrix inversion literature. In multilevel sparse matrix problems the full matrix inverse is not of interest but, instead, inverse matrix sub-blocks matching the non-sparse sections of the original matrix. Arrowhead matrix inversion methods provide much more than is required and, therefore, are overly slow for large sparse two-level matrix problems. In addition, multilevel sparse matrices beyond the two-level case do not have block arrowhead forms. In summary, multilevel sparse matrix problems are, in essence, distinct from arrowhead matrix inversion problems.

Motivated by applications described elsewhere (e.g. Nolan, Menictas & Wand, 2020), our main focus in this article is the provision of full, implementable results for both two-level and three-level sparse matrix problems with statistically relevant matrix inverse sub-blocks. Both general situations and least squares form situations, with QR decomposition enhancement, are covered. We believe that it is best to treat each higher level case separately. If a future application would benefit from the solution to the four-level version of sparse matrix problems treated here then, whilst notationally and algebraically challenging, our two-level and three-level derivations point the way to a solution.

We cover two-level sparse matrix problems in Section 2 and three-level sparse matrix problems in Section 3. Some concluding remarks are made in Section 4.

2 Two-level Sparse Matrix Problems

We begin by defining a two-level sparse matrix problem:

Definition 1**.**

Let $\bm{A}$ by a symmetric and invertible matrix of the form:

[TABLE]

where the dimensions of the sub-blocks of $\bm{A}$ are as follows:

[TABLE]

The two-level sparse matrix problem is defined to be:

1. (I)

solve the linear system

[TABLE] 2. (II)

obtain the sub-blocks of $\bm{A}^{-1}$ corresponding to the positions of the sub-blocks $\bm{A}_{11}$, $\bm{A}_{12,i}$ and $\bm{A}_{22,i}$, $1\leq i\leq m$, in $\bm{A}$.

If $p$ and $q$ are small relative to $m$ then matrices defined by $\bm{A}$ are sparse since, as $m\to\infty$, the leading term of the fraction of non-zero entries of $\bm{A}$ is $\{1+2(p/q)\}m^{-1}$. In motivating statistical applications, $m$ corresponds to sample size, which often is large, whilst $p$ and $q$ correspond to dimensions of model parameter spaces, which are moderate in size. In the area of longitudinal data analysis typical values are $p=q=2$. The dimension variable $m$ matches the number of subjects in the longitudinal study, which may be in the hundreds or even thousands. If $m=1000$ then, for $p=q=2$, only about 0.3% of the entries of $\bm{A}$ are non-zero. Throughout this article we assume that $p$ and $q$ are small relative to $m$. Similar remarks apply to the three-level case treated in Section 3.

Our solution to the two-level sparse matrix problem benefits from the following notation for the sub-matrices of $\bm{A}^{-1}$, $\bm{a}$ and $\bm{x}$:

[TABLE]

where $\times$ generically denotes sub-blocks of $\bm{A}^{-1}$ which are in the same positions as the $\bm{O}$ blocks in $\bm{A}$. The dimensions of the sub-vectors of $\bm{a}$ and $\bm{x}$ are:

[TABLE]

Armed with the notation in (2), part (II) of Definition 1 can be expressed as:

[TABLE]

In applications involving multilevel data analysis the sub-blocks corresponding to the $\times$ symbols are usually not of interest since they correspond to between-group covariances. On the other hand, the sub-blocks of $\bm{A}^{-1}$ which are in the same position as the non-zero sub-blocks of $\bm{A}$ are required for obtaining standard errors of within-group fits. In the case of mean field variational inference, these sub-blocks are sufficient for both coordinate ascent and message passing optimal parameter computation with minimal product restrictions. Details of how these sub-blocks of $\bm{A}^{-1}$ are used in linear mixed model inference are given in Nolan, Menictas & Wand (2020).

Theorem 1 provides a streamlined solution to this problem such that the number of operations is linear in $m$. An analogous expression for $|\bm{A}|$, the determinant of $\bm{A}$, is also provided.

Theorem 1**.**

Consider the two-level sparse matrix problem given by Definition 1 and suppose that all $\bm{A}_{22,i}$, $1\leq i\leq m$, are invertible. The solution to part (II) of Definition 1 is:

[TABLE]

The determinant of $\bm{A}$ is

[TABLE]

The solution to part (I) of Definition 1 is:

[TABLE]

A proof of Theorem 1 is given in Appendix A.1.

2.1 Least Squares Form and QR-decomposition Enhancement

In statistical applications involving linear mixed models, it is common for $\bm{A}$ to admit a least squares form that lends itself to a QR decomposition-based solution. QR decompositions of rectangular matrices are a numerically preferred method for solving least squares problems. A QR-decomposition of a rectangular $n\times p$ ($n\geq p$) matrix $\bm{X}$ involves representing $\bm{X}$ as

[TABLE]

where $\bm{Q}$ is an $n\times n$ orthogonal matrix and $\bm{R}$ is a $p\times p$ upper-triangular matrix. Such a decomposition affords computational stability for least squares problems.

Suppose that $\bm{x}$ is chosen to minimise the least squares criterion

[TABLE]

for matrices

[TABLE]

with sub-matrices and sub-vectors having dimensions:

[TABLE]

such that $\bm{B}$ is full rank. Then it is easily verified that the $\bm{x}$ that minimises (3) is the solution to

[TABLE]

Moreover, $\bm{A}$ is a two-level sparse matrix of the form given by (1). The non-zero sub-blocks of $\bm{A}$ and the sub-vectors of $\bm{a}$ are

[TABLE]

and

[TABLE]

The form of the matrix $\bm{B}$ in (4) arises in statistical models containing both fixed effects and random effects with two-level structure (e.g. Goldstein, 2010). Full details on this motivational connection are given in Nolan et al. (2020).

Theorem 2 extends Theorem 1 by employing a QR decomposition approach for the purpose of numerical stability. Here, and later, we use the following notation for matrices $\bm{M}_{1},\ldots,\bm{M}_{d}$ each having the same number of columns:

[TABLE]

Theorem 2**.**

Suppose that $\bm{A}$ and $\bm{a}$ admit the forms defined by (4) and (5) where $\bm{B}$ is full rank. Then the two-level sparse matrix problem may be solved using the following QR decomposition-based approach:

For $i=1,\ldots,m$:

1. (a)

Decompose \accentset{\mbox{\mbox{{\bullet}}}}{\bm{B}}_{i}=\bm{Q}_{i}\begin{bmatrix}\bm{R}_{i}\\ \bm{0}\end{bmatrix} such that $\bm{Q}_{i}^{-1}=\bm{Q}_{i}^{T}$ and $\bm{R}_{i}$ is upper-triangular. 2. (b)

Then obtain

[TABLE] 2. 2.

Decompose $\operatorname*{stack}\limits_{1\leq i\leq m}(\bm{C}_{2i})=\bm{Q}\begin{bmatrix}\bm{R}\\ \bm{O}\end{bmatrix}$ such that $\bm{Q}^{-1}=\bm{Q}^{T}$ and $\bm{R}$ is upper-triangular and let

[TABLE] 3. 3.

The solutions are

[TABLE]

and, for $1\leq i\leq m$,

[TABLE] 4. 4.

The determinant of $\bm{A}$ is

[TABLE]

A proof of Theorem 2 is in Appendix A.2.

Remarks:

In Theorem 2, Step 1 involves determination of $m$ upper triangular matrices $\bm{R}_{i}$, $1\leq i\leq m$, via QR-decomposition which is a standard procedure within most computing environments. Each of the matrix inversions in Step 3 involve $\bm{R}_{i}^{-1}$, which can be achieved rapidly via back-solving. 2. 2.

Calculations such as $\bm{Q}_{i}^{T}\bm{b}_{i}$ do not require storage of $\bm{Q}_{i}$ and ordinary matrix multiplication. Standard matrix algebraic programming languages are such that information concerning $\bm{Q}_{i}$ is stored in a compact form from which matrices such as $\bm{Q}_{i}^{T}\bm{b}_{i}$ can be efficiently obtained. 3. 3.

Pinheiro & Bates (2000; Section 2.2) make use of this QR decomposition-based approach for fitting two-level linear mixed models. However, their descriptions are restricted to the $\bm{x}_{1}$ and $\bm{x}_{2,i}$ formulae, not those for the sub-blocks of $\bm{A}^{-1}$.

Table 1 summarises the results of a numerical study for which $100$ $\bm{A}$ matrices of the form (1) were randomly generated and the solution to the streamlined two-level sparse matrix problem using Theorem 2 was compared to the naïve solution, where the sparse structure is ignored, for increasingly large versions of the problem. Throughout the study $p=q=2$ and $m$ ranged over the set $\{100,200,400,800,1600\}$. The $n_{i}$ values were generated uniformly on the set $\{30,\ldots,60\}$. To aid maximal speed, both approaches were implemented in the low-level language Fortran 77, with LINPACK subroutines used for numerical linear algebra. The naïve approach involved matrix inversion via Gaussian elimination. The study was run on a MacBook Air laptop computer with a 2.2 gigahertz processor and 8 gigabytes of random access memory. The table lists the average and standard deviation times in seconds across the $100$ replications.

From Table 1 we see that the streamlined solutions are delivered in a small fraction of a second, even for very large versions of the problem. By contrast, naïve computation takes several minutes for the $m=1600$ case and is more than ten thousand times slower.

3 Three-level Sparse Matrix Problems

Three-level sparse matrix problems are such that the $\bm{A}$ matrix is symmetric and invertible but the two-level structure is repeated down the main diagonal. The notation for the general case becomes difficult to digest, so we start with a concrete example of such an $\bm{A}$:

[TABLE]

The general three-level sparse matrix problem is defined by:

Definition 2**.**

Let $\bm{A}$ by a symmetric and invertible matrix with partitioning

[TABLE]

where $\bm{A}_{11}$ is $p\times p$ and

[TABLE]

such that for each $1\leq i\leq m$, $\bm{A}_{12,i}$ is $p\times q_{1}$, and for each $1\leq j\leq n_{i}$, $\bm{A}_{12,ij}$ is $p\times q_{2}$. The lower right block is

[TABLE]

where, for each $1\leq i\leq m$, $\bm{A}_{22,i}$ is $q_{1}\times q_{1}$, and for each $1\leq j\leq n_{i}$, $\bm{A}_{12,\,i,\,j}$ is $q_{1}\times q_{2}$ and $\bm{A}_{22,ij}$ is $q_{2}\times q_{2}$. The three-level sparse matrix problem is defined to be:

1. (I)

solve the linear system

[TABLE] 2. (II)

obtain the sub-blocks of $\bm{A}^{-1}$ corresponding to the positions of the sub-blocks $\bm{A}_{11}$, $\bm{A}_{12,i}$, $\bm{A}_{22,i}$, $1\leq i\leq m$, and $\bm{A}_{12,ij},\ \bm{A}_{12,\,i,\,j},\ \bm{A}_{22,ij}$, $1\leq i\leq m$, $1\leq j\leq n_{i}$, in $\bm{A}$.

For the example three-level sparse matrix given by (6), we have $m=2$, $n_{1}=2$ and $n_{2}=3$. To enhance digestibility we will use these values for $m$, $n_{1}$ and $n_{2}$ throughout our discussion regarding the three-level sparse matrix problem. However, this can be easily generalised to any three-level sparse matrix, where $m$ and $\{n_{i}\}_{1\leq i\leq m}$ are arbitrary.

The solution to the thee-level sparse matrix problem benefits from the following notation for the sub-matrices of $\bm{A}^{-1}$, $\bm{a}$ and $\bm{x}$:

[TABLE]

[TABLE]

The dimensions of the partitioned vectors are:

• •

$\bm{a}_{1}$ and $\bm{x}_{1}$ are $p\times 1$ vectors;

• •

for each $1\leq i\leq m$, $\bm{a}_{2,i}$ and $\bm{x}_{2,i}$ are $q_{1}\times 1$ vectors;

• •

for each $1\leq i\leq m$ and $1\leq j\leq n_{i}$, $\bm{a}_{2,ij}$ and $\bm{x}_{2,ij}$ are $q_{2}\times 1$ vectors.

As in Section 2, $\times$ denotes the blocks of $\bm{A}^{-1}$ that are not of interest. Using the notation exemplified by (7), part (II) of Definition 2 can be expressed as:

[TABLE]

Theorem 3 presents the solution for any matrix that has the same sparsity structure as $\bm{A}$ when $m$ and $\{n_{i}\}_{1\leq i\leq m}$ are arbitrary. An analogous expression for $|\bm{A}|$ is also provided.

Theorem 3**.**

Consider the three-level sparse matrix problem given by Definition 2 and suppose that all $\bm{A}_{22,ij}$, $1\leq i\leq m$, $1\leq j\leq n_{i}$, are invertible. For $1\leq i\leq m$, define

[TABLE]

and suppose that all $\bm{H}_{22,i}$, $1\leq i\leq m$, are invertible. The solution to part (II) of Definition 2 is:

[TABLE]

The determinant of $\bm{A}$ is

[TABLE]

The solution to part (I) of Definition 2 is:

[TABLE]

A proof of Theorem 3 is in Appendix A.3.

3.1 Least Squares Form and QR-decomposition Enhancement

The three-level sparse matrix problem also lends itself to QR-decomposition enhancement. For the special case of $\bm{A}$, with $m=2$, $n_{1}=2$ and $n_{2}=3$, the least squares criterion has the form (3) with

[TABLE]

For general values of $m$ and $\{n_{i}\}_{1\leq i\leq m}$, the forms of $\bm{B}$ and $\bm{b}$ are

[TABLE]

and

[TABLE]

For each $1\leq i\leq m$ and $1\leq j\leq n_{i}$, the dimensions of the sub-blocks of $\bm{B}$ and $\bm{b}$ are:

[TABLE]

We also assume that $\bm{B}$ is full rank. Then the $\bm{x}$ that minimises (3) is the solution to the three-level sparse linear system with

[TABLE]

For general $m$ and $\{n_{i}\}_{1\leq i\leq m}$ the non-zero components of $\bm{A}$ and the sub-vectors of $\bm{a}$ are, for $1\leq i\leq m$,

[TABLE]

and

[TABLE]

The form of the matrix $\bm{B}$ in (9) arises in statistical models containing both fixed effects and random effects with three-level structure (e.g. Goldstein, 2010). Fuller details on this connection are given in Nolan et al. (2020).

Theorem 4 provides a QR decomposition enhancement of Theorem 3 for the least squares forms situation.

Theorem 4**.**

Suppose that $\bm{A}$ and $\bm{a}$ admit the least squares forms defined by (9)–(11) where $\bm{B}$ is full rank. Then the three-level sparse matrix problem may be solved using the following QR decomposition-based approach:

For $i=1,\ldots,m$:

1. (a)

For $j=1,\ldots,n_{i}$:

1. i.

Decompose \accentset{\mbox{\mbox{{\bullet}}\mbox{{\bullet}}}}{\bm{B}}_{ij}=\bm{Q}_{ij}\begin{bmatrix}\bm{R}_{ij}\\ \bm{0}\end{bmatrix} such that $\bm{Q}_{ij}^{-1}=\bm{Q}_{ij}^{T}$ and $\bm{R}_{ij}$ is upper-triangular. 2. ii.

Then obtain

[TABLE] 2. (b)

1. i.

Decompose \operatorname*{stack}\limits_{1\leq\,j\leq n_{i}}(\accentset{\mbox{\mbox{{\bullet}}}}{\bm{D}}_{2ij})=\bm{Q}_{i}\begin{bmatrix}\bm{R}_{i}\\ \bm{0}\end{bmatrix} such that $\bm{Q}_{i}^{-1}=\bm{Q}_{i}^{T}$ and $\bm{R}_{i}$ is upper-triangular. 2. ii.

Then obtain

[TABLE] 2. 2.

Decompose $\operatorname*{stack}\limits_{1\leq i\leq m}(\bm{C}_{2i})=\bm{Q}\begin{bmatrix}\bm{R}\\ \bm{0}\end{bmatrix}$ such that $\bm{Q}^{-1}=\bm{Q}^{T}$ and $\bm{R}$ is upper-triangular and let

[TABLE] 3. 3.

The solutions are, for $1\leq i\leq m$,

[TABLE]

and, for $1\leq i\leq m$, $1\leq j\leq n_{i}$,

[TABLE] 4. 4.

The determinant of $\bm{A}$ is

[TABLE]

Appendix A.4 contains a proof of Theorem 4.

Remarks:

As in Theorem 2, Step 1(a) of Theorem 4 involves determination of $\sum_{i=1}^{m}n_{i}$ upper-triangular matrices $\bm{R}_{ij}$, for $1\leq i\leq m$, $1\leq j\leq n_{i}$, via QR-decomposition. In Step 1(b), $m$ upper-triangular matrices $\bm{R}_{i}$, for $1\leq i\leq m$, are also constructed. A final QR-decomposition is applied in Step 2. Each of the inversions in Step 3 can be solved rapidly via back-solving. 2. 2.

The solutions for $\bm{x}_{1},\ \bm{A}^{11}$ and $\bm{x}_{2,i},\bm{A}^{12,i},\bm{A}^{22,i},1\leq i\leq m,$ have the same forms as in Theorem 4 for two-level sparse matrices. The solutions for $\bm{x}_{2,ij}$, $\bm{A}^{12,ij}$, $\bm{A}^{12,\,i,\,j}$ and $\bm{A}^{22,ij}$ are suggestive of a hierarchical pattern emerging for four-level and higher level classes of the problem.

4 Conclusion

In this short communication we have conveyed the essence of higher level sparse matrix problems as viewed through the prism of fitting and inference for multilevel statistical models. Both time-honoured best linear unbiased prediction and new-fashioned mean field variational Bayes approaches benefit from our four theorems for the two-level and three-level situations, with details given in Nolan, Menictas & Wand (2020). Future extensions to higher level situations are aided by our results and derivations.

Acknowledgements

We are grateful to Gioia Di Credico, Luca Maestrini, Marianne Menictas, Simon Wood and two anonymous referees for their assistance. This research was partially supported by the Australian Research Council Discovery Project DP140100441.

Appendix A Appendix: Proofs of Theorems

A.1 Proof of Theorem 1

In the case of $m=2$ the two-level sparse matrix linear system problem is

[TABLE]

which immediately leads to

[TABLE]

and

[TABLE]

It is clear that the same pattern applies for general $m$, and we have

[TABLE]

and

[TABLE]

Conditions (13) immediately imply that

[TABLE]

Substitution of (14) into (12) then leads to the solution for $\bm{x}_{1}$ as stated in Theorem 1.

For the matrix inverse derivation, we again start with the $m=2$ case and note that

[TABLE]

Observing the pattern from the $m=2$ case and then extending to general $m$ we obtain the system of equations:

[TABLE]

From (17)

[TABLE]

Substitution of (18) into (15) gives

[TABLE]

which implies that

[TABLE]

Substitution of (18) into (16) gives

[TABLE]

implying that

[TABLE]

For the $|\bm{A}|$ result we first prove:

Lemma 1. Let $\bm{M}$ be a symmetric invertible matrix with sub-block partitioning according to the notation

[TABLE]

Then

[TABLE]

Proof of Lemma 1. Lemma 1 is a direct consequence of Theorem 13.3.8 of Harville (2008) concerning the determinant of a matrix with $2\times 2$ sub-block partitioning.

From Lemma 1 we have

[TABLE]

A.2 Proof of Theorem 2

We first note the following simplification:

[TABLE]

where the first equality holds by the orthogonality of $\bm{Q}_{i}$ and the second and third equalities hold by Step 1(b) of Theorem 2. A similar sequence of steps can be used to show that

[TABLE]

These simplifications allow us to represent the non-zero components of $\bm{A}$ and the sub-vectors of $\bm{a}$ as

[TABLE]

and

[TABLE]

The derivation of the inverse matrix problem is as follows:

[TABLE]

where the fifth equality holds by Step 2 of Theorem 2, and we have used the orthogonality of $\bm{Q}$ for the sixth equality. The other components of $\bm{A}^{-1}$ are found by simply substituting the above simplifications into Theorem 1.

The derivation of the solution to the linear system is:

[TABLE]

where we have used Step 2 of Theorem 2 for the fifth and sixth equalities. The other sub-vectors of $\bm{x}$ are found by simply substituting the above simplifications into Theorem 1.

For the $|\bm{A}|$ result, the $|\bm{A}|$ expression from Theorem 1 implies that

[TABLE]

Then

[TABLE]

where we have used the result $|\bm{M}^{T}|=|\bm{M}|$ for any square matrix $\bm{M}$ and the fact that the determinant of an upper-triangular matrix is the product of its diagonal entries (Lemma 13.1.1 of Harville (2008)). Replacement $\bm{R}$ with $\bm{R}_{i}$ in (20) and substitution into (19) leads to the stated result for $|\bm{A}|$.

A.3 Proof of Theorem 3

In the case of $m=2$, $n_{1}=2$ and $n_{2}=3$ the three-level sparse matrix linear system problem is

[TABLE]

For arbitrary values of $m$ and $\{n_{i}\}_{1\leq i\leq m}$, we immediately obtain the following set of equations:

[TABLE]

Conditions (23) imply that

[TABLE]

which is the solution for $\bm{x}_{2,ij}$ as stated in Theorem 3. Substituting this result into conditions (22) leads to

[TABLE]

Solving this equation for $\bm{x}_{2,i}$ and using the definitions for $\bm{H}_{12,i}$ and $\bm{H}_{22,i}$ leads to the solution for $\bm{x}_{2,i}$ as stated in Theorem 3. Finally, substitution of the results for $\bm{x}_{2,i}$ and $\bm{x}_{2,ij}$ into condition (21) gives the solution for $\bm{x}_{1}$ stated in Theorem 3.

For the matrix inverse, we again illustrate the problem with $m=2$, $n_{1}=2$ and $n_{2}=3$:

[TABLE]

Observing the pattern for the $m=2$, $n_{1}=2$ and $n_{2}=3$ case and extending to general $m$ and $\{n_{i}\}_{1\leq i\leq m}$, we obtain the following system of equations:

[TABLE]

Rearranging conditions (29) we obtain

[TABLE]

from which the result for $\bm{A}^{12,ij}$ stated in Theorem 3 quickly follows. Rearranging conditions (28) we obtain

[TABLE]

Substituting the results for each $\bm{A}^{12,ij}$ leads to

[TABLE]

Solving the above set of equations for each $\bm{A}^{12,i}$ and using the definitions of $\bm{H}_{12,i}$ and $\bm{H}_{22,i}$ we obtain the result in Theorem 3. Substituting the results for each $\bm{A}^{12,ij}$ and $\bm{A}^{12,i}$ into (25), using the definition of $\bm{H}_{12,i}$ and solving for $\bm{A}^{11}$ leads to its stated result in Theorem 3. Rearranging conditions (30) we obtain

[TABLE]

from which the result for $\bm{A}^{12,\,i,\,j}$ stated in Theorem 3 follows quickly. Substitution of these results into the corresponding equations of condition (26) leads to

[TABLE]

Solving the above set of equations for each $\bm{A}^{22,i}$ and using the definitions of $\bm{H}_{12,i}$ and $\bm{H}_{22,i}$ in equation (7) of the main article we obtain the result in Theorem 3. Rearrangement of (27) leads to

[TABLE]

From Lemma 1 in the proof of Theorem 1,

[TABLE]

Application of Lemma 1 and Theorem 1 to the two-level sparse matrices in the last-written expression leads to

[TABLE]

and the stated result for $|\bm{A}|$ follows immediately.

A.4 Proof of Theorem 4

We first note the following re-expression:

[TABLE]

where the first equality holds by orthogonality of $\bm{Q}_{ij}$ and the second equality holds by Step 1(a)ii of Theorem 4. A similar sequence of steps can be used to show that

[TABLE]

We also have

[TABLE]

where the first equality holds by orthogonality of $\bm{Q}_{ij}$ and the second equality holds by Step 1(a)i of Theorem 4. A similar sequence of steps can be used to show that

[TABLE]

The above simplifications allow us to represent the non-zero components of $\bm{A}$ and the sub-vectors of $\bm{a}$ as

[TABLE]

and, for $1\leq i\leq m,\ 1\leq j\leq n_{i}$,

[TABLE]

Furthermore, each $\bm{H}_{22,i}$ matrix takes the form

[TABLE]

where the fifth equality holds by Step 1(b)ii of Theorem 4. Similarly,

[TABLE]

We are now in a position to derive the solutions to the inverse matrix problem:

[TABLE]

where the fourth equality holds by orthogonality of each $\bm{Q}_{i}$, the fifth equality holds by Step 1(b)ii, the eighth equality holds by Step 2 and the ninth equality holds by orthogonality of $\bm{Q}$. The other components of $\bm{A}^{-1}$ are found by simply substituting the above simplifications into Theorem 3 and using a similar sequence of steps.

For the linear system solution we have

[TABLE]

where the fourth and fifth equalities hold by Step 1(b)i and Step 1(b)ii, respectively. Also,

[TABLE]

The final task is to derive the expressions in Theorem 4 for the components of $\bm{x}$. For the first block we have:

[TABLE]

where the fourth equality holds by the orthogonality of $\bm{Q}_{i}$, the fifth equality holds by Step 1(b) and the ninth and eleventh equalities hold by Step 2 of Theorem 4. We get the other sub-vectors of $\bm{x}$ by a similar sequence of computations.

From Theorem 3 we have

[TABLE]

From part 3. of Theorem 4 we have $\bm{A}^{11}=\bm{R}^{-1}\bm{R}^{-T}$ so steps similar to those used for simplification of $|\bm{A}|$ in the proof of Theorem 2 lead to

[TABLE]

From the above arguments involving the $\bm{H}_{22,i}$ matrices,

[TABLE]

for $1\leq i\leq m$. Lastly,

[TABLE]

for $1\leq i\leq m,\ 1\leq j\leq n_{i}$. The stated result for $|\bm{A}|$ immediately follows.