Resolvability of Hamming Graphs

TL;DR

This paper characterizes the resolvability of Hamming graphs and hypercubes, introduces efficient methods to determine resolving sets, and demonstrates practical applications in protein sequence analysis using low-dimensional vector representations.

Contribution

It provides a new linear system characterization of Hamming graph resolvability and develops ILP and Gr"obner basis methods for assessing resolving sets.

Findings

Characterized resolvability of Hamming graphs via linear systems.

Developed ILP and Gr"obner basis methods for resolving set determination.

Identified a resolving set of size 77 in octapeptide sequences.

Abstract

A subset of vertices in a graph is called resolving when the geodesic distances to those vertices uniquely distinguish every vertex in the graph. Here, we characterize the resolvability of Hamming graphs in terms of a constrained linear system and deduce a novel but straightforward characterization of resolvability for hypercubes. We propose an integer linear programming method to assess resolvability rapidly, and provide a more costly but definite method based on Gr\"obner bases to determine whether or not a set of vertices resolves an arbitrary Hamming graph. As proof of concept, we identify a resolving set of size 77 in the metric space of all octapeptides (i.e., proteins composed of eight amino acids) with respect to the Hamming distance; in particular, any octamer may be readily represented as a 77-dimensional real-vector. Representing k-mers as low-dimensional numerical vectors may enable new applications of machine learning algorithms to symbolic sequences.

Tables

Table 1. Table 1: Time in seconds required to determine resolvability for each technique. Fifty resolving and fifty non-resolving sets, selected uniformly at random, were considered for each Hamming graph ℍ k , a subscript ℍ 𝑘 𝑎 \mathbb{H}_{k,a} . Means and standard deviations consider five replicates per set.

		Brute Force		Gröbner Basis		ILP
($k,a$)	$a^k$	Mean	SD	Mean	SD	Mean	SD
(2,2)	$4$	$3.88\text{×}{10}^{-05}$	$1.51\text{×}{10}^{-06}$	$6.79\text{×}{10}^{-03}$	$1.06\text{×}{10}^{-03}$	$1.28\text{×}{10}^{-01}$	$3.53\text{×}{10}^{-03}$
(2,4)	$16$	$2.47\text{×}{10}^{-04}$	$6.83\text{×}{10}^{-05}$	$2.25\text{×}{10}^{-02}$	$2.59\text{×}{10}^{-03}$	$1.16\text{×}{10}^{-01}$	$4.84\text{×}{10}^{-03}$
(3,3)	$27$	$5.02\text{×}{10}^{-4}$	$2.45\text{×}{10}^{-4}$	$2.83\text{×}{10}^{-2}$	$7.92\text{×}{10}^{-3}$	$1.21\text{×}{10}^{-01}$	$8.12\text{×}{10}^{-3}$
(5,2)	$32$	$6.61\text{×}{10}^{-04}$	$3.29\text{×}{10}^{-04}$	$3.14\text{×}{10}^{-02}$	$5.27\text{×}{10}^{-03}$	$1.28\text{×}{10}^{-01}$	$3.91\text{×}{10}^{-03}$
(3,5)	$125$	$8.98\text{×}{10}^{-03}$	$5.38\text{×}{10}^{-03}$	$1.12\text{×}{10}^{-01}$	$2.91\text{×}{10}^{-02}$	$1.37\text{×}{10}^{-01}$	$7.02\text{×}{10}^{-03}$
(5,3)	$243$	$2.78\text{×}{10}^{-2}$	$1.96\text{×}{10}^{-2}$	$1.22\text{×}{10}^{-1}$	$7.88\text{×}{10}^{-2}$	$1.20\text{×}{10}^{-1}$	$8.12\text{×}{10}^{-3}$
(8,2)	$256$	$2.85\text{×}{10}^{-02}$	$2.21\text{×}{10}^{-02}$	$9.87\text{×}{10}^{-02}$	$1.96\text{×}{10}^{-02}$	$1.17\text{×}{10}^{-01}$	$1.59\text{×}{10}^{-03}$
(4,4)	$256$	$3.13\text{×}{10}^{-02}$	$1.97\text{×}{10}^{-02}$	$1.27\text{×}{10}^{-01}$	$3.90\text{×}{10}^{-02}$	$1.37\text{×}{10}^{-01}$	$9.58\text{×}{10}^{-03}$
(5,5)	$3125$	$5.19\text{×}{10}^{00}$	$3.17\text{×}{10}^{00}$	$4.00\text{×}{10}^{00}$	$3.54\text{×}{10}^{00}$	$1.35\text{×}{10}^{-01}$	$1.09\text{×}{10}^{-02}$
(12,2)	$4096$	$6.28\text{×}{10}^{00}$	$5.34\text{×}{10}^{00}$	$2.93\text{×}{10}^{-01}$	$7.24\text{×}{10}^{-02}$	$1.24\text{×}{10}^{-01}$	$2.39\text{×}{10}^{-03}$
(6,4)	$4096$	$7.78\text{×}{10}^{00}$	$4.65\text{×}{10}^{00}$	$7.73\text{×}{10}^{-01}$	$3.67\text{×}{10}^{-01}$	$1.52\text{×}{10}^{-01}$	$8.99\text{×}{10}^{-03}$
(8,3)	$6561$	$2.02\text{×}{10}^{01}$	$1.40\text{×}{10}^{01}$	$1.12\text{×}{10}^{01}$	$1.47\text{×}{10}^{01}$	$1.62\text{×}{10}^{-01}$	$1.41\text{×}{10}^{-02}$

Equations

$$\underset{y}{\min }$$

$$\sum _{j\in \mathbb{V}}|d(u,j)-d(v,j)|\cdot y_j\ge 1,\forall u\ne v\in \mathbb{V}$$

$$y_j\in \{0,1\},\forall j\in \mathbb{V}.$$

$${\hbox{vec}}\left(\left[\begin{array}[]{cc}a&b\\ c&d\end{array}\right]\right)=(a,c,b,d).$$

$$A:=\left(\begin{array}[]{c}{\hbox{vec}}(V_{1})\\ \vdots\\ {\hbox{vec}}(V_{n})\end{array}\right).$$

$$A_0=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 1\\ 0& 1& 0& 0& 1& 0\end{array}\right],$$

$$rref(A_1)=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 1\\ 0& 1& 0& 0& 1& 0\\ 0& 0& 1& 0& 0& 1\end{array}\right].$$

$$B:=\left[\begin{array}[]{c}v_{1}-\bar{v_{1}}\\ \vdots\\ v_{n}-\bar{v_{n}}\end{array}\right].$$

$$A=\left[\begin{array}[]{ccccc}\vline&\vline&&\vline&\vline\\ A_{1}&A_{2}&\ldots&A_{2k-1}&A_{2k}\\ \vline&\vline&&\vline&\vline\end{array}\right]$$

$$B$$

$$y$$

$$\Big{(}\llbracket v[1]=1\rrbracket-\llbracket v[1]=0\rrbracket,\ldots,\llbracket v[k]=1\rrbracket-\llbracket v[k]=0\rrbracket\Big{)}=v-\bar{v},$$

$$C:=\left[\begin{array}[]{c}v_{1}\\ \vdots\\ v_{n}\end{array}\right].$$

$$Bz=0\hbox{ if and only if }Cz=\langle 1^{k},z\rangle\left[\begin{array}[]{c}1/2\\ \vdots\\ 1/2\end{array}\right].$$

$$P_1$$

$$P_2$$

$$P_3$$

$$P:=(P_1\cup P_2\cup P_3).$$

$$B_i$$

$$(z_1,\dots ,z_a)⟶(z_{(i-1)a+1},\dots ,z_{ia}).$$

$$G_1=\left\{\sum _{i=1}^a{z}_i\right\}\bigcup _{2\le i\le a}\{z_i^3-z_i\}\bigcup _{2\le i<j\le a}\{z_i^2{z}_j+z_i{z}_j^2\}\bigcup _{2\le i<j<\ell \le a}\{z_i{z}_j{z}_{\ell }\}.$$

