Deciding Subspace reachability problems with application to Skolem’s Problem

Traditionally, work toward determining the decidability of the Orbit Problem and related problems such as Skolem’s Problem has leveraged heavy machinery from various domains of algebraic geometry, mathematical logic, and number theory — particularly algebraic and transcendental number theory. In contrast, our approach is to develop a machinery suitable for application to the Orbit Problem by leveraging basic geometric and dynamical insight.

In essence, we consider not the dynamical system $(\mathbb{R}^{d},A)$ the Orbit Problem is formulated in, but rather the associated induced dynamical system over real projective space (i.e. after projecting onto $S^{d-1}$ and identifying opposite points). Working instead in $\mathbb{R}\mathbb{P}^{d-1}$, we can ask an equivalent version of the Orbit Problem that is decidable if and only if the original formulation is decidable. We find that by compactifying the state space in such a way, the asymptotic behavior of the system becomes far more understandable: in projective space, we have asymptotically stable (converging) orbits, while no-such behavior is expressed in Euclidean space. As a consequence, we can use arguments unavailable in the original formulation. The central contribution of our paper can then be interpreted as a fine-grained examination of the $\omega$-limit sets — or more generally the stable manifolds — of dynamical systems $(\mathbb{R}\mathbb{P}^{d-1},f)$ for $f\in PGL(\mathbb{R}^{d})$, namely linear dynamical systems over real projective space, with this deepened understanding lending itself to deciding the higher-dimensional Orbit Problem.

This reduction is both surprising and noteworthy for a number of reasons. First, we note that although we lose an entire manifold dimension after reducing the problem to projective space, we maintain all the necessary information for deciding the Orbit Problem — indeed, the reduction can be seen as keeping only the “essential” aspects of the Orbit Problem. In addition, this reduction compactifies the state space, and as a consequence enables the application of general tools and intuition from the study of continuous dynamical systems over a compact space, for which there is far more structure and results than in the unbounded case, enabling the application of an entire area of mathematics not used to date.

To the best of our knowledge, no previous work on the Orbit Problem and related questions leverage the techniques introduced here. However, a few works come close in spirit to the asymptotic analysis used here, particularly work on termination of linear loops (see [akshay2022robustness, ouaknine2014positivity, braverman2006termination, hosseiniIntegerTermination, tiwari2004termination]), where the observation that the largest eigenvalues dominate the behavior of the program is used repeatedly. A similar technique is well-known in the study of linear recurrence sequences, whereby restricting the spectral structure of the problem and using the closed-form solution of an LRS, it can be clear how the characteristic roots of maximal modulus dominate the asymptotic behavior. Although this paper leverages the same basic principle, we take it much farther to obtain finer results.

One of the original sources of motivation for studying the Orbit Problem was Skolem’s Problem, due to the fact that Skolem’s Problem reduces to the Orbit Problem when the target space has dimension $d-1$. Although Skolem’s Problem is easy to describe, it is difficulty to prove deep theorems on account of a lack of machinery. Indeed, in terms of lower bounds, it is known that Skolem’s Problem is NP-hard [blondel2002presence] — which translates to the Orbit Problem when restriction is not placed on the dimension of the target subspace. Moreover, deep work has shown that further advances on the Skolem and Orbit Problem for higher dimensions are likely to be hard due to the fact that such advances in related problems would entail major advances in Diophantine approximation [ouaknine2014positivity]. Despite these difficulties, there continues to be notable advances toward Skolem’s problem (see e.g. [lipton2022skolem, biluMFCSskolem]), including recent promising new techniques by way of Universal Skolem Sets [kenison2020skolem, luca2021universal, lucaMFCSuniversal, luca2023skolem]. Nonetheless, despite the intense study on this problem, our understanding remains incomplete.

Another principle source of motivation in studying the Orbit Problem and related questions comes from the problem of program verification. Specifically, enormous effort has been dedicated to solving the “Termination Problem,” chiefly concerned with deciding whether a “linear while loop” will terminate (see e.g. [tiwari2004termination, braverman2006termination, hosseiniIntegerTermination, cook2006termination, cook2006terminator, bozga2012deciding, ben2013linear, bradley2005termination, ouaknine2014termination, karimov2022s, vyalyi2011orbits, luca2022algebraic]). In particular, we note that the Polytope Hitting Problem where the target is an intersection of half spaces, is closely connected to the higher-dimensional Orbit Problem, and more immediately translates to problems of program verification and termination of linear loops [chonev2014polyhedron].

The author is grateful to David Cash for the continued feedback and support.

