# On the computation of a non-parametric estimator by convex optimization

Akshay Seshadri Email<EMAIL_ADDRESS>Department of Physics,
University of Colorado Boulder  Stephen Becker  Department of Applied
Mathematics, University of Colorado Boulder

###### Abstract

Estimation of linear functionals from observed data is an important task in
many subjects. Juditsky & Nemirovski [The Annals of Statistics 37.5A (2009):
2278-2300] propose a framework for non-parametric estimation of linear
functionals in a very general setting, with nearly minimax optimal confidence
intervals. They compute this estimator and the associated confidence interval
by approximating the saddle-point of a function. While this optimization
problem is convex, it is rather difficult to solve using existing off-the-
shelf optimization software. Furthermore, this computation can be expensive
when the estimators live in a high-dimensional space. We propose a different
algorithm to construct this estimator. Our algorithm can be used with existing
optimization software and is much cheaper to implement even when the
estimators are in a high-dimensional space, as long as the Hellinger affinity
(or the Bhattacharyya coefficient) for the chosen parametric distribution can
be efficiently computed given the parameters. We hope that our algorithm will
foster the adoption of this estimation technique to a wider variety of
problems with relative ease.

There are many situations where one wishes to estimate linear functionals of
an unknown state using only observations of quantities determined by the state
(i.e., indirect measurements of the state). Such a scenario is prevalent, for
example, in tomography. Therefore, one would ideally like to incorporate these
different measurements in such a way that not only we find a good estimate of
the linear functional, but also ensure that the associated confidence interval
is tight. Juditsky & Nemirovski [1], extending prior work of Donoho [2],
propose an approach to this problem that constructs an estimator for a
specified linear functional algorithmically, incorporating these indirect
measurements. The (symmetric) confidence interval associated with their
estimator is guaranteed to be nearly minimax optimal. Furthermore, once a
model for the system has been specified, along with the number and type of
measurements that will be recorded, their method can construct the estimator
and the confidence interval even before seeing any data. This fact can be
useful, for example, when one wishes to minimize the measurements that need to
be performed to achieve a desired size of the confidence interval.

To construct this estimator, Juditsky & Nemirovski [1] propose an algorithm
that involves computing the saddle-point of a concave-convex function to a
given precision. While this optimization problem is convex, solving it using
standard optimization algorithms or off-the-shelf optimization software like
CVX [3] is difficult in practice. In some cases, it is possible to extend the
capabilities of these software to handle such saddle-point problems [4].
However, even in such situations, the fact remains that one might need to
perform a high-dimensional optimization for constructing the estimator. These
reasons make the estimation technique hard to use in practice.

In this study, we propose a different algorithm to compute the estimator that
overcomes these difficulties. While our algorithm also approximates the
saddle-point, it can be implemented using off-the-shelf optimization software
currently available. Moreover, our algorithm replaces the potentially high-
dimensional optimization with a much cheaper function evaluation, thus greatly
reducing the computational burden. Furthermore, we write the estimator in a
form that is more amenable to interpretation, which could serve as a starting
point for further theoretical investigations about the estimator. This
algorithm might also be useful in the extensions of Juditsky & Nemirovski’s
framework [1] to several other problems [5, 6].

We begin by reviewing the estimation framework proposed by Juditsky &
Nemirovski [1]. We then present our approach to compute the estimator, along
with the relevant proofs.

##  1 Non-parametric estimation framework

Suppose that we are given a set of “states”
$\mathcal{X}\subseteq\mathbb{R}^{d}$ that is a compact and convex subset of
$\mathbb{R}^{n}$. We wish to estimate the linear functional $g^{T}x$, where
$g\in\mathbb{R}^{n}$ is some fixed vector, while the state $x\in\mathcal{X}$
of the system is unknown to us. While we don’t know $x$, we have access to a
single measurement outcome determined by $x$. Measurements are modeled using
random variables that assign probabilities to the possible outcomes depending
on the state.

To that end, Juditsky & Nemirovski [1] consider a family of random variables
$\mathrm{Z}_{\mu}$ parameterized by $\mu\in\mathcal{M}$, where
$\mathcal{M}\subseteq\mathbb{R}^{m}$ is some set of parameters. These random
variables take values in a separable complete metric space (or Polish space)
$(\Omega,\Sigma)$ equipped with a $\sigma$-finite Borel measure $\mathbb{P}$
[1]. We assume that $\mathrm{Z}_{\mu}$ has a probability density $p_{\mu}$
with respect to this reference measure $\mathbb{P}$. The state
$x\in\mathcal{X}$ determines the random variable $\mathrm{Z}_{A(x)}$ through
an affine function $A\colon\mathcal{X}\to\mathcal{M}$, and we are given one
outcome of this random variable for the purpose of estimation.

Our goal is to construct an estimator for $g^{T}x$ that uses an outcome of the
random variable $\mathrm{Z}_{A(x)}$ to give an estimate. An estimator is a
real-valued Borel measurable function on $\Omega$. The set of estimators
$\mathcal{F}$ we are allowed to work with is any finite-dimensional vector
space comprised of real-valued Borel measurable functions on $\Omega$ as long
as it contains constant functions [1]. The mapping $\mathcal{D}(\mu)=p_{\mu}$
between the parameter $\mu$ and the corresponding probability density
$p_{\mu}$ is called a parametric density family [1]. In order to be able to
choose an appropriate estimator from $\mathcal{F}$ given that the probability
density of the random variable is $p_{A(x)}$, we want the set of estimators
$\mathcal{F}$ and the parametric density family $\mathcal{D}$ to interact well
with each other. This gives rise to the notion of a good pair defined by
Juditsky & Nemirovski [1].

######  Definition 1 (Good pair).

We call a given pair $(\mathcal{D},\mathcal{F})$ of parametric density family
$\mathcal{D}$ and finite-dimensional space $\mathcal{F}$ of Borel functions on
$\Omega$ a good pair if the following conditions hold.

  1. 1.

$\mathcal{M}$ is a relatively open convex set in $\mathbb{R}^{m}$.

  2. 2.

Whenever $\mu\in\mathcal{M}$, we have $p_{\mu}(\omega)>0$ for all
$\omega\in\Omega$.

  3. 3.

Whenever $\mu,\nu\in\mathcal{M}$,
$\phi(\omega)=\ln(p_{\mu}(\omega)/p_{\nu}(\omega))\in\mathcal{F}$.

  4. 4.

Whenever $\phi\in\mathcal{F}$, the function

$F_{\phi}(\mu)=\ln\left(\int_{\Omega}\exp\left(\phi(\omega)\right)p_{\mu}(\omega)\mathbb{P}(d\omega)\right)$

is well-defined and concave in $\mu\in\mathcal{M}$.

Any estimator $\widehat{g}\in\mathcal{F}$ is called an affine estimator (note,
however, that $\widehat{g}$ need not be an affine function in general).

To judge the performance of an arbitrary estimator $\widehat{g}$, we define
the $\epsilon$-risk as follows [1].

######  Definition 2 ($\epsilon$-risk).

Given a confidence level $1-\epsilon\in(0,1)$, we define the $\epsilon$-risk
associated with an estimator $\widehat{g}$ as

$\mathcal{R}(\widehat{g};\epsilon)=\inf\left\\{\delta:\sup_{x\in\mathcal{X}}\operatorname*{Prob}_{\omega\sim
p_{A(x)}}\left\\{\omega:|\widehat{g}(\omega)-g^{T}x|>\delta\right\\}<\epsilon\right\\},$

where $\omega\sim p_{A(x)}$ means that $\omega$ is sampled according to
$p_{A(x)}$. The corresponding minimax optimal risk is defined as

$\mathcal{R}_{*}(\epsilon)=\inf_{\widehat{g}}\mathcal{R}(\widehat{g};\epsilon)$

where the infimum is taken over all Borel functions $\widehat{g}$ on $\Omega$.
Restricting to just the affine estimators, the affine risk is defined as

$\mathcal{R}_{\text{aff}}(\epsilon)=\inf_{\widehat{g}\in\mathcal{F}}\mathcal{R}(\widehat{g};\epsilon).$

Intuitively, the $\epsilon$-risk of an estimator gives the smallest possible
additive error of the estimator for the given confidence level, independent of
the actual state of the system or the observed measurement outcome.

Optimizing over all possible estimators to find the best estimator can be
computationally intractable. For this reason, Juditsky & Nemirovski focus on
constructing an affine estimator. It turns out that we don’t lose much by
restricting our attention to affine estimators. Indeed, Juditsky & Nemirovski
[1] prove that if $(\mathcal{D},\mathcal{F})$ is a good pair, then there is an
estimator $\widehat{g}_{*}\in\mathcal{F}$ with $\epsilon$-risk at most
$\Phi_{*}(\ln(2/\epsilon))$, such that

$\displaystyle\mathcal{R}_{\text{aff}}(\epsilon)$
$\displaystyle\leq\Phi_{*}(\ln(2/\epsilon))\leq\vartheta(\epsilon)\mathcal{R}_{*}(\epsilon)$
$\displaystyle\vartheta(\epsilon)$
$\displaystyle=2+\frac{\ln(64)}{\ln(0.25/\epsilon)}$

for $\epsilon\in(0,0.25)$.

Juditsky & Nemirovski [1] propose the following method to construct the
estimator $\widehat{g}_{*}$ and the risk $\Phi_{*}(\ln(2/\epsilon))$ (to a
precision $2\delta>0$).

  1. 1.

For $r\geq 0$, define the function
$\Phi_{r}\colon(\mathcal{X}\times\mathcal{X})\times(\mathcal{F}\times\mathbb{R}_{+})\to\mathbb{R}$
as

$\displaystyle\Phi_{r}(x,y;\ \phi,\alpha)=g^{T}x-g^{T}y\ +\ $
$\displaystyle\alpha\Bigg{[}\ln\left(\int_{\Omega}\exp(-\phi(\omega)/\alpha)p_{A(x)}(\omega)\mathbb{P}(d\omega)\right)$
$\displaystyle\hskip
7.11317pt+\ln\left(\int_{\Omega}\exp(\phi(\omega)/\alpha)p_{A(y)}(\omega)\mathbb{P}(d\omega)\right)\Bigg{]}+2\alpha
r.$ (1)

Juditsky & Nemirovski [1] show that $\Phi_{r}$ has the following properties.
$\Phi_{r}$ is continuous and concave in $(x,y)$, continuous and convex in
$(\phi,\alpha)$, and we have $\Phi_{r}\geq 0$. Furthermore, $\Phi_{r}$ has a
well-defined saddle-point value.

  2. 2.

Denote the saddle-point value of $\Phi_{r}$ by $2\Phi_{*}(r)$:

$\Phi_{*}(r)=\frac{1}{2}\sup_{x,y\in\mathcal{X}}\inf_{\phi\in\mathcal{F},\alpha>0}\Phi_{r}(x,y;\phi,\alpha)=\frac{1}{2}\inf_{\phi\in\mathcal{F},\alpha>0}\max_{x,y\in\mathcal{X}}\Phi_{r}(x,y;\phi,\alpha).$
(2)

  3. 3.

Given a confidence level $1-\epsilon\in(0.75,1)$ and a precision $\delta>0$,
find $\phi_{*}\in\mathcal{F}$ and $\alpha_{*}>0$ such that

$\max_{x,y\in\mathcal{X}}\Phi_{\ln(2/\epsilon)}(x,y;\phi_{*},\alpha_{*})\leq
2\Phi_{*}(\ln(2/\epsilon))+\delta.$

This is achieved by minimizing the convex function

$\overline{\Phi}_{\ln(2/\epsilon)}(\phi,\alpha)=\max_{x,y\in\mathcal{X}}\Phi_{\ln(2/\epsilon)}(x,y;\phi,\alpha).$

  4. 4.

The estimator $\widehat{g}_{*}\in\mathcal{F}$ is then defined as

$\widehat{g}_{*}=\phi_{*}+c$ (3)

where the constant $c$ is obtained by solving the optimization problem

$\displaystyle c$ $\displaystyle=\frac{1}{2}\max_{\chi\in
X}\left[g^{T}x+\alpha_{*}\ln\left(\int_{\Omega}\exp\left(-\phi_{*}/\alpha_{*}\right)p_{A(x)}\mathbb{P}(d\omega)\right)\right]$
$\displaystyle\hskip 8.5359pt-\frac{1}{2}\max_{y\in
X}\left[-g^{T}y+\alpha_{*}\ln\left(\int_{\Omega}\exp\left(\phi_{*}/\alpha_{*}\right)p_{A(y)}\mathbb{P}(d\omega)\right)\right].$
(4)

Then, the $\epsilon$-risk of the estimator $\widehat{g}_{*}$ is bounded above
by $\Phi_{*}(\ln(2/\epsilon))+2\delta$.

