# Mismatching as a tool to enhance algorithmic performances of Monte Carlo
methods for the planted clique model

Maria Chiara Angelini1,2, Paolo Fachin1, Simone de Feo1

(1 Dipartimento di Fisica, “Sapienza, Università di Roma”, P.le A. Moro 5,
00185, Rome, Italy
2 INFN, Sezione di Roma1, P.le A. Moro 5, 00185, Rome, Italy)

###### Abstract

Over-parametrization was a crucial ingredient for recent developments in
inference and machine-learning fields. However a good theory explaining this
success is still lacking. In this paper we study a very simple case of
mismatched over-parametrized algorithm applied to one of the most studied
inference problem: the planted clique problem. We analyze a Monte Carlo (MC)
algorithm in the same class of the famous Jerrum algorithm. We show how this
MC algorithm is in general suboptimal for the recovery of the planted clique.
We show however how to enhance its performances by adding a (mismatched)
parameter: the temperature; we numerically find that this over-parametrized
version of the algorithm can reach the supposed algorithmic threshold for the
planted clique problem.

The deep learning revolution showed how over-parametrization can be a crucial
ingredient to improve algorithmic performances in machine-learning, while a
complete theory explaining such a success is still lacking [1, 2]. Over-
parametrized models are mismatched models, in the sense that they can express
a different class of functions than the generative model that produced the
data. While the analysis of general cases can be very involved, one could
obtain very helpful insights from the study of simpler models. In this paper
we study a very simple case of mismatched over-parametrized algorithm: We
consider a standard inference problem for which we know the generative model.
In this case, as usually done in Bayesian inference, we can write the
posterior probability distribution of the variables by using Bayes formula.
This posterior probability can be interpreted as a Gibbs-Boltzmann (GB)
measure over an ad-hoc Hamiltonian [3]. This is the Bayes optimal setting,
where many results can be obtained using statistical physics tools [4]. A very
simple way to introduce overparametrization is by adding an additional
parameter, the temperature $T$. To recover the original Bayes posterior one
could just take the value $T=1$, while studying the GB measure at temperatures
different from 1 results in a mismatched model. The introduction of the
temperature modifies the free-energy landscape associated to the algorithms
possibly changing their performances. Algorithms runned at $T\neq 1$ are in
principle sub-optimal if one were able to sample exactly the GB measure.
However sampling algorithms, taking place in a very high-dimensional space and
typically running for a time much shorter than the equilibration time, rarely
allow equilibrium sampling and their performances are not just related to the
equilibrium free-energy. In these situations it is not obvious that their non-
equilibrium sampling is optimal at $T=1$. How do the performances of different
algorithms change with the temperature $T$? If Bayes optimal algorithms behave
worse when $T$ is not equal to 1, we will show that other algorithms, such as
Monte Carlo ones, become more efficient in the mismatched setting, as already
suggested in ref. [5]. We will show that entropic effects will help in finding
the correct solution, showing that the only knowledge of the equilibrium free-
energy is not enough to determine the performances of out-of-equilibrium
algorithms. In a statistical physics language, it is now well-known in
optimization and inference problems that often the choice to work at $T>0$,
that is to minimize the free-energy, leads to better results that working
exactly at $T=0$, trying to directly minimize the energy. In this paper we
show that finding the proper choice of the best $T$ is far from trivial.

In the following, we focus on a well known inference problem: the planted
clique one. The planted clique problem is the following [6]: Consider an
Erdös–Renyi random graph of size $N$, in which each edge is present
independently with probability $1/2$. In such a random graph one plants a
clique (a completely connected subgraph) of size $K$ manually adding any
possible edge between two nodes belonging to the planted clique. Given a
realization of such a graph, we want to identify the planted clique. The
identification of the planted clique is possible with an exhaustive search
algorithm as long as $K>2\log_{2}(N)$, that is the size of the largest clique
present in a random graph [7]. However many known algorithms are proved to
fail in the regime $K/\sqrt{N}\to 0$, including spectral algorithms [8],
Belief-Propagation based algorithms [9], Sum of Squares Hierarchy [10]: the
planted clique model displays one of the most well-studied computational-
statistical gaps in the literature of inference problems.

One could wonder how Monte Carlo (MC) algorithms work for this problem. In
ref. [6], Jerrum introduced the problem, alongside with a MC algorithm to find
the planted clique. The allowed configuration space of Jerrum algorithm is the
one containing perfect cliques of any size. The author proves the failure of
his MC algorithm for $K=o(\sqrt{N})$ but he does not prove its success for
$K\geq\sqrt{N}$. In ref. [5], the authors analyze the performances of a
different type of Monte-Carlo algorithm that spans the space of configurations
at fixed magnetization $m=K$ (more details in the following sections),
including also ”non-clique” configurations. They show that this type of MC is
not able to recover the planted clique below $K\leq N^{2/3}$ in a polynomial
time with $N$. However, they also show that the performance of the MC method
can be salvaged all the way down to $K=\sqrt{N}$, using a mismatched
magnetization $m>K$.

One can also define a third MC method, in the same family of Jerrum one, that
can be constructed directly from the Posterior distribution obtained from the
Bayes formula: in the following we will call it BayesMC. The BayesMC was
introduced in ref. [11], together with an associated Parallel Tempering (PT)
algorithm. The PT version seems to find the planted clique down to
$K=2\log_{2}(N)$ in a polynomial time for relatively small $N$, but a
theoretical analysis of it is still lacking. In this work we will analyze the
performances of the BayesMC: while at $T=1$ it is sub-optimal, increasing $T$
up to an optimal value, it can find the planted solution down to the regime
$K=O(\sqrt{N})$. The results are in agreement with what conjectured in ref.
[5]. Our reformulation of the mismatching in terms of a well defined
statistical mechanics problem leads to an explicit physical intuition behind
the newly introduced mismatched parameter, that in this case is just the
temperature.

##  1 The model and Monte Carlo algorithms

To be concrete, we construct a graph of $N$ nodes with a planted clique
$\mathcal{C}$ of size $K$. A label $v_{i}$ is associated with each node $i$:
$v_{i}=1$ if node $i\in\mathcal{C}$, $v_{i}=0$ if node $i\notin\mathcal{C}$.
We then extract the graph, univocally identified by its adjacency matrix
$\mathbf{A}$. Its element $A_{ij}$ takes the value 1 if an edge is put between
node $i$ and node $j$, and the value 0 otherwise. The event $A_{ij}=1$ happens
with the following probabilities:

$p(A_{ij}=1|\\{{v\\}})=\begin{cases}1&\text{if }v_{i}v_{j}=1\\\
\frac{1}{2}&\text{otherwise}\end{cases}.$ (1)

Given a realization of the graph, we want to estimate the labels $\\{v\\}$. We
will perform this task using a MC algorithm, and we will indicate with
$x_{i}=\\{0,1\\}$ the configuration if the $i$-th node explored by MC.
Following the Bayes formula, the posterior probability for $x_{i}$ given the
graph is

$P(x_{i}|\\{{A\\}})\propto P(A_{ij}|\\{{x\\}})P(x_{i}),$ (2)

where the likelihood $P(A_{ij}|\\{{x\\}})$ is given by eq. (1) and $P(x)$ is
the prior probability. The original problem has a global constraint on the
size of the clique to be recovered (we are treating the case of known $K$):
$P(\\{{x\\}})$ should be zero if $\sum_{i}x_{i}\neq K$. However, following
Jerrum’s idea, we want to construct a Monte-Carlo algorithm, for which the
actual magnetization $m=\sum_{i}x_{i}$ is not fixed. We then replace the
global constraint with a local prior on the single node111We are looking at
the gran-canonical ensemble, in the language of statistical physics:
$P(x)=\left(\frac{K}{N}\right)^{x}\left(1-\frac{K}{N}\right)^{1-x}$. At this
point we introduce an additional parameter, the inverse temperature
$\beta=\frac{1}{T}$, and a new $\beta$-dependent distribution:

$P_{\beta}(\\{x\\}|\\{{A\\}})\propto P^{\beta}(\\{x\\}|\\{{A\\}})$ (3)

that in the case $\beta=1$ reduces to the standard posterior in eq.(2). We
define an Hamiltonian associated to the posterior probability of eq. (3) as:
$P_{\beta}(\\{x\\}|\\{{A\\}})\equiv\frac{1}{\mathcal{N}}e^{-\beta
H(\\{{x\\}})}$, where $\mathcal{N}$ is a normalization factor and the
Hamiltonian $H$ has the form:

$\displaystyle H(\\{{x\\}})=$ $\displaystyle-\sum_{i}\log(P(x_{i}))+$
$\displaystyle-\sum_{ij}\left[\left(1-A_{ij}\right)\log\frac{\left(1-x_{i}x_{j}\right)}{2}+A_{ij}\log\frac{\left(1+x_{i}x_{j}\right)}{2}\right].$
(4)

Let us look carefully to the single terms in the Hamiltonian. The first sum
$\sum_{i}\log(P(x_{i}))$ acts as a local constant field on the single nodes.
Changes in the energy due to this term are of order
$O(\sqrt{N}\log(\sqrt{N}))$ when $K=O(\sqrt{N})$, thus one could naively
conclude that it is subdominant and it can be neglected. However we will show
its effects in the MC acceptance rate, shown in Fig. 1. The second term
$-\sum_{ij}\left[\left(1-A_{ij}\right)\log\frac{\left(1-x_{i}x_{j}\right)}{2}\right]$
implies $H=\infty$ if there exists at least one couple $i,j$ for which
$A_{ij}=0$ and $x_{i}x_{j}=1$, preventing from having configurations without
links between two elements of a clique: it acts as a hard constraint on the
allowed configurations. Finally the last term
$-\sum_{ij}\left[A_{ij}\log\frac{\left(1+x_{i}x_{j}\right)}{2}\right]$ favours
the addition of nodes to the clique whenever it is possible and it can be
rewritten as $-\sum_{(i,j)\in E}\log(2)\left(x_{i}x_{j}-1\right)$, where $E$
is the set of edges of the graph. Changes in the energy due to this term are
extensive in the regime $K=O(\sqrt{N})$. Summarizing, the Hamiltonian
associated to the planted clique can be recast in the Hamiltonian of an Ising
model with an external field plus an hard constraint on the allowed
configurations.

