#  $\Omega$-bounds for the partial sums of some modified Dirichlet characters

Marco Aymone

###### Abstract.

We consider the problem of $\Omega$ bounds for the partial sums of a modified
character, i.e., a completely multiplicative function $f$, such that
$f(p)=\chi(p)$ for all but a finite number of primes $p$, where $\chi$ is a
primitive Dirichlet character. We prove that in some special circumstances,
$\sum_{n\leq x}f(n)=\Omega((\log x)^{|S|})$, where $S$ is the set of primes
$p$ where $f(p)\neq\chi(p)$. This gives credence to a corrected version of a
conjecture of Klurman et al., Trans. Amer. Math. Soc., 374 (11), 2021,
7967–7990. We also compute the Riesz mean of order $k$ for large $k$ of a
modified character, and showed that this is connected with the diophantine
properties of the irrational numbers of the form $\log p/\log q$, for primes
$p$ and $q$.

##  1\. Introduction.

A quite simple statement about a completely multiplicative111A function
$f:\mathds{N}\to\mathds{C}$ is completely multiplicative if $f(nm)=f(n)f(m)$
for all positive integers $n$ and $m$, and hence, such functions are
determined by its values at primes. function, and that was proved only
recently, is that no matter how we choose the values $f(p)$ at primes $p$ on
the unit circle, we will always end up with a completely multiplicative
function whose partial sums are unbounded, i.e.,

$\limsup_{x\to\infty}\left|\sum_{n\leq x}f(n)\right|=\infty.$

This result was proved in 2015 by Tao [8] in the context of the Erdős
discrepancy problem. Dirichlet characters222Dirichlet characters of modulus
$q$, often denoted by $\chi$, are $q$-periodic completely multiplicative
functions such that $\chi(n)=0$ whenever $\gcd(n,q)>1$. play an important role
in the Erdős discrepancy Theory: The non-principal characters $\chi$ are
completely multiplicative and have bounded partial sums. But of course they
are not a counter-example to the result of Tao since they vanish at a finite
subset of primes. Although, they seem to be extremal in discrepancy theory due
to the following fact: Up to this date, the known example of a completely
multiplicative $f$ with $|f|=1$ whose partial sums have the lowest known
fluctuation is obtained by adjusting a non-principal character $\chi$ at the
the primes $p$ where $\chi(p)=0$. For instance (see Borwein - Choi - Coons
[2]) we can take the non-principal character mod $3$, say $\chi_{3}$, and
define $f(3)=\pm 1$ and $f(p)=\chi_{3}(p)$ for all the other primes. The
partial sums up to $x$ of this modification of $\chi_{3}$ are $\ll\log x$.

Going further, one can easily show (see [1] for a proof using a Tauberian
result) that for a modified character, i.e., a function $f$ that is completely
multiplicative and that for some Dirichlet character $\chi$, $f(p)\neq\chi(p)$
only for $p$ in a finite subset of primes $S$, where for each of these primes
$p\in S$, $|f(p)|\leq 1$, we have

(1) $\sum_{n\leq x}f(n)\ll(\log x)^{|S|}.$

A very open question concerns $\Omega$ bounds for the partial sums of such
modifications.

The first treatment to such $\Omega$ bound is due to Borwein - Choi - Coons
[2] where they considered the case $|S|=1$ and showed that for a quadratic
Dirichlet character $\chi$ with prime modulus $q$, if we set $f(p)=\chi(p)$
for all primes except at $q$, where $f(q)=1$, then the partial sums

$\sum_{n\leq x}f(n)=\Omega(\log x).$

Later, Klurman et al. [6], among other results, proved a stronger form
Chudakov’s conjecture [5], and also conjectured that the partial sums of a
modified character $f$ are actually $\Omega((\log x)^{|S|})$. Their Corollary
1.6 states that the partial sums are $\Omega(\log x)$ if the set of
modifications $S$ is finite and for at least one prime $p\in S$, $|f(p)|=1$.
Their proof is correct when the Dirichlet character is primitive, but this
statement is not true for all characters. This is because, if we modify a non-
primitive character $\chi$ at a prime $r$ where $\chi(r)=0$, then there are
circumstances where the modification $f$ might end up in another character
with smaller modulus. Therefore, we restate their conjecture in the following
form.

######  Conjecture 1.1.

Let $\chi$ be a primitive Dirichlet character, and let $f$ be a completely
multiplicative modification of this character such that $f(p)=\chi(p)$ for all
but a finite subset $S$ of primes $p$, where for each prime $p\in S$,
$|f(p)|=1$. Then

$\sum_{n\leq x}f(n)=\Omega((\log x)^{|S|}).$

The first result of this paper gives credence to this conjecture by showing
that it is true in some very special circumstances.

######  Corollary 1.1.

Let $\chi$ be a primitive Dirichlet character with $\chi(-1)=-1$, and let $f$
be a completely multiplicative modification of this character such that
$f(p)=\chi(p)$ for all but a finite subset $S$ of primes $p$, where for each
prime $p\in S$, $f(p)=+1$. Then

$\sum_{n\leq x}f(n)=\Omega((\log x)^{|S|}).$

The result above is a direct consequence of the slightly more general Theorem
below.

######  Theorem 1.1.

Let $\chi$ be a primitive Dirichlet character, and let $f$ be a completely
multiplicative modification of this character such that $f(p)=\chi(p)$ for all
but a finite subset $S$ of primes $p$, where for each prime $p\in S$,
$|f(p)|=1$. Let

$T=\sum_{\begin{subarray}{c}p\in S\\\
f(p)=1\end{subarray}}1-\sum_{\begin{subarray}{c}p\in S\\\
\chi(p)=1\end{subarray}}1.$

Let

$N=\begin{cases}\max\\{0,T\\}$, \mbox{ if } $\chi(-1)=-1,\\\
\max\\{0,T-1\\},\mbox{ if }\chi(-1)=1.\end{cases}$

Then

$\sum_{n\leq x}f(n)=\Omega((\log x)^{N}).$

The main reason behind the result above is that the Dirichlet series of $f$,
say $F(s)$, can be written as

(2) $F(s)=\left(\prod_{p\in
S}\frac{1-\frac{\chi(p)}{p^{s}}}{1-\frac{f(p)}{p^{s}}}\right)L(s,\chi),$

and hence has a pole of order $N$ at $s=0$ that gives the claimed $\Omega$
bound.

###  1.1. Riesz means and irrationality of the numbers $\log p/\log q$

Another way to measure the mean behaviour of a sequence $(f(n))_{n}$ is by its
Riesz means. The Riesz mean of order $1$ is defined as

$R_{1}(x):=\sum_{n\leq x}f(n)\log(x/n).$

After partial summation is not difficulty to see that $R_{1}(x)$ is equal to
$\log x$ times the logarithmic average of $\sum_{n\leq x}f(n)$. Therefore, it
may regarded as a smooth average of $f(n)$.

Going further, we can define (see the book of Montgomery and Vaughan [7]) the
Riesz mean of order $k$:

$R_{k}(x):=\frac{1}{k!}\sum_{n\leq x}f(n)(\log(x/n))^{k}.$

It is not difficult to see that $\Omega$ results for $R_{k}(x)$ transfer to
the classical partial sums, as upper bounds for the classical partial sums
transfer to $R_{k}(x)$.

In our next result we were able to compute with precision the Riesz mean of
order $k$ of a modified character, which allow us to say that the exponent $N$
in Theorem 1.1 is optimal in the Riesz mean context.

######  Theorem 1.2.

Let $\chi$ be a primitive Dirichlet character and $f$ be a completely
multiplicative modification of it such that $f(p)=\chi(p)$ for all but a
finite subset $S$ of primes $p$, where for each prime $p\in S$, $f(p)$ is a
root of unity. Let $N$ be as in Theorem 1.1. Then there exists a constant
$\gamma>0$ such that for all integers $k\geq 10+\gamma(|S|+1)\max_{p\in
S}(\log p)^{2}$, as $x\to\infty$

$\sum_{n\leq x}f(n)(\log x/n)^{k}=P(\log x)+O(1),$

where $P(x)$ is a polynomial of degree $N+k$ with leading coefficient
$a_{N+k}$ given by

$a_{N+k}=c_{\chi}\frac{k!}{(N+k)!}\left(\prod_{\begin{subarray}{c}p\in S\\\
f(p)=1\end{subarray}}\frac{1-\chi(p)}{\log
p}\right)\left(\prod_{\begin{subarray}{c}p\in S\\\
\chi(p)=1\end{subarray}}\frac{\log
p}{1-f(p)}\right)\left(\prod_{\begin{subarray}{c}p\in S\\\ f(p),\,\chi(p)\neq
1\end{subarray}}\frac{1-\chi(p)}{1-f(p)}\right),$

where

$c_{\chi}=\begin{cases}L(0,\chi),\mbox{ if }\chi(-1)=-1,\\\
L^{\prime}(0,\chi),\mbox{ if }\chi(-1)=1.\end{cases}$

It is not difficult to see that if $\sum_{n\leq x}a_{n}\sim(\log x)^{N}$, then
the Riesz mean of order $k$ of this sequence is a polynomial $P(\log x)$ of
degree $N+k$. Therefore, the result above says that Theorem 1.1 is optimal in
the Riesz mean context.

From an analytic point of view, the Riesz mean behaves more and less as the
Césaro mean, and in the context of modified characters, Duke and Nguyen [4]
considered the content of Theorem 1.2 in the case that $|S|=1$. Their results
allowed them to state that there exists a $\pm 1$ sequence with bounded Césaro
“discrepancy”, in contrast with the infinity discrepancy of Tao [8]. Our
result above treats the case $|S|>1$. The new problem that didn’t appeared in
the case $|S|=1$, is that each Euler factor in (2) produces a periodic number
of simple poles at the line $Re(s)=0$ whose contributions to a Perron integral
of the corresponding Riesz mean is roughly at most

$\sum_{n=1}^{\infty}\frac{1}{n^{k}}\prod_{p\in
S\setminus\\{q\\}}\left|1-\exp\left(2\pi in\frac{\log p}{\log
q}\right)\right|^{-1}.$

The convergence of this sum is, therefore, connected to the problem of how
well the numbers $\log p/\log q$ can be approximated by rational numbers.

A parameter of irrationality of an irrational number $\alpha$ is the exponent
$\mu(\alpha)$ defined as the infimum over the real numbers $\eta$ such that
the inequality

$\left|\frac{a}{b}-\alpha\right|\leq\frac{C(\eta)}{b^{\eta}}$

admits only a finite number of rational solutions $a/b$ with $b\geq 1$, where
$C(\eta)$ is a constant that depends only in $\eta$ and the irrational number
$\alpha$. A classical result of Dirichlet is that $\mu(\alpha)\geq 2$ for all
irrational numbers $\alpha$, and this inequality is optimal for almost all
irrational numbers, with respect to the Lebesgue measure.

In our case, the irrationality of the numbers $\log p/\log q$ can be treated
by Baker’s theory of linear forms in logarithms, and we actually have by
Theorem 1.1 of Bugeaud [3] (see also the references therein) the following
Lemma.

######  Lemma 1.1.

Let $p$ and $q$ be distinct prime numbers. Then there exists a constant
$\gamma>0$ that does not depend on $p$ and $q$ such that, for $\mu=\gamma(\log
p)(\log q)$, there is a constant $C=C(p,q)$ such that the inequality

$\left|\frac{a}{b}-\frac{\log p}{\log q}\right|<\frac{C}{b^{\mu+1}}$

is satisfied only for a finite number of rational numbers $a/b$, with $b\geq
1$.

###  1.2. Structure of the paper

In some instances we assume that the reader is familiar with tools from
Analytic Number Theory. In section 2 we state the three main notations used in
this paper, then we quickly proceed with proofs in section 3. We then conclude
with some simulations at the end.

##  2\. Notation

We use the standard notation

  1. (1)

$f(x)\ll g(x)$ or equivalently $f(x)=O(g(x))$;

  2. (2)

$f(x)=o(g(x))$;

  3. (3)

$f(x)=\Omega(g(x))$;

1) is used whenever there exists a constant $C>0$ such that
$|f(x)|\leq|g(x)|$, for all $x$ in a set of numbers. This set of numbers when
not specified is the real interval $[L,\infty]$, for some $L>0$, but also
there are instances where this set can accumulate at the right or at the left
of a given real number, or at complex number. Sometimes we also employ the
notation $\ll_{\epsilon}$ or $O_{\epsilon}$ to indicate that the implied
constant may depends in $\epsilon$.

In case 2), we mean that $\lim_{x}f(x)/g(x)=0$. When not specified, this limit
is as $x\to\infty$ but also can be as $x$ approaches any complex number in a
specific direction.

In the last case 3), we say that $\limsup_{x}|f(x)|/|g(x)|>0$, where the limit
can be taken as in the case 2).

##  3\. proofs

###  3.1. Proof of Theorem 1.1

###### Proof of Theorem 1.1.

Let $f$, $\chi$ and $S$ be as in the statement of Theorem 1.1. Then, the
Dirichlet series of $f$ can be written as follows:

