July 07, 2026
Let \(p\) be a prime and \({\mathsf{R}}={\mathbb{Z}}/p^2{\mathbb{Z}}\) the ring of integers modulo \(p^2\). Any \(A\in{\mathsf{R}}^{n\times n}\) is unimodularly equivalent to its Smith form \[S=\operatorname{diag}\bigl(\underbrace{1,\ldots,1}_{r_0}, \underbrace{p,\ldots,p}_{r_1}, \underbrace{0,\ldots,0}_{r_2}\bigr) \in{\mathsf{R}}^{n\times n},\] i.e., there exist \(U,V\in{\mathsf{R}}^{n\times n}\) such that \(UAV=S\), with \(\det U,\det V\in{\mathsf{R}}^*\) (where \({\mathsf{R}}^*\) is the set of units in \({\mathsf{R}}\), elements not equivalent to \(0\bmod p\)). Our goal in this paper is to determine \(r_0,r_1,r_2\) quickly when \(A\) is sparse or structured. By “sparse” we mean \(A\) is given by a black box such that for any \(v\in{\mathsf{R}}^{n\times 1}\) we can compute \(v\mapsto Av\) with \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n)\) operations in \({\mathsf{R}}\), which captures having few nonzero elements or a multiplicative structure (e.g., Hankel or Toeplitz matrices). We present a randomized algorithm which requires an expected number of \[{O\mskip 1mu\tilde{\,}\mskip 1mu}\left(n^{3-1/(\omega-1)}\right) operations in {\mathsf{R}}\]
to compute the Smith form, where \(\omega\) is the exponent of dense matrix multiplication. Using standard cubic matrix multiplication (\(\omega=3\)) our algorithm thus requires \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n^{2.5})\) operations in \({\mathsf{R}}\), while using the current asymptotically fastest matrix multiplication, with \(\omega<2.371339\), our algorithm requires \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n^{2.270786})\) operations in \({\mathsf{R}}\). Our algorithm is probabilistic of the Monte Carlo type, meaning it fails on any invocation with controllably small probability. We employ iterative block-Wiedemann-style matrix techniques and structured preconditioners. To our knowledge, this is the first algorithm to compute the modular Smith Normal Form modulo \(p^2\) requiring fewer than \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n^\omega)\) operations in \({\mathsf{R}}\), i.e., faster than any dense algorithm.
The Smith normal form (\(\operatorname{SNF}\)), introduced by [1], diagonalizes matrices over \({\mathbb{Z}}\) and has applications in diophantine analysis [2], integer programming [3], combinatorics [4], determining the structure of abelian groups [5] and class groups [6], in control system theory [7], and especially computing simplicial homology [8]. To deal with large instances of these problems we must exploit sparsity in the matrix (see, e.g., [9] in an application to neural data). This remains challenging.
The problem of computing the Smith form of a sparse matrix over a principal ideal ring presents several challenges. The mathematics goes back to [10]. One approach is to simply compute the SNF over \({\mathbb{Z}}\) and then reduce the result modulo the power of the prime ideal. The algorithm of [11] for SNF of a sparse matrix over \({\mathbb{Z}}\) could be used, which reduces the cost using standard matrix arithmetic over \({\mathbb{Z}}\) to \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n^{2.66})\) bit operations, but does not seem directly amenable to asymptotic speedup. The best known algorithm for any matrix over any \({\mathbb{Z}}/(d)\), for any \(d\in{\mathbb{Z}}\), is by [12] and requires \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n^\omega)\) operations mod \(d\), and is not sensitive to sparsity. The “valence” algorithm of [8] uses a hybrid of a “global” method over \({\mathbb{Z}}\) and a (dense) \(p\)-adic elimination method for selected primes \(p\). Our method here may be useful as a subroutine in this approach.
We will address a simple but important case here to attempt to break through this barrier: integer matrices modulo \(p^2\) for a prime \(p\). For matrices over \({\mathsf{R}}={\mathbb{Z}}/p^2{\mathbb{Z}}\), there are only three possible Smith factors, and the SNF has shape \[S=\operatorname{diag}\bigl(\underbrace{1,\ldots,1}_{r_0}, \underbrace{p,\ldots,p}_{r_1}, \underbrace{0,\ldots,0}_{r_2}\bigr) = \operatorname{diag}\bigl(1^{r_0},p^{r_1},0^{r_2}\bigr) \in{\mathsf{R}}^{n\times n}.\]
The first multiplicity is immediate: \[r_0=\operatorname{rank}(A\bmod p),\]
which can be identified by reducing the matrix to the prime field \({\mathbb{F}}_p\) and using well-established algorithms such as [13], which require \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n^2)\) operations in \({\mathsf{R}}\). The difficult part is to compute the number of invariant factors equal to \(p\) without forming dense nullspace bases, dense cokernel bases, or a dense Schur complement, all of which seem expensive to compute, even for sparse matrices.
The goal of this paper is, given a black box matrix over \({\mathsf{R}}\) as above, to find the SNF multiplicities (i.e., \(r_0,r_1,r_2\)) with fewer than \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n^\omega)\) operations in \({\mathsf{R}}\). We assume that two \(n\times n\) matrices can be multiplied with \(O(n^\omega)\) operations, where \(\omega=3\) using classical methods and \(\omega<2.371339\) using the fastest current methods [14].
The key new idea is to split off the unit part modulo \(p\). After fast rank-profile preconditioning, we write an equivalent matrix as \[A'=\begin{bmatrix} B&C\\ D&E \end{bmatrix}, \qquad B\bmod p~~\text{nonsingular},\qquad \dim B=r_0.\]
The Schur complement \[J=E-DB^{-1}C \bmod p^2\] is divisible by \(p\). Thus \(J=pT\) for a matrix \(T\) over the residue field, and \[r_1=\operatorname{rank}(T),\qquad r_2=n-r_0-r_1.\]
The matrix \(T\) is never formed. To evaluate \(TX\), lift \(X\) to \(\widetilde{X}\), solve \(BY=C\widetilde{X}\bmod p^2\), and return \[\frac{E\widetilde{X}-DY}{p}\bmod p,\]
or equivalently the lower block of \[A' \begin{bmatrix}-Y\\ \widetilde{X} \end{bmatrix}\]
divided by \(p\). We show how to perform the solve by Hensel lifting, with residuals computed in the length-two ring before division.
Another important tool we use is efficient sparse blocking. The expensive solves with \(B\bmod p\) are carried out by the sparse-block solve machinery of [11], [15], which provides efficient block projections. To avoid relying on a sparse-projection version of arbitrary-rank block Wiedemann, we compute rank by [16]’s ([16]) binary-search reduction to nonsingularity tests, using the PRECONDSXS/rank-profile conditioners of [17]. Each nonsingularity test is a verified sparse-block solve attempt.
In the nontrivial case \(0<r_0<n\), the sparse-block rank computation will use optimized block size \[t=\min\{r_0,n-r_0,\lceil n^{1/(\omega-1)}\rceil\}.\] The resulting Monte Carlo complexity is \[{O\mskip 1mu\tilde{\,}\mskip 1mu}\left(n^{3-1/(\omega-1)}\right) operations in {\mathsf{R}}.\]
Since \(3-1/(\omega-1)<\omega\) for every \(2<\omega\leq 3\), this is asymptotically below dense \(n^\omega\) arithmetic. The rank computations are Monte Carlo; the sparse-block linear solves inside the algorithm are verified by substitution.
The problem is closely related to earlier algorithms for Smith forms over local rings investigated in [18]. The contribution here is the \(p^2\) specialization in which the matrix \(T=J/p\bmod p\) can be accessed by a carry-safe block oracle, and sparse-block solves amortize the cost of the unit-block inversions enough to obtain a uniform sub-\(n^\omega\) bound for arbitrary sparse input matrices.
The remainder of the paper gives the local algebra, the finite-field sparse-block rank/solve procedures, the carry-safe oracle, the correctness proof, and the complexity analysis.
The base ring for our computations is \({\mathsf{R}}={\mathbb{Z}}/p^2{\mathbb{Z}}\). If the residue field \({\mathsf{R}}/p{\mathsf{R}}\) is too small for the randomized field preconditioners, we extend scalars to the unramified Galois ring \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}=\operatorname{GR}(p^2,d)\), whose residue field is \({\mathsf{F}}= \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}/p \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\cong {\mathbb{F}}_{p^d}\). When no extension is needed, take \(d=1\) and \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}={\mathsf{R}}\), \({\mathsf{F}}={\mathbb{F}}_p\). We will always count operations over \({\mathsf{R}}\) (or \({\mathsf{F}}\), as less expensive than \({\mathsf{R}}\)). Since any extension used has degree \(d=O(\log n + \log(1/\epsilon))\) we will assume these are captured by the polylogarithmic constants in the \({O\mskip 1mu\tilde{\,}\mskip 1mu}(\cdot)\).
The maximal ideal of \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\) is \(p \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\) and \(p^2=0\). Multiplication by \(p\) induces an isomorphism of \({\mathsf{F}}\)-vector spaces \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}/p \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\cong p \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\), \(\bar x\mapsto px\). Thus if \(z\in p \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\), the notation \(z/p\bmod p\) is unambiguous.
For our black box model we assume only that, for any vector \(v\in{\mathsf{R}}^{n\times 1}\), we can evaluate \(v\mapsto Av\) with \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n)\) operations in \({\mathsf{R}}\) (and similarly over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\)). For a block \(V\in{\mathsf{R}}^{n\times s}\), applying this product to each column gives \(V\mapsto AV\) in \({O\mskip 1mu\tilde{\,}\mskip 1mu}(ns)\) operations.
If the input black box is supplied only over the base ring \({\mathsf{R}}\) and the algorithm works over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\), fix an \({\mathsf{R}}\)-basis of the free unramified extension \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}/{\mathsf{R}}\). Since the entries of \(A\) lie in \({\mathsf{R}}\), a product with a vector in \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}^n\) decomposes into \(d\) products with vectors in \({\mathsf{R}}^n\), plus \({O\mskip 1mu\tilde{\,}\mskip 1mu}(nd)\) scalar operations for basis conversion and recombination. In the contexts below \(d=O(\log n+\log(1/\epsilon))\) for fixed or word-sized \(p\), and this overhead is absorbed into the \({O\mskip 1mu\tilde{\,}\mskip 1mu}(\cdot)\) bounds.
Over \({\mathsf{R}}\) or \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\), every Smith form modulo \(p^2\) has the shape \(\operatorname{diag}(1^{r_0},p^{r_1},0^{r_2})\). Our algorithm computes the multiplicities \(r_0,r_1,r_2\). Note that we do not compute the unimodular transformation matrices.
The following shows that moving to an enlarged ring does not change the Smith form.
Lemma 1. Let \(A\in {\mathsf{R}}^{n\times n}\) and extend scalars to \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}=\operatorname{GR}(p^2,d)\). If \[\operatorname{SNF}_{\mathsf{R}}(A)=\operatorname{diag}(1^{r_0},p^{r_1},0^{r_2}),\] then the Smith form of \(A\) over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\) has the same multiplicities.
Proof. Write a Smith decomposition over \({\mathsf{R}}\) as \[A=U\operatorname{diag}(1^{r_0},p^{r_1},0^{r_2})V, \qquad U,V\in\operatorname{GL}_n({\mathsf{R}}).\]
After scalar extension, the same equality holds over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\), and the images of \(U\) and \(V\) remain invertible because their determinants remain units. Thus the displayed diagonal is a Smith diagonal for \(A\) over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\).
The ring \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\) is again a length-two principal local ring with maximal ideal generated by the same element \(p\). Hence its Smith diagonal entries are uniquely determined up to multiplication by units and are represented by \(1\), \(p\), and \(0\). Therefore the multiplicities of these three types are unchanged. Notice that ranks in the extended algorithm are measured over \({\mathsf{F}}={\mathbb{F}}_{p^d}\); they are not multiplied by \(d\) when viewed as \({\mathsf{F}}\)-dimensions. 0◻ ◻
To avoid technicalities, the reader can generally assume that our ring \({\mathsf{R}}\) has a sufficiently large residue field that the algorithms work correctly, but explicit remediation for this is discussed below.
Let \(\bar A=A\bmod p\in {\mathsf{F}}^{n\times n}\). If \[\operatorname{SNF}(A)=\operatorname{diag}(1^{r_0},p^{r_1},0^{r_2}),\]
then \[\bar A\sim \operatorname{diag}(I_{r_0},0,0)\]
over \({\mathsf{F}}\). Therefore \[r_0=\operatorname{rank}_{\mathsf{F}}(\bar A).\] This part is immediate and is computed by a finite-field sparse rank algorithm with \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n^2)\) operations in \({\mathsf{R}}/(p)\).
Let \[a=\operatorname{rank}_{\mathsf{F}}(\bar A),\qquad k=n-a.\] Suppose, after invertible row and column operations over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\), we have \[A'= \begin{bmatrix} B&C\\ D&E \end{bmatrix}, \qquad B\in \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}^{a\times a},\]
with \(\bar B=B\bmod p\) nonsingular.
Lemma 2 (Schur complement reduction). The matrix \(B\) is invertible over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\), and \(A'\) is equivalent over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\) to \[\begin{bmatrix} B&0\\ 0&J \end{bmatrix}, \qquad J=E-DB^{-1}C\bmod p^2.\]
Moreover \(J\equiv 0\bmod p\).
Proof. Since \(\bar B\) is nonsingular over the residue field, \(\det(B)\) is a unit in the local ring \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\), so \(B\) is invertible. Multiplication by the unimodular block matrices \[\begin{bmatrix} I&0\\ -DB^{-1}&I \end{bmatrix} \quad\text{and}\quad \begin{bmatrix} I&-B^{-1}C\\ 0&I \end{bmatrix}\] gives \[\begin{bmatrix} I&0\\ -DB^{-1}&I \end{bmatrix} \begin{bmatrix} B&C\\ D&E \end{bmatrix} \begin{bmatrix} I&-B^{-1}C\\ 0&I \end{bmatrix} = \begin{bmatrix} B&0\\ 0&E-DB^{-1}C \end{bmatrix}.\] Modulo \(p\), the same block elimination is valid over \({\mathsf{F}}\). Since \(\operatorname{rank}(\bar A')=a\) and \(\bar B\) has rank \(a\), the lower-right Schur complement modulo \(p\) must have rank zero. Hence \(J\equiv 0\bmod p\). 0◻ ◻
Thus, there is a unique matrix \(T\in {\mathsf{F}}^{k\times k}\) such that \(J=pT\).
Lemma 3 (The Schur complement divided by \(p\) yields \(r_1\)). With the notation above, \[r_1=\operatorname{rank}_{\mathsf{F}}(T),\qquad r_2=k-\operatorname{rank}_{\mathsf{F}}(T).\]
Proof. The block \(B\) is invertible over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\), so it contributes exactly \(a\) unit invariant factors. It remains to compute the Smith form of \(J=pT\).
Let \(b=\operatorname{rank}_{\mathsf{F}}(T)\). Choose \(U_0,V_0\in\operatorname{GL}_k({\mathsf{F}})\) such that \[U_0 T V_0=\operatorname{diag}(I_b,0).\]
Lift \(U_0,V_0\) to matrices \(\widetilde{U},\widetilde{V}\in\operatorname{GL}_k( \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}})\), and choose any lift \(\widetilde{T}\in \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}^{k\times k}\) of \(T\). Since multiplication by \(p\) kills all terms divisible by \(p\), \[\widetilde{U}\,(p\widetilde{T})\,\widetilde{V} =p\,\widetilde{(U_0TV_0)} =\operatorname{diag}(pI_b,0)\pmod {p^2}.\]
Thus \(pT\) has \(b\) invariant factors equal to \(p\) and \(k-b\) invariant factors equal to \(0\). 0◻ ◻
This section adapts and summarizes the block matrix tools we need for our new algorithm. We make no claim of novelty here, but provide a summary of methods adapted to our particular context. Our finite field routines are based on the algorithms of [11], [15], which we will sometimes refer to as EGGSV. Rank is reduced to the binary-search nonsingularity tests of [16], using the preconditioners of [17], and each nonsingularity test is implemented by an EGGSV sparse-block solve attempt with a random right-hand side.
The only black-box operation in the Smith algorithm is a matrix-vector product with the original matrix \(A\), assumed to cost \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n)\) operations in \({\mathsf{R}}\). Whenever we need to apply \(A\) to an \(n\times s\) block, we apply the black box to the columns, at cost \({O\mskip 1mu\tilde{\,}\mskip 1mu}(ns)\). The finite-field statements below are written for an auxiliary matrix \(M\in {\mathsf{F}}^{m\times m}\) for which one vector product \(v\mapsto Mv\) costs \({O\mskip 1mu\tilde{\,}\mskip 1mu}(L)\) operations. Thus an \(m\times s\) block product by \(M\) costs \({O\mskip 1mu\tilde{\,}\mskip 1mu}(Ls)\) by columnwise application. In the Smith algorithm this will be used with \(L=n\) for leading submatrices of the preconditioned sparse matrix. The divided Schur complement \(T\) is not treated as a new input black box. Section 6 constructs its block oracle from the original black box for \(A\) and from the verified sparse-block solves with \(B_0\).
Throughout this section, \(1\le s\le m\) is a block size. If \(s\nmid m\), pad by fewer than \(s\) dimensions. For nonsingular solves use a direct-sum identity block, and for rank tests use a direct-sum zero block, or an identity block and subtract its known rank. This changes all costs by at most constant factors. Thus, we write \(m=qs\). The sparse-block projection used by EGGSV is \[u=\begin{bmatrix} I_s\\ \vdots\\ I_s \end{bmatrix}\in {\mathsf{F}}^{m\times s},\] so multiplication by \(u\) or \(u^\intercal\) on an \(m\times s\) block costs \(O(ms)\) operations, not \(O(ms^2)\) operations in \({\mathsf{F}}\).
The sparse-block projection and inverse factorization method of [11], [15] is at the foundation of our algorithm and reduced complexity. We summarize its adapted use, features and complexity as follows.
Theorem 1 (Efficient block projections and inverse factorization). Let \(N\in {\mathsf{F}}^{m\times m}\) be nonsingular and \(m=qs\). After fast rank-profile conditioning over a sufficiently large field, the sparse projection \(u\) above has the following properties with controllably high probability.
The right block Krylov matrix \[K^{(r)}=K_q(\widehat N,u)=[u,\widehat N u,\ldots,\widehat N^{q-1}u]\] and the block-row Krylov matrix \[K^{(\ell)}= \begin{bmatrix}u^\intercal\\ u^\intercal\widehat N\\ \vdots\\ u^\intercal\widehat N^{q-1} \end{bmatrix}\] are nonsingular. Here \(\widehat N\) is obtained from \(N\) by fast invertible conditioners and a block diagonal random scaling.
The block-Hankel matrix \[H=K^{(\ell)}\widehat N K^{(r)}\] is nonsingular and has \(s\times s\) Hankel blocks.
There is a factorization \[\widehat N^{-1}=K^{(r)}H^{-1}K^{(\ell)}.\] Undoing the fast conditioners gives the corresponding factorization of \(N^{-1}\).
Moreover, the conditioners apply to blocks in near-linear time, and the failure probability can be made at most \(\epsilon\) by working over an extension field of degree \(O(\log m+\log(1/\epsilon))\) and/or by repetition.
Proof. This is the sparse-block projection and inverse-factorization approach of [11], [15], stated in our context and notation. The sparse projection \(u\) is the projection in [11]. The Krylov nonsingularity statements under leading-minor hypotheses are [11] Theorem 2.1 and Corollary 2.2; the passage to arbitrary nonsingular matrices over sufficiently large fields by fast rank-profile conditioning is [11]. The block-Hankel inverse factorization is [11], and the Las Vegas inversion and inverse-times-matrix consequences are [11]. 0◻ ◻
The following is the parameterized form of the EGGSV inverse-times-matrix algorithm, specialized to one \(m\times s\) right-hand-side block and to the sparse projection \(u=[I_s,\ldots,I_s]^\intercal\).
Theorem 2. Let \(M\in {\mathsf{F}}^{m\times m}\) be nonsingular, and suppose one vector product \(v\mapsto Mv\) costs \({O\mskip 1mu\tilde{\,}\mskip 1mu}(L)\) operations. For \(C\in {\mathsf{F}}^{m\times s}\), there is a Las Vegas algorithm that computes \(M^{-1}C\) with expected cost \[{O\mskip 1mu\tilde{\,}\mskip 1mu}\left(mL+m^2+m s^{\omega-1}\right).\]
Proof. This is [11], parameterized by the block size \(s\) and applied to one \(m\times s\) right-hand-side block. Let \(m=qs\). The EGGSV construction uses the sparse projection \(u=[I_s,\ldots,I_s]^\intercal\), forms the block-Hankel Krylov data \(u^\intercal M^iu\), applies \(K^{(\ell)}\) to the right-hand-side block \(C\), solves the resulting structured block-Hankel system, and evaluates \(K^{(r)}\) on the result.
The cost is assessed as follows. The \(O(q)\) Krylov block products needed for the Hankel data, for \(K^{(\ell)}C\), and for the final \(K^{(r)}\) evaluation each cost \({O\mskip 1mu\tilde{\,}\mskip 1mu}(Ls)\) by \(s\) columnwise vector products. These contributions are \({O\mskip 1mu\tilde{\,}\mskip 1mu}(mL)\). The sparse projections by \(u^\intercal\) over the generated sequences contribute \({O\mskip 1mu\tilde{\,}\mskip 1mu}(m^2)\) field operations. The structured \(q\times q\) block-Hankel solve has \(s\times s\) blocks and costs \({O\mskip 1mu\tilde{\,}\mskip 1mu}(qs^\omega)={O\mskip 1mu\tilde{\,}\mskip 1mu}(ms^{\omega-1})\). Fast conditioning and unconditioning are absorbed in these terms.
The returned candidate is verified by checking \(MX=C\). Since \(M\) is nonsingular, bad internal choices are detected either by a singular Krylov/Hankel construction or by this final verification, and independent resampling gives the stated expected cost. 0◻ ◻
The full inverse algorithm of [11] is stated with access to both \(M\) and \(M^\intercal\). In the present paper we only use the inverse-times-a-block variant. For this task, the product \(K^{(\ell)}C\) is obtained from the Krylov sequence \[C,\; MC,\; M^2C,\;\ldots\]
by applying the sparse projection \(u^\intercal\) to each block. Thus the solve routines used here require only products by \(M\), not by \(M^\intercal\).
We can also adapt the EGGSV framework for testing nonsingularity.
Lemma 4. Let \(N\in {\mathsf{F}}^{r\times r}\), let \(1\leq h\le r\), and suppose a vector product by \(N\) costs \({O\mskip 1mu\tilde{\,}\mskip 1mu}(L_N)\) operations. There is a verified sparse-block solve test which requires \({O\mskip 1mu\tilde{\,}\mskip 1mu}\left(rL_N+r^2+r h^{\omega-1}\right)\) operations in \({\mathsf{F}}\) per trial and satisfies the following:
If \(N\) is nonsingular, a trial reports “nonsingular” with probability at least \(1-\eta\), where \(\eta\) is the chosen EGGSV conditioning failure probability. Repeating independent trials and accepting if any trial succeeds makes this false-negative probability arbitrarily small.
If \(N\) is singular with nullity \(\delta\ge1\), a trial reports “nonsingular” with probability at most \(|{\mathsf{F}}|^{-\delta h}\le |{\mathsf{F}}|^{-h}\). If \(R\) independent trials are used and the algorithm accepts when any trial succeeds, the singular false-positive probability is at most \(R |{\mathsf{F}}|^{-h}\). The extension field and the number of repetitions must be chosen so that \(R|{\mathsf{F}}|^{-h}\) is sufficiently small.
Proof. Choose a random block \(C\in{\mathsf{F}}^{r\times h}\) and run one verified sparse-block solve attempt for \(NX=C\) using Theorem 2 with block size \(h\). Report “nonsingular” exactly when the solve attempt returns a matrix \(X\) and the verification \(NX=C\) succeeds.
We use the same bounded EGGSV solve attempt as an algorithmic procedure on arbitrary input. On nonsingular \(N\), Theorem 2 provides its Las Vegas success behaviour. On singular \(N\), regardless of its internal behaviour, the verification \(NX=C\) can succeed only if every column of \(C\) lies in \(\operatorname{im}(N)\), and the false-positive probability is at most \(|{\mathsf{F}}|^{-\delta h}\).
The cost is the cost of Theorem 2 with dimension \(r\), block size \(h\), and vector-product cost \(L_N\), namely \({O\mskip 1mu\tilde{\,}\mskip 1mu}\left(rL_N+r^2+r h^{\omega-1}\right)\). 0◻ ◻
It is straightforward to employ a Wiedemann-style rank routine, using binary search, singularity testing, and fast solving to determine rank. Using [13]’s ([13]) method may be faster, but some infrastructure around using the block minimal polynomial would have to be developed.
Theorem 3. Let \(M\in {\mathsf{F}}^{m\times m}\) be arbitrary, and suppose one vector product \(v\mapsto Mv\) costs \({O\mskip 1mu\tilde{\,}\mskip 1mu}(L)\) operations. There is a Monte Carlo algorithm that returns \(\operatorname{rank}(M)\) with failure probability at most \(\epsilon\) and cost \[{O\mskip 1mu\tilde{\,}\mskip 1mu}\left(mL+m^2+m s^{\omega-1}\right).\]
The algorithm uses Wiedemann’s rank-by-singularity-reduction approach, the PRECONDSXS rank-profile preconditioner of [17], and the nonsingularity test of Lemma 4 from EGGSV. If the original field is too small, the computation is performed over an extension field of degree \(O(\log m+\log(1/\epsilon))\).
Proof. We recall [16]’s ([16]) rank-by-search reduction. For a candidate integer \(r\), the predicate is \(P(r): r\le \operatorname{rank}(M)\). The PRECONDSXS conditioning of [17], together with their reductions and fast conditioners, supplies fast invertible \(P_r,Q_r\) such that the leading \(r\times r\) block \(N_r=(P_rMQ_r)_{1:r,1:r}\) is nonsingular with high probability when \(r\le\operatorname{rank}(M)\). If \(r>\operatorname{rank}(M)\), every \(r\times r\) submatrix is singular.
Test \(N_r\) with Lemma 4, using \(h=\min(r,s)\). A vector product by \(N_r\) is obtained by embedding into \(m\) coordinates, applying the fast conditioners and \(M\), and restricting; this costs \({O\mskip 1mu\tilde{\,}\mskip 1mu}(L)\), since a product by \(M\) already has output size \(m\). Thus each predicate test costs at most \({O\mskip 1mu\tilde{\,}\mskip 1mu}\left(mL+m^2+m s^{\omega-1}\right)\), including the case \(r<s\). Binary search uses \(O(\log m)\) predicate calls. Repetition reduces false negatives when the predicate is true; when it is false, the only accepting events are the random-right-hand-side false positives of Lemma 4, bounded by \(R|{\mathsf{F}}|^{-h}\) for \(R\) trials. Choosing the extension field and repetition count so that each predicate has error at most \(\epsilon/O(\log m)\) and applying a union bound gives total failure probability at most \(\epsilon\). 0◻ ◻
We will employ the adapted methods of [11], [15], [16], [17], which we summarize as follows.
Corollary 1. Let \(M\in {\mathsf{F}}^{m\times m}\) be a finite-field matrix for which one vector product costs \({O\mskip 1mu\tilde{\,}\mskip 1mu}(L)\) operations. The following operations are available with the displayed costs:
If \(M\) is nonsingular, compute \(M^{-1}C\) for \(C\in {\mathsf{F}}^{m\times s}\) with expected \[{O\mskip 1mu\tilde{\,}\mskip 1mu}\left(mL+m^2+m s^{\omega-1}\right) operations in {\mathsf{F}} (Las Vegas).\]
Compute \(\operatorname{rank}(M)\) in an expected \[{O\mskip 1mu\tilde{\,}\mskip 1mu}\left(mL+m^2+m s^{\omega-1}\right) operations in {\mathsf{F}} (Monte Carlo).\]
by binary search on nonsingularity tests, with failure probability bounded by the chosen parameter.
Later we will apply Corollary 1 to specialize to leading submatrices of the preconditioned sparse matrices. In particular, suppose \(B_0=B\bmod p\in {\mathsf{F}}^{a\times a}\), a leading block of the preconditioned matrix \(A'\bmod p\). A vector product by \(B_0\) is implemented by embedding the vector into the first \(a\) coordinates of an \(n\)-vector, applying \(A'\bmod p\), and restricting to the first \(a\) coordinates. Hence a vector product by \(B_0\) costs \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n)\), and applying \(B_0\) to an \(a\times s\) block costs \({O\mskip 1mu\tilde{\,}\mskip 1mu}(ns)\) by applying this vector operation to the columns. If \(B_0\) is nonsingular, computing \(B_0^{-1}C\) for \(C\in{\mathsf{F}}^{a\times s}\) via Theorem 2 requires \({O\mskip 1mu\tilde{\,}\mskip 1mu}\left(an+a^2+a s^{\omega-1}\right)\) operations in \({\mathsf{F}}\).
The same rank-by-binary-search proof is later applied to the divided Schur complement operator \(T\). We do not introduce a separate black box for \(T\): Section 6 constructs an explicit \(s\)-column oracle for \(T\) from the original \(A\)-oracle and from the verified sparse-block solves with \(B_0\). The rank proof uses \(O(k/s)\) calls to this \(s\)-column black box, together with \({O\mskip 1mu\tilde{\,}\mskip 1mu}(k^2+k s^{\omega-1})\) additional field operations for sparse projections and block-Hankel algebra. This is the counting used in the cost analysis.
We note that the sparse-block projection, PRECONDSXS, and random right-hand-side tests are probabilistic and require a sufficiently large field for precise failure bounds. If \(p\) is too small, choose \({\mathsf{F}}={\mathbb{F}}_{p^d}\), for \(d=O(\log n+\log(1/\epsilon))\), so that \(|{\mathsf{F}}|\) is larger than the polynomial-size sets required by the random conditioners and so that the singular false-positive bound \(R|{\mathsf{F}}|^{-h}\) is negligible for the planned number \(R\) of repetitions. The algorithm is then run over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}=\operatorname{GR}(p^2,d)\). By Lemma 1, this does not change the Smith multiplicities. The extension degree contributes only polylogarithmic factors in the usual algebraic-operation model when \(p\) is fixed or word-sized.
Preconditioners have two distinct roles in our algorithms, which we outline below.
We need fast \(P,Q\in\operatorname{GL}_n( \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}})\) such that \(A'=PAQ\) has a leading \(a\times a\) block \(B\) with \(B\bmod p\) nonsingular, where \(a=\operatorname{rank}(A\bmod p)\). The required properties are:
\(P,Q\) are invertible over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\). This ensures that Smith multiplicities are preserved.
Their reductions \(P_0,Q_0\in\operatorname{GL}_n({\mathsf{F}})\) give generic rank profile to \(P_0(A\bmod p)Q_0\), or at least make its leading \(a\times a\) submatrix nonsingular.
\(P,Q\) apply to \(n\times s\) blocks in \({O\mskip 1mu\tilde{\,}\mskip 1mu}(ns)\) operations.
The following lemma summarizes our global preconditioning requirement.
Lemma 5 (Fast generic rank-profile preconditioning). Let \(M\in {\mathsf{F}}^{n\times n}\) have rank \(a\). Over a sufficiently large finite field, there are randomized fast invertible conditioners \(P_0,Q_0\in\operatorname{GL}_n({\mathsf{F}})\) such that \(P_0MQ_0\) has a nonsingular leading \(a\times a\) block with probability at least \(1-\epsilon\). The conditioners apply to \(n\times s\) blocks in \({O\mskip 1mu\tilde{\,}\mskip 1mu}(ns)\) operations. If the field is too small, the same statement holds after extension degree \(O(\log n+\log(1/\epsilon))\).
Proof. This is the standard generic-rank-profile/PRECONDSXS conditioning theorem used in black-box rank algorithms. [17] explicitly formulate PRECONDSXS as the condition that, for a requested size \(s\le\operatorname{rank}(M)\), the leading \(s\times s\) principal minor of the preconditioned matrix is nonzero; they also state its use in rank computation by binary search [17]. Their Theorem 3.1 reduces PRECONDSXS and generic rank profile to linear-independence preconditioners; Theorem 6.2 gives the butterfly routing network, and Theorem 6.3 gives the associated randomized linear-independence probability bound, yielding fast linear-independence conditioners with near-linear application cost [17] Theorems 3.1, 6.2, and 6.3. [13]’s ([13]) triangular Toeplitz conditioners give another standard fast generic-rank-profile family. EGGSV use the same kind of fast rank-profile conditioning, together with diagonal scalings and butterfly-type conditioners, to make the leading block minors and sparse-block Krylov determinants required by their efficient projection theorem nonzero over sufficiently large fields [11]. The bad choices are zeros of provably nonzero polynomials in the random parameters, so a Schwartz–Zippel bound, or the extension-field choice used in EGGSV, gives failure probability at most \(\epsilon\). 0◻ ◻
In the Smith algorithm, choose such sparse invertible preconditioners \(P_0,Q_0\) for \(M=A\bmod p\) and lift them coefficientwise to \(P,Q\in\operatorname{GL}_n( \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}})\). A lift is invertible over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\) because its determinant is nonzero modulo \(p\), hence a unit in the local ring. Therefore \(A'=PAQ\) is Smith-equivalent to \(A\), and the leading block \(B\bmod p\) is nonsingular whenever the preconditioning event succeeds.
The event that the leading \(a\times a\) block is nonsingular is checked probabilistically by computing \(\operatorname{rank}_{\mathsf{F}}(B\bmod p)\) using Corollary 1. If the reported rank is not \(a\), resample \(P,Q\). If a Monte Carlo rank check falsely reports full rank, this is included in the global Monte Carlo failure probability. The later solve verifications catch many such mistakes, but the proof of correctness does not rely on them as deterministic proofs of nonsingularity.
Inside the finite-field solve algorithms and the solve-based rank tests we use additional fast conditioners and sparse-block projections. These conditioners are not part of the Smith equivalence, but are internal to the finite-field routines. These conditioners must satisfy:
the relevant block Krylov matrices have full rank or the required generic degree profile with high probability;
applying the conditioners and sparse projections to \(m\times t\) blocks costs \({O\mskip 1mu\tilde{\,}\mskip 1mu}(mt)\) per Krylov step;
the total projection overhead over a sequence of length \({O\mskip 1mu\tilde{\,}\mskip 1mu}(m/t)\) is \({O\mskip 1mu\tilde{\,}\mskip 1mu}(m^2)\).
The last item is the key difference from using dense block projections inside the solve/rank tests.
This section presents the main new algorithmic ideas. All operations are done over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\) until an element is known to be divisible by \(p\).
Let \[A'= \begin{bmatrix} B&C\\ D&E \end{bmatrix}, \qquad B\in \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}^{a\times a}, \qquad B_0=B\bmod p\in {\mathsf{F}}^{a\times a}.\] Assume \(B_0\) is nonsingular. Let \[J=E-DB^{-1}C=pT, \qquad T\in {\mathsf{F}}^{k\times k}, \qquad k=n-a.\]
Given \(U\in \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}^{a\times t}\), solve \[BY=U\bmod p^2\]
as follows.
Reduce modulo \(p\): \[U_0=U\bmod p, \qquad B_0=B\bmod p.\]
Solve over \({\mathsf{F}}\): \[B_0Y_0=U_0.\]
Lift \(Y_0\) coefficientwise to \(\widetilde{Y}_0\in \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}^{a\times t}\).
Compute the residual in \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\): \[W=U-B\widetilde{Y}_0\bmod p^2.\]
Because \(B_0Y_0=U_0\), the matrix \(W\) lies in \(p \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}^{a\times t}\).
Divide by \(p\) only now: \[W_1=W/p\bmod p\in {\mathsf{F}}^{a\times t}.\]
Solve over \({\mathsf{F}}\): \[B_0Y_1=W_1.\]
Lift \(Y_1\) to \(\widetilde{Y}_1\in \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}^{a\times t}\) and return \[Y=\widetilde{Y}_0+p\widetilde{Y}_1.\]
Lemma 6 (Correctness of the Hensel solve). The matrix \(Y\) returned above satisfies \(BY=U\bmod p^2\).
Proof. By construction, \[U-B\widetilde{Y}_0=pW_1\bmod p^2.\] Also \(B_0Y_1=W_1\) over \({\mathsf{F}}\), so \[pB\widetilde{Y}_1=pB_0Y_1=pW_1\bmod p^2.\] Therefore \[B(\widetilde{Y}_0+p\widetilde{Y}_1) =B\widetilde{Y}_0+pB\widetilde{Y}_1 =B\widetilde{Y}_0+pW_1 =U\bmod p^2. \tag*{\qed}\] ◻
Suppose we write digit expansions \[B=B_0+pB_1, \qquad U=U_0+pU_1.\]
It is tempting, but incorrect, to set \[W_1=U_1-B_1Y_0.\]
The missing term is the carry from the product \(B_0Y_0\) in representatives. If \[B_0Y_0=U_0+pK\]
in canonical representatives, then the correct residual is \[W_1=U_1-B_1Y_0-K\bmod p.\]
Computing \[W_1=(U-B\widetilde{Y}_0)/p\bmod p\]
after the multiplication in \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\) automatically includes the carry \(K\).
Let \(X\in {\mathsf{F}}^{k\times t}\). The block oracle for \(TX\) is:
Lift \(X\) to \(\widetilde{X}\in \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}^{k\times t}\).
Compute over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\): \[A'\begin{bmatrix}0\\\widetilde{X}\end{bmatrix} = \begin{bmatrix}U\\V\end{bmatrix} = \begin{bmatrix}C\widetilde{X}\\E\widetilde{X}\end{bmatrix} \bmod p^2.\]
Use the Hensel solve above to compute \(Y\) satisfying \[BY=U\bmod p^2.\] This step involves two solves with \(B_0\) over \({\mathsf{F}}\) and one carry-safe residual computation in \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\).
Compute over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\): \[Z=A' \begin{bmatrix} -Y\\\widetilde{X} \end{bmatrix} = \begin{bmatrix} 0\\J\widetilde{X} \end{bmatrix} \bmod p^2.\]
The lower block of \(Z\) is divisible by \(p\). Return \[TX=(\text{lower block of }Z)/p\bmod p.\]
Lemma 7 (Correctness of the \(T\)-oracle). The oracle returns \(TX\), where \(T=J/p\bmod p\).
Proof. Step T2 gives \(U=C\widetilde{X}\). Step T3 computes \(Y=B^{-1}U\bmod p^2\). Therefore \[A'\begin{bmatrix} -Y\\ \widetilde{X} \end{bmatrix} = \begin{bmatrix} -BY+C\widetilde{X}\\ -DY+E\widetilde{X} \end{bmatrix} = \begin{bmatrix} 0\\ (E-DB^{-1}C)\widetilde{X} \end{bmatrix} = \begin{bmatrix} 0\\ J\widetilde{X} \end{bmatrix} \bmod p^2.\] Since \(J=pT\), the lower block is divisible by \(p\), and division by \(p\) gives \(TX\) over \({\mathsf{F}}\). The result is independent of the chosen lift of \(X\): replacing \(\widetilde{X}\) by \(\widetilde{X}+pW\) changes \(J\widetilde{X}\) by \(JpW=p^2TW=0\) modulo \(p^2\). 0◻ ◻
We now present our algorithm for sparse Smith Normal Form modulo \(p^2\).
Input: \(A\in {\mathsf{R}}^{n\times n}\) with fast sparse matrix-vector products.
Output: \((r_0,r_1,r_2)\) such that \(\operatorname{SNF}(A)=\operatorname{diag}(1^{r_0},p^{r_1},0^{r_2})\).
Choose an extension \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}=\operatorname{GR}(p^2,d)\) if needed so that the residue field \({\mathsf{F}}\) is large enough for the sparse-block randomization.
Compute \(a=\operatorname{rank}_{\mathsf{F}}(A\bmod p)\) and set \(r_0=a\) and \(k=n-a\).
If \(k=0\), return \((n,0,0)\).
If \(a=0\), then \(A\equiv 0\bmod p\), so \(A=pT\) for \(T\in {\mathsf{F}}^{n\times n}\). Compute \[r_1=\operatorname{rank}_{\mathsf{F}}(A/p\bmod p),\]
and return \((0,r_1,n-r_1)\).
Choose fast unimodular \(P,Q\in\operatorname{GL}_n( \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}})\) such that, for \[A'=PAQ=\begin{bmatrix}B&C\\D&E\end{bmatrix}, \qquad B\in \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}^{a\times a},\] the matrix \(B_0=B\bmod p\) is nonsingular. Verify this condition probabilistically by computing \(\operatorname{rank}_{\mathsf{F}}(B_0)\). If the check reports failure, resample.
Choose the working block size \[t=\min\left\{a,k,\left\lceil n^{1/(\omega-1)}\right\rceil\right\}.\]
This is the number of simultaneous right hand sides used in all subsequent black-box queries to the Schur-complement oracle. Thus the outer rank computation will query the black box on matrices \(X\in{\mathsf{F}}^{k\times t}\), and each such query reduces internally to solves with coefficient matrix \(B_0\) and \(t\) right hand sides. The bounds \(t\le k\) and \(t\le a\) ensure that the same block size is valid both for the \(k\times k\) matrix \(T\) and for the internal \(a\times a\) solves with \(B_0\).
The cap \(t\le \left\lceil n^{1/(\omega-1)}\right\rceil\) is the value used in the complexity analysis below; it balances the two main block-size dependent terms, namely the repeated sparse/block-access term \(akn/t\) and the dense block-arithmetic term \(ak t^{\omega-2}\).
Use the carry-safe oracle of Section 6 to supply, for every query matrix \(X\in{\mathsf{F}}^{k\times t}\), the block product \[X\longmapsto TX,\qquad T=(E-DB^{-1}C)/p\bmod p.\]
The matrix \(T\) is not formed explicitly.
Compute \(b=\operatorname{rank}_{\mathsf{F}}(T)\) using the rank-by-binary-search part of Corollary 1, with the above block size \(t\) and with access to \(T\) provided by the oracle in Step A7.
Return \((r_0,r_1,r_2)=(a,b,k-b)\).
Correctness is established in the following theorem.
Theorem 4 (Correctness). Assume the rank computations and sparse-block solves return correct results and the preconditioner in Step A5 succeeds. Then the algorithm returns the Smith multiplicities of \(A\) modulo \(p^2\).
Proof. Step A2 computes \(a=\operatorname{rank}(A\bmod p)\), which equals the number of unit invariant factors. Hence \(r_0=a\).
The preconditioners \(P,Q\) are invertible over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\), so \(A\) and \(A'=PAQ\) have the same Smith multiplicities over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\), and therefore the same multiplicities as over \({\mathsf{R}}\). Step A5 ensures that the leading block \(B\) is invertible over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\). By Lemma 2, \[A'\sim \begin{bmatrix}B&0\\0&J\end{bmatrix}, \qquad J=E-DB^{-1}C,\]
and \(J=pT\) for a matrix \(T\in {\mathsf{F}}^{k\times k}\). The block \(B\) contributes exactly \(a\) unit invariant factors. By Lemma 3, the block \(J=pT\) contributes \(\operatorname{rank}_{\mathsf{F}}(T)\) invariant factors equal to \(p\) and \(k-\operatorname{rank}_{\mathsf{F}}(T)\) invariant factors equal to \(0\). Step A8 computes this rank, so the returned multiplicities are correct. 0◻ ◻
The complexity is established as follows. Let \(\omega\) be a feasible exponent for dense matrix multiplication. Let \(\theta=1/(\omega-1)\) and \(\tau=3-\theta=3-1/(\omega-1)\). Note that for \(2<\omega<3\), we have \(2<\tau<\omega\).
A block application \(X\mapsto TX\) with \(X\in {\mathsf{F}}^{k\times t}\) uses:
a constant number of block products by \(A'\) over \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\) with \(t\) columns, costing \({O\mskip 1mu\tilde{\,}\mskip 1mu}(nt)\); this includes the products used to form \(C\widetilde{X}\), to compute the Hensel residual \(U-B\widetilde{Y}_0\), and to form the final vector \(A'[-Y;\widetilde{X}]\);
two solves with \(B_0\in {\mathsf{F}}^{a\times a}\) and \(t\) right hand sides;
sparse-block projection and matrix-polynomial work inside these solves.
A product by \(B_0\) with \(t\) columns costs \({O\mskip 1mu\tilde{\,}\mskip 1mu}(nt)\): embed the \(a\times t\) block into the first \(a\) coordinates of an \(n\times t\) block, apply \(A'\bmod p\), and restrict to the first \(a\) coordinates. Therefore, by Corollary 1, one solve costs \[{O\mskip 1mu\tilde{\,}\mskip 1mu}\left(\frac{a}{t}\,nt+a^2+a t^{\omega-1}\right) = {O\mskip 1mu\tilde{\,}\mskip 1mu}\left(an+a^2+a t^{\omega-1}\right).\] Thus one \(T\)-block application costs \[{O\mskip 1mu\tilde{\,}\mskip 1mu}\left(nt+an+a^2+a t^{\omega-1}\right).\] Since \(a\le n\), the \(a^2\) term is bounded by \(an\), but we keep it visible for the derivation.
The rank computation for \(T\in {\mathsf{F}}^{k\times k}\) is the solve-based binary-search rank computation of Theorem 3. Each nonsingularity test is applied to a leading submatrix of a fast preconditioning of \(T\). A block product by such a tested submatrix is obtained by embedding into the \(k\) coordinates, applying the fast preconditioners, applying the \(T\)-oracle, and restricting. The rank-by-binary-search proof uses \(O(k/t)\) such \(t\)-column oracle calls, plus \({O\mskip 1mu\tilde{\,}\mskip 1mu}(k^2+k t^{\omega-1})\) additional finite-field work. Using the \(T\)-block cost just derived, the cost after the Schur preconditioning is \[\begin{align} C(a,k,t) & = {O\mskip 1mu\tilde{\,}\mskip 1mu}\left( \frac{k}{t}\left(nt+an+a^2+a t^{\omega-1}\right) +k^2+k t^{\omega-1} \right) \\ & = {O\mskip 1mu\tilde{\,}\mskip 1mu}\left( kn+\frac{akn}{t}+\frac{a^2k}{t} +ak t^{\omega-2}+k^2+k t^{\omega-1} \right). \end{align}\] Because \(a\le n\), \[\frac{a^2k}{t}\le \frac{akn}{t}.\] Thus \[C(a,k,t) = {O\mskip 1mu\tilde{\,}\mskip 1mu}\left( kn+\frac{akn}{t}+ak t^{\omega-2}+k^2+k t^{\omega-1} \right).\] The initial rank computation, the verification of \(B_0\), and the setup/profile-check work for the fast preconditioners add \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n^2)\) operations, which is below the final bound. Applications of the internal conditioners during finite-field rank and solve routines are already included in the vector-product and \(m^2\) sparse-projection accounting of Section 4.
Choose \[t=\min\left\{a,k,\left\lceil n^\theta\right\rceil\right\}, \qquad \theta=1/(\omega-1).\] We prove that the cost is \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n^\tau)\) uniformly over all ranks \(a\).
Lemma 8 (Uniform bound). For all \(a,k\ge0\) with \(a+k=n\), and with \(t\) chosen as above when \(a,k>0\), \[C(a,k,t)={O\mskip 1mu\tilde{\,}\mskip 1mu}(n^\tau), \qquad \tau=3-1/(\omega-1).\]
Proof. The terms \(kn\), \(k^2\), and the initial \(n^2\) overhead are all \(O(n^2)\). Also \[k t^{\omega-1}\le {O\mskip 1mu\tilde{\,}\mskip 1mu}\left(n\cdot (n^\theta)^{\omega-1}\right)={O\mskip 1mu\tilde{\,}\mskip 1mu}(n^2).\] It remains to bound \[\frac{akn}{t} \quad\text{and}\quad ak t^{\omega-2}.\]
First suppose \(\min\{a,k\}\ge n^\theta\). Then \(t=\Theta(n^\theta)\), and since \(ak\le n^2\), \[\frac{akn}{t}\le {O\mskip 1mu\tilde{\,}\mskip 1mu}(n^{3-\theta})={O\mskip 1mu\tilde{\,}\mskip 1mu}(n^\tau).\] Similarly, \[ak t^{\omega-2} \le {O\mskip 1mu\tilde{\,}\mskip 1mu}(n^2\cdot n^{\theta(\omega-2)}) ={O\mskip 1mu\tilde{\,}\mskip 1mu}\left(n^{2+(\omega-2)/(\omega-1)}\right) ={O\mskip 1mu\tilde{\,}\mskip 1mu}\left(n^{3-1/(\omega-1)}\right) ={O\mskip 1mu\tilde{\,}\mskip 1mu}(n^\tau).\]
Now suppose \(m=\min\{a,k\}<n^\theta\). Then \(t=m\). Since \(ak/m=\max\{a,k\}\le n\), \[\frac{akn}{t}=\frac{akn}{m}\le n^2.\] Also \[ak t^{\omega-2}=ak m^{\omega-2} =\max\{a,k\}\, m^{\omega-1} \le {O\mskip 1mu\tilde{\,}\mskip 1mu}\left(n\,(n^\theta)^{\omega-1}\right) ={O\mskip 1mu\tilde{\,}\mskip 1mu}(n^2).\] Since \(\tau>2\) for \(\omega>2\), these terms are also \(O(n^\tau)\). 0◻ ◻
The overall cost is summarized as follows.
Theorem 5 (Cost of computing the Smith multiplicities). Using the sparse-block algorithm summarized in Section 4 above and the fast rank-profile preconditioner of Lemma 5, suppose the product \(v\mapsto Av\) costs \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n)\) operations. An \(n\times s\) product is obtained by applying this oracle to each column, and costs \({O\mskip 1mu\tilde{\,}\mskip 1mu}(ns)\).
Then the Smith multiplicities of any \(A\in{\mathsf{R}}^{n\times n}\) can be computed by a randomized algorithm which returns the correct multiplicities with probability at least \(1-\epsilon\) and whose expected number of operations over the working ring \(\mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}\) is \[{O\mskip 1mu\tilde{\,}\mskip 1mu}\left(n^{3-1/(\omega-1)}\right).\] For fixed or word-sized \(p\), and with \([ \mathchoice {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\mathsf{R}}}\vphantom{\mathsf{R}}} {\smash{\tilde{\scriptstyle\mathsf{R}}}\vphantom{\scriptstyle\mathsf{R}}} {\smash{\tilde{\scriptscriptstyle\mathsf{R}}}\vphantom{\scriptscriptstyle\mathsf{R}}}:{\mathsf{R}}]=O(\log n+\log(1/\epsilon))\), this is also \[{O\mskip 1mu\tilde{\,}\mskip 1mu}\left(n^{3-1/(\omega-1)}\right)\] operations over the base ring \({\mathsf{R}}\).
Proof. Correctness follows from Theorem 4. The cost of the main rank computation is bounded by Lemma 8. The initial rank computations, the \(B_0\) check, and the setup/profile-check work for the global preconditioners cost \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n^2)\). Applications of global and internal conditioners inside the finite-field rank and solve routines are part of the vector-product and sparse-projections in Section 4. The field extension and failure-probability amplification add only polylogarithmic factors hidden in the soft-O. 0◻ ◻
The Schur/divide identity used above can be extended to higher powers of \(p\). Let \(A\) be a matrix over \({\mathbb{Z}}/p^e{\mathbb{Z}}\), or over the corresponding unramified Galois ring of length \(e\). After row and column preconditioning, suppose \[A= \begin{bmatrix} B & C\\ D & E \end{bmatrix}, \qquad B\bmod p \text{ nonsingular},\] where the size of \(B\) is \(\operatorname{rank}(A\bmod p)\). Then \(B\) is invertible and \(A\) is equivalent to \[\begin{bmatrix} B & 0\\ 0 & J \end{bmatrix}, \qquad J=E-DB^{-1}C .\] Moreover \(J\equiv 0\bmod p\), so we can divide \(J\) by \(p\) and continue over the shorter ring \({\mathbb{Z}}/p^{e-1}{\mathbb{Z}}\). Thus the same algebra gives a recursive description of the Smith multiplicities. This recursion does not carry over to the complexity result of this paper. For \(e=2\) there is only one divided Schur complement to rank, and its block products are obtained through a single layer of carry-safe solves with the unit block \(B\). For \(e\ge3\), the next divided Schur complement already has products that are themselves defined through lower-level solves. A direct black-box recursion therefore nests the rank and solve procedures, rather than preserving the single-level amortization used above.
Consequently, the straightforward recursive use of the present method does not give a worst-case bound below the dense \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n^\omega)\) benchmark for higher powers of \(p\). The results of this paper are therefore stated only for length-two rings. The corresponding sub-\(n^\omega\) sparse SNF problem for \({\mathbb{Z}}/p^e{\mathbb{Z}}\), \(e\geq 3\), remains open.
We have combined a number of modern iterative block matrix tools with rank-profile Schur splitting over \({\mathsf{R}}={\mathbb{Z}}/p^2{\mathbb{Z}}\) and carry-safe Hensel solves to create a new Monte Carlo algorithm for sparse Smith form modulo \(p^2\) which requires an expected \[{O\mskip 1mu\tilde{\,}\mskip 1mu}(n^{3-\frac{1}{\omega-1}})operations in {\mathsf{R}}.\] When using standard cubic matrix arithmetic with \(\omega=3\) this cost is \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n^{2.5})\) operations in \({\mathsf{R}}\), whereas with the best currently known \(\omega\approx 2.371339\) our algorithm requires \({O\mskip 1mu\tilde{\,}\mskip 1mu}(n^{2.270786})\) operations in \({\mathsf{R}}\), which has exponent well below \(\omega\).
The key points are that the Schur complement is computed only through a block oracle, all divisions by \(p\) occur only after divisibility is established in the ring, and sparse-block projections are used inside the EGGSV solve tests so that the projection cost is \({O\mskip 1mu\tilde{\,}\mskip 1mu}(m^2)\) rather than \({O\mskip 1mu\tilde{\,}\mskip 1mu}(m^2t)\).
For matrices modulo \(p^e\) the same Schur/divide approach is valid, but the simple recursive sparse-block complexity analysis hits an \(n^\omega\) barrier already at \(e=3\). Getting past this appears to require a new multilevel blocking or filtered-module Krylov idea, not merely a black-box recursive invocation of the \(p^2\) algorithm.
The algorithm presented here is of the Monte Carlo type, with a controllably small probability of error. It would be interesting to find a Las Vegas algorithm with this same complexity, perhaps following the certified rank techniques of [11].
The author would like to thank the anonymous referees for their input. The author acknowledges support of the Natural Sciences and Engineering Research Council (NSERC), Canada.