# A Framework for Analyzing Cross-correlators using Price’s Theorem and
Piecewise-Linear Decomposition

Zhili Xiao
Electrical & Systems Engineering
Washington University in Saint Louis
Saint Louis, MO 63130
<EMAIL_ADDRESS>
&Shantanu Chakrabartty
Electrical & Systems Engineering
Washington University in Saint Louis
Saint Louis, MO 63130
<EMAIL_ADDRESS>

###### Abstract

Precise estimation of cross-correlation or similarity between two random
variables lies at the heart of signal detection, hyperdimensional computing,
associative memories, and neural networks. Although a vast literature exists
on different methods for estimating cross-correlations, the question what is
the best and the simplest method to estimate cross-correlations using finite
samples ? is still unclear. In this paper, we first argue that the standard
empirical approach might not be optimal, even though the estimator exhibits
uniform convergence to the true cross-correlation. Instead, we show that there
exists a large class of simple non-linear functions that can be used to
construct cross-correlators with a higher signal-to-noise ratio (SNR). To
demonstrate this, we first present a general mathematical framework using
Price’s Theorem that allows us to analyze arbitrary cross-correlators
constructed using a mixture of piece-wise linear functions. Using this
framework and a high-dimensional mapping, we show that some of the most
promising cross-correlators are based on Huber’s loss functions, margin-
propagation (MP) functions, and the log-sum-exponential (LSE) functions.

††footnotetext: All correspondences related to this manuscript should be
addressed to<EMAIL_ADDRESS>
##  1 Introduction

Estimating cross-correlations between random variables play an important role
in the field of statistics [1], machine learning [2, 3, 4, 5], and signal
detection [6, 7, 8, 9]. This is because the cross-correlation metric measures
some form of similarity between the random variables, revealing how one might
influence the other. With proper normalization, the metric becomes equivalent
to cosine similarity and unitary transforms, both of which are extensively
used in linear algebra [10], natural language processing [11], and computer
vision [12, 13]. In computer vision and signal processing, cross-correlation
is often used for feature extraction [14], where higher precision implies more
information is retained for further data processing and learning. Accurate and
efficient cross-correlation for pattern recognition and template matching [15,
16] also ensures reliable and real-time decision-making in applications like
radar detection or object recognition. In the emerging field of
hyperdimensional computing [17, 18], cross-correlations (or equivalently
inner-products) are used for information retrieval from sparse distributed
memories [19, 20]. Since most of the vectors in higher-dimensions are nearly
orthogonal to each other, precision in cross-correlation is essential to
discriminate between patterns and improve the robustness of learning [21, 22].

In its most general form, cross-correlation
$R:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$ is defined for a pair of
random variables $X\in\mathbb{R}$, $Y\in\mathbb{R}$ as

$R\mathrel{:=}\mathcal{E}[XY]=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}xy\enskip
p(X=x,Y=y)dxdy,$ (1)

where $p:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}^{+}$ denotes the
underlying joint probability distribution from which $x$ and $y$ are drawn
from. The operator $\mathcal{E}[.]$ denotes an expectation under the
probability measure $p$.

In practice, the joint distribution $p$ is not known apriori. Instead, one has
access to $N$ samples independently drawn from the distribution $p$, as
illustrated in Fig. 1c. If we denote the sample vectors as
$\bf{x}\in\mathbb{R}^{N}$ and $\bf{y}\in\mathbb{R}^{N}$, then the cross-
correlation is empirically estimated as

$\hat{R}_{N}=\frac{1}{N}\sum_{n=1}^{N}x_{n}y_{n},$ (2)

where $x_{n}\in\mathbb{R}$ and $y_{n}\in\mathbb{R}$ represent the elements of
the vector $\bf{x}$ and $\bf{y}$. Then, by the law of large numbers (LLN), the
empirical correlation converges uniformly to the true correlation $R$

$\left|\frac{1}{N}\sum_{n=1}^{N}x_{n}y_{n}-R\right|\leq\epsilon\stackrel{{\scriptstyle
N\rightarrow\infty}}{{\longrightarrow}}0,$ (3)

as depicted in Fig. 1.

Figure 1: (a) Statistical definition of cross-correlation $R$ between random
variables $x$ and $y$ under a joint distribution $p(x,y)$; (b) Empirical
cross-correlation $\hat{R}_{N}$ based on samples $(x_{n},y_{n})$; (c) Uniform
convergence of $\hat{R}_{N}$ to $R$; (d) Uniform convergence of empirical non-
linear cross-correlator $\hat{g}_{N}$ to $g(R)$; and (e) Estimation of $R$
using $g^{-1}$.

In this paper, we refer to the above estimation approach using multiplication
of samples as the empirical approach, and we explore an alternate approach
towards estimating cross-correlations using a class of non-linear functions of
samples $f:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$ such that

$\frac{1}{N}\sum_{n=1}^{N}f(x_{n},y_{n})\stackrel{{\scriptstyle
N\rightarrow\infty}}{{\longrightarrow}}\mathcal{E}[f(X,Y)]=g(R),$ (4)

where $g:\mathbb{R}\rightarrow\mathbb{R}$ is a monotonic function.

The uniform convergence of $f$ is illustrated in Fig. 1 where

$\left|\frac{1}{N}\sum_{n=1}^{N}f(x_{n},y_{n})-g(R)\right|\leq\epsilon\stackrel{{\scriptstyle
N\rightarrow\infty}}{{\longrightarrow}}0.$ (5)

The main premise of this paper is that when $x$ and $y$ are drawn from a
stationary distribution, the function $g$ is known apriori or can be estimated
with high accuracy. As a result, for a finite sample size $N$,
$g^{-1}(\frac{1}{N}\sum_{n=1}^{N}f(x_{n},y_{n}))$ is an unbiased cross-
correlation estimator, and its estimation could be closer than $\hat{R}$ to
the true cross-correlation $R$, as illustrated in Fig. 1d.

For the analysis and comparison presented in this paper, we will assume the
following without loss of generality.

  1. 1.