One can thus define a standard Metropolis algorithm from the Hamiltonian in
eq. (4), that we will call BayesMC, where each proposed move is accepted with
probability $\text{min}\left(1,e^{-\beta\Delta E}\right)$, with $\Delta
E=H(\\{{x^{n+1}\\}})-H(\\{{x^{n}\\}})$ [11]. At time $n$ a spin $i$ is chosen
at random among the $N$ ones. If its value is $x_{i}^{n}=0$, it is flipped
with probability:

$P(x_{i}^{n}=0\to x_{i}^{n+1}=1)=\begin{cases}&0\quad\text{ if }\exists
j:x_{j}^{n}=1\text{ and }A_{ij}=0\\\
&\text{min}\left(e^{-\beta\left[\log(1-\frac{K}{N})-\log(\frac{K}{N})+m\log(\frac{1}{2})\right]},1\right)\text{
o.w.}\end{cases},$ (5)

$m=\sum_{j\neq i=1}^{N}x_{j}^{n}$ being the actual magnetization. On the other
hand, if at time $n$ the value of the randomly chosen spin is $x_{i}^{n}=1$,
it is flipped with probability:

$P(x_{i}^{n}=1\to
x_{i}^{n+1}=0)=\text{min}\left(e^{\beta\left[\log(1-\frac{K}{N})-\log(\frac{K}{N})+m\log(\frac{1}{2})\right]},1\right)$
(6)

A single Monte-Carlo step (MCS) is the attempt to change $N$ spins. The
starting condition is the configuration $x_{i}^{0}=0\quad\forall i$. Being the
system fully connected, the computation of the new energy after a flip of a
spin is of order $O(N)$. However, the proposal to flip a spin is accepted only
$O(K)$ times: A single step of the algorithm thus takes $O(K\cdot N)$ time.
The algorithm stops when it finds a clique of size $K$: for the analyzed
values of $K$, there are no random cliques of that size in the graph, thus one
can be sure that when the algorithm stops it just found the planted clique.

One could compare the BayesMC algorithm with the one introduced by Jerrum in
ref. [6], for which:

$P_{\text{Jerrum}}(x_{i}^{n}=0\to x_{i}^{n+1}=1)=\begin{cases}0\quad&\text{ if
}\exists j:x_{j}^{n}=1\text{ and }A_{ij}=0\\\ 1&\text{ otherwise}\end{cases}$
(7) $P_{\text{Jerrum}}(x_{i}^{n}=1\to x_{i}^{n+1}=0)=\lambda^{-1}$

Jerrum suggests the choice $\lambda=N$ for which cliques of size greater than
$(2-\epsilon)\log_{2}(N)$ dominate when the process is in equilibrium: with
this choice the equilibrium measure should be peaked on the planted clique.
While in the Jerrum MC the transition probabilities are fixed by the parameter
$\lambda$, in the BayesMC the transition probabilities vary during the
simulation depending on the actual magnetization $m$: a comparison between the
acceptance rates for the two algorithms is shown in Fig. 1. The more striking
difference is the high acceptance rate for BayesMC w.r.t the Jerrum MC for the
event $(x_{i}^{n}=1\to x_{i}^{n+1}=0)$ at small $m$, that allows to move fast
between cliques of small sizes. 222One could also choose not to follow detail
balance, and to take a mix between Jerrum and BayesMC for the acceptance
rates: just as an example, a particularly good choose could be Jerrum rate for
the event $(x_{i}^{n}=0\to x_{i}^{n+1}=1)$ and BayesMC for the event
$(x_{i}^{n}=1\to x_{i}^{n+1}=0)$. A similar way of reasoning is the one used
to design heuristic algorithms widely used in satisfiability problems, like
WALKSAT or ASAT [12, 13].

Figure 1:  Acceptance rate for the event $(x_{i}^{n}=0\to x_{i}^{n+1}=1)$
(left) and the event $(x_{i}^{n}=1\to x_{i}^{n+1}=0)$ (right) as a function of
the actual magnetization $m$ at time $n$, for the BayesMC, defined in eqs.
(5,6) at different values of $\beta$ and $K$, compared with the Jerrum
algorithm defined in eq. (7), for a graph with $N=2000$.

We will present a numerical study of the Jerrum algorithm in the appendix,
while in the following sections we will focus our attention on the BayesMC,
for which the physical meaning of the parameters is more explicit.

Let us just mention that one could also define a third MC method, that works
at fixed magnetization $m\equiv K$. Defining an energy
$E=\sum_{ij}(1-A_{ij})x_{i}x_{j}$ that associate a unitary cost to couples of
unconnected nodes $i,j$ inserted in the putative planted clique, one can
analyze the Metropolis algorithm associated to it at a given temperature (in
this setting one should find the ground state to discover the planted clique).
Working at fixed magnetization, the possible moves at time $n$ are the ones of
the type $(x_{a}^{n},x_{b}^{n})=(1,0)\to(x_{a}^{n+1},x_{b}^{n+1})=(0,1)$. In
ref. [5], it is proved the impossibility for this class of MC algorithms to
recover the planted clique for $K\leq N^{2/3}$ in a polynomial time with $N$.
However, the authors also show that the performance of this MC algorithm can
be salvaged all the way down to $K=\sqrt{N}$, working with a mismatched fixed
magnetization $m>K$.

##  2 Numerical Thresholds for the Monte Carlo algorithm

As already said, the BayesMC algorithm defined in eqs. (5,6) that aims to
sample the posterior distribution, should work at $\beta=1$. However, we
numerically checked that for $\beta=1$ BayesMC does not succeed to find the
planted clique in a reasonable time ($t<10^{7}$) even for $K>>\sqrt{N}$ at
relatively small $N$ ($N\simeq 10^{3}$). For this reason, we skipped the
analysis in this case and we moved directly to analyze the case of $\beta<1$.
In Fig. 2 we show $t$, that is the average number of MCS needed to find the
planted clique, as a function of $\beta$ for different $K$ at fixed $N=2000$.

Figure 2:  Average number of MCS $t$ needed to find the planted clique with
BayesMC, as a function of $\beta$ for different $K$ at fixed $N=2000$.
Averages are taken over $\sim 100$ samples. For these values of $K$ and $N$,
the planted clique is not recovered in $t\leq 10^{7}$ if $\beta=1$.

Figure 3:  Left: Average magnetization as a function of the number of MCS $t$
of BayesMC for samples with $N=2000$ and $K=50$, at different values of
inverse temperature $\beta$. In this case $\beta_{opt}\simeq 0.3$. Averages
are taken over 200 samples. Please note how for $T=1$, after $t=10^{5}$ steps,
the average magnetization is still far from reaching $m=K$. Right: Average
energy (with changed sign) for the same samples as in the Left Fig.

The minimum of $t$ is located at $\beta_{opt}(K)<1$: entropic effects seem to
help speeding-up the convergence. For a deeper understanding of the function
of $\beta$, in Fig. 3 we plot the average magnetization and energy as a
function of $t$ at fixed value of $N$ and $K$ and for different inverse
temperatures $\beta$. All the curves show the same qualitative behaviour: for
small $t$, there is a first plateau at small magnetization $m_{1}\lesssim
O(\log_{2}(N))$, of the order of the size of the most numerous random cliques.
$m_{1}$ depends on $\beta$ and it is smaller for smaller $\beta$ (energy in
the first plateau is larger for smaller $\beta$). We identify as
$t^{*}(\beta)$ the time needed to escape this first plateau. If $\beta$ is too
large, the time $t^{*}$ grows. In this case it is plausible that random
cliques of larger sizes $m_{1}$ are trapping the dynamics, while working at
smaller $\beta$ allows to explore smaller random cliques without large
barriers between them. For $t>t^{*}(\beta)$, the magnetization grows until a
second plateau at $m_{2}$ that corresponds to the equilibrium magnetization.
For too small $\beta$, $m_{2}<K$, and the planted clique is found only as a
rare fluctuation. Correspondingly, for too small $\beta$, the equilibrium
energy is higher than that found at the optimal $\beta_{opt}(K)$, that
corresponds to the planted ground state of the system. The minimum $t$ needed
to find the planted clique is obtained for a value $\beta_{opt}(K)$ that is a
compromise between these two opposite effects.

In Fig. 4 we plot the value of optimal $t$, extrapolated at $\beta_{opt}(K)$,
as a function of $K$ for graphs of different sizes $N$: times grow lowering
$K$, and their behaviour is well described by a power-law function
$t=\frac{a(N)}{(K-K_{min}(N))^{\nu}}$.

Figure 4:  Left: Average time $t$ needed by BayesMC for the recovery of the
panted clique as a function of $K$ at the optimal temperature $\beta_{opt}(K)$
for graphs of different sizes $N$. The solid line is the best fit with a
power-law function $t=\frac{a}{(K-K_{min}(N))^{\nu}}$. Averages are taken over
$\sim 10^{3}$ samples. Right: Same data as in left. On the $x$ axis, we put
$K/\sqrt{N}$ to compare data of different sizes, while the recovery times are
scaled linearly with the size $N$: using such rescaled variables data at
different sizes $N$ collapse in the whole analyzed region. The orange vertical
line is the threshold $\frac{K}{\sqrt{N}}=\frac{1}{\sqrt{e}}$ for linear-time
local algorithms.

We want to check if the BayesMC working at the optimal temperature is a
polynomial algorithm in the region $K=O(\sqrt{N})$. For this reason in the
right part of Fig. 4 we plot the same data as a function of
$\frac{K}{\sqrt{N}}$. Scaling the recovery times with a factor $N$, data from
different sizes $N$ collapse: the numerical indications suggest that BayesMC
is a polynomial algorithm for $K=O(\sqrt{N})$. It is known that no linear-time
local algorithm can find the planted clique for
$K<K_{BP}\equiv\sqrt{\frac{N}{e}}$ [9, 14]: BayesMC seems to reach this
threshold.

##  3 BP analysis

In this section we want to gain some insights about the possible phase
transitions happening in the landscape identified by the $\beta$-dipendent
distribution in eq. (3). Given a posterior distribution, one can define a
message passing algorithm, called Belief Propagation (BP), that should
converge to the correct marginals for each node variable, if properly
initialized. The BP equations were written and analyzed in ref. [9] in the
case of $\beta=1$. Here we generalize them to $\beta\neq 1$, to extract the
marginal probabilities for each node from eq. (3). We introduce cavity
messages $\psi_{i\rightarrow j}(x_{i})$ that encode the probability that node
$i$ takes value $x_{i}$, conditioned on the absence of edge $(ij)$. The main
assumptions of BP is that the incoming messages are probabilistically
independent. In the present case the property of the graph to be dense makes
single interactions sufficiently weak so that the assumption of independence
of incoming messages is plausible at the leading order in the large size
limit. One can then write iterative equations on the cavity messages, that
read:

$\displaystyle\psi_{i\rightarrow j}(x_{i}=0)=$
$\displaystyle\frac{1}{z_{i\rightarrow
j}}\left(\frac{N-K}{N}\right)^{\beta}\left(\frac{1}{2}\right)^{\beta(N-1)},$
$\displaystyle\psi_{i\rightarrow j}(x_{i}=1)=$
$\displaystyle\frac{1}{z_{i\rightarrow
j}}\left(\frac{K}{N}\right)^{\beta}\left(\frac{1}{2}\right)^{\beta(N-1)}\prod_{k\backslash
j}\left[1+(2^{\beta}A_{ij}-1)\psi_{k\rightarrow i}(1)\right].$

with $z_{i\rightarrow j}$ a normalization factor ensuring that
$\psi_{i\rightarrow j}(0)+\psi_{i\rightarrow j}(1)=1$.

Once the iteration of the cavity messages has reached a fixed point, marginal
probabilities are obtained as:

$\displaystyle\psi_{i}(x_{i}=0)=$
$\displaystyle\frac{1}{z_{i}}\left(\frac{N-K}{N}\right)^{\beta}\left(\frac{1}{2}\right)^{\beta
N},$ $\displaystyle\psi_{i}(x_{i}=1)=$
$\displaystyle\frac{1}{z_{i}}\left(\frac{K}{N}\right)^{\beta}\left(\frac{1}{2}\right)^{\beta
N}\prod_{k}\left[1+(2^{\beta}A_{ij}-1)\psi_{k\rightarrow i}(1)\right].$ (8)

with $z_{i}$ a different normalization factor ensuring that
$\psi_{i}(0)+\psi_{i}(1)=1$. The Bethe free energy associated to the reached
solution is

$f=-\frac{1}{N}\left[\sum_{i}\log(z_{i})-\sum_{ij}\log(z_{ij})\right],$ (9)

with $z_{ij}=\frac{z_{i}}{z_{i\rightarrow j}}$.

At this point one can assign labels to the spins. A first possibility is an
assignment that maximizes the posterior probability [9]:

$\hat{x}^{(1)}_{i}=\operatorname*{argmax}_{x=0,1}\psi_{i}(x)$

However we look also at another assignment $\\{\hat{x}^{(2)}\\}$: supposing to
know the exact $K$, we assign $x_{i}^{(2)}=1$ to the first $K$ nodes ordered
according to their $\psi_{i}(1)$.

Depending on the chosen way to assign the labels, we can define two different
overlaps with the planted solution. If the assignment $\hat{x}^{(1)}$ is
chosen, we define the overlap, following ref [3], as:

$q^{(1)}=\frac{\frac{1}{N}\sum_{i=1}^{N}\delta_{\hat{x}_{i}^{(1)},v_{i}}-q^{(1)}_{ran}}{1-q^{(1)}_{ran}}$
(10)

with $q^{(1)}_{ran}=1-\frac{K}{N}$ being the overlap of an uninformative
assignment $\hat{x}_{i}^{(1)ran}=0\quad\forall i$, obtained when marginals are
just computed according to the prior. On the other hand, if the assignment
$\hat{x}^{(2)}$ is chosen, we define the overlap as:

$q^{(2)}=\frac{\frac{1}{N}\sum_{i=1}^{N}\delta_{\hat{x}_{i}^{(2)},v_{i}}-q^{(2)}_{ran}}{1-q^{(2)}_{ran}}$
(11)

with $q^{(2)}_{ran}=\frac{K^{2}}{N^{2}}+\frac{(N-K)^{2}}{N^{2}}$ being the
overlap of an uninformative assignment for which the $K$ nodes with
$\hat{x}_{i}^{(1)ran}=1$, are chosen at random among the $N$ possible ones.

Let us briefly recap which is the behaviour of the BP algorithm in the
$\beta=1$ case: For $K>K_{BP}=\sqrt{N/e}$, BP always converges to a fixed-
point (FP) that has a high overlap with the planted solution. For
$K_{d}<K<K_{BP}$, BP initialized around the planted solution finds a high-
overlap solution, while it fails to identify the planted solution if randomly
initialized. $K_{d}$ should correspond to $K_{d}=\log_{2}(N)$ in the
$N\to\infty$ limit 333This threshold lacks of an analytical proof. Moreover
$K_{d}$ suffers from huge finite size effects, and it has been numerically
measured to be $K_{d}\simeq 1.3\log_{2}(N)$ for $N\simeq 10^{4}$ in ref.
[11].. One can compute the Bethe free-energy $f$ associated to the two fixed-
points: for $K>K_{s}$ the global minimum of $f$ corresponds to the planted FP
while for $K<K_{s}$ the FP reached from random initialization has lower free-
energy. $K_{s}=2\log_{2}(N)$ in the thermodynamic limit and identifies the
information-theoretical threshold below which it is impossible to recover the
planted solution. For $K<K_{d}$ even BP with planted initialization fails to
find a high-overlap solution: the planted solution becomes locally unstable.
In a statistical mechanics language $K_{d}$ corresponds to the spinodal point
for the existence of the planted state, while $K_{BP}$ is the spinodal point
for the existence of the low-overlap state. In ref. [9], the state evolution
equations are rigorously derived for the evolution of the cavity messages in
the $N\to\infty$ limit: It is proven that
$\Gamma\equiv\log\left(\sqrt{N}\frac{\psi(1)}{\psi(0)}\right)$ at time $t$ is
a Gaussian variable of mean $\mu_{0}(t)$ and variance $\sigma^{2}(t)$ for
nodes that do not belong to the planted clique, and mean $\mu_{1}(t)$ and
variance $\sigma^{2}(t)$ for nodes that do belong to the planted clique. For
$K>K_{BP}$, $\lim_{t\to\infty}\frac{\mu_{1}(t)-\mu_{0}(t)}{\sigma(t)}=\infty$
and the recovery is an easy task. For $K<K_{BP}$ the above limit is finite,
but different from 0 at least in the whole region $K=\kappa\sqrt{N}$, with
$\kappa=O(1)$ 444The region $K=o(\sqrt{N})$ was not studied in ref. [9]: BP
with random initialization does not find a high-overlap solution, however it
reaches a non-trivial inference accuracy; while $q^{(1)}=0$ in this region,
$q^{(2)}>0$: this is the reason why we introduced it. This is a so-called
hybrid-hard phase: a phase of this type has been shown to be present in many
other inference problems [15, 16].

The transitions at $\beta=1$ can be straightforwardly generalized to
$\beta\neq 1$ looking at the fixed points reached by the BP equations in eq.
(8) from random or planted initialization and to the associated free-energies.

First of all we compute the critical temperature $T_{d}(K)$ above which the
planted solution is not locally stable anymore: at a given temperature $T$, we
initialize the cavity messages around the planted solution, taking
$\psi_{i}(1)=1-\delta$ for $i\in\mathcal{C}$ and $\psi_{i}(1)=\delta$ for
$i\notin\mathcal{C}$, with $\delta$ being a random variable uniformly
extracted in $\delta\in[0,10^{-4}]$. We then run the BP iteration following
eq. (8) and check if the reached FP has high overlap with the planted
solution. For $T>T_{d}(K)$ the reached FP has low overlap with the planted
solution.

Then we compute the critical temperature $T_{BP}(K)$ below which the low-
overlap solution is not locally stable anymore: at a given temperature $T$, we
initialize the cavity messages at random, taking $\psi_{i}(1)=\delta$ $\forall
i$ with $\delta$ being a random variable uniformly extracted in
$\delta\in[0,1]$. We then run the BP iteration following eq. (8): For
$T<T_{BP}(K)$ the reached FP always has high overlap with the planted
solution.

Finally, at fixed $T$ we compare the free-energies computed following eq. (9),
associated to the FP reached from planted or random initialization. In this
way we can locate the thermodynamic threshold $T_{s}(K)$ above which the
global minimum of the Bethe free-energy is not associated to the planted
solution anymore.

Figure 5:  Phase diagram for $N=3000$ in the $K-T$ plane. The red line
corresponds to $T_{s}(K)$: for higher $T$ the recovery is impossible.
$T_{d}(K)$ is at higher $T$, outside the range of the plot. The black line
corresponds to $T_{BP}(K)$: for lower $T$ the randomly initialized BP finds a
high-overlap solution. The color scale is associated to the value of the
overlap $q^{(2)}$ between the FP reached from random initialization and the
planted configuration: for $T>T_{BP}(K)$ the high-overlap fixed point is
missed but BP reaches a FP with small $q^{(2)}\neq 0$. The orange dots
identify the optimal $T$ for the MC at fixed $K$ that always stays in the
hybrid-hard region.

In Fig. 5 we show $T_{s}(K)$, $T_{BP}(K)$ and the overlap $q^{(2)}$ for the
solution found from the randomly initialized BP, together with the temperature
$T_{MC}(K)\equiv\beta_{opt}^{-1}(K)$ for which the time needed by MC to find
the planted solution is minimum, for systems of size $N=3000$. We observe that
$T_{MC}(K)$ is always in the hybrid-hard phase and does not corresponds to any
of the transitions extraced by the BP algorithm. This impossibility to locate
$T_{MC}(K)$ just looking at the BP fixed points and free energies can be due
to two separate effects: First of all, in physics language the BP algorithm
assumes Replica Symmetry (RS). Saying it in another way, BP cannot identify
spurious glassy states that can be present in the free-energy landscape. While
Replica Symmetry is garanteed when one considers the posterior probability
obtained having the perfect knowledge of the generative model (in this case
the so-called Nishimori condition holds), and in particular Replica symmetry
holds for $\beta=1$, when $\beta\neq 1$ Replica Symmetry can be broken and
spourious glassy states can appear, that are not detected by RS BP. If this is
the case, one should take Replica symmetry breaking into account.

Secondly, MC is a sampling algorithm that starts out-of-equilibrium. While one
can analyze the stationary, equilibrium distribution, strong out-of-
equilibrium effects can affect its dynamics. The analysis of these effects is
much more complicated and is a key question, also shared with other algorithms
such as (stochastic) Gradient Descent or Langevin dynamics
[mignacco2021stochasticity].