$\displaystyle F(s):=\sum_{n=1}^{\infty}\frac{f(n)}{n^{s}}$
$\displaystyle=\prod_{p\in S}\frac{1}{1-\frac{f(p)}{p^{s}}}\prod_{p\notin
S}\frac{1}{1-\frac{\chi(p)}{p^{s}}}$ $\displaystyle=\left(\prod_{p\in
S}\frac{1-\frac{\chi(p)}{p^{s}}}{1-\frac{f(p)}{p^{s}}}\right)L(s,\chi).$

Now observe that for the primes $p\in S$ such that $f(p)=+1$, we have that
$\chi(p)\neq 1$, and hence, the Euler factor corresponding to this prime in
the Euler product above contributes with a simple pole at $s=0$. Similarly,
for $p\in S$ with $\chi(p)=1$, we have that $f(p)\neq 1$, and hence, the Euler
factor corresponding to this prime in the Euler product above contributes with
simple zero at $s=0$. In the other cases, the Euler factor is a regular and
non-vanishing function at $s=0$. Therefore, as $s\to 0$, we have that

$\prod_{p\in
S}\frac{1-\frac{\chi(p)}{p^{s}}}{1-\frac{f(p)}{p^{s}}}\gg\frac{1}{|s|^{T}},$

where $T$ is as in Theorem 1.1.

Now, recall that we define $N=\max\\{0,T\\}$ if $\chi(-1)=-1$, otherwise, if
$\chi(-1)=1$, $N=\max\\{0,T-1\\}$. If $N=0$, then there is nothing to prove
since the statement of Theorem 1.1 becomes trivial. We therefore assume $N\geq
1$. We begin by recalling the functional equation for $L(s,\chi)$ (see for
instance the book [7]):

(3)
$L(s,\chi)=\varepsilon(\chi)L(1-s,\overline{\chi})2^{s}\pi^{s-1}q^{1/2-s}\Gamma(1-s)\sin\left(\frac{\pi}{2}(s+\kappa)\right),$

where: $|\varepsilon(\chi)|=1$, $\Gamma$ is the classical gamma function and
$\kappa=0$ if $\chi(-1)=1$, $\kappa=1$ if $\chi(-1)=-1$.

By this functional equation, and the fact that $L(1,\chi)\neq 0$ for all
$\chi$, we immediately see that $L(0,\chi)\neq 0$ in the situation that
$\chi(-1)=-1$, and, otherwise, has a simple zero coming from the sine function
at $s=0$ in the situation that $\chi(-1)=1$.

Therefore, as $s\to 0$, we have that $|F(s)|\gg\frac{1}{|s|^{T}}$ in the
situation that $\chi(-1)=-1$, and $|F(s)|\gg\frac{1}{|s|^{T-1}}$ in the
situation that $\chi(-1)=+1$.

Now, in order to complete the proof we will require the following Lemma:

######  Lemma 3.1.

Let $\alpha\geq 0$. Then, there exists a constant that depends only in
$\alpha$, say $c(\alpha)$, such that for all $\sigma>0$:

$I(\sigma,\alpha):=\int_{1}^{\infty}\frac{(\log
x)^{\alpha}}{x^{1+\sigma}}=\frac{c(\alpha)}{\sigma^{1+\alpha}}.$

We continue with the proof of Theorem 1.1 and postpone the proof of the Lemma
above to the end.

We recall that $F(s)$ can be written as the following Mellin transform of
$\sum_{n\leq x}f(n)$ valid for $Re(s)>1$:

$F(s)=s\int_{1}^{\infty}\left(\sum_{n\leq x}f(n)\right)\frac{dx}{x^{1+s}}.$

By the bound (1), this formula actually holds for all $Re(s)>0$. Therefore

$\frac{|F(\sigma)|}{\sigma}\leq\int_{1}^{\infty}\left|\sum_{n\leq
x}f(n)\right|\frac{dx}{x^{1+\sigma}}.$

Now, if $\chi(-1)=-1$, we have that

$\frac{1}{\sigma^{N+1}}\ll\int_{1}^{\infty}\left|\sum_{n\leq
x}f(n)\right|\frac{dx}{x^{1+\sigma}},$

and hence, by Lemma 3.1, the partial sums $\sum_{n\leq x}f(n)$ cannot be
$o((\log x)^{N})$, otherwise the integral would be $o(1/\sigma^{N+1})$.
Similarly, we obtain same conclusions in the case that $\chi(-1)=+1$, and this
completes the proof of the Theorem. ∎

###### Proof of Lemma 3.1.

We begin by making the change $u=\log x$. Then

$I(\sigma,\alpha):=\int_{0}^{\infty}\frac{u^{\alpha}}{e^{u\sigma}}du.$

Now we make another substitution: $v=\sigma u$. This leads to

$I(\sigma,\alpha):=\frac{1}{\sigma^{1+\alpha}}\int_{0}^{\infty}\frac{v^{\alpha}}{e^{v}}dv,$

and this completes the proof. ∎

###  3.2. Proof of Theorem 1.2

###### Proof of Theorem 1.2.

We have that (see the book [7], pg. 143)

$\sum_{n\leq x}f(n)(\log x/n)^{k}=\frac{k!}{2\pi
i}\int_{2-i\infty}^{2+i\infty}\frac{F(s)x^{s}}{s^{k+1}}ds.$

Now, for $X>0$

$\int_{2-i\infty}^{2+i\infty}\frac{F(s)x^{s}}{s^{k+1}}ds=\int_{2-iX}^{2+iX}\frac{F(s)x^{s}}{s^{k+1}}ds+O(x^{2}/X^{k}).$

Let $R_{X}$ be the rectangle with vertices $2-iX,2+iX,-1+iX$ and $-1-iX$. By
formula (2), we see that $F$ has a meromorphic continuation to $\mathds{C}$
with poles only at the line $Re(s)=0$. Actually, has a pole of order $N$ at
$s=0$, and simple poles of the form $2i\pi n/\log p$ where
$n\in\mathbb{Z}\setminus\\{0\\}$ and $p\in S$. Therefore, by the Cauchy
residue Theorem, $1/2\pi i$ times the integral along $\partial R_{X}$,
provided that $F(s)$ is regular at $s=\pm iX$, is equal to

$Res\left(\frac{F(s)x^{s}}{s^{k+1}},s=0\right)+\sum_{p\in
S}\sum_{\begin{subarray}{c}-\frac{X\log p}{2\pi}\leq n\leq\frac{X\log
p}{2\pi}\\\ n\neq
0\end{subarray}}Res\left(\frac{F(s)x^{s}}{s^{k+1}},s=\frac{2\pi in}{\log
p}\right).$

Clearly, $Res\left(\frac{F(s)x^{s}}{s^{k+1}},s=0\right)$ is the claimed
polynomial $P(\log x)$ of degree $N+k$. Further, the exact value of the
leading coefficient $a_{N+k}$ can be easily computed.

Now we are going to show that the contribution of the simple poles at
$Re(s)=0$ is at most $O(1)$ as $x\to\infty$ in the claimed range of $k$. By
formula (2), we see that by denoting $S=\\{p_{1},p_{2},...,p_{|S|}\\}$, then
for a fixed prime $p\in S$

$\displaystyle\sum_{\begin{subarray}{c}-\frac{X\log p}{2\pi}\leq
n\leq\frac{X\log p}{2\pi}\\\ n\neq
0\end{subarray}}Res\left(\frac{F(s)x^{s}}{s^{k+1}},s=\frac{2\pi in}{\log
p}\right)$
$\displaystyle\ll\sum_{n=1}^{\infty}\frac{1}{n^{k+1}}\left|L\left(\frac{2\pi
in}{\log p},\chi\right)\right|\prod_{p_{k}\neq
p}\left|1-f(p_{k})\exp\left(2\pi in\frac{\log p_{k}}{\log
p}\right)\right|^{-1}.$

Since each $f(p_{k})$ is a root of unity, there are coprime positive integers
$a_{k}\leq b_{k}$ such that $f(p_{k})=\exp(2\pi ia_{k}/b_{k})$. Thus, by
letting $\mu$ as in Lemma 1.1, the product above is, except for a finite
number of $n$, $\ll n^{(\mu+1)|S|}$. Hence, by classical estimates for the
Gamma function and for $L(1+it)$ combined with the functional equation (3), we
have that the last series above is absolutely convergent when
$k\geq(\mu+1)|S|+10$.

To complete the argument, by combining the classical estimates for the Gamma
function with (2) and the functional equation (3), the integral

$\int_{-1-iX}^{-1+iX}\frac{F(s)x^{s}}{s^{k+1}}ds\ll\int_{-1-iX}^{-1+iX}\frac{|L(s,\chi)|x^{-1}}{|s|^{k+1}}ds\ll
1/x.$

Further, since the number of simple poles of $F(s)$ up to height $X$ is $\ll
X$, there is a fixed $\delta>0$ and an infinite number of points
$X_{j}\to\infty$ such that the distance between $\pm iX_{j}$ and each one of
these simple poles is at least $\delta/2$. Therefore, along the lines
$I_{\pm}=[-1\pm iX_{j},2\pm iX_{j}]$, we have that

$\int_{I_{\pm}}\frac{F(s)x^{s}}{s^{k+1}}ds\ll\frac{x^{2}}{X_{j}}.$

Then we obtain the claimed result by making $X_{j}\to\infty$. ∎

##  4\. Some simulations

Here we make some simulations with the non-principal Dirichlet character
$\chi$ with modulus $3$, i.e.,

$\chi(n)=\begin{cases}1,\mbox{ if }n\equiv 1\mod 3,\\\ -1,\mbox{ if }n\equiv
2\mod 3,\\\ 0,\mbox{ otherwise. }\end{cases}$

For the first primes we have:

$p$ | $2$ | $3$ | $5$ | $7$ | $11$ | $13$ | $17$ | $19$
---|---|---|---|---|---|---|---|---
$\chi(p)$ | $-1$ | $0$ | $-1$ | $1$ | $-1$ | $1$ | $-1$ | $1$

In Figure 1, we consider four modifications where we turn the values
$\chi(p)=-1$ to $+1$, for four primes $p$. In this case, if $N$ is as in
Theorem 1.1, we have that $N=4$. In Figure 2 we consider four modifications,
but $N=3$, and in Figure 3 we also consider four modifications, but $N=0$.

Figure 1. A case where $N=4$. In blue we plot the partial sums of $f$, where
$f$ is the modification of $\chi$ such that $f(2)=f(3)=f(5)=f(11)=+1$. Figure
2. A case where $N=3$ with four modifications. In blue we plot the partial
sums of $f$, where $f$ is the modification of $\chi$ such that
$f(2)=f(3)=f(5)=f(11)=+1$ and $f(7)=-1$. Figure 3. A case where $N=0$ with
four modifications. In blue we plot the partial sums of $f$, where $f$ is the
modification of $\chi$ such that $f(3)=f(7)=f(13)=f(10)=-1$.

Acknowledgements. I am warmly thankful to Carlos Gustavo Moreira (Gugu) for a
fruitful discussion on Diophantine Theory and for pointing out the reference
of Bugeaud [3], and for Oleksiy Klurman for his comments on a draft version of
this paper. This project was supported by CNPq grant Universal no.
403037/2021-2.

## References

  * [1] M. Aymone, A note on multiplicative functions resembling the Möbius function, J. Number Theory, 212 (2020), pp. 113–121.
  * [2] P. Borwein, S. K. K. Choi, and M. Coons, Completely multiplicative functions taking values in $\\{-1,1\\}$, Trans. Amer. Math. Soc., 362 (2010), pp. 6279–6291.
  * [3] Y. Bugeaud, Effective irrationality measures for quotients of logarithms of rational numbers, Hardy-Ramanujan J., 38 (2015), pp. 45–48.
  * [4] W. Duke and H. N. Nguyen, Riesz means of certain arithmetic functions, J. Number Theory, 210 (2020), pp. 132–141.
  * [5] O. Klurman and A. P. Mangerel, Rigidity theorems for multiplicative functions, Math. Ann., 372 (2018), pp. 651–697.
  * [6] O. Klurman, A. P. Mangerel, C. Pohoata, and J. Teräväinen, Multiplicative functions that are close to their mean, Trans. Amer. Math. Soc., 374 (2021), pp. 7967–7990.
  * [7] H. L. Montgomery and R. C. Vaughan, Multiplicative number theory. I. Classical theory, vol. 97 of Cambridge Studies in Advanced Mathematics, Cambridge University Press, Cambridge, 2007.
  * [8] T. Tao, The Erdős discrepancy problem, Discrete Anal., (2016), pp. Paper No. 1, 29.


# Computing with Residue Numbers in
High-Dimensional Representation

Christopher J. Kymn  Redwood Center for Theoretical Neuroscience, University
of California, Berkeley, CA Denis Kleyko  Centre for Applied Autonomous Sensor
Systems, Örebro University, Sweden  Intelligent Systems Lab, Research
Institutes of Sweden, Kista, Sweden  E. Paxon Frady  Neuromorphic Computing
Lab, Intel, Santa Clara, CA  Connor Bybee  Redwood Center for Theoretical
Neuroscience, University of California, Berkeley, CA Pentti Kanerva  Redwood
Center for Theoretical Neuroscience, University of California, Berkeley, CA
Friedrich T. Sommer  Redwood Center for Theoretical Neuroscience, University
of California, Berkeley, CA Neuromorphic Computing Lab, Intel, Santa Clara, CA
Bruno A. Olshausen  Redwood Center for Theoretical Neuroscience, University of
California, Berkeley, CA