Both the random variables $X$ and $Y$ for which the cross-correlation is being
estimated will be assumed to be zero mean $\mathcal{E}[X]=\mathcal{E}[Y]=0$
and have unit variance $\mathcal{E}[X^{2}]=\mathcal{E}[Y^{2}]=1$. In this
case, the cross-correlation $R$ is equal to Pearson’s correlation coefficient
[23]. Note that if $X,Y$ have nonzero means $\mu_{x}$ and $\mu_{y}$, the
cross-correlation can be expressed as follow,

$\mathcal{E}[XY]=\mathcal{E}[\left(X-\mu_{x}\right)\left(Y-\mu_{y}\right)]+\mu_{x}\mu_{y}.$
(6)

Since $\mu_{x}\mu_{y}$ can be determined apriori, the accuracy of different
cross-correlators is determined by the accuracy of the cross-correlation
between the zero mean variables $X^{\prime}=X-\mu_{x}$ and
$Y^{\prime}=Y-\mu_{y}$.

  2. 2.

The non-linear function $f(x,y)$ used for generating different cross-
correlators is passive and memory-less which implies that its output depends
on the instantaneous values of $X$ and $Y$.

The paper is organized as follows: In section 2, we first propose the
mathematical framework that can be used to analyze cross-correlators for a
general class of non-linear function $f$ and for jointly Gaussian distributed
inputs. We use the framework to analyze different types of estimators, which
include the linear-rectifier cross-correlator, margin-propagation (MP) cross-
correlators [24], Huber-type [25] cross-correlators, and log-sum-exponential
(LSE) [26] estimators. In section 3, we extend the framework to arbitrary
input distributions based on the hyperdimensional mapping using the Walsh-
Hadamard transformation. In section 4, we show experiments evaluating
different correlators and the transformation method. Section 5 discusses the
advantages and disadvantages of different correlators and concludes the paper
with a brief perspective on future directions.

Figure 2: A conceptual scalar demonstration of a) empirical correlators, b)
$L_{1}$ type correlators, and c) an example of correlators using mixtures of
piece-wise linear functions with $L=1,w_{1}=1$, and $\alpha_{1}=0.5$. All
correlators $f(x,y)$ can be viewed as a difference between $h(x+y)$ and
$h(x-y)$. The upper three curves represent $h(x+y)$, the lower three curves
represent $-h(x+y)$, and curves in the middle are $f(x,y)$ for different
values of $y=y_{1},y_{2},y_{3}$.

##  2 Analysis Framework using Price’s Theorem

In this section, we present an analytical framework that can be used to
understand the behavior of different cross-correlators. A cross-correlator can
be viewed as a difference between two functions, and in Fig. 2(a), we
illustrate this for the empirical cross-correlator defined in equation (2),
which can be expressed as

$\frac{1}{N}\sum_{n=1}^{N}x_{n}y_{n}=\frac{1}{4N}\sum_{n=1}^{N}(x_{n}+y_{n})^{2}-(x_{n}-y_{n})^{2}.$
(7)

The symmetric quadratic functions $(x+y)^{2}$ and $-(x-y)^{2}$ are shown in
Fig. 2a, which results in the product $xy$ when summed together. When extended
to $N$ dimensions, the quadratic functions in equation (2) become $L_{2}$
distances $||{\bf x}+{\bf y}||_{2}^{2}$ and $||{\bf x}-{\bf y}||_{2}^{2}$ and
their difference is proportional to the empirical cross-correlation
$\frac{1}{N}\sum_{n=1}^{N}x_{n}y_{n}$. The concept can be generalized to other
norms [27]. As an example, it has been suggested that substituting the
quadratic component with linear rectifiers can yield a comparable cross-
correlator [28, 29] which uses $L_{1}$ distances $|{\bf x}+{\bf y}|$ and
$|{\bf x}-{\bf y}|$. Typically, the linear rectifier cross-correlator is
considered a reasonable approximation of the empirical cross-correlator. In
Fig. 2b, we show the equivalent construction for the $L_{1}$ type cross-
correlators using $L_{1}$ distances $|{\bf x}+{\bf y}|$ and $|{\bf x}-{\bf
y}|$.

Both $L_{2}$ and $L_{1}$ constructions can be viewed as special cases of
mixtures of piece-wise linear functions as shown in Fig. 2c and can be
expressed as

$f(x,y)=h(x+y)-h(x-y),$ (8)

where

$h(x)=\frac{1}{2}\sum_{l=1}^{L}w_{l}\left(\left|x-\alpha_{l}\right|+\left|x+\alpha_{l}\right|\right),$
(9)

with parameters $\alpha_{l}\geq 0$, $w_{l}\geq 0;\sum_{l}w_{l}=1$ and
$\left|\cdot\right|$ is an absolute-value function defined as

$\left|x\right|=\left\\{\begin{array}[]{lcc}x&;&x\geq 0\\\
-x&;&x<0.\end{array}\right.$ (10)

We now state the Lemma that can be used to compute the function
$g(R)=\mathcal{E}[f(x,y)]$.

######  Lemma 2.1.

If the function $f(x,y)$ is given by equation (8) and (9) and the joint
probability distribution of $X\in\mathbb{R}$ and $Y\in\mathbb{R}$ is given by

$p(X=x,Y=y;R)=\frac{1}{2\pi\sqrt{1-R^{2}}}\exp\left[-\frac{x^{2}+y^{2}-2Rxy}{2\left(1-R^{2}\right)}\right],$
(11)

where $R\in[-1,1]$ is the Pearson’s correlation/cross-correlation between $x$
and $y$, then $g(R)=\mathcal{E}[f(X,Y)]$ is given by

$\begin{split}g(R)=\frac{1}{2\sqrt{\pi}}\sum_{l=1}^{L}\int_{0}^{R}\frac{w_{l}}{\sqrt{(1+\rho)}}\exp\left[-\frac{\alpha_{l}^{2}}{4\left(1+\rho\right)}\right]+\frac{w_{l}}{\sqrt{(1-\rho)}}\exp\left[-\frac{\alpha_{l}^{2}}{4\left(1-\rho\right)}\right]d\rho.\end{split}$
(12)

###### Proof.

Since $f$ is a memory-less function with a well-defined Fourier transform and
$X$ and $Y$ are zero-mean, unit variance, jointly distributed Gaussian random
variables, we can apply Price’s theorem [30, 29, 31] which states that

$\frac{\partial\mathcal{E}(f)}{\partial
R}=\mathcal{E}\left[\frac{\partial^{2}f}{\partial x\partial y}\right],$ (13)

where the expectation operator $\mathcal{E}$ is defined as

$\mathcal{E}\left(f\right)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)p(x,y;R)dxdy=g(R).$
(14)