Given an observation $\omega$ of $\mathrm{Z}_{A(x)}$, our estimate for
$g^{T}x$ is given by $\widehat{g}_{*}(\omega)$ with an additive error of
$\Phi_{*}(\ln(2/\epsilon))+2\delta$ and a confidence level of $1-\epsilon$.
Note that computing the estimator $\widehat{g}_{*}$ requires one to perform
optimization and can be time consuming, but once the estimator has been
computed, the estimates can be obtained using $\widehat{g}_{*}$ almost
instantaneously.

So far we have described how to find an estimate for $g^{T}x$ from one outcome
of a single random variable $\mathrm{Z}_{A(x)}$. In practice, we will need to
consider many outcomes of different random variables
$\mathrm{Z}_{A^{(l)}(x)}$, which corresponds to $l=1,\dotsc,L$ different types
of measurement. More precisely, we are given Polish spaces
$(\Omega^{(l)},\Sigma^{(l)})$ equipped with a $\sigma$-finite Borel measure
$\mathbb{P}^{(l)}$ for $l=1,\dotsc,L$. We are also given a set of parameters
$\mathcal{M}^{(l)}$ for $l=1,\dotsc,L$. For each $l=1,\dotsc,L$, we are given
a family of random variables $\mathrm{Z}_{\mu_{l}}$ taking values in
$\Omega^{(l)}$, where $\mu_{l}\in\mathcal{M}^{(l)}$. The random variable
$\mathrm{Z}_{\mu_{l}}$ has a probability density $p^{(l)}_{\mu_{l}}$ with
respect to the reference measure $\mathbb{P}^{(l)}$. As before, we are given
affine mappings $A^{(l)}\colon\mathcal{X}\to\mathcal{M}^{(l)}$ for
$l=1,\dotsc,L$ that map the state $x\in\mathcal{X}$ of the system to a
corresponding parameter. For each $l=1,\dotsc,L$, we can choose estimators for
the $l^{\text{th}}$ measurement from the set $\mathcal{F}^{(l)}$, which is a
finite-dimensional vector space of real-valued Borel measurable functions on
$\Omega^{(l)}$ that contains constant functions. To incorporate the outcomes
of these different random variables, Juditsky & Nemirovski [1] define the
direct product of good pairs, which essentially constructs one large good pair
from many smaller ones.

######  Definition 3 (Direct product of good pairs).

Considering the following quantities for $l=1,\dotsc,L$. Let
$(\Omega^{(l)},\Sigma^{(l)})$ be a Polish space endowed with a Borel
$\sigma$-finite measure $\mathbb{P}^{(l)}$. Let
$\mathcal{D}^{(l)}(\mu_{l})=p^{(l)}_{\mu_{l}}$ be the parametric density
family for $\mu_{l}\in\mathcal{M}^{(l)}$. Let $\mathcal{F}^{(l)}$ be a finite-
dimensional linear space of Borel functions on $\Omega^{(l)}$ containing
constants, such that the pair $(\mathcal{D}^{(l)},\mathcal{F}^{(l)})$ is good.
Then the direct product of these good pairs
$(\mathcal{D},\mathcal{F})=\bigotimes_{l=1}^{L}(\mathcal{D}^{(l)},\mathcal{F}^{(l)})$
is defined as follows.

  1. 1.

The large space is $\Omega=\Omega^{(1)}\times\dotsb\times\Omega^{(L)}$ endowed
with the product measure
$\mathbb{P}=\mathbb{P}^{(1)}\times\dotsb\times\mathbb{P}^{(L)}$.

  2. 2.

The set of parameters is
$\mathcal{M}=\mathcal{M}^{(1)}\times\dotsb\times\mathcal{M}^{(L)}$, and the
associated parametric density family is
$\mathcal{D}(\mu)=p_{\mu}\equiv\prod_{l=1}^{L}p^{(l)}_{\mu_{l}}$ for
$\mu=(\mu_{1},\dotsc,\mu_{L})\in\mathcal{M}$.

  3. 3.

The linear space $\mathcal{F}$ comprises of all functions $\phi$ defined as
$\phi(\omega_{1},\omega_{2},\dotsc,\omega_{L})=\sum_{l=1}^{L}\phi^{(l)}(\omega_{l})$,
where $\phi^{(l)}\in\mathcal{F}^{(l)}$ and $\omega_{l}\in\Omega^{(l)}$ for
$l=1,\dotsc,L$.

It can be verified that the direct product of good pairs is a good pair [1].

To obtain an estimator that accounts for all the given measurement outcomes,
we can apply the procedure outlined for constructing the estimator
$\widehat{g}_{*}$ for the single outcome case to the direct product of good
pairs. Specifically, if we observe $R_{l}$ outcomes of the random variable
$\mathrm{Z}_{l}$ for $l=1,\dotsc,L$, we need to minimize the function

$\overline{\Phi}_{\ln(2/\epsilon)}(\phi,\alpha)=\max_{x,y\in\mathcal{X}}\Phi_{\ln(2/\epsilon)}(x,y;\phi,\alpha)$
(5)

where

$\displaystyle\Phi_{\ln(2/\epsilon)}(x,y;\phi,\alpha)=g^{T}x-g^{T}y\ +\ $
$\displaystyle\alpha\sum_{l=1}^{L}R_{l}\Bigg{[}\ln\left(\int_{\Omega^{(l)}}\exp(-\phi^{(l)}(\omega)/\alpha)p_{A(x)}(\omega)\mathbb{P}^{(l)}(d\omega)\right)$
$\displaystyle\hskip
7.11317pt+\ln\left(\int_{\Omega^{(l)}}\exp(\phi^{(l)}(\omega)/\alpha)p_{A(y)}(\omega)\mathbb{P}^{(l)}(d\omega)\right)\Bigg{]}+2\alpha\ln(2/\epsilon)$
(6)

and $\phi=(\phi^{(1)},\dotsc,\phi^{(L)})$.

###  1.1 Computational hurdles in finding an approximation to the saddle-point

We now point out several difficulties in performing the minimization of
$\overline{\Phi}_{\ln(2/\epsilon)}(\phi,\alpha)$ defined in Eq. (5) over
$\phi\in\mathcal{F}$ and $\alpha>0$. Note that this optimization is necessary
to construct an estimator using Juditsky & Nemirovski’s approach.

  1. 1.

When we have $L$ types of measurements $\mathrm{Z}_{1},\dotsc,\mathrm{Z}_{L}$,
$\phi$ is a vector of the form $\phi=(\phi^{(1)},\dotsc,\phi^{(L)})$. Since
each $\phi^{(l)}$ is a vector by construction, $\phi$ can be a high-
dimensional vector. This can make the minimization
$\inf_{\alpha>0}\inf_{\phi\in\mathcal{F}}\overline{\Phi}_{\ln(2/\epsilon)}(\phi,\alpha)$
costly to implement.

We illustrate through a simple example that the estimator $\phi$ can be a
high-dimensional vector even with only one type of measurement, i.e., $L=1$.
The following is Example 1 from Juditsky & Nemirovski [1].

######  Example 1 (Discrete distributions).

Suppose that the set of parameters
$\mathcal{M}=\\{\mu\in\mathbb{R}^{n}\mid\mu>0,\sum_{i=1}^{n}\mu_{i}=1\\}$
characterizing the observations is the relatively open standard simplex. Let
$\mathcal{D}$ be the mapping $\mathcal{D}(\mu)=\mu$ for any parameter
$\mu\in\mathcal{M}$. In other words, the probability distribution $p_{\mu}$
determined by the parameter $\mu$ is just $p_{\mu}=\mu$. Suppose that the set
of observations is $\Omega=\\{1,\dotsc,n\\}$, and given a parameter
$\mu\in\mathcal{M}$, the element $i\in\Omega$ is observed with a probability
of $\mu_{i}$. We take $\Sigma$ to be the discrete $\sigma$-algebra and
$\mathbb{P}$ to be the counting measure.

Any estimator $\phi\colon\Omega\to\mathbb{R}$ is an $n$-dimensional real
vector. In particular, $\mathcal{F}=\mathbb{R}^{n}$ is a natural choice for
the set of estimators, such that $(\mathcal{D},\mathcal{F})$ is a good pair.
For any system with a large number of observations $n$, the estimator $\phi$
is a high-dimensional vector.

  2. 2.

Since $\overline{\Phi}_{\ln(2/\epsilon)}$ is itself a maximum of the function
$\Phi_{\ln(2/\epsilon)}$, popular gradient based methods can be difficult to
use for minimizing $\overline{\Phi}_{\ln(2/\epsilon)}(\phi,\alpha)$ over
$\phi\in\mathcal{F}$ and $\alpha>0$, even when the function
$\Phi_{\ln(2/\epsilon)}(x,y;\phi,\alpha)$ is smooth in $\phi$ and $\alpha$. We
remark that recent progress has been made to extend the ability of CVX to
handle such problems when the function $\Phi_{\ln(2/\epsilon)}$ possesses
appropriate properties [4]. Nevertheless, solving saddle-point problems on CVX
will still require more computational effort (at least internally) than our
proposed algorithm, which can work with the current capabilities of CVX.

  3. 3.

In some cases, the computation of the function
$\overline{\Phi}_{\ln(2/\epsilon)}(\phi,\alpha)$ can itself be costly. For
example, suppose that $\mathcal{X}\subseteq\mathbb{R}^{d\times d}$ is the set
of positive semidefinite matrices with a trace constraint. Even for a
moderately large $d$, evaluating the function
$\overline{\Phi}_{\ln(2/\epsilon)}(\phi,\alpha)$ can be costly for a given
$\phi$ and $\alpha$. Therefore, we would like to reduce the number of calls to
the function $\overline{\Phi}_{\ln(2/\epsilon)}(\phi,\alpha)$.

Our approach to constructing the estimator circumvents the optimization over
$\phi\in\mathcal{F}$, thereby addressing the above issues to a reasonable
extent. We essentially solve the saddle point system from the other direction
and then construct the estimator, but this requires new justification because
the $\phi$-component of the saddle-point is not unique, and Ref. [1] uses a
specific choice of $\phi$ that is constructed algorithmically.

##  2 A different approach to constructing the estimator

We begin by presenting our algorithm for constructing the estimator. The
justification for the algorithm is given soon after.

  1. 1.

Given a confidence level $1-\epsilon\in(0.75,1)$ and a precision $\delta>0$,
find an $\alpha_{*}>0$, and subsequently $x^{*},y^{*}\in\mathcal{X}$, such
that

$2\Phi_{*}(\ln(2/\epsilon))+\delta\geq
2\alpha_{*}\ln(2/\epsilon)+\left(g^{T}x^{*}-g^{T}y^{*}+2\alpha_{*}\ln(\textnormal{AffH}(A(x),A(y)))\right)$
(7)

by approximately solving the convex optimization problem

$2\Phi_{*}(\ln(2/\epsilon))=\inf_{\alpha>0}\left[2\alpha\ln(2/\epsilon)+\max_{x,y\in
X}\left(g^{T}x-g^{T}y+2\alpha\ln(\textnormal{AffH}(A(x),A(y)))\right)\right]$

where

$\textnormal{AffH}(\mu,\nu)=\int_{\Omega}\sqrt{p_{\mu}(\omega)p_{\nu}(\omega)}\mathbb{P}(d\omega)$
(8)

is the Hellinger affinity between the distributions $p_{\mu}$ and $p_{\nu}$.

  2. 2.

Then, compute the corresponding $\phi$-component of the approximation to the
saddle-point using

$\phi_{*}=\frac{\alpha_{*}}{2}\ln\left(\frac{p_{A(x^{*})}}{p_{A(y^{*})}}\right)$
(9)

  3. 3.

The estimator $\widehat{g}_{*}$ is then readily obtained by setting

$\widehat{g}_{*}=\phi_{*}+\frac{1}{2}\left(g^{T}x^{*}+g^{T}y^{*}\right)$ (10)

The $\epsilon$-risk of $\widehat{g}_{*}$ is bounded above by
$\Phi_{*}(\ln(2/\epsilon))+\delta$.

Observe that our algorithm does not require one to compute the minimum over
$\phi\in\mathcal{F}$, thus avoiding the complications listed in Sec. (1.1).
The inner (continuous) concave maximization over $x,y\in\mathcal{X}$ can often
be performed using standard optimization algorithms or off-the-shelf
optimization software like CVX, assuming that $\mathcal{X}$ is computationally
tractable and the Hellinger affinity given in Eq. (8) can be efficiently
computed. The outer convex minimization is one-dimensional, and therefore, one
can use standard numerical routines (like, for example, those available in
scipy.optimize.minimize_scalar) to solve this problem efficiently.
Furthermore, the estimator given in Eq. (10) is appealing from a theoretical
standpoint because we have a closed-form expression in terms of the saddle-
point approximation $(x^{*},y^{*})$ and $\alpha^{*}$.

