# Permutation Invariant Representations with Applications to Graph Deep
Learning

Radu Balan  ,  Naveed Haghani  and  Maneesh Singh  Department of Mathematics,
University of Maryland, College Park, MD 20742, USA<EMAIL_ADDRESS>Applied
Mathematics and Statistics and Scientific Computation Program, University of
Maryland, College Park, MD 20742, USA<EMAIL_ADDRESS>Verisk Analytics,
Jersey City, NJ 07310, USA<EMAIL_ADDRESS>
###### Abstract.

This paper presents primarily two Euclidean embeddings of the quotient space
generated by matrices that are identified modulo arbitrary row permutations.
The original application is in deep learning on graphs where the learning task
is invariant to node relabeling. Two embedding schemes are introduced, one
based on sorting and the other based on algebras of multivariate polynomials.
While both embeddings exhibit a computational complexity exponential in
problem size, the sorting based embedding is globally bi-Lipschitz and admits
a low dimensional target space. Additionally, an almost everywhere injective
scheme can be implemented with minimal redundancy and low computational cost.
In turn, this proves that almost any classifier can be implemented with an
arbitrary small loss of performance. Numerical experiments are carried out on
two data sets, a chemical compound data set (QM9) and a proteins data set
(PROTEINS_FULL).

##  1\. Introduction

This paper is motivated by a class of problems in graph deep learning, where
the primary task is either graph classification or graph regression. In either
case, the result should be invariant to arbitrary permutations of graph nodes.

As we explain below, the mathematical problem analyzed in this paper is a
special case of the permutation invariance issue described above. To set the
notations consider the vector space ${\mathbb{R}}^{n\times d}$ of $n\times d$
matrices endowed with the Frobenius norm
${\|X\|}=\left(trace(XX^{T})\right)^{1/2}$ and its associated Hilbert-Schmidt
scalar product, ${\langle X,Y\rangle}=trace(XY^{T})$. Let ${\mathcal{S}}_{n}$
denote the symmetric group of $n\times n$ permutation matrices.
${\mathcal{S}}_{n}$ is a finite group of size $|{\mathcal{S}}_{n}|=n!$.

On ${\mathbb{R}}^{n\times d}$ we consider the equivalence relation $\sim$
induced by the symmetric group of permutation matrices ${\mathcal{S}}_{n}$ as
follows. Let $X,Y\in{\mathbb{R}}^{n\times d}$. Then we say $X\sim Y$ if there
is $P\in{\mathcal{S}}_{n}$ so that $Y=PX$. In other words, two matrices are
equivalent if one is a row permutation of the other. The equivalence relation
induces a natural distance on the quotient space
${\widehat{{\mathbb{R}}^{n\times d}}}:={\mathbb{R}}^{n\times d}/\sim$,

(1.1) $d:{\widehat{{\mathbb{R}}^{n\times
d}}}\times{\widehat{{\mathbb{R}}^{n\times
d}}}\rightarrow\mathbb{R}~{}~{},~{}~{}d(\hat{X},\hat{Y})=\min_{\Pi\in{\mathcal{S}}_{n}}{\|X-\Pi
Y\|}$

This makes $({\widehat{{\mathbb{R}}^{n\times d}}},d)$ a complete metric space.

Our main problem can now be stated as follows:

######  Problem 1.1.

Given $n,d\geq 1$ positive integers, find $m$ and a bi-Lipschitz map
$\hat{\alpha}:({\widehat{{\mathbb{R}}^{n\times
d}}},d)\rightarrow(\mathbb{R}^{m},{\|\cdot\|}_{2})$.

Explicitly the problem can be restated as follows. One is asked to construct a
map $\alpha:{\mathbb{R}}^{n\times d}\rightarrow\mathbb{R}^{m}$ that satisfies
the following conditions:

  1. (1)

If $X,Y\in{\mathbb{R}}^{n\times d}$ so that $X\sim Y$ then
$\alpha(X)=\alpha(Y)$

  2. (2)

If $X,Y\in{\mathbb{R}}^{n\times d}$ so that $\alpha(X)=\alpha(Y)$ then $X\sim
Y$

  3. (3)

There are constants $0<a_{0}\leq b_{0}$ so that for any
$X,Y\in{\mathbb{R}}^{n\times d}$,

(1.2) $a_{0}\min_{\Pi\in{\mathcal{S}}_{n}}{\|X-\Pi
Y\|}\leq{\|\alpha(X)-\alpha(Y)\|}_{2}\leq
b_{0}\min_{\Pi\in{\mathcal{S}}_{n}}{\|X-\Pi Y\|}$

Condition (1) allows us to lift $\alpha$ to the quotient space
${\widehat{{\mathbb{R}}^{n\times d}}}$. Thus $\hat{\alpha}(\hat{X})=\alpha(X)$
is well-defined. Condition (2) says that $\hat{\alpha}$ is injective (or, that
$\alpha$ is faithful with respect to the equivalence relation $\sim$).
Condition (3) says that $\hat{\alpha}$ is bi-Lipschitz with constants
$a_{0},b_{0}$. By a slight abuse of notation, when $\alpha$ satisfies (1) we
shall use the same letter to denote the map $\alpha:{\mathbb{R}}^{n\times
d}\rightarrow\mathbb{R}^{m}$ as well as the induced map on the quotient space
$\alpha:{\widehat{{\mathbb{R}}^{n\times d}}}\rightarrow\mathbb{R}^{m}$.

For $X,Y\in{\mathbb{R}}^{n\times d}$, $d(X,Y)$ denotes the same quantity in
(1.1) . In this case $d$ is only a semi-distance on ${\mathbb{R}}^{n\times
d}$, i.e., it is symmetric, non-negative and satisfies the triangle inequality
but fails the positivity condition.

One approach to embedding ${\widehat{{\mathbb{R}}^{n\times d}}}$ is to
consider the convex set of probability measures on $\mathbb{R}^{d}$,
${\mathcal{P}}(\mathbb{R}^{d})$, and the map

(1.3) $\alpha_{\infty}:{\mathbb{R}}^{n\times
d}\rightarrow{\mathcal{P}}(\mathbb{R}^{d})~{}~{},~{}~{}\alpha_{\infty}(X)=\frac{1}{n}\sum_{k=1}^{n}\delta(\cdot-
x_{k})$

where $[x_{1},\ldots,x_{n}]=X^{T}$, i.e., $x_{k}$ is the $k^{th}$ row of $X$
reshaped as a vector, and $\delta$ denotes the Dirac measure. When
${\mathcal{P}}(\mathbb{R}^{d})$ is endowed with the Wasserstein-1 distance
(the Earth Moving Distance), known also as the Kantorovich-Rubinstein metric,

$d_{KR}(p,q)=\inf_{\begin{array}[]{c}\mbox{$\pi\in{\mathcal{P}}(\mathbb{R}^{d}\times\mathbb{R}^{d}):$}\\\
\mbox{$\pi(\cdot,\mathbb{R}^{d})=p$}\\\
\mbox{$\pi(\mathbb{R}^{d},\cdot)=q$}\end{array}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}{\|x-y\|}d\pi(x,y)$

the distance between $a_{\infty}(X)$ and $a_{\infty}(Y)$ becomes

$d_{KR}(a_{\infty}(X),a_{\infty}(Y))=\min_{\Pi\in{\mathcal{S}}_{n}}\sum_{k=1}^{n}{\|x_{k}-(\Pi
Y)_{k}\|}.$

By the Kantorovich-Rubinstein theorem ([10]Theorem 1.14), $d_{KR}$ extends to
a norm on the linear space of bounded signed Borel measures on
$\mathbb{R}^{d}$, ${\mathcal{M}_{b}}(\mathbb{R}^{d})$. It is easy to verify
that

$d(\hat{X},\hat{Y})\leq
d_{KR}(a_{\infty}(X),a_{\infty}(Y))\leq\sqrt{n}d(\hat{X},\hat{Y})$

which proves that $a_{\infty}$ provides an embedding into a normed linear
space. Yet this embedding does not solve the problem since the linear space
${\mathcal{M}_{b}}(\mathbb{R}^{d})$ is infinite dimensional.

Instead of the previous infinite dimensional embedding, we consider two
different classes of embeddings. To illustrate these two constructions,
consider the simplest case $d=1$.

  1. (1)

Algebraic Embedding. For $x\in\mathbb{R}^{n}$, $x=(x_{1},\ldots,x_{n})^{T}$,
construct the polynomial $P_{x}(z)=(z-x_{1})\cdots(z-x_{n})$ and then expand
the product: $P_{x}(z)=z^{n}+c_{1}(x)z^{n-1}+\cdots+c_{n}(x)$. Using Vieta’s
formulas and Newton-Girard identities, an algebraically equivalent description
of $P_{x}$ is given by the symmetric polynomials:

(1.4)
$\alpha:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}~{}~{},~{}~{}\alpha(x)=\left(\sum_{k=1}^{n}x_{k},\sum_{k=1}^{n}x_{k}^{2},\ldots,\sum_{k=1}^{n}x_{k}^{n}\right).$

It is not hard to see that this map satisfies Conditions (1) and (2) and
therefore lifts to an injective continuous map $\hat{\alpha}$ on
$\hat{\mathbb{R}^{n}}$. Yet it is not Lipschitz, let alone bi-Lipschitz. The
approach in [20] can be used to modify $\alpha$ to a Lipschitz continuous map,
but, for the same reason as described in that paper, it cannot be “fixed” to a
bi-Lipschitz embedding. In Section 2 we show how to construct an algebraic
Lipschitz embedding in the case $d>1$.

  2. (2)

Sorting Embedding. For $x\in\mathbb{R}^{n}$, consider the sorting map

(1.5)
$\downarrow:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}~{}~{},~{}~{}\downarrow(x)=(x_{\pi(1)},x_{\pi(2)},\ldots,x_{\pi(n)})^{T}$

where the permutation $\pi$ is so that $x_{\pi(1)}\geq
x_{\pi(2)}\geq\cdots\geq x_{\pi(n)}$. It is obvious that $\downarrow$
satisfies Conditions (1) and (2) and therefore lifts to an injective map on
${\widehat{{\mathbb{R}}^{n\times d}}}$. As we see in Section 3, the map
$\downarrow$ is bi-Lipschitz. In fact it is isometric, and hence produces an
ideal embedding. Our work in Section 3 is to extend such construction to the
more general case $d>1$.

The algebraic embedding is a special case of the more general kernel method
that can be thought of as a projection of the measure $a_{\infty}(X)$ onto a
finite dimensional space, e.g., the space of polynomials spanned by
$\\{X,X^{2},\cdots,X^{n}\\}$. In applications such kernel method is known as a
“Readout Map” [40], based on “Sum Pooling”.

The sorting embedding has been used in applications under the name of “Pooling
Map” [40], based on “Max Pooling”. A naïve extension of the unidimensional map
(1.5) to the case $d>1$ might employ the lexicographic order: order monotone
decreasing the rows according to the first column, and break the tie by going
to the next column. While this gives rise to an injective map, it is easy to
see it is not even continuous, let alone Lipschitz. The main work in this
paper is to extend the sorting embedding to the case $d>1$ using a three-step
procedure, first embed ${\mathbb{R}}^{n\times d}$ into a larger vector space
$\mathbb{R}^{n\times D}$, then apply $\downarrow$ in each column
independently, and then perform a dimension reduction by a linear map into
$\mathbb{R}^{2nd}$. Similar to the phase retrieval problem ([2, 9, 4]), the
redundancy introduced in the first step counterbalances the loss of
information (here, relative order of one column with respect to another) in
the second step.

A summary of main results presented in this paper is contained in the
following result.

######  Theorem 1.2.

Consider the metric space $({\widehat{{\mathbb{R}}^{n\times d}}},d)$.

  1. (1)

(Polynomial Embedding) There exists a Lipschitz injective map

${\hat{\alpha}}:{\widehat{{\mathbb{R}}^{n\times d}}}\rightarrow\mathbb{R}^{m}$

with $m=\left(\begin{array}[]{c}\mbox{$d+n$}\\\ \mbox{$d$}\end{array}\right)$.
Two explicit constructions of this map are given in (2.8) and (2.9).

  2. (2)

(Sorting based Embedding) There exists a class of bi-Lipschitz maps

${\hat{\beta}}_{A,B}:({\widehat{{\mathbb{R}}^{n\times
d}}},d)\rightarrow(\mathbb{R}^{m},{\|\cdot\|})~{},~{}{\hat{\beta}}_{A,B}(\hat{X})=B\left({\hat{\beta}}_{A}(\hat{X})\right)$

with $m=2nd$, where each map ${\hat{\beta}}_{A,B}$ is the composition of two
bi-Lipschitz maps: a full-rank linear operator $B:\mathbb{R}^{n\times
D}\rightarrow\mathbb{R}^{m}$, with the nonlinear bi-Lipschitz map
${\hat{\beta}}_{A}:{\widehat{{\mathbb{R}}^{n\times
d}}}\rightarrow\mathbb{R}^{n\times D}$ parametrized by a matrix
$A\in\mathbb{R}^{d\times D}$ called ”key”. Explicitly,
${\hat{\beta}}(\hat{X})=\downarrow(XA)$, where $\downarrow$ acts column-wise.
These maps are characterized by the following properties:

    1. (a)

For $D=1+(d-1)n!$, any $A\in\mathbb{R}^{d\times(1+(d-1)n!)}$ whose columns
form a full spark frame defines a bi-Lipschitz map ${\hat{\beta}}_{A}$ on
${\widehat{{\mathbb{R}}^{n\times d}}}$. Furthermore, a lower Lipschitz
constant is given by the smallest $d^{th}$ singular value among all $d\times
d$ sub-matrices of $A$, $\min_{J\subset[D],|J|=d}s_{d}(A[J])$.

    2. (b)

For any matrix (“key”) $A\in\mathbb{R}^{d\times D}$ such that the map
${\hat{\beta}}_{A}$ is injective, then
${\hat{\beta}}_{A}:({\widehat{{\mathbb{R}}^{n\times
d}}},d)\rightarrow(\mathbb{R}^{n\times D},{\|\cdot\|})$ is bi-Lipschitz.
Furthermore, an upper Lipschitz constant is given by $s_{1}(A)$, the largest
singular value of $A$.

    3. (c)

Assume $A\in\mathbb{R}^{d\times D}$ is such that the map ${\hat{\beta}}_{A}$
is injective (i.e., a ”universal key”). Then for almost any linear map
$B:\mathbb{R}^{n\times D}\rightarrow\mathbb{R}^{2nd}$ the map
${\hat{\beta}}_{A,B}=B\circ{\hat{\beta}}_{A}$ is bi-Lipschitz.

An immediate consequence of this result is the following corollary whose proof
is included in subsection 3.5:

######  Corollary 1.3.

Let $\beta:\mathbb{R}^{n\times d}\rightarrow\mathbb{R}^{m}$ induce a bi-
Lipschitz embedding ${\hat{\beta}}:{\widehat{{\mathbb{R}}^{n\times
d}}}\rightarrow\mathbb{R}^{m}$ of the metric space
$({\widehat{{\mathbb{R}}^{n\times d}}},d)$ into
$(\mathbb{R}^{m},{\|\cdot\|}_{2})$.

  1. (1)

For any continuous function $f:\mathbb{R}^{n\times d}\rightarrow\mathbb{R}$
invariant to row-permutation (i.e., $f(PX)=f(X)$ for every
$X\in\mathbb{R}^{n\times d}$ and $P\in{\mathcal{S}}_{n}$) there exists a
continuous function $g:\mathbb{R}^{m}\rightarrow\mathbb{R}$ such that
$f=g\circ\beta$. Conversely, for any $g:\mathbb{R}^{m}\rightarrow\mathbb{R}$
continuous function, the function $f=g\circ\beta:\mathbb{R}^{n\times
d}\rightarrow\mathbb{R}$ is continuous and row-permutation invariant.

  2. (2)

For any Lipschitz continuous function $f:\mathbb{R}^{n\times
d}\rightarrow\mathbb{R}$ invariant to row-permutation (i.e., $f(PX)=f(X)$ for
every $X\in\mathbb{R}^{n\times d}$ and $P\in{\mathcal{S}}_{n}$) there exists a
Lipschitz continuous function $g:\mathbb{R}^{m}\rightarrow\mathbb{R}$ such
that $f=g\circ\beta$. Conversely, for any
$g:\mathbb{R}^{m}\rightarrow\mathbb{R}$ Lipschitz continuous function, the
function $f=g\circ\beta:\mathbb{R}^{n\times d}\rightarrow\mathbb{R}$ is
Lipschitz continuous and row-permutation invariant.

The structure of the paper is as follows. Section 2 contains the algebraic
embedding method and encoders $\alpha$ described at part (1) of Theorem 1.2.
Corollary 2.3 contains part (1) of the main result stated above. Section 3
introduces the sorting based embedding procedure and describes the key-based
encoder $\beta$. Necessary and sufficient conditions for key universality are
presented in Proposition 3.8; the injectivity of the encoder described at part
(2.a) of Theorem 1.2 is proved in Theorem 3.9; the bi-Lipschitz property of
any universal key described at part (2.b) of Theorem 1.2 is shown in Theorem
3.10; the dimension reduction statement (2.c) of Theorem 1.2 is included in
Theorem 3.13. Proof of Corollary 1.3 is presented in subsection 3.5. Section 4
contains applications to graph deep learning. These application use Graph
Convolution Networks and the numerical experiments are carried out on two
graph data sets: a chemical compound data set (QM9) and a protein data set
(PROTEINS_FULL).

While the motivation of this analysis is provided by graph deep learning
applications, this is primarily a mathematical paper. Accordingly the formal
theory is presented first, and then is followed by the machine learning
application. Those interested in the application (or motivation) can skip
directly to Section 4.