Let $\mathbb{A}$ denote the field of algebraic numbers; a complex number $\alpha$ is algebraic if it is a root of a single variable polynomial with integer coefficients. The matrix $A$ is taken to have rational entries, so all eigenvalues $\lambda$ of $A$ are in $\mathbb{A}$. This paper requires we effectively perform various operations with the eigenvalues of $A$, which may not lie in $\mathbb{Q}$. Fortunately, there is a sizable literature concerning computation with algebraic numbers (see [AGAlgorithms, cohen2013course], or the appendix of [chonev2016complexity] for a useful review). Below, we summarize the basic properties of the effective manipulation of algebraic numbers pertaining to this paper.

For $\alpha\in\mathbb{A}$, the defining polynomial of $\alpha$, denoted $p_{\alpha}$, is the unique polynomial of least degree vanishing at $\alpha$, where the coefficients do not have common factors. Define the degree and height $H(p)$ of $\alpha$ to be the degree of $p$, and the maximum of the absolute values of $p$’s coefficients, respectively. As such, standard finite encoding of an algebraic number $\alpha$ is as a tuple composed of its defining polynomial, and rational approximations of its real and imaginary parts of precision sufficient to distinguish $\alpha$ from the other roots of $p_{\alpha}$. That is, $\alpha\in\mathbb{A}$ takes the representation $(p_{\alpha},a,b,r)\in\mathbb{Z}[x]\times\mathbb{Q}^{3}$, when $\alpha$ is the unique root of $p_{\alpha}$ inside a circle in $\mathbb{C}$ of radius $r$, centered at $a+bi$. This representation of $\alpha$ is well-defined and computable in part thanks to a useful separation bound of Mignotte [mignotte1982some], which asserts that for distinct roots $\alpha,\beta$ of polynomial $p\in\mathbb{Z}[x]$

where $d=deg(p)$ and $H=H(p)$. As a consequence, when $r$ is less than the root separation bound, we have equality checking between algebraic numbers. Given distinct $\alpha,\beta\in\mathbb{A}$ these are roots of $p_{\alpha}p_{\beta}$, which is of degree at most $deg(\alpha)+deg(\beta)$, and of height at most $H(\alpha)H(\beta)$. Then one may compute $\alpha+\beta$, $\alpha\beta$, $1/\alpha$, $\overline{\alpha}$, $|\alpha|$, decide whether $\alpha>\beta$ for distinct algebraic $\alpha$ and $\beta$, and more.

For the remainder of this paper, we do not make explicit use of these methods, implicitly assuming they are used in the following procedures for deciding instances of the higher-dimensional Orbit Problem. Indeed, it should be clear from the algorithms given in the sequel that the we never leave the field of algebraic numbers. Although, the techniques of this paper do not in fact require the precision the effective manipulation of algebraic numbers provides: all our techniques are quite amenable to the finite-precision setting. So long as precision can be increased as needed, for our purposes a finite approximation of a number paired with an error bound more than suffices. Nonetheless, for the sake of clarity we opt to compute with algebraic numbers as discussed above.

As a consequence of these observations, we see that we have effective techniques for computing the eigenvalues and eigenvectors of rational matrices $A$.

We now consider basic aspects of linear recurrence sequences. For a rich treatment see [everest2003recurrence]. A linear recurrence sequence (LRS) over $\mathbb{Q}$ is an infinite sequence $\{u_{n}\}_{n=0}^{\infty}$ of terms in $\mathbb{Q}$ satisfying the recurrence relation

where $a_{0},...,a_{d-1}\in\mathbb{Q}$, with $a_{0}\neq 0$ and $u_{j}\neq 0$ for at least one $j$ in the range $0\leq j\leq m-1$. We say such an LRS has order $d$. We call the $a_{0},...,a_{d-1}$ the coefficients of the sequence $\{u_{n}\}$ and the initial terms of $\{u_{n}\}$ are $u_{0},...,u_{d-1}$. The characteristic polynomial of the sequence $\{u_{n}\}$ is

with the $\lambda_{1},...,\lambda_{k}$ called the characteristic roots of the sequence $\{u_{n}\}$.

We call the sequence simple if $k=d$, so $m_{1}=\cdots=m_{k}=1$, and non-degenerate if $\lambda_{i}/\lambda_{j}$ is not a root of unity for any distinct $i,j$. The study of arbitrary LRS can be reduced effectively to that of non-degenerate LRS by partitioning the original LRS into finitely many non-degenerate subsequences.

Given an LRS, define

to be the companion matrix of the LRS, and take the vector $x\in\mathbb{Q}^{d}$ as $(u_{d-1},...,u_{0})^{T}$ of initial terms of the LRS. Then iteration of $A$ over $x$ acts as a “shift” on the entries of $x$, shifting in the next term of the LRS and dropping the oldest term. In addition, the characteristic polynomial of the LRS is the characteristic polynomial of $A$, and the characteristic roots of the LRS are the eigenvalues of $A$.