Figure 6:  Time needed by BP algorithm to reach a fixed point starting from a
random configuration for a system with $N=2000$ at different values of $K$ and
inverse temperature $\beta$. The minimum is always at $\beta\simeq 1$ that
corresponds to the Bayes optimal case. A local maximum is present at the
spinodal point $\beta_{BP}(K)<1$ for the existence of a secondary minimum of
the free-energy in addition to the planted one ($\beta_{BP}(50)\simeq 0.86$,
$\beta_{BP}(40)\simeq 0.93$). For $\beta<\beta_{BP}(K)$, the BP algorithm
finds a low-overlap solution, different from the planted one.

One could also wonder how BP works in the mismatched setting $\beta\neq 1$. In
Fig. 6 we show the time $t_{BP}(\beta)$ needed by BP to reach a fixed point as
a function of $\beta$: differently from MC, BP rapidly deteriorates when
$\beta\neq 1$, $t_{BP}(\beta)$ always showing a minimum at $\beta=1$. When one
has the perfect knowledge of the generative model, BP is for sure the
algorithm to prefer because it is able to extract the exact marginal
distribution for each node from the exact Bayes posterior: it is a so called
Bayes optimal algorithm. However, often in real-world problems, one does not
have the perfect knowledge of the model that generated the data. In a
mismatched setting it is possible that different algorithms could outperform
BP. Our work is showing a practical example of a mismatched setting for which
MC can find the optimal solution to the problem in a region of the used
parameters where BP cannot.

##  4 Conclusions

We showed that both Jerrum MC and BayesMC, that is a Metropolis MC based on
the posterior distribution, are suboptimal in finding a planted clique: they
do not succeed in polynomial time even for $K$ much higher than
$K=O(\sqrt{N})$, that corresponds to the conjectured algorithmic threshold.
The BayesMC can be generalized adding an inverse temperature parameter that
takes the value $\beta=1$ in the original problem. We show that letting
$\beta\neq 1$, the generalized BayesMC can reach the performances of Bayes
optimal algorithms such as Belief Propagation, running in a polynomial time
down to $K=\sqrt{\frac{N}{e}}$. An analogous phenomenon was already shown in
ref. [5] for a different class of MC algorithms, for which the assumption of a
mismatched size of the clique $\overline{K}>K$ allowed the MC to reach optimal
recovery. It would be nice if one could apply the same techniques of ref. [5]
to analytically confirm our numerical results, and if one could link the
optimal mismatched $\overline{K}$ for the MC in ref. [5] to the optimal vale
of $\beta$ for the BayesMC.

Our setting of the problem as a standard statistical mechanics problem leads
to a good physical intuition of what is happening: the MC works better for
higher temperature, where the global minimum of the free-energy is still the
planted one, but entropic effects help in moving faster through the phase
space and reach it with a smaller time.

Let us emphasize that the optimal temperature for the MC to reach the planted
solution is on a region where standard Bayes optimal algorithms such as Belief
Propagation are stuck in a low-overlap solution: this is a clear example that
when we do not have the perfect knowledge of the parameters of the model, the
celebrated BP algorithm can work worse than MC.

We do not have a theoretical way to identify the optimal temperature for MC to
work, and we only extracted it numerically. Taking inspiration from recent
literature, the failure of MC for $T\simeq 1$ could be due to the presence of
spurious glassy states that trap the dynamics [17]. Glassy states are not
identified by BP algorithm that is based on a Replica-Symmetric assumption.
This Replica-Symmetry holds at $\beta=1$ (it corresponds to the so called
Nishimori line [4]), but it is broken as long as $\beta\neq 1$. In some recent
papers [18, sarao2020complex] it is shown how the performances of gradient
descent (GD) based algorithms in some planted problems can be affected by
spurious local minima that trap the dynamics. However, at high enough signal
to noise ratio, the spurious minima at higher energy become saddles with a
negative direction towards the planted state. For this reason, GD algorithms
can reach the planted state, surfing on these high energy saddles that lead
them directly on the planted state, even in the presence of stable local
minima uncorrelated with the signal at lower energy. It could be worth looking
whether a similar phenomenon happens in this case, being higher energy states
more probably visited at higher temperatures. In ref. [19], the authors
analyze the topology of the loss function of a neural network, showing as,
although spurious valleys may exist, they are avoided with high probability on
over-parametrized models.

Usually a simple way to identify the dynamical threshold $T_{d}$ for the
existence of glassy states is the following: one plants a glassy state at
$T_{d}$, then run BP initialized inside the state and see if the state remains
stable under iteration or one eventually goes back to the paramagnetic fixed
point. Unfortunately we cannot use this method in the planted clique problem,
because we have seen how the paramagnetic phase is modified by the
introduction of the planted clique, acquiring a small positive overlap with
the planted solution. To identify the temperature for the birth of glassy
states one should indeed run the 1-step Replica Symmetry Breaking equations
[20] in presence of the planted state. We left this as future work.

Finally we would point out that the present analysis on the role of
temperature in the MC behaviour could help to shed light into the behaviour of
more complex MC based algorithms such as Parallel Tempering (PT) and to
explain why PT numerically finds solutions to the planted clique problem down
to the information-theoretical threshold as shown in ref. [11]. We are
currently investigating this point.

We thank David Gamarnik and Federico Ricci-Tersenghi for very interesting
discussions. This research has been supported by the European Research Council
under the European Unions Horizon 2020 research and innovation programme
(grant No. 694925 – Lotglassy, G. Parisi).

## References

  * [1] Ian Goodfellow, Yoshua Bengio, Aaron Courville, and Yoshua Bengio.  Deep learning, volume 1.  MIT press Cambridge, 2016.
  * [2] Giuseppe Carleo, Ignacio Cirac, Kyle Cranmer, Laurent Daudet, Maria Schuld, Naftali Tishby, Leslie Vogt-Maranto, and Lenka Zdeborová.  Machine learning and the physical sciences.  Reviews of Modern Physics, 91(4):045002, 2019.
  * [3] Aurelien Decelle, Florent Krzakala, Cristopher Moore, and Lenka Zdeborová.  Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications.  Physical Review E, 84(6):066106, 2011.
  * [4] Lenka Zdeborová and Florent Krzakala.  Statistical physics of inference: Thresholds and algorithms.  Advances in Physics, 65(5):453–552, 2016.
  * [5] David Gamarnik and Ilias Zadik.  The landscape of the planted clique problem: Dense subgraphs and the overlap gap property.  arXiv preprint arXiv:1904.07174, 2019.
  * [6] Mark Jerrum.  Large cliques elude the metropolis process.  Random Structures & Algorithms, 3(4):347–359, 1992.
  * [7] Geoffrey R Grimmett and Colin JH McDiarmid.  On colouring random graphs.  In Mathematical Proceedings of the Cambridge Philosophical Society, volume 77, pages 313–324. Cambridge University Press, 1975.
  * [8] Noga Alon, Michael Krivelevich, and Benny Sudakov.  Finding a large hidden clique in a random graph.  Random Structures & Algorithms, 13(3-4):457–466, 1998.
  * [9] Yash Deshpande and Andrea Montanari.  Finding hidden cliques of size $\sqrt{N/e}$ in nearly linear time.  Foundations of Computational Mathematics, 15(4):1069–1128, 2015\.
  * [10] Boaz Barak, Samuel Hopkins, Jonathan Kelner, Pravesh K Kothari, Ankur Moitra, and Aaron Potechin.  A nearly tight sum-of-squares lower bound for the planted clique problem.  SIAM Journal on Computing, 48(2):687–735, 2019.
  * [11] Maria Chiara Angelini.  Parallel tempering for the planted clique problem.  Journal of Statistical Mechanics: Theory and Experiment, 2018(7):073404, 2018.
  * [12] John Ardelius and Erik Aurell.  Behavior of heuristics on large and hard satisfiability problems.  Phys. Rev. E, 74:037702, Sep 2006.
  * [13] Sakari Seitz, Mikko Alava, and Pekka Orponen.  Focused local search for random 3-satisfiability.  Journal of Statistical Mechanics: Theory and Experiment, 2005(06):P06006, 2005.
  * [14] Andrea Montanari.  Finding one community in a sparse graph.  Journal of Statistical Physics, 161(2):273–299, 2015.
  * [15] Thibault Lesieur, Florent Krzakala, and Lenka Zdeborová.  Constrained low-rank matrix estimation: Phase transitions, approximate message passing and applications.  Journal of Statistical Mechanics: Theory and Experiment, 2017(7):073403, 2017.
  * [16] Federico Ricci-Tersenghi, Guilhem Semerjian, and Lenka Zdeborová.  Typology of phase transitions in bayesian inference problems.  Physical Review E, 99(4):042109, 2019.
  * [17] Fabrizio Antenucci, Silvio Franz, Pierfrancesco Urbani, and Lenka Zdeborová.  Glassy nature of the hard phase in inference problems.  Physical Review X, 9:011020, 2019.
  * [18] Stefano Sarao Mannelli, Giulio Biroli, Chiara Cammarota, Florent Krzakala, and Lenka Zdeborová.  Who is afraid of big bad minima? analysis of gradient-flow in a spiked matrix-tensor model.  arXiv preprint arXiv:1907.08226, 2019.
  * [19] Luca Venturi, Afonso S Bandeira, and Joan Bruna.  Spurious valleys in two-layer neural network optimization landscapes.  arXiv preprint arXiv:1802.06384, 2018.
  * [20] Marc Mézard and Giorgio Parisi.  The bethe lattice spin glass revisited.  Eur. Phys. J. B, 20(2):217–233, 2001.
  * [21] LA Fernandez, Víctor Martin-Mayor, G Parisi, and B Seoane.  Temperature chaos in 3d ising spin glasses is driven by rare events.  EPL (Europhysics Letters), 103(6):67003, 2013.

##  Appendix A Numerical analysis of the Jerrum algorithm

We present here a numerical analysis of the Jerrum Monte-Carlo algorithm [6]
for the recovery of a planted clique of size $K$ inside a graph of size $N$.
The acceptance probability for a move for Jerrum algorithm is defined in the
main text in eq. (7). The natural quantity to measure to describe the
algorithmic performances would be the average time needed to find the planted
clique as a function of $K$. However we notice that, decreasing $K$, the
average time over a large number of samples is dominated by worst cases: in
Fig. 7 we show the average over the fastest $n\%$ samples, for different
values of $n$: the growing rate, decreasing $K$, strongly depend on $n$. The
same phenomenon is evident when plotting the cumulative probability $P(t>x)$
that the time to solution $t$ for a given sample is greater than a certain
value $x$, as shown in Fig. 8: while the median smoothly increases decreasing
$K$, alongside the distributions develop heavy tails at very large times. When
this happens, the time averaged over all the samples is no longer concentrated
around the typical case but is dominated by rare slow events. This situation
of hardness driven by rare events, is quite common in the study of disordered
systems at finite sizes [21].