The partial derivatives of the sub-functions $h(x+y),h(x-y)$ in (9) are

$\frac{\partial^{2}h(x+y)}{\partial x\partial
y}=\frac{1}{2}\sum_{l=1}^{L}w_{l}\left[\delta\left(x+y-\alpha_{l}\right)+\delta\left(x+y+\alpha_{l}\right)\right],$
(15a) $-\frac{\partial^{2}h(x-y)}{\partial x\partial
y}=\frac{1}{2}\sum_{l=1}^{L}w_{l}\left[\delta\left(x-y-\alpha_{l}\right)+\delta\left(x-y+\alpha_{l}\right)\right],$
(15b)

where $\delta(.)$ denotes the Dirac-delta function. The expectation operators
can be computed as follow

$\displaystyle\mathcal{E}\left[\frac{\partial^{2}h(x+y)}{\partial x\partial
y}\right]$
$\displaystyle=\frac{1}{2}\sum_{l=1}^{L}w_{l}\int_{-\infty}^{\infty}p(x,-x+\alpha_{l})+p(x,-x-\alpha_{l})dx$
$\displaystyle=\frac{1}{\sqrt{\pi(1+R)}}\sum_{l=1}^{L}w_{l}\exp\left[-\frac{\alpha_{l}^{2}}{4\left(1+R\right)}\right]$
(16a) $\displaystyle-\mathcal{E}\left[\frac{\partial^{2}h(x-y)}{\partial
x\partial y}\right]$
$\displaystyle=\frac{1}{2}\sum_{l=1}^{L}w_{l}\int_{-\infty}^{\infty}p(x,x-\alpha_{l})+p(x,x+\alpha_{l})dx$
$\displaystyle=\frac{1}{\sqrt{\pi(1-R)}}\sum_{l=1}^{L}w_{l}\exp\left[-\frac{\alpha_{l}^{2}}{4\left(1-R\right)}\right].$
(16b)

Substituting the results into (8) and (13) will get

$\begin{split}\frac{\partial g}{\partial
R}=\frac{1}{\sqrt{\pi(1+R)}}\sum_{l=1}^{L}w_{l}\exp\left[-\frac{\alpha_{l}^{2}}{4\left(1+R\right)}\right]+\frac{1}{\sqrt{\pi(1-R)}}\sum_{l=1}^{L}w_{l}\exp\left[-\frac{\alpha_{l}^{2}}{4\left(1-R\right)}\right],\end{split}$
(17)

which leads to the expression for $g(R)$,

$\begin{split}g(R)=\int_{0}^{R}\frac{1}{\sqrt{\pi(1+\rho)}}\sum_{l=1}^{L}w_{l}\exp\left[-\frac{\alpha_{l}^{2}}{4\left(1+\rho\right)}\right]d\rho+\int_{0}^{R}\frac{1}{\sqrt{\pi(1-\rho)}}\sum_{l=1}^{L}w_{l}\exp\left[-\frac{\alpha_{l}^{2}}{4\left(1-\rho\right)}\right]d\rho.\end{split}$
(18)

∎

Example 1: When $L=1,w_{1}=1,\alpha_{1}=0$, the function $f$ is reduced to

$f(x,y)=|x+y|-|x-y|,$ (19)

which is the well-studied linear rectifier correlator. In this case, the
relation (12) can be evaluated in closed form and is given by

$g_{L1}(R)=\frac{2}{\sqrt{\pi}}(\sqrt{1+R}-\sqrt{1-R}).$ (20)

Example 2: When $w_{l}=1/L,\alpha_{l}=c/L,l=1,..,L$ and
$c,L\rightarrow\infty$, the function $f$ is reduced to

$f(x,y)=\frac{1}{2c}(\left(x+y\right)^{2}-\left(x-y\right)^{2}),$ (21)

where $c>|x|$ is the range of inputs. So the function $f$ becomes the
empirical correlator. In this case, the summation in the relation (12) can be
replaced by integrals in the limit $L\rightarrow\infty$ in which case

$\displaystyle\frac{\partial g_{L2}}{\partial
R}=\frac{1}{c\sqrt{\pi(1+R)}}\int_{0}^{\infty}\exp\left[-\frac{x^{2}}{4\left(1+R\right)}\right]dx+\frac{1}{c\sqrt{\pi(1-R)}}\int_{0}^{\infty}\exp\left[-\frac{x^{2}}{4\left(1-R\right)}\right]dx=\frac{2}{c}.$
(22)

Therefore, $g_{L2}(R)=\frac{2}{c}R$ matches the result for a scaled empirical
cross-correlation.

Figure 3: a):$\frac{\partial g}{\partial R}$ for $L=1,w_{1}=1$, and
$\alpha_{1}=\\{0,0.2,0.4,0.6,0.8,1.0\\}$; b) The $g(R)$ of the empirical and
linear rectifier correlators and the $g(R)$ derived by integrating
$\frac{\partial g}{\partial R}$ for $L=1,w_{1}=1$, and $\alpha_{1}=0.5,0.725$
and 1. It demonstrates how the analysis framework can be used to construct
varying correlators with different output characteristics.

The above examples show how specific cross-correlators can be constructed
using Lemma 2.1. Even though the equation (12) may not be solved in closed-
form, it is analytical and hence can be used to visualize the form of $g(R)$
for specific choices of $\alpha_{l}$ and $w_{l}$. Fig. 3a visualizes the
partial derivative $\frac{\partial g}{\partial R}$ in (17) when $L=1,w_{1}=1$,
and $\alpha_{1}=\\{0,0.2,0.4,0.6,0.8,1\\}$. The zoomed-in window demonstrates
that the partial derivatives stay almost constant around zero correlation,
which implies the linearity of $g(R)$ around zero correlation for any cross-
correlators. Fig. 3b displays the $g(R)$ of both empirical and linear
rectifier correlators and the $g(R)$ for different values of a single offset
$\alpha_{1}=0.5,0.725,1$ computed by integration. The results will be verified
by the Monte Carlo experiment for MP correlators in section 4.