(November 7, 2023)

###### Abstract

We introduce Residue Hyperdimensional Computing, a computing framework that
unifies residue number systems with an algebra defined over random, high-
dimensional vectors. We show how residue numbers can be represented as high-
dimensional vectors in a manner that allows algebraic operations to be
performed with component-wise, parallelizable operations on the vector
elements. The resulting framework, when combined with an efficient method for
factorizing high-dimensional vectors, can represent and operate on numerical
values over a large dynamic range using vastly fewer resources than previous
methods, and it exhibits impressive robustness to noise. We demonstrate the
potential for this framework to solve computationally difficult problems in
visual perception and combinatorial optimization, showing improvement over
baseline methods. More broadly, the framework provides a possible account for
the computational operations of grid cells in the brain, and it suggests new
machine learning architectures for representing and manipulating numerical
data.

##  1 Introduction

The problem of representing and computing on high-dimensional representations
of numerical values—such as position, velocity, and color—is central to both
machine learning and computational neuroscience. In machine learning, vector
representations of numbers are useful for defining position or function
encodings in neural networks [1, 2, 3], improving robustness to adversarial
examples [4], and generating efficient classifiers [5]. In neuroscience,
experimentalists seek to understand how populations of neurons in the brain
represent and transform perceptual or cognitive variables, and so numerous
theorists have constructed models for how these variables could be encoded in
and decoded from high-dimensional vector encodings (e.g., [6, 7, 8]).

A particularly salient example of high-dimensional representation in
neuroscience is the ‘grid cell’ encoding of spatial position in the medial
entorhinal cortex [9]. Grid cells have multiple peaks in their firing rates
that correlate with spatial positions arranged in a hexagonal lattice. While
somewhat perplexing at first glance, the usefulness of such a coding scheme
becomes apparent from how it functions as a population code. In comparison to
a population of neurons with traditional unimodal encoding functions whose
coding resolution increases linearly with the number of neurons, a grid cell
population possesses a coding resolution that grows exponentially in the
number of neurons [10]. In particular, Fiete and colleagues have emphasized
that this advantage of grid cell encoding utilizes properties of residue
numbers (see Section 2.2 below) [11].

Inspired by this observation, we propose a comprehensive algebraic framework
for distributed neural computation based on residue number systems. Our novel
framework builds on an existing algebraic framework for computing with high-
dimensional random vectors. The idea was originally called computing in
Holographic Reduced Representation [12] and is now referred to also as Vector
Symbolic Architectures (VSA) [13] and Hyperdimensional Computing [14]. We call
our new framework Residue Hyperdimensional Computing (RHC) and demonstrate
that it inherits the computational advantages of both standard residue number
systems and hyperdimensional computing. This enables fault-tolerant
computations that can efficiently represent numbers and search over a large
dynamic range with greatly reduced memory requirements. Furthermore, as we
shall see, the new framework provides a useful formalism for understanding
computations in grid cells.

To summarize, we list the four key coding properties we achieve with Residue
Hyperdimensional Computing: (1) Algebraic structure: Simple operations on
vectors perform addition and multiplication on encoded values, (2)
Expressivity: Feasible encoding range scales better than linearly with
dimension, (3) Efficient decoding: Required resources to decode scale better
than linearly with encoding range, and (4) Robustness to noise. Although a
number of previously proposed models achieve some of these properties (Table
S.1, Supplement A), Residue Hyperdimensional Computing is the first, to our
knowledge, to achieve all four of these desiderata, as we shall now show.

##  2 Results

We first define the key concepts on which the Residue Hyperdimensional
Computing framework is based (Section 2.1) and then describe the framework
fully in Section 2.2. We then demonstrate its favorable encoding, decoding,
and robustness properties (Section 2.3), as well as how it can be extended to
multiple dimensions (Section 2.4) and sub-integer encodings (Section 2.5). Of
particular note, we construct a hexagonal residue encoding system, analogous
to grid cell coordinates, that provides higher spatial resolution than square
lattices. Finally, we describe how the framework can be applied to problems in
visual scene analysis and combinatorial optimization (Section 2.6).

###  2.1 Preliminary definitions

Definition 2.1.1. A Residue Number System [15] encodes an integer
$x\in\mathbb{Z}$ by its value modulo $\\{m_{1},m_{2},\dots,m_{K}\\}$, where
the $m_{k}$ are the moduli of the system. For example, relative to moduli
$\\{3,5,7\\}$, $x=20$ would be encoded by the residue $[2,0,6]$—i.e.,
$[20\bmod 3,\;20\bmod 5,\;20\bmod 7]$. The Chinese Remainder Theorem states
that if the moduli are pairwise co-prime, then for any $x$ such that $0\leq
x<M\vcentcolon=\prod_{k}m_{k}$, the integer is uniquely encoded by its residue
[16]. From here on, we will assume that the pairwise co-prime condition is
fulfilled.

Definition 2.1.2. Fractional Power Encoding (FPE) [12] defines a randomized
mapping from an integer $x$ to a high-dimensional vector $\mathbf{z}(x)$. Let
$D$ be the dimension of the vector, with $D$ typically in the range
$10^{2}\leq D\leq 10^{4}$. The encoding scheme involves two steps. First, draw
a random base vector, $\mathbf{z}$, defined as follows:

$\mathbf{z}=[e^{i\phi_{1}},e^{i\phi_{2}},\dots,e^{i\phi_{D}}]$ (1)

where each element $e^{i\phi_{j}}$ is a complex number with unit magnitude (a
phasor), and each $\phi_{j}$ is a random sample from a specified probability
distribution. Second, define a function from $x$ to $\mathbb{C}^{D}$ via
component-wise exponentiation of the random vector:

$\mathbf{z}(x)=\mathbf{z}^{x}=[e^{i\phi_{1}x},e^{i\phi_{2}x},\dots,e^{i\phi_{D}x}]$
(2)

Definition 2.1.3. A kernel, $K(x_{1},x_{2})$, is a function
${\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}\crossproduct{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}\to\mathbb{R}$
that measures the similarity between two objects in a set $\textstyle\chi$
(e.g., vectors in $\mathbb{R}^{n}$). Notably, FPE implements kernel
approximation [17], which is widely used in machine learning [18]. More
specifically, we can induce a kernel based on the inner products of FPEs:

$K(x_{1},x_{2})=\frac{1}{D}\ \mathfrak{R}\\{\mathbf{z}(x_{1})^{T}\
\overline{\mathbf{z}(x_{2})}\\}$ (3)

where $\overline{\mathbf{z}(x_{2})}$ is the complex conjugate of
$\mathbf{z}(x_{2})$. This defines a translation-invariant kernel $K(\Delta x)$
(where $\Delta x=x_{1}-x_{2}$), which converges to a particular $K^{*}(\Delta
x)$ as $D\to\infty$, where the shape of the kernel is determined by the
probability distribution used to draw $\mathbf{z}$ [19, 20].

###  2.2 Residue Hyperdimensional Computing

We now introduce how FPE can implement a residue number system. As a first
step, we explain how FPE can implement congruence (representing a remainder,
modulo $m$).

Definition 2.2.1. Fractional Power Encoding, modulo $m$: For $\mathbf{z}_{m}$,
let the support of our probability distribution for the $\phi_{j}$ be the
$m$-th roots of unity. In other words, each $\phi_{j}$ must be a multiple of
$2\pi/m$. Then congruent values are mapped to the same vector:

$\displaystyle\mathbf{z}_{m}(x+m)$
$\displaystyle=[e^{i\phi_{1}(x+m)},e^{i\phi_{2}(x+m)},\dots,e^{i\phi_{D}(x+m)}]$
$\displaystyle=[e^{i\phi_{1}x}\cdot e^{i\phi_{1}m},e^{i\phi_{2}x}\cdot
e^{i\phi_{2}m},\dots,e^{i\phi_{D}x}\cdot e^{i\phi_{D}m}]$
$\displaystyle=[e^{i\phi_{1}x},e^{i\phi_{2}x},\dots,e^{i\phi_{D}x}]$
$\displaystyle=\mathbf{z}_{m}(x)$

because $e^{i\phi_{j}m}=e^{2\pi\cdot i\cdot k}$ for some integer $k$, and
$e^{2\pi\cdot i\cdot k}=1$ for any integer $k$. Put another way,
$\mathbf{z}_{m}(x)$ is a representation (in the abstract algebraic sense) of
the additive group of integers modulo $m$.

The kernel induced by $\mathbf{z}_{m}(\Delta x)$ is $1$ if $\Delta x=0$ (mod
$m$), and $\approx 0$ otherwise, as shown in Figure 1a. This is highly useful,
because it implies that distinct integers behave as quasi-orthogonal vectors,
just like symbols in hyperdimensional computing. Unlike symbols, however, we
can perform algebraic manipulations transforming one integer into another.

Definition 2.2.2. Residue Hyperdimensional Computing: Let
$\mathbf{z}_{m_{1}}$, $\mathbf{z}_{m_{2}}$, $\dots$, $\mathbf{z}_{m_{K}}$
denote FPE vectors with moduli $m_{1},m_{2},\dots,m_{K}$ respectively. Let
$\odot$ denote component-wise multiplication (a.k.a. Hadamard product). Then,
we encode an integer $x$ by combining our modulo representations via the
Hadamard product:

$\mathbf{z}(x)=\bigodot_{k=1}^{K}\mathbf{z}_{m_{k}}(x)$ (4)

The above encoding represents the remainder of $x$, because each
$\mathbf{z}_{m_{k}}$ represents its value modulo $m_{k}$. The code is fully
distributed, as every element of the vector $\mathbf{z}$ contains information
about each encoding vector $\mathbf{z}_{m_{k}}(x)$. By contrast, typical
implementations of residue number systems compartmentalize information about
individual moduli [21].

The kernel induced by $\mathbf{z}(\Delta x)$ is $1$ if $\Delta x=0$ (mod
$\prod_{k}m_{k}=M$), and $\approx 0$ for other integer intervals $\Delta x$
(Figure 1b). This means that the kernel maps different remainders of our
residue number system to quasi-orthogonal directions in high-dimensional
vector space, and enables computing in superposition over these variables.
Examples of possible applications enabled by such a scheme are presented in
Section 2.6.

The hallmark of a residue number system is carry-free arithmetic; that is,
addition, subtraction, and multiplication can be performed component-wise on
the remainders. This enables residue number systems to be highly parallel,
avoiding carryover operations required in binary number systems. Residue
Hyperdimensional Computing implements arithmetic with component-wise
operations, thus inheriting the premier computational property of residue
number systems.

Addition is defined as the Hadamard product between vectors, that is
$\mathbf{z}(x_{1}+x_{2})=\mathbf{z}(x_{1})\odot\mathbf{z}(x_{2})$ (Methods
4.1.1). This follows from the fact that component-wise multiplication of
phasors corresponds to phase addition, and that component-wise multiplication
is commutative. Subtraction is defined by addition of the additive inverse.

Next, we define a second binding operation that implements multiplication,
denoted as $\star$: $\mathbf{z}(x_{1}\cdot
x_{2})=\mathbf{z}(x_{1})\star\mathbf{z}(x_{2})$. Just as variable addition is
implemented by element-wise multiplication, variable multiplication is
implemented by another element-wise operation, this one involving
exponentiation (Methods 4.1.2). Here, it is crucial that our encoding function
is restricted to integers as its domain, because multiplication is commutative
and integer powers commute (that is,
$(c^{x_{1}})^{x_{2}}=c^{x_{1}x_{2}}=(c^{x_{2}})^{x_{1}}$ for $c\in\mathbb{C}$
and integers $x_{1},x_{2}$). We show how this definition for integer
multiplication can be generalized to multiplication for vector encodings
$\mathbf{z}(x_{1})$ and $\mathbf{z}(x_{2})$, without invoking the costly step
of first decoding the integers back from the vectors.

Division is not well-defined for residue number systems, because integers are
not closed under division. Still, when a modulus $m_{k}$ is prime,
multiplications by non-zero integers are invertible, because each non-zero
integer has a unique multiplicative inverse with respect to $m_{k}$.

The existence of two distinct binding operators for addition and
multiplication is a new contribution to hyperdimensional computing. Previous
formulations only supported addition, or only multiplication via addition
after taking logarithmic transformations [22]. Consequently, we can now
formulate a fully distributed residue number system that inherits the benefits
of computing with high-dimensional vectors.

Figure 1: Residue Hyperdimensional Computing defines a kernel separating
different remainder values, and it enables algebraic operations. a, For
fractional power encoding, modulo $m=5$, inner products between vectors
reflect the similarity of points with the same remainder value and are quasi-
orthogonal elsewhere. The light blue curve shows the kernel shape when $\Delta
x$ is a real-valued scalar; integers occur approximately at zero crossings. A
further derivation is provided in Supplement B. b, The kernel induced by an
RHC vector (green) is the product kernel of the moduli used to form it
(orange, blue). c, Demonstration of addition and multiplication. Blue and
orange show encodings of 2 and 3, respectively. Teal shows the decoded value
of 2+3 (i.e., 5); purple shows decoded value of $2\times 3$ (i.e., 6).