$$f(z):=\Big{(}\sum_{j=1}^{ka}z_{j}^{2}\Big{)},$$

$$\left\{\begin{array}[]{rcl}z_{1}+z_{6}&=&0;\\ z_{2}+z_{5}&=&0.\end{array}\right.$$

$$\left\{\begin{array}[]{ccl}0&=&z_{1}^{3}-z_{1};\\ 0&=&z_{2}^{3}-z_{2};\\ 0&=&z_{3}^{3}-z_{3};\\ 0&=&z_{1}+z_{2}+z_{3};\\ 0&=&(2-z_{1}^{2}-z_{2}^{2}-z_{3}^{2})\cdot(z_{1}^{2}+z_{2}^{2}+z_{3}^{2});\\ \hline\cr 0&=&z_{4}^{3}-z_{4};\\ 0&=&z_{5}^{3}-z_{5};\\ 0&=&z_{6}^{3}-z_{6};\\ 0&=&z_{4}+z_{5}+z_{6};\\ 0&=&(2-z_{4}^{2}-z_{5}^{2}-z_{6}^{2})\cdot(z_{4}^{2}+z_{5}^{2}+z_{6}^{2});\end{array}\right.$$

$$f(z)=z_1^2+z_2^2+z_3^2+z_4^2+z_5^2+z_6^2.$$

$$R\text{ does text@underlinenot resolve }{\mathbb{H}}_{k,a}⟺\exists z\ne 0\text{such that}z\in ker(A)\cap \{P=0\}.$$

$$\mathcal{I}:=\bigcap _{i=1}^k\left\{z\in {\mathbb{Z}}^{ak}such\; that\sum _{j=(i-1)a+1}^{ia}{z}_j=0\; and\sum _{j=(i-1)a+1}^{ia}|z_j|\le 2\right\}.$$

$${\mathcal{J}}_1$$

$${\mathcal{J}}_2$$

$$\left\{z\in {\mathbb{R}}^{ak}such\; that\sum _{j=1}^a|z_j|\le 2\right\}=\bigcap _{w\in \{-1,1{\}}^{ak}}\left\{z\in {\mathbb{R}}^{ak}such\; that\sum _{j=1}^a{w}_j{z}_j\le 2\right\},$$

$${\mathcal{I}}_i:=\left\{z\in {\mathbb{Z}}^{ak}such\; that\sum _{j=(i-1)a+1}^{ia}{z}_j=0\; and\sum _{j=(i-1)a+1}^{ia}|z_j|\le 2\right\}.$$

$$\underset{z\in {\mathbb{R}}^{ak}}{\min }$$

$$Az=0\text{and}z\in \mathcal{I},$$

$$z\in {\mathbb{R}}^{ak}$$

Full text

{AMS}

05C12, 05C50, 05C62, 68R10, 90C35, 92C40

1 Introduction

In what follows, $k\geq 1$ and $a\geq 2$ are fixed integers, and we refer to elements in the set $\mathbb{V}:=\{0,\ldots,a-1\}^{k}$ as $k$-mers, which we represent either as strings or row vectors depending on the context.

The Hamming distance between two $k$-mers $u$ and $v$, from now on denoted as $d(u,v)$, is the number of coordinates where the $k$-mers differ, and is a valid metric. The Hamming graph $\mathbb{H}_{k,a}$ has $\mathbb{V}$ as its vertex set, and two $k$-mers $u$ and $v$ are adjacent (i.e. connected by an undirected edge) if and only if $d(u,v)=1$, i.e. $u$ and $v$ differ at exactly one coordinate. As a result, the (geodesic) distance between two vertices in $\mathbb{H}_{k,a}$ is precisely their Hamming distance (see Figure 1). The literature refers to the Hamming graph with $a=2$ as the ($k$-dimensional) hypercube.

A non-empty set $R\subseteq\mathbb{V}$ is called resolving when for all $u,v\in\mathbb{V}$, with $u\neq v$, there exists $r\in R$ such that $d(u,r)\neq d(v,r)$. In other words, $R$ multilaterates $\mathbb{V}$. For instance, $\mathbb{V}$ resolves $\mathbb{H}_{k,a}$ because $d(u,v)=0$ if and only if $u=v$. Equivalently, $R\subseteq\mathbb{V}$ is resolving if and only if the transformation $\Phi:\mathbb{V}\to\mathbb{R}^{|R|}$ defined as $\Phi(v):=(d(v,r))_{r\in R}$ is one-to-one. In particular, the smaller a resolving set of $\mathbb{H}_{k,a}$, the lower the dimension needed to represent $k$-mers as points in a Euclidean space, which may be handy e.g. to represent symbolic data numerically for machine learning tasks [29].

The metric dimension of $\mathbb{H}_{k,a}$, which we denote $\beta(\mathbb{H}_{k,a})$, is defined as the size of a minimal resolving set in this graph [17, 27]. For instance, $\beta(\mathbb{H}_{1,a})=(a-1)$ because $\mathbb{H}_{1,a}$ is isomorphic to $K_{a-1}$, the complete graph on $(a-1)$ vertices [5]. Unfortunately, computing the metric dimension of an arbitrary graph is a well-known NP-complete problem [6, 13, 20], and it remains unknown if this complexity persists when restricted to Hamming graphs. In fact, the metric dimension of hypercubes is only known up to dimension $k=10$ [17], and values have been conjectured only up to dimension $k=17$ [24]—see OEIS sequence A303735 for further details [25].

Integer linear programming (ILP) formulations have been used to search for minimal resolving sets [5, 9]. In the context of Hamming graphs, a potential resolving set $R$ is encoded by a binary vector $y$ of dimension $a^{k}$ such that $y_{j}=1$ if $j\in R$ and $y_{j}=0$ if $j\in\mathbb{V}\setminus R$. One can then search for a minimal resolving set for $\mathbb{H}_{k,a}$ by solving the ILP [5]:

[TABLE]

The first constraint ensures that for all pairs of different vertices $u$ and $v$, there is some $j\in R$ such that $|d(u,j)-d(v,j)|>0$, hence $R$ resolves $\mathbb{H}_{k,a}$. The objective penalizes the size of the resolving set. (A variant due to [9] is similar but stores $a^{k}$ copies of a binary version of the distance matrix of the graph.) One downside of this formulation is that forming the distance matrix of $\mathbb{H}_{k,a}$ requires $\mathcal{O}(a^{2k})$ storage, as well as significant computation. Moreover, standard approaches to reduce the computation below $\mathcal{O}(a^{2k})$, such as fast multipole methods [15] and kd-trees [2], do not obviously apply. Even if one could compute all pairwise distances between nodes, simply storing the distance matrix is impractical. To fix ideas, the graph $\mathbb{H}_{8,20}$—which is associated with octapeptides (see Section 6)—has $20^{8}$ nodes, so storing the distance matrix with $\log_{2}(8)=3$ bits per entry and taking advantage of symmetry would require $3{20^{8}\choose 2}$ bits, or approximately a prohibitive 123 exabytes.

Due to the above difficulties, other efforts have focused on finding small resolving sets rather than minimal ones. When $a^{k}$ is small, resolving sets for $\mathbb{H}_{k,a}$ may be determined using the so-called Information Content Heuristic (ICH) algorithm [18], or a variable neighborhood search algorithm [24]. Both approaches quickly become intractable with increasing $k$. However, the highly symmetric nature of Hamming graphs can be taken advantage of to overcome this problem. Indeed, recent work [29] has shown that $\beta(\mathbb{H}_{k,a})\leq\beta(\mathbb{H}_{k-1,a})+\lfloor a/2\rfloor$; in particular, $\beta(\mathbb{H}_{k,a})\leq(k-1)\lfloor a/2\rfloor+(a-1)$ i.e., just $\mathcal{O}(k)$ nodes are enough to resolve all the $a^{k}$ nodes in $\mathbb{H}_{k,a}$. Moreover, one can find a resolving set of size $\mathcal{O}(k)$ in only $\mathcal{O}(ak^{2})$ time [29].

This manuscript is based on the recent Bachelor’s thesis [21], and has two overarching goals. First, it aims to develop practical methods for certifying the resolvability, or lack thereof, of subsets of nodes in arbitrary Hamming graphs. So far, this has been addressed for hypercubes in the literature [1] but remains unexamined for arbitrary values of the parameter $a$. While our work does not directly address the problem of searching for minimal resolving sets, verifying resolvability is a key component of any such search and may shed new light on the precise metric dimension of $\mathbb{H}_{k,a}$ in future investigations. Second, this paper aims also to exploit said characterization to remove unnecessary nodes—if any—in known resolving sets. This problem, which is infeasible by brute force when $a^{k}$ is large, has not received any attention in the literature despite being crucial for the embedding of $k$-mers into the Euclidean space of a lowest possible dimension.

The paper is organized as follows. Our main theoretical results are presented first in Section 2. Theorem 2.3 provides the foundation from which we address the problem of verifying resolvability in Hamming graphs and implies a new characterization of resolvability of hypercubes (Corollary 2.7). An illustrative example shows the utility of Theorem 2.3 but raises several practical challenges in its implementation on large Hamming graphs. Section 3 describes a computationally demanding verification method based on Gröbner bases that is nevertheless more efficient than the brute force approach and determines with certainty whether or not a set of nodes in $\mathbb{H}_{k,a}$ is resolving. Computational issues are addressed in Section 4 with a novel ILP formulation of the problem. This approach is fast but stochastic and hence has the potential to produce false positives or false negatives. Section 5 compares the run time of these methods against a brute force approach across small Hamming graphs. Combining the techniques from sections 3 and 4, Section 6 presents a simple approach to discovering and removing redundant nodes in a given resolving set. This approach allows us to improve on previous bounds on the metric dimension of the Hamming graph $\mathbb{H}_{8,20}$. Finally, two appendices provide background information about Gröbner bases and linear programming.

All code used in this manuscript is available on GitHub (https://github.com/hamming-graph-resolvability/Hamming_Resolvability).

2 Main results

In what follows ${\hbox{Tr}}(A)$ denotes the trace of a square matrix $A$, $B^{\prime}$ the transpose of a matrix or vector $B$, and ${\hbox{vec}}(C)$ the column-major ordering of a matrix $C$ i.e. the row vector obtained by appending from left to right the entries in each column of $C$. For instance:

[TABLE]

In addition, $\bar{D}$ denotes the flip of the entries in a binary matrix (or vector) $D$, that is 0 is mapped to 1, and vice versa.

The one-hot encoding of a $k$-mer $v$ is defined as the binary matrix $V$ of dimension $(a\times k)$ such that $V[i,j]=1$ if and only if $(i-1)=v[j]$ (the offset in $i$ is needed since the reference alphabet is $\{0,...,a-1\}$ instead of $\{1,\ldots,a\}$). Here, $V[i,j]$ denotes the entry in row-$i$ and column-$j$ of the matrix $V$, and similarly $v[j]$ denotes the $j$-th coordinate of the vector $v$. We also follow the convention of capitalizing $k$-mer names to denote their one-hot encodings.

Our first result links one-hot encodings of $k$-mers with their Hamming distance. Note this result applies to any alphabet size, not just binary.

Lemma 2.1**.**

*If $u,v$ are $k$-mers with one-hot encodings $U,V$, respectively, then $d(u,v)=k-{\hbox{Tr}}(U^{\prime}V)$; in particular, $d(u,v)={\hbox{Tr}}(U^{\prime}\bar{V})$. *

Proof 2.2**.**

*Let $U_{i}$ and $V_{i}$ be the $i$-th column of $U$ and $V$, respectively. Clearly, if $u[i]=v[i]$ then $\langle U_{i},V_{i}\rangle=1$, and if $u[i]\neq v[i]$ then $\langle U_{i},V_{i}\rangle=0$, because all but one of the entries in $U_{i}$ and $V_{i}$ vanish and the non-vanishing entries are equal to 1. As a result, ${\hbox{Tr}}(U^{\prime}V)=\sum_{i=1}^{k}\langle U_{i},V_{i}\rangle$ counts the number of positions where $u$ and $v$ are equal; in particular, $d(u,v)=k-{\hbox{Tr}}(U^{\prime}V)$. Finally, observe that if $1^{a\times k}$ denotes the $(a\times k)$ matrix with all entries equal to 1 then ${\hbox{Tr}}(U^{\prime}1^{a\times k})=k$ because every row of $U^{\prime}$ has exactly one 1 and all other entries vanish. As a result, $d(u,v)={\hbox{Tr}}(U^{\prime}(1^{a\times k}-V))={\hbox{Tr}}(U^{\prime}\bar{V})$, as claimed. *

We can now give a necessary and sufficient condition for a subset of nodes in an arbitrary Hamming graph to be resolving.

Theorem 2.3**.**

Let $v_{1},\ldots,v_{n}$ be $n\geq 1$ $k$-mers and $V_{1},\ldots,V_{n}$ their one-hot encodings, respectively, and define the $(n\times ak)$ matrix with rows

[TABLE]

*Then $R:=\{v_{1},\ldots,v_{n}\}$ resolves $\mathbb{H}_{k,a}$ if and only if [math] is the only solution to the linear system $Az=0$, with $z$ a column vector of dimension $ak$, satisfying the following constraints: if $z$ is parsed into $k$ consecutive but non-overlapping subvectors of dimension $a$, namely $z=((z_{1},\ldots,z_{a}),(z_{a+1},\ldots,z_{2a}),...,(z_{(k-1)a+1},\ldots,z_{ka}))^{\prime}$, then each subvector is the difference of two canonical vectors. *

Proof 2.4**.**

Before showing the theorem observe that, for any pair of matrices $A$ and $B$ of the same dimension, ${\hbox{Tr}}(A^{\prime}B)=\langle{\hbox{vec}}(A),{\hbox{vec}}(B)\rangle$, where $\langle\cdot,\cdot\rangle$ is the usual inner product of real vectors.

*Consider $k$-mers $x$ and $y$, and let $X$ and $Y$ be their one-hot encodings, respectively. Due to Lemma 2.1, $d(v_{i},x)=d(v_{i},y)$ if and only if ${\hbox{Tr}}(V_{i}^{\prime}(X-Y))=0$ i.e. $\langle{\hbox{vec}}(V_{i}),{\hbox{vec}}(X-Y)\rangle=0$. As a result, the set $R$ does not resolve $\mathbb{H}_{k,a}$ if and only if there are $k$-mers $x$ and $y$ such that $Az\neq 0$, where $z:={\hbox{vec}}(X)-{\hbox{vec}}(Y)$. Note however that each column $X$ and $Y$ equals a canonical vector in $\mathbb{R}^{a}$; in particular, if we parse ${\hbox{vec}}(X)$ and ${\hbox{vec}}(Y)$ into $k$ subvectors of dimension $a$ as follows: ${\hbox{vec}}(X)=(x_{1},\ldots,x_{k})$ and ${\hbox{vec}}(Y)=(y_{1},\ldots,y_{k})$, then $z=(x_{1}-y_{1},\ldots,x_{k}-y_{k})$ with $x_{i}^{\prime}$ and $y_{i}^{\prime}$ canonical vectors in $\mathbb{R}^{a}$. This shows the theorem. *

2.1 Illustrative Example

In $H_{2,3}$ consider the set of nodes $R_{0}=\{02,11\}$. From Theorem 2.3, $R_{0}$ resolves $H_{2,3}$ if and only if $A_{0}z=0$, with

[TABLE]

has no non-trivial solution $z$ which satisfies the other constraints in the theorem when writing z=\big{(}(z_{1},z_{2},z_{3}),(z_{4},z_{5},z_{6})\big{)}^{\prime}. Something useful to note about this decomposition is that if a subvector of $z$ has two identical entries, then all the entries in that subvector must vanish.

Note that $A_{0}$ is already in its reduced row echelon form [26], and has two pivots: $z_{1}=-z_{6}$ and $z_{2}=-z_{5}$. Seeking non-trivial solutions to the constrained linear system, we examine permissible values for $z_{5}$ and $z_{6}$:

• (a)

If $z_{5}=-1$ then we must have $(z_{4},z_{6})\in\{(0,1),(1,0)\}$. Furthermore, if $z_{6}=1$ then $(z_{1},z_{2},z_{3})=(-1,1,0)$, but if $z_{6}=0$ then $(z_{1},z_{2},z_{3})=(0,1,-1)$. Consequently, $z=(-1,1,0,0,-1,1)$ and $z=(0,1,-1,1,-1,0)$ solve the constrained system.

• (b)

Similarly, we find that $z=(-1,0,1,-1,0,1)$ and $z=(1,0,-1,1,0,-1)$ solve the constrained system when we assume that $z_{5}=0$.

• (c)

Finally, $z=(1,-1,0,0,1,-1)$ and $z=(0,-1,1,-1,1,0)$ are also found to solve the constrained system when we impose that $z_{5}=1$.

Having found at least one non-trivial solution to the constrained linear system, we conclude that $R_{0}$ does not resolve $H_{2,3}$. (The found $z$’s are in fact the only non-trivial solutions.)

From the proof of Theorem 2.3, we can also determine pairs of vertices in $H_{2,3}$ which are not resolved by $R_{0}$. Indeed, using the non-trivial solutions found above we find that $12$ and $01$, $21$ and $10$, and $00$ and $22$ are the only pairs of nodes in $H_{2,3}$ which are unresolved by $R_{0}$. In particular, because the distances between the nodes in each pair and $22$ are different, $R_{1}:=R_{0}\cup\{22\}$ resolves $H_{3,2}$.

We can double-check this last assertion noticing that the reduced echelon form of the matrix $A_{1}$ associated with $R_{1}$ is

[TABLE]

In particular, $z_{1}=-z_{6}$, $z_{2}=-z_{5}$, and $z_{3}=-z_{6}$. The first and third identity imply that $z_{1}=z_{2}$, hence $(z_{1},z_{2},z_{3})=(0,0,0)$. This together with the first and second identity now imply that $(z_{4},z_{5},z_{6})=(0,0,0)$. So, as anticipated, $z=0$ is the only solution to the constrained linear system $A_{1}z=0$.

In general, if the reduced row echelon form of the matrix given by Theorem 2.3 has $j$ free variables, then there could be up to $3^{j}$ possible solutions to the associated linear system, each of which would have to be checked for the additional constraints. This exhaustive search could be very time consuming if not impossible. Handling the linear system constraints more systematically and efficiently is the motivation for sections 3 and 4.

2.2 Specializations to Hypercubes

In [1] a necessary and sufficient condition for the resolvability of hypercubes is provided exploiting that $d(u,v)=\|u-v\|_{2}^{2}$ when $u$ and $v$ are binary $k$-mers. Next, we reproduce this result using our framework of one-hot encodings instead.

Corollary 2.5**.**

[1, Theorem 2.2]** Let $R=\{v_{1},\ldots,v_{n}\}$ be a set of $n\geq 1$ binary $k$-mers, and define the $(n\times k)$ matrix with rows

[TABLE]

*Then, $R$ resolves $H_{k,2}$ if and only if $\hbox{ker}(B)\cap\{0,\pm 1\}^{k}=\{0\}$.

Proof 2.6**.**

Let

[TABLE]

be the $(n\times 2k)$ matrix with columns $A_{1},\ldots,A_{2k}$ given by Theorem 2.3 for $R$. It follows that $R$ resolves $H_{k,2}$ if and only if $Az=0$, with $z=((x_{1},y_{1}),\ldots,(x_{k},y_{k}))^{\prime}\in\{0,\pm 1\}^{2k}$ and $(x_{i}+y_{i})=0$ for each $i=1,\ldots,k$, has only a trivial solution. Note however that $Az=By$, where

[TABLE]

Therefore $R$ is resolving if and only if $By=0$, with $y\in\{0,\pm 1\}^{k}$, has only a trivial solution. But recall from Theorem 2.3 that the rows of $A$ are the column-major orderings of the one-hot encodings of the binary $k$-mers in $R$. In particular, using $\llbracket\cdot\rrbracket$ to denote Iverson brackets, we find that the row in $B$ associated with $v\in R$ is:

[TABLE]

*from which the corollary follows. *

We can provide an even simpler characterization of sets of $k$-mers that resolve the hypercube, provided that $1^{k}:=(1,\ldots,1)$ is one of them. This seemly major assumption is only superficial. Indeed, hypercubes are transitive; that is, given any two binary $k$-mers there is an automorphism (i.e., a distance preserving bijection $\sigma:\{0,1\}^{k}\to\{0,1\}^{k}$) that maps one into the other [29, §3.1]. Hence, given any set $R$ of binary $k$-mers there is an automphism $\sigma$ such that $1^{k}\in\sigma(R)$. In particular, because $R$ is resolving if and only if $\sigma(R)$ is resolving, one can assume without any loss of generality that $1^{k}$ is an element of $R$.

Corollary 2.7**.**

Let $R=\{v_{1},\ldots,v_{n}\}$ be a set of $n$ binary $k$-mers such that $1^{k}\in R$, and define the $(n\times k)$ matrix with rows

[TABLE]

*Then, $R$ resolves $H_{k,2}$ if and only if $\hbox{ker}(C)\cap\{0,\pm 1\}^{k}=\{0\}$.

Proof 2.8**.**

Note that for all binary $k$-mers $v$: $(v+\bar{v})=1^{k}$; in particular, $(v-\bar{v})=(2v-1^{k})$. Hence, if $B$ is as given in Corollary 2.7 and $C$ as defined above then

[TABLE]

*But, because $1^{k}\in R$ and $2\cdot 1^{k}-1^{k}=1^{k}$, one of the entries in $Bz$ equals $\langle 1^{k},z\rangle$. Since the entries in $Cz$ are proportional to $\langle 1^{k},z\rangle$, $Bz=0$ if and only if $Cz=0$, from which the corollary follows. *

3 Polynomial Roots Formulation

In this section, we express the constraints of the linear system in Theorem 2.3 as roots of a multi-variable polynomial system, and we reveal various properties about this system which can drastically reduce the complexity of determining whether a subset of nodes resolves or not $\mathbb{H}_{k,a}$.

In what follows, for any given non-empty set $P$ of polynomials in a possibly multi-variable $z$, $\{P=0\}$ denotes the set of $z$’s such that $p(z)=0$, for each $p\in P$. Unless otherwise stated, we assume that $z$ has dimension $ka$, i.e. $z=(z_{1},\ldots,z_{ka})$.

Consider the polynomial sets

[TABLE]

Our first result characterizes the constraints of the linear system in Theorem 2.3 in terms of the roots of the above polynomials. Ahead, unless otherwise stated:

[TABLE]

Lemma 3.1**.**

$z\in\{P=0\}$* if and only if when parsing $z$ into $k$ consecutive but non-overlapping subvectors of dimension $a$, each subvector is the difference of two canonical vectors.

Proof 3.2**.**

*The polynomials in $P_{1}$ enforce that each entry in $z$ must be a ${-1}$, [math], or $1$, while the polynomials in $P_{2}$ enforce that there is a $(-1)$ for every $1$ in each subvector of $z$. Finally, the polynomials in $P_{3}$ enforce that each subvector of $z$ has exactly two non-zero entries or no non-zero entries. Altogether, $z\in P$ if and only if each subvector is identically zero, or it has exactly one 1 and one $(-1)$ entry and all other entries vanish, i.e. each subvector of $z$ is the difference of two canonical vectors in $\mathbb{R}^{a}$. *

The following is now an immediate consequence of this lemma and Theorem 2.3.

Corollary 3.3**.**

*Let $R$ be a set of nodes in $\mathbb{H}_{k,a}$ and $A$ the matrix given by equation Eq. 2. Then, $R$ resolves $\mathbb{H}_{k,a}$ if and only if $\hbox{ker}(A)\cap\{P=0\}=\{0\}$.

Our focus in what remains of this section is to better characterize the non-trivial roots of the polynomial system $\{P=0\}$. To do so, we rely on the concepts of polynomial ideals and (reduced) Gröbner bases, and the following fundamental result from algebraic geometry. For a primer to these and other concepts on which our results rely see Appendix A.

Theorem 3.4**.**

*(Hilbert’s Weak Nullstellensatz [8, §4.1].)

For any non-empty finite set of polynials $P$, $\{P=0\}=\emptyset$ if and only if $\{1\}$ is the reduced Gröbner basis of $I(P)$, the ideal generated by $P$. *

Define for each $i=1,\ldots,k$:

[TABLE]

Observe that $B_{i}$ is a set of polynomials in $(z_{(i-1)a+1},\ldots,z_{ia})$, i.e. the $i$-th subvector of $z$; in particular, each of these polynomials may be regarded as a function of $z$, and $B_{1},\ldots,B_{k}$ partition $P$, i.e. $P=\sqcup_{i=1}^{k}B_{i}$. Accordingly, we call $B_{i}$ the $i$-th block of $P$, and denote the reduced Gröbner basis of $B_{i}$ as $G_{i}$. The computational advantage of these observations is revealed by the following results.

Lemma 3.5**.**

$G=\cup_{i=1}^{k}G_{i}$* is the reduced Gröbner bases of $P$ in equation Eq. 7. Furthermore, $G_{i}$ may be obtained from $G_{1}$ using the change of variables:*

[TABLE]

Proof 3.6**.**

*The case with $k=2$ follows from [8, Proposition 2] due to the fact that no variable and hence no polynomial is shared between the blocks of $P$. A straightforward inductive argument in $k\geq 2$ then shows that $\cup_{i=1}^{k}G_{i}$ is the reduced Gröbner bases of $P$. Finally, note that $B_{1}$ is up to the change of variables in equation (9) identical to $B_{i}$; in particular, since Buchberger’s algorithm (Algorithm 1) and the Gröbner basis reduction algorithm (Algorithm 2) build upon polynomial division, the reduced Gröbner bases of $B_{i}$ may be obtained from that of $B_{1}$ using the same change of variables. *

Lemma 3.7**.**

The reduced Gröbner bases of $B_{1}$ under the lexicographic ordering is

[TABLE]

Proof 3.8**.**

Let $G$ be the set of polynomials on the right-hand side above. Since $G$ depends on $a$ but not on the parameter $k$ of $\mathbb{H}_{k,a}$, and the identity for $a\in\{2,3,4\}$ can be checked using algorithms 1 and 2, without loss of generality we assume in what follows that $a>5$.

Since reduced Gröbner basis are unique, it suffices to show that (i) $I(G)=I(B_{1})$; and that for all $f,g\in G$: (ii) the reduction of $\hbox{Spoly}(f,g)$ by $G$ is 0; (iii) $LC(f)=1$; and (iv) if $f\in G\setminus\{g\}$ then no monomial of $f$ is divisible by $LM(g)$. We omit the tedious but otherwise straightforward verification of properties (ii) and (iv). Since property (iii) is trivially satisfied, it only remains to verify property (i).

*To prove $I(G)=I(B_{1})$, it suffices to show that $\{G=0\}=\{B_{1}=0\}$. Indeed, the polynomials of the form $z_{i}z_{j}z_{\ell}$ imply that if $z\in\{G=0\}$ then $(z_{2},\ldots,z_{a-1})$ has at most two non-zero coordinates. In the case two of these coordinates are non-zero, say $z_{i}$ and $z_{j}$, the polynomials $z_{i}^{3}-z_{i}=z_{i}(z_{i}-1)(z_{i}+1)$, $z_{j}^{3}-z_{j}=z_{j}(z_{j}-1)(z_{j}+1)$, and $z_{i}^{2}z_{j}+z_{i}z_{j}^{2}=z_{i}z_{j}(z_{i}+z_{j})$ imply that $(z_{i},z_{j})=(1,-1)$ or $(z_{i},z_{j})=(-1,1)$; in particular, because we must have $\sum_{\ell=1}^{a}z_{\ell}=0$, $z_{1}=0$. Instead, if exactly one of the coordinates in $(z_{2},\ldots,z_{a-1})$ is non-zero, say $z_{j}$, then the polynomial $\sum_{\ell=1}^{a}z_{\ell}$ together with $z_{j}^{3}-z_{j}$ imply that $(z_{1},z_{j})=(1,-1)$ or $(z_{1},z_{j})=(-1,1)$. Finally, if $(z_{2},\ldots,z_{a-1})=0$ then the polynomial $\sum_{\ell=1}^{a}z_{\ell}$ implies that $z_{1}=0$. In all of these three exhaustive cases, it follows that $(z_{1},\ldots,z_{a})$ is identically zero, or it has exactly one 1 and one (-1) coordinate and all other coordinates vanish; in other words, $(z_{1},\ldots,z_{a})$ is a difference of two canonical vectors in $\mathbb{R}^{a}$. Since this is precisely the constraint imposed on this subvector of $z$ by the polynomials in $B_{1}$, we obtain that $\{G=0\}=\{B_{1}=0\}$ i.e. $I(G)=I(B_{1})$. *

A minor issue for using the Weak Nullstellensatz in our setting is that the polynomials in $P$ have no constant terms; in particular, $0\in\{P=0\}$. To exclude this trivial root, observe that if $z\in P$ then $\sum_{j=(i-1)a+1}^{ia}z_{j}^{2}\in\{0,2\}$, for each $i=1,\ldots,k$. As a result, if $z$ is a non-trivial root of $\{P=0\}$ then $\sum_{j=1}^{ka}z_{j}^{2}=2i$ for some $i$. This motivates to introduce the auxiliary polynomial:

[TABLE]

so that $R$ resolves $\mathbb{H}_{k,a}$ if and only if $\hbox{ker}(A)\cap\{P=0\}\cap\{f-2i=0\}=\emptyset$ for all $i=1,\ldots,k$.

Lemma 3.9**.**

*Consider a (finite) reduced Gröbner basis $G\neq\{1\}$ and a polynomial $f$. If $f\xrightarrow{G}r$ then, for each $c\in\mathbb{R}$, $(f+c)\xrightarrow{G}(r+c)$. *

Proof 3.10**.**

*Let $G=\{g_{1},\ldots,g_{n}\}$. Without loss of generality fix a constant $c\neq 0$. Note that $G$ contains no constant polynomial (except for [math]) because $G\neq\{1\}$ hence $1\notin G$. As a result, the leading monomial of each $g_{i}$ does not divide $c$, hence $c\xrightarrow{G}c$. Since $f\xrightarrow{G}r$, and reductions by a Gröbner basis are unique, $(f+c)\xrightarrow{G}(r+c)$ as claimed. *

The following is now a direct consequence of the lemma.

Corollary 3.11**.**

*Let $G$ be the reduced Gröbner basis of $P$ in equation (7). If $f$ is as defined in Eq. 11 and $f\xrightarrow{G}r$ then, for each $i=1,2,\ldots,k$, $(f-2i)\xrightarrow{G}(r-2i)$.

The results from this section allow for a computational method for checking resolvability on $\mathbb{H}_{k,a}$. Lemmas 9 and 10 are used to construct the reduced Gröbner basis $G$ directly, and Corollary 3.11 efficiently removes the trivial solution from consideration in the criteria provided by Theorem 3.4. Altogether these results significantly reduce the number of polynomial reductions required to assess the resolvability of a set of nodes on $\mathbb{H}_{k,a}$.

3.1 Illustrative Example (Continuation)

We saw in Section 2.1 that $R_{0}=\{02,11\}$ does not resolve $H_{2,3}$ whereas $R_{1}=R_{0}\cup\{22\}=\{02,11,22\}$ does. We can double-check these assertions using Corollary 3.3 as follows.

First, recall that for $H_{2,3}$ the variable $z$ is 6-dimensional and should be decomposed in the form z=\big{(}(z_{1},z_{2},z_{3}),(z_{4},z_{5},z_{6})\big{)}. Next, the kernel of the matrix given by the corollary for $R_{0}$ (denoted as $A_{0}$, see Eq. Eq. 3) is described by the linear system:

[TABLE]

On the other hand, the roots in $\{P=0\}$ given by Corollary 3.3 correspond to the polynomial system:

[TABLE]

where the horizontal line distinguishes between the first and second block of $P$ (see Eq. (3)). Finally, recall the auxiliary polynomial given by equation Eq. 11:

[TABLE]

Assuming the lexicographic order over the monomials, one can determine that the reduced Gröbner basis of $\{A_{0}z\}\cup P\cup\{f-2\}$ is $\{1\}$; in particular, $\hbox{ker}(A_{0})\cap\{P=0\}\cap\{f=2\}=\emptyset$. On the other hand, because the reduced Gröbner basis of $\{A_{0}z\}\cup P\cup\{f-2\}$ is $\{z_{1}+z_{6},z_{2}+z_{5},z_{3}-z_{5}-z_{6},z_{4}+z_{5}+z_{6},z_{5}^{2}+z_{5}z_{6}+z_{6}^{2}-1,z_{6}^{3}-z_{6}\}$, it follows that $\hbox{ker}(A_{0})\cap\{P=0\}\cap\{f=4\}\neq\emptyset$ i.e. $\hbox{ker}(A_{0})\cap\{P=0\}$ has a non-trivial solution. Consequently, $R_{0}$ does not resolve $H_{2,3}$.

To confirm that $R_{1}=R_{0}\cup\{22\}$ does resolve $H_{2,3}$, note that we only need to add the equation $z_{3}+z_{6}=0$ to the previous linear system (the full linear system is now described by the matrix $A_{1}$, see Eq. Eq. 4). Using our code, we find that $\hbox{ker}(A_{1})\cap\{P=0\}\cap\{f-2\}=\emptyset$ and also that $\hbox{ker}(A_{1})\cap\{P=0\}\cap\{f-4\}=\emptyset$ because the associated reduced Gröbner basis are both equal to $\{1\}$. As a result, $\hbox{ker}(A_{1})\cap\{P=0\}$ has no non-trivial solution, i.e. $R_{1}$ resolves $H_{2,3}$.

4 Novel Integer Linear Programming Formulation

For some background about Integer Linear Programming (ILP), see Appendix B.

In contrast to the ILP approaches of [5, 9], our ILP formulation checks the resolvability of a given set rather than searching for minimal resolving sets. Furthermore, it does not pre-compute the distance matrix of a Hamming graph. As before, fix $\mathbb{H}_{k,a}$ and a subset of vertices $R$. Letting $z=(z_{1},\ldots,z_{ka})$ and using the polynomial set $P$ from equation (7), we leverage Lemma 3.1 (with $A$ as in equation (2), each row corresponding to a vertex in $R$) to reformulate Theorem 2.3 as follows:

[TABLE]

To formulate this as an ILP, we use the following result.

Lemma 4.1**.**

Define

[TABLE]

*Then $\mathcal{I}$ is the intersection of a closed convex polyhedron with the integer lattice $\mathbb{Z}^{ak}$, and $z\in\{P=0\}$ if and only if $z\in\mathcal{I}$.

Proof 4.2**.**

Since the intersection of convex sets is convex, and the intersection of a finite number of polyhedra is a polyhedron, it follows from standard arguments that

[TABLE]

are convex subsets of $\mathbb{R}^{ak}$, and $\mathcal{J}_{1}$ is a polyhedron. We claim that $\mathcal{J}_{2}$ is also a polyhedron, for which it suffices to check that each set in the intersection that defines it is a polyhedron. Without loss of generality, we do so only for the case with $i=1$. Indeed, because $\{z\in\mathbb{R}^{ak}\hbox{ such that }\sum_{j=1}^{a}|z_{j}|\leq 2\}$ is invariant under arbitrary coordinate sign flips, we have that

[TABLE]

which implies that $\mathcal{J}_{2}$ is also a polyhedron. Since $\mathcal{I}=(\mathcal{J}_{1}\cap\mathcal{J}_{2}\cap\mathbb{Z}^{ak})$, the first part of the lemma follows.

From the proof of Lemma 3.1 it is immediate that $\{P=0\}\subset\mathcal{I}$. To show the converse inclusion, observe that $\{P=0\}=\cap_{i=1}^{k}\{B_{i}=0\}$ where the $B_{i}$’s are as defined in equation (3). To complete the proof, it suffices therefore to show that $\mathcal{I}_{i}\subset\{B_{i}=0\}$, where

[TABLE]

*Indeed, if $z\in\mathcal{I}_{1}$ then, because the coordinates of $z$ are integers, the condition $\sum_{j=1}^{a}|z_{j}|\leq 2$ implies that $|z_{j}|\in\{0,1,2\}$ for $j=1,\ldots,a$. If $|z_{j}|=2$ for some $j$ then $\sum_{j=1}^{a}z_{j}=\pm 2$, which is not possible. Thus $z_{j}\in\{0,\pm 1\}$ for $j=1,\ldots,a$; in particular, $z_{j}^{3}-z_{j}=0$. On the other hand, the condition $\sum_{j=1}^{a}z_{j}=0$ implies that the number of 1’s and (-1)’s in $(z_{1},\ldots,z_{a})$ balance out; in particular, since $\sum_{j=1}^{a}|z_{j}|\leq 2$, either $(z_{1},\ldots,z_{a})$ vanishes, or it has exactly one 1 and one (-1) entry and all other entries vanish; in particular, $(2-\sum_{j=1}^{a}z_{j}^{2})\cdot\sum_{j=1}^{a}z_{j}^{2}=0$. Thus, $z\in\{B_{1}=0\}$. The case for $i>1$ is of course the same. *

Remark 4.3**.**

*With current ILP solvers, one can impose that $z\in\{0,\pm 1\}^{ak}$ simply as $|z_{i}|\leq 1$ for $i=1,\ldots,ak$. On the other hand, while a constraint like $\sum_{j=1}^{a}|z_{j}|\leq 2$ is clearly polyhedral, it is not in the form of an affine equality or inequality suitable for ILP solvers. Nevertheless, standard reformulation techniques can convert this into a set of affine equalities and inequalities in a higher dimensional space. For example, in the product space with variables $(\tilde{z},w)$, we can write the constraint as $\sum_{j=1}^{a}w_{j}\leq 2$ and $|\tilde{z}_{j}|\leq w_{j}$ (i.e., $\tilde{z}_{j}\leq w_{j}$ and $-\tilde{z}_{j}\leq w_{j}$), which leads to an equivalent formulation of the original ILP. One may handle such reformulations automatically using the Matlab package CVX [14]. *

It only remains to encode the fact that we look for a nonzero root in $\{P=0\}$, which we do via the ILP in the following theorem:

Theorem 4.4**.**

A subset of vertices $R$ is not resolving on $\mathbb{H}_{k,a}$ if and only if the solution to the following ILP is less than zero:

[TABLE]

*where $A$ is defined in equation Eq. 2. *

Proof 4.5**.**

*Using equation Eq. 12 and Lemma 4.1, it remains to show that the objective function is less than zero if and only if there is a non-zero feasible $z$. Suppose there is not a non-zero feasible $z$. Clearly $z=0$ is feasible, hence it is the only feasible point for the ILP, and the objective value is zero. Now suppose there is some non-zero feasible $z$. Let $j^{\prime}$ be the largest non-zero coordinate. Then because $\sum_{j=1}^{j^{\prime}-1}2^{j}<2^{j^{\prime}}$, and because each entry is bounded $|z_{j}|\leq 1$, the objective value at this $z$ is non-zero. If the objective value is negative, this proves the value of the ILP is negative; if the objective value is positive, then observe that $(-z)$ is also feasible and has a negative objective value, and hence the value of the ILP is negative. *

Remark 4.6**.**

*If the solution to the ILP is less than zero and hence $R$ is not a resolving set, then each optimal vector $z$ is the difference of the column-major ordering of the one-hot encodings of two $k$-mers which are not resolved by $R$; in particular, a vector that resolves these $k$-mers needs to be added to $R$ to resolve $\mathbb{H}_{k,a}$. *

4.1 Practical formulations and roundoff error

When $ak$ is small, it is feasible to directly solve the ILP in equation Theorem 4.4. One issue with larger values of $ak$, besides an obvious increase in run-time, is that the values of $2^{j}$ in the objective function quickly lead to numerical overflow. A simple fix is to replace each coefficient $c_{j}=2^{j}$ with an independently drawn realization of a standard normal random variable $\mathcal{N}(0,1)$. Since these new coefficients are independent of the feasible set, if the latter is truly larger than $\{0\}$, the probability that the entire feasible set is in the null-space of the linear function $\sum_{j=1}^{ak}c_{j}z_{j}$ is zero. Of course, again due to finite machine precision, this otherwise almost surely exact method may only be approximate. Admittedly, when running the ILP with the random coefficients $c_{j}$’s, finding an undoubtedly negative solution to the ILP would certify that the set $R$ is not resolving. However, if the solution is just slightly negative or vanishes within machine precision, the assessment about $R$ should be taken with a grain of salt. In this case, one should draw a new set of random coefficients and re-run the ILP to reassess the resolvability of $R$.

Another consideration is that the ILP solver wastes time finding a feasible point with the smallest possible objective, when we only care if there is a feasible point with objective smaller than [math]. Thus we could solve the feasiblity problem

[TABLE]

where $c_{j}=2^{j}$ or $c_{j}\sim\mathcal{N}(0,1)$ as discussed above. (Feasibility problems can be encoded in software by minimizing the [math] function.) Unfortunately this is not an ILP because $\{c\mid\langle c,z\rangle<0\}$ is not a closed set. We can partially ameliorate this by solving

[TABLE]

where $\delta>0$ is a small number (our code uses $\delta=10^{-3}$). Finding a feasible point $z$ is then proof that the set $R$ does not resolve $\mathbb{H}_{k,a}$. If the solver says the above problem is infeasible, it could be that $\delta$ was too large and hence the computation was inconclusive. In this case, one could run the slower program Theorem 4.4.

5 Computational Complexity Experiments

The theoretical framework and algorithms proposed in this paper provide a novel way of approaching resolvability on Hamming graphs. To show the computational feasibility and practicality of our methods, we compare the average run-time of both the ILP and Gröbner basis algorithms against the brute force approach for checking resolvability. Our experiments use Python 3.7.3 and SymPy version 1.1.1 [23], and the commercial ILP solver gurobi ver. 7.5.2 [16].

In the Table 1, we present the average run-time and standard deviation of the algorithms on reference test sets for Hamming graphs of increasing sizes. Fig. 2 displays the mean run-times as a function of the graph size, and best linear fit for each method. As seen in the table and figure, the brute force approach is faster on only the smallest Hamming graphs (with fewer than $1000$ nodes) whereas the ILP solution is exceptionally fast even as the Hamming graph grows to more than $6000$ nodes. For small problems, the time taken to solve the ILP is likely dominated by the overhead cost of using CVX to recast the ILP into standard form. The run-time results show a promising improvement in computational time over the brute force approach which will only become more pronounced on massive Hamming graphs. Additionally, the brute force approach is infeasible on these larger graphs due to significant memory costs.

The ILP algorithm is exceptionally quick, beating all other methods for Hamming graphs with more than 1000 nodes, but it cannot guarantee that a set is resolving. The Gröbner basis algorithm by contrast is much slower on average but is a deterministic method of showing resolvability. ILP can be used to quickly determine possible resolving sets which are then verified by the Gröbner basis algorithm. In this way, the two methods are symbiotic and cover each other’s weaknesses. We illustrate this in the next section.

6 Low-dimensional Protein Representations

Symbolic information pervades modern data science. With the advent and popularization of high-throughput sequencing assays, this is particularly true in the field of computational biology where large volumes of biological sequence data have become critical for studying and understanding the behavior of cells. Analysis of these sequences, however, presents significant challenges. One major issue is that many powerful analysis techniques deal with numeric vectors, not arbitrary symbols. As a result, biological sequence data is typically mapped to a real space before such methods are applied. Two of the most common mappings use K-mer count [22] and one-hot encodings (also called binary vectors) [4]. K-mer count vectors represent symbolic sequences by their counts of each possible K-mer.

Resolving sets can be used to define low-dimensional mappings as well. To fix ideas we focus on octapeptides, that is proteins composed of 8 amino acids. With a total of 20 possible amino acids (which we represent as a,r,n,d,c,q,e,g,h,i,l,k,m,f, p,s,t,w,y,v) and imposing the Hamming distance across these sequences, we have the Hamming graph $\mathbb{H}_{8,20}$. This graph is massive. It has $25.6$ billion vertices and roughly $1.9$ trillion edges rendering most methods of discovering small resolving sets, including the ICH algorithm, useless. Utilizing a constructive algorithm, a resolving set of size 82, which we call $R$, was discovered for $\mathbb{H}_{8,20}$ in [29]. However, it is not known whether $R$ contains a proper subset that is still resolving. Here, we address this problem applying the results of sections 3 and 4.

Starting with lower and upper bounds $L=1$ and $U=82$ respectively, we implement a binary search for $\beta(\mathbb{H}_{8,20})$. With $s=\frac{L+U}{2}$ as the current subset size to check, up to 1000 subsets of $R$ are selected at random. The ILP approach (Section 4) then provides an efficient method for testing the feasibility problem outlined in Theorem 2.3 for these subsets. If any subset passes this test, the upper bound is set to $s$. Otherwise, $s$ becomes the lower bound. This process is repeated until $L=(U-1)$. Following this procedure, we found the following set of size $77$:

[TABLE]

Since the ILP formulation does not guarantee that this set is resolving, we verified the result using a parallelized version of the Polynomial Roots Formulation (Section 3) so that the Gröbner bases of multiple auxiliary polynomials (Eq. Eq. 11) may be determined simultaneously. Thus, we have found a set $r\subset R$ of size 77 that resolves $\mathbb{H}_{8,20}$; in particular, $\beta(\mathbb{H}_{8,20})\leq 77$, which improves the bound of [29], and all $25.6$ billion octapeptides may be uniquely represented with only 77 dimensions. In contrast, a $2$-mer count vector representation would require 400 dimensions and a one-hot encoding 160 dimensions.

Remark 6.1**.**

*We replicated the verification of $r$ as a resolving set of $H_{8,20}$ using our Polynomial Roots Formulation 10 times across 32 computer cores. Overall, a maximum of approximately 380 megabytes of memory per core (SD $\sim 0.5$ MB) and 6 hours and 20 minutes (SD $\sim 142$ s) were required to demonstrate the resolvability of $r$. Memory usage was determined using the Slurm workload manager sacct command and maxRSS field, while time was measured using Python’s time module. *

Appendix A Gröbner Basis

In what follows, $z=(z_{1},\ldots,z_{d})$ is a $d$-dimensional variable and, unless otherwise stated, polynomials are functions of $z$ with real coefficients.

A.1 Polynomial Ideals

A polynomial ideal $I$ is a non-empty set of polynomials with the property that if $f,g\in I$ and $c\in\mathbb{R}$ then $cf,f+g,fg\in I$.

The polynomial ideal associated with a non-empty and finite set $P=\{p_{1},\ldots,p_{m}\}$ of polynomials in the variable $z$ is the set defined as

[TABLE]

Polynomial ideals are useful to characterize the complex numbers $z$ such that $p_{i}(z)=0$ for all $i=1,\ldots,m$. Indeed, $z\in\{P=0\}$ if and only if $z\in\{I(P)=0\}$.

A.2 Monomial Orderings

A monomial (in the variable $z$) is any product of the form $z_{1}^{a_{1}}\cdots z_{d}^{a_{d}}$, where $a_{1},\ldots,a_{d}\geq 0$ are integers. This product is often written $z^{a}$ with $a=(a_{1},\ldots,a_{d})$.

A monomial ordering is a total ordering of the monomials such that (i) if $z^{a}<z^{b}$ then, for any monomial $z^{c}$, $z^{a+c}<z^{b+c}$; and (ii) $1<z^{a}$ when $z^{a}\neq 1$.

A common example of monomial ordering is the so-called lexicographic order. Under this ordering, $z^{a}<z^{b}$ if there is an index $i$ such that $a_{j}=b_{j}$ for all $1\leq j<i$ but $a_{i}<b_{i}$. Another example is the graded lexicographic order under which $z^{a}<z^{b}$ if either (i) $\sum_{i=1}^{d}a_{i}<\sum_{i=1}^{d}b_{i}$; or (ii) $\sum_{i=1}^{d}a_{i}=\sum_{i=1}^{d}b_{i}$ but $z^{a}$ is smaller than $z^{b}$ under the lexicographic order. Both of these orderings can be reversed giving the reversed lexicographic order and the graded reverse lexicographic order, respectively.

In what remains of Appendix A, a fixed monomial ordering is assumed. In particular, each non-zero polynomial $p$ may be uniquely written in the form $p=\sum_{i=1}^{t}c_{i}m_{i}$, where $t\geq 1$ is an integer, $c_{1},\ldots,c_{t}$ are real coefficients, and $m_{1}>\cdots>m_{t}$ are monomials in descending order. This allows us to define $\hbox{LM}(p):=m_{1}$ (the leading monomial of $p$), $\hbox{LC}(p):=c_{1}$ (the leading coefficient of $p$), and $\hbox{LT}(p):=c_{1}m_{1}$ (the leading term of $p$).

A.3 Polynomial Reductions

For a given non-empty set $P$ of polynomials, every polynomial $f$ can be represented in the form $f=(r+g)$, with $r,g$ polynomials such that $g\in I(P)$. Representing $f$ in this form is called reducing $f$ by $P$. The term $r$ is called the reduction of $f$ by $P$; which is expressed in writing as $f\stackrel{{\scriptstyle P}}{{\rightarrow}}r$. Observe that reductions are typically not unique. This is because, if $f=(r+g)$, with $g\in I(P)$, then $f=(r-h)+(g+h)$ for any polynomial $h$, however, $(g+h)\in I(P)$ when $h\in I(P)$.

Reductions can be computed using multivariate long-division as follows. First, set $r=0$ and $g=f$, and look for the smallest index $i$ such that $LM(p_{i})$ divides $LM(g)$. If such an index exists, set $g=g-LT(g)\cdot p_{i}/LT(p_{i})$ so that the $LT(g)$ and the $LT(p_{i})$ cancel. Otherwise, set $g=g-LT(g)$ and $r=r+LT(g)$. Continue this process until $g=0$. This will produce a remainder $r$ where no monomial of $r$ is divisible by any $LM(p_{i})$.

In general, reductions of the form $f\stackrel{{\scriptstyle P}}{{\rightarrow}}r$ with $r\notin I(P)$ are also not unique. This is because the polynomial $r$ obtained by the long-division depends on the order in which the polynomials in $P$ are indexed. This lack of uniqueness is the primary motivation for Gröbner bases.

A.4 Buchberger’s Criterion

In this section we give a characterization of Gröbner bases due to Buchberger [7].

The least common multiple between $z^{a}$ and $z^{b}$ is the monomial $\hbox{LCM}(z^{a},z^{b}):=z^{c}$ with $c=(\max\{a_{1},b_{1}\},\ldots,\max\{a_{d},b_{d}\})$.

Given two polynomials $p$ and $q$ such that $\hbox{LCM}(\hbox{LM}(p),\hbox{LM}(q))=z^{c}$, their S-polynomial is defined as

[TABLE]

A non-empty set $G=\{g_{1},\ldots,g_{n}\}$ of polynomials is called a Gröbner basis for a polynomial ideal $I$ if (i) $I(G)=I$, and (ii) $\hbox{Spoly}(g_{i},g_{j})\stackrel{{\scriptstyle G}}{{\rightarrow}}0$ for all $g_{i},g_{j}\in G$.

Gröbner bases have the following property with regards to reductions: $f\stackrel{{\scriptstyle G}}{{\rightarrow}}0$ if and only if $f\in I(G)$, otherwise $f\stackrel{{\scriptstyle G}}{{\rightarrow}}r$ for some $r\notin I(G)$.

A Gröbner basis $G$ is called reduced if for all $g_{i}\in G$, (iii) $LC(g_{i})=1$, and (iv) for all $g_{j}\in G\setminus\{g_{i}\}$, no monomial of $g_{j}$ is divisible by $LM(g_{i})$.

Unlike Gröbner bases, the reduced Gröbner basis of a polynomial ideal is unique.

A.5 Buchberger’s Algorithm

This is a method for generating a Gröbner basis for a polynomial ideal $I(P)$ based on Buchberger’s criterion (see Algorithm 1). The key idea of the algorithm is to add to the initially empty Gröbner basis S-polynomials of pairs in $P$ which do not reduce to [math]. This by construction satisfies Buchberger’s criterion and hence produces a Gröbner basis.

Buchberger’s Algorithm, however, does not necessarily produce a reduced Gröbner basis. Such a reduction can be computed using Algorithm 2. This algorithm is the simplest but also least efficient for computing Gröbner bases.

The computation of Gröbner bases is an active field of research with many different approaches. There are matrix reduction based algorithms, such as Faugére’s F4 algorithm [11], as well as signature based algorithms, like Faugére’s F5 algorithm and its variants F5C and F5B [12, 10, 28].

Appendix B Linear Programming

A linear program (LP) is any optimization problem that minimizes or maximizes a linear objective function over a (possibly unbounded) closed convex polyhedron in Euclidean space, i.e., a finite number of affine equalities and inequalities (but not strict inequalities). In the combinatorial model of complexity, LP’s are known to be solvable in polynomial time [19]. If the variables are constrained to be integers, the LP is an integer linear program (ILP), or more generally, if some of the variables are constrained to be integers and others are not, it is a mixed integer linear program (MILP); one usually does not make the distinction between ILP and MILP, since all ILP solvers also solve MILP.

Unlike LPs, ILP can encode NP-Hard decision problems, and thus in the worst case they are intractable for large problems. However, because ILPs usually admit relaxation (e.g. dropping some constrains), one can find upper and lower bounds on the value of an ILP, and thus, when the bounds meet, produce a certificate of optimality. Standard techniques such as branch-and-bound produce a large tree of different relaxations of the ILP, each node a LP and thus solvable. Once all leaf nodes are visited, one has a guarantee of the optimal solution, but there may be combinatorially many leaf nodes. Modern ILP solvers use clever heuristics to determine the order in which to traverse the tree, and with a good choice, entire branches of the tree can quickly be pruned. Thus ILP solvers combine the speed of heuristics with provable optimality. Furthermore, because ILPs are common in industry, there are many high-quality commercial solvers with excellent implementations such as gurobi [16].

Because in the worst-case ILPs can take exponentially long to solve, they are sometimes overlooked in the mathematician’s toolbox. Nevertheless, because they are general, they are the focus of much research and thus the software to solve them continues to improve rapidly. To quote Dimitris Bertsekas, “in the last twenty-five years (1990–2014), algorithmic advances in integer optimization combined with hardware improvements have resulted in an astonishing 200 billion factor speedup in solving MIO [mixed-integer optimization] problems” [3]. Thus the main challenge is recognizing when a problem with a combinatorial flavor can be recast as an ILP.

Acknowledgments

This research was partially funded by NSF grant No. 1836914. The authors acknowledge the BioFrontiers Computing Core at the University of Colorado Boulder for providing High-Performance Computing resources (funded by NIH grant No. 1S10OD012300), supported by BioFrontiers IT group.