Let’s denote $\mathcal{E}\left(f\right)$ or $g(R)$ by $y$, we can derive the
following from $g_{L1}(R)$ in (20)

$R^{2}=\frac{\pi}{4}y^{2}-\frac{\pi^{2}}{64}y^{4}.$ (23)

This suggests that $g^{-1}$ in Fig. 1 can be robustly estimated using a
polynomial expansion with a relatively low degree. This is important since the
closed-form solution for equation (18) may not be computed for different
choices of $w_{l},\alpha_{l}$. As a result, $g^{-1}$ has to be
learned/estimated by drawing samples with known apriori cross-correlation,
which is then used to estimate $R$ according to Fig. 1. As we will show later
in section 3, this calibration procedure and procedure to estimate $R$ can be
agnostic to the input distribution.

We now apply the calibration procedure to three other types of functions of
the type given by expression (8). The first is the margin-propagation (MP)
function given by

$\displaystyle(x-z)_{+}+(-x-z)_{+}$ $\displaystyle=\gamma,$ (24)
$\displaystyle h(x)$ $\displaystyle=z,$ (25)

where $(\cdot)_{+}$ is the ReLU function, and $\gamma\geq 0$ is the hyper-
parameter. It can be easily verified that the MP function in (24) is
equivalent to (9) with $L=1$ and $\alpha_{1}=0.5\gamma$.

The second function is the Huber function which requires an infinite number of
splines and is given by

$h(x)=\begin{cases}&0.5x^{2}/\delta,\enskip x<\delta\\\
&|x|-\frac{1}{2}\delta,\enskip x\geq\delta,\end{cases}$ (26)

where $\delta>0$ is a threshold parameter.

The third function is a log-sum-exp (LSE) function which also requires an
infinite number of splines and is given by