###  2.3 The resonator network enables efficient decoding of residue
numbers

Given the vector representation of a residue number, $\mathbf{z}(x)$, how do
we decode it to recover its value, $x$? One method for decoding commonly used
in hyperdimensional computing is codebook decoding [23], which involves taking
the inner product of $\mathbf{z}(x)$ with a set of $M$ codebook vectors with
known values of $x$. However, this procedure requires $\mathcal{O}(M*D)$
storage and $M$ inner product evaluations.

Fortunately, we can improve the situation by utilizing the fact that residue
numbers break the dynamic range of a number into a set of smaller numbers,
each with lower overall dynamic range. For example, a number with dynamic
range of 105, when represented modulo $[3,\ 5,\ 7]$, consumes a total dynamic
range of $3+5+7=15$. This in turn reduces both the memory and computation
resources for decoding by a factor of $105/15=7$. To make this work though
requires that we invert Equation 4—that is, we must factorize $\mathbf{z}(x)$
into the set of constituent vectors
$\\{\mathbf{z}_{m_{1}}(x),\mathbf{z}_{m_{2}}(x),\dots,\mathbf{z}_{m_{K}}(x)\\}$
representing $x$ modulo $m_{k}$, from which $x$ can be easily recovered. For
this we can use a resonator network [24, 25], a recently discovered method for
efficiently factorizing vectors in hyperdimensional computing. Figure 2a shows
that for a range of $M$ values, the resonator network can recover vectors over
an order of magnitude faster than standard codebook decoding. Two parameters
that contribute to this are the vector dimension ($D$) and number of moduli
($K$).

To evaluate the dependence of resonator decoding on vector dimension, we fix
the number of moduli ($K=2$) and calculate the empirical accuracy (Figure 2b)
of the resonator on the effective range $M$. We find that for a fixed $D$, the
accuracy remains almost perfect up to a certain range of $M$, after which
accuracy rapidly decays. To evaluate scaling with $D$, we define the capacity
$C$ of a dimension to be the largest $M$ up to which empirical accuracy is at
least 95 percent. We find that the scaling of $C(D)$ is well-fit with a
quadratic polynomial (Figure 2c), consistent with previous scaling laws
studied for a resonator network with two states per component [25]. Further
tests with higher dimensions would help confirm quadratic scaling, but even
the linear scaling has high slope (and note that $C(4096)>2\times 10^{6}$).

To evaluate dependence on the number of moduli, we fix $D=1024$ and vary $K$.
We find that resonator capacity decreases as $K$ increases (Figure 2d), also
consistent with prior work [25]. Still, we emphasize that resonators with
higher $K$ have two advantages: decreased computation per decoding (Figure
2a), and decreased memory requirements. The resonator requires only
$\sum_{k}m_{k}=b$ codebook vectors, rather than $\prod_{k}m_{k}=M$. This means
that increasing $K$ can increase the effective range by several orders of
magnitude given a fixed codebook budget (Figure 2e). Remarkably, the maximal
$M$ for a given $b$ is given by Landau’s function $g(b)$, which scales as
$g(b)=e^{(1+o(1))\sqrt{b\ln b}}$ [26]. This implies an exponential scaling
between storage required and effective range, if we proffer sufficient $K$ to
achieve it.

Finally, we evaluate how robust resonator network decoding is to noise. We
draw phase noise from a von Mises distribution with mean 0 and concentration
$\kappa$; higher $\kappa$ indicates less noise. In Figure 2f, we observe that
performance degrades gradually as a function of noise, yet capacity remains
remarkably high even at high noise levels.

Figure 2: The resonator network can efficiently recover moduli of different
numbers. a, The resonator network outperforms codebook decoding by over an
order of magnitude when the effective range is within resonator network
capacity. b, Demonstration of resonator network capacity. For a fixed vector
dimension, accuracy remains high up to a given range, before gradually falling
off. c, Scaling of resonator network capacity, $C$, as a function of
dimension, $D$, is well-described by a quadratic fit (orange dashed line; cf.
linear fit in blue dash-dotted line). Quadratic fit coefficients are
$(1.3\times 10^{-1})D^{2}+(2.4\times 10^{1})D-3.7\times 10^{3}$, linear fit
coefficients are $(7.6\times 10^{2})D-7.6\times 10^{5}.$ d, Resonator network
performance is slightly worse for a higher number of moduli, $K$, but e, an
advantage of lower moduli is a much higher effective range given encoding
resources. f, The capacity of the resonator network remains high even in the
presence of large amounts of phase noise.

###  2.4 Generalization to multiple dimensions

####  2.4.1 Cartesian representations of $\mathbf{Z}^{n}$

Next, we generalize Residue Hyperdimensional Computing from scalars to multi-
dimensional variables, showing that the core operations and principles still
apply. Let $\mathbf{x}\in\mathbb{Z}^{n}$ be a low-dimensional vector. Let
$x_{1},x_{2},\dots,x_{n}$ denote the components of $\mathbf{x}$. To encode
$\mathbf{x}$ with a single high-dimensional vector
$\mathbf{z}\in\mathbb{V}^{D}$ ($D>>n)$ we form the encoding $\mathbf{z}(x)$ by
taking the Hadamard product of the distributed representations of individual
components:

$\mathbf{z}(\mathbf{x})=\mathbf{z}_{1}(x_{1})\odot\mathbf{z}_{2}(x_{2})\odot\dots\odot\mathbf{z}_{n}(x_{n})$
(5)

where each $\mathbf{z}_{i}$ is a random vector as generated for the residue
representation of a single number in Equation (4). (Each $\mathbf{z}_{i}$ is
created from binding multiple vectors for each of the moduli, distinct for
each dimension $i$.) Since the binding operation is commutative, we can
rearrange our terms to show the following useful properties hold:

  1. (a)

The Hadamard product ($\odot$) performs vector addition:
$\mathbf{z}(\mathbf{x})\odot\mathbf{z}(\mathbf{x^{\prime}})=\mathbf{z}(\mathbf{x+x^{\prime}})$.

  2. (b)

The multiplicative binding operation ($\star$) performs component-wise
multiplication of $\mathbf{x}$ and $\mathbf{x^{\prime}}$:
$\mathbf{z}(\mathbf{x})\star\mathbf{z}(\mathbf{x^{\prime}})=\mathbf{z}(\mathbf{x}\odot\mathbf{x^{\prime}})$.

  3. (c)

The kernel induced by $\mathbf{z}$ is the product of the kernels of the
individual components: $K_{z}(\mathbf{\Delta
x})=\prod_{i=1}^{n}K_{z_{i}}(\Delta x_{i})$, where $\mathbf{\Delta
x}=\mathbf{x}-\mathbf{x^{\prime}}$.

While (a) and (c) are general properties of FPE, (b) is once again unique to
Residue Hyperdimensional Computing. In addition, the savings in decoding
resources and computation required (Figure 2) also scale in higher dimensions,
as utilized in Section 2.6.

####  2.4.2 Hexagonal coordinate systems

When working in a multi-dimensional space, there are multiple alternatives to
a Cartesian coordinate system. For example, grid cells in medial entorhinal
cortex encode spatial location with a hexagonal coordinate system; research in
theoretical neuroscience suggests that this is because such a tiling of space
has the highest resolution (Fisher information) in 2D space [27]. As an
illustrative example, we show that Residue Hyperdimensional Computing can also
implement hexagonal coordinate systems, and that such a hexagonal lattice
retains coding advantages over square lattices. Of particular note, we
formulate a self-consistent encoding that extends a residue number system to
non-negative hexagonal coordinates.

Figure 3: Definition of a residue number system in non-negative hexagonal
coordinates. a and c show discrete phase distributions chosen for a 3D
coordinate system of a 2D space with moduli 3 and 5, respectively. In addition
to requiring that phases are drawn from the $m$-th roots of unity, we enforce
that the three phases sum to $0\ (\text{mod}\ 2\pi)$. b and d show the
respective kernels generated by these phase distributions. e, An example of
the Voronoi tesselation of different states composed of a hexagonal coordinate
system with modulus 5. Each color corresponds to a different state
representation in the vector space of integers. f, The kernel induced by a
hexagonal residue HD vector (period of $3\cdot 5=15$). g, Compared to square
encodings of space, hexagonal encodings approximately triple the effective
range of encodable states with only a 50 percent increase in required storage
space. h, The Shannon entropy of the hexagonal code is higher than that of the
square code with the same modulus, reflecting the benefit of hexagonal
packing.

To encode a two-dimensional position into a three-coordinate frame requires
two steps. First, we project the 2D vector $\mathbf{x}$ into a 3D vector
$\mathbf{y}$ with unit vectors whose angles each differ by $\frac{2\pi}{3}$
(Methods 4.3). This coordinate representation is known as the ‘Mercedes-Benz’
frame in $\mathbb{R}^{2}$; it is well-studied in signal processing [28] and
has attracted recent interest in hyperdimensional computing (e.g., [29, 19]).
For our purposes, this step is necessary, but not sufficient, because
projections can result in negative values. Second, to rectify values, we
encode $\mathbf{y}$ with the method described in Section 2.4.1, but with the
additional constraint that $\mathbf{z}([1,1,1])=\mathbf{z}([0,0,0])$, ensuring
that every state with negative coordinates has an equivalent representation to
one with non-negative coordinates (Figures 3a and 3c, Methods 4.3). This
constraint also reflects the fact that equal movement in every direction
cancels out, and thus it enforces that different paths to the same 2D position
result in the same high-dimensional encoding. The kernels induced by vectors
of individual moduli (Figures 3b and 3d) and by the residue vector (Figure 3f)
exhibit the six-fold symmetry characteristic of hexagonal lattices and grid
cells [9].

We can, therefore, represent the hexagonal coordinate system with a Voronoi
tesselation (Figure 3e) in which different regions of space are mapped to
their nearest integer-valued 3D coordinate. A hexagonal system with modulus
$m$ has $3m^{2}-3m+1$ distinct states and requires $3m$ codebook vectors,
whereas a square lattice has $m^{2}$ distinct states and requires $2m$
codebook vectors (Figure 3g). Thus, the hexagonal system achieves better
spatial resolution (it is a higher entropy code with regards to space) than a
square lattice (Figure 3h) for the same number of resources.

###  2.5 Extensions to sub-integer decoding resolution

Figure 4: The resonator network enables retrieval of fractional (sub-integer)
values. a, Example of the resonator network converging for $\mathbf{z}(x)$,
$x=40.4$. b, Inner product values decoded by the resonator network are
predicted by fractional offsets from a Dirac comb convolved with a sinc
function. c, Sub-integer encoding accuracy under a low noise regime for a RHC
vector ($r$ denotes number of sub-integer partitions). d, Bits per vector for
different fractional decodings. e, f, same as c and d, respectively, but under
higher noise conditions.

In previous sections, we worked exclusively with integer states and residue
number systems implementing them. Intriguingly, however, we can extend our
definition of FPE to rational numbers (Methods 4.4.1), and the resonator
converges to FPE encodings of non-integers, even when codebooks contain only
encodings of integers (Figures 4a and b). Strictly speaking, such extensions
beyond integers are not residue number systems, and multiplicative binding is
no longer well-defined. However, extensions towards sub-integer resolution
have been considered in theoretical analyses of grid cells, e.g., in [11, 30],
and we show that the resonator dynamics achieve this sub-integer resolution.

An efficient procedure for decoding with sub-integer precision is suggested by
Figure 4b. The inner products between codebook states and the final resonator
network state are well-described by evaluations of a Dirac comb convolved with
a sinc function (Supplement B). For integer encodings, this function would
evaluate to 0 for all features not near the peak, but for non-integers this is
no longer the case. Still, we can find the sub-integer offset that best
matches the resonator network state in order to decode the sub-integer value
(Methods 4.4.2).

Phase noise is the limiting factor in decoding sub-integer states. To quantify
this more rigorously, we evaluate a resonator network with varying effective
ranges $M$ under different noise regimes ($\kappa=16$ and $1$, respectively).
We then split each unit interval into $r$ partitions, so that there are
$M\cdot r$ distinct numbers represented. Figures 4c and 4e show the accuracy
of decoding for a different number of partitions with $\kappa=16$ and
$\kappa=1$, respectively. To account for accuracy and the number of different
states distinguished, we also report the ‘bits per vector’ metric [31] in
Figures 4d and 4f. This metric validates that with lower noise, we can more
reliably decode a higher number of states.

###  2.6 Applications

####  2.6.1 Efficient disentangling of object shape and pose from images

Figure 5: Residue Hyperdimensional Computing enables efficient disentangling
of images. a, Processing pipeline: an image is first represented in terms of
its local shape features via convolutional sparse coding, then converted to a
high-dimensional vector $\mathbf{z}$ by superimposing the residue number
encodings of the positions of each feature in the image, which is finally
factorized into object identity and position via a resonator network. b,
Simulations of the resonator network on visual scenes without the residue
number encoding (Standard) and with the residue number encoding (Residue). c,
With a residue number system, the resonator network requires less memory
overhead (40 vs. 220 codebook vectors) and less total computation to converge
to the correct solution (792.4 vs. 2085.6 average codebook evaluations). Error
bars show the 25th and 75th percentile of the number of evaluations.