Notations. For an integer $d\geq 1$, $[d]=\\{1,2,\ldots,d\\}$. For a matrix
$X\in{\mathbb{R}}^{n\times d}$, $x_{1},\ldots x_{d}\in\mathbb{R}^{n}$ denote
its columns, $X=[x_{1}|\cdots|x_{d}]$. All norms are Euclidean; for a matrix
$X$, ${\|X\|}=\sqrt{trace(X^{T}X)}=\sqrt{\sum_{k,j}|X_{k,j}|^{2}}$ denotes the
Frobenius norm; for vectors $x$,
${\|x\|}={\|x\|}_{2}=\sqrt{\sum_{j}|x_{j}|^{2}}$.

###  1.1. Prior Works

Several methods for representing orbits of vector spaces under the action of
permutation (sub)groups have been studied in literature. Here we describe some
of these results, without claiming an exhaustive literature survey.

A rich body of literature emanated from the early works on symmetric
polynomials and group invariant representations of Hilbert, Noether, Klein and
Frobenius. They are part of standard commutative algebra and finite group
representation theory.

Prior works on permutation invariant mappings have predominantly employed some
form of summing procedure, though some have alternatively employed some form
of sorting procedure.

The idea of summing over the output nodes of an equivariant network has been
well studied. The algebraic invariant theory goes back to Hilbert and Noether
(for finite groups) and then continuing with the continuous invariant function
theory of Weyl and Wigner (for compact groups), who posited that a generator
function $\psi:X\rightarrow\mathbb{R}$ gives rise to a function
$E:X\rightarrow\mathbb{R}$ invariant to the action of a finite group $G$ on
$X$, $(g,x)\mapsto g.x$, via the averaging formula
$E(x)=\frac{1}{|G|}\sum_{g\in G}\psi(g.x)$.

More recently, this approach provided the framework for universal
approximation results of $G$-invariant functions. [27] showed that invariant
or equivariant networks must satisfy a fixed point condition. The equivariant
condition is naturally realized by GNNs. The invariance condition is realized
by GNNs when followed by summation on the output layer, as was further shown
in [21], [28] and [30]. Subsequently, [39] proved universal approximation
results over compact sets for continuous functions invariant to the action of
finite or continuous groups. In [16], the authors obtained bounds on the
separation power of GNNs in terms of the Weisfeiler-Leman (WL) tests by
tensorizing the input-output mapping. [35] studied approximations of
equivariant maps, while [11] showed that if a GNN with sufficient expressivity
is well trained, it can solve the graph isomorphism problem.

The authors of [36] designed an algorithm for processing sets with no natural
orderings. The algorithm applies an attention mechanism to achieve permutation
invariance with the attention keys being generated by a Long-Short Term Memory
(LSTM) network. Attention mechanisms amount to a weighted summing and
therefore can be considered to fall within the domain of summing based
procedures.

In [24], the authors designed a permutation invariant mapping for graph
embeddings. The mapping employs two separate neural networks, both applied
over the feature set for each node. One neural network produces a set of new
embeddings, the other serves as an attention mechanism to produce a weighed
sum of those new embeddings.

Sorting based procedures for producing permutation invariant mappings over
single dimensional inputs have been addressed and used by [40], notably in
their max pooling procedure.

The authors of [31] developed a permutation invariant mapping $pointnet$ for
point sets that is based on a $max$ function. The mapping takes in a set of
vectors, processes each vector through a neural network followed by an scalar
output function, and takes the maximum of the resultant set of scalars.

The paper [41] introduced SortPooling. SortPooling orders the latent
embeddings of a graph according to the values in a specific, predetermined
column. All rows of the latent embeddings are sorted according to the values
in that column. While this gives rise to an injective map, it is easy to see
it is not even continuous, let alone Lipschitz. The same issue arises with any
lexicographic ordering, including the well-known Weisfeiler-Leman embedding
[37]. Our paper introduces a novel method that bypasses this issue.

As shown in [28], the sum pooling-based GNNs provides universal approximations
for of any permutation invariant continuous function but only on _compacts_.
Our sorting based embedding removes the compactness restriction as well as it
extends to all Lipschitz maps.

While this paper is primarily mathematical in nature, methods developed here
are applied to two graph data sets, QM9 and PROTEINS_FULL. Researchers have
applied various graph deep learning techniques to both data sets. In
particular, [17] studied extensively the QM9 data set, and compared their
method with many other algorithms proposed by that time.

##  2\. Algebraic Embeddings

The algebraic embedding presented in this section can be thought of a special
kernel to project equation (1.3) onto.

###  2.1. Kernel Methods

The kernel method employs a family of continuous kernels (test) functions,
$\\{K(x;y)~{};~{}x\in\mathbb{R}^{d}~{},~{}y\in Y\\}$ parametrized/indexed by a
set $Y$. The measure representation $\mu=a_{\infty}(X)$ in (1.3) yields a
nonlinear map

$\alpha:\mathbb{R}^{n\times d}\rightarrow C(Y)~{}~{},~{}~{}X\mapsto
F(y)=\int_{R^{d}}K(x;y)d\mu$

given by

$\alpha(X)(y)=\frac{1}{n}\sum_{k=1}^{n}K(x_{k};y)$

The embedding problem 1.1) can be restated as follows. One is asked to find a
finite family of kernels $\\{K(x;y)~{};~{}x\in\mathbb{R}^{d}~{},~{}y\in Y\\}$,
$m=|Y|$ so that

(2.1) ${\hat{\alpha}}:({\widehat{{\mathbb{R}}^{n\times d}}},d)\rightarrow
l^{2}(Y)\sim(\mathbb{R}^{m},{\|\cdot\|}_{2})~{}~{},~{}~{}({\hat{\alpha}}(\hat{X}))_{y}=\frac{1}{n}\sum_{k=1}^{n}K(x_{k};y)$

is injective, Lipschitz or bi-Lipschitz.

Two natural choices for the kernel $K$ are the Gaussian kernel and the complex
exponential (or, the Fourier) kernel:

$K_{G}(x,y)=e^{-{\|x-y\|}^{2}/\sigma^{2}}~{}~{},K_{F}(x,y)=e^{2\pi i{\langle
x,y\rangle}}$

where in both cases $Y\subset\mathbb{R}^{d}$. In this paper we analyze a
different kernel, namely the polynomial kernel
$K_{P}(x,y)=x_{1}^{y_{1}}x_{2}^{y_{2}}\cdots x_{d}^{y_{d}}$,
$Y\subset\\{0,1,2,\ldots,n\\}^{d}$.

###  2.2. The Polynomial Embedding

Since the polynomial representation is intimately related to the Hilbert-
Noether algebraic invariants theory [18] and the Hilbert-Weyl theorem, it is
advantageous to start our construction from a different perspective.

The linear space ${\mathbb{R}}^{n\times d}$ is isomorphic to $\mathbb{R}^{nd}$
by stacking the columns one on top of each other. In this case, the action of
the permutation group $S_{n}$ can be recast as the action of the subgroup
$I_{d}\otimes S_{n}$ of the bigger group $S_{nd}$ on $\mathbb{R}^{nd}$.
Specifically, let us denote by $\sim_{G}$ the equivalence relation

$x,y\in\mathbb{R}^{nd}~{}~{},~{}~{}x\sim_{G}y\Longleftrightarrow y=\Pi
x~{},~{}{\rm for~{}some}~{}\Pi\in G$

induced by a subgroup $G$ of $S_{nd}$. In the case $G=I_{d}\otimes
S_{n}=\\{diag_{d}(P)~{},~{}P\in S_{n}\\}$ of block diagonal permutation
obtained by repeating $d$ times the same $P\in S_{n}$ permutation along the
main diagonal, two vectors $x,y\in\mathbb{R}^{nd}$ are $\sim_{G}$ equivalent
iff there is a permutation matrix $P\in S_{n}$ so that
$y(1+(k-1)n:kn)=Px(1+(k-1)n:kn)$ for each $1\leq k\leq d$. In other words,
each disjoint $n$-subvectors in $y$ and $x$ are related by the same
permutation. In this framework, the Hilbert-Weyl theorem (Theorem 4.2, Chapter
XII, in [25]) states that the ring of invariant polynomials is finitely
generated. The Göbel’s algorithm (Section 3.10.2 in [18]) provides a recipe to
find a complete set of invariant polynomials. In the following we provide a
direct approach to construct a complete set of polynomial invariants.

Let $\mathbb{R}[{\bf x}_{1},{\bf x}_{2},...,{\bf x}_{d}]$ denote the algebra
of polynomials in $d$-variables with real coefficients. Let us denote
$X\in{\mathbb{R}}^{n\times d}$ a generic data matrix. Each row of this matrix
defines a linear form over ${\bf x}_{1},...{\bf x}_{d}$,
$\lambda_{k}=X_{k,1}{\bf x}_{1}+\cdots+X_{k,d}{\bf x}_{d}$. Let us denote by
$\mathbb{R}[{\bf x}_{1},\ldots,{\bf x}_{d}][{\bf t}]$ the algebra of
polynomials in variable ${\bf t}$ with coefficients in the ring
$\mathbb{R}[{\bf x}_{1},\ldots,{\bf x}_{d}]$. Notice $\mathbb{R}[{\bf
x}_{1},{\bf x}_{2},\ldots,{\bf x}_{d}][{\bf t}]=\mathbb{R}[{\bf t},{\bf
x}_{1},{\bf x}_{2},\ldots,{\bf x}_{d}]$ by rearranging the terms according to
degree in ${\bf t}$. Thus $\lambda_{k}\in\mathbb{R}[{\bf x}_{1},\ldots,{\bf
x}_{d}]\subset\mathbb{R}[{\bf x}_{1},\ldots,{\bf x}_{d}][{\bf t}]$ can be
encoded as zeros of a polynomial $P_{X}$ of degree $n$ in variable ${\bf t}$
with coefficients in $\mathbb{R}[{\bf x}_{1},\ldots,{\bf x}_{d}]$:

(2.2) $P_{X}({\bf t},{\bf x}_{1},\ldots,{\bf x}_{d})=\prod_{k=1}^{n}({\bf
t}-\lambda_{k}({\bf x}_{1},\ldots,{\bf x}_{d}))=\prod_{k=1}^{n}({\bf
t}-X_{k,1}{\bf x}_{1}-\ldots-X_{k,d}{\bf x}_{d})$

Due to identification $\mathbb{R}[{\bf x}_{1},{\bf x}_{2},\ldots,{\bf
x}_{d}][{\bf t}]=\mathbb{R}[{\bf t},{\bf x}_{1},{\bf x}_{2},\ldots,{\bf
x}_{d}]$, we obtain that
$P_{X}\in\mathbb{R}[{\bf t},{\bf x}_{1},{\bf x}_{2},\ldots,{\bf x}_{d}]$ is a
homogeneous polynomial of degree $n$ in $d+1$ variables. Let
$\mathbb{R}_{n}[{\bf t},{\bf x}_{1},\ldots,{\bf x}_{d}]$ denote the vector
space of homogeneous polynomials in $d+1$ variables of degree $n$ with real
coefficients. Notice the real dimension of this vector space is

(2.3) $\dim_{\mathbb{R}}\mathbb{R}_{n}[{\bf t},{\bf x}_{1},\ldots,{\bf
x}_{d}]=\left(\begin{array}[]{c}n+d\\\
d\end{array}\right)=\left(\begin{array}[]{c}n+d\\\ n\end{array}\right).$

By noting that $P_{X}$ is monic in ${\bf t}$ (the coefficient of ${\bf t}^{n}$
is always 1) we obtain an injective embedding of
${\widehat{{\mathbb{R}}^{n\times d}}}$ into $\mathbb{R}^{m}$ with
$m=\dim_{\mathbb{R}}\mathbb{R}_{n}[{\bf t},{\bf x}_{1},\ldots,{\bf x}_{d}]-1$
via the coefficients of $P_{X}$ similar to (1.4). This is summarized in the
following theorem:

######  Theorem 2.1.

The map $\alpha_{0}:{\mathbb{R}}^{n\times d}\rightarrow\mathbb{R}^{m-1}$ with
$m=\left(\begin{array}[]{c}n+d\\\ d\end{array}\right)$ given by the (non-
trivial) coefficients of polynomial $P_{X}\in\mathbb{R}_{n}[{\bf t},{\bf
x}_{1},\ldots,{\bf x}_{d}]$ lifts to an analytic embedding
${\hat{\alpha}}_{0}$ of $({\widehat{{\mathbb{R}}^{n\times d}}},d)$ into
$\mathbb{R}^{m}$. Specifically, for $X\in{\mathbb{R}}^{n\times d}$ expand the
polynomial

(2.4) $P_{X}({\bf t},{\bf x}_{1},\ldots,{\bf x}_{d})=\prod_{k=1}^{n}({\bf
t}-X_{k,1}{\bf x}_{1}-\ldots-X_{k,d}{\bf x}_{d})={\bf
t}^{n}+\sum_{\begin{array}[]{c}\mbox{$p_{0},p_{1},...,p_{d}\geq 0$}\\\
\mbox{$p_{0}+\cdots+p_{d}=n$}\\\
\mbox{$p_{0}<n$}\end{array}}c_{p_{0},p_{1},\ldots,p_{d}}{\bf t}^{p_{0}}{\bf
x}_{1}^{p_{1}}\cdots{\bf x}_{d}^{p_{d}}$

Then

(2.5) $X\in{\mathbb{R}}^{n\times
d}\mapsto\alpha_{0}(X)=(c_{p_{0},p_{1},\ldots,p_{d}})_{(p_{0},p_{1},\ldots,p_{d})\in
I_{n,d}}$

where the index set is given by

(2.6) $I_{n,d}=\\{(p_{0},p_{1},\ldots,p_{d})~{},~{}0\leq
p_{0},p_{1},\ldots,p_{d}~{},~{}p_{0}<n~{},~{}p_{0}+p_{1}+\cdots+p_{d}=n\\}$

and is of cardinal $|I_{n,d}|=m-1$. The map
${\hat{\alpha}}_{0}:{\widehat{{\mathbb{R}}^{n\times
d}}}\rightarrow\mathbb{R}^{m-1}$ is the lifting of $\alpha_{0}$ to the
quotient space.

Proof

Since for any permutation $\pi$ with associated permutation matrix
$\Pi\in{\mathcal{S}}_{n}$,

$P_{\Pi X}({\bf t},{\bf x}_{1},\cdots,{\bf x}_{d})=\prod_{k=1}^{n}({\bf
t}-X_{\pi(k),1}{\bf x}_{1}-\ldots-X_{\pi(k),d}{\bf x}_{d})=P_{X}({\bf t},{\bf
x}_{1},\ldots,{\bf x}_{d}),$

it follows that $\alpha_{0}$ is invariant to the action of
${\mathcal{S}}_{n}$, $\alpha_{0}(X)=\alpha_{0}(\Pi X)$. Thus $\alpha_{0}$
lifts to a map ${\hat{\alpha}}_{0}$ on ${\widehat{{\mathbb{R}}^{n\times d}}}$.

The coefficients of polynomial $P_{X}$ depend analytically on its roots
(Vieta’s formulas), hence on entries of matrix $X$.

The only remaining claim is that if $X,Y\in{\mathbb{R}}^{n\times d}$ so that
$\alpha_{0}(X)=\alpha_{0}(Y)$ then there is $\Pi\in{\mathcal{S}}_{n}$ so that
$Y=\Pi X$. Assume $P_{X}=P_{Y}$. For each choice $({\bf x}_{1},{\bf
x}_{2},\ldots,{\bf x}_{d})=(f(1),\ldots,f(d))$ in $\mathbb{R}^{d}$, the $n$
real zeros of the two polynomials in ${\bf t}$, $P_{X}({\bf
t},f(1),\ldots,f(d))$ and $P_{Y}({\bf t},f(1),\ldots,f(d))$ coincide.
Therefore $Xf\sim Yf$ for each $f\in\mathbb{R}^{d}$, Let $D=1+(d-1)n!$ and
choose $F\in\mathbb{R}^{d\times D}$ so that each subset of $d$ columns are
linearly independent, in other words, the set
${\mathcal{F}}=\\{f_{1},f_{2},\ldots,f_{D}\\}$ formed by the $D$ columns of
$F$ is a full spark frame in $\mathbb{R}^{d}$, see [1]. As proved in [1],
almost every such set is a full spark frame. Then for each $1\leq k\leq D$
there is a permutation $\Pi_{k}\in{\mathcal{S}}_{n}$ so that
$Xf_{k}=\Pi_{k}Yf_{k}$. By the pigeonhole principle, since
$|{\mathcal{S}}_{n}|=n!$, there are $1\leq k_{1}<k_{2}<\cdots<k_{d}\leq D$ so
that $\Pi_{k_{1}}=\Pi_{k_{2}}=\cdots=\Pi_{k_{d}}$. Then
$(X-\Pi_{k_{1}}Y)f_{k_{j}}=0$ for every $1\leq j\leq d$. Since
$\\{f_{k_{1}},\ldots,f_{k_{d}}\\}$ is linearly independent it follows that
$X-\Pi_{k_{1}}Y=0$. Thus $X\sim Y$ which ends the proof of this result. ∎.

######  Remark 2.2.

The set of invariants produced by map $\alpha_{0}$ are proportional to those
produced by the Göbel’s algorithm in [18], §3.10.2. Indeed, the $nd$ primary
invariants are given by

$\\{c_{p,n-p,0,\ldots,0}~{},~{}0\leq p\leq
n-1\\}\cup\cdots\cup\\{c_{p,0,\ldots,0,n-p}~{},~{}0\leq p\leq n-1\\}$

corresponding to the elementary symmetric polynomials in entries of each
column. The secondary invariants correspond to the remaining coefficients that
have at least 2 nonzero indices among $p_{1},\ldots,p_{d}$.

The embedding provided by $\alpha_{0}$ is analytic and injective but is not
globally Lipschitz because of the polynomial growth rate. Next we show how a
simple modification of this map will make it Lipschitz. First, let us denote
by $L_{0}$ the Lipschitz constant of $\alpha_{0}$ when restricted to the
closed unit ball $B_{1}({\mathbb{R}}^{n\times d}):\\{X\in{\mathbb{R}}^{n\times
d}~{},~{}{\|X\|}\leq 1\\}$ of ${\mathbb{R}}^{n\times d}$, i.e.
${\|\alpha_{0}(X)-\alpha_{0}(Y)\|}\leq L_{0}{\|X-Y\|}$ for any
$X,Y\in{\mathbb{R}}^{n\times d}$ with ${\|X\|},{\|Y\|}\leq 1$. Second, let