$h(x)=\frac{1}{a}(\log(\exp\left[ax\right]+\exp\left[-ax\right]),$ (27)

where $a>0$ is a scaling factor.

As we will show in section 4, the Huber and LSE correlators fall between the
empirical and linear rectifier correlators, hence the inverse cross-
correlation function $g^{-1}$ can be approximated by a polynomial of degree
lower or equal to the degree needed for $g_{L1}^{-1}$. In practice, we found
that fourth order polynomial is sufficient for calibration of $g_{L1}^{-1}$.
For calibrating MP correlators, higher degree polynomials were needed as the
value of $\gamma$ increases as its $|\frac{\partial g}{\partial R}|$ is not
monotonically increasing as correlation increases. As such, a fifth-order
polynomial was used to learn the inverse cross-correlation function $g^{-1}$.

##  3 Extension to non-Gaussian distributions

The theoretical results presented in section 2 assumed random variables with
joint Gaussian distributions. In this section, we extend the previous results
for non-Gaussian distributions. To achieve this, we use results from the
hyperdimensional computing literature, which state that variances and cross-
correlations are preserved when random variables are mapped into high-
dimensional space using unitary random matrices.

######  Lemma 3.1.

Let $X$ and $Y$ be zero-mean random variables with unit variance and with a
cross-correlation $R$. Let $\Phi:\mathbb{R}^{N}\rightarrow\mathbb{R}^{M}$
denote a high-dimensional embedding using a Unitary transform such that
$\mathcal{E}[\Phi({\bf x})]=\mathbf{0}$. Then, as $N\rightarrow\infty$,
$\frac{1}{N}\left<\Phi({\bf x}),\Phi({\bf y})\right>\rightarrow R$, where
$\left<\cdot,\cdot\right>$ is the inner product.

###### Proof.

Suppose ${\bf x}$ and ${\bf y}$ are $\mathbb{R}^{N}$-valued random vector, and
each entry $x_{n}$, $y_{n}$ are independently identically distributed (i.i.d)
variables with the joint probability density function $p(X=x_{n},Y=y_{n};R)$

$\displaystyle\frac{1}{N}\left<\Phi({\bf x}),\Phi({\bf
y})\right>=\frac{1}{N}\left<{\bf x},{\bf
y}\right>=\frac{1}{N}\sum_{n=1}^{N}x_{n}y_{n}\stackrel{{\scriptstyle
N\rightarrow\infty}}{{\longrightarrow}}R.$ (28)

∎

The Walsh-Hadamard-Transform is one such unitary transform
$\mathcal{H}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ and it can be
represented by a $N\times N$ Hadamard matrix. The transformation generally
requires the input zero-padded to a power of two. An example $4\times 4$ WHT
matrix is shown below

$H_{4}=\frac{1}{2}\begin{pmatrix}1&1&1&1\\\ 1&-1&1&-1\\\ 1&1&-1&-1\\\
1&-1&-1&1\end{pmatrix}.$ (29)

The Hadamard matrices are orthogonal and symmetric matrices composed of +1 and
-1 with a normalization factor $1/\sqrt{N}$, which makes it easy for
implementations and computations. Besides keeping the covariance between
random variables, it can also be shown that the transformed zero-mean
variables converge to the joint Gaussian distribution with the same variance
and covariance.

######  Lemma 3.2.

Let $X$ and $Y$ be zero-mean random variables with finite variances and the
cross-correlation $R$. Let
$\mathcal{H}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ denote the Walsh-
Hadamard-Transform. Suppose the entries of the vector equation
$\bf{x^{\prime}}=\mathcal{H}(\bf{x})$ are given by
$x_{n}^{\prime}=h_{n}(x_{1},...,x_{n})$, and $h_{n}(\cdot)$ is therefore

$h_{n}(x_{1},...,x_{n})=\begin{cases}&\frac{1}{\sqrt{N}}\sum_{n=1}^{N}x_{n},\enskip\enskip\enskip
n=1,\\\
&\frac{1}{\sqrt{N}}(\sum_{n=1}^{N/2}x_{n}-\sum_{n=1}^{N/2}x_{n}),\enskip n\neq
1,\end{cases}$ (30)

where the $x_{n}$ are independently identically distributed (i.i.d.) samples
from $X$. Then, as $N\rightarrow\infty$, the joint probability distribution
$p(x_{n}^{\prime},y_{n}^{\prime})$ converges to a bivariate Gaussian
distribution with zero-mean, same variances $\sigma^{2}_{x}$,
$\sigma^{2}_{y}$, and covariance $R$.

###### Proof.

Suppose $X$ and $Y$ are zero-mean random variables with finite variances and
covariance $R$. According to the multivariate Central Limit Theorem (CLT)
[32], as $N\rightarrow\infty$ the joint distribution
$p(\sqrt{N}\bar{x}_{N},\sqrt{N}\bar{y}_{N})$ converges to bivariate Gaussian
distribution with zero mean and the same variances and covariance R, where
$\bar{x}_{N}=\frac{1}{N}\sum_{n=1}^{N}x_{n}$ is the average of $N$
independently identically distributed samples of $X$.

Notice that the transformed entry after WHT
$x_{n}^{\prime}=h_{n}(x_{1},...,x_{n})$ can be expressed as

$x_{n}^{\prime}=\begin{cases}&\sqrt{N}\bar{x}_{N},\enskip\enskip\enskip
n=1,\\\
&\frac{1}{\sqrt{2}}\sqrt{\frac{N}{2}}(\bar{x}_{N/2}-\bar{y}_{N/2}),\enskip
n\neq 1.\end{cases}$ (31)

For $n=1$, the CLT can be applied to the transformed input $x_{1}^{\prime}$
and $y_{1}^{\prime}$. For the case of $n\neq 1$, note that the CLT also
applies to $\pm\sqrt{\frac{N}{2}}\bar{x}_{N/2}$ and
$\pm\sqrt{\frac{N}{2}}\bar{y}_{N/2}$, so they become bivariate Gaussian with
zero means, and same variances and covariance, and so is their sum divided by
$\sqrt{2}$. ∎

As such, using the proof in section 2, it can be shown that for correlators
that can be expressed by equations (8) and (9), the expected output
$\mathcal{E}[f(\mathcal{H}({\bf x}),\mathcal{H}({\bf y}))]$ is equal to
$g(R)$. In other words, for non-Gaussian distributed variables with zeros mean
and finite covariance $R$, we can first transform the inputs to jointly
Gaussian distribution, which preserves the covariance, and then use the cross-
correlator in section 2 to estimate the cross-correlation $R$ using the
transformed data and the same $g^{-1}$.

##  4 Experiments Results and Analysis

In the first part of the experiments, we validate the analytical results in
section 2 when the random variables are drawn from joint Gaussian
distribution. For this specific case, a closed-form expression of the
Cramér–Rao bound can be computed. This bound can then be used to evaluate the
effectiveness of any unbiased cross-correlation estimators. We will then
validate our analytical results for arbitrary input distributions.

###  4.1 Cross-correlation Cramér–Rao bound

According to the Cramér–Rao bound (CRB) [33], the variance or standard
deviation of any unbiased estimator (satisfying regularity conditions) is
lower bounded by the reciprocal of the Fisher information [34]. The following
lemma presents the Cramér–Rao bound of the correlation estimator under the
standard bivariate normal distribution.

######  Lemma 4.1.

For a collection of $N$ independent and identically distributed (iid)
bivariate variables $X\in\mathbb{R}$ and $Y\in\mathbb{R}$ drawn from the
jointly Gaussian distribution defined by (11) the Cramér–Rao bound of the
correlation estimator is given by

$\sigma^{2}(\hat{R})\geq\frac{(1-R^{2})^{2}}{N(1+R^{2})}.$ (32)

###### Proof.

For the distribution (11), the Fisher information for the correlation
coefficient of a single pair of samples $(x,y)$ can be computed as,

$I(R)=-\mathcal{E}\left[\frac{\partial^{2}(l(x,y;R))}{\partial R^{2}}\right],$
(33)

where $l(x,y;R)$ represents the natural logarithm likelihood function of a
single sample pair $(x,y)$ for the joint density function $p(X=x,Y=y;R)$

$l(x,y;R)=-\log(2\pi)-\log(\sqrt{1-R^{2}})-\frac{x^{2}+y^{2}-2Rxy}{2\left(1-R^{2}\right)}.$
(34)

From (34) and using $\mathcal{E}[X^{2}]=\mathcal{E}[Y^{2}]=1$ and
$\mathcal{E}[XY]=R$ leads to the Fisher’s information metric

$I(R)=\frac{1+R^{2}}{(1-R^{2})^{2}}.$ (35)

For $N$ independent and identically distributed (iid) samples, the total
Fisher information is the sum of information from each individual sample.
Therefore, the Cramér–Rao bound is given by

$\sigma^{2}(\hat{R})\geq\frac{1}{I(R)}=\frac{(1-R^{2})^{2}}{N(1+R^{2})}.$ (36)

∎

In other words, the accuracy of all cross-correlators is upper bounded by the
Fisher information. In section 4, we use the standard deviation of estimation
errors $\sigma$ to assess the performance of cross-correlators.

###  4.2 Monte-Carlo Experiments for Jointly Gaussian Inputs

This section presents the results using different correlators to estimate the
covariance for standard jointly normal distribution. To study and compare
their performance, vectors of different lengths are randomly sampled from the
zero mean and unit variance Gaussian distribution. Each vector pair
$\bf{s_{1}},\bf{s_{2}}\in\mathbb{R}^{N}$ is mixed in the following way to
generate a bivariate Gaussian distribution $\bf{x}=(\bf{x_{1}},\bf{x_{2}})$
with different correlation $R$ to learn and test the inverse cross-correlation
function $g(R)$,

$\displaystyle\bf{x_{1}}$ $\displaystyle=\bf{s_{1}},$ (37)
$\displaystyle\bf{x_{2}}$
$\displaystyle=R\bf{s_{1}}+\sqrt{1-R^{2}}\bf{s_{2}},$ (38)
$\displaystyle\bf{x}$ $\displaystyle=(\bf{x_{1}},\bf{x_{2}})\sim
N(\vec{0},\Sigma),\Sigma=\begin{pmatrix}1&R\\\ R&1\end{pmatrix}.$ (39)

Figure 4: The estimated $\hat{g}(R)$ for the MP, Huber, and LSE functions for
the standard bivariate normal distribution from Monte Carlo experiments. The
expected outputs of the empirical and linear rectifier correlators are added
for comparison. (a) The $\hat{g}(R)$ of MP correlators with
$\gamma=\\{0,1,1.45,2\\}$; (b) The $\hat{g}(R)$ of Huber correlators with
$\delta=\\{0.5,1.4,2\\}$; (c) The $\hat{g}(R)$ of LSE correlators with
$a=\\{0.1,1,3\\}$.

Fig. 4 shows the average output of cross-correlation function $\hat{g}(R)$
corresponding to the MP, Huber, and LSE functions with different parameters,
which is used as an approximation of the correlation function $g(R)$ to get
the calibration function $g^{-1}$. The output of the empirical correlator is
shown for comparison. Note that the linear rectifier correlator is a special
case for the MP, Huber, and LSE functions, which is included and labeled as
"MP $\gamma=0$," "Huber $\delta=0$," and "LSE $a\rightarrow\infty$." In fact,
the normalized $g(R)$ for the Huber function and LSE functions are bounded
above and below by the normalized $g_{L1}(R)$ and $g_{L2}(R)$. This is easy to
see for Huber functions, as it’s a combination of the quadratic function
($L_{2}$) and absolute value function ($L_{1}$). For LSE functions, as the
scaling factor $a$ increases, the $h(x)$ for LSE functions in expression (27)
can be simplified to $\frac{1}{a}log(exp\left[a|x|\right]=|x|$ as the negative
part will go to zero exponentially fast. On the other hand, as $a$ decreases,
we can Taylor expand the exponential function at zero, which leads to

$\displaystyle
h(x)=\frac{1}{a}\log(\sum_{n=0}^{\infty}\frac{(ax)^{n}}{n!}+\frac{(-ax)^{n}}{n!}).$
(40)

The odd-degree terms cancel each other out and higher-order terms decay fast,
which leads to

$\displaystyle h(x)\approx\frac{1}{a}\log(2+(ax)^{2}).$ (41)

Applying the same trick to the logarithm but expanding at 2, we have

$\displaystyle
h(x)\approx\frac{1}{a}(\log(2)+\frac{1}{2}(ax)^{2}+...)\approx\frac{1}{a}+\frac{1}{2}ax^{2}.$
(42)

Therefore, the LSE correlator approaches the empirical correlator as $a$
decreases.

Figure 5: The standard deviation of the cross-correlation estimation error of
MP, Huber, and LSE correlators for different correlation levels and the CRB at
the dimension of 256, with the SNR of empirical and linear rectifier added for
reference. (a) The error plot of MP correlators with $\gamma=\\{1,1.45,2\\}$;
(b) The error plot of Huber correlators with $\sigma=\\{0.5,1.4,2\\}$; (c) The
error plot of LSE correlators with $a=\\{0.1,1,3\\}$. Figure 6: The SNR plots
of MP, Huber, and LSE correlators from correlator length of 16 to 65536, with
the SNR of empirical and linear rectifier correlators as a reference. a) The
SNR plot of MP correlators with $\gamma=\\{1,1.45,2\\}$; b) The SNR plot of
Huber correlators with $\sigma=\\{0.5,1.4,2\\}$; c) The SNR plot of LSE
correlators with $a=\\{0.1,1,3\\}$.