Here we study the disentangling problem in vision—that is, the task of
recovering the underlying components of a scene given only the image pixel
values. Such problems abound in visual perception; examples include inferring
structure-from-motion or separating spectral reflectance properties from
shading due to 3D shape and lighting. How brains disentangle these factors of
variation in general is unknown, and it is computationally challenging due to
the combinatorial explosion in how factors combine to create any given scene
[32].

Here we demonstrate how Residue Hyperdimensional Computing can efficiently
tackle the combinatorial complexity inherent in visual perception by
considering a simple case of images containing three factors of variation:
object shape, horizontal position, and vertical position. Let $O,H,V$ be
finite sets listing the possible vectors for each factor. The goal is to infer
the $O_{a}\in O,H_{b}\in H,V_{c}\in V$ provided an image, $I$. In this setup,
the search space is $|O|\cdot|H|\cdot|V|$, and in our example (Figure 5),
$|O|=10$, and $|H|=|V|=105$, giving a search space of $\sim 10^{5}$.

We solve this problem in two stages. First, we form a latent, feature-based
representation of the image via convolutional sparse coding (Methods 4.5).
This step mirrors the neural representation in primary visual cortex, which is
hypothesized to describe image content in terms of a small number of image
features [33]. We observe that this step is useful as it helps to decorrelate
image patterns, thus achieving higher accuracy and faster convergence for the
resonator network.

Second, we encode the latent feature representation into a high-dimensional
vector that can be subsequently factorized into its components ($O_{a}$,
$H_{b}$, $V_{c}$) via a resonator network. This is accomplished, following
[34], by superimposing the residue number encodings of the position of each
image feature into a single scene vector, $\mathbf{s}$ (Methods 4.5, equation
10). The resulting vector, $\mathbf{s}$, can be expressed equivalently as a
product of vectors representing object shape and position, and thus the
problem of disentangling these factors essentially amounts to a vector
factorization problem.

The standard way to factorize the scene vector (e.g., as in [34]) would be to
use three codebooks corresponding to shape, horizontal position and vertical
position, for a total of $10+105*2=220$ codebook vectors. By contrast, a
residue number system with moduli $\\{3,5,7\\}$ uses 7 factors but only
$10+(3+5+7)*2=40$ vectors. Example runs of both problem setups are shown in
Figure 5b.

Figure 5c demonstrates the two main advantages of the residue resonator
compared to the standard resonator baseline: a reduction in both memory
requirements (as just described) and the required number of iterations.
Whereas the standard resonator takes over $\approx 2{,}000$ codebook
evaluations on average in our simulations, the residue resonator averages only
$\approx 800$ codebook evaluations. (Both dramatically improve over the brute
force search, which requires $110{,}250$ codebook evaluations). The key lesson
is that this simple change to a resonator network leads to a multiplicative
decrease in the number of computations required.

####  2.6.2 Generating exact solutions to the subset sum problem

Figure 6: The resonator network, with Residue Hyperdimensional Computing,
enables successful searches for solutions to the subset sum problem. a,
Demonstration of the resonator network on the subset sum problem, with
$S=\\{18,4,5,10,2,23\\}$. The resonator network converges to the correct
solution after a few iterations. b, Performance of the resonator network on
randomly selected subset sum problems with a fixed set size. c, The
probability of success scales as expected with independent trials on each
start. d, Success of the resonator network on the subset sum problem depends
on the range of items indexed in the subset sum problem (larger ranges are
harder to factorize). e, Performance of the resonator network in terms of the
expected number of iterations compared to chance. The expected number of
iterations for the resonator network scales favorably compared to brute force
search and improves with higher encoding dimension. f, Comparison of average
compute time of the resonator network versus an exact algorithm. The
initialization cost of setting up the resonator network is higher; even so,
for large set sizes, the resonator network is faster.

Here we apply Residue Hyperdimensional Computing to the subset sum problem.
Formally, the problem asks if a multiset of integers, $S$, contains a subset
$S^{*}$ that sums to a target integer $T$. A further demand is to return
$S^{*}$ if it exists. When all integers are positive, the subset-sum problem
is NP-complete [35]. It is a useful case to consider because there are well-
known polynomial-time reductions from other NP-complete problems (e.g., 3-SAT)
to subset sum [36].

To find solutions to the subset sum problem, we encode $T$ as a vector,
$\mathbf{z}(T)$, and use $|S|$ factors. Each factor, $F_{k}$, has a codebook
of two items: an identity vector $\mathbf{z}(0)$ and
$\mathbf{z}(S_{k})$—reflecting the binary decision to include that item in the
sum or not. Figure 6a demonstrates the resonator network successfully finding
the solution when $T=21$ and $|S|=6$.

In order to use a residue number system, we need to choose $M$ so that
$M>\sum_{s\in S}s$. This means that we need $\lceil\text{log}\ M\rceil$ bits
per unit. This is an improvement from previous work using resonator networks
to solve factorization problems [37], which requires floating point precision
to perform semi-prime factorization.

To understand the scaling capacity of the residue number system, we evaluate
the performance of the resonator network as the set size increases. We observe
that the resonator network finds exact solutions to the subset sum problem for
large sets, and that performance improves with higher vector dimension (Figure
6b). Figure 6c illustrates that the success probability after up to 10 trials
matches what is expected from 10 independent runs of the 1-trial accuracy.
This finding suggests that the resonator network constitutes a ‘Las Vegas’
algorithm [38], in which each run has a success probability $p$, $p$ is
independent across runs, and so the algorithm requires $1/p$ iterations on
average. Accuracy also depends on the integer range searched over, even for
the same set size (Figure 6d), perhaps because larger integer ranges reduce
the probability of multiple subsets matching the target.

Finally, we compare our subset sum algorithm to brute force search and an
exponential-time algorithm that solves the decision problem. The average
number of iterations required by the resonator network to find a solution is
drastically less than the exponentially increasing cost of a brute force
search (Figure 6e) and also improves with higher dimension. We find that on a
CPU, the resonator network has faster clock time than the exponential time
algorithm for $|S|>28$ (Figure 6f, Methods 4.6); most of the compute time is
spent on generating the vector representations, rather than the resonator
network dynamics itself. More significantly, whereas the baseline algorithm
required $\mathcal{O}(2^{|S|})$ memory to keep candidate subsets in memory,
the resonator network only requires $\mathcal{O}(D\cdot|S|)$ memory, since it
never needs to explicitly represent every subset. We emphasize that the CPU
implementation of the resonator network is primarily intended as a proof of
concept, and that further performance gains would likely result from
implementing the resonator network on emerging computing platforms, as in [34,
39].

##  3 Discussion

Our study provides the first definition of residue number systems with
hyperdimensional computing. The framework inherits the benefits of both
systems: the carry-free arithmetic and Chinese Remainder theoretic guarantees
from residue number systems, along with the robustness and computing-in-
superposition properties of hyperdimensional computing [40]. The framework
provides a favorable way to encode, transform, and decode variables that is
robust to noise. Taken together, these properties make Residue
Hyperdimensional Computing an appealing framework for efficiently solving
difficult computational problems, especially those involving combinatorial
optimization. It also has implications for modeling population codes in the
brain, in particular grid cells.

Prior work in computational neuroscience [11, 30] has emphasized that residue
number systems endow grid cells with useful computational properties,
including high spatial resolution, modular updates, and error correction. We
demonstrate that Residue Hyperdimensional Computing successfully achieves each
of these coding properties (Sections 2.2 & 2.3) in a working neural
implementation. Reciprocally, our algebraic framework makes two contributions
to computational neuroscience. First, we show how to extend a residue number
system to a self-consistent, non-negative hexagonal coordinate system. Second,
we provide a new algorithm for collectively coupling spatial position to grid
cell modules via the resonator network. The core prediction our framework
makes for systems neuroscience is that each grid cell module corresponds to a
factor estimate in the resonator network. More specifically, each module
implements a toroidal attractor network, and multiplicative couplings with
hippocampus and other grid cell modules enable error correction. This
prediction is consistent both with recent experimental analysis supporting the
existence of continuous attractor networks in single grid cell modules [41],
and with recent theories of joint attractor dynamics in hippocampus and medial
entorhinal cortex [42].

The framework also has implications for how neural populations can solve
difficult optimization problems, such as disentangling of visual scenes.
Recent work has emphasized the promise of hyperdimensional computing as an
abstraction for neuromorphic computing [43, 34, 40]. Residue Hyperdimensional
Computing substantially reduces the storage and average number of operations
needed for solving decoding problems and combinatorial optimization,
contributing a simple yet powerful improvement. In addition, the phasor
representations suggested by our framework directly map onto $Q$-state phasor
networks [44], suggesting promising implementations in spiking neural networks
[45] and strategies for solving combinatorial optimization problems [46].

Finally, the performance of our framework on the subset sum problem suggests a
new route for solving optimization problems with distributed representations
and unconventional hardware. The subset sum problem is a particularly good fit
for our framework because it is easily implemented by the Hadamard product
operation on high-dimensional vectors. Since other hard problems, such as
3-SAT, can be efficiently mapped to subset sum, our results potentially point
the way to a new class of parallel algorithms for efficiently solving NP-hard
problems.

##  4 Methods

###  4.1 Definitions of algebraic operations

####  4.1.1 Definition of additive binding operation

We implement addition within Residue Hyperdimensional Computing by using the
Hadamard product operation (component-wise multiplication, $\odot$) : i.e.,
$\mathbf{z}(x_{1}+x_{2})=\mathbf{z}(x_{1})\odot\mathbf{z}(x_{2})$. The
Hadamard product correctly implements addition because it is commutative. This
can be seen as follows:

$\displaystyle\mathbf{z}(x_{1})\odot\mathbf{z}(x_{2})$
$\displaystyle=(\mathbf{z}_{m_{1}}(x_{1})\odot\mathbf{z}_{m_{2}}(x_{1})\odot\dots\odot\mathbf{z}_{m_{K}}(x_{1}))\odot(\mathbf{z}_{m_{1}}(x_{2})\odot\mathbf{z}_{m_{2}}(x_{2})\odot\dots\odot\mathbf{z}_{m_{K}}(x_{2}))$
$\displaystyle=(\mathbf{z}_{m_{1}}(x_{1})\odot\mathbf{z}_{m_{1}}(x_{2}))\odot(\mathbf{z}_{m_{2}}(x_{1})\odot\mathbf{z}_{m_{2}}(x_{2}))\odot\dots\odot(\mathbf{z}_{m_{K}}(x_{1})\odot\mathbf{z}_{m_{K}}(x_{2}))$
$\displaystyle=\mathbf{z}_{m_{1}}(x_{1}+x_{2})\odot\mathbf{z}_{m_{2}}(x_{1}+x_{2})\odot\dots\odot\mathbf{z}_{m_{K}}(x_{1}+x_{2})$
$\displaystyle=\mathbf{z}(x_{1}+x_{2}).$

####  4.1.2 Definition of multiplicative binding operation

To implement a second binding operation ($\star$), such that
$\mathbf{z}(x_{1}\cdot x_{2})$ = $\mathbf{z}(x_{1})\star\mathbf{z}(x_{2})$,
every component of the vector $\mathbf{z}(x_{1})$ must be multiplied by
$x_{2}$. If we had the value of $x_{2}$ explicitly, then we could directly
implement $\mathbf{z}(x_{1}\cdot x_{2})$ by component-wise exponentiation of
$\mathbf{z}(x_{1})$ by $x_{2}$. However, decoding incurs additional
computational costs, and we show here that multiplications can be computed
without this intermediate step.

We require a few simplifying assumptions to define our multiplication
operation. First, we assume that we have access to the individual base vectors
for each moduli (e.g., $\mathbf{z}_{m_{1}}(x_{1})$). If we do not, then we can
use the resonator network to recover them. The key observation is that if $x$
is an integer, then each component of $\mathbf{z}_{m_{k}}(x)$ is itself a
$m_{k}$-th root of unity. More specifically, it equals
$e^{i\frac{2\pi}{m_{k}}r_{j}}$, for some integer
$r_{j}=\frac{m_{k}}{2\pi}\phi_{j}x\ (\text{mod}\ m_{k})$.

Therefore, we define an operation, $f$, that can multiply two discrete phases
when they are both drawn from the $m_{k}$-th roots of unity:
$f\left(e^{i\frac{2\pi}{m_{k}}r},e^{i\frac{2\pi}{m_{k}}s}\right)=e^{i\frac{2\pi}{m_{k}}rs}$.
When $f$ is applied to two vectors of the same dimension, the multiplication
is applied component-wise. Supposing that $r=\frac{m_{k}}{2\pi}\phi x_{1}$,
and $s=\frac{m_{k}}{2\pi}\phi x_{2}$, we obtain
$rs=\frac{m_{k}}{2\pi}\phi^{2}x_{1}x_{2}$, which is off from our desired
result by a multiplicative factor of $\phi$. This motivates a final step of
cancelling out this extra factor.