(2.7)
$\varphi_{0}:\mathbb{R}\rightarrow[0,1]~{},~{}\varphi_{0}(x)=min(1,\frac{1}{x})=\left\\{\begin{array}[]{rcl}1&if&\mbox{$x\leq
1$}\\\ \mbox{$\frac{1}{x}$}&if&\mbox{$x>1$}\end{array}\right.$

be a Lipschitz monotone decreasing function with Lipschitz constant 1.

######  Corollary 2.3.

Consider the map:

(2.8) $\alpha_{1}:{\mathbb{R}}^{n\times
d}\rightarrow\mathbb{R}^{m}~{}~{},~{}~{}\alpha_{1}(X)=\left(\begin{array}[]{c}\mbox{$\alpha_{0}\bigg{(}\varphi_{0}({\|X\|})X\bigg{)}$}\\\
\mbox{${\|X\|}$}\end{array}\right),$

with $m=\left(\begin{array}[]{c}n+d\\\ d\end{array}\right)$. The map
$\alpha_{1}$ lifts to an injective and globally Lipschitz map
${\hat{\alpha}}_{1}:{\widehat{{\mathbb{R}}^{n\times
d}}}\rightarrow\mathbb{R}^{m}$ with Lipschitz constant
$Lip({\hat{\alpha}}_{1})\leq\sqrt{1+L_{0}^{2}}$.

Proof

Clearly $\alpha_{1}(\Pi X)=\alpha_{1}(X)$ for any $\Pi\in{\mathcal{S}}_{n}$
and $X\in{\mathbb{R}}^{n\times d}$. Assume now that
$\alpha_{1}(X)=\alpha_{1}(Y)$. Then ${\|X\|}={\|Y\|}$ and since
${\hat{\alpha}}_{0}$ is injective on ${\widehat{{\mathbb{R}}^{n\times d}}}$ it
follows $\varphi({\|X\|})X=\Pi\varphi({\|Y\|})Y$ for some
$\Pi\in{\mathcal{S}}_{n}$. Thus $X\sim Y$ which proves $\alpha_{1}$ lifts to
an injective map on ${\widehat{{\mathbb{R}}^{n\times d}}}$.

Now we show ${\hat{\alpha}}_{1}$ is Lipschitz on
$({\widehat{{\mathbb{R}}^{n\times d}}},d)$ of appropriate Lipschitz constant.
Let $X,Y^{\prime}\in{\mathbb{R}}^{n\times d}$ and
$\Pi_{0}\in{\mathcal{S}}_{n}$ so that
$d(\hat{X},\hat{Y^{\prime}})={\|X-\Pi_{0}Y^{\prime}\|}$. Let
$Y=\Pi_{0}Y^{\prime}$ so that $d(\hat{X},\hat{Y})={\|X-Y\|}$.

Choose two matrices $X,Y\in{\mathbb{R}}^{n\times d}$. We claim
${\|\alpha_{1}(X)-\alpha_{1}(Y)\|}\leq\sqrt{1+L_{0}^{2}}{\|X-Y\|}$. This
follows from two observations:

(i) The map

$X\mapsto\rho(X):=\varphi_{0}({\|X\|})X$

is the nearest-point map to (or, the metric projection map onto) the convex
closed set $B_{1}({\mathbb{R}}^{n\times d})$. This means
${\|\varphi_{0}({\|X\|})X-Z\|}\leq{\|X-Z\|}$ for any $Z\in
B_{1}({\mathbb{R}}^{n\times d})$.

(ii) The nearest-point map to a convex closed subset of a Hilbert space is
Lipschitz with constant 1, i.e. it shrinks distances, see [29].

These two observations yield:

${\|\alpha_{1}(X)-\alpha_{1}(Y)\|}^{2}={\|\alpha_{0}(\rho(Y))-\alpha_{0}(\rho(Y))\|}^{2}+|{\|X\|}-{\|Y\|}|^{2}\\\
\leq
L_{0}^{2}{\|\rho(X)-\rho(Y)\|}^{2}+{\|X-Y\|}^{2}\leq(1+L_{0}^{2}){\|X-Y\|}^{2}.$

This concludes the proof of this result. ∎

A simple modification of $\phi_{0}$ can produce a $C^{\infty}$ map by
smoothing it out around $x=1$.

On the other hand the lower Lipschitz constant of ${\hat{\alpha}}_{1}$ is 0
due to terms of the form $X_{i,j}^{k}$ with $k\geq 2$. In [20], the authors
built a Lipschitz map by a retraction to the unit sphere instead of unit ball.
Inspired by their construction, a modification of $\alpha_{0}$ in their spirit
reads:

(2.9) $\alpha_{2}:{\mathbb{R}}^{n\times
d}\rightarrow\mathbb{R}^{m}~{}~{},~{}~{}\alpha_{2}(X)=\left(\begin{array}[]{c}\mbox{${\|X\|}\alpha_{0}\bigg{(}\frac{X}{{\|X\|}}\bigg{)}$}\\\
\mbox{${\|X\|}$}\end{array}\right)~{},~{}if~{}X\neq
0~{}~{},~{}and~{}~{}\alpha_{2}(0)=0.$

It is easy to see that $\alpha_{2}$ satisfies the non-parallel property in
[20] and is Lipschitz with a slightly better constant than $\alpha_{1}$ (the
constant is determined by the tangential derivatives of $\alpha_{0}$). But,
for the same reasons as in [20] this map is not bi-Lipschitz.

###  2.3. Dimension reduction in the case $d=2$ and consequences

In this subsection we analyze the case $d=2$. The embedding dimension for
$\alpha_{0}$ is $\left(\begin{array}[]{c}n\\\
2\end{array}\right)-1=\frac{n(n-1)}{2}-1$. On the other hand, consider the
following approach. Each row of $X$ defines a complex number
$z_{1}=X_{1,1}+i\,X_{1,2}$, … , $z_{n}=X_{n,1}+i\,X_{n,2}$ that can be encoded
by one polynomial of degree $n$ with complex coefficients
$Q\in\mathbb{C}_{n}[t]$,

$Q({\bf t})=\prod_{k=1}^{n}({\bf t}-z_{k})={\bf t}^{n}+\sum_{k=0}^{n-1}{\bf
t}^{k}q_{k}$

The coefficients of $Q$ provide a $2n$-dimensional real embedding $\zeta_{0}$,

$\zeta_{0}:\mathbb{R}^{n\times
2}\rightarrow\mathbb{R}^{2n}~{}~{},~{}~{}\zeta_{0}(X)=(Re(q_{n-1}),Im(q_{n-1}),\ldots,Re(q_{0}),Im(q_{0}))$

with properties similar to those of $\alpha_{0}$. One can similarly modify
this embedding to obtain a globally Lipschitz embedding $\hat{\zeta}_{1}$ of
$\hat{R_{n,2}}$ into $\mathbb{R}^{2n+1}$.

It is instructive to recast this embedding in the framework of commutative
algebras. Indeed, let $\langle{\bf x}_{1}-1,{\bf x}_{2}^{2}+1\rangle$ denote
the ideal generated by polynomials ${\bf x}_{1}-1$ and ${\bf x}_{2}^{2}+1$ in
the algebra $\mathbb{R}[{\bf t},{\bf x}_{1},{\bf x}_{2}]$. Consider the
quotient space $\mathbb{R}[{\bf t},{\bf x}_{1},{\bf x}_{2}]/\langle{\bf
x}_{1}-1,{\bf x}_{2}^{2}+1\rangle$ and the quotient map
$\sigma:\mathbb{R}[{\bf t},{\bf x}_{1},{\bf x}_{2}]\mapsto\mathbb{R}[{\bf
t},{\bf x}_{1},{\bf x}_{2}]/\langle{\bf x}_{1}-1,{\bf x}_{2}^{2}+1\rangle$. In
particular, let $S=\sigma(\mathbb{R}_{n}[{\bf t},{\bf x}_{1},{\bf x}_{2}])$
denote the vector space projected through this quotient map. Then a basis for
$S$ is given by $\\{1,{\bf t},\ldots,{\bf t}^{n},{\bf x}_{2},{\bf x}_{2}{\bf
t},\ldots,{\bf x}_{2}{\bf t}^{n-1},{\bf x}_{2}{\bf t}^{n}\\}$. Thus $\dim
S=2n+2$. Let $\mathfrak{S}=\\{P_{X}~{},~{}X\in\mathbb{R}^{n\times
2}\\}\subset\mathbb{R}_{2}[{\bf t},{\bf x}_{1},{\bf x}_{2}]$ denote the set of
polynomials realizable as in (2.4). Then the fact that
$\hat{\zeta}_{0}:\mathbb{R}^{n\times 2}\rightarrow\mathbb{R}^{2n}$ is
injective is equivalent to the fact that
$\sigma{|}_{\mathfrak{S}}:\mathfrak{S}\rightarrow S$ is injective. On the
other hand note

$\sigma(\mathfrak{S})\subset{\bf t}^{n}+span_{\mathbb{R}}\\{1,{\bf
t},\ldots,{\bf t}^{n-1},{\bf x}_{2},{\bf x}_{2}{\bf t},\ldots,{\bf x}_{2}{\bf
t}^{n-1}\\}$

where the last linear subspace is of dimension $2n$.

In the case $d=2$ we obtain the identification $\mathbb{R}[{\bf t},{\bf
x}_{1},{\bf x}_{2}]/\langle{\bf x}_{1}-1,{\bf
x}_{2}^{2}+1\rangle=\mathbb{C}[{\bf t}]$ due to uniqueness of polynomial
factorization.

This observation raises the following open problem:

For $d>2$, is there a non-trivial ideal $I=\langle
Q_{1},\ldots,Q_{r}\rangle\subset\mathbb{R}[{\bf t},{\bf x}_{1},\ldots,{\bf
x}_{d}]$ so that the restriction $\sigma{|}_{\mathfrak{S}}$ of the quotient
map $\sigma:\mathbb{R}[{\bf t},{\bf x}_{1},\ldots,{\bf
x}_{d}]\rightarrow\mathbb{R}[{\bf t},{\bf x}_{1},\ldots,{\bf x}_{d}]/I$ is
injective? Here $\mathfrak{S}$ denote the set of polynomials in
$\mathbb{R}_{n}[{\bf t},{\bf x}_{1},\ldots,{\bf x}_{d}]$ realizable via (2.4).

######  Remark 2.4.

One may ask the question whether the quaternions can be utilized in the case
$d=4$. While the quaternions form an associative division algebra,
unfortunately polynomials have in general an infinite number of factorization.
This prevents an immediate extension of the previous construction to the case
$d=4$.

######  Remark 2.5.

Similar to the construction in [20], a linear dimension reduction technique
may be applicable here (which, in fact, may answer the open problem above)
which would reduce the embedding dimension to $m=2nd+1$ (twice the intrinsec
dimension plus one for the homogenization variable). However we did not
explore this approach since, even if possible, it would not produce a bi-
Lipschitz embedding. Instead we analyze the linear dimension reduction
technique in the next section in the context of sorting based embeddings.

##  3\. Sorting based Embedding

In this section we present the extension of the sorting embedding (1.5) to the
case $d>1$.

The embedding is performed by a linear-nonlinear transformation that resembles
the phase retrieval problem. Consider a matrix $A\in\mathbb{R}^{d\times D}$
and the induced nonlinear transformation:

(3.1) $\beta_{A}:{\mathbb{R}}^{n\times d}\rightarrow\mathbb{R}^{n\times
D}~{}~{},~{}~{}\beta_{A}(X)=\downarrow(XA)$

where $\downarrow$ is the monotone decreasing sorting operator acting in each
column independently. Specifically, let $Y=XA\in\mathbb{R}^{n\times D}$ and
note its column vectors $Y=[y_{1},y_{2},\ldots,y_{D}]$. Then

$\beta_{A}(X)=\left[\begin{array}[]{cccc}\mbox{$\Pi_{1}y_{1}$}&\mbox{$\Pi_{2}y_{1}$}&\cdots&\mbox{$\Pi_{D}y_{D}$}\end{array}\right]$

for some $\Pi_{1},\Pi_{2},\ldots,\Pi_{D}\in{\mathcal{S}}_{n}$ so that each
column is sorted monotonically decreasing:

$(\Pi_{k}y_{k})_{1}\geq(\Pi_{k}y_{k})_{2}\geq\cdots\geq(\Pi_{k}y_{k})_{n}.$

Note the obvious invariance $\beta_{A}(\Pi X)=\beta_{A}(X)$ for any
$\Pi\in{\mathcal{S}}_{n}$ and $X\in{\mathbb{R}}^{n\times d}$. Hence
$\beta_{A}$ lifts to a map $\hat{\beta_{A}}$ on
${\widehat{{\mathbb{R}}^{n\times d}}}$.

######  Remark 3.1.

Notice the similarity to the phase retrieval problem, e.g., [4], where the
data is obtained via a linear transformation of the input signal followed by
the nonlinear operation of taking the absolute value of the frame
coefficients. Here the nonlinear transformation is implemented by sorting the
coefficients. In both cases it represents the action of a particular subgroup
of the unitary group.

In this section we analyze necessary and sufficient conditions so that maps of
type (3.1) are injective, or injective almost everywhere. First a few
definitions.

######  Definition 3.2.

A matrix $A\in\mathbb{R}^{d\times D}$ is called a _universal key_ (for
${\mathbb{R}}^{n\times d}$) if $\hat{\beta_{A}}$ is injective on
${\widehat{{\mathbb{R}}^{n\times d}}}$.

In general we refer to $A$ as a key for encoder $\beta_{A}$.

######  Definition 3.3.

Fix a matrix $X\in{\mathbb{R}}^{n\times d}$. A matrix $A\in\mathbb{R}^{d\times
D}$ is said _admissible_ (or an admissible key) for $X$ if for any
$Y\in{\mathbb{R}}^{n\times d}$ so that $\beta_{A}(X)=\beta_{A}(Y)$ then $Y=\Pi
X$ for some $\Pi\in{\mathcal{S}}_{n}$.

In other words,
$\hat{\beta_{A}}^{-1}(\hat{\beta_{A}}(\hat{X}))=\\{\hat{X}\\}$. We let
${\mathcal{A}}_{D}(X)$, or simply ${\mathcal{A}}(X)$, denote the set of
admissible keys for $X$.

######  Definition 3.4.

Fix $A\in\mathbb{R}^{d\times D}$. A matrix $X\in{\mathbb{R}}^{n\times d}$ is
said to be separated by $A$ if $A\in{\mathcal{A}}(X)$.

For a key $A$, we let $\mathfrak{S}_{n}(A)$, or simply $\mathfrak{S}(A)$,
denote the set of matrices separated by $A$. Thus a matrix
$X\in\mathfrak{S}_{n}(A)$ if and only if, for any matrix
$Y\in\mathbb{R}^{n\times d}$, if $\beta_{A}(X)=\beta_{A}(Y)$ then $X\sim Y$.

Thus a key $A$ is universal if and only if
$\mathfrak{S}_{n}(A)={\mathbb{R}}^{n\times d}$.

Our goal is to produce keys that are admissible for all matrices in
${\mathbb{R}}^{n\times d}$, or at least for almost every data matrix. As we
show in Proposition 3.6 below this requires that $D\geq d$ and $A$ is full
rank. In particular this means that the columns of $A$ form a frame for
$\mathbb{R}^{d}$.

###  3.1. Characterizations of ${\mathcal{A}}(X)$ and $\mathfrak{S}(A)$

We start off with simple linear manipulations of sets of admissible keys and
separated data matrices.

######  Proposition 3.5.

Fix $A\in\mathbb{R}^{d\times D}$ and $X\in{\mathbb{R}}^{n\times d}$.

  1. (1)

For an invertible $d\times d$ matrix $T\in\mathbb{R}^{d\times d}$,

(3.2) $\mathfrak{S}_{n}(TA)=\mathfrak{S}_{n}(A)T^{-1}.$

In other words, if $X$ is separated by $A$ then $XT^{-1}$ is separated by
$TA$.

  2. (2)

For any permutation matrix $L\in\SS_{D}$ and diagonal invertible matrix
$\Lambda\in\mathbb{R}^{D\times D}$,

(3.3) $\mathfrak{S}_{n}(AL\Lambda)=\mathfrak{S}_{n}(A\Lambda
L)=\mathfrak{S}_{n}(A).$

In other words, if $X$ is separated by $A$ then $X$ is separated also by
$AL\Lambda$ as well as by $A\Lambda L$.

  3. (3)

Assume $T\in\mathbb{R}^{d\times d}$ is a $d\times d$ invertible matrix. Then

(3.4) ${\mathcal{A}}_{D}(XT)=T^{-1}{\mathcal{A}}_{D}(X).$

In other words, if $A$ is an admissible key for $X$ then $T^{-1}A$ is an
admissible key for $XT$.

Proof

The proof is immediate, but we include it here for convenience of the reader.

(1) Denote $B=TA$. Let $Y\in\mathbb{R}^{n\times d}$. Then

$\beta_{B}(Y)=\beta_{B}(X)\Longleftrightarrow\downarrow(XB)=\downarrow(YB)\Longleftrightarrow\downarrow(XTA)=\downarrow(YTA)\Longleftrightarrow\beta_{A}(XT)=\beta_{A}(YT).$

Thus, if $X\in\mathfrak{S}_{n}(A)$ and $Y^{\prime}\in{\mathbb{R}}^{n\times d}$
so that $\beta_{B}(Y^{\prime})=\beta_{B}(X^{\prime})$ with
$X^{\prime}=XT^{-1}$, then $\beta_{A}(Y^{\prime}T)=\beta_{A}(X)$. Therefore
there exists $\Pi\in{\mathcal{S}}_{n}$ so that $Y^{\prime}T=\Pi X$. Thus
$Y^{\prime}\sim X^{\prime}$. Hence $X^{\prime}\in\mathfrak{S}_{n}(B)$. This
shows $\mathfrak{S}_{n}(A)T^{-1}\subset\mathfrak{S}_{n}(TA)$. The reverse
include follows by replacing $A$ with $TA$ and $T$ with $T^{-1}$. Together
they prove (3.2).

(2) Let $Y\in{\mathbb{R}}^{n\times d}$ such that
$\beta_{AL\Lambda}(X)=\beta_{AL\Lambda}(Y)$. For every $1\leq j\leq D$ let
$k\in[D]$ be so that $L_{jk}=1$.

If $\Lambda_{kk}>0$ then $\downarrow((XA)_{j})=\downarrow((YA)_{j})$.

If $\Lambda_{kk}<0$ then $\downarrow(-(XA)_{j})=\downarrow(-(YA)_{j})$. But
this implies also $\downarrow((XA)_{j})=\downarrow((YA)_{j})$ since
$\downarrow(-z)=L_{0}\downarrow(z)$ where $L_{0}$ is the permutation matrix
that has 1 on its main antidiagonal.

Either way, $\downarrow((XA)_{j})=\downarrow((YA)_{j})$. Hence
$\downarrow(XA)=\downarrow(YA)$. Therefore $X\sim Y$ and thus
$X\in\mathfrak{S}_{n}(AL\Lambda)$. This shows
$\mathfrak{S}_{n}(A)\subset\mathfrak{S}_{n}(AL\Lambda)$. the reverse inclusion
follows by a similar argument. Finally, notice $\\{L\Lambda\\}$ forms a group
since $L^{-1}\Lambda L$ is also a diagonal matrix. This shows
$\mathfrak{S}_{n}(A\Lambda L)=\mathfrak{S}(AL\Lambda^{\prime})$ for some
diagonal matrix $\Lambda^{\prime}$, and the conclusion (3.3) then follows.

(3) The relation (3.4) follows from noticing
$\beta_{T^{-1}A}(Y)=\beta_{A}(YT)$. ∎

Relation (3.3) shows that, since $A$ is assumed full rank, without loss of
generality we can assume the first $d$ columns are linearly independent. Let
$V$ denote the first $d$ columns of $A$ so that

(3.5)
$A=V\left[\begin{array}[]{ccc}\mbox{$I$}&\mbox{$|$}&\mbox{$\tilde{A}$}\end{array}\right]$

where $\tilde{A}\in\mathbb{R}^{d\times(D-d)}$. The following result shows
that, unsurprisingly, when $D=d>1$, almost every matrix $X$ is not separated
by $A$. By Proposition 3.5 we can reduce the analysis to the case $A=I$ by a
change of coordinates.

######  Proposition 3.6.

Assume $D=d>1$, $n>1$. Then

  1. (1)

The set of data matrices not separated by $I_{d}$ includes:

(3.6) $\mathbb{B}:=\\{X\in\mathbb{R}^{n\times d}~{},~{}\exists
i,j,k,l~{},~{}1\leq i<j\leq n,1\leq k<l\leq d~{}\Rightarrow~{}X_{i,k}\neq
X_{j,k}~{}\&~{}X_{i,l}\neq X_{j,l}\\}\subset\mathfrak{S}_{n}(I_{d})^{c}.$

  2. (2)

The set $\mathbb{B}$ is generic with respect to Zariski topology, i.e., open
and dense. Specifically, its complement is the zero set of the polynomial

$P(X)=\sum_{1\leq i<j\leq n}~{}\sum_{1\leq k<l\leq
d}(X_{i,k}-X_{j,k})^{2}(X_{i,l}-X_{j,l})^{2}.$

  3. (3)

For an invertible matrix $A\in\mathbb{R}^{d\times d}$,

$\mathbb{B}\cdot A^{-1}\subset\mathfrak{S}_{n}(A)^{c}.$

Hence almost every matrix (w.r.t. Lebesgue measure) $X\in{\mathbb{R}}^{n\times
d}$ is _not_ separated by $A$.

Proof

(1) We need to show that any matrix $X$ that on some columns $k$ and $l$ has
distinct elements on same row positions is not separated by $I_{d}$. Indeed if
$X$ is such a matrix, let $Y$ denote a copy of $X$ except on those 4 entries
where we set

$Y_{i,k}=X_{j,k}~{}~{},~{}~{}Y_{j,k}=X_{i,k}~{}~{},~{}~{}Y_{i,l}=X_{i,l}~{}~{},~{}~{}Y_{j,l}=X_{j,l}.$

Note $X\not\sim Y$ yet $\downarrow(X)=\downarrow(Y)$. Hence such matrices are
not separated by $I_{d}$.

(2) By negation, the complement of $\mathbb{B}$ is given by

$\mathbb{B}^{c}=\\{X\in\mathbb{R}^{n\times d}~{},~{}\forall
i,j,k,l~{},~{}1\leq i<j\leq n,1\leq k<l\leq
d\&(X_{i,k}=X_{j,k}~{}or~{}X_{i,l}=X_{j,l})\\}$

This shows $\mathbb{B}^{c}$ is the zero set of polynomial $P$ as claimed. Thus
$\mathbb{B}^{c}$ is a closed Zariski set. Its complement is generic with
respect to the Zariski topology since $\mathbb{B}^{c}\neq\mathbb{R}^{n\times
d}$.

(3) The inclusion is immediate. Density claim follows from this inclusion.
$\Box$

On the other hand, extending the identity matrix by only one column produces
an almost universal key:

######  Proposition 3.7.

Assume $d\geq 2$ and $n\geq 3$.

Let $a\in\mathbb{R}^{d}$ be a vector with non-zero entries, i.e.,
$\prod_{i=1}^{d}a_{i}\neq 0$. Let
$A=\left[\begin{array}[]{ccc}I_{d}&|&a\end{array}\right]\in\mathbb{R}^{d\times(d+1)}$
be a key. Then $\mathfrak{S}_{n}(A)$ is generic with respect to the Zariski
topology (i.e., open and dense), however
$\mathfrak{S}_{n}(A)\neq\mathbb{R}^{n\times d}$. In particular, its complement
$\mathfrak{S}_{n}(A)^{c}:={\mathbb{R}}^{n\times
d}\setminus\mathfrak{S}_{n}(A)$ is non-empty but has Lebesgue measure zero.
Thus almost every matrix $X\in\mathbb{R}^{n\times d}$ is separated by $A$.

Proof

First we show that $\mathfrak{S}_{n}(A)\neq\mathbb{R}^{n\times d}$. Consider
the matrices $X,Y\in\mathbb{R}^{n\times d}$ full of zeros except for the 3x2
top left corner where:

$X_{1,1}=Y_{1,1}=\frac{1}{a_{1}}~{}~{},~{}~{}X_{2,1}=Y_{2,1}=-\frac{1}{a_{1}}~{}~{},~{}~{}X_{3,1}=Y_{3,1}=0$
$X_{1,2}=Y_{3,2}=-\frac{1}{a_{2}}~{}~{},~{}~{}X_{2,2}=Y_{1,2}=0~{}~{},~{}~{}X_{3,2}=Y_{2,2}=\frac{1}{a_{2}}$

Clearly $\beta_{A}(X)=\beta_{A}(Y)$ (the two left columns and the last column
contain $1$, $0$ repeated $n-2$ times and $-1$) and yet $X\not\sim Y$.

Next we show that $\mathfrak{S}_{n}(A)^{c}$ is included in a finite union of
linear spaces each of positive codimension. This proves the clam.

To simplify notation we introduce the following two operators. Let
$\Pi,\Pi_{0},\Pi_{1},\cdots,\Pi_{d}\in{\mathcal{S}}_{n}$ denote permutation
matrices of size n. For $X\in\mathbb{R}^{n\times d}$ denote by
$x_{1},\ldots,x_{d}$ its columns. Thus
$X=\left[x_{1}|x_{2}|\cdots|x_{d}\right]$.

$L_{\Pi_{0},\Pi_{1},\ldots,\Pi_{d}}:\mathbb{R}^{n\times
d}\rightarrow\mathbb{R}^{d}~{}~{},~{}~{}L_{\Pi_{0},\Pi_{1},\ldots,\Pi_{d}}X=\Pi_{0}Xa-(a_{1}\Pi_{1}x_{1}+\cdots
a_{d}\Pi_{d}x_{d})$

and

$M_{\Pi,\Pi_{1},\ldots,\Pi_{d}}:\mathbb{R}^{n\times
d}\rightarrow\mathbb{R}^{n\times
d}~{}~{},~{}~{}M_{\Pi,\Pi_{1},\ldots,\Pi_{d}}X=\Pi
X-\left[\begin{array}[]{ccccc}\Pi_{1}x_{1}&|&\cdots&|&\Pi_{d}x_{d}\end{array}\right]$

A matrix $X\in{\mathbb{R}}^{n\times d}$, $X=[x_{1}|\cdots|x_{d}]$ is not
separated by $A=[I_{d}|a]\in\mathbb{R}^{d\times(d+1)}$, i.e.,
$X\in\mathfrak{S}_{n}(A)^{c}$ if there are permutation matrices
$\Pi_{1},\ldots,\Pi_{d}\in{\mathcal{S}}_{n}$ such that the matrix
$Y=[\Pi_{1}x_{1}|\cdots|\Pi_{d}x_{d}]$ satisfies:

$Y\not\sim X~{}~{}and~{}~{}\downarrow(X\cdot a)=\downarrow(Y\cdot a)$

This is equivalent to say:

$\exists\Pi_{0}\in{\mathcal{S}}_{n}~{}~{},~{}~{}\Pi_{0}Xa-(a_{1}\Pi_{1}x_{1}+\cdots+a_{d}\Pi_{d}x_{d})=0$
$\forall\Pi\in{\mathcal{S}}_{n}\exists
k\in[d]~{}~{},~{}~{}(\Pi-\Pi_{k})x_{k}\neq 0$

Hence

$\mathfrak{S}_{n}(A)^{c}=\bigcup_{(\Pi_{0},\Pi_{1},\ldots,\Pi_{d})\in{\mathcal{S}}_{n}^{d+1}}\left(\ker\,L_{\Pi_{0},\Pi_{1},\ldots,\Pi_{d}}\setminus\left(\bigcup_{\Pi\in{\mathcal{S}}_{n}}\ker\,M_{\Pi,\Pi_{1},\ldots,\Pi_{d}}\right)\right)$

Let $\Delta$ denote the diagonal in ${\mathcal{S}}_{n}^{d+1}$,

$\Delta=\\{(\Pi_{0},\Pi_{1},\Pi_{2},\ldots,\Pi_{d})\in{\mathcal{S}}_{n}^{d+1}~{}~{},~{}~{}\Pi_{1}=\Pi_{2}=\cdots=\Pi_{d}\\}$

parametrized by the first two permutation matrices. For any
$(\Pi_{0},\Pi_{1},\ldots,\Pi_{d})\in\Delta$, we have
$\ker\,M_{\Pi_{1},\Pi_{1},\ldots,\Pi_{d}}={\mathbb{R}}^{n\times d}$. Thus

$\ker\,L_{\Pi_{0},\Pi_{1},\ldots,\Pi_{d}}\setminus\left(\bigcup_{\Pi\in{\mathcal{S}}_{n}}\ker\,M_{\Pi,\Pi_{1},\ldots,\Pi_{d}}\right)=\emptyset$

It follows:

$\mathfrak{S}_{n}(A)^{c}=\bigcup_{(\Pi_{0},\Pi_{1},\ldots,\Pi_{d})\in{\mathcal{S}}_{n}^{d+1}\setminus\Delta}\left(\ker\,L_{\Pi_{0},\Pi_{1},\ldots,\Pi_{d}}\setminus\left(\bigcup_{\Pi\in{\mathcal{S}}_{n}}\ker\,M_{\Pi,\Pi_{1},\ldots,\Pi_{d}}\right)\right)$

Consider now
$(\Pi_{0},\Pi_{1},\ldots,\Pi_{d})\in{\mathcal{S}}_{n}^{d+1}\setminus\Delta$.
Then

$\ker\,L_{\Pi_{0},\Pi_{1},\ldots,\Pi_{d}}=\ker\,L_{I,\Pi_{0}^{-1}\Pi_{1},\ldots,\Pi_{0}^{-1}\Pi_{d}}$

Hence there is $k\in[d]$ so that $\Pi_{0}^{-1}\Pi_{k}\neq I$. Choose
$x_{k}\in\mathbb{R}^{n}$ so that $(\Pi_{0}^{-1}\Pi_{k})x_{k}\neq x_{k}$. Set
$x_{j}=0$ for $j\in[d]$, $j\neq k$ and consider the matrix
$X=[x_{1}|\cdots|x_{d}]$. Then
$L_{\Pi_{0},\Pi_{1},\ldots,\Pi_{d}}X=a_{k}(\Pi_{0}-\Pi_{k})x_{k}\neq 0$. This
show that $\ker\,L_{\Pi_{0},\Pi_{1},\ldots,\Pi_{d}}\neq{\mathbb{R}}^{n\times
d}$ and hence it is a subspace of positive codimension. We obtain:

$\mathfrak{S}_{n}(A)^{c}\subset\bigcup_{(\Pi_{0},\Pi_{1},\ldots,\Pi_{d})\in{\mathcal{S}}_{n}^{d+1}\setminus\Delta}\ker\,L_{\Pi_{0},\Pi_{1},\ldots,\Pi_{d}}$

This shows that $\mathfrak{S}_{n}(A)^{c}$ is included in a finite union of
proper subspaces of ${\mathbb{R}}^{n\times d}$ which in turn is a closed set
with respect to the Zariski topology of empty interior. This ends the proof of
this result. $\Box$

The next result provides a characterization of the set $\mathfrak{S}_{n}(A)$.
To do so we need to introduce additional notation that extends the operators
$L_{\Pi_{0},\ldots,\Pi_{d}}$ and $M_{\Pi,\ldots,\Pi_{d}}$ defined in the proof
of Proposition 3.7. For $E_{1},E_{2},\ldots,E_{d}\in\mathbb{R}^{n\times n}$
and $b\in\mathbb{R}^{d}$, with $b=\left(b_{1},b_{2},\cdots,b_{d}\right)^{T}$
define

(3.7) $L_{E_{1},E_{2},\ldots,E_{d};b}:\mathbb{R}^{n\times
d}\rightarrow\mathbb{R}^{n}~{}~{},~{}~{}X=\left[x_{1}|x_{2}|\cdots|x_{d}\right]\mapsto
L_{E_{1},E_{2},\ldots,E_{d};b}(X)=b_{1}E_{1}x_{1}+\cdots+b_{d}E_{d}x_{d}.$

######  Proposition 3.8.

Fix $a_{1},\ldots,a_{D-d}\in\mathbb{R}^{d}$ and consider the key
$A=\left[I_{d}|a_{1}|\cdots|a_{D-d}\right]\in\mathbb{R}^{d\times D}$. Let
$X=\left[x_{1}|x_{2}|\cdots|x_{d}\right]$.

  1. (1)

$X\in\mathfrak{S}_{n}(A)^{c}:={\mathbb{R}}^{n\times
d}\setminus\mathfrak{S}_{n}(A)$ if and only if there are
$\Pi_{1},\Pi_{2},\ldots,\Pi_{d},\Xi_{1},\ldots,\Xi_{D-d}\in{\mathcal{S}}_{n}$
such that:

    1. (a)

$\forall j\in[D-d]$,
$\left[(\Xi_{j}-\Pi_{1})x_{1},\ldots,(\Xi_{j})-\Pi_{d})x_{d}\right]a_{j}=0$

    2. (b)

