Solving p-adic polynomial systems via iterative eigenvector algorithms

TL;DR

This paper presents a method for solving 0-dimensional p-adic polynomial systems using iterative eigenvector algorithms, along with an improved eigenvalue/eigenvector computation for p-adic matrices.

Contribution

It introduces an implementation of a polynomial system solver for p-adic systems and enhances an existing eigenvalue algorithm for p-adic matrices.

Findings

Effective approximate solutions for p-adic polynomial systems.

Improved eigenvalue and eigenvector computation for p-adic matrices.

Enhanced algorithm performance in finite precision p-adic arithmetic.

Abstract

In this article, we describe an implementation of a polynomial system solver to compute the approximate solutions of a 0-dimensional polynomial system with finite precision p-adic arithmetic. We also describe an improvement to an algorithm of Caruso, Roe, and Vaccon for calculating the eigenvalues and eigenvectors of a p-adic matrix.

Equations

$$\Vert v{\Vert }_p:=\underset{1\le i\le n}{\sup }\{{\left|v_i\right|}_p\}.$$

$$\Vert A{\Vert }_p:=\underset{v\in {\mathbb{Q}}_p^m\backslash \{0\}}{\sup }\frac{\Vert Av{\Vert }_p}{\Vert v{\Vert }_p}.$$

$$a=a_{-r}{p}^{-r}+\dots +a_0+p{a}_1+p^2{a}_2+\dots +a_{N-1}{p}^{N-1}+O(p^N)$$

$$\frac{(1+p^{99}+O(p^{100}))-(1+O(p^{100}))}{p^{100}+O(p^{200})}=p^{-1}+O(1)$$

$$A:=\left[\begin{array}{cc}{p}^3& 1\\ 0& -p^3\end{array}\right],B:=\left[\begin{array}{cc}{p}^3& 1\\ p^6& -p^3\end{array}\right],$$

$$v_A:=\left[\begin{array}{c}1\\ 0\end{array}\right],v_B:=\left[\begin{array}{c}1\\ p^3\sqrt{2}\end{array}\right].$$

$$\left|{\lambda }_1\dots {\lambda }_{n-m}\right|={\left|{\sigma }_1\dots {\sigma }_{n-m}\right|}_p=1\text{ and }{\left|{\lambda }_{n-m}\right|}_p>{\left|{\lambda }_{n-m+1}\right|}_p.$$

$$(B^m{)}^N=\left[\begin{array}{cc}{B}_{11}^{Nm}& \ast \\ 0& 0\end{array}\right]+O(p^N).$$

$$(A-\lambda I{)}^{Nm}(AM)=A(A-\lambda I{)}^{Nm}V=0$$

$$A\mathbf{V}=\mathbf{V}\left[\begin{array}{cccc}{X}_1& & & \\ & X_2& & \\ & & \ddots & \\ & & & X_r\end{array}\right]+O(p^N).$$

$$O\left(\rho m^3+\ell (n^3+M(n){\log }_2(mN))\right).$$

$$H:=I-\frac{2}{v^{T}v}v{v}^T$$

$$p^{2r}{v}^{T}v$$

$$=2(x^{T}x-(x_i+O(p^{r+1}))x_1)$$

$$=2(x_i^2-x_i{x}_1)+O(p^{2r+1})$$

$$=2x_i(x_i-x_1)+O(p^{2r+1})$$

$$HH$$

$$=I-\frac{4}{v^{T}v}v{v}^T+\frac{4}{(v^{T}v{)}^2}(v{v}^T)(v{v}^T)$$

$$=I-\frac{4}{v^{T}v}v{v}^T+\frac{4}{(v^{T}v)}(v{v}^T)$$

$$=I.$$

$$0\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}$$