Because each phase $\phi$ is drawn from the $m_{k}$-th roots of unity, it can
be written as $\frac{2\pi}{m_{k}}u$, where $u\in\mathbb{Z}\ (\text{mod}\
m_{k})$. When $m_{k}$ is prime, then any non-zero integer $u$ has a unique
modular multiplicative inverse $v\in\mathbb{Z}\ (\text{mod}\ m_{k})$, such
that $u\times v=1$ (mod $m_{k}$). For example, the modular multiplicative
inverse of 3 (mod 5) is 2. Consequently,
$f\left(e^{i\frac{2\pi}{m_{k}}u^{2}},e^{i\frac{2\pi}{m_{k}}v}\right)=e^{i\frac{2\pi}{m_{k}}u}$.
Motivated by these reasons, we therefore assume that whenever multiplicative
binding is used, every moduli $m_{k}$ is prime. This assumption allows us to
define an “anti-base” vector, $\mathbf{y}_{m_{k}}$, whose components are
defined by the modular multiplicative inverses of $\mathbf{z}_{m_{k}}$. That
is, if the $j$-th component of $\mathbf{z}_{m_{k}}$ is
$e^{i\frac{2\pi}{m_{k}}u}$, then the $j$-th component of $\mathbf{y}_{m_{k}}$
is $e^{i\frac{2\pi}{m_{k}}v}$.

These assumptions motivate the following definition of the multiplicative
operation, which we show successfully performs the multiplication of the
arguments along with necessary cancellations:

$\displaystyle\mathbf{z}(x_{1})\star\mathbf{z}(x_{2})$
$\displaystyle:=f(f(\mathbf{z}_{m_{1}}(x_{1}),\mathbf{z}_{m_{1}}(x_{2})),\mathbf{y}_{m_{1}})\odot\dots\odot
f(f(\mathbf{z}_{m_{K}}(x_{1}),\mathbf{z}_{m_{K}}(x_{2})),\mathbf{y}_{m_{K}})$
$\displaystyle=[e^{i\phi_{1,1}^{2}x_{1}x_{2}\phi_{1,1}^{-1}},\dots,e^{i\phi_{1,D}^{2}x_{1}x_{2}\phi_{1,D}^{-1}}]\odot\dots\odot[e^{i\phi_{K,1}^{2}x_{1}x_{2}\phi_{K,1}^{-1}},\dots,e^{i\phi_{K,D}^{2}x_{1}x_{2}\phi_{K,D}^{-1}}]$
$\displaystyle=[e^{i\phi_{1,1}x_{1}x_{2}},\dots,e^{i\phi_{1,D}x_{1}x_{2}}]\odot\dots\odot[e^{i\phi_{K,1}x_{1}x_{2}},\dots,e^{i\phi_{K,D}x_{1}x_{2}}]$
$\displaystyle=\mathbf{z}(x_{1}\times x_{2})$

We implement $f$ by extracting $s$ by taking the angle of the phasor
$e^{i\frac{2\pi}{m_{k}}s}$, multiplying the angle by $\frac{m_{k}}{2\pi}$ and
exponentiating $e^{i\frac{2\pi}{m_{k}}r}$ by the result. We compute modular
multiplicative inverses via the built-in pow function in Python. However, we
note that both functions can also be implemented by lookup tables, and
precomputing all input-output pairs may be optimal when many computations are
re-used and lookups are inexpensive.

###  4.2 Decoding methods

In the context of high-dimensional distributed representations, the decoding
problem is to recover a variable $x$ from a distributed representation
$\mathbf{z}(x)$. In all of our decoding experiments, $x$ is either an integer
or rational number. A survey of decoding methods, applied to symbolic
hyperdimensional computing models, can be found in [23].

####  4.2.1 Codebook decoding

Codebook decoding estimates $x$ by taking inner products between
$\mathbf{z}(x)$ and a precomputed set of reference vectors:
$\hat{x}=\underset{x_{k}}{\operatorname*{arg\,max}}\langle\mathbf{z}(x),\mathbf{z}(x_{k})\rangle$.

####  4.2.2 Resonator network details

The resonator network [24, 25] is an algorithm for factoring an input vector,
$\mathbf{z}$, into the primitives
$\\{\mathbf{z}_{1},\mathbf{z}_{2},\dots,\mathbf{z}_{K}\\}$ that compose it via
Hadamard product binding:
$\mathbf{z}=\mathbf{z}_{1}\odot\mathbf{z}_{2}\odot\dots\odot\mathbf{z}_{K}$.
Each $\mathbf{z}_{j}$ is specified to come from a set of candidate vectors
concatenated in a codebook, $\mathbf{Z}_{j}$, and therefore the search space
grows combinatorially in terms of the product of codebook sizes.

The resonator network functions as a dynamical system with the following
update equations:

$\hat{\mathbf{z}}_{j}(t+1)=g(\mathbf{Z}_{j}\mathbf{Z}_{j}^{\dagger}(\mathbf{z}\odot\prod_{i\neq
j}\mathbf{z}_{i}^{\dagger}(t)))$ (6)

where $g$ is a non-linearity preserving phase and discarding angles of each
complex component. In the case where we are decoding the representation of a
residue number, the codebook for each factor (or moduli) $\mathbf{Z_{j}}$ is
composed of $m_{j}$ entries, which are the vector encodings for the residues
$[\mathbf{z}_{m_{j}}(0),\mathbf{z}_{m_{j}}(1),\dots,\mathbf{z}_{m_{j}}(m_{j}-1)]$.
In all experiments, we use an asynchronous update rule in which at each time
step, only one factor estimate is updated, and every set of time steps, each
vector is updated once. The algorithm runs either until convergence or until a
maximum number of iterations has been reached. We consider the resonator
converged when the normalized cosine similarity between two successive states
exceeds a threshold, $\alpha$ (for all experiments, $\alpha=0.95$).

####  4.2.3 Evaluation of resonator network decoding accuracy, capacity and
robustness to noise

We evaluate resonator network accuracy as a function of vector dimension
($D$), effective range ($M$), number of moduli ($K$), and noise level
(dependent on $\kappa$). We add noise only in experiments shown in Figure 2f,
and $K=2$ unless stated otherwise. $D=1024$ in Figure 2d, and $D=512$ in
Figure 2f. To compute data points for curves that are a function of $M$, we
generate a list of ascending primes, and select $K$ consecutive primes as
moduli. The effective range, $M$, is the product of these moduli. We continue
experiments for a fixed $D$ and increasingly large $M$ until empirical
accuracy falls below a given threshold (0.95 for Figures 2a and 2c, and 0.05
otherwise). To report the required number of comparisons for Figure 2a, we
normalize the average number of inner product iterations by the accuracy, and
visualize curves only in the high-accuracy regime (above 95 percent).

###  4.3 Hexagonal residue encodings

To project a $2D$ vector $\mathbf{x}$ to a $3D$ hexagonal coordinate
$\mathbf{y}$, we multiply it by a matrix $\Psi$:

$\Psi=\begin{bmatrix}-\sqrt{3}/2&-1/2\\\ \sqrt{3}/2&-1/2\\\ 0&1\end{bmatrix}$

The resulting vector $\mathbf{y}=\Psi\,\mathbf{x}$ is encoded as a high-
dimensional vector using the generalization to multiple dimensions specified
in Equation 5. As an additional constraint, we require that
$\mathbf{z}([1,1,1])=\mathbf{z}([0,0,0])$. This condition implements the self-
cancellation property (i.e, that moving equally along the three equiangular
directions cancels out). It also converts possible negative values arising
from the projection step to an equivalent non-negative coordinate encoding.
Somewhat fortuitously, this constraint is naturally enforced in RHC by
ensuring that for each component of $\mathbf{z}$, the three phases
corresponding to the three directions sum to 0 (mod $2\pi$). This is achieved
by constraining the joint probability distribution over triplets of phases so
that this requirement is met (Figures 3a and 3c).

###  4.4 Decoding with sub-integer precision

####  4.4.1 Extension of encoding scheme to rational numbers

For a rational number $q\in\mathbb{Q}$, we define $\mathbf{z}_{m}(q)$ based on
Fractional Power Encoding, modulo $m$ as:

$\mathbf{z}_{m}(q)=[e^{i\phi_{1}q},e^{i\phi_{2}q},\dots,e^{i\phi_{D}q}]$ (7)

We then form our representation of the remainders (modulo $m_{k}$) via the
same process described in Equation 4. If $q$ is an integer, then this
procedure matches that of Definition 2.2.1. But in general,
$\mathbf{z}_{m}(q)\neq(\mathbf{z}_{m})^{q}$. This is significant because while
we can still evaluate similarity via inner products and perform addition
operations, multiplication operations are no longer well defined.

####  4.4.2 Sub-integer decoding with the resonator network

Sub-integer decoding with the resonator network proceeds in three steps.
First, we let the resonator update its factor estimates until convergence to a
fixed point. We emphasize that sub-integer encodings are also fixed points of
the resonator, even when the codebooks of the resonator contain only integers
(Section 2.5). Second, we find the nearest integer codebook for each moduli,
and generate the nearest codebooks for the fractional values within range 1 of
that decoded integer ($r$ in total). Third, we use codebook decoding over
these vectors encoding fractional values to return our result.

####  4.4.3 Evaluation of sub-integer decoding with noise

We fix $D=512$, and let $\kappa=\\{16.0,1.0\\}$. We run the resonator network
until a maximum number of iterations or convergence, and evaluate if both the
nearest integer and nearest fractional state are correct. If so, we regard the
solution as correct, reporting accuracy and bits per vector.

####  4.4.4 Measuring bits per vector

To measure the total amount of information decoded, we account for the
accuracy of decoding and the number of states distinguished. The amount of
information decoded for a single number (denoted as $I_{\mathrm{num}}$) is
calculated using the corresponding accuracy ($a$) and size ($P$) of the total
search space as

$\begin{split}I_{\mathrm{num}}(a,P)=&a\log_{2}(Pa)+(1-a)\log_{2}\left(\frac{P}{P-1}(1-a)\right).\end{split}$
(8)

For a detailed derivation of this equation, please refer to Section 2.2.3 of
[31]. According to this metric, the amount of decoded information is $0$ when
the accuracy is at chance ($1/P$).

###  4.5 Visual scene factorization experiments

Convolutional sparse coding learns a dictionary of basis functions,
$\\{\phi_{j}(x,y)\\}$, and infers a set of sparse latent representations,
$\\{A_{j}(x,y)\\}$, for each image, $I(x,y)$, by minimizing the following
energy function, E:

$E=\frac{1}{2}||I-\sum_{j=1}^{n}\phi_{j}\ast
A_{j}||^{2}_{2}+\lambda\sum_{j=1}^{n}||A_{j}||_{1}$ (9)

where $\ast$ denotes convolution, and $\lambda$ is a hyperparameter weighting
the tradeoff between reconstruction error vs. sparsity. We use the SPORCO
implementation of convolutional sparse coding introduced by [47] to learn the
$\\{\phi_{j}(x,y)\\}$ for an ensemble of MNIST digits, and to infer the sparse
representation $A$ for each image $I$.

A useful feature of convolutional sparse coding is its equivariance to 2D
translation; that is, 2D translation in the image domain results in 2D
translation of the sparse representations, $\\{A_{j}(x,y)\\}$. We can thus
convert the set of sparse feature maps $\\{A_{j}(x,y)\\}$ to a high-
dimensional vector as follows:

$\mathbf{s}=\sum_{j,x,y}\textbf{h}(x)\odot\textbf{v}(y)\odot\textbf{d}_{j}\cdot
A_{j}(x,y)$ (10)

Here $\textbf{h}(x)$ and $\textbf{v}(y)$ denote the RHC encodings of
horizontal ($x$) and vertical ($y$) position. $\mathbf{d}_{j}$ is a random
vector generated i.i.d. that represents the identity of each basis function
$\phi_{j}$. By expectation, most values of each $A_{j}(x,y)$ will be zero,
because the energy function for sparse coding penalizes non-zero coefficients.
Thus, the scene vector, $\mathbf{s}$, can be seen as a sparse superposition of
position encodings of features ($\mathbf{\phi_{i}}$) contained in the image.

Now, we can separately define the vector encoding of each object $i$ to be
recognized as

$\mathbf{O}^{(i)}=\sum_{j,x,y}\textbf{h}(x)\odot\textbf{v}(y)\odot\textbf{d}_{j}\cdot
o^{(i)}_{j}(x,y)$ (11)

where $o^{(i)}(x,y)$ is the sparse representation of the image of object $i$
within a canonical reference frame. If we were to place object $i$ at position
$(x^{\prime},y^{\prime})$ within an image, the resulting scene vector computed
according to equation 10 will be given as

$\mathbf{s}=\mathbf{h}(x^{\prime})\odot\mathbf{v}(y^{\prime})\odot\mathbf{O}^{(i)}$
(12)

Therefore, our scene analysis problem amounts to one of factorizing
$\mathbf{s}$ into its constituent vectors $\mathbf{h}(x^{\prime})$,
$\mathbf{v}(y^{\prime})$, and $\mathbf{O}^{(i)}$.