$\forall\Pi\in{\mathcal{S}}_{n}~{}\exists k\in[d]$ so that
$(\Pi_{k}-\Pi)x_{k}\neq 0$.

  2. (2)

The following hold true:

(3.8)
$\mathfrak{S}_{n}(A)=\bigcap_{\begin{array}[]{c}\Pi_{1},\ldots,\Pi_{d}\in{\mathcal{S}}_{n}\\\
\Xi_{1},\ldots,\Xi_{D-d}\in{\mathcal{S}}_{n}\end{array}}\left[\bigcup_{j=1}^{D-d}\left(\ker\,\
L_{\Xi_{j}-\Pi_{1},\ldots,\Xi_{j}-\Pi_{d};a_{j}}\right)^{c}\bigcup\bigcup_{\Pi\in{\mathcal{S}}_{n}}\bigcap_{k=1}^{d}\ker\,L_{\Pi_{1}-\Pi,\ldots,\Pi_{d}-\Pi;\delta_{k}}\right]$

and

(3.9)
$\mathfrak{S}_{n}(A)^{c}=\bigcup_{\begin{array}[]{c}\Pi_{1},\ldots,\Pi_{d}\in{\mathcal{S}}_{n}\\\
\Xi_{1},\ldots,\Xi_{D-d}\in{\mathcal{S}}_{n}\end{array}}\left(\bigcap_{j=1}^{D-d}\ker\,\
L_{\Xi_{j}-\Pi_{1},\ldots,\Xi_{j}-\Pi_{d};a_{j}}\right)\bigcap\left(\bigcup_{\Pi\in{\mathcal{S}}_{n}}\bigcap_{k=1}^{d}\ker\,L_{\Pi_{1}-\Pi,\ldots,\Pi_{d}-\Pi;\delta_{k}}\right)^{c}$

where $\delta_{k}=(0,\ldots,0,1,0,\ldots,0)^{T}$ has only one 1 on the
$k^{th}$ position.

Proof

The proof is a consequence of linear algebra analysis of map $\beta_{A}$.

(1) Assume $X$ is not separated by $A$. Then there is
$Y\in{\mathbb{R}}^{n\times d}$ so that $\beta_{A}(X)=\beta_{A}(Y)$ yet
$X\not\sim Y$.

Let $Y=[y_{1}|\cdots|y_{d}]$. Then $\beta_{A}(X)=\beta_{A}(Y)$ implies that
there are permutation matrices
$\Pi_{1},\ldots,\Pi_{d},\Xi_{1},\ldots,\Xi_{D-d}\in{\mathcal{S}}_{n}$ so that:

$y_{1}=\Pi_{1}x_{1},\ldots,y_{d}=\Pi_{d}x_{d},Ya_{1}=\Xi_{1}Xa_{1},\ldots,Ya_{D-d}=\Xi_{D-d}Xa_{D-d}$

Substituting the expressions for $y_{1},\ldots,y_{d}$ provided by the first
$d$ equations into the latter $D-d$ equations, we obtain part 1.(a).

For same $Y$, the condition $X\not\sim Y$ implies that for every
$\Pi\in{\mathcal{S}}_{n}$, $Y-\Pi X\neq 0$. Thus part 1(b) is proved.

(3) Equation (3.9) is a transcription of part 1.

(2) Equation (3.8) follows from (3.9) by taking the complement. $\Box$

###  3.2. Construction of universal keys

In this subsection we construct universal keys. Proposition 3.8 provides us
with an algorithm to check whether a key $A$ is universal. Unfortunately the
algorithm has an exponential complexity in data size.

If the key $A\in\mathbb{R}^{d\times D}$ is universal then $A$ must have full
rank. Therefore there are permutation matrix $L\in\SS_{D}$ and invertible
$T\in GL(d,\mathbb{R})$ so that $A=T\left[I_{d}~{}\tilde{A}\right]L$, with
$\tilde{A}\in\mathbb{R}^{d\times(D-d)}$. Proposition 3.5 shows that $A$ is a
universal key if and only if $\left[I_{d}~{}\tilde{A}\right]$ is a universal
key. This observation allows us to prove the main result of this subsection
stated earlier as part b of Theorem 2.1. Recall a set of vectors
$\\{f_{1},\ldots,f_{m}\\}$ in a linear space $V$ of finite dimension $n\leq m$
is called a full spark frame if any subset of $n$ vectors is linearly
independent. See [1, 26] for more information and explicit constructions of
full spark frames.

######  Theorem 3.9.

Consider the metric space $({\widehat{{\mathbb{R}}^{n\times d}}},d)$. Set
$D=1+(d-1)n!$ and let $A\in\mathbb{R}^{d\times D}$ be a matrix whose columns
form a full spark frame, i.e., any subset of $d$ columns is linearly
independent. Then the key $A$ is universal and the induced map
${\hat{\beta}}_{A}:{\widehat{{\mathbb{R}}^{n\times
d}}}\rightarrow\mathbb{R}^{n\times D}$, $X\mapsto\beta_{A}(X)=\downarrow(XA)$
is injective. Furthermore, ${\hat{\beta}}_{A}$ is bi-Lipschitz, with estimates
of the bi-Lipschitz constants $a_{0}=\min_{J\subset[D],|J|=d}s_{d}(A[J])$ and
$b_{0}=s_{1}(A)$, where $s_{1}(A)$ denotes the largest singular value of $A$,
$A[J]$ denotes the submatrix of $A$ formed by columns indexed by $J$, and
$s_{d}(A[J])$ denotes the $d^{th}$ singular value (in this case, the smallest)
of $A[J]$. Specifically, for any $X,Y\in{\mathbb{R}}^{n\times d}$,

(3.10) $a_{0}d(\hat{X},\hat{Y})\leq{\|\beta_{A}(X)-\beta_{A}(Y)\|}\leq
b_{0}d(\hat{X},\hat{Y})$

where all norms are Frobenius norms.

Proof

Let $a_{1},\ldots,a_{D}$ denote the columns of $A$, $A=[a_{1}|\cdots|a_{D}]$.

Fix $X,Y\in{\mathbb{R}}^{n\times d}$ two matrices. Then there are permutation
matrices
$P_{0},\Pi_{1},\ldots,\Pi_{D},\Xi_{1},\ldots,\Xi_{D}\in{\mathcal{S}}_{n}$ so
that $d(\hat{X},\hat{Y})={\|P_{0}X-Y\|}$ and

$\beta_{A}(X)=\left[\begin{array}[]{ccccc}\Pi_{1}Xa_{1}&|&\cdots&|&\Pi_{D}Xa_{D}\end{array}\right]~{},~{}\beta_{A}(Y)=\left[\begin{array}[]{ccccc}\Xi_{1}Ya_{1}&|&\cdots&|&\Xi_{D}Ya_{D}\end{array}\right]$

Thus

(3.11)
${\|\beta_{A}(X)-\beta_{A}(Y)\|}^{2}=\sum_{k=1}^{D}{\|(\Pi_{k}X-\Xi_{k}Y)a_{k}\|}_{2}^{2}=\sum_{k=1}^{D}{\|(\Xi_{k}^{T}\Pi_{k}X-Y)a_{k}\|}_{2}^{2}$

Permutations $\Pi_{k}$ and $\Xi_{k}$ satisfy the optimality condition:
${\|\Pi_{k}Xa_{k}-\Xi_{k}Ya_{k}\|}_{2}=\min_{P\in{\mathcal{S}}_{n}}{\|PXa_{k}-Ya_{k}\|}_{2}$.
Hence ${\|\Pi_{k}Xa_{k}-\Xi_{k}Ya_{k}\|}_{2}\leq{\|P_{0}Xa_{k}-Ya_{k}\|}_{2}$.
Therefore:

(3.12)
${\|\beta_{A}(X)-\beta_{A}(Y)\|}^{2}\leq\sum_{k=1}^{D}{\|P_{0}Xa_{k}-Ya_{k}\|}_{2}^{2}={\|(PX-Y)A\|}^{2}\leq
s_{1}(A)^{2}{\|P_{0}X-Y\|}^{2}$

from where we obtain the upper bound in (3.10).

The lower bound in (3.10) follows from the pigeonhole principle similar to the
one employed in the proof of Theorem 2.1. In equation (3.11) there are
$D=1+(d-1)n!$ terms. Since only $n!$ permutations are distinct, there is a
permutation $Q$ that repeats at least $d$ times. Say
$J=\\{j_{1},j_{2},\ldots,j_{d}\\}\subset[D]$ is a set of indices so that
$\Xi_{j_{1}}^{T}\Pi_{j_{1}}=\cdots=\Xi_{j_{d}}^{T}\Pi_{j_{d}}=Q$. Then

$\displaystyle{\|\beta_{A}(X)-\beta_{A}(Y)\|}^{2}$
$\displaystyle\geq\sum_{k=1}^{d}{\|(\Xi_{j_{k}}^{T}\Pi_{j_{k}}X-Y)a_{j_{k}}\|}_{2}^{2}={\|(QX-Y)A[J]\|}^{2}$
$\displaystyle\geq s_{d}(A[J])^{2}{\|QX-Y\|}^{2}\geq
s_{d}(A[J])^{2}{\|P_{0}X-Y\|}^{2}\geq a_{0}^{2}d(\hat{X},\hat{Y})^{2}.$

The lower bound in (3.10) implies that
${\hat{\beta}}_{A}:{\widehat{{\mathbb{R}}^{n\times
d}}}\rightarrow\mathbb{R}^{n\times D}$ is injective and hence $A$ is a
universal key. This ends the proof of Theorem 3.9. $\Box$

###  3.3. Bi-Lipschitz properties of universal keys

In this subsection we prove that any universal key defines a bi-Lipschitz
encoding map, regardless of $D$.

######  Theorem 3.10.

Assume the key $A\in\mathbb{R}^{d\times D}$ is universal, i.e., the induced
map ${\hat{\beta}}_{A}:{\widehat{{\mathbb{R}}^{n\times
d}}}\rightarrow\mathbb{R}^{n\times D}$, $X\mapsto\beta_{A}(X)=\downarrow(XA)$
is injective. Then ${\hat{\beta}}_{A}$ is bi-Lipschitz, that is, there are
constants $a_{0}>0$ and $b_{0}>0$ so that for all $X,Y\in{\mathbb{R}}^{n\times
d}$,

(3.13) $a_{0}\,d(\hat{X},\hat{Y})\leq{\|\beta_{A}(X)-\beta_{A}(Y)\|}\leq
b_{0}\,d(\hat{X},\hat{Y})$

where all are Frobenius norms. Furthermore, an estimate for $b_{0}$ is
provided by the largest singular value of $A$, $b_{0}=s_{1}(A)$.

Proof

The upper bound in (3.13) follows as in the proof of Theorem 3.9, from
equations (3.11) and (3.12). Notice that no property is assumed in order to
obtain the upper Lipschitz bound.

The lower bound in (3.13) is more difficult. It is shown by contradiction
following the strategy utilized in the Complex Phase Retrieval problem [6].

Assume $\inf_{X\not\sim
Y}\frac{{\|\beta_{A}(X)-\beta_{A}(Y)\|}_{2}^{2}}{d(\hat{X},\hat{Y})^{2}}=0$.

Step 1: Reduction to local analysis. Since
$d(\hat{tX},\hat{tY})=t\,d(\hat{X},\hat{Y})$ for all $t>0$, the quotient
$\frac{{\|\beta_{A}(X)-\beta_{A}(Y)\|}_{2}}{d(\hat{X},\hat{Y})}$ is scale
invariant. Therefore, there are sequences $(X^{t})_{t},(Y^{t})_{t}$ with
${\|Y^{t}\|}\leq{\|X^{t}\|}=1$ and $d(\hat{X^{t}},\hat{Y^{t}})>0$ so that
$\lim_{t\rightarrow\infty}\frac{{\|\beta_{A}(X^{t})-\beta_{A}(Y^{t})\|}_{2}}{d(\hat{X^{t}},\hat{Y^{t}})}=0$.
By compactness of the closed unit ball, one can extract convergence
subsequences. For easiness of notation, assume $(X^{t})_{t},(Y^{t})_{t}$ are
these subsequences. Let ${X^{\infty}}=\lim_{t}X^{t}$ and
${Y^{\infty}}=\lim_{t}Y^{t}$ denote their limits. Notice
$\lim_{t}{\|\beta_{A}(X^{t})-\beta_{A}(Y^{t})\|}_{2}=0$. This implies
${\|\beta_{A}({X^{\infty}})-\beta_{A}({Y^{\infty}})\|}=0$ and thus
$\beta_{A}({X^{\infty}})=\beta_{A}({Y^{\infty}})$. Since $\widehat{\beta_{A}}$
is assumed injective, it follows that
$\widehat{{X^{\infty}}}=\widehat{{Y^{\infty}}}$.

This means that, if the lower Lipschitz bound vanishes, then this is achieved
by vanishing of a local lower Lipschitz bound. To follow the terminology in
[6], the type I local lower Lipschitz bound vanishes at some
$Z_{0}\in{\mathbb{R}}^{n\times d}$, with ${\|Z_{0}\|}=1$:

(3.14) ${A}(Z_{0}):=\lim_{r\rightarrow
0}\inf_{\begin{array}[]{c}\hat{X}\neq\hat{Y}\\\ d(\hat{X},\hat{Z_{0}})<r\\\
d(\hat{Y},\hat{Z_{0}})<r\end{array}}\frac{{\|\beta_{A}(X)-\beta_{A}(Y)\|}_{2}^{2}}{d(\hat{X},\hat{Y})^{2}}=0.$

Note that, in general, the infimum of the type I local lower Lipschitz bound
over the unit sphere may be strictly larger than the global lower Lipschitz
bound (see Theorems 2.1 and Theorem 2.2 in [6] and Theorem 4.3 in [5]). The
compactness argument forces the local lower Lipschitz bound to vanish when the
global lower bound vanishes.

Step 2. Local Linearization. The following stability subgroups of
${\mathcal{S}}_{n}$ play an important role:

$G=\\{P\in{\mathcal{S}}_{n}~{}:~{}PZ_{0}=Z_{0}\\}~{}~{},~{}~{}H_{j}=\\{P\in{\mathcal{S}}_{n}~{}:~{}PZ_{0}a_{j}=Z_{0}a_{j}\\}~{},~{}1\leq
j\leq D.$

Obviously $I_{n}\in G\subset H_{j}\subset{\mathcal{S}}_{n}$, for every
$j\in[D]$. The group $G$ is the stabilizer of $Z_{0}$, whereas $H_{j}$ is the
stabilizer of $Z_{0}a_{j}$. Let
$\delta_{0}=\min_{P\in{\mathcal{S}}_{n}\setminus G}{\|(I_{n}-P)Z_{0}\|}$
denote the smallest variation of $Z_{0}$ under row permutations. Note
$\delta_{0}>0$ by the definition of $G$.

Consider $X=Z_{0}+U$ and $Y=Z_{0}+V$ where $U,V\in{\mathbb{R}}^{n\times d}$
are “aligned” in the sense that $d(\hat{X},\hat{Y})={\|U-V\|}$. This property
requires that ${\|U-V\|}\leq{\|PX-Y\|}$, for every $P\in{\mathcal{S}}_{n}$.
Next result replaces equivalently this condition by requirements involving
$(U,V)$ and the group $G$ only.

######  Lemma 3.11.

Assume ${\|U\|},{\|V\|}<\frac{1}{4}\delta_{0}$, where
$\delta_{0}=\min_{P\in{\mathcal{S}}_{n}\setminus G}{\|(I_{n}-P)Z_{0}\|}$. Let
$X=Z_{0}+U$, $Y=Z_{0}+V$. Then:

  1. (1)

$d(\hat{X},\hat{Z_{0}})={\|U\|}$ and $d(\hat{Y},\hat{Z_{0}})={\|V\|}$.

  2. (2)

$d(\hat{X},\hat{Y})=\min_{P\in G}{\|U-PV\|}=\min_{P\in G}{\|PU-V\|}$

  3. (3)

The following are equivalent:

    1. (a)

$d(\hat{X},\hat{Y})={\|U-V\|}$.

    2. (b)

For every $P\in G$, ${\|U-V\|}\leq{\|PU-V\|}$.

    3. (c)

For every $P\in G$, ${\langle U,V\rangle}\geq{\langle PU,V\rangle}$.

Proof of Lemma 3.11. (1)

Note that is $U=0$ then the claim follows. Assume $U\neq 0$. Then

$d(\hat{X},\hat{Z_{0}})=\min_{P\in{\mathcal{S}}_{n}}{\|X-PZ_{0}\|}=\min_{P\in{\mathcal{S}}_{n}}{\|(I_{n}-P)Z_{0}+U\|}\leq{\|U\|}$

On the other hand, assume the minimum is achieved for a permutation
$P_{0}\in{\mathcal{S}}_{n}$. If $P_{0}\in G$ then
$d(\hat{X},\hat{Z_{0}})={\|(I_{n}-P_{0})Z_{0}+U\|}={\|U\|}$. If $P_{0}\not\in
G$ then

$d(\hat{X},\hat{Z_{0}})\geq{\|(I_{n}-P_{0})Z_{0}\|}-{\|U\|}>\frac{3\delta_{0}}{4}>{\|U\|}\geq
d(\hat{X},\hat{Z_{0}})$

which yields a contradiction. Hence $d(\hat{X},\hat{Z_{0}})={\|U\|}$.
Similarly, one shows $d(\hat{X},\hat{Z_{0}})={\|V\|}$.

(2) Clearly

$d(\hat{X},\hat{Y})=\min_{P\in{\mathcal{S}}_{n}}{\|PX-Y\|}\leq\min_{P\in
G}{\|PX-Y\|}=\min_{P\in G}{\|PU-V\|}$

On the other hand, for $P\in{\mathcal{S}}_{n}\setminus G$ and $Q\in G$,

${\|PX-Y\|}={\|(P-I_{n})Z_{0}+PU-V\|}\geq{\|(I_{n}-P)Z_{0}\|}-{\|U\|}-{\|V\|}\geq$
$\geq\delta_{0}-2{\|U\|}-2{\|V\|}+{\|QU-V\|}\geq\min_{Q\in G}{\|QU-V\|}\geq
d(\hat{X},\hat{Y}).$

(3)

(a)$\Rightarrow$(b).

If $d(\hat{X},\hat{Y})={\|U-V\|}$ then

${\|U-V\|}\leq{\|PX-Y\|}={\|(P-I_{n})Z_{0}+PU-V\|}~{}~{},~{}~{}\forall
P\in{\mathcal{S}}_{n}.$

In particular, for $P\in G$, $(P-I_{n})Z_{0}=0$ and the above inequality
reduces to (b).

(b)$\Rightarrow$(a).

Assume (b). For $P\in G$,

${\|U-V\|}={\|X-Y\|}\leq{\|PU-V\|}={\|PX-Y\|}.$

For $P\in{\mathcal{S}}_{n}\setminus G$,

${\|PX-Y\|}={\|(P-I_{n})Z_{0}+PU-V\|}\geq{\|(I_{n}-P)Z_{0}\|}-{\|U\|}-{\|V\|}\geq$
$\geq\delta_{0}-2{\|U\|}-2{\|V\|}+{\|U-V\|}\geq{\|U-V\|}={\|X-Y\|}.$

This shows $d(\hat{X},\hat{Y})={\|X-Y\|}={\|U-V\|}$.

(b)$\Longleftrightarrow$(c). This is immediate after squaring (b) and
simplifying the terms.

$\Box$

Consider now sequences $(\hat{X^{t}})_{t},(\hat{Y^{t}})_{t}$ that converge to
$\hat{Z_{0}}$ and achieve lower bound 0 as in (3.14). Choose representatives
$X_{t}$ and $Y_{t}$ in their equivalence classes that satisfy the hypothesis
of Lemma 3.11 so that $X_{t}=Z_{0}+U_{t}$, $Y_{t}=Z_{0}+V_{t}$,
${\|U_{t}\|},{\|V_{y}\|}<\frac{1}{4}\delta_{0}$,
$d(\hat{X_{t}},\hat{Z_{0}})={\|U_{t}\|}$,
$d(\hat{Y_{t}},\hat{Z_{0}})={\|V_{t}\|}$ and
$d(\hat{X_{t}},\hat{Y_{t}})={\|U_{t}-V_{t}\|}>0$. With
$A=[a_{1}|\cdots|a_{D}]$ we obtain:

${\|\beta_{A}(X_{t})-\beta_{A}(Y_{t})\|}_{2}^{2}=\sum_{j=1}^{D}{\|\downarrow(X_{t}a_{j})-\downarrow(Y_{t}a_{j})\|}_{2}^{2}=\sum_{j=1}^{D}{\|(Z_{0}+U_{t})a_{j}-\Pi_{j,t}(Z_{0}+V_{t})a_{j}\|}_{2}^{2}$