A potential drawback of our algorithm is that the Hellinger affinity between
$p_{\mu}$ and $p_{\nu}$ must be efficiently computable given $\mu$ and $\nu$.
For commonly encountered distributions, like the discrete distribution,
Poisson distribution, and the Gaussian distribution, one can find closed-form
expressions for the Hellinger affinity, and therefore,
$\textnormal{AffH}(\mu,\nu)$ can be efficiently computed for appropriate
choice of the parameters $\mu$ and $\nu$. Note that a similar problem exists
in Juditsky & Nemirovski’s [1] approach, wherein

$\int_{\Omega}e^{-\phi(\omega)/\alpha}p_{\mu}(\omega)\mathbb{P}(d\omega)$

must be efficiently computable given the parameter $\mu$. In a sense, the
requirement of our algorithm might be easier to satisfy because the Hellinger
affinity is a well-studied quantity, and furthermore, it is independent of the
set $\mathcal{F}$ of affine estimators one chooses for the problem.

We now provide justification of our algorithm. Our algorithm differs from that
of Juditsky & Nemirovski [1] on two fronts. One, we approximate the saddle-
point by solving a different optimization problem, which is based on the
expression for saddle-point value obtained in Prop. (1). Two, we make a choice
for the $\phi$-component of the saddle-point which, as noted in the proof of
Lemma (1), is possible because the $\phi$-component of the saddle-point is not
unique. We show in Prop. (2) that this approach still gives an estimator with
the same guarantees on the $\epsilon$-risk as the estimator constructed by
Juditsky & Nemirovski’s [1] algorithm.

We begin by proving the existence of saddle-point in $(x,y)$ and $\phi$
components for a fixed $\alpha>0$.

######  Lemma 1.

Let $\alpha>0$ be fixed. Then, for $r\geq 0$, the function

$\displaystyle\Phi^{\alpha}_{r}(x,y;\phi)=g^{T}x-g^{T}y\ +\ $
$\displaystyle\alpha\Bigg{[}\ln\left(\int_{\Omega}\exp(-\phi(\omega)/\alpha)p_{A(x)}(\omega)\mathbb{P}(d\omega)\right)$
$\displaystyle\hskip
7.11317pt+\ln\left(\int_{\Omega}\exp(\phi(\omega)/\alpha)p_{A(y)}(\omega)\mathbb{P}(d\omega)\right)\Bigg{]}+2\alpha
r$ (11)

defined on $(\mathcal{X}\times\mathcal{X})\times\mathcal{F}$ has a saddle-
point (maximum in $(x,y)$ and minimum in $\phi$). Furthermore, if
$(x^{*},y^{*})$ is the $(x,y)$ component of the saddle-point, then the
$\phi$-component of the saddle-point can be expressed as

$\phi_{*}=\frac{\alpha}{2}\ln\left(\frac{p_{A(x^{*})}}{p_{A(y^{*})}}\right)$

###### Proof.

We essentially follow the proof of Thm. 2.1 in Goldenshluger et al. [7] to
show this result. We begin by noting that $\Phi^{\alpha}_{r}$ is continuous
and concave in $(x,y)\in\mathcal{X}\times\mathcal{X}$ and continuous and
convex in $\phi\in\mathcal{F}$. This follows from the corresponding properties
of $\Phi_{r}$ proved in Ref. [1], noting that
$\Phi^{\alpha}_{r}(x,y;\phi)=\Phi_{r}(x,y;\phi,\alpha)$. Recall that
$\mathcal{X}$ is a compact set and $\mathcal{F}$ is a finite-dimensional real
vector space. Then, it follows from Sion-Kakutani theorem [8] that
$\Phi^{\alpha}_{r}$ has a well-defined saddle-point value, that is

$\inf_{\phi\in\mathcal{F}}\max_{x,y\in\mathcal{X}}\Phi^{\alpha}_{r}(x,y;\phi)=\sup_{x,y\in\mathcal{X}}\inf_{\phi\in\mathcal{F}}\Phi^{\alpha}_{r}(x,y;\phi);$

see Ref. [1] for details. Since $\Phi^{\alpha}_{r}(x,y;\phi)$ is continuous in
$\phi\in\mathcal{F}$ for each $x,y\in\mathcal{X}$, we have that
$\inf_{\phi\in\mathcal{F}}\Phi^{\alpha}(x,y;\phi)$ is upper semi-continuous in
$(x,y)\in\mathcal{X}\times\mathcal{X}$. Then, since $\mathcal{X}$ is compact,
the maximum of $\inf_{\phi\in\mathcal{F}}\Phi^{\alpha}(x,y;\phi)$ is attained
in $\mathcal{X}\times\mathcal{X}$. Therefore, to show the existence of a
saddle-point, it suffices to show that the minimum
$\inf_{\phi\in\mathcal{F}}\Phi^{\alpha}(x,y;\phi)$ is attained in
$\mathcal{F}$ for all $x,y\in\mathcal{X}$.

We show this by proving that $\Phi^{\alpha}_{r}(x,y;\phi)$ is coercive in
$\phi\in\mathcal{F}$ for all $x,y\in\mathcal{X}$, since $\mathcal{F}$ is a
finite-dimensional vector space and $\Phi^{\alpha}_{r}(x,y;\phi)$ is
continuous in $\phi\in\mathcal{X}$ for all $x,y\in\mathcal{X}$. For
convenience, denote $\Phi^{x,y,\alpha}_{r}(\phi)=\Phi^{\alpha}_{r}(x,y;\phi)$.
Recall that $\Phi^{x,y,\alpha}_{r}(\phi)$ is said to be coercive if
$\lim_{\left\lVert\phi\right\rVert\to\infty}\Phi^{x,y,\alpha}_{r}(\phi)=\infty$.
Note that Goldenshluger et al. [7] prove in Thm. (2.1) that the function

$\Theta^{x,y}(\phi)=\ln\left(\int_{\Omega}\exp(-\phi(\omega))p_{A(x)}(\omega)\mathbb{P}(d\omega)\right)+\ln\left(\int_{\Omega}\exp(\phi(\omega))p_{A(y)}(\omega)\mathbb{P}(d\omega)\right)$

is coercive in $\phi\in\mathcal{F}$ for all $x,y\in\mathcal{X}$. Since
$\Phi^{x,y,\alpha}_{r}(\phi)=g^{T}x-g^{T}y+2\alpha
r+\alpha\Theta^{x,y}(\phi/\alpha)$ for fixed $x,y\in\mathcal{X}$ and
$\alpha>0$, it follows that $\Phi^{x,y,\alpha}_{r}(\phi)$ is coercive in
$\phi\in\mathcal{F}$. Therefore, $\Phi^{\alpha}_{r}$ has a saddle-point.

Now, suppose that $(x^{*},y^{*})$ is the $(x,y)$ component of the saddle-point
of $\Phi^{\alpha}_{r}$. Then, if $\phi_{*}$ is the $\phi$-component of the
saddle-point, it minimizes the function $\Phi^{\alpha}_{r}(x^{*},y^{*};\phi)$,
and therefore, $\phi_{*}/\alpha$ minimizes the function
$\Theta^{x^{*},y^{*}}(\phi/\alpha)$. In the proof of Thm (2.1) of
Goldenshluger et al. [7] (Remark (A.1) in particular), it is shown that every
minimum of the function $\Theta^{x,y}(\phi)$ is of the form
$1/2\ln(p_{A(x)}/p_{A(y)})+s$, where $s\in\mathbb{R}$ is a constant. Observe
that $\Phi^{\alpha}_{r}$ is invariant under translations of the form
$\phi\to\phi+s$ for $s\in\mathbb{R}$. Therefore, $\phi_{*}$ can be expressed
as

$\frac{\phi_{*}}{\alpha_{*}}=\frac{1}{2}\ln\left(\frac{p_{A(x^{*})}}{p_{A(y^{*})}}\right)$

by translating the $\phi$-component of the saddle-point with a constant, if
necessary. ∎

Using this result, we express the saddle-point value of $\Phi_{r}$ in terms of
$(x,y)$ and $\alpha$ components.

######  Proposition 1.

The saddle-point value of the function $\Phi_{r}$ defined in Eq. (1) can be
expressed as

$2\Phi_{*}(r)=\inf_{\alpha>0}\left[2\alpha
r+\max_{x,y\in\mathcal{X}}\left(g^{T}x-g^{T}y+2\alpha\ln(\textnormal{AffH}(A(x),A(y)))\right)\right]$

where AffH is the Hellinger affinity defined in Eq. (8). The quantity
$\ln(\textnormal{AffH}(\mu,\nu))$ is well-defined for every
$\mu,\nu\in\mathcal{M}$.

###### Proof.

The function $\Phi_{r}$ can be written as
$\Phi_{r}(x,y;\phi,\alpha)=\Phi^{\alpha}_{r}(x,y;\phi)$, where the function
$\Phi^{\alpha}_{r}$ is defined in Eq. (11). Then, by Lemma (1), we know that
$\Phi^{\alpha}_{r}$ has a well-defined saddle-point for each $\alpha>0$.
Moreover, from the proof of Prop. (3.1) of Ref. [1] (or Thm. (2.1) of Ref.
[7]), we know that

$\inf_{\phi\in\mathcal{F}}\Theta^{x,y}(\phi)=2\ln(\text{AffH}(A(x),A(y)))$

where

$\Theta^{x,y}(\phi)=\ln\left(\int_{\Omega}\exp(-\phi(\omega))p_{A(x)}(\omega)\mathbb{P}(d\omega)\right)+\ln\left(\int_{\Omega}\exp(\phi(\omega))p_{A(y)}(\omega)\mathbb{P}(d\omega)\right).$

Noting that $\Phi_{r}(x,y;\phi)=2\alpha
r+g^{T}x-g^{T}y+\alpha\Theta^{x,y}(\phi/\alpha)$, we have

$\displaystyle 2\Phi_{*}(r)$
$\displaystyle=\inf_{\phi\in\mathcal{F},\alpha>0}\max_{x,y\in\mathcal{X}}\Phi_{r}(x,y;\phi,\alpha)$
$\displaystyle=\inf_{\alpha>0}\left[2\alpha
r+\inf_{\phi\in\mathcal{F}}\max_{x,y\in\mathcal{X}}\left(g^{T}x-g^{T}y+\alpha\Theta^{x,y}(\phi/\alpha)\right)\right]$
$\displaystyle=\inf_{\alpha>0}\left[2\alpha
r+\max_{x,y\in\mathcal{X}}\left(g^{T}x-g^{T}y+\alpha\inf_{\phi^{\prime}\equiv\phi/\alpha\in\mathcal{F}}\Theta^{x,y}(\phi^{\prime})\right)\right]$
$\displaystyle=\inf_{\alpha>0}\left[2\alpha
r+\max_{x,y\in\mathcal{X}}\left(g^{T}(x-y)+2\alpha\ln\left(\textnormal{AffH}(A(x),A(y))\right)\right)\right]$

giving the desired result. We note that $\ln(\textnormal{AffH}(A(x),A(y)))$ is
continuous and concave on $\mathcal{X}\times\mathcal{X}$, as shown in Prop.
(3.1) of Ref. [1]. Observe that the maximum over $x,y\in\mathcal{X}$ is always
non-negative because we obtain a value of zero when $x=y$.

Note that $\textnormal{AffH}(\mu,\nu)>0$ for all $\mu,\nu\in\mathcal{M}$. This
is because
$\textnormal{AffH}(\mu,\nu)=\int_{\Omega}\sqrt{p_{\mu}(\omega)p_{\nu}(\omega)}\mathbb{P}(d\omega)$,
and since $p_{\mu}(\omega)>0$ for every $\omega\in\Omega$ and
$\mu\in\mathcal{M}$ by definition of a good pair, the integrand is positive.
Consequently, $\ln(\textnormal{AffH}(\mu,\nu))$, and therefore
$\ln(\textnormal{AffH}(A(x),A(y)))$, is well-defined. ∎

Therefore, we can use the above result to compute the $(x,y)$ and $\alpha$
components of the saddle-point of $\Phi_{r}$ to any precision $\delta>0$. We
can then construct an estimator as described in Eq. (10). In the following
proposition, we show that the $\epsilon$-risk of this estimator is bounded
above by $\Phi_{*}(\ln(2/\epsilon))+\delta$.

######  Proposition 2.

Let $1-\epsilon\in(0,1)$ specify the confidence level, and let $\delta>0$ be
the given precision. Let $\alpha_{*}>0$ be chosen such that

$2\Phi_{*}(\ln(2/\epsilon))+\delta\geq
2\alpha_{*}\ln(2/\epsilon)+\max_{x,y\in\mathcal{X}}\left(g^{T}x-g^{T}y+2\alpha_{*}\ln(\textnormal{AffH}(A(x),A(y)))\right)$