We can factorize $\mathbf{s}$ using a resonator network with three codebooks,
$\mathbf{O},\mathbf{H},$ and $\mathbf{V}$. Each element
$\mathbf{O}^{(i)}\in\mathbf{O}$ consists of an encoding of each object as
above, and $\mathbf{H}$ and $\mathbf{V}$ contain RHC encodings of horizontal
and vertical position.

For our object examples, we use 10 images from the MNIST dataset. Sparse
coding dictionary elements are optimized over a subset of the MNIST dataset.
After inferring a sparse code for each image, we encode it as a high-
dimensional vector ($D=10,000$). We use a residue number system with bases
$\\{3,5,7\\}$ for both horizontal and vertical dimension and then either
enumerate all 105 codebooks for a single factor (Standard) or use 3 factors
with 3, 5, and 7 codebooks, respectively (Residue). In either case, we ran the
resonator network until convergence to a vector matching the scene
representation (including reinitialization, if it did not converge after a
fixed number of iterations or became stuck in a local minima), and record the
average number of iterations multiplied by the average number of codebook
evaluations (which is smaller for the residue encoding).

###  4.6 Subset sum experiments

We use a residue number system with 3 moduli, $\\{m-1,m,m+1\\}$, where $m$ is
a positive integer, ensuring that our moduli are co-prime. To generate random
subset sum problems, we first define a maximum sum range to be $M/2$. For
Figures 6b and 6c, $m=200$, $M\approx 200^{3}$. Then, we draw random variables
from a uniform distribution (scaled between 0, and half of the maximum sum
over the largest set size tested). We then select a random subset of the set
(all subsets are equally likely) and compute the sum. This sum forms the input
to the resonator network, and we treat its solution is correct if it converged
to the same sum. If the resonator network returns the wrong output, we restart
it from a different random initialization, up to a maximum number of trials.
We vary both the vector dimension ($D$) and set size ($|S|$), reporting
accuracy after multiple simulations. For Figure 6d, $D=400$.

To compare the number of evaluations relative to brute force (Figure 6e), we
record the average number of evaluations on each set size. We divide the
number of inner product comparisons required for brute force evaluation by the
number of comparisons per resonator network iteration. Further, we normalize
the number of resonator iterations by the accuracy to ensure a fair
comparison. In comparing our algorithm to a solver, we implement an exact
subset-sum algorithm as a baseline [48]. We let $m=1{,}000$,
$D=\\{10{,}000,20{,}000\\}$, and draw integers uniformly from the range
$[0,5000]$.

## References

  * [1] Ashish Vaswani et al. “Attention is all you need” In _Advances in Neural Information Processing Systems_ 30, 2017
  * [2] Vincent Sitzmann et al. “Implicit neural representations with periodic activation functions” In _Advances in neural information processing systems_ 33, 2020, pp. 7462–7473
  * [3] Matthew Tancik et al. “Fourier features let networks learn high frequency functions in low dimensional domains” In _Advances in Neural Information Processing Systems_ 33, 2020, pp. 7537–7547
  * [4] J. Buckman, A. Roy, C. Raffel and I. Goodfellow “Thermometer Encoding: One Hot Way to Resist Adversarial Examples” In _International Conference on Learning Representations (ICLR)_ , 2018, pp. 1–22
  * [5] C. Diao, D. Kleyko, J.. Rabaey and B.. Olshausen “Generalized Learning Vector Quantization for Classification in Randomized Neural Networks and Hyperdimensional Computing” In _International Joint Conference on Neural Networks_ , 2021, pp. 1–9
  * [6] Alexandre Pouget, Peter Dayan and Richard Zemel “Information processing with population codes” In _Nature Reviews Neuroscience_ 1.2 Nature Publishing Group, 2000, pp. 125–132
  * [7] Nikolaus Kriegeskorte, Marieke Mur and Peter A Bandettini “Representational similarity analysis-connecting the branches of systems neuroscience” In _Frontiers in Systems Neuroscience_ Frontiers, 2008, pp. 4
  * [8] Blake Bordelon and Cengiz Pehlevan “Population codes enable learning from few examples by shaping inductive bias” In _Elife_ 11 eLife Sciences Publications Limited, 2022, pp. e78606
  * [9] Torkel Hafting et al. “Microstructure of a spatial map in the entorhinal cortex” In _Nature_ 436.7052 Nature Publishing Group, 2005, pp. 801–806
  * [10] Alexander Mathis, Andreas VM Herz and Martin B Stemmler “Resolution of nested neuronal representations can be exponential in the number of neurons” In _Physical review letters_ 109.1 APS, 2012, pp. 018103
  * [11] Ila R Fiete, Yoram Burak and Ted Brookings “What grid cells convey about rat location” In _Journal of Neuroscience_ 28.27 Soc Neuroscience, 2008, pp. 6858–6871
  * [12] T.. Plate “Distributed Representations and Nested Compositional Structure” University of Toronto, PhD Thesis, 1994
  * [13] R.. Gayler “Vector Symbolic Architectures Answer Jackendoff’s Challenges for Cognitive Neuroscience” In _Joint International Conference on Cognitive Science (ICCS/ASCS)_ , 2003, pp. 133–138
  * [14] Pentti Kanerva “Hyperdimensional computing: An introduction to computing in distributed representation with high-dimensional random vectors” In _Cognitive computation_ 1 Springer, 2009, pp. 139–159
  * [15] Harvey L Garner “The residue number system” In _Western Joint Computer Conference_ , 1959, pp. 146–153
  * [16] Oded Goldreich, Dana Ron and Madhu Sudan “Chinese remaindering with errors” In _Thirty-first Annual ACM Symposium on Theory of Computing_ , 1999, pp. 225–234
  * [17] Tony A Plate “Holographic Reduced Representation: Distributed representation for cognitive structures” CSLI Publications Stanford, 2003
  * [18] Ali Rahimi and Benjamin Recht “Random Features for Large-Scale Kernel Machines” In _Advances in Neural Information Processing Systems_ 20, 2007, pp. 1–8
  * [19] E.. Frady et al. “Computing on Functions Using Randomized Vector Representations” In _arXiv:2109.03429_ , 2021, pp. 1–33
  * [20] E.. Frady et al. “Computing on Functions Using Randomized Vector Representations (in brief)” In _Neuro-Inspired Computational Elements Conference_ , 2022, pp. 115–122
  * [21] Amos R Omondi and A Benjamin Premkumar “Residue number systems: theory and implementation” World Scientific, 2007
  * [22] Denis Kleyko et al. “Vector symbolic architectures as a computing framework for emerging hardware” In _Proceedings of the IEEE_ 110.10 IEEE, 2022, pp. 1538–1571
  * [23] Denis Kleyko et al. “Efficient decoding of compositional structure in holistic representations” In _Neural Computation_ 35.7 MIT Press, 2023, pp. 1159–1186
  * [24] E.. Frady, S.. Kent, B.. Olshausen and F.. Sommer “Resonator Networks, 1: An Efficient Solution for Factoring High-Dimensional, Distributed Representations of Data Structures” In _Neural Computation_ 32.12, 2020, pp. 2311–2331
  * [25] Spencer J Kent, E Paxon Frady, Friedrich T Sommer and Bruno A Olshausen “Resonator networks, 2: Factorization performance and capacity compared to optimization-based methods” In _Neural computation_ 32.12 MIT Press One Rogers Street, Cambridge, MA 02142-1209, USA journals-info …, 2020, pp. 2332–2388
  * [26] Edmund Landau “Über die maximalordnung der permutationen gegebenen grades” In _Archiv der Math. und Phys_ 3, 1903, pp. 92–103
  * [27] Alexander Mathis, Martin B Stemmler and Andreas VM Herz “Probable nature of higher-dimensional symmetries underlying mammalian grid-cell activity patterns” In _Elife_ 4 eLife Sciences Publications Limited, 2015, pp. e05979
  * [28] Vassili Nikolaevich Malozemov and Aleksandr Borisovich Pevnyi “Equiangular tight frames” In _Journal of Mathematical Sciences_ 157.6 Springer, 2009, pp. 789–815
  * [29] Brent Komer “Biologically Inspired Spatial Representation”, 2020
  * [30] Sameet Sreenivasan and Ila Fiete “Grid cells generate an analog error-correcting code for singularly precise neural computation” In _Nature Neuroscience_ 14.10 Nature Publishing Group, 2011, pp. 1330–1337
  * [31] E Paxon Frady, Denis Kleyko and Friedrich T Sommer “A theory of sequence indexing and working memory in recurrent neural networks” In _Neural Computation_ 30.6 MIT Press One Rogers Street, Cambridge, MA 02142-1209, USA journals-info …, 2018, pp. 1449–1513
  * [32] Bruno A. Olshausen “Perception as an Inference Problem” In _The Cognitive Neurosciences_ MIT Press, 2014, pp. 295–304
  * [33] Bruno A Olshausen and David J Field “Emergence of simple-cell receptive field properties by learning a sparse code for natural images” In _Nature_ 381.6583 Nature Publishing Group, 1996, pp. 607–609
  * [34] Alpha Renner et al. “Neuromorphic visual scene understanding with resonator networks” In _arXiv:2208.12880_ , 2022
  * [35] Jon Kleinberg and Eva Tardos “Algorithm design” Pearson Education India, 2006
  * [36] Richard Karp “chapter Reducibility among Combinatorial Problems” In _Complexity of Computer Computations_ Plenum, 1972
  * [37] Denis Kleyko et al. “Integer factorization with compositional distributed representations” In _Neuro-Inspired Computational Elements Conference_ , 2022, pp. 73–80
  * [38] László Babai “Monte-Carlo algorithms in graph isomorphism testing” In _Université tde Montréal Technical Report, DMS_ , 1979
  * [39] Jovin Langenegger et al. “In-memory factorization of holographic perceptual representations” In _Nature Nanotechnology_ 18.5 Nature Publishing Group UK London, 2023, pp. 479–485
  * [40] Denis Kleyko et al. “Vector symbolic architectures as a computing framework for emerging hardware” In _Proceedings of the IEEE_ 110.10 IEEE, 2022, pp. 1538–1571
  * [41] Richard J Gardner et al. “Toroidal topology of population activity in grid cells” In _Nature_ 602.7895 Nature Publishing Group UK London, 2022, pp. 123–128
  * [42] Haggai Agmon and Yoram Burak “A theory of joint attractor dynamics in the hippocampus and the entorhinal cortex accounts for artificial remapping and grid cell field-to-field variability” In _Elife_ 9 eLife Sciences Publications, Ltd, 2020, pp. e56894
  * [43] E Paxon Frady and Friedrich T Sommer “Robust computation with rhythmic spike patterns” In _Proceedings of the National Academy of Sciences_ 116.36 National Acad Sciences, 2019, pp. 18050–18059
  * [44] André J Noest “Discrete-state phasor neural networks” In _Physical Review A_ 38.4 APS, 1988, pp. 2196
  * [45] Connor Bybee and Friedrich Sommer “Optimal Oscillator Memory Networks” In _Neuro-Inspired Computational Elements Conference_ , 2022, pp. 81–83
  * [46] Tianshi Wang and Jaijeet Roychowdhury “OIM: Oscillator-based Ising machines for solving combinatorial optimisation problems” In _Unconventional Computation and Natural Computation: 18th International Conference, UCNC 2019, Tokyo, Japan, June 3–7, 2019, Proceedings 18_ , 2019, pp. 232–256 Springer
  * [47] Brendt Wohlberg “SPORCO: A Python package for standard and convolutional sparse representations” In _Python in Science Conference_ 15, 2017, pp. 1–8
  * [48] Soumendra Nanda “Subset Sum Problem”, 2005 URL: https://cs.dartmouth.edu/~ac/Teach/CS105-Winter05/Notes/nanda-scribe-3.pdf
  * [49] D.. Rachkovskij “Linear Classifiers based on Binary Distributed Representations” In _Information Theories and Applications_ 14.3, 2007, pp. 270–274
  * [50] O. Räsänen and J. Saarinen “Sequence Prediction with Sparse Distributed Hyperdimensional Coding Applied to the Analysis of Mobile Phone Use Patterns” In _IEEE Transactions on Neural Networks and Learning Systems_ 27.9, 2016, pp. 1878–1889
  * [51] D. Kleyko, E. Osipov, N. Papakonstantinou and V. Vyatkin “Hyperdimensional Computing in Industrial Systems: The Use-Case of Distributed Fault Isolation in a Power Plant” In _IEEE Access_ 6, 2018, pp. 30766–30777
  * [52] A. Rahimi, P. Kanerva, L. Benini and J.. Rabaey “Efficient Biosignal Processing Using Hyperdimensional Computing: Network Templates for Combined Learning and Classification of ExG Signals” In _Proceedings of the IEEE_ 107.1, 2019, pp. 123–143
  * [53] Kaspar A Schindler and Abbas Rahimi “A Primer on Hyperdimensional Computing for iEEG Seizure Detection” In _Frontiers in Neurology_ Frontiers, 2021, pp. 1–12
  * [54] D. Kleyko, D.. Rachkovskij, E. Osipov and A. Rahimi “A Survey on Hyperdimensional Computing aka Vector Symbolic Architectures, Part II: Applications, Cognitive Models, and Challenges” In _ACM Computing Surveys_ 55.9, 2023, pp. 1–52
  * [55] D. Kleyko, D.. Rachkovskij, E. Osipov and A. Rahimi “A Survey on Hyperdimensional Computing aka Vector Symbolic Architectures, Part I: Models and Data Transformations” In _ACM Computing Surveys_ 55.6, 2022, pp. 1–40
  * [56] P.. Penz “The closeness code: An integer to binary vector transformation suitable for neural network algorithms” In _IEEE First Annual International Conference on Neural Networks_ , 1987, pp. 515–522
  * [57] A.. Goltsev “An assembly neural network for texture segmentation” In _Neural Networks_ 4.9, 1996, pp. 643–653
  * [58] D. Smith and P. Stanford “A Random Walk in Hamming Space” In _International Joint Conference on Neural Networks_ 2, 1990, pp. 465–470
  * [59] D.. Rachkovskij, S.. Slipchenko, E.. Kussul and T.. Baidyk “Sparse Binary Distributed Encoding of Scalars” In _Journal of Automation and Information Sciences_ 37.6, 2005, pp. 12–23
  * [60] D. Kleyko et al. “Density Encoding Enables Resource-Efficient Randomly Connected Neural Networks” In _IEEE Transactions on Neural Networks and Learning Systems_ 32.8, 2020, pp. 3777–3783
  * [61] Denis Kleyko et al. “Classification and recall with binary hyperdimensional computing: Tradeoffs in choice of density and mapping characteristics” In _IEEE Transactions on Neural Networks and Learning Systems_ 29.12, 2018, pp. 5880–5898