Figure 7:  Average time needed by the Jerrum algorithm to reach the planted
solution for a system with $N=3000$ at different values of $K$: times $t_{n}$
are averaged over the fastest $n\%$ of the samples, curves are shown for
different $n$. For $K\lesssim 500$ the average is dominated by slowest samples
and curves at different values of $n$ grow differently decreasing $K$.

Figure 8:  Left: Cumulative probability $P(t>x)$ that the time $t$ needed by
the Jerrum algorithm to reach the planted solution for a given system with
$N=3000$ would be greater than $x$, at different values of $K$: while the
median smoothly increases decreasing $K$, alongside the distributions develop
heavy tails at very large times. When this happens, the averaged time is no
longer concentrated around the typical case. Right: Cumulative probability
$P(t>x)$ for the BayesMC for systems with $N=3000$ at the optimal inverse
temperature $\beta_{opt}(K)$. The heavy tails phenomenon is not present in
this case.

The heavy tails make the quantification of average times an impossible task.
For this reason we choose to look at the average time over the fastest $50\%$
of the samples, $t_{50}$, that is the quantity shown in Fig. 10. With this
choice we disregard the heavy tail phenomenon that could be strongly related
to finite size effects. In practice we are underestimating the minimum clique
size $K_{min}$ that could be found by the Jerrum algorithm. However, all we
want to show is that the scaling of such threshold is higher than
$O(\sqrt{N})$. We are aware that it is very hard to deduce the exact form of
the scaling of the threshold just by numerical experiments, specially for a
problem like the planted clique with very large finite size effects and it is
beyond the scope of this paper. For this reason, we focus on $t_{50}$, that is
not the perfect choice but it leads to a lower bound on the estimation of
$K_{min}$.

Figure 9:  Time to solution $t_{50}$ averaged over the fastest $50\%$ of the
samples as a function of the size of the clique for different values of $N$.
Solid lines are the best fit of the form in eq. (12).

$t_{50}$ diverges lowering $K$. We fit its behaviour with a function of the
type:

$t_{50}(K)=\frac{a}{(K-K_{min})^{\nu}}.$ (12)

with $a$ and $K_{min}$ depending on $N$. In this way we can extrapolate the
threshold $K_{min}(N)$ that we plot as a function of $N$ in Fig. 10. Assuming
a behaviour of the type:

$K_{min}(N)=bN^{\alpha},$ (13)

we obtain an estimate $\alpha=0.91(4)$ that is much higher than the lower
bound computed analytically by Jerrum $\alpha_{L}=0.5$. The Jerrum MC is thus
suboptimal, as already conjectured in ref. [5]. As already said, numerical
estimation of the threshold exponent $\alpha$ is affected by strong finite
size effects. The natural methodological framework for its analytical
computation would be the one already used in ref. [5] for a different class of
MC algorithms.

Figure 10:  Minimum clique size $K_{min}$ that can be found by the Jerrum
algorithm as a function of the graph size $N$. The red solid line is the best
fit of the form in eq. (13) while the blue line is the analytical lower bound
computed by Jerrum [6].


# Contact Information Flow and Design of Compliance

Kevin Haninger, Marcel Radke, Richard Hartisch, Jörg Krüger Authors affiliated
with the Department of Automation at Fraunhofer IPK, Berlin, Germany.
Corresponding email<EMAIL_ADDRESS>This project has
received funding from the European Union’s Horizon 2020 research and
innovation programme under grant agreement No 820689 — SHERLOCK.

###### Abstract

Identifying changes in contact during contact-rich manipulation can detect
task state or errors, enabling improved robustness and autonomy. The ability
to detect contact is affected by the mechatronic design of the robot,
especially its physical compliance. Established methods can design physical
compliance for many aspects of contact performance (e.g. peak contact force,
motion/force control bandwidth), but are based on time-invariant dynamic
models. A change in contact mode is a discrete change in coupled robot-
environment dynamics, not easily considered in existing design methods.

Towards designing robots which can robustly detect changes in contact mode
online, this paper investigates how mechatronic design can improve contact
estimation, with a focus on the impact of the location and degree of
compliance. A design metric of information gain is proposed which measures how
much position/force measurements reduce uncertainty in the contact mode
estimate. This information gain is developed for fully- and partially-observed
systems, as partial observability can arise from joint flexibility in the
robot or environmental inertia. Hardware experiments with various compliant
setups validate that information gain predicts the speed and certainty with
which contact is detected in (i) monitoring of contact-rich assembly and (ii)
collision detection.

##  1 Introduction

Manipulation aims to achieve a relative arrangement of objects, often with a
desired set of contacts between them. In structured environments, this desired
contact can be achieved with a certain robot pose, but as robots are applied
in less-structured environments, the required robot pose may vary. A robot
which can robustly detect contact mode can improve autonomy, for example using
online contact detection to switch the controller [1] or to provide a reward
signal to a learning process [2].

Figure 1: Contact-rich manipulation task, where a discrete collection of
contact dynamics occur, corresponding to different task phases. Detecting
these changes can help error monitoring like in-hand slip, and the ability to
detect these changes is a function of sensor placement and physical dynamics.

We call a contact mode a constraint over generalized robot and environment
position $q$, i.e. $\phi_{i}(q)=0$ for the $i^{\mathrm{th}}$ mode. This is
often considered an infinitely stiff kinematic constraint in robotics, either
eliminated or enforced with impulses in time-stepping approaches [3]. Here, we
consider a compliant contact model where forces are measured over position $q$
[4, 5]. Contact can also be described with an inequality $\phi_{i}(q)\geq 0$
to describe non-penetration of rigid bodies [6, 3], but we separate the
contact conditions into discrete contact modes with equality constraints [5].

Contact mode can be used for task representation and planning. As the contact
depends on surface geometry, contact modelling is more common in locomotion
where flat surfaces simplify the expression of $\phi(q)$ [6]. In manipulation,
when constraints are simple and known in advance, a task can be represented as
a sequence of contact modes [7], which can then be used for motion planning
[8, 9]. Contact geometry can also be demonstrated with special instruments
[10]. Data-based methods have been used, learning contact geometry on
simplified geometries [11], but not yet demonstrated on robotic manipulation.

Offline contact models are not robust to variation in environmental objects
[12], which can be improved with online detection of contact. Estimation of
environment dynamics has been established, using a Kalman filter to estimate
environment stiffness [13] or impedance [5]. In locomotion, foot contact mode
estimation has been proposed in a stochastic framework [1], shown to improve
robustness to unexpected changes in floor height. In manipulation, discrete
contact mode estimation has been proposed through recursive Bayesian
estimation [14, 15, 16] and support vector machines [17].

The ability to sense contact mode is affected by the hardware design. Hardware
design also affects peak collision force [18, 19, 20], and contact stability
[21, 22]. The location, degree and structure of compliance can be designed to
improve energy efficiency [23, 24], reduce the need for information [25], or
reduce the problem dimension [26]. Compliance can have negative effects,
reducing the feasible feedback control gain [27] and the effectiveness of
feedforward control [28], limiting motion control performance.

The contribution of this paper is a method to design compliance to improve
estimation of the contact mode. The estimation of contact mode is established
[1, 14, 15, 16], but the impact of mechatronic design on this estimation is
not yet investigated. Analytical rules for ‘distinguishability’ of contact
modes have been proposed [29], but the approach here considers compliance and
sensor noise. This paper first introduces the compliant contact models, then
the information gain from position and force measurements in the fully- and
partially-observed case. The information gain is written in closed form for a
linear inertial system. Finally, experiments validate that information gain
predicts the certainty and speed of assembly monitoring and collision
detection, where the location and degree of compliance is varied.

#### Notation

Random variables are typeset in boldface ${\bf n}\sim p(n)$. The normal
distribution is denoted $\mathcal{N}(\mu,\Sigma)$, and
$\mathcal{N}(x|\mu,\Sigma)$ denotes an evaluation of the Gaussian probability
distribution function at $x$. Proportionality $\propto$ indicates a
distribution is written without a normalization constant for compactness. The
matrix determinant $\mathrm{det}$, and can be ignored when the argument is
scalar. A sequence of variables is denoted $q_{1:t}=[q_{1},\dots,q_{t}]$.

##  2 Contact and robot models

This section presents the dynamic model for compliant contact and observation
models for position and force.

###  2.1 Compliant robot model

Let the robot and environment be described by generalized position $q\in
Q=\mathbb{R}^{D}$ with velocity $\dot{q}$. Compliance (in the joint, gripper,
or environment) can then be modelled by adding to the standard manipulator
equations of motion as [30]

$\displaystyle M(q)\ddot{q}+\tau(q,\dot{q})-K(q)=B(q)\tau_{m},$ (1)

where $M(q)\in\mathbb{R}^{D\times D}$ is the inertia matrix,
$\tau(q,\dot{q})\in\mathbb{R}^{D}$ contains the gravitational, Coriolis, and
friction forces, spring-like forces
$K:\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}$, and input matrix
$B(q)\in\mathbb{R}^{D\times K}$ modifies motor torque input
$\tau_{m}\in\mathbb{R}^{K}$.

This model will be reduced to a two-mass model in the direction of contact
motion, but provides the general statement for contact state and observations.

###  2.2 Contact State

Contact is modelled as an equality constraint over total configuration,
$\phi(q)=0$ [3]. A compliant contact model is used here [5, 4] instead of
impulse-based [3] – it is assumed that sufficient compliance is present that
collision forces can be represented as a function over configuration $q$.

While inequality constraints $\phi(q)\geq 0$ can capture non-penetration
constraints, thereby modelling both free-space and constrained dynamics in a
unified way [6], the approach here instead uses a collection of equality
constraints, only one of which is active at once. Therefore, the contact state
is then a discrete variable $n\in[1,\dots,N]$, where each contact state $n$
has a corresponding constraint $\phi_{n}(q)=0$.

###  2.3 Observation model

Denote a position measurement at time step $t$ which measures some of the
robot/environment positions $\theta_{t}\in\Theta\subset Q$, representing, e.g.
that motor positions are measured, but the position of links or environmental
inertias are not.