As discussed in section 4.1, the standard deviation of estimation can be used
to evaluate the performance of correlators. Since the expected value of
$g^{-1}$ is equal to $R$, the standard deviation of estimation
$\sigma(\bar{R})$ is equivalent to the standard deviation of estimation error
$\sigma(R-g^{-1}(f))$. In Fig. 5, we display the standard deviation of cross-
correlation estimation errors and its CRB lower bound, denoted by $\sigma$,
made by the MP, Huber, and LSE correlators with different parameters at
different levels of cross-correlation using the learned inverse cross-
correlation function $g^{-1}$. It is observed that the linear rectifier
correlator is more accurate when the signal of interest is highly correlated,
while the empirical correlator makes less error in the other case. The error
profile of the Huber and LSE correlators becomes more similar to the empirical
correlator when $\sigma$ is high and $a$ is small and approaches the linear
rectifier correlator otherwise. The MP correlator in (24) is equivalent to the
linear rectifier correlator when $\gamma=0$. As $\gamma$ increases, the
performance degrades and performs even worse than the empirical correlator.

In Fig. 6, we plot the signal-to-noise-ratio (SNR) for these cross-correlators
for different correlator length $N$, where the SNR is defined as:

$\text{SNR}=20\log_{10}(\frac{1}{\bar{\sigma}}),$ (43)

where $\bar{\sigma}$ is the average standard deviation of the error across
multiple Monte-Carlo runs. The result shows that the linear rectifier
correlator has the best SNR in joint Gaussian input distribution. Also, the
SNR increases by 3dB when the correlation length $N$ is doubled, which can be
attributed to the reduction in the estimation error due to simple averaging.

Figure 7: Scatter plots showing the random variables $x$ and $y$ drawn from
different probability distributions with zero-mean, unit variance and a
correlation R=0.8. The data points denoted by the blue cross are drawn from
the original distribution and the data points denoted by the red circles are
the inputs after applying WHT. After applying WHT, the resulting distribution
converges towards a standard bivariate normal distribution with a correlation
$R$=0.8.

###  4.3 Monte-Carlo Experiments with WHT and Non-Gaussian Inputs

In this section, the WHT method in section 3 was tested on non-Gaussian
distributions to verify that the function $g(R)$ remains unchanged for non-
Gaussian inputs after being transformed by the WHT. Fig. 7 shows varying joint
input distributions before and after the WHT, which are denoted by blue
crosses and red circles, respectively. The non-Gaussian distributions before
the transformation are combinations of the uniform, Gaussian, and gamma
distributions, and the inputs have a correlation of 0.8. It can be seen that
the vectors become jointly Gaussian distributed after the transformation, and
their correlation is retained. Monte-Carlo experiments show that the expected
correlator output $\mathcal{E}[f(\mathcal{H}(x),\mathcal{H}(y))]$ is the same
as the $g(R)$ for jointly Gaussian inputs. So no calibration is needed to
learn $g^{-1}$ for different distributions.

