Dyadic Factorization and Efficient Inversion of Sparse Positive Definite Matrices

TL;DR

The paper addresses efficient inversion of large sparse positive definite matrices by exploiting a hidden dyadic sparsity via nested dissection principles. It introduces a two-step framework: pack the matrix into a dyadic form and perform a sparse Gram–Schmidt orthogonalization, with a formal algebraic treatment of dyadic structures and a stable, level-by-level implementation. The key contributions are the dyadic factorization algorithm with complexity $O\left(dk^2\log^2\left({d}/{k}\right)\right)$ (and $O\left(dk^2\log\left({d}/{k}\right)\right)$ in block-tridiagonal cases) and a novel packing method based on the $l_1$ norm that uses distance matrices and multidimensional scaling to identify near-diagonal separators. The approach is implemented in the R package DyadiCarma, and simulations show robust performance for permuted band, block-tridiagonal, and dyadic matrices, including iterative strategies to recover dyadic structures from irregular patterns. Overall, the work provides a principled, scalable pathway to sparse matrix inversion by revealing and exploiting dyadic structure, with practical implications for high-dimensional statistics and numerical linear algebra.

Abstract

In inverting large sparse matrices, the key difficulty lies in effectively exploiting sparsity during the inversion process. One well-established strategy is the nested dissection, which seeks the so-called sparse Cholesky factorization. We argue that the matrices for which such factors can be found are characterized by a hidden dyadic sparsity structure. This paper builds on that idea by proposing an efficient approach for inverting such matrices. The method consists of two independent steps: the first packs the matrix into a dyadic form, while the second performs a sparse (dyadic) Gram-Schmidt orthogonalization of the packed matrix. The novel packing procedure works by recovering block-tridiagonal structures, focusing on aggregating terms near the diagonal using the $l_1$-norm, which contrasts with traditional methods that prioritize minimizing bandwidth, i.e. the $l_\infty$-norm. The algorithm performs particularly well for matrices that can be packed into banded or dyadic forms which are moderately dense. Due to the properties of $l_1$-norm, the packing step can be applied iteratively to reconstruct the hidden dyadic structure, which corresponds to the detection of separators in the nested dissection method. We explore the algebraic properties of dyadic-structured matrices and present an algebraic framework that allows for a unified mathematical treatment of both sparse factorization and efficient inversion of factors. For matrices with a dyadic structure, we introduce an optimal inversion algorithm and evaluate its computational complexity. The proposed inversion algorithm and core algebraic operations for dyadic matrices are implemented in the R package DyadiCarma, utilizing Rcpp and RcppArmadillo for high-performance computing. An independent R-based matrix packing module, supported by C++ code, is also provided.

Figures (6)

• Figure 1: Effect of increasing the neighborhood size for finding the optimal permutation for the matrix $\boldsymbol{\Sigma}$. The first row shows $\boldsymbol{\Sigma}(s)$, $s=1,2,3,4$ for a $100 \times 100$ full band matrix $\boldsymbol{\Sigma}$ with a half-bandwidth $\lambda = 20$. The second row shows the matrices recovered by applying the permutations obtained from the packing algorithm to the randomly permuted matrices in the first row.
• Figure 2: The average half-width $\Vert \boldsymbol l^{(\pi)} \Vert_1/d$ (the first two columns) and $\Vert \boldsymbol l^{(\pi)} \Vert_\infty$ (the remaining two columns) as a function of $p$. Rows correspond to different values of parameter $s \in \{1,2,3\}$ and columns correspond to different values of the parameter $\lambda \in \{10,40\}$. The blue lines in each plot represent values calculated over the initial (before permutation) matrices, while the red lines correspond to averages calculated over the final matrices obtained by the packing algorithm.
• Figure 3: The first plot presents the average $F(s)$ obtained for band matrices, $\lambda =10$. The second and the third plots present the average proportion of nonzero elements within the tridiagonal block structure returned by the packing algorithm applied to randomly permuted tridiagonal matrices with widths $k=10$ and $k=40$. The graphs are functions of the level of filling $p$ for $s \in \{1,2,3\}$.
• Figure 4: The average half-width $\Vert \boldsymbol l^{(\pi)} \Vert_1/d$ for permuted tridiagonal matrices shown as a function of $p$ with $N =4$. Rows correspond to $s \in \{1,2,3\}$, and columns to $k \in \{10,40\}$. The blue lines represent values calculated over the initial (not permuted) matrices, while the red lines correspond to averages calculated over the final matrices obtained by the packing algorithm.
• Figure 5: The performance of the packing algorithm for a banded dyadic matrix, $\lambda=60$, (the first two columns) and a dyadic matrix (the two remaining columns), ($N = 5, k = 10, p = 0.5$). The top graphs display the original (non-permuted) matrices along with their second power. The bottom graphs show the matrices recovered by the algorithm.

Theorems & Definitions (63)

• Definition 2.1
• Remark 2.2
• Definition 2.3
• Remark 2.4
• Example 1
• Theorem 2.5
• Remark 2.6
• Proposition 2.7
• Theorem 2.8
• Corollary 2.9