Let the measured torques/forces $f_{t}\in\mathbb{R}^{L}$ be a subset of the
forces on compliant elements,

${\bf f}_{t}=K^{n}_{f}(q_{t})+{\bf\epsilon}_{t}$ (2)

where ${\bf\epsilon}_{t}\sim\mathcal{N}(0,\Sigma_{f})$ is i.i.d. noise
capturing electrical or quantisation noise, and
$K^{n}_{f}:Q\rightarrow\mathbb{R}^{L}$ the effective sensor stiffness for the
mode $n$. Therefore

$p(f_{t}|n,q_{t})=\mathcal{N}(K^{n}_{f}(q_{t}),\Sigma_{f}).$ (3)

As force/torque sensors measure displacement - either strain or bulk
displacement [31] \- practical force measurements can be described with (2).

##  3 Contact Information Gain

Bayesian inference is used to recursively infer contact mode given a prior
belief of contact mode, position and force measurements. The prior
distribution can reflect a priori assumptions about the environment, e.g.,
that a contact occurs at a certain location, or it could be a belief informed
by other sensors.

###  3.1 Fully Observed

When the system is fully observed, i.e. $\Theta=Q$, a prior belief $p({n})$ is
updated by a measurement $f,\,\theta$ as

$\displaystyle
p({n}|f,\theta)=\frac{p(f|n,\theta)p({n})}{\sum_{m=1}^{N}p(f|m,\theta)p(m)},$
(4)

where the denominator is a normalization constant.

In the fully observed case, multiple measurements are independent as
$\epsilon_{t}$ in (2) is i.i.d., and the full posterior is

$\displaystyle p(n|f_{1:t},\theta_{1:t})\propto
p(n)\prod_{i=1}^{t}p(f_{i}|n,\theta_{i}).$ (5)

This can be updated recursively as a hidden Markov model [32].

###  3.2 Partially observed

In the partially observed case, some positions are not directly measured. This
can arise from joint flexibility (motor position is measured, link position is
not), or from inertia in the environment. In this case, there are force
contributions from the unmeasured inertias which will appear with the forces
resulting from the contact mode.

We denote measured position $\theta_{t}$ and unmeasured position
$\overline{q}_{t}$, so $q_{t}=[\theta_{t}^{T},\overline{q}_{t}^{T}]^{T}$,
making the force observation model $p(f_{t}|\overline{q}_{t},\theta_{t},n)$. A
dynamic model such as a discretized version of (1) can be used for inference,
where we assume the dynamics
$p(\overline{q}_{t+1}|n,\overline{q}_{t},\theta_{t})$ are known (these are
derived for linear inertial system in Section 4.2). With these dynamics, a
standard recursive estimate can be written [32, Ch. 2.4]

$\begin{split}p(\overline{q}_{t}|n,\theta_{1:t-1},f_{1:t-1}&)\\\ \propto\int
p(\overline{q}_{1})\prod_{i=1}^{t-1}&p(\overline{q}_{i\texttt{+}1}|n,\overline{q}_{i},\theta_{i})p(f_{i}|n,\overline{q}_{i},\theta_{i})\overline{q}_{1:t\texttt{-}1},\end{split}$
(6)

where $p(\overline{q}_{1})$ is a distribution of initial position
${\bf\overline{q}}_{1}$. To calculate (6), a Kalman filter can be used in the
linear case [32, Ch. 3.2], and in the nonlinear flexible-joint robot case,
various observers of link-side position/velocity are established [33].

The observation model (3) then needs to be marginalized over
${\bf\overline{q}}_{t}$ as

$\begin{split}p(f_{t}|n,f_{1:t\texttt{-}1},\theta_{1:t})\\\ \texttt{=}\int
p(f_{t}|n,\overline{q}_{t},&\theta_{t})p(\overline{q}_{t}|f_{1:t\texttt{-}1},\theta_{1:t\texttt{-}1},n)d\overline{q}_{t}.\end{split}$
(7)

We assume the force observations (2) are linear, i.e.
$K^{n}_{f}(q_{t})=K^{n}_{f}(q_{t}-q_{0}^{n})$ where
$K^{n}_{f}\in\mathbb{R}^{L\times D}$ and $q_{0}^{n}$ is the rest position of
the spring. Denote the posterior belief of
$p(q_{t}|n,f_{1:t-1},\theta_{1:t-1})=\mathcal{N}(\mu^{n}_{q},\Sigma^{n}_{q})$,
then the force measurement is distributed as

$p(f_{t}|n,f_{1:t\texttt{-}1},\theta_{1:t\texttt{-}1})=\mathcal{N}(K_{f}^{n}(\mu^{n}_{q}\mathtt{-}q^{n}_{0}),K_{f}^{n}\Sigma_{q}^{n}K_{f}^{n,T}\mathtt{+}\Sigma_{f}).$
(8)

Note that the covariance in ${\bf f}_{t}$ now depends on the uncertainty in
the position $\Sigma_{q}^{n}$, the stiffness $K_{f}^{n}$, and force sensor
noise $\Sigma_{f}$.

###  3.3 Information gain

The information gain is proposed as the reduction in uncertainty from $p(n)$
to $p(n|f,\theta)$, i.e. how much a measurement reduces belief entropy. Recall
that entropy measures the uncertainty of a discrete random variable $H({\bf
n})=-\sum_{n=1}^{N}p(n)\ln(p(n))$ and differential entropy over continuous
random variables $h({\bf x})=-\int p(x)\ln p(x)dx$ [34, Ch. 2], we define the
information gain as

$\displaystyle\Delta_{I}=H({\bf n})-H({\bf n}|{\bf f},\theta).$ (9)

A larger $\Delta_{I}$ means that the measurement brings more certainty.

Denoting a belief vector $b\in\mathbb{R}^{N}$, where $p(n=m)=b_{m}$,
$b_{m}>0$, $\|b\|_{1}=1$, and the distribution of $\bf f$ per mode as

$p({f}|n,\theta)=\mathcal{N}({f}|\mu_{n},\Sigma_{n}),$ (10)

the following relations hold:

$\displaystyle H({\bf n}|{\bf f},\theta)=$ $\displaystyle H({\bf n})-h({\bf
f}|\theta)+h({\bf f}|{\bf n},\theta),$ (11) $\displaystyle h({\bf f}|{\bf
n},\theta)=$ $\displaystyle\frac{L}{2}\ln 2\pi
e+\sum_{n=1}^{N}\frac{b_{n}}{2}\ln\det\Sigma_{n},$ (12) $\displaystyle h({\bf
f}|\theta)\geq$
$\displaystyle-\sum_{n=1}^{N}b_{n}\ln\left(\sum_{m=1}^{N}b_{m}\mathcal{N}\left(\mu_{m}|\mu_{n},\Sigma_{n}\mathtt{+}\Sigma_{m}\right)\right),\,\,\,\,\,\,$
(13)

where (11) results from writing $I({\bf f},{\bf n}|\theta)$ two ways, $h({\bf
f}|\theta)-h({\bf f}|{\bf n},\theta)=H({\bf n}|\theta)-H({\bf n}|{\bf
f},\theta)$ [35], (12) is standard multivariate Gaussian entropy with the
expectation taken over $\bf n$ [34], and (13) is directly Theorem 2 of [36].
Then, $\Delta_{I}$ can be lower bounded as
$\Delta_{I}\geq\underline{\Delta}_{I}$ where

$\displaystyle\begin{split}\underline{\Delta}_{I}=&-\sum_{n=1}^{N}b_{n}\ln\sum_{m=1}^{N}b_{m}\mathcal{N}(\mu_{n}|\mu_{m},\Sigma_{n}+\Sigma_{m})\\\
&-\frac{L}{2}\ln 2\pi
e-\sum_{n=1}^{N}\frac{b_{n}}{2}\ln\det\Sigma_{n}.\end{split}$ (14)

To demonstrate, the information gain bound $\underline{\Delta}_{i}$ with two
modes $N=2$, one dimension $L=1$, and a flat prior belief $b_{n}=0.5$ is shown
in Figure 2, where $\mu_{1}=0$, $\Sigma_{1}=1$, and $\mu_{2},\,\Sigma_{2}$ are
varied. The minimum is reached when $\mu_{2}=0$, $\Sigma_{2}=1$, where the two
distributions are identical. As $\Sigma_{2}$ decreases, the information gain
increases, as well as when $|\mu_{2}-\mu_{1}|$ increases. Note that the
information gain saturates - the maximum reduction in uncertainty is the prior
uncertainty $H({\bf n})$.

Interpreting this for the contact mode, when the mean force expected from two
contacts vary widely, a force measurement should provide more information gain
and the mode uncertainty should decrease quickly. If the mean of the force
expected in the two modes is the same, there may be some information gain as
long as the covariance is different. However, if the mean and covariance of
two contact mode forces are the same, there is no reduction in mode
uncertainty.

In the partially-observed case, $\underline{\Delta}_{I}=H({\bf n})-H({\bf
n}|{\bf f_{t}},f_{1:t-1},\theta_{1:t})$, modifying the force distribution
(10).

Figure 2: Information gain for single observation with two modes, where
$\mu_{1}=0$, $\Sigma_{1}=1$. As the differences between the two distributions
increase, the information gain increases, and the mode can be more accurately
inferred.

##  4 Information gain in inertial systems

Figure 3: Linear inertial model, where $q^{1}$ is measured robot or motor
position, and $q^{2}$ possibly unmeasured.

In this section, information gain is considered on a linear, inertial system
as seen in Figure 3. We consider that force is measured as
$f_{t}=K_{1}(q^{2}_{t}-q^{1}_{t})+\epsilon_{t}$, where
$\epsilon_{t}\sim\mathcal{N}(0,\Sigma_{f})$.

Consider two modes, free space where $K_{1}=0$, and contact with $M_{2}$ where
$K_{1}>0$. This gives

$\begin{split}p(f|n=1)=&\,\,\mathcal{N}(0,\Sigma_{f})\\\
p(f|n=2,q^{1},q^{2})=&\,\,\mathcal{N}(K_{1}(q^{1}-q^{2}),\Sigma_{f}).\end{split}$
(15)

###  4.1 Fully observed

In the fully observed case, the positions of $q^{1},q^{2}$ are known and there
are no unmeasured inertias. In other words, the robot is stiff, the
environment has no inertias, and the contact force model (3) can be directly
evaluated. If we consider the information gain at a flat prior $b_{n}=0.5$,
(14) can be used with (15) to find