for some $\Pi_{j,t}\in{\mathcal{S}}_{n}$. In fact $\Pi_{j,t}\in argmin_{\Pi\in
H_{j}}{\|U_{t}-\Pi V_{t})a_{j}\|}_{2}$. Pass to sub-sequences (that will be
indexed by $t$ for an easier notation) so that $\Pi_{j,t}=\Pi_{j}$ for some
$\Pi_{j}\in{\mathcal{S}}_{n}$. Thus

${\|\beta_{A}(X_{t})-\beta_{A}(Y_{t})\|}_{2}^{2}=\sum_{j=1}^{D}{\|(I_{n}-\Pi_{j})Z_{0}a_{j}+(U_{t}-\Pi_{j}V_{t})a_{j}\|}_{2}^{2}$

Since the above sequence must converge to $0$ as $t\rightarrow\infty$, while
$U_{t},V_{t}\rightarrow 0$, it follows that necessarily $\Pi_{j}\in H_{j}$ and
the expressions simplify to

${\|\beta_{A}(X_{t})-\beta_{A}(Y_{t})\|}_{2}^{2}=\sum_{j=1}^{D}{\|(U_{t}-\Pi_{j}V_{t})a_{j}\|}_{2}^{2}$

Thus equation (3.14) implies that for every $j\in[D]$,

(3.15)
$\lim_{t\rightarrow\infty}\frac{{\|(U_{t}-\Pi_{j}V_{t})a_{j}\|}_{2}^{2}}{{\|U_{t}-V_{t}\|}^{2}}=0$

where $\Pi_{j}\in H_{j}$, ${\|U_{t}\|},{\|V_{t}\|}\rightarrow 0$, and
$U_{t},V_{t}$ are aligned so that ${\langle U_{t},V_{t}\rangle}\geq{\langle
PU_{t},V_{t}\rangle}$ for every $P\in G$. Equivalently, relation (3.14) can be
restated as:

(3.16) $\inf_{\begin{array}[]{c}U,V\in{\mathbb{R}}^{n\times d}\\\ s.t.\\\
U\neq V\\\ {\langle U,V\rangle}\geq{\langle PU,V\rangle},\forall P\in
G\end{array}}\frac{\sum_{j=1}^{D}{\|(U-\Pi_{j}V)a_{j}\|}_{2}^{2}}{{\|U-V\|}^{2}}=0$

for some permutations $\Pi_{j}\in H_{j}$, $j\in[D]$. By Lemma 3.11 the
constraint in the optimization problem above implies ${\|U-V\|}=\min_{P\in
G}{\|U-PV\|}$. Hence (3.16) implies:

(3.17) $\inf_{\begin{array}[]{c}U,V\in{\mathbb{R}}^{n\times d}\\\ s.t.\\\
U\neq PV,\forall P\in G\end{array}}\max_{P\in
G}\frac{\sum_{j=1}^{D}{\|(U-\Pi_{j}V)a_{j}\|}_{2}^{2}}{{\|U-PV\|}^{2}}=0$

for same permutation matrices $\Pi_{j}$’s. While the above optimization
problem seems a relaxation of (3.16), in fact (3.17) implies (3.16) with a
possibly change of permutation matrices $\Pi_{j}$, but remaining still in
$H_{j}$.

Step 3. Existence of a Minimizer.

The optimization problem (3.16) is a Quadratically Constrained Ratio of
Quadratics (QCRQ) optimization problem. A significant number of papers have
been published on this topic [7, 8]. In particular, [3] presents a formal
setup for analysis of QCRQ problems. Our interest is to utilize some of these
techniques in order to establish the existence of a minimizer for (3.16) or
(3.17). Specifically we show:

######  Lemma 3.12.

Assume the key $A$ has linearly independent rows (equivalently, the columns of
$A$ form a frame for $\mathbb{R}^{d}$) and the lower Lipschitz bound of
${\hat{\beta}}_{A}$ is $0$. Then there are
$\tilde{U},\tilde{V}\in{\mathbb{R}}^{n\times d}$ so that:

  1. (1)

$\tilde{U}\neq P\tilde{V}$, for every $P\in G$;

  2. (2)

For every $j\in[D]$, $(\tilde{U}-\Pi_{j}\tilde{V})a_{j}=0$.

Proof of Lemma 3.12

We start with the formulation (3.17). Therefore there are sequences
$(U_{t},V_{t})_{t\geq 1}$ so that $U_{t}\neq PV_{t}$ for any $P\in G,t\geq 1$,
and yet for any $P\in G$,

$\lim_{t\rightarrow\infty}\frac{\sum_{j=1}^{D}{\|(U_{t}-\Pi_{j}V_{t})a_{j}\|}_{2}^{2}}{{\|U_{t}-PV_{t}\|}^{2}}=0.$

Let $E=\\{(U,V)\in{\mathbb{R}}^{n\times d}\times{\mathbb{R}}^{n\times
d}~{},~{}(U-\Pi_{j})V)a_{j}=0~{},~{}\forall j\in[D]\\}$ denote the null space
of the linear operator

$T:{\mathbb{R}}^{n\times d}\times{\mathbb{R}}^{n\times
d}\rightarrow\mathbb{R}^{D}~{},~{}(U,V)\mapsto\left[\begin{array}[]{ccccc}(U-\Pi_{1}V)a_{1}&|&\cdots&|&(U-\Pi_{D}V)a_{D}\end{array}\right],$

associated to the numerator of the above quotient. Let
$F_{P}=\\{(U,V)\in{\mathbb{R}}^{n\times d}\times{\mathbb{R}}^{n\times
d}~{},~{}U-PV=0\\}$ be the null space of the linear operator

$R_{P}:{\mathbb{R}}^{n\times d}\times{\mathbb{R}}^{n\times
d}\rightarrow{\mathbb{R}}^{n\times d}~{},~{}(U,V)\mapsto U-PV.$

A consequence of (3.17) is that for every $P\in G$, $E\setminus
F_{P}\neq\emptyset$. In particular, $F_{p}\cap E$ is a subspace of $E$ of
positive codimension. Using the Baire category theorem (or more elementary
linear algebra arguments), we conclude that

$E\setminus\left(\cup_{P\in G}F_{P}\right)\neq\emptyset.$

Let $(\tilde{U},\tilde{V})\in E\setminus\left(\cup_{P\in G}F_{P}\right)$. This
pair satisfies the conclusions of Lemma 3.12.

$\Box$

Step 4. Contradiction with the universality property of the key.

So far we obtained that if the lower Lipschitz bound of ${\hat{\beta}}_{A}$
vanishes than there are $Z_{0},\tilde{U},\tilde{V}\in{\mathbb{R}}^{n\times d}$
with $Z_{0}\neq 0$ and $\tilde{U}\neq P\tilde{V}$, for all $P\in G$ that
satisfy the conclusions of Lemma 3.12. Notice ${\langle
Z_{0},Z_{0}\rangle}={\langle PZ_{0},Z_{0}\rangle}$ for all $P\in G$ and
$(Z_{0}-\Pi_{j}Z_{0})a_{j}=0$ for all $j\in[D]$. Choose $s>0$ but small enough
so that $s{\|\tilde{U}\|},s{\|\tilde{V}\|}<\frac{1}{4}\delta_{0}$ with
$\delta_{0}=\min_{P\in{\mathcal{S}}_{n}\setminus G}{\|(I_{n}-P)Z_{0}\|}$. Let
$X=Z_{0}+s\tilde{U}$ and $Y=Z_{0}+s\tilde{V}$. Then Lemma 3.11 implies
$d(\hat{X},\hat{Y})=\min_{P\in G}{\|\tilde{U}-P\tilde{V}\|}>0$. Hence
$\hat{X}\neq\hat{Y}$. On the other hand, for every $j\in[D]$,
$Xa_{j}=\Pi_{j}Ya_{j}$. Thus
${\hat{\beta}}_{A}(\hat{X})={\hat{\beta}}_{A}(\hat{Y})$. Contradiction with
the assumption that ${\hat{\beta}}_{A}$ is injective.

This ends the proof of Theorem 3.10.

$\Box$

###  3.4. Dimension Reduction

Theorem 3.9 provides an Euclidean bi-Lipschitz embedding of very high
dimension, $D=1+(d-1)n!$. On the other hand, Theorem 3.10 shows that any
universal key $A\in\mathbb{R}^{d\times D}$ for
${\widehat{{\mathbb{R}}^{n\times d}}}$, and hence any injective map
$\hat{\beta}_{A}$ is bi-Lipschitz. In this subsection we show that any bi-
Lipschitz Euclidean embedding $\hat{\beta}_{A}:{\widehat{{\mathbb{R}}^{n\times
d}}}\rightarrow\mathbb{R}^{n\times D}$ with $D>2d$ can be further compressed
to a smaller dimension space $\mathbb{R}^{m}$ with $m=2nd$ thus yielding bi-
Lipschitz Euclidean embeddings of redundancy 2. This is shown in the next
result.

######  Theorem 3.13.

Assume $A\in\mathbb{R}^{d\times D}$ is a universal key for
${\widehat{{\mathbb{R}}^{n\times d}}}$ with $D\geq 2d$. Then, for $m\geq 2nd$,
a generic linear operator $B:\mathbb{R}^{n\times D}\rightarrow\mathbb{R}^{m}$
with respect to Zariski topology on $\mathbb{R}^{n\times D\times m}$, the map

(3.18) $\hat{\beta}_{A,B}:{\widehat{{\mathbb{R}}^{n\times
d}}}\rightarrow\mathbb{R}^{2nd}~{},~{}\hat{\beta}_{A,B}(\hat{X})=B\left(\hat{\beta}_{A}(\hat{X})\right)$

is bi-Lipschitz. In particular, almost every full-rank linear operator
$B:\mathbb{R}^{n\times D}\rightarrow\mathbb{R}^{2nd}$ produces such a bi-
Lipschitz map.

######  Remark 3.14.

The proof shows that, in fact, the complement set of linear operators $B$ that
produce bi-Lipschitz embeddings is included in the zero-set of a polynomial.

######  Remark 3.15.

Putting together Theorems 3.9, 3.10, 3.13 we obtain that the metric space
${\widehat{{\mathbb{R}}^{n\times d}}}$ admits a global bi-Lipschitz embedding
in the Euclidean space $\mathbb{R}^{2nd}$. This result is compatible with a
Whitney embedding theorem (see §1.3 in [19]) with the important caveat that
the Whitney embedding result applies to smooth manifolds, whereas here
${\widehat{{\mathbb{R}}^{n\times d}}}$ is merely a non-smooth algebraic
variety.

######  Remark 3.16.

These three theorems are summarized in part two of the Theorem 2.1 presented
in the first section.

######  Remark 3.17.

While the embedding dimension grows linearly in $nd$, in fact $m=2nd$, the
computational complexity of constructing ${\hat{\beta}}_{A,B}$ is NP due to
the $1+(d-1)n!$ intermediary dimension.

######  Remark 3.18.

As the proofs show, for $D\geq 1+(d-1)n!$, a generic $(A,B)$ with respect to
Zariski topology, $A\in\mathbb{R}^{d\times D}$ and linear map
$B:\mathbb{R}^{n\times D}\rightarrow\mathbb{R}^{2nd}$, produces a bi-Lipschitz
embedding $({\hat{\beta}}_{A,B},d)$ of ${\widehat{{\mathbb{R}}^{n\times d}}}$
into $(\mathbb{R}^{2nd},{\|\cdot\|}_{2})$.

Proof of Theorem 3.13

The proof follows a similar approach as in Theorem 3 of [20]. See also [13].

Without loss of generality, assume $m<nD$.

Notice $\beta_{A}:\mathbb{R}^{n\times d}\rightarrow\mathbb{R}^{n\times D}$ is
already homogeneous of degree 1 (with respect to positive scalars). Let
$\Delta:\mathbb{R}^{n\times d}\times\mathbb{R}^{n\times
d}\rightarrow\mathbb{R}^{n\times D}$ be defined by
$\Delta(X,Y)=\beta_{A}(X)-\beta_{A}(Y)$. Denote
$E=Ran(\Delta)=\\{\beta_{A}(X)-\beta_{A}(Y)~{},~{}X,Y\in\mathbb{R}^{n\times
d}\\}$.

Recall $A=[a_{1}|\cdots|a_{D}]$ is a notation for the columns of key $A$.
Notice that

$\Delta(X,Y)=\left[P_{1}Xa_{1}-Q_{1}Ya_{1}|\cdots|P_{D}Xa_{D}-Q_{D}Ya_{D}\right]$

for some $P_{1},\ldots,P_{D},Q_{1},\ldots,Q_{D}\in{\mathcal{S}}_{n}$, so that
for each $k\in[D]$, $P_{k},Q_{k}$ are permutations that sort monotone
decreasingly vectors $Xa_{k}$ and $Ya_{k}$, respectively. In particular,

$E\subset
F:=\bigcup_{\gamma\in({\mathcal{S}}_{n})^{2D}}F_{\gamma}~{}~{},~{}~{}F_{\gamma}:=Ran(L_{\gamma})$

where the $(n!)^{2D}$ linear operators $L_{\gamma}:\mathbb{R}^{n\times
d}\times\mathbb{R}^{n\times d}\rightarrow\mathbb{R}^{n\times D}$, are defined
by

$L_{\gamma}(X,Y)=\left[P_{1}Xa_{1}-Q_{1}Ya_{1}|\cdots|P_{D}Xa_{D}-Q_{D}Ya_{D}\right]$

when
$\gamma=(P_{1},\ldots,P_{D},Q_{1},\ldots,Q_{D})\in({\mathcal{S}}_{n})^{2D}$.

Claim: We claim that, for $m\geq 2nd$ and a generic linear operator
$B:\mathbb{R}^{n\times D}\rightarrow\mathbb{R}^{m}$ we have $\ker(B)\cap
F=\\{0\\}$. Such a generic linear operator has the kernel of dimension
$\dim(\ker(B))=nD-m\leq n(D-2d)$. It is therefore sufficient to show that, for
a generic subspace $V\subset\mathbb{R}^{n\times D}$ of dimension $r\leq
n(D-2d)$, for every $\gamma\in({\mathcal{S}}_{n})^{2D}$, $V\cap
F_{\gamma}=\\{0\\}$. This last claim follows from the observation
$\dim(F_{\gamma})\leq 2nd$.

We now show how this claim proves the Theorem. Let $B$ be such a linear map,
and let $\beta_{A,B}:\mathbb{R}^{n\times d}\rightarrow\mathbb{R}^{m}$ be the
map $\beta_{A,B}(X)=B(\downarrow(XA))$. Then $\beta_{A,B}(X)=\beta_{A,B}(Y)$
implies $\Delta(X,Y)=\beta_{A}(X)-\beta_{A}(Y)\in\ker(B)$. Thus
$\Delta(X,Y)=0$ which implies $\beta_{A}(X)=\beta_{A}(Y)$. Since
${\hat{\beta}}_{A}$ is injective on ${\widehat{{\mathbb{R}}^{n\times d}}}$ it
follows $\hat{X}=\hat{Y}$. Thus ${\hat{\beta}}_{A,B}$ is injective. On the
other hand, for each
$\gamma=(P_{1},\ldots,P_{D},Q_{1},\ldots,Q_{D})\in{\mathcal{S}}_{n}^{2D}$, the
restriction of $B$ to the linear space $Ran(L_{\gamma})$ is injective, and
thus bounded below as a linear map: there is $a_{\gamma}>0$ so that for every
$X,Y\in\mathbb{R}^{n\times d}$, ${\|B(L_{\gamma}(X,Y))\|}\geq
a_{\gamma}{\|L_{\gamma}(X,Y)\|}$. Let $a_{\infty}=\min_{\gamma}a_{\gamma}>0$.
Thus

${\|\beta_{A,B}(X)-\beta_{A,B}(Y)\|}={\|B(L_{\gamma_{0}}(X,Y))\|}\geq
a_{\infty}{\|L_{\gamma_{0}}(X,Y)\|}=a_{\infty}{\|\beta_{A}(X)-\beta_{A}(Y)\|}$

where $\gamma_{0}\in(S_{n})^{2D}$ is a particular $2D$-tuple of permutations.
This shows that $B{|}_{\beta_{A}(\mathbb{R}^{n\times
d})}:\beta_{A}(\mathbb{R}^{n\times d})\rightarrow\mathbb{R}^{m}$ is bi-
Lipschitz. By Theorem 3.10, the map ${\hat{\beta}}_{A}$ is bi-Lipschitz.
Therefore we get ${\hat{\beta}}_{A,B}$ is bi-Lipschitz as well. $\Box$

###  3.5. Proof of Corollary 1.3

(1) It is clear that any continuous $f$ induces a continuous
$\varphi:\beta(\mathbb{R}^{n\times d})\rightarrow\mathbb{R}$ via
$\varphi(\beta(X))=f(X)$. Furthermore, $F:=\beta(\mathbb{R}^{n\times
d})={\hat{\beta}}({\widehat{{\mathbb{R}}^{n\times d}}})$ is a closed subset of
$\mathbb{R}^{m}$ since ${\hat{\beta}}$ is bi-Lipschitz. Then a consequence of
Tietze extension theorem (see problem 8 in §12.1 of [33]) implies that
$\varphi$ admits a continuous extension
$g:\mathbb{R}^{m}\rightarrow\mathbb{R}$. Thus $g(\beta(X))=f(X)$ for all
$X\in\mathbb{R}^{n\times d}$. The converse is trivial.

(2) As at part (1), the Lipschitz continuous function $f$ induces a Lipschitz
continuous function $\varphi:F\rightarrow\mathbb{R}$. Since
$F\subset\mathbb{R}^{m}$ is a subset of a Hilbert space, by Kirszbraun
extension theorem (see [38]), $\varphi$ admits a Lipschitz continuous
extension (even with the same Lipschitz constant!)
$g:\mathbb{R}^{m}\rightarrow\mathbb{R}$ so that $g(\beta(X))=f(X)$ for every
$X\in\mathbb{R}^{n\times d}$. The converse is trivial. $\Box$