## Acknowledgments

We thank Anthony Thomas, Eric Weiss, Josh Cynamon, and Amir Khosrowshahi for
insightful discussions and feedback. The work of CJK was supported by the
Department of Defense (DoD) through the National Defense Science & Engineering
Graduate (NDSEG) Fellowship Program. The work of DK, CB, FTS, and BAO was
supported in part by Intel’s THWAI program. The work of CJK, CB, PK, and BAO
was supported by the Center for the Co-Design of Cognitive Systems (CoCoSys),
one of seven centers in JUMP 2.0, a Semiconductor Research Corporation (SRC)
program sponsored by DARPA. DK has received funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Sklodowska-
Curie grant agreement No 839179. The work of FTS was supported in part by NIH
under Grant R01-EB026955, and in part by NSF under Grant IIS-118991.

## Corresponding authors

Correspondence to Christopher J. Kymn or Bruno A. Olshausen.

Supplemental material

##  Appendix A A brief survey of distributed coding schemes

In order to process vector representations of numbers, such as in machine
learning settings [49, 50, 51, 52, 53, 54], previous work combined
hyperdimensional computing with different kinds of locality-preserving
encodings for representing numeric data with vectors. The requirement to be
locality-preserving is that inner products between vectors encode similarity
of the underlying data. Here we briefly review some locality-preserving
encoding schemes that have been used in the past (see also [55]), assessing
their kernel properties.

Encoding scheme: | Algebra | Expressivity |  Efficient decoding |  Robust to noise
---|---|---|---|---
One-hot | ✗ | ✗ | ✓ | ✗
Thermometer [56], Supplement A.1 | ✗ | ✗ | ✓ | $\sim$
Float [57], Supplement A.2 Gaussian Population Codes [6] | ✗ | ✗ | ✓ | $\sim$
Scatter [58], Supplement A.3 | ✗ | ✗ | ✗ | ✓
Fractional Power Encoding [12], Section 2.1 | $\sim$ | ✓ | ✗ | ✓
Residue Hyperdimensional Computing | ✓ | ✓ | ✓ | ✓

Table S.1: Existing high-dimensional vector-based schemes for encoding numbers
(first five rows) in comparison to our proposed framework (last row).

###  A.1 The thermometer code

Figure S.1:  Similarity kernel of the thermometer code shown for several
levels; $D$ was set to $50$.

The thermometer code [56, 59, 4, 60] is a simple and structured way to form a
locality-preserving encoding for a range of discrete levels $s$, $s\in[0,D]$.
The first code $\mathbf{z}(0)$ consists of all -1s. For other levels, the
components of $\mathbf{z}(s)$ are determined as:

$z_{i}(s)=\begin{cases}+1,&i\leq s\\\ -1,&\text{otherwise}\end{cases}$ (S.1)

Thus, the last code $\mathbf{z}(D)$ consists of all +1s, and, in total, the
thermometer code can represent $D+1$ levels. Figure S.1 shows how cosine
similarity appears for several different levels when $D=50$. Thermometer codes
produce a translation-invariant kernel that is triangular and has a width of
$2D+1$ levels. It is a nonlocal kernel, in the sense that there are no two
points in the encoding range that have a similarity of zero. In practice,
thermometer codes are used commonly when applying hyperdimensional computing
to classification problems.

###  A.2 The float code

Figure S.2:  Similarity kernel of the float code shown for several levels; $D$
was set to $60$ while $w$ was set to $10$.

The float code, also known as the sliding code, [57, 59] addresses the issue
of the thermometer code, i.e., that the similarity decay is not local. This is
done by using $w$ consecutive +1 components (“float”) where the size of $w$
regulates similarity characteristics of the code. For the binary case, the
similarity kernel of the float code is the triangular kernel of width $2w+1$
levels. To encode the lowest value $\mathbf{z}(0)$, the first $w$ components
of the vector are set to +1s while the rest of the components are 0s. In
general, the components of $\mathbf{z}(s)$ are determined as:

$z_{i}(s)=\begin{cases}+1,&s\leq i<s+w\\\ 0,&\text{otherwise}\end{cases}$
(S.2)

Figure S.2 depicts how similarity (inner product normalized by $w$) decays for
several levels in the float code for $D=60$, $w=10$. The float code also
produces a triangular kernel, but in contrast to the thermometer code, it
allows controlling the width of the triangular kernel. The number of levels it
could encode is still limited and equals to $n-w+1$.

###  A.3 The scatter code

Figure S.3:  Similarity kernel of a scatter code; $D$ was set to $1000$, $p$
was $0.05$. The values of similarities were averaged over $50$ random
initializations of the code.

Scatter codes [58, 59, 61] are another alternative to form a locality-
preserving encoding where similarity decays nonlinearly. In scatter codes, the
code for the first level $\mathbf{z}(0)$ is chosen randomly while each
subsequent code is obtained from the previous one by randomly swapping its
components with some probability $p$:

$z_{i}(s)=\begin{cases}-z_{i}(s-1),&r_{i}\leq p\\\
z_{i}(s-1),&\text{otherwise}\end{cases}$ (S.3)

where $r_{i}$ is a random value for the $i$-th component of $\mathbf{z}(s)$
chosen from the uniform distribution. Note that potentially there is no
limitation on how many levels can be created with the scatter codes.

Figure S.3 shows how the cosine similarity looks for several different levels
formed with the scatter code. Interestingly, the kernels are ‘bell-shaped”
with the exact shape depending on the parameter settings. To better figure out
which standard kernel will correspond to this similarity, we have empirically
fitted three kernels: exponentiated triangular:

$K(s_{1},s_{2})=(1-\gamma|s_{2}-s_{1}|)^{\alpha};$ (S.4)

squared exponential:

$K(s_{1},s_{2})=e^{-\frac{(s_{2}-s_{1})^{2}}{2l^{2}}};$ (S.5)

and rational quadratic:

$K(s_{1},s_{2})=\left(1+\frac{(s_{2}-s_{1})^{2}}{2\alpha
l^{2}}\right)^{-\alpha};$ (S.6)

The parameters of the kernels were chosen using the mean squared error as the
fit criterion.

##  Appendix B Kernel properties of Residue Hyperdimensional
Computing

To show that Fractional Power Encoding (modulo $m$) results in approximation
of a particular periodic kernel with period $m$, we observe that our
probability distribution can be written as a Dirac comb function pointwise
multiplied by a rect function. This fact becomes useful when we see that our
kernel approximates a Fourier integral. Letting $x=x_{1}-x_{2}$, we take the
following steps to show convergence in the infinite-dimensional limit:

$\displaystyle K_{m}^{*}(x_{1},x_{2})$
$\displaystyle=\lim_{D\to\infty}\frac{1}{D}\sum_{d=1}^{D}e^{i\phi_{d}x_{1}}\overline{e^{i\phi_{d}x_{2}}}$
$\displaystyle=\lim_{D\to\infty}\frac{1}{D}\sum_{d=1}^{D}e^{i\phi_{d}(x_{1}-x_{2})}$
$\displaystyle=\int e^{i\phi(x_{1}-x_{2})}p(\phi)d\phi$
$\displaystyle=\mathscr{F}^{-1}[p(\phi)](x)$
$\displaystyle=\mathscr{F}^{-1}[\left(\frac{1}{m}(\sum_{s\in\mathbb{Z}}\delta(\omega-\frac{2\pi}{m}s)\right)\cdot\left(\text{rect}(\frac{\omega}{2\pi})\right)](x)$
$\displaystyle=[\sum_{s\in\mathbb{Z}}\delta(x-ms)]\circledast\text{sinc}(x)$
$\displaystyle=\sum_{s\in\mathbb{Z}}\text{sinc}(x-ms)$

Thus, our kernel is a ‘sinc comb’ function: a sum of sinc functions spaced
with a period of $m$. This result is particularly notable because sinc
evaluates to 0 for integers that are not a multiple of $m$ and means that
distinct integers (and remainders) are orthogonal in the high-dimensional
space.

To simplify the equation even further, we can derive a sum-less expression for
an infinite number of sinc functions. We need to consider two cases: 1) $x=ms$
for some $s\in\mathbb{Z}$, and 2) $x\neq ms$ for all $s\in\mathbb{Z}$.

Case 1 is straightforward. Without loss of generality, let $x$ be the value
for which $x-ms=0$. Then we have:

$\displaystyle K_{m}(x)$
$\displaystyle=\sum_{s\in\mathbb{Z}}\text{sinc}(x-ms)$
$\displaystyle=1+\sum_{n\in\mathbb{N}}\text{sinc}(-mn)+\text{sinc}(mn)$
$\displaystyle=1\,\text{(since sinc evaluates to 0 for non-zero integers)}$

For case 2, without loss of generality, let $0<|x|\leq m/2$. The answer
differs subtly depending on whether $m$ is even or odd. Let us derive the even
case first:

$\displaystyle K_{m}(x)$
$\displaystyle=\sum_{s\in\mathbb{Z}}\text{sinc}(x-ms)$
$\displaystyle=\sum_{s\in\mathbb{Z}}\frac{\text{sin}(\pi(x+ms))}{\pi(x+ms)}$
$\displaystyle=\frac{1}{\pi}\sum_{s\in\mathbb{Z}}\frac{\text{sin}(\pi
x)}{x+ms}$ $\displaystyle=\frac{\text{sin}(\pi
x)}{\pi}\sum_{s\in\mathbb{Z}}\frac{1}{x+ms}$
$\displaystyle=\frac{\text{sin}(\pi x)}{\pi
m}\sum_{s\in\mathbb{Z}}\frac{1}{\frac{x}{m}+s}$
$\displaystyle=\frac{\text{sin}(\pi x)}{\pi m}\pi\text{cot}(\frac{\pi x}{m})$
$\displaystyle=\frac{1}{m}\text{sin}(\pi x)\text{cot}(\frac{\pi x}{m})$

where the infinite sum is replaced due to an identity via the Herglotz trick:
$\pi\text{cot}(\pi x)=\sum_{s\in\mathbf{Z}}\frac{1}{x+s}$. We can also use a
second Herglotz identity, that $\pi\text{csc}(\pi
x)=\sum_{n\in\mathbf{N}}\frac{(-1)^{n}}{x+n}$, to solve for the odd case:

$\displaystyle K_{m}(x)$
$\displaystyle=\sum_{s\in\mathbb{Z}}\text{sinc}(x-ms)$
$\displaystyle=\sum_{s\in\mathbb{Z}}\frac{\text{sin}(\pi(x-ms))}{\pi(x-ms)}$
$\displaystyle=\frac{\text{sin}(\pi
x)}{\pi}\sum_{s\in\mathbb{Z}}\frac{(-1)^{n}}{x-ms}$
$\displaystyle=\frac{\text{sin}(\pi
x)}{\pi}\sum_{s\in\mathbb{Z}}\frac{(-1)^{n}}{x+ms}$
$\displaystyle=\frac{\text{sin}(\pi x)}{\pi
m}\sum_{s\in\mathbb{Z}}\frac{(-1)^{n}}{\frac{x}{m}+s}$
$\displaystyle=\frac{1}{m}\text{sin}(\pi x)\text{csc}(\frac{\pi x}{m})$

We confirm our result by comparing our analytic values for both cases to the
kernel induced by a high-dimensional vector (Figure S.4). In the limit of
$m\to\infty$, both kernels converge to the sinc function:
$K(x)=\frac{\text{sin}(\pi x)}{\pi x}$.

Figure S.4: The analytic kernel expected by dashed lines matches the
approximate kernel generated by a random vector of sufficiently high dimension
($D$=50,000). a, match for an odd modulus ($m=5$), b, match for an even
modulus ($m=6$).