The algorithm presented in this paper requires the computation of angles between Euclidean subspaces. The subject of computing angles between Euclidean subspaces — also known as computing angles between flats — is a well-studied area and there are many effective and efficient techniques for doing so. See [BFb0062107, gunawan2005formula, bjorck1973numerical, jiang1996angles] for readable introductions to the subject, and a variety of effective techniques for computing angles between subspaces.

Let $U,W\subset\mathbb{R}^{d}$ be $k$ and $l$ dimensional subspaces, respectively, presented by orthogonal rational basis $\{u_{1},...,u_{k}\}$ and $\{w_{1},...,w_{l}\}$. Note any linearly independent set of rational vectors can be effectively transformed into an orthogonal set of rational vectors via the Gram-Schmidt Process (before normalization). If instead it is desired that the basis be orthonormal, then the entries of the vectors will be algebraic after the renormalization step.

In this paper we are primarily concerned with determining the minimal angle $0\leq\theta_{1}\leq\pi/2$ between two subspaces, defined as the least angle between any pair of non-zero vectors from $U$ and $W$. This is the first principle angle between subspaces, whose variational characterization is

where $\langle\cdot,\cdot\rangle$ is the standard inner (dot) product. We refer to the quantity expressed above in Equation 1 as the minimal angle between subspaces $U$ and $W$. We always assume $U$ and $W$ have trivial intersection $U\cap W=\{0\}$, so that the minimal angle is well-defined and non-zero.

It is well-known that the variational characterization of singular values implies a variational characterization of the angles between subspaces [golub2013matrix]. Indeed, as a consequence of the variational characterization of the singular value decomposition of a $m\times n$ matrix $A$, we see that for $x\in\mathbb{C}^{m}$, $y\in\mathbb{C}^{n}$, then

where $\sigma_{1}$ is the first singular value of $A$. Dropping the $\arccos$ term in Equation 1, we obtain a new problem, namely maximizing the quantity

Note then that $\sigma\leq 1$. The pair of vectors $u,v$ maximizing the above quantity can be used to obtain the minimal angle between the subspaces $U$ and $W$.

Combining the above observations, take matrices $X\in\mathbb{R}^{n\times k}$ and $Y\in\mathbb{R}^{n\times l}$ to have columns consisting of orthonormal bases $\{u_{1},...,u_{k}\}$ and $\{w_{1},...,w_{l}\}$ of $U$ and $W$, respectively. The optimization problem in Equation 2 can then be written as

the solution to which is the largest singular value of $X^{\top}Y$, by the variational characterization of the singular value decomposition. For more details and proof of the fact that cosines of principle angles come from the singular value decomposition see [knyazev2002principal], or refer to Chapter 5 of [meyer2023matrix] for more details on the variational characterization of singular value decomposition, and additional effective methods for computing the minimal angle between subspaces.

The value $\theta=\arccos\sigma$ is the minimal angle between subspaces $U,W$. However, we need not learn $\theta$ — in fact $\theta$ need not be an algebraic number due to the $\arccos$. Rather, it is sufficient for the purposes of this paper to instead work only with the cosine of the minimal angle, namely $\sigma$. The algorithm presented in this paper only requires that we compare the minimal angle between subspaces, which we can accomplish equipped just with $\sigma$: supposing $\sigma^{\prime}<\sigma$, we know $\arccos(\sigma^{\prime})>\arccos(\sigma)$, and hence the angle corresponding to $\sigma$ is smaller than the angle corresponding to $\sigma^{\prime}$. In fact it is sufficient not to work with $\sigma$, but $\sigma^{2}$, which is the largest of the positive eigenvalues obtained in the singular value decomposition.

The singular value $\sigma$ is merely the square root of an eigenvalue computed in the singular value decomposition (or just the eigenvalue in the $\sigma^{2}$ case), and is hence algebraic, admitting a finite representation as per Section 2.1. As a consequence, we may effectively compare the minimal angles between different pairs of Euclidean subspaces described by rational basis vectors.

For the remainder of this paper, when we speak of angles between subspaces, we are implicitly speaking of the singular value $\sigma$ expressed above. Define the function

to be the function computing the minimal angle between subspaces $U$ and $W$ as above. Whether $\Gamma(U,W)$ returns $\sigma$ or $\sigma^{2}$ does not matter for the purposes of this paper — intuitively $\Gamma(U,W)$ computes the minimal angle between the subspaces ($\arccos\Gamma$ is the actual minimal cosine angle), and everything remains computable: the columns of $U$ and $W$ will always be composed of rational (or possibly algebraic) numbers, and $\Gamma$ is effectively computable with such input.