$\displaystyle\underline{\Delta}_{I}(q^{1},q^{2})=$
$\displaystyle-\ln\left(\mathcal{N}(0|K_{1}(q^{1}-q^{2}),2\Sigma_{f})\right)$
$\displaystyle-\frac{1}{2}\ln 2\pi e-\frac{1}{2}\ln\Sigma_{f}$ (16)
$\displaystyle=$
$\displaystyle\frac{1}{2}K_{1}^{2}\left(q^{1}-q^{2}\right)^{2}(2\Sigma_{f})^{-1}-\frac{1}{2}(1-\ln
2),$ (17)

using
$\mathcal{N}(\mu_{n}|\mu_{m},\Sigma_{n}+\Sigma_{m})=\mathcal{N}(\mu_{m}|\mu_{n},\Sigma_{n}+\Sigma_{m})$.

Here, we note that the information gain increases as $K_{1}$ does,
$|q^{1}-q^{2}|$ grows, or as sensor noise $\Sigma_{f}$ decreases. When we have
a stiff environment, deeper contact, or a less noisy sensor, the ability to
estimate contact mode increases.

###  4.2 Partially observed

Here we consider when $q^{2}$ is not observed, e.g. $q^{1}$ is motor position
and $q^{2}$ is the position of an environment or link inertia. The contact
force model (3) cannot be directly evaluated, first $q^{2}$ must be estimated,
increasing uncertainty. As time is relevant in the partially observed case, we
denote $q^{1}_{t}$ and $q^{2}_{t}$ as the values at time step $t$, and
consider that $q^{1}_{t}$ is measured, along with force as in (15). The first-
order discretized system dynamics can be written as

$\begin{split}\left[\begin{array}[]{c}q^{2}_{t+1}\\\
\dot{q}^{2}_{t+1}\end{array}\right]&\approx\left[\begin{array}[]{cc}1&T_{s}\\\
-T_{s}\frac{K_{2}+K_{1}}{M_{2}}&1-T_{s}\frac{B_{2}}{M_{2}}\end{array}\right]\left[\begin{array}[]{c}q^{2}_{t}\\\
\dot{q}^{2}_{t}\end{array}\right]\\\ &+\left[\begin{array}[]{c}0\\\
T_{s}M_{2}^{-1}K_{1}\end{array}\right]q^{1}_{t}+\left[\begin{array}[]{c}0\\\
M_{2}^{-1}\end{array}\right]{\bf w}_{t}\\\
&=A\left[\begin{array}[]{c}q^{2}_{t+1}\\\
\dot{q}^{2}_{t+1}\end{array}\right]+Bq^{1}_{t}+B_{w}{\bf w_{t}}\end{split}$
(18)

where ${\bf w}_{t}\sim\mathcal{N}(0,\Sigma_{w})$ is i.i.d. process noise which
is assumed to be a force disturbance, and thus multiplied by $M_{2}^{-1}$.

A Kalman estimator can be written, where, for simplicity, we consider the
covariance of $p(q^{2}_{t}|f_{1:t},q^{1}_{1:t})$ as $t\rightarrow\infty$,
which is given by the steady-state covariance of $P$, which satisfies the
discrete algebraic Riccati equation of

$P=APA^{T}-APC^{T}(CPC^{T}+\Sigma_{f})^{-1}CPA^{T}+B_{w}\Sigma_{w}B_{w}^{T},$
(19)

where $A$, $B$ and $B_{w}$ are from (18) and $C=[K_{1},0]$ [32].

While the steady-state covariance $P$ is not necessarily the covariance during
transient phases (i.e. when starting from an uncertain initial condition), it
has a simple closed-form solution which is differentiable [37].

With this assumption,
$p(q^{2}_{t}|f_{1:t-1},q^{1}_{1:t-1})\sim\mathcal{N}(\mu^{2}_{t},P^{+})$,
where $P^{+}=APA^{T}+B_{w}\Sigma_{w}B_{w}^{T}$ and the posterior distribution
on ${\bf f}_{t}$ can be written as:

$\begin{split}p(f_{t}|n=2,q^{1}_{1:t},f_{1:t-1})=\qquad\qquad\qquad\qquad\qquad\\\
\mathcal{N}(K_{1}(\mu^{2}_{t}-q^{1}_{t}),K_{1}P^{+}K_{1}^{T}+\Sigma_{f}),\end{split}$
(20)

updating the contact force distribution in (15).

##  5 Experimental Results

We test information gain in experiments with industrial robots (7 and 60 Kg
payload) and varying compliant conditions. An F/T sensor is integrated at the
flange and measured at $1250$ Hz, and a standard admittance control and
collision detector implemented, see [38] for details. When collision is
detected based on a force threshold, the position commanded to the robot is
held to the last position before collision was detected. The F/T sensors used
in all setups have similar noise properties, and $\Sigma_{f}=1.25$ is used for
all conditions.

The experimental data and analysis code is available at
https://owncloud.fraunhofer.de/index.php/s/ovtOlGRb7o7b9Fb.

###  5.1 Monitoring Assembly

We consider an admittance-controlled gear assembly task from [16], where the
gear assembly has multiple stages, as seen in Figure 1. One failure mode is
in-hand slip, which may occur during the gear mating process. We modify the
setup by introducing additional foam at the fingertips, between the gripper
and gear, and compare the ability to infer the task mode. The parameters of
the estimator are available with the experimental data in the above mentioned
cloud folder.

A typical mode detection from a trial with in-hand slip can be seen in Figure
4, where the belief entropy (i.e. mode uncertainty) increases substantially
during the transitions between the modes. Comparing the total entropy over a
number of experiment iterations, we find that decreasing the gripper stiffness
increases the total uncertainty, as seen in Figure 5.

The information gain is calculated using the fully-observed information gain
(17) with the gripper stiffness as $K_{1}$, and the average taken over
$|q^{1}-q^{z}2|\in[0,0.5]$ mm with a flat prior $b_{n}=1/N$. It can be seen
that the higher stiffness has higher information gain, corresponding with
lower total belief entropy.

Figure 4: Detail of transition in mode estimate, with estimated contact state
over time (top), force/position (middle) and belief entropy at each time step
(bottom).  Figure 5: Total belief entropy over compliance at the gripper. For
both successful and in hand slip trials, the uncertainty increases as the
compliance decreases.

###  5.2 Collision detection - High Payload Robot

A 16 Kg plate is mounted to a large industrial robot as seen in Figure 7, and
brought into contact with a compliant magazine, where compliance is introduced
in one of three locations: (i) in the joints of the magazine, (ii) on the
contact surface, and (iii) under the feet of the magazine. The robot was
driven into contact in using admittance control with a desired contact force
of 30 N. Pure admittance control is used to collect data for model
identification, then contact detection is activated, where a simple force
threshold of $6$ N is used (gravitational forces of the plate are
compensated).

Condition | $K_{1}$ | $M_{2}$ | $B_{2}$ | $K_{2}$ | $\underline{\Delta}_{I}$
---|---|---|---|---|---
Flex Joints | 17.4e4 | 20.1 | 305 | 2630 | 0.131
Compliant Surface | 1.3e4 | 63.9 | 1870 | 7350 | 0.0743
Compliant Feet | 2.61e4 | 69.3 | 1080 | 1.81e4 | 0.348
Flex Joints Gradient | 7.7e-6 | 7.6e-3 | 2.2e-6 | -4.8e-10 | -
Comp Surf Gradient | 1.5e-5 | 1.8e-3 | 1.7e-6 | 1.3e-10 | -
Comp Feet Gradient | 2.5e-5 | 0.8e-3 | 1.3e-6 | 4.8e-10 | -
Table 1: Identified dynamic parameters for the magazine contact  Figure 6:
Contact force, position, and velocity in experiments with the contact
detection

As seen in Figure 6, the fastest contact detection occurred with compliance on
the feet (bottom), where the hard contact and larger effective magazine
inertia results in an initial peak in force. Compliance in the joints results
in the 2nd fastest contact detection – note that in this case contact was lost
as the magazine springs back with the lower $K_{2}$. The slowest contact
detection was with compliance at the contact surface, where force increases
slowly.

A model for these three compliance cases is then identified following the
inertial model in Figure 3, where the measured robot position is $q^{1}$,
unmeasured magazine position is $q^{2}$ and the measured force is that on
$K_{1}$. These parameters are identified by doing a least-squares fit
$\min_{K_{1},M_{2},B_{2},K_{2}}\sum_{t}\|f^{m}_{t}-f_{t}\|$, where $f^{m}_{t}$
is the model output given the observed $q^{1}_{t}$ as input, and can be seen
in Table 1, with time plots in Figure 8.

From these models, the $\underline{\Delta}_{I}$ can be calculated in the
partially observed case with $q_{t}$ at time $t=0.002$. It can be seen in
Table 1 that the compliant feet has the highest information gain, followed by
the compliant joints and compliant surface, corresponding to the contact
detection performance. Because the $\underline{\Delta}_{I}$ calculation is
composed of differentiable operations, the gradient can be numerically
calculated with an automatic differentiation toolbox [39], also seen in Table
1. It can be seen that increasing $K_{1}$ and $M_{2}$ make the largest
increase in $\underline{\Delta}_{I}$, while $B_{2}$ and $K_{2}$ play a more
minor role. The gradient validates an experimental observation: a higher
contact inertia $M_{2}$ helps improve contact detection.

Figure 7: Contact magazine with the various locations of compliance, using
strips of viscoelastic foam, or 3d printed flexure joints.

Figure 8: Identified model and real response, for flex joint magazine (top)
and compliance in the feet (bottom). Contact is made with force control, with
desired force of 30 N. Note that admittance parameters are set poorly for the
bottom contact, but the identified model is still accurate.

###  5.3 Collision Detection - Medium Payload Robot

Figure 9: Compliant table with 7 Kg payload robot, where compliance can be
adjusted in the feet of the table or between the failsafe block and table.

A compliant table, seen in Figure 9, was used to investigate the impact of
location and quantity of compliance on the ability to detect contact. By
varying the compliance in the feet and contact surface, as well as the inertia
of the table, a range of contact responses were achieved as seen in Figure 10.
A simple contact detection threshold of $6$ N was applied, and the performance
can be seen in the attached video. Model identification was done as above, but
a two-inertia model was found to provide better model fit, requiring a
straightforward extension to (18) which is omitted due to space constraints
but in the published code.