##  4\. Applications to Graph Deep Learning

In this section we take an empirical look at the permutation invariant
mappings presented in this paper. We focus on the problems of graph
classification, for which we employ the PROTEINS_FULL dataset [12], and graph
regression, for which we employ the quantum chemistry QM9 dataset [32]. In
both problems we want to estimate a function $F:(A,Z)\rightarrow p$, where
$(A,Z)$ characterizes a graph where $A\in\mathbb{R}^{n\times n}$ is an
adjacency matrix and $Z\in\mathbb{R}^{n\times r}$ is an associated feature
matrix where the $i^{th}$ row encodes an array of $r$ features associated with
the $i^{th}$ node. $p$ is a scalar output where we have $p\in\\{0,1\\}$ for
binary classification and $p\in\mathbb{R}_{+}$ for regression.

We estimate $F$ using a deep network that is trained in a supervised manor.
The network is comprised of three successive components applied in series:
$\Gamma$, $\phi$, and $\zeta$. $\Gamma$ represents a graph deep network [23],
which produces a set of embeddings $X\in\mathbb{R}^{N\times d}$ across the
nodes in the graph. Here $N\geq n$ is chosen to accommodate the graph with the
largest number of nodes. In this case, the last $N-n$ rows of $Y$ are filled
with 0’s. $\phi:\mathbb{R}^{N\times d}\rightarrow\mathbb{R}^{m}$ represents a
permutation invariant mapping such as those proposed in this paper.
$\zeta:\mathbb{R}^{m}\rightarrow\mathbb{R}$ is a fully connected neural
network. The entire end-to-end network is shown in Figure 1.

In this paper, we model $\Gamma$ using a Graph Convolutional Network (GCN)
outlined in [23]. Let ${\bf D}\in\mathbb{R}^{n\times n}$ be the associated
degree matrix for our graph $\mathcal{G}$. Also let $\tilde{A}$ be the
associated adjacency matrix of $\mathcal{G}$ with added self connection:
$\tilde{A}=I+A$, where $I$ is the $n\times n$ identity matrix, and
$\tilde{{\bf D}}={\bf D}+I$. Finally, we define the modified adjacency matrix
$\hat{A}=\tilde{{\bf D}}^{-1/2}\tilde{A}\tilde{{\bf D}}^{-1/2}$. A GCN layer
is defined as $H^{(l+1)}=\sigma(\hat{A}H^{(l-1)}W^{(l)})$. Here $H^{(l-1)}$
represents the GCN state coming into the $l^{th}$ layer, $\sigma$ represents a
chosen nonlinear element-by-element operation such as ReLU, and $W^{(l)}$
represents a matrix of trainable weights assigned to the $l^{th}$ layer whose
number of rows match the number of columns in $H^{l}$ and number of columns is
set to the size of the embeddings at the (l+1)’th layer. The initial state
$H^{(0)}$ of the network is set to the feature set of the nodes of the graph
$H^{(0)}=Z$.

For $\phi$ we employ seven (7) different methods that are described next.

  1. (1)

ordering: For the ordering method, we set $D=d+1$,
$\phi_{ordering}(X)=\beta_{A}(X)=\downarrow(XA)$ with $A=[I~{}1]$ the identity
matrix followed by a column of ones. The ordering and identity-based mappings
have the notable disadvantage of not producing the same output embedding size
for different sized graphs. To accommodate this and have consistently sized
inputs for $\eta$, we choose to zero-pad $\phi(X)$ for these methods to
produce a vector in $\mathbb{R}^{m}$, where $m=ND=N(d+1)$ and $N$ is the size
of the largest graph in the dataset.

  2. (2)

kernels: For the kernels method,

$(\phi_{kernel}(X))_{j}=\sum_{k=1}^{n}K_{G}(x_{k},a_{j})=\sum_{k=1}^{n}exp(-{\|x_{k}-a_{j}\|}^{2}),~{}~{}j\in[m],$

for $X=[x_{1}|\cdots|x_{n}]^{T}$, where kernel vectors
$a_{1},\ldots,a_{m}\in\mathbb{R}^{d}$ are generated randomly, each element of
each vector is drawn from a standard normal distribution. Each resultant
vector is then normalized to produce a kernel vector of magnitude one. When
inputting the embedding $X$ to the kernels mapping, we first normalized the
embedding for each respective node.

  3. (3)

identity: In this case $\phi_{id}(X)=X$, which is obviously not a permutation
invariant map.

  4. (4)

data augmentation: In this case $\phi_{data\;augment}(X)=X$ but data
augmentation is used. Our data augmentation scheme works as follows. We take
the training set and create multiple permutations of the adjacency and
associated feature matrix for each graph in the training set. We add each
permuted graph to the training set to be included with the original graphs. In
our experiments we use four added permutations for each graph when employing
data augmentation.

  5. (5)

sum pooling: The sum pooling method sums the feature values across the set of
nodes: $\phi_{sum\;pooling}(X)=\mathbf{1}_{n\times 1}^{T}X$.

  6. (6)

sort pooling: The sort pooling method flips entire rows of $X$ so that the
last column is ordered descendingly, $\phi_{sort\;pool}(X)=\Pi X$ where
$\Pi\in{\mathcal{S}}_{n}$ so that $\Pi\,X(:,d)=\downarrow(X(:,d))$.

  7. (7)

set-2-set: This method employs a recurrent neural network that achieves
permutation invariance through attention-based weighted summations. It has
been introduced in [36].

For our deep neural network $\eta$ we use a simple multilayer perceptron of
size described below.

Size parameters related to $\Gamma$ and $\zeta$ components are largely held
constant across the different implementations. However the network parameters
are trained independently for each method.

Figure 1. .

###  4.1. Graph Classification

####  4.1.1. Methodology

For our experiments in graph classification we consider the PROTEINS_FULL
dataset obtained from [22] and originally introduced in [12]. The dataset
consists of 1113 proteins falling into one of two classes: those that function
as enzymes and those that do not. Across the dataset there are 450 enzymes in
total. The graph for each protein is constructed such that the nodes represent
amino acids and the edges represent the bonds between them. The number of
amino acids (nodes) vary from around 20 to a maximum of 620 per protein with
an average of 39.06. Each protein comes with a set of features for each node.
The features represent characteristics of the associated amino acid
represented by the node. The number of features is $r=29$. We run the end-to-
end model with three GCN layers in $\Gamma$, each with 50 hidden units.
$\zeta$ consists of three dense multi-layer perceptron layers, each with 150
hidden units. We set d equal to 1, 10, 50 and 100.

For each method and embedding size we train for 300 epochs. Note though that
the data augmentation method will have experienced five times as many training
steps due to the increased size of its training set. We use a batch size of
128 graphs. The loss function minimized during training is the binary cross
entropy loss (BCE) defined as

(4.1)
$BCE=-\frac{1}{B}\sum_{t=1}^{B}p_{t}log(\sigma(\eta(\phi(X^{(t)}))))+(1-p_{t})log(1-\sigma(\eta(\phi(X^{(t)}))))$

where $B=128$ is the batch size, $p_{t}=1$ when the $t^{th}$ graph (protein)
is an enzyme and $p_{t}=0$ otherwise, $\sigma(x)=\frac{1}{1+e^{-x}}$ is the
sigmoid function that maps the output $\eta(\phi(X^{(t)})$ of the 3-layer
fully connected network $\eta$ to $[0,1]$. Three performance metrics were
computed: accuracy (ACC), area under the receiver operating characteristic
curve (AUC), and average precision (AP) as area under the precision-recall
curve from precision scores. These measure are defined as follows (see
sklearn.metrics module documentation in pytorch, or [15]).

For a threshold $\tau\in[0,1]$, the classification decision
$\hat{p}_{t}(\tau)$ is given by:

(4.2)
$\hat{p}_{t}(\tau)=\left\\{\begin{array}[]{rcl}1&\mbox{if}&\mbox{$\sigma(\eta(\phi(X^{(t)}))\geq\tau$}\\\
0&\mbox{if}&\mbox{otherwise}\end{array}\right..$

By default $\tau=\frac{1}{2}$. For a given threshold, one computes the four
scores, true positive (TP), false positive (FP), true negative (TN) and false
negative (FN):

(4.3)
$TP(\tau)=\frac{1}{B_{1}}\sum_{t=1}^{B}1_{\hat{p}_{t}(\tau)=1}1_{p_{t}=1}~{}~{},~{}~{}TN(\tau)=\frac{1}{B_{0}}1_{\hat{p}_{t}(\tau)=0}1_{p_{t}=0}$
(4.4)
$FP(\tau)=\frac{1}{B_{0}}\sum_{t=1}^{B}1_{\hat{p}_{t}(\tau)=1}1_{p_{t}=0}=1-TN(\tau)~{}~{},~{}~{}FN(\tau)=\frac{1}{B_{1}}\sum_{t=1}^{B}1_{\hat{p}_{t}(\tau)=0}1_{p_{t}=1}=1-TP(\tau)$

where $B_{0}=\sum_{t=1}^{B}1_{p_{t}=0}$ and
$B_{1}=\sum_{t=1}^{B}1_{p_{t}=1}=B-B_{0}$.

These four statistics predict Precision $P(\tau)$, Recall $R(\tau)$ (also
known as sensitivity or true positive rate), and Specificity $S(\tau)$ (also
known as true negative rate)

(4.5)
$P(\tau)=\frac{TP(\tau)}{TP(\tau)+FP(\tau)}~{}~{},~{}~{}R(\tau)=\frac{TP(\tau)}{TP(\tau)+FN(\tau)}~{}~{},~{}~{}S(\tau)=\frac{TN(\tau)}{TN(\tau)+FP(\tau)}$

Accuracy (ACC) is defined as the fraction of correct classification for
default threshold $\tau=\frac{1}{2}$ over the set of batch samples:

(4.6)
$ACC=\frac{1}{B}\sum_{t=1}^{B}1_{p_{t}=\hat{p}_{t}(\frac{1}{2})}=\frac{B_{0}}{B}TN(\frac{1}{2})+\frac{B_{1}}{B}TP(\frac{1}{2})$

Area under the receiver operating characteristic curve (AUC) is computed from
prediction scores as the area under true positive rate (TPR) vs. false
positive rate (FPR) curve, i.e. the recall vs. 1-specificity curve

(4.7)
$AUC=\frac{1}{2}\sum_{k=1}^{K}(S(\tau_{k-1})-S(\tau_{k}))(R(\tau_{k-1})+R(\tau_{k}))$

where $K$ is the number of thresholds. Average precision (AP) summarizes a
precision-recall curve as the weighted mean of precision achieved at each
threshold, with the increase in recall from the previous thresholds used as
the weight:

(4.8) $AP=\sum_{k=1}^{K}(R(\tau_{k})-R(\tau_{k-1}))P(\tau_{k}).$

We track the binary cross entropy (BCE) through training and we compute it on
the holdout set and a random node permutation of the holdout set (see Figures
2 and 3). The lower the value the better.

We look at the three performance metrics on the training set, the holdout set,
and a random node permutation of the holdout set: see Figures 4, and 5 for
accuracy (ACC); see Figures 6, and 7 for area under the receiver operating
characteristic curve (AUC); and see Figures 8, and 9 for average precision
(AP). For all these performance metrics, the higher the score the better.

####  4.1.2. Discussion

Tables 1-12 list values of the three performance metrics (ACC, AUC, AP) at the
end of training (after 300 epochs). Performances over the course of training
are plotted in Figures 2 through 9.

The authors of [22] utilized a Support Vector Machine (1-layer perceptron) for
classification and obtained an accuracy (ACC) of 77% on the entire data set
using 52 features, and an accuracy of 80% on a smaller set of 36 features. By
comparison, our data augmentation method for $d=100$ achieved an accuracy of
97.5% on training data set, but dropped dramatically to 73% on holdout data,
and 72% on holdout data set with randomly permuted nodes. On the other hand,
both the kernels method and the sum-pooling method with $d=50$ achieved an
accuracy of around 79% on training data set, while dropping accuracy
performance by only 2% to around 77% on holdout data (as well as holdout data
with nodes permuted).

For $d=1$, data augmentation performed the best on the training set with an
area under the receiver operating characteristic (AUC) of 0.896, followed
closely by the identity method with an AUC of 0.886. On the permuted holdout
set however, sort-pooling performed the best with an AUC of 0.803.

For $d=10$, sum-pooling, ordering, and kernels performed well on the permuted
holdout set with AUC’s of 0.821, 0.820, and 0.818 respectively. The high
performance of the identity method, data augmentation, and sort-pooling on the
training set did not translate to the permuted holdout set at $d=10$. By
$d=100$, sum-pooling still performed the best on the permuted holdout set with
an AUC of 0.817. This was followed by the kernels method which achieved an AUC
of 0.801 on the permuted holdout set.

For experiments where $d>1$, the identity method and data augmentation show a
notable drop in performance from the training set to the holdout set. This
trend is also, to a lesser extent, visible in the sort pooling and ordering
methods. In the holdout permuted set we see significant oscillations in the
performance of both the identity and data augmentation methods.

###  4.2. Graph Regression

####  4.2.1. Methodology

For our experiments in graph regression we consider the qm9 dataset [32]. This
dataset consists of 134 thousand molecules represented as graphs, where the
nodes represent atoms and edges represent the bonds between them.

Each graph has between 3 and 29 nodes, $3\leq n\leq 29$. Each node has 11
features, $r=11$. We hold out 20 thousand of these molecules for evaluation
purposes. The dataset includes 19 quantitative features for each molecule.

For the purposes of our study, we focus on electron energy gap (units $eV$),
which is $\Delta\varepsilon$ in [14] whose chemical accuracy is $0.043eV$ and
whose prediction performance of any machine learning technique is worse than
any other feature. The best existing estimator for this feature is
enn-s2s-ens5 from [17] and has a mean absolute error (MAE) of $0.0529eV$ which
is $1.23$ larger than the chemical accuracy. We run the end to end model with
three GCN layers in $\Gamma$, each with 50 hidden units. $\eta$ consists of
three multi-layer perceptron layers, each with 150 hidden units. We use
rectified linear units as our nonlinear activation function. Finally, we vary
$d$, the size of the node embeddings that are outputted by $\Gamma$. We set
$d$ equal to 1, 10, 50 and 100.

For each method and embedding size we train for 300 epochs. Note though that
the data augmentation method will have experienced five times as many training
steps due to the increased size of its training set. We use a batch size of
128 graphs. The loss function minimized during training is the mean square
error (MSE) between the ground truth and the network output (see Figures 10,
11)

(4.9)
$MSE=\frac{1}{B}\sum_{t=1}^{B}|\Delta\varepsilon_{t}-\eta(\phi(X^{(t)}))))|^{2}$

where $B=128$ is the batch size of 128 graphs and $\Delta\varepsilon_{t}$ is
the electron energy gap of the $t^{th}$ graph (molecule). The performance
metric is Mean Absolute Error (MAE)

(4.10)
$MAE=\frac{1}{B}\sum_{t=1}^{B}|\Delta\varepsilon_{t}-\eta(\phi(X^{(t)}))))|.$

We track the mean absolute error through the course of training. We look at
this performance metric on the training set, the holdout set, and a random
node permutation of the holdout set (see Figures 12, and 13).

####  4.2.2. Discussion

Numerical results at the end of training (after 300 epochs) are included in
Tables 13, 14, 15 and 16. From the results we see that the ordering method
performed best for $d=100$ followed closely by the data augmentation method,
while both the ordering method and the kernels method performed well for
$d=10$, though both fell slightly short of data augmentation which performed
marginally better on both the training data and the holdout data, though with
significantly more training iterations. For $d=1$, the kernels method failed
to train adequately. The identity mapping performed relatively well on
training data (for $d=100$ it achieved the smallest MAE among all methods and
all parameters) and even the holdout data, however it lost its performance on
the permuted holdout data. The identity mapping’s failure to generalize across
permutations of the holdout set is likely exacerbated by the fact that the QM9
data as presented to the network comes ordered in its node positions from
heaviest atom to lightest. Data augmentation notably kept its performance
despite this due to training on many permutations of the data.

For $d=100$, our ordering method achieved a MAE of $0.155eV$ on training data
set and $0.187eV$ on holdout data set, which are $3.6$ and $4.35$ times larger
than the chemical accuracy ($0.043eV$), respectively. This is worse than the
enn-s2s-ens5 method in [17] (current best method) that achieved a MAE $0.0529$
(eV), $1.23$ larger than the chemical accuracy, but better than the Coulomb
Matrix (CM) representation in [34] that achieved a MAE $5.32$ larger than the
chemical accuracy whose features were optimized for this task.