Real projective space $\mathbb{R}\mathbb{P}^{d-1}$ is defined to be the quotient of $\mathbb{R}^{d}\setminus\{0\}$ by the equivalence relation $x\sim\alpha x$ where $\alpha\neq 0$ is real and $x\in\mathbb{R}^{d}$. Clearly, it suffices to consider only vectors $x$ with Euclidean norm $||x||=1$. Hence, geometrically, projective space is the space of lines through the origin, or alternatively the space obtained by identifying opposite points on the unit sphere $S^{d-1}$. The projective linear group $PGL(\mathbb{R}^{d})=GL(\mathbb{R}^{d})/\alpha\cdot\text{Id}$, $\alpha\neq 0$, is the induced action of the general linear group on projective space. The induced action of non-invertible square matrices over elements of projective space is similarly defined.

Write $\mathbb{P}:\mathbb{R}^{d}\setminus\{0\}\rightarrow\mathbb{R}\mathbb{P}^{d-1}$ for the projection of elements from Euclidean space to real projective space, denoting elements of $\mathbb{R}\mathbb{P}^{d-1}$ by $p=\mathbb{P}x$, where $0\neq x\in\mathbb{R}^{d}$. A metric on $\mathbb{R}\mathbb{P}^{d-1}$ is given by

In the case that $W$ is a subspace of $\mathbb{R}^{d}$, define

Observe that these metrics are intimately related to the function $\Gamma$ defined in the previous subsection, as expressed in the following trivial, albiet useful lemma.

Given $x\in\mathbb{R}^{d}$, let $X$ denote the one-dimensional subspace (line) spanned by $x$. Then, for a matrix $A\in\mathbb{R}^{d\times d}$ and subspace $U$, we see that there exists $n\in\mathbb{N}$ such that $A^{n}x\in U$ if and only if $A^{n}X\subseteq U$, where $A^{n}X$ is the subspace spanned by the vector $A^{n}x$. Hence, given an instance $(A,x,U)$ of the Orbit Problem, it is immediate that there is an $n\in\mathbb{N}$ such that $A^{n}x\in U$ if and only if $\mathbb{P}A^{n}x\in\mathbb{P}U$.

These considerations lead us to the simple yet powerful conclusion that by projecting the Orbit Problem into projective space we still retain the needed information to decide the Orbit Problem — indeed the formulations are equivalent — and yet in projective space the asymptotic behavior of orbits is far more manageable than in Euclidean space. In projective space, we have asymptotically stable (converging) orbits, while no-such behavior is expressed in Euclidean space. We abuse this structure in the sequel, applying the basic insight that if an orbit of a dynamical system monotonically converges to some set, and the attracting set is separated from the target set by some $\epsilon>0$, then by running the system for finite time we can decide whether the orbit intersects the target set. See [colonius2014dynamical] for details on linear dynamical systems over real projective space.

We begin by stating the following well-known result about Jordan blocks, easily proved by induction.

Note that that when $\lambda=\alpha\pm i\beta$ is complex $|\lambda|=\sqrt{\alpha^{2}+\beta^{2}}$, hence with $\theta\in[0,2\pi)$ determined by $\cos\theta=\alpha/\sqrt{\alpha^{2}+\beta^{2}}$, we can rewrite the matrix $R$ above as

Hence $R$ describes a rotation by the angle $\theta$ followed by multiplication by $|\lambda|$. Letting $H$ denote the $2D_{i}$ dimensional nilpotent matrix of ones just above the principle diagonal, for Jordan block $J_{i}$ associated with a complex conjugate eigenvalue pair, we have

where $\tilde{R}$ is the block diagonal matrix with blocks $R$, and $n\geq D_{i}-2$.

The following lemmas establish the asymptotic behavior of a non-zero vector under iteration of Jordan blocks. Tiwari and Braverman in [tiwari2004termination, braverman2006termination] touch upon a similar idea, although the outcome and use is different.

We also note that the following Lemmas 5, 6, and 7 are essentially folklore in dynamical systems; the asymptotic behavior of real projective transformations is well understood — see e.g. [colonius2014dynamical, ayala2014topological, KUIPER197613, he2017topological] which all contain relevant background and results related to the Lemmas below. However, we write the following Lemmas in a way that corresponds their computable nature, thereby translating them into a form suitable for the problem this paper is attacking.

With Lemma 5 in hand, we obtain the following set of lemmas giving asymptotic behavior of orbits when there are multiple Jordan blocks.

The following lemma states that when there are many Jordan blocks all associated with eigenvectors of the same modulus, the largest block(s) dictate the asymptotic behavior of orbits.

The following lemma establishes the case when the Jordan blocks composing a block-diagonal Jordan matrix $J$ correspond to eigenvalues of different magnitude.

With these lemmas in hand, we move to prove the theorems stated in the introduction of this paper.