On the other hand, it should be noticed that the standard deviation for
different cross-correlation estimators is not guaranteed to be the same for
different input distributions, even after the WHT transformation. To see this,
notice that the WHT process will not change the output for the empirical
correlator because the WHT transformation is unitary. However, the standard
deviation of the empirical correlator, which is given by

$\begin{split}\mathcal{E}(XY-\mathcal{E}[XY])^{2}=\mathcal{E}[X^{2}Y^{2}]-g^{2}(R)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}x^{2}y^{2}\enskip
p(x,y)dxdy-g^{2}(R),\end{split}$ (44)

will change as the joint probability density function $p(x,y)$ changes. In
particular, for joint Gaussian distribution, it can be verified that
$\mathcal{E}(xy-\mathcal{E}[xy])^{2}=1+R^{2}$, whereas it is
$1-\frac{1}{5}R^{2}$ in the case of jointly uniform distribution. As an
example, Fig. 8 shows the standard deviation of estimation error plots for the
joint uniform distribution as shown in Fig. 7a. It is observed that the
contour of the error plot changes for all correlators except for the linear
rectifier correlator, which is still the best-performing cross-correlation
estimator in this specific non-Gaussian input distribution.

Figure 8: The standard deviation error plot of MP, Huber, and LSE correlators
for the symmetric input distribution test when dimension is 256. a) The error
plot of MP correlators with $\gamma=\\{1,1.45,2\\}$; b) The error plot of
Huber correlators with $\sigma=0.5,1.4,2\\}$; c) The error plot of LSE
correlators with $a=\\{0.1,1,3\\}$.

##  5 Discussions and Conclusions

In this paper, we present a mathematical framework for analyzing different
types of cross-correlators. The closed-form solution facilitates the
comparison of different cross-correlators and allows us to understand
performance trade-offs. The analysis framework has been verified by Monte
Carlo simulation for different input distributions. The error analysis reveals
that the shape of the error profile exhibits a trade-off. Cross-correlation
estimators that exhibit high errors near $R\approx 0$ make fewer errors near
$|R|\approx 1$. The hyperparameters of the Huber estimators, MP estimators,
and LSE estimators can be adapted to achieve different error profiles.

However, the complexity of implementing these different online estimators on
hardware could be significantly different. The empirical and Huber cross-
correlator relies on the quadratic function and hence may be difficult to
implement on hardware. On the other hand, the MP and LSE correlators can be
easily implemented on analog hardware [35]. Another benefit of employing other
correlators is the potential for improved computational efficiency and wider
dynamic range when implementing them on digital systems. It is evident that
the linear rectifier correlator and the MP correlator (without additional
offsets) are more cost-effective than the empirical correlator since each
multiplication is replaced by three and five addition operations [36]. The
computation of the linear rectifier can be further simplified to condition and
shift operations, and the addition operation is resistant to underflow on a
fixed-point system. Moreover, the LSE correlator exhibits superior numerical
stability at the expense of computational complexity, a common technique
employed in machine learning to address the issue of gradient updates.

If an accurate estimate of the covariance or Pearson’s correlation is
necessary, then the calibration is required to learn the inverse function
$g^{-1}$ for inputs with diverse distributions and variances. To this end, the
WHT is proposed as a pre-processing technique to convert input distributions
to joint normal distributions, allowing the calibration process to be
indifferent to input distributions. The limitation is that the transformation
requires the dimension of inputs to be a power of two. However, calibration is
still needed for varying variances, which inevitably introduces extra
computational and resource expenses. On the other hand, the raw correlator
output $f(\cdot)$ would be sufficient if the objective of the task is to
obtain only a similarity score for different signals. For instance, in
hyperdimensional computing, whether some $\bf{y}\in\mathbb{R}^{N}$ is
contained in a set $\mathcal{S}\in\mathbb{R}^{N}$ is checked by if the dot
product $\left<\phi(y),\phi(\mathcal{S})\right>$ is above a certain threshold,
where $\phi(\cdot)$ is the hyperdimensional representation. The WHT can be
incorporated into the hyperdimensional mapping process $\phi$, and alternative
correlation functions, such as those discussed, can be used instead of
resource-intensive and computationally inefficient inner products.

Regarding the accuracy for estimating the cross-correlation of jointly
Gaussian distributed inputs, it appears that the empirical method may not be
the most effective correlator. The linear rectifier correlator is superior in
estimating the covariance for highly correlated signals and in terms of
overall SNR. Its standard deviation of error plot has a similar trend with the
Cramér–Rao bound (CRB). The performance gained for high correlation can be
explained by the shape of $g(R)$. The variance of output $f(x,y)$ is
relatively small with respect to the gradient of $g(R)$ for high correlation,
which results in a larger confidence interval for estimations. The Huber and
LSE correlator’s performance is bounded by the linear rectifier and empirical
correlator. The MP correlator in (24) is equivalent to the linear rectifier
correlator when $\gamma=0$. As $\gamma$ increases, its performance for higher
$\gamma$ values can be potentially improved by introducing offsets. Of course,
the above observations are not guaranteed to hold for other input
distributions and are left for future research.