Figure 10: Traces of contact force (solid line) vs identified model response
(dotted). Detailed responses can be better seen on the online repository,
these responses show the range of force profiles and model fit.

The information gain, calculated with (20) at the state $t=0.002$ after
contact, is compared with the time taken to detect the contact in Figure 11.
It can be seen that as information gain increases, time to detect contact
decreases.

The configurations where contact compliance was removed, and where additional
weight was added to the table, had higher information gain and a higher
contact velocity could be safely achieved - further detailed plots are
available online with the experimental data.

Figure 11: Information gain vs time to detect contact, colors are consistent
with those used in Figure 10. As information gain increases, the time required
to detect contact decreases. Hard contact and higher table inertia gave the
fastest detection.

##  6 Conclusion

The impact of mechatronic design, especially compliance, on the ability to
detect contact conditions was presented, and formalized in a metric of the
information gain. This metric is validated to predict the certainty with a
complex assembly task can be monitored, and the speed with which collision can
be detected.

In the case of introducing compliance to the environment, it is shown that
this compliance is best introduced not at the point of contact, but deeper
into the kinematic structure. This allowed faster contact detection in both
the magazine and table contact tasks.

## References

  * [1] G. Bledt, P. M. Wensing, S. Ingersoll, and S. Kim, “Contact Model Fusion for Event-Based Locomotion in Unstructured Terrains,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , May 2018, pp. 4399–4406.
  * [2] M. Vecerik, O. Sushkov, D. Barker, T. Rothörl, T. Hester, and J. Scholz, “A Practical Approach to Insertion with Variable Socket Position Using Deep Reinforcement Learning,” _arXiv:1810.01531 [cs]_ , Oct. 2018\.
  * [3] D. E. Stewart, “Rigid-Body Dynamics with Friction and Impact,” _SIAM Review_ , vol. 42, no. 1, pp. 3–39, Jan. 2000.
  * [4] A. M. Castro, A. Qu, N. Kuppuswamy, A. Alspach, and M. Sherman, “A Transition-Aware Method for the Simulation of Compliant Contact With Regularized Friction,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 1859–1866, Apr. 2020.
  * [5] M. V. Minniti, R. Grandia, K. Fäh, F. Farshidian, and M. Hutter, “Model predictive robot-environment interaction control for mobile manipulation tasks,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 1651–1657.
  * [6] M. Posa, C. Cantu, and R. Tedrake, “A direct method for trajectory optimization of rigid bodies through contact,” _The International Journal of Robotics Research_ , vol. 33, no. 1, pp. 69–81, 2014.
  * [7] Feng Pan and J. M. Schimmels, “Robust procedures for obtaining assembly contact state extremal configurations,” in _IEEE International Conference on Robotics and Automation, 2004. Proceedings. ICRA ’04. 2004_ , vol. 1, Apr. 2004, pp. 972–977 Vol.1.
  * [8] J. De Schutter, T. De Laet, J. Rutgeerts, W. Decré, R. Smits, E. Aertbeliën, K. Claes, and H. Bruyninckx, “Constraint-based task specification and estimation for sensor-based robot systems in the presence of geometric uncertainty,” _The International Journal of Robotics Research_ , vol. 26, no. 5, pp. 433–455, 2007.
  * [9] M. P. Polverini, A. M. Zanchettin, and P. Rocco, “A constraint-based programming approach for robotic assembly skills implementation,” _Robotics and Computer-Integrated Manufacturing_ , vol. 59, pp. 69–81, 2019\.
  * [10] W. Meeussen, J. Rutgeerts, K. Gadeyne, H. Bruyninckx, and J. De Schutter, “Contact-state segmentation using particle filters for programming by human demonstration in compliant-motion tasks,” _IEEE Transactions on Robotics_ , vol. 23, no. 2, pp. 218–231, 2007.
  * [11] S. Pfrommer, M. Halm, and M. Posa, “ContactNets: Learning of Discontinuous Contact Dynamics with Smooth, Implicit Representations,” _arXiv preprint arXiv:2009.11193_ , 2020.
  * [12] J. Rosell, L. Basanez, and R. Suárez, “Contact Identification for robotic assembly tasks with uncertainty,” _IFAC Proceedings Volumes_ , vol. 33, no. 27, pp. 141–146, 2000.
  * [13] L. Roveda, “Sensorless environment stiffness and interaction force estimation for impedance control tuning in robotized interaction tasks,” _Autonomous Robots_ , p. 18.
  * [14] I. F. Jasim and P. W. Plapper, “Contact-state recognition of compliant motion robots using expectation maximization-based Gaussian mixtures,” in _ISR/Robotik 2014; 41st International Symposium on Robotics; Proceedings Of_. VDE, 2014, pp. 1–8.
  * [15] T. Ren, Y. Dong, D. Wu, and K. Chen, “Collision detection and identification for robot manipulators based on extended state observer,” _Control Engineering Practice_ , vol. 79, pp. 144–153, 2018.
  * [16] K. Haninger and D. Surdilovic, “Multimodal Environment Dynamics for Interactive Robots: Towards Fault Detection and Task Representation,” in _Proc. IEEE/RSJ Intl Conf on Intelligent Robots and Systems (IROS)_ , 2018, pp. 6932–6937.
  * [17] S. Golz, C. Osendorfer, and S. Haddadin, “Using tactile sensation for learning contact knowledge: Discriminate collision from physical interaction,” in _2015 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2015, pp. 3788–3794.
  * [18] P. M. Wensing, A. Wang, S. Seok, D. Otten, J. Lang, and S. Kim, “Proprioceptive Actuator Design in the MIT Cheetah: Impact Mitigation and High-Bandwidth Physical Interaction for Dynamic Legged Robots,” _IEEE Transactions on Robotics_ , vol. 33, no. 3, pp. 509–522, Jun. 2017.
  * [19] A. Bicchi and G. Tonietti, “Fast and ”Soft-Arm” Tactics,” _IEEE Robotics & Automation Magazine_, vol. 11, no. 2, pp. 22–33, Jun. 2004.
  * [20] K. Haninger and D. Surdilovic, “Bounded Collision Force by the Sobolev Norm: Compliance and Control for Interactive Robots,” in _2019 IEEE International Conference on Robotics and Automation (ICRA)_ , 2019, pp. 8259–8535.
  * [21] S. D. Eppinger and W. P. Seering, “Three dynamic problems in robot force control,” _IEEE Transactions on Robotics and Automation_ , vol. 8, no. 6, pp. 751–758, 1992.
  * [22] G. Pratt and M. Williamson, “Series elastic actuators,” in _1995 IEEE/RSJ International Conference on Intelligent Robots and Systems 95. ’Human Robot Interaction and Cooperative Robots’, Proceedings_ , vol. 1, Aug. 1995, pp. 399–406 vol.1.
  * [23] E. A. Bolıvar-Nieto, G. C. Thomas, E. Rouse, and R. D. Gregg, “Convex Optimization for Spring Design in Series Elastic Actuators: From Theory to Practice,” 2021.
  * [24] S. Haddadin, K. Krieger, A. Albu-Schäffer, and T. Lilge, “Exploiting elastic energy storage for “blind” cyclic manipulation: Modeling, stability analysis, control, and experiments for dribbling,” _IEEE Transactions on Robotics_ , vol. 34, no. 1, pp. 91–112, 2018.
  * [25] K. Haninger, “Minimum directed information: A design principle for compliant robots,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2021.
  * [26] M. Hamaya, R. Lee, K. Tanaka, F. von Drigalski, C. Nakashima, Y. Shibata, and Y. Ijiri, “Learning robotic assembly tasks with lower dimensional systems by leveraging physical softness and environmental constraints,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 7747–7753.
  * [27] P. Tomei, “A simple PD controller for robots with elastic joints,” _Automatic Control, IEEE Transactions on_ , vol. 36, no. 10, pp. 1208–1213, 1991.
  * [28] C. Della Santina, M. Bianchi, G. Grioli, F. Angelini, M. Catalano, M. Garabini, and A. Bicchi, “Controlling Soft Robots: Balancing Feedback and Feedforward Elements,” _IEEE Robotics Automation Magazine_ , vol. 24, no. 3, pp. 75–83, Sep. 2017.
  * [29] T. J. Debus, P. E. Dupont, and R. D. Howe, “Distinguishability and identifiability testing of contact state models,” _Advanced Robotics_ , vol. 19, no. 5, pp. 545–566, 2005.
  * [30] M. W. Spong, “Modeling and control of elastic joint robots,” _Journal of dynamic systems, measurement, and control_ , vol. 109, no. 4, pp. 310–318, 1987\.
  * [31] N. Kashiri, J. Malzahn, and N. G. Tsagarakis, “On the Sensor Design of Torque Controlled Actuators: A Comparison Study of Strain Gauge and Encoder-Based Principles,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 2, pp. 1186–1194, Apr. 2017.
  * [32] S. Thrun, “Probabilistic robotics,” _Communications of the ACM_ , vol. 45, no. 3, pp. 52–57, 2002.
  * [33] S. Haddadin, A. De Luca, and A. Albu-Schäffer, “Robot collisions: A survey on detection, isolation, and identification,” _IEEE Transactions on Robotics_ , vol. 33, no. 6, pp. 1292–1312, 2017.
  * [34] T. M. Cover, _Elements of Information Theory_. John Wiley & Sons, 1999.
  * [35] A. Kolchinsky and B. D. Tracey, “Estimating mixture entropy with pairwise distances,” _Entropy_ , vol. 19, no. 7, p. 361, 2017.
  * [36] M. F. Huber, T. Bailey, H. Durrant-Whyte, and U. D. Hanebeck, “On entropy approximation for Gaussian mixture random vectors,” in _2008 IEEE International Conference on Multisensor Fusion and Integration for Intelligent Systems_. IEEE, 2008, pp. 181–188.
  * [37] T.-C. Kao and G. Hennequin, “Automatic differentiation of Sylvester, Lyapunov, and algebraic Riccati equations,” _arXiv:2011.11430 [cs, math]_ , Nov. 2020.
  * [38] K. Haninger, M. Radke, A. Vick, and J. Kruger, “Towards High-Payload Admittance Control for Manual Guidance with Environmental Contact,” _IEEE Robotics and Automation Letters_ , 2022.
  * [39] J. A. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, “CasADi: A software framework for nonlinear optimization and optimal control,” _Mathematical Programming Computation_ , vol. 11, no. 1, pp. 1–36, 2019.