Let $x^{*},y^{*}\in\mathcal{X}$ be any points that attain the maximum in
$\max_{x,y\in\mathcal{X}}\left(g^{T}x-g^{T}y+2\alpha_{*}\ln(\textnormal{AffH}(A(x),A(y)))\right)$,
and correspondingly, define

$\phi_{*}=\frac{\alpha_{*}}{2}\ln\left(\frac{p_{A(x^{*})}}{p_{A(y^{*})}}\right)$

Then the estimator

$\widehat{g}_{*}=\phi_{*}+\frac{1}{2}\left(g^{T}x^{*}+g^{T}y^{*}\right)$

has $\epsilon$-risk bounded above by $\Phi_{*}(\ln(2/\epsilon))+\delta$.

###### Proof.

Since $x^{*},y^{*}\in\mathcal{X}$ correspond to the maximum, we can write

$\max_{x,y\in\mathcal{X}}\left(g^{T}x-g^{T}y+2\alpha_{*}\ln(\textnormal{AffH}(A(x),A(y)))\right)=g^{T}x^{*}-g^{T}y^{*}+2\alpha_{*}\ln(\textnormal{AffH}(A(x^{*}),A(y^{*})))$

For the choice

$\phi_{*}=\frac{\alpha_{*}}{2}\ln\left(\frac{p_{A(x^{*})}}{p_{A(y^{*})}}\right)$

we have