## References

  * [1] B. Alexeev, J. Cahill, and Dustin G. Mixon, _Full spark frames_ , J. Fourier Anal. Appl 18 (2012), 1167–1194.
  * [2] A.S. Bandeira, J. Cahill, D. Mixon, A.A. Nelson, _Saving phase: Injectivity and stability for phase retrieval_ , Appl. Comp. Harm. Anal. 37 (2014), no. 1, 106–125.
  * [3] A. Auslender and Teboulle.M., _Asymptotic cones and functions in optimization and variational inequalities_ , Springer, 2003.
  * [4] R. Balan, _Frames and phaseless reconstruction_ , vol. Finite Frame Theory: A Complete Introduction to Overcompleteness, Proceedings of Symposia in Applied Mathematics, no. 73, pp. 175–199, AMS Short Course at the Joint Mathematics Meetings, San Antonio, January 2015 (Ed. K.Okoudjou), 2016.
  * [5] R. Balan and Y. Wang, _Invertibility and robustness of phaseless reconstruction_ , Applied and Comput. Harmon. Analysis 38 (2015), no. 3, 469–488.
  * [6] R. Balan and D. Zou, _On lipschitz analysis and lipschitz synthesis for the phase retrieval problem_ , Linear Algebra and Applications 496 (2016), 152–181.
  * [7] A. Beck, A. Ben-Tal, and M. Teboulle, _Finding a global optimal solution for a quadratically constrained fractional quadratic problem with applications to the regularized total least squares_ , SIAM J. Mat. Anal. Appl. 28 (2006), no. 2, 425–445.
  * [8] A. Beck and M. Teboulle, _On minimizing quadratically constrained ratio of two quadratic functions_ , J. Convx. Anal. 17 (2010), no. 3,4, 789–804.
  * [9] B.G. Bodmann, N. Hammen, _Stable phase retrieval with low-redundancy frames_ , Adv. Comput. Math. 41 (2015), 317–331.
  * [10] Cédric Villani, _Topics in optimal transportation_ , American Mathematical Society, 2003.
  * [11] Zhengdao Chen, Soledad Villar, Lei Chen, and Joan Bruna, _On the equivalence between graph isomorphism testing and function approximation with gnns_ , Advances in Neural Information Processing Systems (H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, E. Fox, and R. Garnett, eds.), vol. 32, Curran Associates, Inc., 2019.
  * [12] Paul D. Dobson and Andrew J. Doig, _Distinguishing Enzyme Structures from Non-enzymes Without Alignments_ , Journal of Molecular Biology 330 (2003), no. 4, 771–783.
  * [13] Emilie Dufresne, _Separating invariants and finite reflection groups_ , Advances in Mathematics 221 (2009), no. 6, 1979–1989.
  * [14] Felix A. Faber, Bing Hutchison, Luke Huang, Justin Gilmer, Samuel S. Schoenholz, George E. Dahl, Oriol Vinyals, Steven Kearnes, Patrick F. Riley, and O. Anatole von Lilienfeld, _Machine learning prediction errors better than dft accuracy_ , J. Chem. Theory Comput. (2017), no. 13, 5255–5264.
  * [15] Tom Fawcett, _An introduction to roc analysis_ , Pattern Recognition Letters 27 (2006), 861–874.
  * [16] Floris Geerts and Juan L Reutter, _Expressiveness and approximation properties of graph neural networks_ , International Conference on Learning Representations, 2022.
  * [17] Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl, _Neural message passing for quantum chemistry_ , Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, JMLR.org, 2017, p. 1263–1272.
  * [18] G. Kemper H. Derksen, _Computational invariant theory_ , Springer, 2002.
  * [19] Morris Hirsch, _Differential topology_ , Springer, 1994.
  * [20] A.C. Hip J. Cahill, A. Contreras, _Complete set of translation invariant measurements with lipschitz bounds_ , Appl. Comput. Harm. Anal. 49 (2020), no. 2, 521–539.
  * [21] Nicolas Keriven and Gabriel Peyré, _Universal invariant and equivariant graph neural networks_ , Advances in Neural Information Processing Systems 32 (2019), 7092–7101.
  * [22] Kristian Kersting, Nils M. Kriege, Christopher Morris, Petra Mutzel, and Marion Neumann, _Benchmark data sets for graph kernels_ , 2016, http://graphkernels.cs.tu-dortmund.de.
  * [23] Thomas N. Kipf and Max Welling, _Semi-supervised classification with graph convolutional networks_ , International Conference on Learning Representations (ICLR), 2017.
  * [24] Yujia Li, Daniel Tarlow, Marc Brockschmidt, and Richard Zemel, _Gated Graph Sequence Neural Networks_ , arXiv e-prints (2015), arXiv:1511.05493.
  * [25] D.G. Schaeffer M. Golubitsky, I. Stewart, _Singularities and groups in bifurcation theory_ , vol. 2, Springer, 1988.
  * [26] Romanos Diogenes Malikiosis and Vignon Oussa, _Full spark frames in the orbit of a representation_ , Applied and Computational Harmonic Analysis 49 (2020), no. 3, 791–814.
  * [27] Haggai Maron, Heli Ben-Hamu, Nadav Shamir, and Yaron Lipman, _Invariant and equivariant graph networks_ , International Conference on Learning Representations, 2019.
  * [28] Haggai Maron, Ethan Fetaya, Nimrod Segol, and Yaron Lipman, _On the universality of invariant networks_ , Proceedings of the 36th International Conference on Machine Learning (Kamalika Chaudhuri and Ruslan Salakhutdinov, eds.), Proceedings of Machine Learning Research, vol. 97, PMLR, 09–15 Jun 2019, pp. 4363–4371.
  * [29] R.R. Phelps, _Convex sets and nearest points_ , Proc. Amer. Math. Soc. 8 (1957), 790 – 797.
  * [30] Omri Puny, Matan Atzmon, Edward J. Smith, Ishan Mishra, Aditya Grover, Heli Ben-Hamu, and Yaron Lipman, _Frame averaging for invariant and equivariant network design_ , International Conference on Learning Representations, 2022.
  * [31] Charles R Qi, Hao Su, Kaichun Mo, and Leonidas J Guibas, _Pointnet: Deep learning on point sets for 3d classification and segmentation_ , Proceedings of the IEEE conference on computer vision and pattern recognition, 2017, pp. 652–660.
  * [32] Raghunathan Ramakrishnan, Pavlo O Dral, Matthias Rupp, and O Anatole Von Lilienfeld, _Quantum chemistry structures and properties of 134 kilo molecules_ , Scientific data 1 (2014), no. 1, 1–7.
  * [33] H.L. Royden and P.M. Fitzpatrick, _Real analysis_ , 4th ed., Pearson Education, Inc., 2010.
  * [34] Matthias Rupp, Alexandre Tkatchenko, Klaus-Robert Müller, and O. Anatole von Lilienfeld, _Fast and accurate modeling of molecular atomization energies with machine learning_ , Phys. Rev. Lett. 108 (2012), 058301\.
  * [35] Akiyoshi Sannai, Yuuki Takai, and Matthieu Cordonnier, _Universal approximations of permutation invariant/equivariant functions by deep neural networks_ , 2020.
  * [36] Oriol Vinyals, Samy Bengio, and Manjunath Kudlur, _Order Matters: Sequence to sequence for sets_ , International Conference on Learning Representations, 2016.
  * [37] B.Yu Weisfeiler and A.A. Leman, _The reduction of a graph to canonical form and the algebra which appears therein_ , Nauchno-Technicheskaya Informatsia 2 (1968), no. 9, 12–16, English translation by G. Ryabov is available at https://www.iti.zcu. cz/wl2018/pdf/wl_paper_translation.pdf.
  * [38] J.H. Wells and L.R Williams, _Embeddings and extensions in analysis_ , Springer-Verlag, 1975, Ergebnisse der Mathematik und ihrer Grenzgebiete Band 84\.
  * [39] Dmitry Yarotsky, _Universal approximations of invariant maps by neural networks_ , Constructive Approximation (2021), 1–68.
  * [40] Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Russ R Salakhutdinov, and Alexander J Smola, _Deep sets_ , Advances in Neural Information Processing Systems (I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, eds.), vol. 30, Curran Associates, Inc., 2017.
  * [41] Muhan Zhang, Zhicheng Cui, Marion Neumann, and Yixin Chen, _An end-to-end deep learning architecture for graph classification_ , Thirty-Second AAAI Conference on Artificial Intelligence, 2018.

##  Appendix A Results for the PROTEINS_FULL dataset

d = 1 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  76 |  72 |  80 |  81.6 |  76.2 |  78 |  72.4
Holdout |  74 |  74 |  72.5 |  76.5 |  70.5 |  74.5 |  72
Holdout Perm |  74 |  74 |  67.5 |  75 |  70.5 |  74.5 |  72

Table 1. Accuracy ACC(%) for enzyme/non-enzyme classification of the seven
algorithms on PROTEINS_FULL dataset after 300 epochs for embedding dimension
$d=1$

d = 10 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  84.5 |  78.2 |  87 |  90.6 |  77.8 |  85.2 |  72.5
Holdout |  74 |  75.5 |  73 |  76 |  75 |  71 |  74.5
Holdout Perm |  74 |  75.5 |  62.5 |  73.5 |  75 |  71 |  74.5

Table 2. Accuracy ACC(%) for enzyme/non-enzyme classification of the seven
algorithms on PROTEINS_FULL dataset after 300 epochs for embedding dimension
$d=10$

d = 50 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  83.1 |  78.8 |  91 |  96 |  79.2 |  83.7 |  76.7
Holdout |  71.5 |  76.5 |  72.5 |  71 |  77 |  71 |  76
Holdout Perm |  71.5 |  76.5 |  69.5 |  72 |  77 |  71 |  76

Table 3. Accuracy ACC(%) for enzyme/non-enzyme classification of the seven
algorithms on PROTEINS_FULL dataset after 300 epochs for embedding dimension
$d=50$

d = 100 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  88 |  77 |  97.5 |  97.5 |  78.1 |  87.3 |  76.6
Holdout |  71 |  74.5 |  72.5 |  73 |  75.5 |  69.5 |  74.5
Holdout Perm |  71 |  74.5 |  68.5 |  72 |  75.5 |  69.5 |  74.5

Table 4. Accuracy ACC(%) for enzyme/non-enzyme classification of the seven
algorithms on PROTEINS_FULL dataset after 300 epochs for embedding dimension
$d=100$

d = 1 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  0.846 |  0.758 |  0.886 |  0.896 |  0.818 |  0.858 |  0.778
Holdout |  0.794 |  0.775 |  0.766 |  0.796 |  0.777 |  0.803 |  0.788
Holdout Perm |  0.794 |  0.775 |  0.747 |  0.785 |  0.777 |  0.803 |  0.788

Table 5. Area under the receiver operating characteristic curve (AUC) for
enzyme/non-enzyme classification of the seven algorithms on PROTEINS_FULL
dataset after 300 epochs for embedding dimension $d=1$

d = 10 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  0.913 |  0.849 |  0.941 |  0.970 |  0.842 |  0.930 |  0.787
Holdout |  0.820 |  0.817 |  0.782 |  0.796 |  0.821 |  0.798 |  0.779
Holdout Perm |  0.820 |  0.817 |  0.668 |  0.784 |  0.821 |  0.798 |  0.779

Table 6. Area under the receiver operating characteristic curve (AUC) for
enzyme/non-enzyme classification of the seven algorithms on PROTEINS_FULL
dataset after 300 epochs for embedding dimension $d=10$

d = 50 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  0.922 |  0.847 |  0.965 |  0.994 |  0.856 |  0.920 |  0.820
Holdout |  0.791 |  0.818 |  0.775 |  0.768 |  0.821 |  0.791 |  0.777
Holdout Perm |  0.791 |  0.818 |  0.716 |  0.768 |  0.821 |  0.791 |  0.777

Table 7. Area under the receiver operating characteristic curve (AUC) for
enzyme/non-enzyme classification of the seven algorithms on PROTEINS_FULL
dataset after 300 epochs for embedding dimension $d=50$

d = 100 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  0.949 |  0.832 |  0.997 |  0.997 |  0.849 |  0.948 |  0.842
Holdout |  0.754 |  0.801 |  0.766 |  0.775 |  0.817 |  0.784 |  0.776
Holdout Perm |  0.754 |  0.801 |  0.708 |  0.775 |  0.817 |  0.784 |  0.776

Table 8. Area under the receiver operating characteristic curve (AUC) for
enzyme/non-enzyme classification of the seven algorithms on PROTEINS_FULL
dataset after 300 epochs for embedding dimension $d=100$

d = 1 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  0.788 |  0.709 |  0.844 |  0.857 |  0.754 |  0.811 |  0.707
Holdout |  0.720 |  0.698 |  0.692 |  0.725 |  0.636 |  0.680 |  0.708
Holdout Perm |  0.720 |  0.698 |  0.622 |  0.710 |  0.636 |  0.680 |  0.708

Table 9. Average precision (AP) for enzyme/non-enzyme classification of the
seven algorithms on PROTEINS_FULL dataset after 300 epochs for embedding
dimension $d=1$

d = 10 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  0.890 |  0.804 |  0.922 |  0.961 |  0.797 |  0.904 |  0.722
Holdout |  0.738 |  0.749 |  0.631 |  0.646 |  0.753 |  0.693 |  0.693
Holdout Perm |  0.738 |  0.749 |  0.497 |  0.664 |  0.753 |  0.693 |  0.693

Table 10. Average precision (AP) for enzyme/non-enzyme classification of the
seven algorithms on PROTEINS_FULL dataset after 300 epochs for embedding
dimension $d=10$

d = 50 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  0.899 |  0.797 |  0.950 |  0.991 |  0.814 |  0.891 |  0.757
Holdout |  0.700 |  0.738 |  0.627 |  0.589 |  0.750 |  0.676 |  0.666
Holdout Perm |  0.700 |  0.738 |  0.520 |  0.600 |  0.750 |  0.676 |  0.666

Table 11. Average precision (AP) for enzyme/non-enzyme classification of the
seven algorithms on PROTEINS_FULL dataset after 300 epochs for embedding
dimension $d=50$

d = 100 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  0.933 |  0.777 |  0.995 |  0.995 |  0.806 |  0.927 |  0.782
Holdout |  0.627 |  0.729 |  0.601 |  0.622 |  0.747 |  0.637 |  0.704
Holdout Perm |  0.627 |  0.729 |  0.529 |  0.656 |  0.747 |  0.637 |  0.704

Table 12. Average precision (AP) for enzyme/non-enzyme classification of the
seven algorithms on PROTEINS_FULL dataset after 300 epochs for embedding
dimension $d=100$

Figure 2. Binary Cross Entropy (BCE) for enzyme/non-enzyme classification on
PROTEINS_FULL dataset using; $d=1$ left column; $d=10$ right column.

Figure 3. Binary Cross Entropy (BCE) for enzyme/non-enzyme classification on
PROTEINS_FULL dataset using; $d=50$ left column; $d=100$ right column.

Figure 4. Accuracy ACC(%) for enzyme/non-enzyme classification on
PROTEINS_FULL dataset using; $d=1$ left column; $d=10$ right column.

Figure 5. Accuracy ACC(%) for enzyme/non-enzyme classification on
PROTEINS_FULL dataset using; $d=50$ left column; $d=100$ right column.

Figure 6. Area under the ROC curve (AUC) for enzyme/non-enzyme classification
on PROTEINS_FULL dataset using; $d=1$ left column; $d=10$ right column.

Figure 7. Area under ROC curve (AUC) for enzyme/non-enzyme classification on
PROTEINS_FULL dataset using; $d=50$ left column; $d=100$ right column.

Figure 8. Average precision (AP) for enzyme/non-enzyme classification on
PROTEINS_FULL dataset using; $d=1$ left column; $d=10$ right column.

Figure 9. Average precision (AP) for enzyme/non-enzyme classification on
PROTEINS_FULL dataset using; $d=50$ left column; $d=100$ right column.

##  Appendix B Results for the QM9 dataset

d = 1 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  0.302 |  0.867 |  0.320 |  0.281 |  0.349 |  0.309 |  0.389
Holdout |  0.304 |  0.868 |  0.331 |  0.285 |  0.344 |  0.313 |  0.385
Holdout Perm |  0.304 |  0.868 |  2.433 |  0.298 |  0.344 |  0.313 |  0.385

Table 13. Mean Absolute Error (MAE) for regression of the electron energy gap
$\Delta\varepsilon=LUMO-HOMO$ (eV) of the seven algorithms on QM9 dataset
after 300 epochs for embedding dimension $d=1$

d = 10 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  0.220 |  0.219 |  0.182 |  0.175 |  0.214 |  0.226 |  0.282
Holdout |  0.232 |  0.222 |  0.244 |  0.208 |  0.223 |  0.278 |  0.287
Holdout Perm |  0.232 |  0.222 |  1.099 |  0.216 |  0.223 |  0.278 |  0.287

Table 14. Mean Absolute Error (MAE) for regression of the electron energy gap
$\Delta\varepsilon=LUMO-HOMO$ (eV) of the seven algorithms on QM9 dataset
after 300 epochs for embedding dimension $d=10$

d = 50 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  0.163 |  0.257 |  0.163 |  0.172 |  0.182 |  0.166 |  0.196
Holdout |  0.191 |  0.258 |  0.234 |  0.212 |  0.204 |  0.227 |  0.211
Holdout Perm |  0.191 |  0.258 |  1.607 |  0.219 |  0.204 |  0.277 |  0.211

Table 15. Mean Absolute Error (MAE) for regression of the electron energy gap
$\Delta\varepsilon=LUMO-HOMO$ (eV) of the seven algorithms on QM9 dataset
after 300 epochs for embedding dimension $d=50$

d = 100 |  ordering |  kernels |  identity |  data augment |  sum-pooling |  sort-pooling |  set-2-set
---|---|---|---|---|---|---|---
Training |  0.155 |  0.269 |  0.139 |  0.164 |  0.178 |  0.199 |  0.173
Holdout |  0.187 |  0.267 |  0.227 |  0.206 |  0.201 |  0.239 |  0.201
Holdout Perm |  0.187 |  0.267 |  1.086 |  0.213 |  0.201 |  0.239 |  0.201

Table 16. Mean Absolute Error (MAE) for regression of the electron energy gap
$\Delta\varepsilon=LUMO-HOMO$ (eV) of the seven algorithms on QM9 dataset
after 300 epochs for embedding dimension $d=100$

Figure 10. Mean-Square Error (MSE) Results for regression of the electron gap
energy $\Delta\varepsilon=LUMO-HOMO$ (eV) on QM9 dataset using the seven
algorithms; $d=1$ left column; $d=10$ right column.

Figure 11. Mean-Square Error (MSE) Results for regression of the electron gap
energy $\Delta\varepsilon=LUMO-HOMO$ (eV) on QM9 dataset using the seven
algorithms; $d=50$ left column; $d=100$ right column.

Figure 12. Mean-Absolute Error (MAE) Results for regression of the electron
gap energy $\Delta\varepsilon=LUMO-HOMO$ (eV) on QM9 dataset using the seven
algorithms; $d=1$ left column; $d=10$ right column.

Figure 13. Mean-Absolute Error (MAE) Results for regression of the electron
gap energy $\Delta\varepsilon=LUMO-HOMO$ (eV) on QM9 dataset using the seven
algorithms; $d=1$ left column; $d=10$ right column.