## References

  * [1] D.S. Moore, G.P. McCabe, and B.A. Craig.  Introduction to the Practice of Statistics.  W. H. Freeman, 2014.
  * [2] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [3] K.I. Diamantaras and Sun-Yuan Kung.  Cross-correlation neural network models.  IEEE Transactions on Signal Processing, 42(11):3218–3223, 1994.
  * [4] M. Betke and N.C. Makris.  Fast object recognition in noisy images using simulated annealing.  In Proceedings of IEEE International Conference on Computer Vision, pages 523–530, 1995.
  * [5] Richard O. Duda and Peter E. Hart.  Pattern classification and scene analysis.  In A Wiley-Interscience publication, 1974.
  * [6] Jae-Chern Yoo and Tae Hee Han.  Fast normalized cross-correlation.  Circuits, systems and signal processing, 28:819–843, 2009.
  * [7] S. Adrián-Martínez, M. Ardid, M. Bou-Cabo, I. Felis, C. Llorens, J. A. Martínez-Mora, and M. Saldaña.  Acoustic signal detection through the cross-correlation method in experiments with different signal to noise ratio and reverberation conditions, 2015.
  * [8] Hsueh-Jyh Li, Yung-Deh Wang, and Long-Huai Wang.  Matching score properties between range profiles of high-resolution radar targets.  IEEE Transactions on Antennas and Propagation, 44(4):444–452, 1996.
  * [9] H.-J. Li and S.-H. Yang.  Using range profiles as feature vectors to identify aerospace objects.  IEEE Transactions on Antennas and Propagation, 41(3):261–268, 1993.
  * [10] Gilbert Strang.  Linear algebra and its applications.  Thomson, Brooks/Cole, Belmont, CA, 2006.
  * [11] Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg S Corrado, and Jeff Dean.  Distributed representations of words and phrases and their compositionality.  In C.J. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 26. Curran Associates, Inc., 2013.
  * [12] Chen Wang, Wenshan Wang, Yuheng Qiu, Yafei Hu, and Sebastian A. Scherer.  Visual memorability for robotic interestingness via unsupervised online learning.  CoRR, abs/2005.08829, 2020.
  * [13] Chen Wang, Le Zhang, Lihua Xie, and Junsong Yuan.  Kernel cross-correlator.  Proceedings of the AAAI Conference on Artificial Intelligence, 32(1), Apr. 2018.
  * [14] Zitong Yu, Xiaobai Li, Pichao Wang, and Guoying Zhao.  Transrppg: Remote photoplethysmography transformer for 3d mask face presentation attack detection.  CoRR, abs/2104.07419, 2021.
  * [15] Elhanan Elboher and Michael Werman.  Asymmetric correlation: A noise robust similarity measure for template matching.  IEEE Transactions on Image Processing, 22(8):3062–3073, 2013.
  * [16] J.P. Lewis.  Fast normalized cross-correlation.  Ind. Light Magic, 10, 10 2001.
  * [17] Pentti Kanerva.  Hyperdimensional computing: An introduction to computing in distributed representation with high-dimensional random vectors.  Cognitive Computation, 1:139–159, 2009.
  * [18] Anthony Thomas, Sanjoy Dasgupta, and Tajana Rosing.  A theoretical perspective on hyperdimensional computing.  J. Artif. Int. Res., 72:215–249, jan 2022.
  * [19] Mohamad H. Hassoun.  Associative neural memories : theory and implementation.  1993\.
  * [20] Trenton Bricken and Cengiz Pehlevan.  Attention approximates sparse distributed memory.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 15301–15315. Curran Associates, Inc., 2021.
  * [21] Mohsen Imani, Zhuowen Zou, Samuel Bosch, Sanjay Anantha Rao, Sahand Salamat, Venkatesh Kumar, Yeseong Kim, and Tajana Rosing.  Revisiting hyperdimensional learning for fpga and low-power architectures.  In 2021 IEEE International Symposium on High-Performance Computer Architecture (HPCA), pages 221–234, Feb 2021.
  * [22] Alejandro Hernández-Cano, Namiko Matsumoto, Eric Ping, and Mohsen Imani.  Onlinehd: Robust, efficient, and single-pass online learning using hyperdimensional system.  In 2021 Design, Automation And Test in Europe Conference and Exhibition (DATE), pages 56–61, 2021.
  * [23] Karl Pearson.  Note on regression and inheritance in the case of two parents.  Proceedings of the Royal Society of London, 58:240–242, 1895.
  * [24] Shantanu Chakrabartty and Gert Cauwenberghs.  Margin propagation and forward decoding in analog VLSI.  2003\.
  * [25] Peter J. Huber.  Robust Estimation of a Location Parameter.  The Annals of Mathematical Statistics, 35(1):73 – 101, 1964.
  * [26] Giuseppe C. Calafiore, Stephane Gaubert, and Corrado Possieri.  A universal approximation result for difference of log-sum-exp neural networks.  IEEE Transactions on Neural Networks and Learning Systems, 31(12):5603–5612, 2020.
  * [27] Hanting Chen, Yunhe Wang, Chunjing Xu, Boxin Shi, Chao Xu, Qi Tian, and Chang Xu.  Addernet: Do we really need multiplications in deep learning?  CVPR, 2020.
  * [28] Jr. Faran, James J.  A simple electronic correlator.  The Journal of the Acoustical Society of America, 24(4):452–452, 06 2005.
  * [29] E. McMahon.  An extension of price’s theorem (corresp.).  IEEE Transactions on Information Theory, 10(2):168–168, 1964.
  * [30] R. Price.  A useful theorem for nonlinear devices having gaussian inputs.  IRE Transactions on Information Theory, 4(2):69–72, 1958.
  * [31] A. Papoulis.  Comments on ’an extension of price’s theorem’ by mcmahon, e. l.  IEEE Transactions on Information Theory, 11(1):154–154, 1965.
  * [32] Achim Klenke.  Probability Theory: A Comprehensive Course.  Springer, 2007.
  * [33] C. adhakrishna Rao.  Information and accuracy attainable in the estimation of statistical parameters.  Bulletin of the Calcutta Mathematical Society, 37(3):81–91, 1945.
  * [34] Rory A. Fisher.  On the mathematical foundations of theoretical statistics.  Philosophical Transactions of the Royal Society A, 222:309–368.
  * [35] Ming Gu and Shantanu Chakrabartty.  Synthesis of bias-scalable cmos analog computational circuits using margin propagation.  IEEE Transactions on Circuits and Systems I: Regular Papers, 59(2):243–254, 2012.
  * [36] Abhishek Ramdas Nair, Pallab Kumar Nath, Shantanu Chakrabartty, and Chetan Singh Thakur.  Multiplierless mp-kernel machine for energy-efficient edge devices.  IEEE Transactions on Very Large Scale Integration (VLSI) Systems, 30(11):1601–1614, 2022.