$2\ln(\textnormal{AffH}\left(A(x^{*}),A(y^{*})\right)=\ln\left(\int_{\Omega}\exp(-\phi_{*}(\omega)/\alpha_{*})p_{A(x^{*})}(\omega)\mathbb{P}(d\omega)\right)+\ln\left(\int_{\Omega}\exp(\phi_{*}(\omega)/\alpha_{*})p_{A(y^{*})}(\omega)\mathbb{P}(d\omega)\right)$

Therefore, we can write

$2\alpha_{*}\ln(2/\epsilon)+\max_{x,y\in\mathcal{X}}\left(g^{T}x-g^{T}y+2\alpha_{*}\ln(\textnormal{AffH}(A(x),A(y)))\right)=\Phi^{\alpha_{*}}_{\ln(2/\epsilon)}(x^{*},y^{*};\phi_{*})$

where $\Phi^{\alpha_{*}}_{\ln(2/\epsilon)}$ is obtained using Eq. (11). Now,
observe that we have

$\max_{x,y\in\mathcal{X}}\inf_{\phi\in\mathcal{F}}\Phi^{\alpha_{*}}_{\ln(2/\epsilon)}(x,y;\phi)=2\alpha_{*}\ln(2/\epsilon)+\max_{x,y\in\mathcal{X}}\left(g^{T}x-g^{T}y+2\alpha_{*}\ln(\textnormal{AffH}(A(x),A(y)))\right)$

which follows from the proof of Prop. (1). Therefore, $(x^{*},y^{*})$ and
$\phi_{*}$ correspond to $(x,y)$ and $\phi$ components of the saddle-point of
$\Phi^{\alpha_{*}}_{\ln(2/\epsilon)}$, respectively. Consequently, the points
$x^{*},y^{*}\in\mathcal{X}$ achieve the maximum in
$\max_{x,y\in\mathcal{X}}\Phi^{\alpha_{*}}_{\ln(2/\epsilon)}(x,y;\phi_{*})$.
In particular, we have

$\displaystyle\Phi^{\alpha_{*}}(x,y^{*};\phi_{*})\leq\Phi^{\alpha^{*}}_{\ln(2/\epsilon)}(x^{*},y^{*};\phi_{*})\leq
2\Phi_{*}(\ln(2/\epsilon))+\delta\quad(\forall x\in\mathcal{X})$
$\displaystyle\Phi^{\alpha_{*}}(x^{*},y;\phi_{*})\leq\Phi^{\alpha^{*}}_{\ln(2/\epsilon)}(x^{*},y^{*};\phi_{*})\leq
2\Phi_{*}(\ln(2/\epsilon))+\delta\quad(\forall y\in\mathcal{X})$ (12)

where the last inequality follows from the definition of $\alpha_{*}$.

Next, we rewrite the expression for the constant in the estimator
$\widehat{g}_{*}$ in a convenient form. Since

$\int_{\Omega}\exp(-\phi_{*}(\omega)/\alpha_{*})p_{A(x^{*})}(\omega)\mathbb{P}(d\omega)=\int_{\Omega}\exp(\phi_{*}(\omega)/\alpha_{*})p_{A(y^{*})}(\omega)\mathbb{P}(d\omega)$

holds for our choice of $\phi_{*}$, we can write

$\displaystyle c$
$\displaystyle\equiv\frac{1}{2}(g^{T}x^{*}+g^{T}y^{*})=\frac{1}{2}\left[g^{T}x^{*}+\alpha_{*}\ln\left(\int_{\Omega}\exp(-\phi_{*}(\omega)/\alpha_{*})p_{A(x^{*})}(\omega)P(d\omega)\right)+\alpha_{*}\ln(2/\epsilon)\right]$
$\displaystyle\hskip
85.35826pt-\frac{1}{2}\left[-g^{T}y^{*}+\alpha_{*}\ln\left(\int_{\Omega}\exp(\phi_{*}(\omega)/\alpha_{*})p_{A(y^{*})}(\omega)P(d\omega)\right)+\alpha_{*}\ln(2/\epsilon)\right]$
$\displaystyle=\frac{1}{2}\Phi^{\alpha_{*}}_{\ln(2/\epsilon)}(x^{*},y^{*};\phi_{*})-\left[-g^{T}y^{*}+\alpha_{*}\ln\left(\int_{\Omega}\exp(\phi_{*}(\omega)/\alpha_{*})p_{A(y^{*})}(\omega)P(d\omega)\right)+\alpha_{*}\ln(2/\epsilon)\right]$
(13)
$\displaystyle=\left[g^{T}x^{*}+\alpha_{*}\ln\left(\int_{\Omega}\exp(-\phi_{*}(\omega)/\alpha_{*})p_{A(x^{*})}(\omega)P(d\omega)\right)+\alpha_{*}\ln(2/\epsilon)\right]-\frac{1}{2}\Phi^{\alpha_{*}}_{\ln(2/\epsilon)}(x^{*},y^{*};\phi_{*})$
(14)

Note that our estimator is given as
$\widehat{g}_{*}=\phi_{*}+c\in\mathcal{F}$. We define the quantity
$R=\Phi_{*}(\ln(2/\epsilon))+\delta$. Then, for any $x\in X$, we have

$\displaystyle
g^{T}x+\alpha_{*}\ln\left(\int_{\Omega}\exp(-\widehat{g}_{*}/\alpha_{*})p_{A(x)}(\omega)P(d\omega)\right)+\alpha_{*}\ln(2/\epsilon)$
$\displaystyle\qquad=g^{T}x+\alpha_{*}\ln\left(\int_{\Omega}\exp(-\phi_{*}(\omega)/\alpha_{*})p_{A(x)}(\omega)P(d\omega)\right)+\alpha_{*}\ln(2/\epsilon)-c$
$\displaystyle\qquad=\Phi^{\alpha_{*}}_{\ln(2/\epsilon)}(x,y^{*};\phi_{*})-\frac{1}{2}\Phi^{\alpha_{*}}_{\ln(2/\epsilon)}(x^{*},y^{*};\phi_{*})$
$\displaystyle\qquad\leq\Phi_{*}(\ln(2/\epsilon))+\frac{\delta}{2}$
$\displaystyle\qquad=R-\frac{\delta}{2}$

where we used Eq. (13) and Eq. (12). Similarly, using Eq. (14) and Eq. (12),
we find that

$-g^{T}y+\alpha_{*}\ln\left(\int_{\Omega}\exp(\widehat{g}_{*}/\alpha_{*})p_{A(y)}(\omega)P(d\omega)\right)+\alpha_{*}\ln(2/\epsilon)\leq
R-\frac{\delta}{2}$

for all $y\in\mathcal{X}$. Then, dividing by $\alpha_{*}>0$ and rearranging
the terms, we find that

$\displaystyle\ln\left(\int_{\Omega}\exp((g^{T}x-\widehat{g}_{*}-R)/\alpha_{*})p_{A(x)}(\omega)P(d\omega)\right)$
$\displaystyle\leq\ln\left(\frac{\epsilon}{2}\right)-\frac{\delta}{2\alpha_{*}}\equiv\ln\left(\frac{\epsilon^{\prime}}{2}\right)\quad(\forall
x\in\mathcal{X})$
$\displaystyle\ln\left(\int_{\Omega}\exp((-g^{T}y+\widehat{g}_{*}-R)/\alpha_{*})p_{A(y)}(\omega)P(d\omega)\right)$
$\displaystyle\leq\ln\left(\frac{\epsilon}{2}\right)-\frac{\delta}{2\alpha_{*}}\equiv\ln\left(\frac{\epsilon^{\prime}}{2}\right)\quad(\forall
y\in\mathcal{X})$

where $\epsilon^{\prime}=\epsilon e^{-\delta/2\alpha_{*}}<\epsilon$. Then,
from Markov’s inequality, it follows that

$\displaystyle\operatorname*{Prob}_{\omega\sim
p_{A(x)}}\left\\{g^{T}x-\widehat{g}_{*}-R\geq
0\right\\}\leq\mathbb{E}[\exp((g^{T}x-\widehat{g}_{*}-R)/\alpha_{*})]$
$\displaystyle\leq\frac{\epsilon^{\prime}}{2}\quad(\forall x\in\mathcal{X})$
$\displaystyle\operatorname*{Prob}_{\omega\sim
p_{A(y)}}\left\\{-g^{T}y+\widehat{g}_{*}-R\geq
0\right\\}\leq\mathbb{E}[\exp((-g^{T}y+\widehat{g}_{*}-R)/\alpha_{*})]$
$\displaystyle\leq\frac{\epsilon^{\prime}}{2}\quad(\forall y\in\mathcal{X}).$

Therefore, using the union bound, we can conclude that

$\operatorname*{Prob}_{\omega\sim
p_{A(x)}}\left\\{|\widehat{g}_{*}(\omega)-g^{T}x|\geq
R\right\\}\leq\epsilon^{\prime}<\epsilon\quad(\forall x\in\mathcal{X}).$

Therefore, the $\epsilon$-risk of the estimator $\widehat{g}_{*}$ is bounded
above by $R=\Phi_{*}(\ln(2/\epsilon))+\delta$. ∎

## References

  * [1] Anatoli Juditsky and Arkadi Nemirovski.  Nonparametric estimation by convex programming.  The Annals of Statistics, 37(5A):2278–2300, 2009.
  * [2] D. Donoho.  Statistical estimation and optimal recovery.  Annals of Statistics, 22(1):238–270, 1994.
  * [3] M. Grant and S. Boyd.  CVX: Matlab Software for Disciplined Convex Programming, version 2.2, 2020.  Available at http://cvxr.com/cvx.
  * [4] Anatoli Juditsky and Arkadi Nemirovski.  On well-structured convex-concave saddle point problems and variational inequalities with monotone operators.  arXiv preprint arXiv:2102.01002, 2021.
  * [5] Anatoli Juditsky and Arkadi Nemirovski.  Near-optimal recovery of linear and n-convex functions on unions of convex sets.  Information and Inference: A Journal of the IMA, 9(2):423–453, 2020\.
  * [6] Anatoli Juditsky and Arkadi Nemirovski.  Statistical Inference via Convex Optimization.  Princeton University Press, 2020.
  * [7] Alexander Goldenshluger, Anatoli Juditsky, and Arkadi Nemirovski.  Hypothesis testing by convex optimization.  Electronic journal of statistics, 9(2):1645–1712, 2015.
  * [8] Maurice Sion.  On general minimax theorems.  Pacific Journal of mathematics, 8(1):171–176, 1958.


abstract=false, titlepage=false, twocolumn=true, title subtitle disposition
date author publishers captionlabel USenglishFigureFig. *0.5em

# Elastic shakedown and roughness evolution in repeated elastic-plastic
contact

Lucas Frérot<EMAIL_ADDRESS>https://orcid.org/0000-0002-4138-1052 Lars Pastewka
Department of Microsystems Engineering, University of Freiburg,
Georges-Köhler-Allee 103, 79110 Freiburg, Germany
https://orcid.org/0000-0001-8351-7336

Surface roughness emerges naturally during mechanical removal of material,
fracture, chemical deposition, plastic deformation, indentation, and other
processes. Here, we use continuum simulations to show how roughness which is
neither Gaussian nor self-affine emerges from repeated elastic-plastic contact
of rough and rigid surfaces on a flat elastic-plastic substrate. Roughness
profiles change with each contact cycle, but appear to approach a steady-state
long before the substrate stops deforming plastically and has hence “shaken-
down” elastically. We propose a simple dynamic collapse for the emerging
power-spectral density, which shows that the multi-scale nature of the
roughness is encoded in the first few indentations. In contrast to macroscopic
roughness parameters, roughness at small scales and the skewness of the height
distribution of the resulting roughness do not show a steady-state, with the
latter vanishing asymptotically with contact cycle.

###### keywords:

surface; roughness; contact; elastoplasticity; tribology; elastic shakedown

##  1 Introduction

Friction, wear, thermal transfer and other tribological phenomena depend on
the roughness of surfaces in contact: Its presence means that contact only
occurs on a limited subset of the apparent area of contact, the _true contact
area_ $A$, whose magnitude and geometry are of prime importance for tribology
[1]. Understanding the origins of roughness, with its multiscale nature [2,
3], and making predictions on the true contact area are therefore objectives
of the study of rough interfaces between solids. They are also both intimately
linked to material properties.

Early models for contact of rough surfaces assumed that deformation is purely
plastic [4, 5, 6, 7, 8]. Assuming a surface hardness $p_{\mathrm{m}}$, this
means that for a normal force $F_{\mathrm{n}}$ the true area of contact is
given by $A=F_{\mathrm{n}}/p_{\mathrm{m}}$, because the mean pressure within
contacts equals $p_{\mathrm{m}}$. However, researchers in the 1950s and 60s
discussed the idea of an “elastic shakedown” [9]: It seemed unlikely, that a
contact would continue to behave plastically during repeated loading.
Eventually, the contact should “shake down” to a point where it only responds
elastically. This eventually lead to the development of elastic contact
models, such as the famous 1966 model of Greenwood and Williamson [10]. The
idea of elastic shakedown is probably one of the oldest ideas in the rough
contact literature. Here, we present numerical calculations of the contact of
rough surfaces with an elastic-plastic substrate to gain insights into the
deformation processes during shakedown.

It seems reasonable to assume that a signature of shakedown is hidden in
surface roughness: Self-affine surface roughness is typically associated with
irreversible deformation, such as fracture [11, 12] or wear [2, 13, 14, 15,
16, 17], but plastic deformation can also contribute to the emergence of
roughness [18, 19, 20]. In this work, we focus on the roughness created
through quasi-static repeated indentation (i.e. without sliding), at a small
contact area ratio, of an initially flat elastic-plastic half-space with a
number of rough, rigid, periodic, isotropic counter-faces, each different. Our
work stands apart from prior theoretical work on elastic-plastic contact [21,
22], which has investigated the contact of a rough deformable half-space with
a rigid but flat plane. While those calculations show smoothening of
asperities, we find that the initially flat counterbody continuously roughens
during our calculations. However, the statistical properties of the deformed
body differ significantly from the indenting surface: Profiles are non-
Gaussian and do not appear to be self-affine, even if the indenting surface
has ideal self-affine properties. Even though simple macroscopic properties,
such as the RMS height, appear to saturate at a steady-state value after a few
contact cycles, our calculations reveal that surface topography and subsurface
plastic deformation continues to evolve beyond this apparent “shake down” of
the contact.

##  2 Methods

###  2.1 Elastoplastic Contact

Figure 1: Two dimensional schematic view of the repeated indentation procedure
we employ in this work. Indentation is made with a different surface at each
iteration to reflect that in practice, contacts would be unlikely to occur
twice with the same alignment of surfaces. $\mathcal{B}$ is a semi-infinite
elastic-plastic solid obeying a $J_{2}$ isotropic linear hardening plasticity
model. The rigid, rough, isotropic counterface has a period $L$.

Figure 1 illustrates the repeated indentation procedure described in the
introduction. We simulate the contact between an elastic-plastic half-space
$\mathcal{B}$ and a periodic rigid rough surface using a volume-integral
equation approach [23, 24, 25]. In these repeated contact simulations, without
horizontal movement, we denote $h_{k}$ the rigid surface of the $k$-th
indentation. Unlike conventional boundary-integral equation methods used in
elastic and elastic-plastic rough contacts [26, 27, 28, 21, 29], our approach
fully considers plastic residual strains in the subsurface region of the half-
space. The subsurface (volumetric) plastic strain field contributes to
displacements and stress [30]. Equilibrium equations are solved with Green’s
functions applied in the Fourier domain [25].

The plasticity model assumes an additive strain decomposition and a von Mises
($J_{2}$) yield criterion, $f_{\mathrm{y}}$:

$\displaystyle{\bm{\sigma}}$
$\displaystyle=\mathcal{C}:({\bm{\varepsilon}}-{\bm{\varepsilon}}^{p}),$ (1)
$\displaystyle f_{\mathrm{y}}({\bm{\sigma}})$
$\displaystyle=\sqrt{\frac{3}{2}}\|{\bm{s}}\|,\quad\text{with
}{\bm{s}}={\bm{\sigma}}-\frac{1}{3}\mathrm{tr}({\bm{\sigma}}){\bm{I}}.$ (2)

We have defined ${\bm{\sigma}}$ the Cauchy stress tensor, ${\bm{\varepsilon}}$
the total strain, $\mathcal{C}$ the small-strain elasticity tensor,
${\bm{\varepsilon}}^{p}$ the plastic strain and ${\bm{s}}$ the deviatoric
stress. The colon operator here introduces contraction on the two right-most
indices of a tensor. This type of model, also referred to as $J_{2}$
plasticity [31], is the most widely-used canonical plasticity model. It
assumes plastic deformation progresses as volume-conserving laminar flow
beyond a yield shear (von Mises) stress, but not under hydrostatic conditions.
The model is empirical and material-unspecific, but captures the basic
phenomenology of plastic flow at macroscopic scales.

Since our simulations are quasi-static, and the pseudo-time variable is given
by the number of different surfaces used in indentation, $k$, we describe
plastic strain rates as increments [31]. The equivalent cumulated plastic
strain is defined as:

$\displaystyle
e^{\mathrm{p}}=\sqrt{\frac{2}{3}}\sum_{i=1}^{k}{\|\Delta{{\bm{\varepsilon}}}^{\mathrm{p}}_{i}\|}.$
(3)

The plasticity conditions can be written as, in terms of $e^{\mathrm{p}}$ and
$\sigma$:

$\displaystyle f_{\mathrm{y}}({\bm{\sigma}})-f_{\mathrm{h}}(e^{\mathrm{p}})$
$\displaystyle\leq 0,$ (4) $\displaystyle\Delta{e}^{\mathrm{p}}$
$\displaystyle\geq 0,$ (5)
$\displaystyle\left(f_{\mathrm{y}}({\bm{\sigma}})-f_{\mathrm{h}}(e^{\mathrm{p}})\right)\Delta{e}^{\mathrm{p}}$
$\displaystyle=0,$ (6)

where $f_{\mathrm{h}}(e^{\mathrm{p}})=\sigma_{0}+E_{\mathrm{h}}e^{\mathrm{p}}$
is the linear hardening function with initial yield stress $\sigma_{0}$ and
hardening modulus $E_{\mathrm{h}}$. These conditions express that the von
Mises stress should never exceed the yield stress and that plastic deformation
is irreversible. We assume an associated flow rule, so the plastic strain
increment is expressed as:

$\displaystyle\Delta{\bm{\varepsilon}}^{\mathrm{p}}=\frac{3\Delta{e}^{\mathrm{p}}}{2f_{\mathrm{y}}({\bm{\sigma}}+\mathcal{C}:\Delta{\bm{\varepsilon}})}{\bm{s}}({\bm{\sigma}}+\mathcal{C}:\Delta{\bm{\varepsilon}}).$
(7)

The usual return-mapping algorithm [31] is used to compute $\Delta
e^{\mathrm{p}}$. Solving the coupled elastic-plastic contact problem is done
with a fixed-point iterative approach [24, 25] accelerated with an Anderson
mixing procedure [32, 33]. The code, available as a supplementary material, is
built on top of the open-source high-performance contact library Tamaas [34].

For reference, we also used the surface plasticity approach of [21] and a
rigid-plastic, as known as a bearing area model. The former solves an elastic
contact problem with the constraint that the normal surface pressure be less
than or equal to the indentation hardness $p_{\mathrm{m}}$ [35], correcting
the counterface so that the constraint is satisfied. We use this correction as
the residual plastic displacement [21, 29]. The choice of the hardness
$p_{\mathrm{m}}$ is a debated question: While $p_{\mathrm{m}}\approx
3\sigma_{\mathrm{y}}$ is a common choice stemming from spherical indentation
[35], finite element simulations of sinusoidal [36, 37] and rough surfaces
[38] report values in the order of $6\sigma_{\mathrm{y}}$, and a quantitative
agreement with the $J_{2}$ model requires a surface-specific adjustment of
$p_{m}$ [30]. Here we nonetheless set
$p_{\mathrm{m}}={{3}}\sigma_{\mathrm{y}}$. The model is hereafter referred to
as the saturated plasticity model. Finally, the rigid-plastic, or bearing area
approach, consists in solving the contact entirely without elastic
deformation, by finding a plane intersecting the rough counterface such that
the intersection area equals the normal force, $F_{\mathrm{n}}$, divided by
$p_{\mathrm{m}}$. Note that these two models corresponds roughly to historic
plastic contact models by [8].

###  2.2 Roughness Statistics

Gaussian, self-affine rough surfaces with a period $L$, are generated using a
random-phase algorithm [39, 40, 3, 41] for each indentation step. They follow
the isotropic power-spectral density (PSD)

$\displaystyle\left|\mathcal{F}[h_{k}]\right|^{2}$
$\displaystyle=\begin{cases}C&\quad\text{if }q_{0}\leq q\leq q_{1},\\\
C\left(\frac{q}{q_{1}}\right)^{-2(H+1)}&\quad\text{if }q_{1}\leq q\leq
q_{2},\\\ 0&\quad\text{otherwise,}\end{cases}$ (8)

where $\mathcal{F}$ is the Fourier transform, $q_{i}$ are angular wavenumbers
for the cutoffs in the surface spectrum, $H$ is the Hurst exponent, and $C$ an
unspecified constant, adjusted such that the root-mean-square of surface
slopes has the desired expected value.

The residual displacement resulting from plastic deformation at step $k$ is
noted ${\bm{u}}^{\mathrm{p}}_{k}$ and is computed using a Green’s function
[42], which accounts for “spring-back” effects,

$\displaystyle{\bm{u}}^{\mathrm{p}}_{k}(x)=\int_{\mathcal{B}}{{\bm{\nabla}}_{y}{\bm{U}}(x,y):\mathcal{C}:{\bm{\varepsilon}}^{\mathrm{p}}_{k}\,\mathrm{d}V_{y}},$
(9)

where ${\bm{U}}$ is the Mindlin tensor [43, 44, 25]. We note
$h^{\mathrm{p}}_{k}$ the trace on $\partial\mathcal{B}$, the boundary of
$\mathcal{B}$, of the component of ${\bm{u}}^{\mathrm{p}}_{k}$ normal to the
surface: It is the emerging surface roughness due to successive indentations.

Be it from plastic indentation or from synthetic generation, we derive from a
rough surface $h$ the following statistical properties [41],

$\displaystyle\left(h^{(\alpha)}_{\mathrm{rms}}\right)^{2}$
$\displaystyle=\frac{1}{{{L^{2}}}}\int_{{{A_{L}}}}{|\nabla^{\alpha}h|^{2}\,\mathrm{d}A},$
(10) $\displaystyle h_{\mathrm{rms}}$
$\displaystyle=h^{(0)}_{\mathrm{rms}},\quad
h^{\prime}_{\mathrm{rms}}=h^{(1)}_{\mathrm{rms}},\quad
h^{\prime\prime}_{\mathrm{rms}}=\frac{1}{2}h^{(2)}_{\mathrm{rms}},$ (11)
$\displaystyle s$
$\displaystyle=\frac{1}{h_{\mathrm{rms}}^{3}{{L^{2}}}}\int_{{{A_{L}}}}{(h-\langle
h\rangle)^{3}\,\mathrm{d}A},$ (12)

where $A_{L}=[0,L]^{2}$, $\nabla^{\alpha}$ is the fractional Riesz derivative,
computed in the Fourier domain, which links $h^{(\alpha)}_{\mathrm{rms}}$ to
Nayak’s definition of PSD moments [45]. The parameters $h_{\mathrm{rms}}$,
$h^{\prime}_{\mathrm{rms}}$ and $h^{\prime\prime}_{\mathrm{rms}}$ are the
root-mean-square (RMS) of heights, RMS of slopes and RMS of curvatures,
respectively. The parameter $s$ is the skewness of the surface height
distribution.

Rather than reporting the curvature directly, we report it in form of Nayak’s
parameter [45],

$\alpha=4(h_{\mathrm{rms}}h_{\mathrm{rms}}^{\prime\prime}/h_{\mathrm{rms}}^{\prime
2})^{2}.$ (13)

The advantage of using Nayak’s parameter is that it removes the influence of
changes in roughness amplitude.

Since rough surfaces are multi-scale objects, we use a multi-scale
quantitative roughness analysis in addition to the scalar measures defined
above. As one of the most common statistical descriptors, we report the
isotropic PSD $\phi(q)=|\mathcal{F}[h]|^{2}$. We also apply a variable
bandwidth method (VBM) in form of detrended fluctuation analysis [46, 47, 48,
41]. Specifically, the surface with side length $L$ is divided into square
patches of side length $\ell$ (see inset in fig. 7 for an illustration). Each
patch is then independently detrended by subtracting the plane that minimizes
$h_{\mathrm{rms}}$ in that patch. Evaluating the above roughness parameters on
the detrended patches for decreasing $\ell$ gives length-scale dependent
metrics. Both PSD and VBM allow the study of a surface’s potential self-
affinity, with self-affine surfaces following $\phi(q)\propto q^{-2-2H}$ or
$h_{\mathrm{rms}}(\ell)\propto\ell^{H}$ over the self-affine bandwidth.

We also compute the dissipated plastic energy $d^{\mathrm{p}}_{k}$ during
successive indentations, given by

$d^{\mathrm{p}}_{k}=\sum_{i=1}^{k}\int_{{{\mathcal{B}_{L}}}}{\left({\bm{\sigma}}_{i}:\Delta{{\bm{\varepsilon}}}^{\mathrm{p}}_{i}-E_{\mathrm{h}}e^{\mathrm{p}}_{i}\Delta
e^{\mathrm{p}}_{i}\right)\mathrm{d}V},$ (14)

for the $J_{2}$ model ($\mathcal{B}_{L}=A_{L}\times[0,+\infty[$), and the
irrecoverable energy, which includes the dissipated plastic energy and the
stored elastic energy due to residual stresses, and can be computed directly
on the surface for both models:

$\delta^{\mathrm{p}}_{k}=\sum_{i=1}^{k}\int_{{{A_{L}}}}{p_{i}\Delta{h}^{\mathrm{p}}_{i}\mathrm{d}A}.\\\
$ (15)

Note that the saturated model lacks a description of eigenstresses [30], thus
$\delta^{\mathrm{p}}_{k}\equiv d^{\mathrm{p}}_{k}$ is the dissipated plastic
energy for this model.

##  3 Results

We simulate the repeated indentation of 90 different surfaces with the same
statistical properties and the following PSD parameters: $q_{0}L/2\pi=2$,
$q_{1}L/2\pi=2$, $q_{2}L/2\pi=128$, $H\in\\{0.5,0.8\\}$,
$h^{\prime}_{\mathrm{rms}}=0.4$, where $L$ is the horizontal period of the
surfaces. Since the computational cost of the saturated model is significantly
lower than the $J_{2}$ model, we simulate it in 270 repeated indentation
steps. Contact is made against a half-space with material properties
$\nu={{0.3}}$, $\sigma_{0}=0.1E$, $E_{h}=0.1E$, with $E$ the Young modulus.
All properties are nondimensionalized by $E$ and by the horizontal period $L$
and only those parameters control the outcome of the simulations.

Figure 2 shows a rendered view of the emerging surface roughness due to the
accumulation of plastic residual deformation, for a normal force of
$F_{\mathrm{n}}/(Eh^{\prime}_{\mathrm{rms}}{{L^{2}}})=0.025$ (which
corresponds to a contact area ratio of 7%), and a discretization size of
729$\times$729$\times$55 points for a domain size of $L\times L\times 0.075L$
(total of 175M degrees-of-freedom). The highlighted areas in fig. 2 show where
plastic energy is dissipated on the surface for the current indentation step.
We also show select one-dimensional cross-sectional profiles of these
topography maps below them. After enough steps, no remaining “undeformed” area
of the surface can be identified by eye, and the surface appears entirely
rough from plasticity.

Figure 2: Evolution of surface roughness through repeated elastic-plastic
indentation. Highlighted areas on the three-dimensional views show where
plastic energy is dissipated at the specific step that is shown. The panels
show the topography after (a) $k=1$, (b) $10$ and (c) $90$ repeated contacts.
The bottom row shows the evolution of one-dimensional cross-sectional profiles
of the above surfaces. After sufficient number of iterations, no undeformed
(flat) areas can be identified in the surface. Figure 3: Evolution of the
normalized root-mean-square (RMS) of (a) heights, (b) slopes, (c) dissipated
and irrecoverable energy, and (d) Nayak’s parameter with indentation step.
While, for the $J_{2}$ model, the RMS of heights and slopes show similar
evolution, with a steady-state after $\sim 40$ contact cycles, Nayak’s
parameter converges earlier and the dissipated energy does not reach a steady-
state value within $90$ contacts. The scaling of the saturated results with
$k$ is similar to the $J_{2}$ model, with a steady state reached after a
larger number of iterations. Nayak’s parameter (d) differs qualitatively
between calculations with $J_{2}$ and saturated plasticity models due to slope
discontinuities in the residual displacement.

We first evaluate macroscopic scalar roughness parameters on the emerging
roughness. Figure 3 shows the evolution, for the indented roughness, of the
root-mean-square of heights $h^{\mathrm{p}}_{k,\mathrm{rms}}$ (panel a), the
RMS of slopes $h^{\mathrm{p}\prime}_{k,\mathrm{rms}}$ (b), the dissipated
energy $d^{\mathrm{p}}_{k}$ and the irrecoverable energy
$\delta^{\mathrm{p}}_{k}$ normalized by the elastic conforming energy of the
counterface $w^{e}=[h^{(1/2)}_{\mathrm{rms}}]^{2}E/(4(1-\nu^{2}))$ (c), and
Nayak’s parameter, $\alpha^{\mathrm{p}}$ (d). The dynamics of the first three
parameters is similar in the two plasticity models, $J_{2}$ (von Mises), and
saturated. Both RMS metrics initially scale as a power-law with an exponent
close to $1/2$ before plateauing at steady-state values, with the $J_{2}$
model transitioning earlier to lower values than the saturated model.
Interestingly, the saturated model, which does not have hardening, also
reaches a plateau. The overall scaling of RMS metrics is independent of
$p_{\mathrm{m}}$ (not shown), and $p_{\mathrm{m}}$ only controls the contact
area at each indentation: A lower contact ratio delays the appearance of a
steady-state roughness. While energy dissipation in the saturated models
appears to be constant per indentation step, the $J_{2}$ model shows lower
energy dissipation rates after the first $\sim 20$ contact cycles. The
irrecoverable energy shows a trend similar to the dissipated energy. It has a
larger value, which indicates a build-up of residual stresses with successive
indentations. Residual stresses require work for their formation, and while
the stored elastic energy is in principle reversible, it is not recoverable in
practice. The ratio $\delta^{\mathrm{p}}_{k}/d^{\mathrm{p}}_{k}$ increases
from 1.22 to 1.34 with $k$, reflecting an increase in the proportion of
elastic stored energy in $\delta^{\mathrm{p}}_{k}$, consistent with hardening
in the $J_{2}$ model.

As a geometric measure, Nayak’s parameter evolves differently in the $J_{2}$
and the saturated plasticity model. In the saturated plasticity model,
residual displacements are non-zero only in areas where the contact pressure
$p=p_{\mathrm{m}}$. This criterion creates a slope discontinuity at the edge
of zones where $p=p_{\mathrm{m}}$ which dominates the high-frequency part of
the PSD. They are reflected in large values of
$h^{\prime\prime}_{\mathrm{rms}}$, contributing to the sustained increase in
$\alpha_{k}^{\mathrm{p}}$ in the saturated model.

The influence of the Hurst exponent $H$ disappears when we nondimensionalize
lengths, slopes and energies by RMS heights, slopes and conforming elastic
energy, respectively, of the counterface. However, normalization does not
collapse the steady-state value of RMS heights and Nayak’s parameter,
indicating that details of the counterface self-affine structure non-trivially
affect the long-wavelength part of the emerging roughness spectrum.

Figure 4: Evolution of the probability density function (PDF) of residual
roughness heights for (a) the $J_{2}$, (b) the saturated plasticity model and
(c) the rigid-plastic model. The peak corresponds to the residual deformation
in untouched areas. As more areas come in contact, the peak broadens and
eventually disappears. The height distribution then tends towards a non-
Gaussian, asymmetric form, which resembles common extreme value distributions.

Having established the evolution of macroscopic parameters with contact cycle,
we now detail the evolution of the distribution of heights. Figure 4 shows,
for the $J_{2}$ (a), saturated (b), and rigid-plastic (c) models, the height
probability density function of the emerging roughness through the indentation
steps, for the counterface with $H=0.8$.

While the distributions are always non-Gaussian, there is a clear evolution in
the shape of the distribution with $k$: For low indentation count, the
distribution obtained from the $J_{2}$ model is dominated by a well-defined
peak, similar to a rectified distribution, with an exponential tail. This peak
occurs for a height which corresponds to the residual displacement in the
zones that have not been in contact yet. The peak has a finite width due to
the non-local influence of subsurface residual strains. As more surfaces come
in contact with increasing $k$, the “undeformed” area of the surface shrinks:
the peak height decreases and its breadth increases due to a larger variation
of residual displacements, caused by ever-increasing, in number and size,
plastic zones. The peak eventually disappears, and the height distribution
settles on a more regular shape.

Unlike the $J_{2}$ model, the saturated model gives a true rectified
distribution, where the peak in fig. 4a becomes a Dirac distribution in fig.
4b, whose magnitude is given by the remaining (truly) undeformed area. The
tail of the distribution at large (positive) heights is similar for both
models, although the saturated model has not reached a steady-state RMS of
heights at the last iteration shown in fig. 4b. Similarly to the saturated,
the rigid-plastic model (fig. 4c) also shows a rectified distribution, except
that the residual height range is larger by a factor of approximately 2. This
is due to the initial elastic deformation needed to overcome the plastic yield
stress, absent in the rigid-plastic approach.

Figure 4 shows fits of the simulation data with two extreme value
distributions: the Weibull111cumulative distribution
$\mathrm{cdf}(x)=1-\exp((-x/\lambda)^{\kappa})$, $\lambda=0.3$, $\kappa=1.9$
distribution, which [49] used in an elastic contact analysis of non-Gaussian
surfaces, and the Gumbel distribution222cumulative distribution
$\mathrm{cdf}(x)=\exp(-\exp(-x/\beta))$, $\beta=1/9$. We use extreme value
distributions because, owing to the relatively low contact area and to long-
range elastic interactions—which makes contact near large asperities unlikely,
only the largest asperities of the counterface come into contact, leading to a
biased sampling of the counterface roughness by the deformable elastic-plastic
manifold, which retains deformation from the roughness it “saw” in contact.
Neither of the two distribution functions appears to be a perfect fit, but the
Gumbel distribution captures the largest values better than Weibull. This
clearly shows that the overall distribution has non-Gaussian tails. We note
that the saturated plasticity model also has non-Gaussian tails, as it samples
the peaks of the counterface. In contrast, the rigid-plastic (bearing area)
model simply reflects the Gaussian distribution of the counterface.

Figure 5: Evolution of the skewness of the height distribution with
indentation step $k$ for the $J_{2}$ plasticity model. The skewness follows a
power-law $k^{z}$, whose dynamics exponent $z$ is independent of the
counterface’s Hurst exponent. This exponent describes the dynamics of the
roughness evolution.

From the prior results is it clear that the height distribution is not
symmetric, which we now characterize from its skewness. Figure 5 shows, for
the $J_{2}$ model only, that the skewness of the height distribution,
$s^{\mathrm{p}}_{k}$, decreases with $k$ as a power-law that is independent of
the counterface’s Hurst exponent. Unlike the RMS of heights or slopes, there
is no observable transition from a power-law to a constant regime. We obtain a
power-law exponent of $z=-0.61(1)$ from a least-squares fit to the logarithm
of the data. Although $-1<z<0$ indicates a slow evolution towards a symmetric
height distribution, the skewness may converge to a finite value if the
surface hardens enough to prevent further plastic dissipation, which would
indicate elastic shakedown. Nonetheless, as we show next, the exponent $z$ is
representative of the evolution dynamics of the roughness at all scales.

Figure 6: (a) Evolution of the power-spectral density (PSD) with indentation
step for $H=0.8$, and (b) final power-spectral density for both $H=0.5$ and
$H=0.8$. Scaling the PSD with $k^{z}$ collapses all curves to a single master
curve. The master curve is constant at low wavenumbers (large wavelengths) and
has two distinct power-law regimes at intermediate and small wavelengths. The
intermediate scaling regime follows self-affine scaling of the counterface
only for $H=0.8$. The topographies obtained after indentation with $H=0.5$ and
$H=0.8$ show universal $q^{-8}$ scaling beyond the counterface’s short
wavelength cutoff.

Figure 6 (a) shows, for each indentation step $k$ and for $H=0.8$, the power-
spectral density (PSD) of the emerging roughness scaled by $k^{z}$: individual
PSDs then collapse onto a singular master curve above a wavenumber of $\sim$
10\. The scaled PSD shows a constant regime up to that wavenumber, where it
transitions to a regime with a Hurst exponent of $H^{\mathrm{p}}\approx 0.8$,
which corresponds to the counterface. Beyond the counterface’s short
wavelength cutoff (at wavenumber 128), the PSD scales as $q^{-8}$, which
corresponds to $H^{\mathrm{p}}=3$. This scaling regime is due to the
mechanical response at small scales: fig. 6b shows the PSDs of the last
indentation for counterfaces with $H=0.5$ and $H=0.8$, and both have identical
scaling beyond the short wavelength cutoff. The self-affine regime does not
exist in the simulations we conducted in the case of a counterface with
$H=0.5$. The transition from a constant (white noise) to power-law (self-
affine or elastic response) occurs on a range of scales too large to be
captured within the counterface’s spectral bandwidth ($q_{2}/q_{1}=64$). We
discuss below whether both resulting surfaces are in fact self-affine in the
regime where the counterface is self-affine.

Note that only “vertical” scaling of the PSD was required to obtain a
collapse, indicating that the final surface’s scaling behavior is already
encoded in the first indentations. This is corroborated by the behavior of
Nayak’s parameter in fig. 3 which is independent of the PSD’s absolute
magnitude, and shows little variation with $k$.

Figure 7: Variable bandwidth (VBM) analysis of (a) RMS heights and (b)
skewness, for the simulations with a counterface with $H=0.8$. The inset in
the left plot illustrates the successive surface subdivision which the VBM
employs.

In fig. 7, we show a multi-scale VBM analysis—illustrated by the schematic
inset in (a)—of two roughness parameters: the scale-dependent RMS of heights
(a), $h^{\mathrm{p}}_{k,\mathrm{rms}}(\ell)$, and skewness (b),
$s^{\mathrm{p}}_{k}(\ell)$, where $\ell$ is the window size. An infinite self-
affine surface shows a power-law behavior of
$h_{\mathrm{rms}}(\ell)\propto\ell^{H}$. Due to the restricted bandwidth of
the spectrum, the counterface follows this scaling approximately (but not
perfectly). In contrast to this, the emerging roughness is clearly not self-
affine.

The skewnews, shown in fig. 7b, is normalized by $k^{-z}$. The collapse of
skewness for all scales of curves with different $k$ corroborates the
observation that $z$ describes the roughness evolution dynamics at all scales.
As magnification increases, the skewness decreases, i.e. the height
distributions become more symmetric at smaller scales.

Figure 8: Height distributions after $90$ repeated contacts for magnifications
$L/\ell=1$ and $L/\ell=128$. The plot shows the results obtained for
counterfaces with $H=0.5$ and $H=0.8$. Normalization by
$h^{\mathrm{p}}_{\mathrm{rms}}(\ell)$ shows that the Hurst exponent does not
affect the distribution shape. At small scales, the height distribution
appears symmetric but non-Gaussian with exponential tails.

This behavior is directly shown in fig. 8, which depicts the height
distribution at magnifications $L/\ell=1$ and $128$ for both values of the
Hurst exponent at the last step $k=90$. The plot shows that the Hurst exponent
has no influence on the shape of the height distribution. At large
magnification, corresponding to a small wavelength cutoff of the counterface,
the height distribution is symmetric, although non-Gaussian, with
approximately exponential tails. This is roughly consistent with the Gumbel
tails observed at the length scale $L$ of the overall simulations.

##  4 Discussion

Our simulation show that some roughness properties, such as the RMS height or
slope, reach a steady-state long before the surface actually “shakes downs”,
in the sense that it stops deforming plastically and reacts only elastically
and reversibly. Indeed, our simulations never reach a true shakedown, but our
analysis of the dissipated energy indicates that the rate of plastic
deformation decreases with contact cycle, but only for the $J_{2}$ model that
explicitly considers hardening. Shakedown is never reached in the saturated
plasticity model without hardening, even though the roughness may appear to
have reached a state-state.

In fig. 3a, the RMS of heights for both models initially scales as $\sqrt{k}$.
Despite memory (from plasticity), and long-range correlations (from
elasticity), the RMS of heights behaves similarly to a simple ballistic
deposition model at short times [50], indicating that early iterations are
essentially uncorrelated. This can also be observed in the dissipated energy,
as early indentations are unaffected by hardening, while later indentations
occur on already hardened areas. However, since a perfect-plasticity
hypothesis also yields a steady-state, we conclude that hardening only plays a
minor role in the number of iterations required to observe a steady-state.

Another interesting aspect of our simulations is the emergence of non-Gaussian
height distributions. Non-Gaussianity is manifested in skewness as well as
exponential tails. The skewness has a well-defined monotonous dynamics
described by a simple power-law $s\propto k^{z}$ with $-1<z<0$. This type of
dynamics implies that the emerging rough geometry slowly becomes symmetric.

More insights are obtained by looking at the skewness as a function of scales.
Figure 7b clearly shows that the small-scale roughness is more symmetric than
the larger scales, and it appears that a factor of $k^{-z}$ collapses the
dynamic evolution of the skewness even at small scales. This seems to indicate
that during repeated indentation, the small-scale roughness diffuses towards
larger scales, or in other words that the small-scales are representative of
height-distribution at larger scales obtained at later indentation cycles. The
height distributions hence become symmetric in the asymptotic steady-state.

However, the small-scale roughness does not follow a Gaussian distribution but
has exponential tails (see fig. 8). The above argument seems to indicate, that
even in the steady-state where height distributions are symmetric they are
still not Gaussian. It is noteworthy that height distributions with
exponential tails are typically assumed in independent-asperity models of
rough contacts to make them analytically tractable [10, 51]. The height
distribution has also been shown to have implications in the way the true
contact area depends on the applied normal force [49], which could have
consequences in contact sealing [52, 53] and wear particle formation [54, 55,
56], although some of these applications would also see the development of
anisotropic roughness, due to sliding for example [13].

##  5 Summary & Conclusions

In summary, we have shown how surface roughness of an initially flat elastic-
plastic substrate evolves during repeated contact with a rigid self-affine
counterface. The resulting rough topography is neither self-affine nor
Gaussian, even if the rigid body has those properties. Common scalar roughness
measures, such as the RMS height, converge quickly to a steady-state value.
Yet, subsurface plastic deformation still occurs long after the topography
appears to have shaken-down geometrically. The topography continues to evolve,
but this can only be seen at intermediate scales obtained from a scale-
dependent geometric analyses, such as PSD.

The question of elastic shakedown has preoccupied tribologists for more than
half a century [10, 57, 1]. Shakedown occurs during complex phenomena such
friction [58], fretting [59] and rolling [60]. Our simplified simulations are
directly representative of processes such as sheet metal forming, where a die
is brought into repeated contact with a counterface [61]. However, they
certainly also hold insights for roughness evolution during sliding but
clearly cannot explain the development of anisotropic topographies [13]. Our
results contribute to the recent discussion on the emergence of surface
roughness on solid but deformable bodies [14, 19, 16].

#### Acknowledgments

The authors acknowledge the financial support of the Deutsche
Forschungsgemeinschaft (DFG Grant No. 461911253, AWEARNESS). Computational
resources from NEMO, at the University of Freiburg (DFG Grant No. INST
39/963-1 FUGG) are also acknowledged. The library Tamaas [34] was used for
contact simulation, roughness post-processing was done with Tamaas and
SurfaceTopography [62]. Blender and Matplotlib [63] were used for
visualization.

#### Author Contributions

L.F. and L.P. designed the study, discussed results and edited the final
manuscript. L.F. carried out calculations, analyzed the data and wrote the
first manuscript draft.

#### Data Availability

Simulation and figure codes are available on Zenodo [64].

#### Competing Interests

The authors declare no competing interests.

## References

  * [1] A.. Vakis et al.  “Modeling and Simulation in Tribology across Scales: An Overview”  In _Tribology International_ 125, 2018, pp. 169–199  DOI: 10.1016/j.triboint.2018.02.005
  * [2] François Renard, Thibault Candela and Elisabeth Bouchaud  “Constant Dimensionality of Fault Roughness from the Scale of Micro-Fractures to the Scale of Continents”  In _Geophysical Research Letters_ 40.1, 2013, pp. 83–87  DOI: 10.1029/2012GL054143
  * [3] B… Persson, O. Albohr, U. Tartaglino, A.. Volokitin and E. Tosatti  “On the Nature of Surface Roughness with Application to Contact Mechanics, Sealing, Rubber Friction and Adhesion”  In _Journal of Physics: Condensed Matter_ 17.1, 2005, pp. R1  DOI: 10.1088/0953-8984/17/1/R01
  * [4] Ludwig Binder  “Der Widerstand von Kontakten”  In _Elektrotechnik und Maschinenbau_ 30, 1912, pp. 781–788
  * [5] Ragnar Holm  “Eine Bestimmung der wirklichen Berührungsfläche eines Bürstenkontakes”  In _Wissenschaftliche Veröffentlichungen aus den Siemens-Werken_ 17.4, 1938, pp. 405–409
  * [6] Ragnar Holm  “Über die auf die wirkliche Berührungsfläche bezogene Reibkraft”  In _Wissenschaftliche Veröffentlichungen aus den Siemens-Werken_ 17.4, 1938, pp. 400–404
  * [7] Ragnar Holm  “Beitrag zur Kenntnis der Reibung”  In _Wissenschaftliche Veröffentlichungen aus den Siemens-Werken_ 20.1, 1941, pp. 68–84
  * [8] Ragnar Holm  “Electric Contacts: Theory and Applications”  Berlin ; New York: Springer, 2000
  * [9] K.. Johnson  “A Review of the Theory of Rolling Contact Stresses”  In _Wear_ 9.1, 1966, pp. 4–19  DOI: 10.1016/0043-1648(66)90010-X
  * [10] J.. Greenwood and J… Williamson  “Contact of Nominally Flat Surfaces”  In _Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences_ 295.1442, 1966, pp. 300–319  DOI: 10.1098/rspa.1966.0242
  * [11] Benoit B. Mandelbrot, Dann E. Passoja and Alvin J. Paullay  “Fractal Character of Fracture Surfaces of Metals”  In _Nature_ 308.5961, 1984, pp. 721–722  DOI: 10.1038/308721a0
  * [12] Laurent Ponson  “Statistical Aspects in Crack Growth Phenomena: How the Fluctuations Reveal the Failure Mechanisms”  In _International Journal of Fracture_ 201.1, 2016, pp. 11–27  DOI: 10.1007/s10704-016-0117-7
  * [13] Thibault Candela and Emily E. Brodsky  “The Minimum Scale of Grooving on Faults”  In _Geology_ 44.8, 2016, pp. 603–606  DOI: 10.1130/G37934.1
  * [14] Enrico Milanese, Tobias Brink, Ramin Aghababaei and Jean-François Molinari  “Emergence of Self-Affine Surfaces during Adhesive Wear”  In _Nature Communications_ 10.1, 2019, pp. 1116  DOI: 10.1038/s41467-019-09127-8
  * [15] Son Pham-Ba and Jean-François Molinari  “Creation and Evolution of Roughness on Silica under Unlubricated Wear”  In _Wear_ 472–473, 2021, pp. 203648  DOI: 10.1016/j.wear.2021.203648
  * [16] Ramin Aghababaei, Emily E. Brodsky, Jean-François Molinari and Srinivasan Chandrasekar  “How Roughness Emerges on Natural and Engineered Surfaces”  In _MRS Bulletin_ 47.12, 2022, pp. 1229–1236  DOI: 10.1557/s43577-022-00469-1
  * [17] Joaquin Garcia-Suarez, Tobias Brink and Jean-François Molinari  “Roughness Evolution Induced by Third-Body Wear”, 2023  DOI: 10.48550/arXiv.2306.08993
  * [18] P.. Thomson and P.. Nayak  “The Effect of Plastic Deformation on the Roughening of Free Surfaces of Sheet Metal”  In _International Journal of Machine Tool Design and Research_ 20.1, 1980, pp. 73–86  DOI: 10.1016/0020-7357(80)90020-7
  * [19] Adam R. Hinkle, Wolfram G. Nöhring, Richard Leute, Till Junge and Lars Pastewka  “The Emergence of Small-Scale Self-Affine Surface Roughness from Deformation”  In _Science Advances_ 6.7 American Association for the Advancement of Science, 2020, pp. eaax0847  DOI: 10.1126/sciadv.aax0847
  * [20] Wolfram G. Nöhring, Adam R. Hinkle and Lars Pastewka  “Nonequilibrium Plastic Roughening of Metallic Glasses Yields Self-Affine Topographies with Strain-Rate and Temperature-Dependent Scaling Exponents”  In _Physical Review Materials_ 6.7 American Physical Society, 2022, pp. 075603  DOI: 10.1103/PhysRevMaterials.6.075603
  * [21] A. Almqvist, F. Sahlin, R. Larsson and S. Glavatskih  “On the Dry Elasto-Plastic Contact of Nominally Flat Surfaces”  In _Tribology International_ 40.4, 2007, pp. 574–579  DOI: 10.1016/j.triboint.2005.11.008
  * [22] Avinash Tiwari, Anle Wang, Martin H. Müser and B… Persson  “Contact Mechanics for Solids with Randomly Rough Surfaces and Plasticity”  In _Lubricants_ 7.10 Multidisciplinary Digital Publishing Institute, 2019, pp. 90  DOI: 10.3390/lubricants7100090
  * [23] Marc Bonnet and Subrata Mukherjee  “Implicit BEM Formulations for Usual and Sensitivity Problems in Elasto-Plasticity Using the Consistent Tangent Operator Concept”  In _International Journal of Solids and Structures_ 33.30, 1996, pp. 4461–4480  DOI: 10.1016/0020-7683(95)00279-0
  * [24] C. Jacq, D. Nélias, G. Lormand and D. Girodin  “Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code”  In _Journal of Tribology_ 124.4, 2002, pp. 653  DOI: 10.1115/1.1467920
  * [25] Lucas Frérot, Marc Bonnet, Jean-François Molinari and Guillaume Anciaux  “A Fourier-accelerated Volume Integral Method for Elastoplastic Contact”  In _Computer Methods in Applied Mechanics and Engineering_ 351, 2019, pp. 951–976  DOI: 10.1016/j.cma.2019.04.006
  * [26] H.. Stanley and T. Kato  “An FFT-Based Method for Rough Surface Contact”  In _Journal of Tribology_ 119.3, 1997, pp. 481–485  DOI: 10.1115/1.2833523
  * [27] Carlos Campañá, Martin H. Müser and Mark O. Robbins  “Elastic Contact between Self-Affine Surfaces: Comparison of Numerical Stress and Contact Correlation Functions with Analytic Predictions”  In _Journal of Physics: Condensed Matter_ 20.35, 2008, pp. 354013  DOI: 10.1088/0953-8984/20/35/354013
  * [28] Vladislav A. Yastrebov, Guillaume Anciaux and Jean-François Molinari  “Contact between Representative Rough Surfaces”  In _Physical Review E_ 86.3, 2012
  * [29] B. Weber, T. Suhina, T. Junge, L. Pastewka, A.. Brouwer and D. Bonn  “Molecular Probes Reveal Deviations from Amontons’ Law in Multi-Asperity Frictional Contacts”  In _Nature Communications_ 9.1, 2018, pp. 888  DOI: 10.1038/s41467-018-02981-y
  * [30] Lucas Frérot, Guillaume Anciaux and Jean-François Molinari  “Crack Nucleation in the Adhesive Wear of an Elastic-Plastic Half-Space”  In _Journal of the Mechanics and Physics of Solids_ 145, 2020, pp. 104100  DOI: 10.1016/j.jmps.2020.104100
  * [31] J.. Simo and Thomas J.. Hughes  “Computational Inelasticity”, Interdisciplinary Applied Mathematics v. 7  New York: Springer, 1998
  * [32] Donald G. Anderson  “Iterative Procedures for Nonlinear Integral Equations”  In _Journal of the ACM_ 12.4, 1965, pp. 547–560  DOI: 10.1145/321296.321305
  * [33] V. Eyert  “A Comparative Study on Methods for Convergence Acceleration of Iterative Vector Sequences”  In _Journal of Computational Physics_ 124.2, 1996, pp. 271–285  DOI: 10.1006/jcph.1996.0059
  * [34] Lucas Frérot, Guillaume Anciaux, Valentine Rey, Son Pham-Ba and Jean-François Molinari  “Tamaas: A Library for Elastic-Plastic Contact of Periodic Rough Surfaces”  In _Journal of Open Source Software_ 5.51, 2020, pp. 2121  DOI: 10.21105/joss.02121
  * [35] David Tabor  “The Hardness of Metals”, Monographs on the Physics and Chemistry of Materials  Oxford: Clarendon Press, 1951
  * [36] Y.. Gao, A.. Bower, K.-S. Kim, L. Lev and Y.. Cheng  “The Behavior of an Elastic–Perfectly Plastic Sinusoidal Surface under Contact Loading”  In _Wear_ 261.2, 2006, pp. 145–154  DOI: 10.1016/j.wear.2005.09.016
  * [37] Hamid Ghaednia, Xianzhang Wang, Swarna Saha, Yang Xu, Aman Sharma and Robert L. Jackson  “A Review of Elastic–Plastic Contact Mechanics”  In _Applied Mechanics Reviews_ 69.6, 2017, pp. 060804-060804–30  DOI: 10.1115/1.4038187
  * [38] L. Pei, S. Hyun, J.-F. Molinari and Mark O. Robbins  “Finite Element Modeling of Elasto-Plastic Contact between Rough Surfaces”  In _Journal of the Mechanics and Physics of Solids_ 53.11, 2005, pp. 2385–2409  DOI: 10.1016/j.jmps.2005.06.008
  * [39] Jiunn-Jong Wu  “Simulation of Rough Surfaces with FFT”  In _Tribology International_ 33.1, 2000, pp. 47–58  DOI: 10.1016/S0301-679X(00)00016-5
  * [40] Srinivasa B. Ramisetti, Carlos Campañá, Guillaume Anciaux, Jean-Francois Molinari, Martin H. Müser and Mark O. Robbins  “The Autocorrelation Function for Island Areas on Self-Affine Surfaces”  In _Journal of Physics: Condensed Matter_ 23.21, 2011, pp. 215004  DOI: 10.1088/0953-8984/23/21/215004
  * [41] Tevis D.. Jacobs, Till Junge and Lars Pastewka  “Quantitative Characterization of Surface Topography Using Spectral Analysis”  In _Surface Topography: Metrology and Properties_ 5.1, 2017, pp. 013001  DOI: 10.1088/2051-672X/aa51f8
  * [42] H.. Bui  “Some Remarks about the Formulation of Three-Dimensional Thermoelastoplastic Problems by Integral Equations”  In _International Journal of Solids and Structures_ 14.11, 1978, pp. 935–939  DOI: 10.1016/0020-7683(78)90069-0
  * [43] Raymond D. Mindlin  “Force at a Point in the Interior of a Semi-Infinite Solid”  In _Journal of Applied Physics_ 7.5, 1936, pp. 195–202  DOI: 10.1063/1.1745385
  * [44] Lucas Frérot  “The Mindlin Fundamental Solution - A Fourier Approach”, 2018  Zenodo  DOI: 10.5281/zenodo.1492149
  * [45] P. Nayak  “Random Process Model of Rough Surfaces”  In _Journal of Lubrication Technology_ 93.3, 1971, pp. 398–407  DOI: 10.1115/1.3451608
  * [46] Michael J. Cannon, Donald B. Percival, David C. Caccia, Gary M. Raymond and James B. Bassingthwaighte  “Evaluating Scaled Windowed Variance Methods for Estimating the Hurst Coefficient of Time Series”  In _Physica A: Statistical Mechanics and its Applications_ 241.3, 1997, pp. 606–626  DOI: 10.1016/S0378-4371(97)00252-5
  * [47] C.-K. Peng, S.. Buldyrev, S. Havlin, M. Simons, H.. Stanley and A.. Goldberger  “Mosaic Organization of DNA Nucleotides”  In _Physical Review E_ 49.2 American Physical Society, 1994, pp. 1685–1689  DOI: 10.1103/PhysRevE.49.1685
  * [48] C.-K. Peng, Shlomo Havlin, H. Stanley and Ary L. Goldberger  “Quantification of Scaling Exponents and Crossover Phenomena in Nonstationary Heartbeat Time Series”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 5.1, 1995, pp. 82–87  DOI: 10.1063/1.166141
  * [49] T. Silva Sabino, A.. Couto Carneiro, R. Pinto Carvalho and F.. Andrade Pires  “Evolution of the Real Contact Area of Self-Affine Non-Gaussian Surfaces”  In _International Journal of Solids and Structures_ 268, 2023, pp. 112173  DOI: 10.1016/j.ijsolstr.2023.112173
  * [50] A.- L. Barabási and H.. Stanley  “Fractal Concepts in Surface Growth”  Cambridge University Press, 1995  DOI: 10.1017/CBO9780511599798
  * [51] Bo N.. Persson  “Sliding Friction”, NanoScience and Technology  Berlin, Heidelberg: Springer Berlin Heidelberg, 2000  DOI: 10.1007/978-3-662-04283-0
  * [52] Wolf B. Dapp, Andreas Lücke, Bo N.. Persson and Martin H. Müser  “Self-Affine Elastic Contacts: Percolation and Leakage”  In _Physical Review Letters_ 108.24, 2012, pp. 244301  DOI: 10.1103/PhysRevLett.108.244301
  * [53] A.. Shvarts and V.. Yastrebov  “Trapped Fluid in Contact Interface”  In _Journal of the Mechanics and Physics of Solids_ 119, 2018, pp. 140–162  DOI: 10.1016/j.jmps.2018.06.016
  * [54] Lucas Frérot, Ramin Aghababaei and Jean-François Molinari  “A Mechanistic Understanding of the Wear Coefficient: From Single to Multiple Asperities Contact”  In _Journal of the Mechanics and Physics of Solids_ 114, 2018, pp. 172–184  DOI: 10.1016/j.jmps.2018.02.015
  * [55] Valentin L. Popov and Roman Pohrt  “Adhesive Wear and Particle Emission: Numerical Approach Based on Asperity-Free Formulation of Rabinowicz Criterion”  In _Friction_ , 2018  DOI: 10.1007/s40544-018-0236-4
  * [56] Tobias Brink, Lucas Frérot and Jean-François Molinari  “A Parameter-Free Mechanistic Model of the Adhesive Wear Process of Rough Surfaces in Sliding Contact”  In _Journal of the Mechanics and Physics of Solids_ 147, 2021, pp. 104238  DOI: 10.1016/j.jmps.2020.104238
  * [57] A. Kapoor, J.. Williams and K.. Johnson  “The Steady State Sliding of Rough Surfaces”  In _Wear_ 175.1, 1994, pp. 81–92  DOI: 10.1016/0043-1648(94)90171-6
  * [58] R.. Flicek, D.. Hills, J.. Barber and D. Dini  “Determination of the Shakedown Limit for Large, Discrete Frictional Systems”  In _European Journal of Mechanics - A/Solids_ 49, 2015, pp. 242–250  DOI: 10.1016/j.euromechsol.2014.08.001
  * [59] S. Fouvry, Ph. Kapsa and L. Vincent  “An Elastic–Plastic Shakedown Analysis of Fretting Wear”  In _Wear_ 247.1, 2001, pp. 41–54  DOI: 10.1016/S0043-1648(00)00508-1
  * [60] L. Berthe, P. Sainsot, A.. Lubrecht and M.. Baietto  “Plastic Deformation of Rough Rolling Contact: An Experimental and Numerical Investigation”  In _Wear_ 312.1–2, 2014, pp. 51–57  DOI: 10.1016/j.wear.2014.01.017
  * [61] Sy-Wei Lo and Tzu-Chern Horng  “Surface Roughening and Contact Behavior in Forming of Aluminum Sheet”  In _Journal of Tribology_ 121.2, 1999, pp. 224–233  DOI: 10.1115/1.2833925
  * [62] Michael C. Röttger, Antoine Sanner, Luke A. Thimons, Till Junge, Abhijeet Gujrati, Joseph M. Monti, Wolfram G. Nöhring, Tevis D.. Jacobs and Lars Pastewka  “Contact.Engineering—Create, Analyze and Publish Digital Surface Twins from Topography Measurements across Many Scales”  In _Surface Topography: Metrology and Properties_ 10.3 IOP Publishing, 2022, pp. 035032  DOI: 10.1088/2051-672X/ac860a
  * [63] John D. Hunter  “Matplotlib: A 2D Graphics Environment”  In _Computing in Science & Engineering_ 9.3, 2007, pp. 90–95  DOI: 10.1109/MCSE.2007.55
  * [64] Lucas Frérot and Lars Pastewka  “Supplementary Codes and Data to "Elastic Shakedown and Roughness Evolution in Repeated Elastic-Plastic Contact"”, 2023  Zenodo  DOI: 10.5281/zenodo.8280362

Figure 9: Graphical Abstract