$$\lx@xy@svg{\hbox{\raise 0.0pt\hbox{\kern 5.5pt\hbox{\ignorespaces\ignorespaces\ignorespaces\hbox{\vtop{\kern 0.0pt\offinterlineskip\halign{\entry@#!@&&\entry@@#!@\cr&&&&\crcr}}}\ignorespaces{\hbox{\kern-5.5pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{0\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{\hbox{\kern 29.5pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 29.5pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{I\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{\hbox{\kern 64.68054pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 64.68054pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{R\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 89.01181pt\raise 4.50694pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern 0.0pt\raise-1.50694pt\hbox{$\scriptstyle{\pi}$}}}\kern 3.0pt}}}}}}\ignorespaces{\hbox{\kern 102.35068pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 102.35068pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{R/I\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{\hbox{\kern 154.64572pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 154.64572pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{0}$}}}}}}}\ignorespaces}}}}\ignorespaces.$$

$$B\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}$$

$$0\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}$$

$$B\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}$$

Full text

1. Introduction

Let $k$ be a field and let $f_{1},\ldots,f_{m}\in k[x_{1},\ldots,x_{n}]$ be polynomials such that the variety defined by $f_{1},\ldots,f_{m}$ has dimension [math]. It is often of interest to compute the solutions of such a system, either exactly or approximately. A significant bottleneck in computing the exact solutions of such a polynomial system is the complexity of the field extension required to write down all of the solutions. A popular method to compute the solutions exactly is to compute a triangular decomposition for the ideal $\langle f_{1},\ldots,f_{m}\rangle$ and solve for the coordinates via back-substitution. A difficulty with this approach is the explosion of the size of the coefficients, which in conjunction with taking complicated field extensions, usually renders a computer algebra system unresponsive.

Computations in an exact field (such as $\mathbb{Q}$) can often be approximated by computations in an inexact field (such as $\mathbb{R},\mathbb{C}$), where by inexact, we mean that the computer representation of an element is only an approximation up to some precision. It is often far faster to compute an approximation to the solutions, though the trade-off is that the solutions are not known perfectly. In numerous applications, a real or complex approximation to the solutions is sufficient.

In Linear Algebra, $p$-adic methods go back to at least Dixon [dixon1982exact]. One advantage pointed out by Dixon is that a linear system can be ill-conditioned with respect to an archimedean norm, but well-conditioned $p$-adically; general polynomial systems exhibit this behavior as well. Recently, new developments in number theory and tropical geometry highlight that non-archimedean data is not just a route to integral solutions, but is significant in and of itself.

Our polynomial system solver is based on the truncated normal form solver of [TMV2018truncated]. To transplant their implementation to the $p$-adic setting, we rely on finite precision linear algebra over the $p$-adic field, which in analogy to the case over $\mathbb{R}$ we refer to as pnumerical linear algebra. It is essential for us to be able to compute the eigenvalues and eigenvectors of an approximate $p$-adic matrix quickly and stably. The naive algorithm to compute the eigenvalues of a matrix relies on the calculation and factorization of the characteristic polynomial. In the numerical setting (i.e, over $\mathbb{R}$) it is desirable to avoid this step due to the instability of solving for the roots of a polynomial. There is a similar loss of precision in solving for the roots of a $p$-adic polynomial when the roots of the polynomial are $p$-adically close together. Unlike the archimedean setting, the accurate calculation of the characteristic polynomial is more difficult since a division free algorithm must be used [CRV2017characteristic, Introduction]. Our idea to improve on the algorithm of [CRV2017characteristic] is to adapt the ideas from classical numerical linear algebra to the pnumerical setting; specifically to use an iterative scheme to compute the eigenvectors or eigenvalues. To our knowledge, this is the first appearance of a $p$-adic numerical algorithm based on iterating matrix multiplication to solve for eigenvalues and eigenvectors.

The significant deviation from the classical numerical setting is the preconditioning strategy. A general matrix over $\mathbb{R}$ tends to have eigenvalues of distinct absolute values, so aspects of the iteration are dominated by the unique largest eigenvalue. In the $p$-adic setting, it tends to be the case that all of the eigenvalues have identical absolute value, and this is true for most shifts of the form $A\mapsto A-\mu I$ as well. Our strategy is to first approximate the eigenvalues by computing the characteristic polynomial over the residue field, then refine the initial approximation via iteration. The advantage of this approach is that the computation of the characteristic polynomial of a matrix in $M_{n}(\mathbb{F}_{q})$ is significantly faster than computing the characteristic polynomial over $\mathbb{Q}_{p}$ [CRV2017characteristic, Storjohann2001frobenius]. The basic Hensel-lifting based strategy applied to the system of equations $Av=\lambda v,{v_{i}=1}$, with $v_{i}$ some fixed entry of $v$, would require the computation of the inverse of the Jacobian matrix of the map $F\colon(v,\lambda)\mapsto(Av-\lambda v,v_{i}-1)$ at each step in the iteration. Our cost per iteration step is much smaller. Unfortunately, matrices of the form $\mu I+N$, with $N\in M_{n}(\mathbb{Z}_{p})$ topologically nilpotent, do not allow us to apply our deflation strategy as the characteristic polynomial over the residue field is unhelpful; we hope that future improvements to our strategy can be made.

The structure of the article is as follows. In Section 2 we discuss some background and terminology regarding the subject of finite precision computation over $\mathbb{Q}_{p}$. We review some specific results on linear algebra over the $p$-adic field, especially the $p$-adic analogues of condition-stable algorithms based on matrix factorizations. We describe our iterative strategy for the computation of the eigenvectors of a matrix over $\mathbb{Q}_{p}$ and compare it to the existing alternatives. In Section 3, we describe our adaptations of the algorithm of [TMV2018truncated] to solve polynomial systems over $\mathbb{Q}_{p}$. Our algorithms are available as Julia packages, and can be obtained from:

[TABLE]

(For technical reasons related to the Julia package system, we choose to divide our implementation between two packages.)

2. Linear algebra

2.1. Precision, norms, and condition numbers

We state some basic definitions for our discourse. We direct the reader to [CRV2015linear, caruso2017computationswithpadics, Kedlaya2010differential] for more details. For a matrix $A\in M_{n}(\mathbb{Z}_{p})$, we will denote by $\chi_{A}$ the characteristic polynomial and denote $\chi_{A,p}:=\chi_{A}\pmod{p}$. We denote the standard basis of $\mathbb{Q}_{p}^{n}$ (or $\mathbb{Z}_{p}^{n}$) by $\{e_{1},\ldots,e_{n}\}$.

Definition 2.1**.**

Let $v:=(v_{1},\ldots,v_{n})\in\mathbb{Q}_{p}^{n}$ be a vector. The norm of $v$ is

[TABLE]

Definition 2.2**.**

Let $A\in\mathbb{Q}_{p}^{n\times m}$ be a matrix. The operator norm of $A$ with respect to $\lVert\cdot\rVert_{p}$ is

[TABLE]

The condition number for an invertible square matrix is $\kappa(A):=\lVert A\rVert_{p}\cdot\lVert A^{-1}\rVert_{p}$. The condition number for a singular matrix is $\infty$.

We can identify a subgroup of $\operatorname{GL}_{n}(\mathbb{Q}_{p})$ where every matrix is well-conditioned, serving the analogous role to $\operatorname{O}_{n}(\mathbb{R})$ in the real setting.

Lemma 2.3**.**

Let $A\in\operatorname{GL}_{n}(\mathbb{Z}_{p})$. Then $\kappa(A)=1$.

Proof.

See the comments following [Kedlaya2010differential, Definition 4.3.3]. ∎

We now come to the discussion of $p$-adic precision. There are many possible ways to represent a $p$-adic element $a\in\mathbb{Q}_{p}$ in a computer system [caruso2017computationswithpadics]. In the Nemo/Hecke systems [nemohecke], an element is represented by a series

[TABLE]

where the $O(p^{N})$ is the $p$-adic ball representing the uncertainty of the remaining digits. The relative precision of $a$ is the quantity $N+r$, and the absolute precision is the number $N$. In the terminology of [caruso2017computationswithpadics], we consider a system with the zealous (i.e, interval) implementation of arithmetic. The operations $-,+$ preserve the minimum of the absolute precision of the operands, and $\times,\div$ preserve the minimum relative precision of the operands. If $u\in\mathbb{Z}_{p}^{\times}$, $a\in\mathbb{Z}_{p}$, and $N\leq N^{\prime}$, then we have that $(u+O(p^{N^{\prime}}))(a+O(p^{N}))=ua+O(p^{N})$. Multiplication by $p$ preserves the relative precision and increases the absolute precision by $1$. The worst operation when it comes to absolute $p$-adic precision is dividing a small number by $p$. For example, the expression

[TABLE]

begins with $3$ numbers with an absolute and relative precision of at least $100$, and ends with a result where not even the constant term is known.

In Nemo/Hecke, elements of a matrix store their own precision. In the terminology of [CRV2017characteristic], the associated precision structure is the jagged precision lattice. However, we will not keep track of this finer precision here. Instead, in the terminology of [CRV2017characteristic], we will look at the flat precision of the matrix:

Definition 2.4**.**

Let $A,B\in M_{n}(\mathbb{Z}_{p})$ be matrices such that $a_{i,j}=b_{i,j}+O(p^{N_{i,j}})$. Then we write $A=B+O(p^{N})$, where $N:=\min_{i,j}N_{i,j}$.

To refer to a matrix $A\in M_{n}(\mathbb{Q}_{p})$ whose elements are known at an absolute precision at least $N$, we will simply write $A+O(p^{N})$.

2.1.1. Design choices for the implementation

We discuss some of the design choices of the computer algebra package accompanying this article. The user we had in mind while designing our system is one who is interested in data that is inherently $p$-adic. More specifically, we have assumed the following hypotheses:

Assumptions**.**

• •

Input and output are inherently approximate.

• •

The desired output precision is roughly the input precision.

• •

Element-wise gains in absolute precision after $a\mapsto pa$ or $a\mapsto a^{p}$ are usually lost in a subsequent step.

We believe that our linear algebra package will be used in the middle of an ongoing $p$-adic computation, such as solving polynomial systems of equations. In such a circumstance, it does not make sense for our algorithm to “request” extra digits of precision from the input data since the input itself is approximate. We also assume that if a user inputs data at precision $N$, they expect output at precision roughly $N$, or at least $N$ minus the condition number for the problem. Computations using relaxed $p$-adic numbers [caruso2017computationswithpadics] allow the desired output precision to be specified. Unfortunately, relaxed arithmetic does not seem well-suited for use in an iterative scheme due to the extra memory costs (see [caruso2017computationswithpadics, Section 2.4]). We also assume for a dense matrix that the entries are all given at the same flat absolute precision. We believe that this situation is fairly common in practical applications, but of course not ubiquitous. That said, the implementation of $p$-adic matrices in [nemohecke] uses an element-wise precision model, so it is possible for a computation using our software to validly return a result with more accuracy than we indicate here.

In our analysis, we do not account for the extra precision on an element $a\in\mathbb{Q}_{p}$ gained after $a\mapsto pa$ or $a\mapsto a^{p}$. Experimentally, this seems to be a reasonable assumption for randomly selected dense matrices over a $p$-adic field, with $p\geq 40$. For a more structured scenario this assumption might not be reasonable. A nice feature of capped precision is that we do not have to account for the growth of arithmetic costs due to increased precision.

With these assumptions made, it is often useful to think of the computer as performing arithmetic in a $\mathbb{Z}/p^{N}\mathbb{Z}$-module for a computation with no divisions by $p$. It is useful to consider the image of a matrix $A\in M_{n}(\mathbb{Z}_{p})$ in the ring $M_{n}(\mathbb{Z}/p^{N}\mathbb{Z})$. It is an essential difficulty of finite precision $p$-adic linear algebra that $\mathbb{Z}/p^{N}\mathbb{Z}$ is a local principal ideal ring, but not a domain when $N>1$.

Remark 2.5**.**

Our eigenvector algorithms are focused on finding only the eigenvectors defined over $\mathbb{Q}_{p}$. We made this choice for the following two reasons. First, in the last section of the article, we focus on finding the $\mathbb{Q}_{p}$-solutions of a [math]-dimensional system of polynomial equations over $\mathbb{Q}_{p}$. Secondly, our present implementation is written using components of the OSCAR system (specifically [nemohecke]) that are currently under development. At the time of writing the first version of our implementation, the “$q$-adics” interface or an interface to handle ramifed extensions of $\mathbb{Q}_{p}$ did not exist. Thus, our article is written with $\mathbb{Q}_{p}$ in mind to more closely reflect the actual software. Rapid progress is being made on the system, so we expect our implementation to evolve over time as well.

2.2. Matrix factorizations

In this section, we very briefly review some matrix invariants and factorizations. Matrix factorizations relating to the topological structure of $\mathbb{Q}_{p}$ are unsurprisingly related to the algebraic structure of the $\mathbb{Z}_{p}$-module spanned by the columns (or rows) of the matrix. We invite the reader to consider [Kedlaya2010differential, Chapter 4] for further details.

Proposition 2.6** (Iwasawa decomposition).**

Let $k$ be either a finite extension of $\mathbb{Q}_{p}$, $\mathbb{R}$, or $\mathbb{C}$ and let $G\lhook\joinrel\relbar\joinrel\rightarrow\operatorname{GL}_{n}(k)$ be a linear algebraic group, let $K$ be a maximal compact subgroup of $G$, and let $B$ be a Borel subgroup. For any $A\in G$, there exists a $Q\in K$ and an $R\in B$ such that $A=QR$.

Note that a maximal compact subgroup of $\operatorname{GL}_{n}(\mathbb{R})$ is the orthogonal group $\operatorname{O}_{n}(\mathbb{R})$, and the subgroup of invertible real upper triangular matrices is a Borel subgroup. We have as a consequence:

Corollary 2.7** ($QR$-factorization).**

Let $A\in\mathbb{R}^{n\times m}$ be a matrix. Then there exists a $Q\in\operatorname{O}_{n}(\mathbb{R})$ and an upper triangular matrix $R\in\mathbb{R}^{n\times m}$ such that $A=QR$.

Of course, $QR$-factorization is really just a consequence of Gram-Schmidt. From the point of view of Proposition 2.6, we also have a $p$-adic analogue of $QR$-factorization. In $\operatorname{GL}_{n}(\mathbb{Q}_{p})$, a maximal compact subgroup is $\operatorname{GL}_{n}(\mathbb{Z}_{p})$, and the subgroup of invertible upper triangular matrices is a Borel subgroup.

Corollary 2.8** ($p$-adic $QR$-factorization).**

Let $A\in\mathbb{Z}_{p}^{n\times m}$ be a matrix. Then there exists a $Q\in\operatorname{GL}_{n}(\mathbb{Z}_{p})$ and an upper triangular matrix $R\in\mathbb{Z}_{p}^{n\times m}$ such that $A=QR$.

We decided to refer to the factorization above as a $QR$-factorization since many applications are similar to the real case. However, one important difference is that $Q$ is generally not an element of $\operatorname{O}_{n}(\mathbb{Q}_{p})$. As in the case over $\mathbb{R}$, an elementary proof of Corollary 2.8 can be given by exhibiting an algorithm to compute a $QR$-decomposition. The algorithm is well-known, and in fact it is just the algorithm to compute a $PLU$-decomposition, though with each row-pivot chosen by taking the vector with the largest $p$-adic norm [Kedlaya2010differential, Chapter 4]. We also note since $Q\in\operatorname{GL}_{n}(\mathbb{Z}_{p})$, that $R$ can be chosen to be the Hermite normal form of $A$. We summarize the discussion as:

Remark 2.9**.**

The $p$-adic $QR$-decomposition of $A$ is computed using a $PLU$-decomposition, with pivots chosen by $p$-adic norm.

An essential property of the $QR$-factorization of a real matrix $A$ is that $Q$ has a well-conditioned inverse. The same is true $p$-adically. The columns of $Q$ in the $p$-adic case are also “orthogonal” in the sense of [Schikhof1984ultrametric, Section 50]. Unfortunately, this notion of orthogonality does not help to compute the inverse of $Q$ like in the real case. As we might expect, $QR$-factorizations give us a way to compute the singular value decomposition of a matrix in the $p$-adic setting.

Proposition 2.10**.**

If $A\in M_{n}(\mathbb{Q}_{p})$, then there exist $U,V\in\operatorname{GL}_{n}(\mathbb{Z}_{p})$ and a diagonal matrix $\Sigma\in M_{n}(\mathbb{Q}_{p})$ such that $A=U\Sigma V^{T}$.

Proof.

Let $A=URP^{\prime}$ be a $QR$-decomposition with column pivoting. Now let $R^{T}=V\Sigma$ be a $QR$-decomposition. The result follows. ∎

The singular value decomposition is the Smith normal form of the $\mathbb{Z}_{p}$-module generated by columns of $A$ [Kedlaya2010differential, Chapter 4]. Since the matrices $U,V$ are well-conditioned, we can use a singular value decomposition to stabilize pnumerical computations (as in the real setting).

2.3. The eigenvector problem: semantics

Let $A\in M_{n}(\mathbb{Z}_{p})$ be a matrix. In the setting of a finite $p$-adic precision computation, we may formulate two versions of the eigenvector problem:

Problem 1** (Eigenvector problem, Version I).**

Let $(\lambda,v)$ be an exact eigenpair for $A$ over $\mathbb{Z}_{p}$. Given an approximation $\widetilde{A}$ to $A$ such that $A=\widetilde{A}+O(p^{N})$, compute a pair $(\widetilde{\lambda},\widetilde{v})$ such that $\lVert\widetilde{v}\rVert=1$ and $\widetilde{A}\widetilde{v}=\widetilde{\lambda}\widetilde{v}+O(p^{N})$.

Problem 2** (Eigenvector problem, Version II).**

Let $(\lambda,v)$ be an exact eigenpair for $A$ over $\mathbb{Z}_{p}$. Given an approximation $\widetilde{A}$ to $A$ such that $A=\widetilde{A}+O(p^{N})$, compute a triple $(\widetilde{\lambda},\widetilde{v},M)\in\mathbb{Z}_{p}\times\mathbb{Z}_{p}^{n}\times\mathbb{Z}$ such that:

1. (i)

$\lVert\widetilde{v}\rVert=1$, $\lambda=\widetilde{\lambda}+O(p^{M})$ and $v=\widetilde{v}+O(p^{M})$, 2. (ii)

for any $B=0+O(p^{N})$, we have that $(\tilde{A}+B)\widetilde{v}=\widetilde{\lambda}\widetilde{v}+O(p^{M})$, 3. (iii)

$M$ is maximal among such triples.

In Version II, the second criteria ensures that the problem is well-posed.

Lemma 2.11**.**

Let $A\in M_{n}(\mathbb{Z}_{p})$, let $\widetilde{A}$ be an approximation such that $A=\widetilde{A}+O(p^{N})$, and let $(\lambda,v,M)$ satisfy the first two criteria for a solution to Version II. Then $M\leq N$.

Proof.

If $\lVert(A-\widetilde{A})\widetilde{v}\rVert_{p}\geq p^{-N}$ then clearly $M\leq N$. Otherwise, choose any $B\in M_{n}(\mathbb{Z}_{p})$ such that $\lVert B\widetilde{v}\rVert_{p}=1$. Then $\lVert(A-\widetilde{A}+p^{N}B)\widetilde{v}\rVert_{p}=p^{-N}\lVert B\widetilde{v}\rVert_{p}=p^{-N}$. ∎

The following example demonstrates the difference between Version I and Version II of the problem.

Example 2.12**.**

Let $p=7$, and consider the matrices

[TABLE]

Both matrices are rounded to the same result with precision $N:=6$. However, the respective eigenvectors are

[TABLE]

Note in $\mathbb{Z}_{7}$ that $\sqrt{2}=\pm 4+O(p)$. We see that the equation $Av_{A}+O(p^{6})=p^{3}v_{A}+O(p^{6})=Bv_{A}+O(p^{6})$. However, $v_{A}\neq v_{B}+O(p^{6})$.

Version II exhibits pathological behavior whenever the characteristic polynomial of $A$ is not square-free modulo $p^{M}$. For instance, as $X$ ranges over all $\mathbb{Z}_{p}^{n}$ matrices, we see that the eigenvectors of $I+p^{M}X$ could be anything. In other words, any solution to Version II of the problem for $I+p^{M}X$ has $M\leq 0$. A similar difficulty arises in the eigenvector problem for real matrices. Because of this pathology, we restrict our attention to solving Version I.

Of course, we can give analogous formulations for computing a Jordan form of $A$ or for computing a block Schur form for $A$. Generally, we will work with Version I of these problems as well.

2.4. Commentary on the classical algorithm

The naive algorithm to compute the eigenvalues of a matrix relies on the calculation and factorization of the characteristic polynomial. In the numerical setting over $\mathbb{R}$ it is desirable to avoid this step due to the instability of solving for the roots of a polynomial. Also $p$-adically, there is also a loss of precision in solving for the roots of a polynomial. Thus, in order to accurately compute the nullspace of $A-\lambda I$ the roots of the characteristic polynomial, and thus the characteristic polynomial itself, must be known to as much precision as possible. Presently, there are two algorithms to compute the characteristic polynomial at the maximum possible precision, both of which use more than $O(n^{3})$ arithmetic operations (asymptotically):

1. (1)

Use a division-free calculation. The asymptotic run-time is at least $O(n^{4})$. 2. (2)

Use the algorithm of [CRV2017characteristic]. Compute the characteristic polynomial of a matrix with inflated precision using division, then truncate to the optimal precision computed from the precision lattice analysis. The runtime is $O(\rho n^{3})$, where $\rho$ represents the factor of performing arithmetic at the higher precision. Practically, the computation of the comatrix and the precision increase required for the precision analysis appears to impose a large cost [CRV2017characteristic, Remark 5.3].

In the next subsection, we discuss how to avoid an expensive computation of the characteristic polynomial for most matrices in $M_{n}(\mathbb{Z}_{p})$.

Remark 2.13**.**

To clarify the term “most”, if $A+O(p^{N})$ is an $n\times n$-matrix whose entries are chosen with the uniform probability distribution on $[0,\ldots,p^{N}-1]$, then the limit as $n\rightarrow\infty$ of the probability that $\chi_{A}$ is square-free is at least $\frac{1-p^{-5}}{1+p^{-3}}$ [Fulman2002random].

2.5. The power iteration algorithm

Algorithm 2.14

Eigvecs($A$)

0:  $A+O(p^{N})$, an $n\times n$-matrix over $\mathbb{Z}_{p}$.

0:  The eigenpairs $(\lambda_{i},v_{i})$ of the matrix $A$.

1:  if $A$ is diagonal  then

2:     return  $\{(a_{11},e_{1}),\ldots,(a_{nn},e_{n})\}$.

3:  if $A\mod p\equiv 0$ then

4:     $\nu=\min_{i,j}\{\operatorname{ord}x_{ij}\}$

5:      $V$ = Eigvecs($p^{-v}X$)

6:      Reset precision of $V$

7:     return  $V$

8:  Compute $\chi_{A,p}$

9:  if $\chi_{A,p}$ has no linear factors then

10:     return  Nothing.

11:  if  $\chi_{A,p}\not\equiv(x-\lambda)^{n}\pmod{p}$ then

12:     $\{(X_{1},V_{1}),\ldots,(X_{r},V_{r})\}$ := PowerIterationDecomposition($A$, $\chi_{A,p}$)

13:     return  Eigvecs($X_{i}$) for $i=1,\ldots,r$

14:  else

15:     return   ClassicalAlgorithm($A$)

Note in step 6, we may reset the precision since the input matrix on the previous step $X$ was a multiple of $p^{\nu}$. If $p$ is moderately large, then for most matrices the size of the Jordan blocks modulo $p$ will be small (see Remark 2.13). The square matrix in step $13$ is expected to be small, so we apply the costly method to compute the characteristic polynomial of this block.

Algorithm 2.15

PowerIterationDecomposition($A$, $\chi_{A,p}$)

0:   $A+O(p^{N})$, an $n\times n$-matrix over $\mathbb{Z}_{p}$. $\chi_{A,p}$, the characteristic polynomial of $A\pmod{p}$.

0:  A list of pairs of matrices $\{(X_{i},V_{i})\}_{i}$ such that $AV_{i}=V_{i}X_{i}+O(p^{N})$.

1:  Compute eigenpairs $(\lambda_{i},V_{i})$ mod $p$

2:  Set $m$ to be the largest multiplicity of a linear factor of $\chi_{A,p}$.

3:  for $V=V_{1},\ldots,V_{r}$  do

4:      Lift the eigenpair to an approximate eigenpair $(\lambda_{i},V_{i})$ in $\mathbb{Z}_{p}$

5:     Set $B=(A-\lambda_{i}I)$.

6:     Iterate for $\log_{2}(mN)$ times, $B=B^{2}$

7:     Set $V$ := nullspace($B+O(p^{N})$).

8:      Solve $VX_{i}=AV$.

9:  return  $\{(X_{1},V_{1}),\ldots,(X_{r},V_{r})\}$.

To clarify, the lifting in step 4 is just the trivial operation of interpreting the elements $\{0,\ldots,p-1\}$ as elements of $\mathbb{Z}_{p}$.

2.6. Proof of correctness

We give a quick proof of correctness of the power iteration algorithm.

Lemma 2.16**.**

Let $\lambda$ be an approximate eigenvalue of $A\pmod{p}$, let $m$ be the multiplicity of $\lambda$ as a root of $\chi_{A,p}$, and let $M$ be the nullspace of $(A-\lambda I)^{Nm}$. Then $M$ is a free $\mathbb{Z}/p^{N}\mathbb{Z}$-module with trivial annihilator. Letting $V$ be a matrix whose columns are the generators of $M$, there exists an $X\in M_{N}(\mathbb{Z}_{p})$ with such that $AV=VX+O(p^{N})$.

Proof.

First, note that $(A-\lambda I)\pmod{p}$ has [math] as an eigenvalue with multiplicity $m$. Ordering the singular values and eigenvalues by size, we have that $(A-\lambda I)$ satisfies

[TABLE]

Thus, by the Hodge-Newton decomposition [Kedlaya2010differential, Theorem 4.3.11], there is a matrix $U\in\operatorname{GL}_{n}(\mathbb{Z}_{p})$ such that $B:=U^{-1}(A-\lambda I)U$ is block upper triangular, with the lower right block $B_{22}$ accounting for the $m$ small eigenvalues and the upper left block $B_{11}$ accounting for the unit eigenvalues. Now, we see that $B_{22}^{m}=0\pmod{p}$, so

[TABLE]

But as $U\in\operatorname{GL}_{n}(\mathbb{Z}_{p})$, we see that the Smith normal form of $(A-\lambda I)^{Nm}$ is $\begin{bmatrix}I&0\\ 0&0\end{bmatrix}.$ Thus, the kernel of $(A-\lambda I)^{Nm}$ is free as a $\mathbb{Z}/p^{N}\mathbb{Z}$-module with trivial annihilator in $\mathbb{Z}/p^{N}\mathbb{Z}$. Next, note that we have

[TABLE]

so $AM\subseteq M$. Thus $A$ restricted to $M$ defines an endomorphism of $M$, represented by a square matrix. ∎

As an immediate corollary, we have:

Corollary 2.17**.**

Let $(V_{1},X_{1}),\ldots,(V_{r},X_{r})$ be the blocks returned from power iteration. With $\mathbf{V}:=[V_{1},\ldots,V_{r}]$ the horizontal join, we have

[TABLE]

We hope that step 15 can be improved, as we believe it is helpful in solving polynomial systems whose solutions have large valuations. Thus, we raise the question:

Question 2.18**.**

Given a topologically nilpotent matrix $B\in M_{n}(\mathbb{Z}_{p})$, is there an efficient iterative strategy to determine the eigenvectors of $B+O(p^{N})$?

2.7. The Schur form algorithm

Solving for the nullspace at the end of the iteration is the dominant cost in the power iteration algorithm. In the classical numerical setting, we can avoid this drawback by using the $QR$-iteration algorithm. We also have a $p$-adic version of the $QR$-algorithm. With the $p$-adic $QR$ step replacing the usual $QR$ factorization, our algorithm is simply the original $LR$-algorithm [rutishauser1958eigenvalue] that inspired the $QR$-algorithm for real matrices. The $LR$-algorithm fell out of favor some time in the 1960’s since it is numerically less stable than the $QR$-algorithm, but $p$-adically this situation is reversed!

We now give the algorithm.

Algorithm 2.19

QR-Iteration($A$, $\chi_{A,p}$)

0:   $A+O(p^{N})$, an $n\times n$-matrix over $\mathbb{Z}_{p}$. $\chi_{A,p}$, the characteristic polynomial of $A\pmod{p}$.

0:  A (block) triangular form $T$ for $A$, and matrix $V$ such that $AV=VT+O(p^{N})$.

1:  Set $\lambda_{1},\ldots,\lambda_{\ell}$ to be the roots of $\chi_{A,p}$ in $\mathbb{F}_{p}$, lifted to $\mathbb{Z}_{p}$.

2:  Set $m_{1},\ldots,m_{\ell}$ to be the multiplicities of the roots of $\chi_{A,p}$.

3:  Compute $B,V$ such that $AV=VB$ and $B$ is in Hessenberg form.

4:  for $i=1,\ldots,\ell$ do

5:     for $j=1,\ldots,m_{i}N$ do

6:        Factor $(B-\lambda_{i}I)=QR$

7:        Set $B:=RQ+\lambda_{i}I$

8:        Set $V:=Q^{-1}V$

9:  return  $B,V$.

Algorithm 2.20

BlockSchurForm($A$)

0:  $A+O(p^{N})$, an $n\times n$-matrix over $\mathbb{Z}_{p}$.

0:  A block upper triangular matrix $T$ with $\ell+1$ distinct blocks and matrix $V$ such that $AV=VT+O(p^{N})$. Here, $\ell$ is the number of linear factors of $\chi_{A,p}=\chi_{A}\pmod{p}$.

1:  if $A$ is diagonal  then

2:     return  $A$, $\{e_{1},\ldots,e_{n}\}$.

3:  if $A\mod p\equiv 0$ then

4:     $\nu=\min_{i,j}\{\operatorname{ord}_{p}x_{ij}\}$

5:      $V$ = BlockSchurForm($p^{-v}X$)

6:      Reset precision of $V$

7:     return  $V$

8:  Compute $\chi_{A,p}:=\chi_{A}\pmod{p}$

9:  if $\chi_{A,p}$ has no linear factors then

10:     return  $A$, $\{e_{1},\ldots,e_{n}\}$

11:  if  $\chi_{A,p}\neq(x-\lambda)^{n}$ then

12:     return  $B,V$ := QR-Iteration($A$, $\chi_{A,p}$)

13:  else

14:     $B,V$ := $A$, $\{e_{1},\ldots,e_{n}\}$

15:  for $j:=1\ldots\ell$ do

16:     if  $\chi_{B_{ii},p}\equiv(x-\lambda)^{m}$ for $m>1$  then

17:         $C_{j},X_{j}$ := ClassicalAlgorithm($B_{j}$)

18:        Update $B$ and $V$

19:  return  $A,V$

We briefly comment on some aspects of this algorithm analogous to the archimedean setting (the proofs for the various statements are analogous as well). Note that each $p$-adic $QR$-iteration preserves the Hessenberg form, so the cost of applying the $QR$-step is $O(n^{2})$. From the discussion of the $LR$-algorithm in [wilkinson1965convergence, Sections 5-9], we see that the subdiagonal entry $(i,i+1)$ converges to [math] at a rate of $O\left(\frac{\left\lvert\lambda_{i+1}\right\rvert_{p}^{s}}{\left\lvert\lambda_{i}\right\rvert_{p}^{s}}\right)$, provided that $\left\lvert\lambda_{i+1}\right\rvert_{p}<\left\lvert\lambda_{i}\right\rvert_{p}$. For $B\in M_{n}(\mathbb{Z}_{p})$, if $\left\lvert\lambda_{i+1}\right\rvert_{p}<\left\lvert\lambda_{i}\right\rvert_{p}=1$, then $\frac{\left\lvert\lambda_{i+1}\right\rvert_{p}}{\left\lvert\lambda_{i}\right\rvert_{p}}<p^{-\frac{1}{m}}$, where $m$ is the largest multiplicity of an eigenvalue of $B\pmod{p}$. Additionally, if some subdiagonal entries are identified to be zero during the iteration, the standard deflation strategies can be used to accelerate the algorithm.

Remark 2.21**.**

In our version of the iteration we use a constant shift factor, and so the iteration converges in at most $mN$ steps [wilkinson1965convergence]. We believe that a $p$-adic analogue of the Rayleigh quotient strategy to refine the eigenvalue approximation will lead to a drastic speed-up in convergence.

Remark 2.22**.**

It is possible to use the $QR$-iteration to compute the valuations of the eigenvalues iteratively. After applying a single round of the $QR$-iteration to a Hessenberg matrix $A$ with no shift (i.e, $\lambda_{i}=0$), the valuations of the eigenvalues can be read from the diagonal blocks in the resulting matrix. If $B_{ii}$ is one of the diagonal blocks with no subdiagonal entry equal to [math], then all of the eigenvalues of $B_{ii}$ must have the same valuation, as otherwise some of the subdiagonal entries of $B_{ii}$ would have converged to [math]. In fact, the valuations of the eigenvalues of $A$ can be computed this way even if some of the eigenvalues are not defined over $\mathbb{Q}_{p}$.

2.8. Complexity analysis

We denote by $M(n)$ the cost of matrix multiplication.

Proposition 2.23**.**

Let $N$ be the precision and let $m$ be the largest size of a Jordan block of $A\pmod{p}$. Then the number of $\mathbb{Q}_{p}$-arithmetic operations performed by Algorithm 2.14 is as most

[TABLE]

Proof.

We tabulate the costs in a table.

[TABLE]

∎

The $\rho m^{3}$ term corresponds to blocks which cannot be decomposed using power iteration, i.e, where the classical algorithm is invoked. Note that the cost of computing $\chi_{A}\pmod{p}$ is negligible, since over a finite field the characteristic polynomial can be computed in $O(M(n)(\log n)(\log\log n))$ finite field operations [Storjohann2001frobenius].

If the linear factors occur with small multiplicities mod $p$, our algorithm significantly outperforms the naive strategy. If all linear factors mod $p$ occur with multiplicity $1$ the classical algorithm is never called, so the asymptotic run-time simplifies to ${O(n^{3}+\ell M(n)\log_{2}(N))}$. The worst-case scenario of the algorithm is if $\chi_{A}\pmod{p}\equiv(x-\alpha)^{n}$, in which case the default algorithm is called immediately; the only wasted time in this case is the inexpensive computation of $\chi_{A}\pmod{p}$.

We now study the complexity of the block Schur form algorithm.

Proposition 2.24**.**

Let $N$ be the precision, let $m$ be the largest size of a Jordan block of $A\pmod{p}$, and let $\ell$ be the number of linear factors of $\chi_{A}\pmod{p}$. Then the block Schur form can be computed in $O(n^{3}+\ell mNn^{2}+\rho m^{3})$ $\mathbb{Q}_{p}$-arithmetic operations.

Proof.

We tabulate the costs in a table.

[TABLE]

∎

2.9. $p$-adic Householder reflections

We close the section with a quick observation about Householder reflections.

Lemma 2.25**.**

Let $p$ be an odd prime. Let $x\in\mathbb{Z}_{p}^{n}$ be a vector with exactly one coordinate with minimal valuation $r$. We further assume this coordinate is not $x_{1}$. Let $e_{1}$ be the first standard basis vector. Then:

1. (a)

$x^{T}x$* is a non-zero square in $\mathbb{Z}_{p}$.* 2. (b)

Let $\alpha:=\sqrt{x^{T}x}\in\mathbb{Z}_{p}$ be one of the square roots and let $v^{\prime}:=x-\alpha e_{1}$. Choose (the unique) $v\in\mathbb{Z}_{p}^{n}$ such that $v^{\prime}=p^{r}v$. Then $\left\lvert v^{T}v\right\rvert_{p}=1$, and the Householder transformation

[TABLE]

is a well-defined reflection. In particular $H\in\operatorname{O}_{n}(\mathbb{Z}_{p})$. 3. (c)

$Hx=\alpha e_{1}$.

Proof.

1. (a)

Assume $x_{i}$ is the coordinate with minimal valuation. Write $x=p^{r}y$. We have that $x^{T}x=p^{2r}(y^{T}y)=p^{2r}(y_{i}^{2}+O(p))$. As $y_{i}$ is a unit, we see by Hensel’s lemma that this is a square. 2. (b)

Assume again $x_{i}$ is the coordinate with minimal valuation. Note that $\alpha^{2}=\sum_{j=1}^{n}x_{j}^{2}$, so as $x_{i}$ is the unique coordinate with minimal valuation, we have $\alpha=x_{i}+O(p^{r+1})$. Now

[TABLE]

By our assumption on $x$, we have $x_{i}-x_{1}=x_{i}+O(p^{r+1})$. Thus, $v^{T}v$ is a unit and $H$ is well-defined. That it is a reflection follows from the usual calculation, which is

[TABLE]

Finally, it is clear that all of the entries of $H$ are integral. 3. (c)

The result is clear from direct calculation.

∎

We now offer a Householder version of the Hessenberg algorithm. Our version of the Hessenberg algorithm is more a theoretical curiosity; for our practical implementations, we use [CRV2017characteristic, Algorithm 1] as it appears to be more efficient.

Algorithm 2.26

HessenbergForm($A$)

1:  for  $j=1,\ldots,(n-1)$ do

2:     if  The $j$-th subdiagonal is 0 then

3:        continue

4:     Swap the last (i.e, $n$-th) row with row $i$ to ensure $a_{nj}$ has the minimal valuation of the sub-diagonal elements

5:     Add multiples of the last ($n$-th row) to the $i$-th row, for $j<i<n$, to ensure $a_{nj}$ is the unique sub-diagonal element with minimal valuation

6:     Apply the corresponding column operations to preserve similarity. None of these operations affect columns $1,\ldots,j$

7:     Apply the Householder reflection to set the $j$-th subcolumn to a multiple of $e_{1}$.

3. Polynomial system solving

In this section, we solve 0-dimensional systems of polynomial equations over $\mathbb{Q}_{p}$. The general method is based on the truncated normal form solver of [TMV2018truncated]. We begin by reviewing some of their terminology. Throughout, we let $R:=\mathbb{Q}_{p}[x_{1},\ldots,x_{n}]$ and let $I$ be an ideal such that $d=\dim_{\mathbb{Q}_{p}}R/I<\infty$.

Definition 3.1**.**

A normal form on $R$ with respect to $I$ is a linear map $\mathcal{N}:R\rightarrow B$ where $B\subseteq R$ is a $\mathbb{Q}_{p}$-vector subspace of dimension $d$ such that the sequence

[TABLE]

is exact, and $\mathcal{N}|_{B}=\mathrm{id}_{B}$.

Consider the canonical exact sequence

[TABLE]

Let $s\colon R/I\rightarrow R$ be a section of $\mathbb{Q}_{p}$-vector spaces and let $B$ be the image of $R/I$. We can construct a normal form by setting $\mathcal{N}:=\pi\circ s$, giving the commutative diagram

[TABLE]

Vice-versa, any normal form $\mathcal{N}$ defines a section of $\pi$.

Following [TMV2018truncated], there is no need to compute an entire normal form for our purposes, and we can often restrict to a finite dimensional subspace.

Definition 3.2**.**

Let $B\subseteq V\subseteq R$ with $B,V$ finite dimensional $\mathbb{Q}_{p}$-vector subspaces such that $x_{i}B\subseteq V$ for all $i=1,\ldots,n$, and $\dim_{\mathbb{Q}_{p}}B=\dim_{\mathbb{Q}_{p}}(R/I)$. A Truncated Normal Form (TNF) on $V$ with respect to $I$ is a linear map $\mathcal{N}:V\rightarrow B$ such that $N$ is the restriction to $V$ of a normal form with respect to $I$. That is, the sequence

[TABLE]

is exact, and $\mathcal{N}|_{B}=\mathrm{id}_{B}$.

Similar to normal forms, from a commutative diagram of $\mathbb{Q}_{p}$-vector spaces

[TABLE]

we obtain a truncated normal form by $\mathcal{N}=s\circ\pi$, where $s$ is a section of $\pi$ and $x_{i}B\subseteq V$ for all $i$.

One way to construct a truncated normal form is via resultant matrices.

Definition 3.3**.**

Let $k$ be a field and let $f_{1},\ldots,f_{r}\in k[x_{1},\ldots,x_{n}]$ be polynomials. Let $D$ be a positive integer and let $d_{i}:=\deg f_{i}$. A resultant matrix is a $k$-linear map of the form

[TABLE]

where $V_{d}$ is the $k$-vector space of polynomials of degree at most $d$. (In particular, $V_{-n}=\{0\}$ for $n$ a positive integer.) We denote the corresponding matrix as $[\mathrm{Res}_{ij}]$.

Note that the image of the resultant map is contained in $V_{D}\cap I$. If $D$ is sufficiently large we obtain a truncated normal form. (The result below is quite old, and originally due to Macaulay.)

Proposition 3.4**.**

If $D=\sum_{i=1}^{r}(d_{i}-1)+1$, then $\operatorname{Im}\mathrm{Res}=I\cap V_{D}$. In particular, the canonical quotient $\pi\colon V_{D}\rightarrow V_{D}/(V_{D}\cap I)$ is represented by the matrix $(\ker[\mathrm{Res}_{ij}]^{T})^{T}$.

Proof.

See [TMV2018truncated, Section 4]. ∎

3.1. The algorithm

We now summarize our modifications of [TMV2018truncated, Algorithm 1]. Our implementation is available at:

\url

https://github.com/a-kulkarn/pAdicSolver

The critical theoretical underpinning of [TMV2018truncated, Algorithm 1] is Stickelberger’s Theorem [coxlittleoshea2015ideals, page 621].

Theorem 3.5**.**

Let $k$ be a field, let $R:=k[x_{1},\ldots,x_{n}]$, and let $I$ be an ideal such that $\dim_{k}R/I$ is finite. Let $p_{1},\ldots,p_{\ell}$ be the geometric points of the affine scheme $\operatorname{Spec}(R/I)$ with multiplicities $m_{1},\ldots,m_{\ell}$ respectively. Let $[f]\colon R/I\rightarrow R/I$ be the endomorphism of $R/I$ induced by multiplication by $f\in R$. Then the eigenvalues of $[f]$ and their multiplicities are $(f(p_{1}),m_{1}),\ldots,(f(p_{\ell}),m_{\ell})$. The $i$-th invariant subspace of $[f]$ is $\{g\in R/I:g(p_{j})=0\text{ for }j\neq i\}$.

Proof.

Let $k^{\textrm{al}}$ be an algebraic closure for $k$. We may assume that $R/I\otimes_{k}k^{\textrm{al}}$ is local by the Chinese Remainder Theorem. Note $R/I\otimes_{k}k^{\textrm{al}}$ is Artinian, so for $\mathfrak{m}$ the maximal ideal we have $\mathfrak{m}(\mathfrak{m}^{r})=\mathfrak{m}^{r}$ for some $r>0$. By Nakayama’s lemma we have $\mathfrak{m}^{r}=0$. Thus, for any $f\in\mathfrak{m}$, $[f]$ is nilpotent and so has [math] as an eigenvalue with multiplicity $\dim_{k^{\textrm{al}}}(R/I\otimes_{k}k^{\mathrm{al}})$. Since $R/I\otimes_{k}k^{\textrm{al}}=k^{\textrm{al}}\oplus\mathfrak{m}$ as a $k^{\textrm{al}}$-vector space, we now need only check the result for constant functions. For $f=1$, the result is clear. ∎

Note that if $k$ is not algebraically closed, the eigenvalues and invariant subspaces may only be defined over a finite extension of $k$; these correspond to points whose coordinates lie in a finite extension of $k$. Critically, the eigenvalues of $[x_{i}]$ are the $i$-th coordinates of the geometric points. We now give the $p$-adic solver algorithm:

Algorithm 3.6

SolveQpSystem($f_{1},\ldots,f_{m}$)

0:  Polynomials $f_{1},\ldots,f_{m}\in\mathbb{Q}_{p}[x_{1},\ldots,x_{n}]$ defining a [math]-dimensional polynomial system.

0:  A set of approximate solutions to the distinct $\mathbb{Q}_{p}$-solutions.

1:  Construct the resultant map $[\mathrm{Res}_{ij}]$ with $D$ sufficiently large

2:   Compute $[\pi_{ij}]=(\ker[\mathrm{Res}_{ij}]^{T})^{T}$

3:  Extract the submatrix $M$ of $[\pi_{ij}]$ whose columns correspond to monomials $m$ such that $\deg{m}<D$.

4:   Compute a square subblock of $M$ that is pnumerically well conditioned, using the $p$-adic $QR$-factorization. Let $\mathbf{b}$ be the basis corresponding to this sub-block.

5:  Construct the multiplication by $x_{i}$ matrices $[x_{1}]_{\mathbf{b}},\ldots,[x_{n}]_{\mathbf{b}}$ from the columns of $[\pi_{ij}]_{\mathbf{b}}$.

6:   Compute the eigenvectors of a random linear combination via iteration.

7:  return  The eigenvalues for $[x_{1}]_{\mathbf{b}},\ldots,[x_{n}]_{\mathbf{b}}$ indexed by the invariant subspaces.

By $[x_{i}]_{\mathbf{b}}$, we mean the matrix corresponding to the multiplication-by-$x_{i}$ operator in $\mathbf{b}$-coordinates. The choice of $B:=\mathrm{Span}(\mathbf{b})$ in step 4 defines a subspace of $V_{D}$ of the same dimension as $R/I$, hence a section $s\colon R/I\rightarrow V_{D}$ and thus a truncated normal form. Since $[x_{i}]B\subseteq V_{D}$, we can easily read off the action of the operator $[x_{i}]\colon B\rightarrow B$ from the columns of $[\pi_{ij}]$. For explicit details see [TMV2018truncated, Theorem 4 and following text].

The differences between the $\mathbb{Q}_{p}$ and $\mathbb{R}$ versions of the algorithm are precisely the adaptations in the linear algebra, i.e, steps 2, 4, and 6. Note that if $k$ is any finite precision field for which these steps can be implemented, the general recipe of [TMV2018truncated, Algorithm 1] allows us to compute the approximate solutions. Examples include ${k=\mathbb{R},\mathbb{C},\mathbb{Q}_{p},\mathbb{F}_{q}(\!(t)\!)}$.

Finally, instead of computing a truncated normal form in steps $1$ and $2$, we may simply use a Gröbner basis instead. For systems where a Gröbner basis is already known or where the polynomials in the basis have small coefficients, the numerous optimizations of Gröbner basis algorithms are likely preferable from a practical standpoint. The major advantage of our implementation is the fact that the coefficient size is capped at $p^{N}$.

4. Acknowledgements

I would like to thank Yue Ren, David Roe, and Simon Telen for helpful discussions about the article. I would also like to thank Marta Panizzut, Emre Sertöz, and Bernd Sturmfels for their helpful comments.