# Efficient Exploration of Multi-Modal Posterior Distributions

Yi-Ming Hu<EMAIL_ADDRESS>Martin Hendry  Ik Siong Heng  School of
Physics and Astronomy, University of Glasgow, Glasgow, G12 8QQ, UK

###### Abstract

The Markov Chain Monte Carlo (MCMC) algorithm is a widely recognised as an
efficient method for sampling a specified posterior distribution. However,
when the posterior is multi-modal, conventional MCMC algorithms either tend to
become stuck in one local mode, become non-Markovian or require an excessively
long time to explore the global properties of the distribution. We propose a
novel variant of MCMC, _mixed MCMC_ , which exploits a specially designed
proposal density to allow the generation candidate points from any of a number
of different modes. This new method is efficient by design, and is strictly
Markovian. We present our method and apply it to a toy model inference problem
to demonstrate its validity.

††preprint: APS/123-QED

##  I Introduction

Bayesian inference methods have been applied to an increasingly wide range of
data analysis problems in the physical sciences – particularly problems
requiring _parameter estimation_ and _model selection_ (see, for example,
Gregory (50 X)Sivia (back)Rover (2007)Loredo (2013)). This growing popularity
of the Bayesian approach has been driven by the ready availability of
increasingly powerful computational facilities, allied with significant
advances in the data analysis methodology and algorithms that have been
developed. These twin drivers have permitted the application of Bayesian
methods to data sets of a size, complexity and dimensionality that until
recently would have rendered them intractable.

One of the main targets of Bayesian inference is to estimate the posterior
distribution of desired parameters. To estimate posterior distributions, one
naive solution is to use an exhaustive algorithm to calculate the posterior
over a dense grid of points in the parameter space. Such brute-force methods
will have little or no practical value when dealing with medium-to-high
dimensional problems since the computational burden will be prohibitively
high. For such problems the ability to concentrate sampling in regions where
the posterior probability is high is very important if we are to implement
Bayesian inference methods efficiently.

Methods such as _Markov Chain Monte Carlo (MCMC)_ and _Nested Sampling_ are
well tailored to explore the posterior distribution over high dimensional
parameter spaces. While the computational cost of brute-force methods
increases exponentially with the dimension, MCMC usually only grows linearly
with dimension Allison and Dunkley (2013)Veitch and Vecchio (2010).

Generally, the method of MCMC works well so long as the posterior surface is
sufficiently smooth. However, when the posterior distribution has a
complicated structure, MCMC will become inefficient. For example, MCMC
samplers are known to get “caught” in a local mode of the posterior, and
unable to jump out and explore any other isolated modes in the parameter space
Kirkpatrick  _et al._ (1983)Neal (1996). So a lot of methods have been
proposed to make the MCMC sampling more efficient (e.g. Tierney and Mira
(1999)Trias  _et al._ (2009)Mira (2001))

In this work, we propose a novel method called _mixed MCMC_ to deal with such
issues. The conventional MCMC algorithm is robust for exploring the detailed
structure of the posterior surface, and we want to retain that property while
enabling some global “communication” between different regions of the
parameter space so that the sampler can make jumps between those regions
without requiring a very long exploration time.

This paper is structured as follows. In section II, we introduced the general
realisation of MCMC and discussed how the MCMC can be used as a Bayesian tool.
In section III, we discuss the main difference of mixed MCMC and conventional
MCMC and give the pseudo code for its realisation. We apply the method to a
toy model as illustrated in section IV. Finally in section V, we summarise the
motivation and properties of mixed MCMC, and discussed possible extension to
it.

##  II Method

The basic principles of an MCMC algorithm are simply stated. The algorithm
sets out to sample a chain of points in the parameter space and at the
$i^{th}$ iteration (i.e. after $i-1$ points have already been sampled) a
_candidate_ point $\bm{\theta^{*}}$ is randomly sampled from some specified
proposal distribution, based solely on the position of the _previous point_ in
the chain $\bm{\theta^{(i-1)}}$. The corresponding posterior for this
candidate point is calculated, and compared with the posterior at
$\bm{\theta^{(i-1)}}$. If the value of the posterior at the candidate point is
larger than that of the previous point, the candidate point is accepted as the
next point in the chain. Otherwise, the candidate is accepted only with a
certain _acceptance probability_ (see next section). One finds therefore that
the sampling will generally proceed “uphill” – i.e. to regions of the
parameter space where the value of the posterior is larger – while sometimes
it can also go “downhill” to regions where the posterior takes on lower
values. The precise form of the acceptance probability achieves what is termed
_detailed balance_ , which ensures that the chain of sampled points is indeed
a random sample drawn from the desired posterior distribution. This method can
thus be used to efficiently explore the posterior distribution, avoiding the
need for a global optimization via an exhaustive grid search.Gregory (50
X)Sivia (back)Rover (2007)Loredo (2013)

###  II.1 Markov Chain Monte Carlo

Interested reader is refered to Gregory (50 X) and Sivia (back) for detailed
discussion of Bayesian Inference. Hereafter, we define the posterior
$f(\bm{\theta})=p(\bm{\theta}|D,I)$, the prior
$\pi(\bm{\theta})=p(\bm{\theta},I)$ and likelihood
$\ell(\bm{\theta})=p(D|\bm{\theta},I)$, where $\bm{\theta}$ is the parameter
set, $D$ is the data and $I$ is the information.

The simplest form of Markov Chain Monte Carlo (MCMC) is known as the
Metropolis algorithm, which can be achieved by the following steps Gregory (50
X)Sivia (back)Rover (2007)Press  _et al._ (2007).

  1. 1.

Arbitrarily choose a starting point $\bm{\theta}^{(0)}$ that satisfies
$f(\bm{\theta}^{(0)})>0$, and a symmetric proposal distribution
$J(\bm{\theta}_{a}|\bm{\theta}_{b})$. Set step index _i_ =0.

  2. 2.

Increment _i_ by 1.

  3. 3.

Randomly propose a new parameter set $\bm{\theta}^{*}$ by sampling from
$J(\cdot|\bm{\theta}^{(i-1)})$.

  4. 4.

Calculate the _Metropolis ratio_ given by

$r=\frac{f(\bm{\theta}^{*})}{f(\bm{\theta}^{(i-1)})}$ (1)

  5. 5.

Accept the proposed parameter set $\bm{\theta}^{*}$ with _acceptance
probability_

$\alpha(\bm{\theta}^{(i-1)},\bm{\theta}^{*})\triangleq\min(1,r)$ (2)

If $r\geq 1$, then the candidate is accepted, so the new point is
$\bm{\theta}^{(i)}=\bm{\theta}^{*}$.

If $r<1$, draw a random number $rand$ from a uniform distribution $U[0,1]$,
and if $rand<r$, then set $\bm{\theta}^{(i)}=\bm{\theta}^{*}$; otherwise set
$\bm{\theta}^{(i)}=\bm{\theta}^{(i-1)}$.

Step 2-5 are repeated until a large enough number of points have been sampled.
This termination could be controlled by a preset number, or by monitoring the
samples’ distribution and check if it’s sufficiently stable.Gelman and Rubin
(1992) The beginning peorid, which is generaly called as “burn-in” stage, is
discarded to prevent the influence of the arbitrary choice of starting point
$\bm{\theta}^{(0)}$.

The Metropolis-Hastings (M-H) algorithm is a more general form of the
Metropolis algorithm. In the Metropolis algorithm, the proposal distribution
is symmetric, that is
$J(\bm{\theta}_{a}|\bm{\theta}_{b})=J(\bm{\theta}_{b}|\bm{\theta}_{a})$, but
this condition is not necessary. In the M-H algorithm we relax this symmetric
condition, so that equation (1) should be modified as follows

$r=\frac{f(\bm{\theta}^{*})J(\bm{\theta}^{(i-1)}|\bm{\theta}^{*})}{f(\bm{\theta}^{(i-1)})J(\bm{\theta}^{*}|\bm{\theta}^{(i-1)})}.$
(3)

It is clear that when the proposal distribution is symmetric, equation (3) is
identical to equation (1).

It can be shown that the number density of the sampled points will represent a
sample from the posterior distribution Gregory (50 X). Thus estimation of the
parameter(s) that characterise the posterior distribution becomes possible
with a sufficiently large number of sampling points.

###  II.2 Convergence of MCMC

We can estimate the posterior distribution from a histogram of the values
sampled by our MCMC chain, and the mean of the parameters can be estimated
trivially as

$\frac{1}{N}\sum_{i=1}^{N}\bm{\theta}^{(i)}$

Any realization of MCMC is guaranteed to be converged if it satisfies the
requirement of _detailed balance_.Gregory (50 X)

$f(\bm{\theta}^{(i)})J(\bm{\theta}^{(i-1)}|\bm{\theta}_{t})=f(\bm{\theta}^{(i-1)})J(\bm{\theta}^{(i)}|\bm{\theta}^{(i-1)})$
(4)

The concept of detailed balance in thermodynamics can help us to understand
this requirement for the convergence of the MCMC chain. In thermodynamics, we
can define the probability of a particle to be in state $\bm{\theta}^{(i)}$ as
$f(\bm{\theta}^{(i)})$, and the probability of the particle to jump to state
$\bm{\theta}^{(i-1)}$ as $J(\bm{\theta}^{(i-1)}|\bm{\theta}^{(i)})$. Detailed
balance requires that, after a sufficiently long timescale, the probability
for a particle to jump from state $\bm{\theta}^{(i)}$ to state
$\bm{\theta}^{(i-1)}$ should be exactly the same as the probability to jump
from state $\bm{\theta}^{(i-1)}$ to state $\bm{\theta}^{(i)}$.

We need to notice that detailed balance is a stronger requirement than
convergence, in the sense that a Markovian Chain that is not in detailed
balanced may still converge to the target distribution.Foreman-Mackey  _et
al._ (2013)

##  III mixed MCMC

If the starting point and/or the proposal density is not properly chosen, the
MCMC sampler might become stuck in a local mode, and will not be able to
appropriately explore the whole parameter space. This might introduce a
statistical bias in the parameter estimation carried out by MCMC, particularly
when the target distribution is multi-modal. Thus motivates the realisation of
mixed MCMC as a really Markovain realisation of MCMC that can sample posterior
efficiently.Neal (1996)Gelman and Rubin (1992)Foreman-Mackey  _et al._ (2013)

Here we propose a novel method which we term _mixed MCMC_ to perform Bayesian
inference on multi-modal posterior distributions. This method can allow the
sampler to communicate between different local maxima, so that the sampler
will be able to represent local peaks, as well as to explore the global
structure. As noted previously, our method requires some limited information
about the location of the multiple modes before sampling. In many cases,
however, we will have at least some rough prior knowledge about the posterior,
and we can use this information to guide the sampler. Even in the absence of
such prior knowledge, other existing global sampling methods can be tailored
for this purpose speed up this process.Neal (1996)Feroz and Skilling
(2013)Crowder  _et al._ (2006)

###  III.1 Algorithm

The main difference between the algorithm for mixed MCMC and the conventional
MCMC algorithms simply roots in the use of a novel form of proposal density.
The sampler should be able to generate candidates from different sub-regions,
while proper choice of Metropolis ratio will ensure that the sampling between
those different sub-regions satisfies detailed balance.

Suppose, as a result of existing prior knowledge, or with the help of some
other global sampling method, we have some information about the posterior
distribution that is sufficient to identify the existence and the rough
location of the several modes in posterior distribution, where the location of
the $t^{th}$ mode is labeled as $\bm{\theta^{0}_{t}}$. We can then divide the
parameter space into several distinct _sub-regions_ each of which we assume
contains a single mode of the posterior Gregory (50 X)Gregory  _et al._
(2007).

We should bare in mind that this method is designed for multi-modal posterior,
thus the proposal density should be designed in a way that it can propose new
candidates in all posterior modes. Thus we assign to the $t^{th}$ sub-region
what we term a _picking up probability_ , $p_{t}$, which determines the
probability to get a new candidate in the $t^{th}$ sub-region. Ideally, this
probability should be the same as the marginal likelihood (also known as the
evidence) within the sub-region – i.e. the probability that the candidate
point lies within that sub-region. Note also that the picking up probability
should satisfy the normalisation requirement _ $\sum\limits_{t}p_{t}=1$_. At
the same time it will maximise the efficiency of our approach if
$p_{t}\propto\int_{V_{t}}f(\bm{\theta})d\bm{\theta}$, where $V_{t}$ is the
volume of the $t^{th}$ sub-region of the parameter space.

Suppose we decide to generate a candidate point in the $t^{th}$ sub-region,
while the current (i.e. most recently updated) point $\bm{\theta^{(i-1)}}$ is
located in the $s^{th}$ sub-region. Then a normalised multivariate
distribution (most conveniently taken to be a Gaussian) centering around the
point $\bm{\theta^{0}_{t}}-\bm{\theta^{0}_{s}}+\bm{\theta^{(i-1)}}$ is used as
proposal density, and a candidate is drawn from this distribution. After
calculating the value of the posterior at this candidate point, and then
computing the Metropolis ratio, $r$, in the usual way, we can decide to accept
the candidate point with the acceptance probability $\alpha$ as before.

In more detail our mixed MCMC algorithm can be illustrated with the following
pseudo-code.

  1. 1.

Obtain some rough approximation to the posterior distribution using other
methods. Identify _m_ modes in the parameter space, and estimate their central
locations given by $\bm{\theta^{0}_{t}}$.

  2. 2.

Set $p_{t}$ to be the picking up probability, defined as proportional to the
volume of the $t^{th}$ sub-region, with $\sum_{t}p_{t}=1$. Set step label
$i=0$

  3. 3.

Randomly pick a starting point, $\bm{\theta^{(i)}}$.

  4. 4.

_while_(not converged)

    1. (a)

Set i=i+1

    2. (b)

Randomly pick a sub-region number $t$ with probability $p_{t}$ and assign $s$
to be the current sub-region index. $s,~{}t\in\\{1,...,m\\}$ where m is the
number of all sub-regions.

    3. (c)

Generate the candidate point
$\bm{\theta^{*}}=\bm{\theta^{(i-1)}}+\bm{\theta^{0}_{t}}-\bm{\theta^{0}_{s}}+\bm{\delta\theta}$
drawn from the proposal density $\sim J(\cdot|\bm{\theta^{(i)}})$

    4. (d)

Calculate the Metropolis ratio $r$ based on the candidate and the previous
point, $r=\rm\frac{f(\bm{\theta^{*}})p_{t}}{f(\bm{\theta^{(i-1)}})p_{s}}$.

    5. (e)

Generate a random number ${\rm rand}\sim U[0,1]$.

    6. (f)

Accept the proposed parameter set $\bm{\theta}^{*}$ with _acceptance
probability_ $\alpha=\min(1,r)$ as follows:

_if_($r>{\rm rand}$), update, $\bm{\theta^{(i)}}=\bm{\theta^{*}}$
_else_ $\bm{\theta^{(i)}}=\bm{\theta^{(i-1)}}$

The mixed MCMC algorithm set out above is strictly Markovian, and detailed
balance is achieved by construction. Thus the number of points sampled in
given sub-region should provide an estimate of the local evidence. As noted
above, in order to maximise the efficiency of the algorithm the picking up
probability $p_{t}$ should better be proportional to the local evidence.

Also, we can notice that when the proposed point and the previous point are
located in the same sub-region, then the algorithm reduces to the conventional
M-H algorithm, which further verifies its validity.

##  IV Toy Model

We demonstrate our mixed MCMC algorithm using a simple toy model. On a two
dimensional $x-y$ parameter space, we considered a posterior distribution is
the sum of a pair of well-separated bivariate normal distributions. Where the
parameters are $x$ and $y$, and the two artificial posterior modes locate in
$(\mu^{x}_{1},\mu^{y}_{1})$ and $(\mu^{x}_{2},\mu^{y}_{2})$, each mode can be
described by a bivariate normal distribution with a diagonal covariance
matrix, where standard deviation in each direction is
$\sigma^{x}_{1},\sigma^{y}_{1}$ for the first mode and
$\sigma^{x}_{2},\sigma^{y}_{2}$ for the second mode. The form of this
posterior is, therefore:

$\begin{split}f(x,y)=C_{1}\exp[-\frac{(x-\mu^{x}_{1})^{2}}{2(\sigma^{x}_{1})^{2}}-\frac{(y-\mu^{y}_{1})^{2}}{2(\sigma^{y}_{1})^{2}}]\\\
+C_{2}\exp[-\frac{(x-\mu^{x}_{2})^{2}}{2(\sigma^{x}_{2})^{2}}-\frac{(y-\mu^{y}_{2})^{2}}{2(\sigma^{y}_{2})^{2}}].\end{split}$
(5)

The coefficients $C_{1}$ and $C_{2}$ allow the two modes to differ in height,
and when integrated over the entire parameter space the normalization
condition implies that

$\int_{x,y}\mathrm{d}x~{}\mathrm{d}y~{}f(x,y)=1.$ (6)

For simplicity, we chose $\mu_{1}^{x}=-\mu_{2}^{x}=-3$,
$\sigma_{1}^{x}=\sigma_{2}^{x}=0.1$, $\mu_{y}^{1}=\mu_{y}^{2}=0$,
$\sigma_{1}^{y}=\sigma_{2}^{y}=0.1$ and the ratio of two coefficients
$C_{1}:C_{2}$ is kept as $1:3$.

In this toy model test, we only concentrate on the validity of the mixed MCMC
method, and do not consider in detail other factors such as its efficiency or
generality. Thus, we assume prior knowledge of the separated structure of the
posterior distribution. Given this assumption, it is possible to analytically
calculate the $\Delta\chi^{2}$ value that corresponds to the contour within
which a certain fraction of the entire volume of the posterior is located,
thus providing us with an exact theoretical reference result with which to
compare. For a one-dimensional Gaussian distribution, the $1\sigma$, $2\sigma$
and $3\sigma$ credible regions correspond to $68.27\%,~{}95.45\%$ and
$99.73\%$ of the cumulative probability function (CDF) respectively. Thus, it
is convenient to consider for our toy model posterior the $\Delta\chi^{2}$
values that correspond to the $68.27\%,~{}95.45\%$ and $99.73\%$ of the CDF,
and compare it with the sample estimates obtained from application of our
mixed MCMC algorithm.

Under the Gaussian assumption, $\Delta\chi^{2}=-\Delta 2\log(\cal L)$, so in
the general case we will compare the value of $2\Delta\log(\cal L)$ with its
theoretical evaluation. Details of this theoretical calculation can be found
in the Appendix. The posterior of the toy model can be taken as two
independent bivariant Gaussians, each with a diagonal covariance matrix.

We generated a chain with $10^{5}$ points. The proposal density was set to be
a bivariate Gaussian distribution in addition to the probabilistic shift
between sub-regions, with covariance matrix equal to the identity matrix
multiplied by $\sigma=0.1$.

For this particular toy model, theoretically the corresponding
$68.27\%,~{}95.45\%$ and $99.73\%$ credible regions should have values of
$2\Delta\log(\cal L)$ equal to $3.11,~{}6.99$ and $12.64$. A typical
realisation gives result as $3.14,~{}7.04$ and $12.24$, and the numbers of
points sampled in the two sub-regions are 2559 and 7441, which is consistent
with the $1:3$ ratio assumed for the coefficients $C_{1}$ and $C_{2}$.

In figure 1 an example of the sampling results is shown, with blue, green and
red colour points representing the highest (i.e. largest value of the
posterior) $68.27\%,~{}95.45\%~{}\rm{and}~{}99.73\%$ fraction of the samples,
after sorting the posterior values in descending order.

Figure 1: Example realisation of mixed MCMC applied to the toy model posterior
distribution described in the text. Here blue, green and red colour points
represent subsets of the sampled points with the highest
$68.27\%,~{}95.45\%~{}\rm{and}~{}99.73\%$ fractions respectively of the
posterior value among all samples. The different areas of the two different
modes reflect the fact that the two modes have different weights. Like all
MCMC results, the density of the sampling points reflect the posterior.

##  V Discussion

A novel method, which we call _mixed MCMC_ , has been proposed. In the
situations when the multi-modal characteristics of the posterior distribution
are already roughly known, the parameter space can be split into several sub-
regions, each of which hosts a single mode, and our mixed MCMC method can be
applied. The proposal density can generate candidates in different sub-regions
by adding a shift from the current sub-region to the proposed new sub-region.
In this way, a comparison between different sub-regions can be done globally,
which improves efficiency. This algorithm is strictly Markovian, so the
detailed balance requirement is fulfilled. The concept of mixed MCMC is
realised by enabling proposed candidate points to be generated from different
modes of the posterior. Admittedly, the mixed MCMC approach must rely on other
methods to first identify the multiple modes of the posterior distribution.
However, since that identification will generally require only rough
information, we can expect this initial stage to be rapid. Moreover, other
existing methods already provide some solutions to the problem of identifying
multiple modesFeroz  _et al._ (2009)Feroz and Skilling (2013).

From another perspective, if we view separate parameter subspaces as different
models, this mixed MCMC algorithm can be viewed as an special form of
reversible jump MCMCGreen (1995)Green and Mira (2001), which can sample from
different models even when they have different dimensionality, and thus
provides the Bayesian odds ratio of two models.

This method is a novel realisation of MCMC which can achieve high efficiency
in analysing multi-modal posterior distributions by virtue of its unique form
of proposal density. It relies not only on local information, but also on the
global structure through swapping between different sub-chains. In particular
the candidate point is accepted with an _acceptance probability_
$a=\min(1,r)$, where _r_ is the _Metropolis ratio_ , which takes into account
the global information about the multiple modes of the posterior.

So far, we have not discussed in detail how to obtain rough information about
the posterior modes. We note that methods such as MultiNestFeroz  _et al._
(2009)Feroz and Skilling (2013) aim to solve similar problems, so their
approach could be directly applied here. Some other methods like parallel
tempering MCMCNeal (1996) or k-means MacQueen (1967) can also be modified and
applied here.

We leave the detailed comparison with other methods, like parallel tempering
MCMC to future work. However, by not throwing away points in parallel chains,
and the design of the proposal density to have a relatively short
autocorrelation length we expect the mixed MCMC algorithm to be quite
efficient Peskun (1973).

In this work, we also applied our method to a simple toy model, with two
distinct well-separated modes, to demonstrate its efficacy. With $10^{5}$
samples, our mixed MCMC was able to both find the picking up probability,
which represents the bulk distribution of the posterior (i.e. the probability
of belonging to each mode) and also the Bayesian credible regions for the
posterior as a whole – each of which show excellent agreement with the exact,
theoretically computed values for our toy model.

This toy model investigation shows that the idea of mixed MCMC is
theoretically sound and practically useful. Further investigation can be done
in future work by investigating ways of optimising the local proposal
densities individually so that the efficiency of the method can be further
improved.

##  Appendix A Analytical Evaluation of $2\Delta\log(\cal L)$

In this section we present the calculation of $2\Delta\log(\cal L)$, which is
defined as

$2\Delta\log({\cal L})=2\log({\cal L}_{max})-2\log({\cal L}),$

for our toy model posterior.

In the toy model, there are two well separated modes, and we can simply assume
they are fully independent. The two modes have an evidence ratio of
$w_{1}:w_{2}$, where $w_{1}+w_{2}=1$ is the normalisation requirement.

We define the two independent parts of the posterior as

$P_{1}(x,y)=\frac{1}{2\pi\sigma^{x}_{1}\sigma^{y}_{1}}\exp[-\frac{(x-\mu^{x}_{1})^{2}}{2(\sigma^{x}_{1})^{2}}-\frac{(y-\mu^{y}_{1})^{2}}{2(\sigma^{y}_{1})^{2}}]$
$P_{2}(x,y)=\frac{1}{2\pi\sigma^{x}_{2}\sigma^{y}_{2}}\exp[-\frac{(x-\mu^{x}_{2})^{2}}{2(\sigma^{x}_{2})^{2}}-\frac{(y-\mu^{y}_{2})^{2}}{2(\sigma^{y}_{2})^{2}}]$

and the posterior can be written as

$f(x,y)=w_{1}P_{1}(x,y)+w_{2}P_{2}(x,y)$ (7)

The peak values of the posterior for its two modes are
$C_{1}\triangleq\frac{w_{1}}{2\pi\sigma_{1}^{x}\sigma_{1}^{y}}$ and
$C_{2}\triangleq\frac{w_{2}}{2\pi\sigma_{2}^{x}\sigma_{2}^{y}}$ respectively.

For simplicity, we replace
$\frac{(x-\mu^{x}_{1})^{2}}{(\sigma^{x}_{1})^{2}}+\frac{(y-\mu^{y}_{1})^{2}}{(\sigma^{y}_{1})^{2}}=r_{1}^{2}$
and
$\frac{(x-\mu^{x}_{2})^{2}}{(\sigma^{x}_{2})^{2}}+\frac{(y-\mu^{y}_{2})^{2}}{(\sigma^{y}_{2})^{2}}=r_{2}^{2}$
and rewrite the posterior as

$f(x,y)=C_{1}\exp(-r_{1}^{2}/2)+C_{2}\exp(-r_{2}^{2}/2)$ (8)

Without losing generality, we assume $C_{1}>C_{2}$, and so the highest
posterior value $f_{max}=C_{1}$, and highest posterior value for the secondary
peak is $C_{2}$. We define $r_{0}$ as $C_{1}\exp(-r_{0}^{2}/2)=C_{2}$,
equivalently, $\exp(-r_{0}^{2}/2)=\frac{C_{2}}{C_{1}}$.

Our aim is to find the expression for $\Delta\chi^{2}(C)$, so that given $C$,
we have

$\int\limits_{f>\exp(-\frac{\Delta\chi^{2}}{2})}f(x,y)~{}\mathrm{d}x~{}\mathrm{d}y=C$

When $f>\exp(-r_{0}^{2}/2)$,

$\displaystyle C$ $\displaystyle=$
$\displaystyle\int\limits_{f>\exp(-\Delta\chi^{2}/2)}C_{1}\exp(-\frac{r^{2}}{2})\mathrm{d}x~{}\mathrm{d}y$
(9) $\displaystyle=$ $\displaystyle 2\pi
C_{1}\sigma_{1}^{x}\sigma_{1}^{y}\int_{0}^{\Delta\chi^{2}}\exp(-\frac{r^{2}}{2})\mathrm{d}r^{2}/2$
$\displaystyle=$ $\displaystyle w_{1}[1-\exp(-\Delta\chi^{2}/2)].$

This expression is valid so long as $C<C_{0}\triangleq
w_{1}(1-\frac{C_{1}}{C_{2}})=w_{1}-w_{2}\frac{\sigma_{1}^{x}\sigma_{1}^{y}}{\sigma_{2}^{x}\sigma_{2}^{y}}$,
$=w_{1}-w_{2}$; if, however, $C$ is bigger, than we have to include the
contribution from the secondary mode.

$\begin{split}&C=C_{0}+\int_{r_{0}}^{r_{1}}w_{1}\exp(-\frac{r^{2}}{2})r\mathrm{d}r+\int_{0}^{r_{2}}w_{2}\exp(-\frac{r^{2}}{2})r\mathrm{d}r\\\
&C-C_{0}=w_{1}[\frac{C_{2}}{C_{1}}-\exp(-r_{1}^{2}/2)]+w_{2}[1-\exp(-r_{2}^{2}/2)]\\\
&C-w_{1}+w_{2}\frac{\sigma_{1}^{x}\sigma_{1}^{y}}{\sigma_{2}^{x}\sigma_{2}^{y}}=w_{2}\frac{\sigma_{1}^{x}\sigma_{1}^{y}}{\sigma_{2}^{x}\sigma_{2}^{y}}-w_{1}\exp(-r_{1}^{2}/2)\\\
&~{}~{}~{}~{}~{}~{}+w_{2}-w_{1}\frac{\sigma_{2}^{x}\sigma_{2}^{y}}{\sigma_{1}^{x}\sigma_{1}^{y}}\exp(-r_{1}^{2}/2)\\\
&C=1-w_{1}(1+\frac{\sigma_{2}^{x}\sigma_{2}^{y}}{\sigma_{1}^{x}\sigma_{1}^{y}})\exp(-r_{1}^{2}/2)\end{split}$
(10)

In the third line we used the relation that $w_{1}P(r_{1})=w_{2}P(r_{2}))$.

Furthermore, we have

$\begin{split}\exp(-\frac{r_{1}^{2}}{2})&=\frac{1-C}{w_{1}(1+\frac{\sigma_{2}^{x}\sigma_{2}^{y}}{\sigma_{1}^{x}\sigma_{1}^{y}})}\\\
-r_{1}^{2}/2&=\log(1-C)-\log(w_{1})-\log(1+\frac{\sigma_{2}^{x}\sigma_{2}^{y}}{\sigma_{1}^{x}\sigma_{1}^{y}})\\\
\Delta\chi^{2}=r_{1}^{2}&=-2[\log(1-C)-\log(w_{1})\\\
&~{}~{}~{}~{}~{}~{}~{}-\log(1+\frac{\sigma_{2}^{x}\sigma_{2}^{y}}{\sigma_{1}^{x}\sigma_{1}^{y}})]\end{split}$
(11)

We determined $w_{1}=\frac{3}{4}$ and $w_{2}=\frac{1}{4}$, while keeping
$\sigma_{1}^{x}=\sigma_{2}^{x}$ and $\sigma_{1}^{y}=\sigma_{2}^{y}$, and
choosing the C value as $68.27\%,~{}95.45\%~{}\rm{and}~{}99.73\%$. This yields
the corresponding $\Delta\chi^{2}$ values as 3.11, 6.99 and 12.64.

## References

  * Gregory (50 X) P. C. Gregory,  _Bayesian Logical Data Analysis for the Physical Sciences: A Comparative Approach with Mathematica Support_ (Cambridge, Cambridge, 2005 (ISBN: 0-521-84150-X)).
  * Sivia (back) D. S. Sivia,  _Data Analysis: A Bayesian Tutorial_ (Clarendon (Oxford Univ. Press), Oxford, 1996 (ISBN: 0-19-851762-9 or 0-19-851889-7 in paperback)).
  * Rover (2007) C. Rover,  _Bayesian inference on astrophysical binary inspirals based on gravitational-wave measurements_ , Ph.D. thesis, The University of Auckland (2007).
  * Loredo (2013) T. Loredo, Astrostatistical Challenges for the New Astronomy c (2013).
  * Allison and Dunkley (2013) R. Allison and J. Dunkley, Instrumentation and Methods for Astrophysics 12, 13 (2013), 1308.2675 .
  * Veitch and Vecchio (2010) J. Veitch and A. Vecchio, Physical Review D 81, 062003 (2010).
  * Kirkpatrick  _et al._ (1983) S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, SCIENCE 220, 671 (1983).
  * Neal (1996) R. Neal, Statistics and computing 6, 353 (1996).
  * Tierney and Mira (1999) L. Tierney and a. Mira, Statistics in medicine 18, 2507 (1999).
  * Trias  _et al._ (2009) M. Trias, A. Vecchio, and J. Veitch, arXiv preprint arXiv:0904.2207 (2009).
  * Mira (2001) A. Mira, Metron (2001).
  * Press  _et al._ (2007) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical Recipes 3rd Edition: The Art of Scientific Computing_ , 3rd ed. (Cambridge University Press, New York, NY, USA, 2007).
  * Gelman and Rubin (1992) A. Gelman and D. Rubin, Statistical science 7, 457 (1992).
  * Foreman-Mackey  _et al._ (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publications of the Astronomical Society of the Pacific 125, 306 (2013), 1202.3665 .
  * Feroz and Skilling (2013) F. Feroz and J. Skilling (2013) pp. 106–113.
  * Crowder  _et al._ (2006) J. Crowder, N. Cornish, and J. Reddinger, Physical Review D 73, 063011 (2006).
  * Gregory  _et al._ (2007) P. C. Gregory, K. H. Knuth, A. Caticha, J. L. Center, A. Giffin, and C. C. Rodríguez, AIP Conference Proceedings 954, 307 (2007).
  * Feroz  _et al._ (2009) F. Feroz, M. P. Hobson, and M. Bridges, Monthly Notices of the Royal Astronomical Society 398, 1601 (2009).
  * Green (1995) P. Green, Biometrika 82, 711 (1995).
  * Green and Mira (2001) P. Green and A. Mira, Biometrika 88, 1035 (2001).
  * MacQueen (1967) J. B. MacQueen, in  _Proc. of the fifth Berkeley Symposium on Mathematical Statistics and Probability_ , Vol. 1, edited by L. M. L. Cam and J. Neyman (University of California Press, 1967) pp. 281–297.
  * Peskun (1973) P. H. Peskun, Biometrika 60, 607 (1973).


11institutetext: L. T. K. Hien and D. Papadimitriou 22institutetext: 3NLab,
Huawei Belgium Research Center, Leuven 3001, Belgium
22email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# An inertial ADMM for a class of nonconvex composite optimization with
nonlinear coupling constraints

Le Thi Khanh Hien  Dimitri Papadimitriou

###### Abstract

In this paper, we propose an inertial alternating direction method of
multipliers for solving a class of non-convex multi-block optimization
problems with _nonlinear coupling constraints_. Distinctive features of our
proposed method, when compared with other alternating direction methods of
multipliers for solving non-convex problems with nonlinear coupling
constraints, include: (i) we apply the inertial technique to the update of
primal variables and (ii) we apply a non-standard update rule for the
multiplier by scaling the multiplier by a factor before moving along the
ascent direction where a relaxation parameter is allowed. Subsequential
convergence and global convergence are presented for the proposed algorithm.

###### Keywords:

alternating direction methods of multipliers nonlinear coupling constraints
logistic matrix factorization multiblock nonconvex optimization majorization
minimization

##  1 Introduction

A nonconvex nonsmooth composite optimization problem with nonlinear coupling
constraints has the following form

$\begin{split}\operatorname*{minimize}_{x=(x_{1},\ldots,x_{s})\in\mathbb{R}^{n},y\in\mathbb{R}^{m}}&\quad\Theta(x,y):=F(x)+\sum_{i=1}^{s}f_{i}(x_{i})+G(y)\\\
\mbox{subject to}&\quad h(x)+g(y)=0,\end{split}$ (1)

where $h(x)=(h_{1}(x),\ldots,h_{q}(x))$ is a mapping from $\mathbb{R}^{n}$ to
$\mathbb{R}^{q}$, $h_{i}(x)$ for $i=1,\ldots,q$ are continuously
differentiable functions, $g(y)$ is a mapping from $\mathbb{R}^{m}$ to
$\mathbb{R}^{q}$, $F:\mathbb{R}^{n}\to\mathbb{R}$ is a lower semi-continuous
function, $f_{i}:\mathbb{R}^{n}\to\mathbb{R}\cup\\{+\infty\\}$ are proper
lower semi-continuous functions, and $G$ is $L_{G}$-smooth (that is, $\nabla
G$ is $L_{G}$-Lipschitz continuous). Alternating direction methods of
multipliers (ADMM) for solving Problem (1) with _linear_ coupling constraints
(that is when both $h(x)$ and $g(y)$ are affine mappings) have been deeply
studied in the literature, see e.g., Bot2020 ; Hien_iADMM ; HongADMM2020 ;
Wang2019  and the references therein. However, ADMM for solving Problem (1)
with _nonlinear_ coupling constraints have only appeared in Bolte2018_ADMM ,
Cohen2021  and mADMM2022 . The authors in Bolte2018_ADMM  consider Problem (1)
with $s=1$ and propose a universal framework to study global convergence
analysis of _Lagrangian sequences_ (see (Bolte2018_ADMM, , Section 3.3)). To
establish the global convergence of the Lagrangian sequences, the notion of
_information zone_ was introduced in Bolte2018_ADMM  and the boundedness of
the multiplier sequence is the key assumption to fulfill the role of the
information zone. The authors in mADMM2022  develop the idea of information
zone of Bolte2018_ADMM  to propose mADMM - a multiblock alternating direction
method of multipliers - for solving Problem (1) with any $s\geq 1$. Although
the convergence analysis presented in mADMM2022  does not use the property of
the Lagrangian sequence proposed in Bolte2018_ADMM  to establish the global
convergence of its generated sequence, it still relies on the boundedness
assumption of the multiplier sequence (as a consequence of using the
information zone). It is important to note that the boundedness assumption of
the multiplier sequence is difficult to verify in practice. This assumption
can be omitted when $g(y)$ is a linear mapping (but $h(x)$ is not necessary to
be linear), see (mADMM2022, , Section 4). The authors in Cohen2021  consider
Problem (1) with $g(y)$ being a linear mapping and $s=1$, they design a
proximal linearized alternating direction method of multipliers that requires
a backtracking procedure to generate the proximal parameters. The convergence
analysis of Cohen2021  neither requires the boundedness of the multiplier
sequence but the backtracking procedure used in Cohen2021  relies on the
boundedness assumption of the generated sequence (see, (Cohen2021, , Lemma
5.2)) to guarantee the boundedness of the proximal parameters; however,
Cohen2021  does not present a sufficient condition to guarantee the
boundedness of the generated sequence.

We notice that _inertial techniques_ have not been used in Bolte2018_ADMM ,
Cohen2021  and mADMM2022 . Let us provide a very brief review for inertial
techniques. As far as we know, POLYAK1964  is the first paper using an
inertial technique that adds an inertial force, also known as a “momentum”,
$\zeta^{k}(x^{k}-x^{k-1})$ (where $x^{k}$ is the current iterate, $x^{k-1}$ is
the previous iterate, and $\zeta^{k}$ is an extrapolation parameter) to the
gradient direction to accelerate the gradient descent method. Later, Nesterov
proposed his well-known accelerated fast gradient methods in the series of
works Nesterov1983 ; Nesterov1998 ; Nesterov2004 ; Nesterov2005 . Since the
appearance of Nesterov’s accelerated gradient methods, inertial techniques
have been widely applied in convex as well as non-convex problems to
accelerate the convergence of first-order methods, see e.g., Hien_iADMM ;
Hien_ICML2020 ; Titan2020 ; Ochs2019 ; Ochs2014 ; Pock2016 ; Xu2013 ; Xu2017 ;
Zavriev1993  and the references therein.

In this paper, we employ the inertial technique proposed in Titan2020  to
propose iADMMn - an inertial alternating direction method of multipliers for
solving the following composite problem with nonlinear coupling constraints

$\begin{split}\operatorname*{minimize}_{x=(x_{1},\ldots,x_{s})\in\mathbb{R}^{n},y\in\mathbb{R}^{m}}&\quad\Theta(x,y):=F(x)+\sum_{i=1}^{s}f_{i}(x_{i})+G(y)\\\
\mbox{subject to}&\quad h(x)+\mathcal{B}y=0,\end{split}$ (2)

where $\mathcal{B}$ is a linear mapping from $\mathbb{R}^{m}$ to
$\mathbb{R}^{q}$ and we assume the following Assumption 1. This assumption is
useful for us to establish the properties that are necessary for the upcoming
convergence proof when the inertial technique is applied to the update of
block $x_{i}$. In Remark 2, we will also consider a situation when this
assumption can be removed.

###### Assumption 1

For $i=1,\ldots,s$, $x_{i}\mapsto\frac{1}{2}\|h(x)\|^{2}$ is
$l_{i}(x_{1},\ldots,x_{i-1},x_{i+1},\ldots,x_{s})$-smooth.

Let us provide some examples that satisfy Assumption 1.

  * •

If $h(x)=\sum_{i=1}^{s}\mathcal{A}_{i}x_{i}$ ( $h$ is linear), then
$x_{i}\mapsto\frac{1}{2}\|h(x_{i},y_{\neq i})\|^{2}$ is
$\|\mathcal{A}_{i}^{*}\mathcal{A}_{i}\|$-smooth.

  * •

If $h(x_{1},x_{2})=x_{1}x_{2}$ (bi-linear function), then
$x_{1}\mapsto\frac{1}{2}\|h(x_{1},x_{2})\|^{2}$ is
$\|x_{2}x_{2}^{\top}\|$-smooth and
$x_{2}\mapsto\frac{1}{2}\|h(x_{1},x_{2})\|^{2}$ is
$\|x_{1}^{\top}x_{1}\|$-smooth.

  * •

If $h$ is a multilinear mapping then $h$ satisfies Assumption 1.

We also assume the following assumption throughout the paper.

###### Assumption 2

We assume that
$\sigma_{\mathcal{B}}:=\lambda_{\min}(\mathcal{B}\mathcal{B}^{*})>0$ and
$\Theta(x,y)$ is bounded from below by $\nu$.

Problem (2) is Problem (1) with $g(y)=\mathcal{B}y$. Similarly to (mADMM2022,
, Section 4) and Cohen2021 , by considering the linear case of $g$, we would
be able to avoid the boundedness assumption of the multiplier sequence. When
$\mathcal{B}=-\mathcal{I}$, where $\mathcal{I}$ is an identity mapping,
Problem (1) is equivalent to

$\operatorname*{minimize}_{x\in\mathbb{R}^{n}}F(x)+\sum_{i=1}^{s}f_{i}(x_{i})+G(h(x)).$
(3)

An example of (3) is the following logistic matrix factorization problem
Liu2016_PLOS

$\min_{U,V}\sum_{i=1}^{m}\sum_{j=1}^{n}(1+cy_{ij}-y_{ij})\log(1+\exp(u_{i}v_{j}^{\top}))-cy_{ij}u_{i}v_{j}^{\top}+\frac{\lambda_{d}}{2}\|U\|_{F}^{2}+\frac{\lambda_{t}}{2}\|V\|_{F}^{2},$
(4)

where $Y=\mathbb{R}^{m\times n}$ is a given data set with each element
$y_{ij}\in\\{0,1\\}$, $\lambda_{d}$ and $\lambda_{t}$ are regularization
parameters, and $c$ is some given constant. Problem (4) has the form of
Problem (3) with $x=(U,V)$,
$F(x)=\frac{\lambda_{d}}{2}\|U\|_{F}^{2}+\frac{\lambda_{t}}{2}\|V\|_{F}^{2}$,
$f_{i}=0$, $h(U,V)=UV$, and
$G(W)=\sum_{i,j}(1+cy_{ij}-y_{ij})\log(1+\exp(W_{ij}))-cy_{ij}W_{ij}$. A few
other examples of Problem (3) are the PDE-constrained inverse problem
Roosta2014 , the risk parity portfolio selection problem Maillard2010 , the
robust phase retrieval problem Duchi2018 , the nonlinear regression problem
Dutter1981 ; GoodBengCour16  ($h$ represents the model to train, $G$ is a loss
function, and $F(x)+\sum_{i=1}^{s}f_{i}(x_{i})$ is a regularizer), and the
generative adversarial networks mao2017squares .

It is worth noting that inertial techniques have also been applied in
Hien_iADMM  to accelerate the convergence of ADMM for solving Problem (1) with
both $h$ and $g$ being linear mappings. We allow $h$ to be nonlinear in this
paper. On the other hand, iADMMn uses a non-standard update rule for the
multiplier. More specifically, the multiplier is scaled by a factor before
moving along the ascent direction, see (10). This update rule for the
multipliers is inspired by recent papers SunSun2021 ; YANG2022110551  which
give new perspectives on the multiplier update of primal-dual methods in the
non-convex setting. Specifically, classical primal-dual methods are
interpreted as methods that alternatively update the primal variable by
minimizing the primal problem (primal descent) and update the dual (the
multiplier) by maximizing the dual problem (dual ascent), see (Bertsekas, ,
Chapter 7). However, when the primal problem is highly non-convex and may not
be done in closed form, the classic dual ascent step may lose its valid
interpretation. Hence, it makes sense to consider some modifications for the
dual update. For example, SunSun2021  proposes a scaled dual _descent_ update
and YANG2022110551  imposes a discounting factor to the multiplier before
moving along the ascent direction. The multiplier update (10) of our iADMMn is
similar to the multiplier update proposed in YANG2022110551  in the sense that
it scales the multiplier by a factor $\tau_{1}$ before moving along the ascent
direction, but iADMMn also allows a relaxation parameter $\tau_{2}$ in the
step-size of updating the multiplier.

The paper is organized as follows. In the next section, we provide some
preliminaries on block surrogate functions, inertial techniques, and the
augmented Lagrangian function. In Section 3, we describe iADMMn and analyze
its convergence properties. We conclude the paper in Section 4.

#### Notations

We denote $[s]=\\{1,\ldots,s\\}$. For $y\in\mathbb{R}^{s}$, we use $y_{\neq
i}$ to denote $(y_{1},\ldots,y_{i-1},y_{i+1},\ldots,y_{s})$. For a given
mapping $h:\mathbb{R}^{n}\to\mathbb{R}^{s}$, we use $\nabla
h(x)\in\mathbb{R}^{s\times n}$ to denote the Jacobian of $h$ at $x$, that is
$\nabla h(x)=[\nabla h_{1}(x)\ldots\nabla h_{s}(x)]^{\top}$, and we use
$\nabla_{x_{i}}h(x)\in\mathbb{R}^{s\times n_{i}}$ to denote the partial
Jacobian of $h$ with respect to $x_{i}$, that is
$\nabla_{x_{i}}h(x)=[\nabla_{x_{i}}h_{1}(x)\ldots\nabla_{x_{i}}h_{s}(x)]^{\top}$.
For a given sequence $\\{y^{k}\\}$, we denote $\Delta y^{k}=y^{k}-y^{k-1}$.

##  2 Preliminaries

We refer the readers to (Hien_iADMM, , Appendix 1) for some preliminaries of
non-convex optimization such as the definition of limiting subdifferential and
the definition of the KŁ property. In the following, we give some
preliminaries of block surrogate functions, inertial technique, augmented
Lagrangian function, and $\varepsilon$-stationary point.

###  2.1 Block surrogate and nearly sufficiently decrease property

######  Definition 1 (Block surrogate function)

Let $\mathcal{X}_{i}\subseteq\mathbb{R}^{n_{i}}$, for $i=1,\ldots s$, and
$\mathcal{X}\subseteq\mathbb{R}^{n}$. A continuous function
$u_{i}:\mathcal{X}_{i}\times\mathcal{X}\to\mathbb{R}$ is called a block
$x_{i}$ surrogate function of $\varphi:\mathbb{R}^{n}\to\mathbb{R}$ on
$\mathcal{X}=\mathcal{X}_{1}\times\ldots\times\mathcal{X}_{s}$ if the
following conditions are satisfied:

  * (a)

$u_{i}(z_{i},z)=\varphi(z)$ for all $z\in\mathcal{X}$,

  * (b)

$u_{i}(x_{i},z)\geq\varphi(x_{i},z_{\neq i})$ for all
$x_{i}\in\mathcal{X}_{i}$ and $z\in\mathcal{X}$,

where $(x_{i},z_{\neq i})$ denotes
$(z_{1},\ldots,z_{i-1},x_{i},z_{i+1},\ldots,z_{s}).$ The block approximation
error is defined as $e_{i}(x_{i},z):=u_{i}(x_{i},z)-\varphi(x_{i},z_{\neq
i}).$

Given $\varphi$, $u_{i}(x_{i},z)=\varphi(x_{i},z_{\neq
i})+\frac{\rho_{i}}{2}\|x_{i}-z_{i}\|^{2}$, where $\rho_{i}>0$, is a simple
example of a block surrogate function. More examples such as a Lipschitz
gradient surrogate, a quadratic surrogate, a DC programming surrogate, a
saddle point surrogate, and a Jensen surrogate can be found in Hien_iADMM ;
Titan2020 ; Mairal_ICML13 .

The authors in Titan2020  propose TITAN, an inertial block majorization
minimization framework for solving _unconstrained_ non-convex optimization
problems. TITAN updates one block of variables at a time by minimizing an
inertial block surrogate function that is formed by adding an inertial force
to a block surrogate function of the objective. The key property to establish
the convergence of TITAN is the following nearly sufficiently decreasing
property (NSDP), which holds when an inertial block surrogate is minimized
(that is when (5) is performed).

###### Proposition 1

(Titan2020, , Theorem 3) Let $\Phi:\mathbb{R}^{n}\to\mathbb{R}$ be a lower
semi-continuous function and
$f_{i}:\mathbb{R}^{n}\to\mathbb{R}\cup\\{+\infty\\}$ be proper lower semi-
continuous functions. Denote $x^{k,0}=x^{k}$,
$x^{k,i}=(x^{k+1}_{1},\ldots,x^{k+1}_{i},x^{k}_{i+1},\ldots,x^{k}_{s})$ for
$i\in[s]$, $x^{k+1}=x^{k,s}$, and
$\Psi(x):=\Phi(x)+\sum_{i=1}^{s}f_{i}(x_{i})$. Suppose
$\mathcal{G}^{k}_{i}:\mathbb{R}^{n_{i}}\times\mathbb{R}^{n_{i}}\to\mathbb{R}^{n_{i}}$
be some extrapolation operator that satisfies
$\mathcal{G}^{k}_{i}(x^{k}_{i},x^{k-1}_{i})\leq
a_{i}^{k}\|x^{k}_{i}-x^{k-1}_{i}\|$ for some $a_{i}^{k}$. Let $u_{i}(x_{i},z)$
be a block $x_{i}$ surrogate function of $\Phi(x)$ and

$x_{i}^{k+1}=\operatorname*{argmin}_{x_{i}}u_{i}(x_{i},x^{k,i-1})+f_{i}(x_{i})-\langle\mathcal{G}^{k}_{i}(x^{k}_{i},x^{k-1}_{i}),x_{i}\rangle.$
(5)

We assume one of the following conditions holds:

  * •

$x_{i}\mapsto u_{i}(x_{i},z)+f_{i}(x_{i})$ is $\rho_{i}$-strongly convex,

  * •

the approximation error $e_{i}(x_{i},z):=u_{i}(x_{i},z)-\Phi(x_{i},z_{\neq
i})$ satisfying

$e_{i}(x_{i},z)\geq\frac{\rho_{i}(z)}{2}\|x_{i}-z_{i}\|^{2}\mbox{for
all}\,\,x_{i}.$

(Note that $\rho_{i}(z)$ may depend on $z$.) Then we have the following NSDP

$\Psi(x^{k,i-1})+\gamma_{i}^{k}\|x_{i}^{k}-x_{i}^{k-1}\|^{2}\geq\Psi(x^{k,i})+\eta_{i}^{k}\|x_{i}^{k+1}-x_{i}^{k}\|^{2},$
(6)

where

$\begin{array}[]{ll}\gamma^{k}_{i}=\frac{(a^{k}_{i})^{2}}{2\nu\rho_{i}^{k}},\qquad\eta^{k}_{i}=\frac{(1-\nu)\rho_{i}^{k}}{2},\end{array}$

$\rho_{i}^{k}=\rho_{i}(x^{k,i-1})$, and $\nu\in(0,1)$ is a constant. If we do
not apply extrapolation, that is $a_{i}^{k}=0$, then (6) is satisfied with
$\gamma_{i}^{k}=0$ and $\eta_{i}^{k}=\rho^{k}_{i}/2$.

The following proposition, which is derived from (Hien_ICML2020, , Remark 3)
and (Xu2013, , Lemma 2.1), provides an NSDP when $\Phi$ is multi-convex,
Lipschitz gradient suggorates are used for $\Phi$, and $f_{i}$, $i\in[s]$, are
also convex. We then obtain better values of $\gamma_{i}^{k}$ for the NSDP and
larger extrapolation parameters could be used.

###### Proposition 2

Let $\Phi$, $f_{i}$ and $\Psi$ be defined as in Proposition 1. For $i\in[s]$,
suppose $x_{i}\mapsto\Phi(x)$ is an $L_{i}(x_{\neq i})$-smooth convex function
(the Lipschitz constant may depend on $x_{\neq i}$, the values of the other
blocks) and $f_{i}(x_{i})$ is convex. Define
$\hat{x}_{i}^{k}=x_{i}^{k}+\alpha_{i}^{k}(x_{i}^{k}-x_{i}^{k-1})$,
$\bar{x}_{i}^{k}=x_{i}^{k}+\zeta_{i}^{k}(x_{i}^{k}-x_{i}^{k-1})$, and
$\bar{x}^{k,i-1}=(x^{k+1}_{1},\ldots,x^{k+1}_{i-1},\bar{x}^{k}_{i},x^{k}_{i+1},\ldots,x^{k}_{s})$,
where $\alpha_{i}^{k}$ and $\zeta_{i}^{k}$ are extrapolation parameters.
Suppose

$x_{i}^{k+1}=\operatorname*{argmin}_{x_{i}}\langle\nabla_{i}\Phi(\bar{x}^{k,i-1}),x_{i}\rangle+f_{i}(x_{i})+\frac{L_{i}^{k}}{2}\|x_{i}-\hat{x}_{i}^{k}\|^{2},$
(7)

where $L_{i}^{k}=L_{i}(x^{k,i-1}_{\neq i})$. Note that (7) is equivalent to
(5) with

$u_{i}(x_{i},x^{k,i-1})=\Phi(x^{k,i-1})+\langle\nabla_{i}\Phi(x^{k,i-1}),x_{i}-x_{i}^{k}\rangle+\frac{L_{i}^{k}}{2}\|x_{i}-x_{i}^{k}\|^{2},$

and
$\mathcal{G}_{i}^{k}(x_{i}^{k},x_{i}^{k-1})=\nabla_{i}\Phi(\bar{x}^{k,i-1})-\nabla_{i}\Phi(x^{k,i-1})+L_{i}^{k}(x_{i}^{k}-\hat{x}_{i}^{k})$.
Then we have

$\Phi(x^{k,i-1})+f_{i}(x^{k}_{i})+\gamma_{i}^{k}\|x_{i}^{k}-x_{i}^{k-1}\|^{2}\geq\Phi(x^{k,i})+f_{i}(x^{k+1}_{i})+\eta_{i}^{k}\|x^{k+1}_{i}-x_{i}^{k}\|^{2},$

where

$\gamma^{k}_{i}=\frac{L^{k}_{i}}{2}\big{(}(\zeta_{i}^{k})^{2}+\frac{(\gamma_{i}^{k}-\alpha_{i}^{k})^{2}}{\nu}\big{)},\qquad\eta^{k}_{i}=\frac{(1-\nu)L^{k}_{i}}{2}.$

It implies that the NSDP (6) is also satisfied.

If $\alpha_{i}^{k}=\zeta_{i}^{k}$ then we have Inequality (6) is satisfied
with

$\gamma^{k}_{i}=\frac{L^{k}_{i}}{2}(\zeta_{i}^{k})^{2},\qquad\eta^{k}_{i}=\frac{L^{k}_{i}}{2}.$

###  2.2 Augmented Lagrangian and stationary point

Considering Problem (2), its augmented Lagrangian is

$\mathcal{L}_{\beta}(x,y,\omega)=\Theta(x,y)+\langle
h(x)+\mathcal{B}y,\omega\rangle+\frac{\beta}{2}\|h(x)+\mathcal{B}y\|^{2}.$ (8)

###### Definition 2

We call $(x^{*},y^{*})$ a stationary point of Problem (2) if there exists
$\omega^{*}$ such that the following conditions are satisfied.

$0\in\partial_{x_{i}}(f_{i}(x_{i}^{*})+F(x^{*}))+\nabla_{x_{i}}h(x^{*})^{\top}\omega^{*},\,\,\nabla
G(y^{*})+\mathcal{B}^{*}\omega^{*}=0,\,\,h(x^{*})+\mathcal{B}y^{*}=0.$ (9)

We know that finding a stationary point of (2) is equivalent to finding a
critical point of the augmented Lagrangian $\mathcal{L}_{\beta}$, see
(mADMM2022, , Section 2.2). We are also interested in an
$\varepsilon$-stationary point of (2), which is defined as follows.

###### Definition 3

We call $(x^{*},y^{*})$ an $\varepsilon$-stationary point of Problem (2) if
there exists $\omega^{*}$ and
$\chi_{i}\in\partial_{x_{i}}(f_{i}(x_{i}^{*})+F(x^{*}))$ such that

$\|\chi_{i}+\nabla_{x_{i}}h(x^{*})^{\top}\omega^{*}\|\leq\varepsilon,\|\nabla
G(y^{*})+\mathcal{B}^{*}\omega^{*}\|\leq\varepsilon,\|h(x^{*})+\mathcal{B}y^{*}\|\leq\varepsilon.\vspace{-0.05in}$

##  3 Algorithm description and convergence analysis

Before describing iADMMn and presenting its convergence analysis, let us have
a discussion on the proof methodology for a convergence guarantee of an
alternating direction method of multipliers for solving a _non-convex_
optimization problem. An informal general proof recipe that describes the main
steps to prove the global convergence of a first-order method to a critical
point of the objective function of an _unconstrained non-convex_ optimization
problem was introduced in Attouch2013 ; Bolte2014 . The three main steps of
the recipe are (i) sufficient decrease property of the objective sequence,
(ii) a subgradient lower bound for the iterates gap, and (iii) using the KL
property, see (Bolte2014, , Section 3.2) for more details. This proof
methodology has been used wisely to prove the global convergence of the
proximal point algorithm, the forward-backward splitting algorithm, the
regularized Gauss–Seidel method Attouch2009 ; Attouch2013 , the proximal
alternating linearized minimization algorithm Bolte2014 , and the multi-block
Bregman proximal alternating linearized minimization algorithm
ahookhosh2019multi . By replacing the sufficient decrease property of the
objective sequence in Step (i) with the sufficient decrease property of the
augmented Lagrangian sequence or the regularized augmented Lagrangian
sequence, the proof methodology of Attouch2013 ; Bolte2014  is developed to
prove the global convergence of ADMM for solving non-convex problems with
linear/nonlinear coupling constraints, see (Wang2019, , Section 3),
(Bolte2018_ADMM, , Section 3.3), (Bot2020, , Theorem 1), (mADMM2022, , Theorem
2), (Yashitini2022, , Section 3). The augmented Lagrangian function or its
regularized/auxiliary form plays the role of a Lyapunov function (sometimes it
is called a potential function) of an original Lyapunov methodology. The
sufficient decrease property of the Lyapunov function, the subgradient lower
bound for the iterate gaps, together with the KL property of the Lyapunov
function are sufficient to prove the global convergence of the generated
sequence to a critical point of the Lyapunov function. This general proof
recipe is also applied when an inertial technique is used to accelerate the
convergence, see Hien_iADMM ; Hien_ICML2020 ; Titan2020 ; Hien_BMME2022 ;
Ochs2014 ; Xu2013 . Choosing the appropriate Lyapunov function and
establishing its sufficient decrease property is the cornerstone of proving
global convergence. In this paper, we will also use this general proof recipe
to prove the global convergence of iADMMn. To that end, in the following, we
sequentially describe the update of $x_{i}$, $y$, and $\omega$, and establish
the necessary NSDPs. These NSDPs will be used to prove the sufficient decrease
property of an appropriate Lyapunov function that is defined later in (41).

###  3.1 Algorithm description and the NSDPs

Let us remind the notations

$x^{k,0}=x^{k},\,\,x^{k,i}=(x^{k+1}_{1},\ldots,x^{k+1}_{i},x^{k}_{i+1},\ldots,x^{k}_{s})\,\text{
for }i\in[s],\,\text{ and }\,\,x^{k+1}=x^{k,s}.$

We have

$\mathcal{L}_{\beta}(x,y,\omega)=\sum_{i=1}^{s}f_{i}(x_{i})+\varphi(x,y,\omega),$

where

$\begin{split}\varphi(x,y,\omega)&=F(x)+G(y)+\langle
h(x)+\mathcal{B}y,\omega\rangle+\frac{\beta}{2}\|h(x)+\mathcal{B}y\|^{2}\\\
&=F(x)+G(y)+\langle
h(x)+\mathcal{B}y,w\rangle+\frac{\beta}{2}\|h(x)\|^{2}+\frac{\beta}{2}\|\mathcal{B}y\|^{2}+\beta\langle
h(x),\mathcal{B}y\rangle.\end{split}$

Throughout the paper we assume the following assumption for
$\mathcal{L}_{\beta}$.

###### Assumption 3

We have

$\partial\mathcal{L}_{\beta}(x,y,\omega)=\partial_{x_{1}}\mathcal{L}_{\beta}(x,y,\omega)\times\ldots\times\partial_{x_{s}}\mathcal{L}_{\beta}(x,y,\omega)\times\partial_{y}\mathcal{L}_{\beta}(x,y,\omega)\times\partial_{\omega}\mathcal{L}_{\beta}(x,y,\omega).$

This assumption is naturally satisfied if $F$ is a sum of a continuously
differentiable function and a block separable function, see (Attouch2010, ,
Proposition 2.1).

We summarize the description of iADMMn in Algorithm 1. The detailed update
rules of $x_{i}$ and $y$ are elaborated in the following.

Suppose $\hat{u}_{i}$ is a block surrogate function of $F$ with respect to
block $x_{i}$ and $\hat{u}_{i}$ satisfies Assumption 4. Parameters of the
algorithm are chosen as in (29).

Choose initial points $x^{0}$, $y^{0}$, $\omega^{0}$.

for $k=0,\ldots$ do

for $i=1,\ldots,s$ do

Compute the extrapolation point
$\bar{x}_{i}^{k}=x_{i}^{k}+\alpha_{i}^{k}(x_{i}^{k}-x_{i}^{k-1})$.

Under Assumption 1, we update $x_{i}$ as in (15).

If together with Assumption 1 we also have $x_{i}\mapsto F(x_{i},y_{\neq i})$
is $L_{i}(y_{\neq i})$-smooth then we update $x_{i}$ as in (21), and if
Assumption 1 is not satisfied but $x_{i}\mapsto F(x_{i},y_{\neq i})$ is
$L_{i}(y_{\neq i})$-smooth then we update $x_{i}$ as in (25).

end for

Update $y$ as in (27).

Update $\omega$ as
$\omega^{k+1}=\tau_{1}\omega^{k}+\tau_{2}\beta(h(x^{k+1})+\mathcal{B}y^{k+1}).$
(10)

end for

Algorithm 1 iADMMn for solving Problem (2)

#### Update rule for block $x_{i}$ and NSDP when updating $x_{i}$

Following the inertial technique proposed in Titan2020 , we conduct three
steps to update block $x_{i}$: (1) form a block $x_{i}$ surrogate function of
$\varphi(x,y,\omega)$, (2) add an inertial term
$\langle\mathcal{G}_{i}^{k}(x_{i}^{k},x_{i}^{k-1}),x_{i}\rangle$ to this
surrogate function, and (3) update $x_{i}$ by minimizing the sum of
$f_{i}(x_{i})$ and the inertial surrogate function. To form a block $x_{i}$
surrogate function of $\varphi(x,y,\omega)$, we note that if $u_{i}(x_{i},z)$
is a block $x_{i}$ surrogate function of $F(x)+\frac{\beta}{2}\|h(x)\|^{2}$
then

$u_{i}(x_{i},z)+G(y)+\langle h(x_{i},z_{\neq
i})+\mathcal{B}y,w\rangle+\frac{\beta}{2}\|\mathcal{B}y\|^{2}+\beta\langle
h(x_{i},z_{\neq i}),\mathcal{B}y\rangle$ (11)

is a block $x_{i}$ surrogate function of $\varphi(x,y,\omega)$. A general
update rule for $x_{i}$ can be described as follows.

$\begin{split}&x_{i}^{k+1}\in\operatorname*{argmin}_{x_{i}}\Big{\\{}f_{i}(x_{i})+\langle
h(x_{i},x^{k,i-1}_{\neq
i}),\omega^{k}+\beta\mathcal{B}y^{k}\rangle+u_{i}(x_{i},x^{k,i-1})\\\
&\qquad\qquad\qquad\qquad\qquad-\langle\mathcal{G}_{i}^{k}(x_{i}^{k},x_{i}^{k-1}),x_{i}\rangle\Big{\\}}.\end{split}$
(12)

Denote $h_{2}(x)=\frac{1}{2}\|h(x)\|^{2}$. As we assume $x_{i}\mapsto
h_{2}(x_{i},y_{\neq i})$ is $l_{i}(y_{\neq i})$-smooth (see Assumption 1),
then we take the Lipschitz gradient surrogate for $\frac{1}{2}\|h(x)\|^{2}$
and formulate $u_{i}$ of (11) as follows

$u_{i}(x_{i},z)=\hat{u}_{i}(x_{i},z)+\beta
h_{2}(z)+\beta\langle\nabla_{i}h_{2}(z),x_{i}-z_{i}\rangle+\frac{\beta\kappa_{i}(z_{\neq
i})}{2}\|x_{i}-z_{i}\|^{2},$ (13)

where $\kappa_{i}(z_{\neq i})\geq l_{i}(z_{\neq i})$ and $\hat{u}_{i}$ is a
block surrogate function of $F$ with respect to $x_{i}$. We assume the
following assumption for $\hat{u}_{i}$.

###### Assumption 4

We assume that $x_{i}\mapsto\hat{u}_{i}(x_{i},z)-F(x_{i},z_{\neq i})$ is
continuously differentiable and $\nabla_{x_{i}}(\hat{u}_{i}(z_{i},z)-F(z))=0$

After forming the block surrogate for $\varphi$, the second step is to take
the inertial term

$\mathcal{G}_{i}^{k}(x_{i}^{k},x_{i}^{k-1})=\beta\nabla_{i}h_{2}(x_{i}^{k},x^{k,i-1}_{\neq
i})-\beta\nabla_{i}h_{2}(\bar{x}_{i}^{k},x^{k,i-1}_{\neq
i})+\beta\kappa_{i}^{k}(\bar{x}_{i}^{k}-x_{i}^{k}),$ (14)

where $\kappa_{i}^{k}\geq l_{i}(x^{k,i-1}_{\neq i})=l_{i}^{k}$ and
$\bar{x}_{i}^{k}=x_{i}^{k}+\alpha_{i}^{k}(x_{i}^{k}-x_{i}^{k-1})$, here
$\alpha_{i}^{k}$ is an extrapolation parameter.

Then the update (12) of $x_{i}$ is rewritten as follows.

$\begin{split}&x_{i}^{k+1}\in\operatorname*{argmin}_{x_{i}}\Big{\\{}f_{i}(x_{i})+\hat{u}_{i}(x_{i},x^{k,i-1}_{\neq
i})+\langle h(x_{i},x^{k,i-1}_{\neq
i}),\omega^{k}+\beta\mathcal{B}y^{k}\rangle\\\
&\,+\big{\langle}\beta\nabla_{i}h(\bar{x}_{i}^{k},x^{k,i-1}_{\neq
i})^{\top}h(\bar{x}_{i}^{k},x^{k,i-1}_{\neq
i}),x_{i}\big{\rangle}+\frac{\beta\kappa_{i}^{k}}{2}\|x_{i}-\bar{x}_{i}^{k}\|^{2}\Big{\\}}.\end{split}$
(15)

Let us establish the NSDP for $\mathcal{L}_{\beta}$ when the update in (15) is
used. We have

$\|\mathcal{G}_{i}^{k}(x_{i}^{k},x_{i}^{k-1})\|\leq\beta(l_{i}^{k}+\kappa_{i}^{k})\alpha_{i}^{k}\|\Delta
x_{i}^{k}\|.$

The approximation error between $\varphi(x,y,\omega)$ and its block $x_{i}$
surrogate in (11) satisfies that

$\begin{split}&\hat{u}_{i}(x_{i},z)-F(x_{i},z_{\neq i})+\beta
h_{2}(z)+\beta\langle\nabla_{i}h_{2}(z),x_{i}-z_{i}\rangle\\\
&\qquad\qquad\qquad\qquad\qquad+\frac{\kappa_{i}(z_{\neq
i})\beta}{2}\|x_{i}-z_{i}\|^{2}-\beta h_{2}(x_{i},z_{\neq i})\\\
&\geq\theta_{i}(x_{i},z):=\beta h_{2}(z)-\beta h_{2}(x_{i},z_{\neq
i})+\beta\langle\nabla_{i}h_{2}(z),x_{i}-z_{i}\rangle+\frac{\kappa_{i}(z_{\neq
i})\beta}{2}\|x_{i}-z_{i}\|^{2}.\end{split}$

Note that $\nabla^{2}_{x_{i}}\theta_{i}(x_{i},z)=\kappa_{i}(z_{\neq
i})\beta\mathcal{I}-\beta\nabla^{2}_{x_{i}}h_{2}(x_{i},z_{\neq i})$. This
implies $x_{i}\mapsto\theta_{i}(x_{i},z)$ is $\beta(\kappa_{i}(z_{\neq
i})-l_{i}(z_{\neq i}))$-strongly convex (since
$\nabla^{2}_{x_{i}}h_{2}(x_{i},z_{\neq i})\preceq l_{i}(z_{\neq
i})\mathcal{I}$). On the other hand, $\nabla_{x_{i}}\theta_{i}(z_{i},z)=0$. It
follows from (Hien_iADMM, , Lemma 1) that

$\theta_{i}(x_{i},z)\geq\frac{1}{2}\beta(\kappa_{i}(z_{\neq i})-l_{i}(z_{\neq
i}))\|x_{i}-z_{i}\|^{2}.$

We then apply Proposition 1 to obtain the following NSDP.

$\mathcal{L}(x^{k,i},y^{k},\omega^{k})+\eta^{k}_{i}\|\Delta
x^{k+1}_{i}\|^{2}\leq\mathcal{L}(x^{k,i-1},y^{k},\omega^{k})+\gamma_{i}^{k}\|\Delta
x^{k}_{i}\|^{2},$ (16)

where $\kappa_{i}^{k}>l_{i}^{k}$ and

$\eta^{k}_{i}=\frac{(1-\nu_{i})(\kappa^{k}_{i}-l_{i}^{k})\beta}{2},\quad\gamma_{i}^{k}=\frac{(a_{i}^{k})^{2}}{2\nu_{i}(\kappa^{k}_{i}-l_{i}^{k})\beta},\quad
a_{i}^{k}=\beta(l_{i}^{k}+\kappa_{i}^{k})\alpha_{i}^{k}.$ (17)

If $x_{i}\mapsto f_{i}(x_{i})+\hat{u}_{i}(x_{i},x^{k,i-1}_{\neq i})$ is convex
and $x_{i}\mapsto h(x)$ is linear then we take $\kappa_{i}^{k}=l_{i}^{k}$ and
the NSDP in (16) can be tighter. Indeed, as $x_{i}\mapsto h_{2}(x)$ and
$x_{i}\mapsto f_{i}(x_{i})+\hat{u}_{i}(x_{i},x^{k,i-1}_{\neq i})+\langle
h(x_{i},x^{k,i-1}_{\neq i}),\omega^{k}+\beta\mathcal{B}y^{k}\rangle$ are
convex, we apply Proposition 2 to obtain

$\begin{split}&h_{2}(x^{k,i})+f_{i}(x_{i}^{k})+\hat{u}_{i}(x_{i}^{k},x^{k,i-1}_{\neq
i})+\langle h(x_{i}^{k},x^{k,i-1}_{\neq
i}),\omega^{k}+\beta\mathcal{B}y^{k}\rangle+\eta_{i}^{k}\|\Delta
x_{i}^{k+1}\|^{2}\\\ &\leq
h_{2}(x^{k,i-1})+f_{i}(x_{i}^{k-1})+\hat{u}_{i}(x_{i}^{k-1},x^{k,i-1}_{\neq
i})+\langle h(x_{i}^{k-1},x^{k,i-1}_{\neq
i}),\omega^{k}+\beta\mathcal{B}y^{k}\rangle\\\
&\qquad\qquad+\gamma_{i}^{k}\|\Delta x_{i}^{k}\|^{2},\end{split}$ (18)

where $\gamma_{i}^{k}=\frac{1}{2}\beta l_{i}^{k}(\alpha_{i}^{k})^{2}$ and
$\eta_{i}^{k}=\frac{1}{2}\beta l_{i}^{k}$. Furthermore, note that
$\hat{u}_{i}(x_{i}^{k-1},x^{k,i-1}_{\neq i})=F(x^{k,i-1})$ and
$\hat{u}_{i}(x_{i}^{k},x^{k,i-1}_{\neq i})\geq F(x_{i}^{k},x^{k,i-1}_{\neq
i})=F(x^{k,i})$. Therefore, (18) implies that (16) is satisfied with
$\gamma_{i}^{k}=\frac{1}{2}\beta l_{i}^{k}(\alpha_{i}^{k})^{2}$ and
$\eta_{i}^{k}=\frac{1}{2}\beta l_{i}^{k}$. In the following proposition, we
summarize the NSDP when using (15) to update $x_{i}$ .

######  Proposition 3 (NSDP when updating $x_{i}$)

Suppose Assumption 1 is satisfied and $x_{i}$ is updated as in (15). We choose
$\kappa_{i}^{k}>l_{i}^{k}$, then the NSDP (16) is satisfied with
$\eta_{i}^{k}$ and $\gamma_{i}^{k}$ defined in (17). If $x_{i}\mapsto
f_{i}(x_{i})+\hat{u}_{i}(x_{i},x^{k,i-1}_{\neq i})$ is convex and
$x_{i}\mapsto h(x)$ is linear then we take $\kappa_{i}^{k}=l_{i}^{k}$ and then
the NSDP (16) is satisfied with $\gamma_{i}^{k}=\frac{1}{2}\beta
l_{i}^{k}(\alpha_{i}^{k})^{2}$ and $\eta_{i}^{k}=\frac{1}{2}\beta l_{i}^{k}$.

We discuss two other situations in the following remarks.

###### Remark 1

Consider the case that together with Assumption 1 we also have $x_{i}\mapsto
F(x_{i},y_{\neq i})$ is an $L_{i}(y_{\neq i})$-smooth function. Then
$x_{i}\mapsto\hat{h}(x):=F(x)+\frac{\beta}{2}\|h(x)\|^{2}$ is
$\mathbf{l}_{i}(y_{\neq i})=L_{i}(y_{\neq i})+\beta l_{i}(y_{\neq i})$-smooth.
We take the Lipschitz gradient surrogate for
$F(x)+\frac{\beta}{2}\|h(x)\|^{2}$, that is, $u_{i}$ in (11) is

$u_{i}(x_{i},z)=\hat{h}(z)+\langle\nabla_{i}\hat{h}(z),x_{i}-z_{i}\rangle+\frac{\kappa_{i}(z_{\neq
i})}{2}\|x_{i}-z_{i}\|^{2},$ (19)

where $\kappa_{i}(z_{\neq i})\geq\mathbf{l}_{i}(z_{\neq i})$. And we take the
inertial term

$\mathcal{G}_{i}^{k}(x_{i}^{k},x_{i}^{k-1})=\nabla_{i}\hat{h}(x_{i}^{k},x^{k,i-1}_{\neq
i})-\nabla_{i}\hat{h}(\bar{x}_{i}^{k},x^{k,i-1}_{\neq
i})+\kappa_{i}^{k}(\bar{x}_{i}^{k}-x_{i}^{k}),$ (20)

where $\kappa_{i}^{k}\geq\mathbf{l}_{i}^{k}=\mathbf{l}_{i}(x^{k,i-1}_{\neq
i})$. The update (12) of $x_{i}$ is rewritten as follows

$\begin{split}&x_{i}^{k+1}\in\operatorname*{argmin}_{x_{i}}\Big{\\{}f_{i}(x_{i})+\langle
h(x_{i},x^{k,i-1}_{\neq i}),\omega^{k}+\beta\mathcal{B}y^{k}\rangle\\\
&+\big{\langle}\nabla_{i}F(\bar{x}_{i}^{k},x^{k,i-1}_{\neq
i})+\beta\nabla_{i}h(\bar{x}_{i}^{k},x^{k,i-1}_{\neq
i})^{\top}h(\bar{x}_{i}^{k},x^{k,i-1}_{\neq
i}),x_{i}\big{\rangle}+\frac{\kappa_{i}^{k}}{2}\|x_{i}-\bar{x}_{i}^{k}\|^{2}\Big{\\}}.\end{split}$
(21)

Similarly to the above reasoning of Proposition 3, we take
$\kappa_{i}^{k}>\mathbf{l}_{i}^{k}$ and apply Proposition 1 to obtain the NSDP
in (16) with

$\eta^{k}_{i}=\frac{(1-\nu_{i})(\kappa_{i}^{k}-\mathbf{l}_{i}^{k})\beta}{2},\quad\gamma_{i}^{k}=\frac{(a_{i}^{k})^{2}}{2\nu_{i}(\kappa_{i}^{k}-\mathbf{l}_{i}^{k})\beta},\quad
a_{i}^{k}=\beta(\mathbf{l}_{i}^{k}+\kappa_{i}^{k})\alpha_{i}^{k}.$ (22)

If together with the smoothness of $x_{i}\mapsto\hat{h}(x)$ we assume that
$x_{i}\mapsto f_{i}(x_{i})$ and $x_{i}\mapsto F(x)$ are convex and
$x_{i}\mapsto h(x)$ is linear, then we take
$\kappa_{i}^{k}=\mathbf{l}_{i}^{k}$ and (16) holds with
$\gamma_{i}^{k}=\frac{1}{2}\beta\mathbf{l}_{i}^{k}(\alpha_{i}^{k})^{2}$ and
$\eta_{i}^{k}=\frac{1}{2}\beta\mathbf{l}_{i}^{k}$.

###### Remark 2

In this remark, we discuss a case when Assumption 1 can be removed: we assume
$x_{i}\mapsto F(x_{i},y_{\neq i})$ is an $L_{i}(y_{\neq i})$-smooth function
(but $x_{i}\mapsto h_{2}(x_{i},y_{\neq i})$ can be non-smooth), then we take

$u_{i}(x_{i},z)=\beta h_{2}(x_{i},z_{\neq
i})+F(z)+\langle\nabla_{i}F(z),x_{i}-z_{i}\rangle+\frac{\kappa_{i}(z_{\neq
i})}{2}\|x_{i}-z_{i}\|^{2},$ (23)

where $\kappa_{i}(z_{\neq i})\geq L_{i}(z_{\neq i})$. And we take the inertial
term

$\mathcal{G}_{i}^{k}(x_{i}^{k},x_{i}^{k-1})=\nabla_{i}F(x_{i}^{k},x^{k,i-1}_{\neq
i})-\nabla_{i}F(\bar{x}_{i}^{k},x^{k,i-1}_{\neq
i})+\kappa_{i}^{k}(\bar{x}_{i}^{k}-x_{i}^{k}),$ (24)

where $\kappa_{i}^{k}\geq L_{i}(x^{k,i-1}_{\neq i})=L_{i}^{k}$. The update
(12) of $x_{i}$ is rewritten as follows

$\begin{split}&x_{i}^{k+1}\in\operatorname*{argmin}_{x_{i}}\Big{\\{}f_{i}(x_{i})+\langle
h(x_{i},x^{k,i-1}_{\neq
i}),\omega^{k}+\mathcal{B}y^{k}\rangle+\frac{\beta}{2}\|h(x_{i},x^{k,i-1}_{\neq
i})\|^{2}\\\ &\,+\big{\langle}\nabla_{i}F(\bar{x}_{i}^{k},x^{k,i-1}_{\neq
i}),x_{i}\big{\rangle}+\frac{\kappa_{i}^{k}}{2}\|x_{i}-\bar{x}_{i}^{k}\|^{2}\Big{\\}}.\end{split}$
(25)

Similarly to the reasoning of Proposition 3, we take
$\kappa_{i}^{k}>L_{i}^{k}$ and apply Proposition 1 to obtain the NSDP in (16)
with

$\eta^{k}_{i}=\frac{(1-\nu_{i})(\kappa_{i}^{k}-L_{i}^{k})\beta}{2},\quad\gamma_{i}^{k}=\frac{(a_{i}^{k})^{2}}{2\nu_{i}(\kappa_{i}^{k}-L_{i}^{k})\beta},\quad
a_{i}^{k}=\beta(L_{i}^{k}+\kappa_{i}^{k})\alpha_{i}^{k}.$ (26)

If together with the smoothness of $x_{i}\mapsto F(x)$ we assume that
$x_{i}\mapsto f_{i}(x_{i})$ and $x_{i}\mapsto F(x)$ are convex and
$x_{i}\mapsto h(x)$ is linear, then we take $\kappa_{i}^{k}=L_{i}^{k}$ and
(16) holds with $\gamma_{i}^{k}=\frac{1}{2}\beta
L_{i}^{k}(\alpha_{i}^{k})^{2}$ and $\eta_{i}^{k}=\frac{1}{2}\beta L_{i}^{k}$.

#### Update $y$

As $G$ is $L_{G}$-smooth, we form a block $y$ surrogate of $\varphi$ by
summing the Lipschitz gradient surrogate of $G$ and the remaining part of
$\varphi$:

$\begin{split}\hat{\varphi}(y,x,y^{\prime},\omega)&=G(y^{\prime})+\langle\nabla
G(y^{\prime}),y-y^{\prime}\rangle+\frac{L_{G}}{2}\|y-y^{\prime}\|^{2}\\\
&\qquad+F(x)+\langle
h(x)+\mathcal{B}y,\omega\rangle+\frac{\beta}{2}\|h(x)+\mathcal{B}y\|^{2}.\end{split}$

Block $y$ is updated as follows

$\begin{split}y^{k+1}&\in\operatorname*{argmin}_{y}\hat{\varphi}(y,x^{k+1},y^{k},\omega^{k})\\\
&=\operatorname*{argmin}_{y}\langle\mathcal{B}^{*}\omega^{k}+\nabla
G(y^{k}),y\rangle+\frac{\beta}{2}\|h(x^{k+1})+\mathcal{B}y\|^{2}+\frac{L_{G}}{2}\|y-y^{k}\|^{2},\\\
&=(\beta\mathcal{B}^{*}\mathcal{B}+L_{G}\mathbf{I})^{-1}\big{(}L_{G}y^{k}-\nabla
G(y^{k})-\mathcal{B}^{*}(\omega^{k}+\beta h(x^{k+1}))\big{)}.\end{split}$ (27)

######  Proposition 4 (Sufficient decrease when updating $y$)

The update in (27) guarantees a sufficient decrease

$\mathcal{L}(x^{k+1},y^{k+1},\omega^{k})+\frac{\delta}{2}\|y^{k+1}-y^{k}\|^{2}\leq\mathcal{L}(x^{k+1},y^{k},\omega^{k}),$
(28)

where $\delta=L_{G}+\beta\lambda_{\min}(\mathcal{B}^{*}\mathcal{B})$. If $G$
is convex, then (28) is satisfied with $\delta=L_{G}$.

###### Proof

We note that $y\mapsto\hat{\varphi}(y,x,y^{\prime},\omega)$ is
$L_{G}+\beta\lambda_{\min}(\mathcal{B}^{*}\mathcal{B})$-strongly convex.
Applying Proposition 1 (note that we do not use inertial term for updating
$y$), we obtain (28). If $G$ is convex, then we apply Proposition 2 (we simply
take $\alpha_{i}^{k}=\beta_{i}^{k}=0$ in Proposition 2) to obtain the result.

###  3.2 Convergence analysis

Given $\tau_{1}$, $\tau_{2}$, we denote
$C_{1}=\frac{(\tau_{1}+1)|\tau_{1}-\tau_{2}|}{2\sigma_{\mathcal{B}}\tau_{2}\beta(1-|\tau_{1}-\tau_{2}|)}$,
$C_{2}=\frac{(\tau_{1}+1)\tau_{2}/\tau_{1}}{2\sigma_{\mathcal{B}}\beta(1-|\tau_{1}-\tau_{2}|)(1-|1-\tau_{2}/\tau_{1}|)}$,
$C_{3}=\frac{\delta}{2}-2C_{2}L_{G}^{2}$, where $\delta$ is defined in
Proposition 4. We remind that $\gamma_{i}^{k}$ and $\eta_{i}^{k}$ are the
coefficients in the NSDP(16). Their formulas are defined in Proposition 3 if
we assume Assumption 1 and use the update rule in (15). If we use the update
rule in (21) then the formulas of $\gamma_{i}^{k}$ and $\eta_{i}^{k}$ are
defined in Remark 1. And if we use the update rule in (25) then the formulas
of $\gamma_{i}^{k}$ and $\eta_{i}^{k}$ are defined in Remark 2.

We need to choose the parameters such that they satisfy the following
conditions.

$\begin{split}&\tau_{1}\in(0,1],\quad\tau_{2}/\tau_{1}\in(0,2),\quad|\tau_{1}-\tau_{2}|<1\\\
&\gamma_{i}^{k}\leq B_{1}\eta_{i}^{k},\quad 8C_{2}L_{G}^{2}\leq
B_{2}C_{3},\end{split}$ (29)

where $B_{1},B_{2}\in(0,1)$ are some constants. In this section, we will
establish the subsequential convergence and the global convergence for iADMMn.
The following proposition provides a recursive inequality for
$\\{\mathcal{L}(x^{k},y^{k};\omega^{k})\\}$.

###### Proposition 5

The following inequality holds

$\begin{split}&\mathcal{L}(x^{k+1},y^{k+1},\omega^{k+1})-\frac{1-\tau_{1}}{2\tau_{2}\beta}\|\omega^{k+1}\|^{2}+\sum_{i=1}^{s}\eta_{i}^{k}\|\Delta
x_{i}^{k+1}\|^{2}+C_{3}\|\Delta y^{k+1}\|^{2}\\\
&\leq\mathcal{L}(x^{k},y^{k},\omega^{k})-\frac{1-\tau_{1}}{2\tau_{2}\beta}\|\omega^{k}\|^{2}+\sum_{i=1}^{s}\gamma_{i}^{k}\|\Delta
x_{i}^{k}\|^{2}\\\
&\qquad+C_{1}\big{(}\|\mathcal{B}^{*}\Delta\omega^{k}\|^{2}-\|\mathcal{B}^{*}\Delta\omega^{k+1}\|^{2}\big{)}+8C_{2}L_{G}^{2}\|\Delta
y^{k}\|^{2}.\end{split}$

###### Proof

Optimality condition of (27) gives us

$\mathcal{B}^{*}(\omega^{k}+\beta(h(x^{k+1})+\mathcal{B}y^{k+1}))+\nabla
G(y^{k})+L_{G}(y^{k+1}-y^{k})=0$ (30)

Denote $t^{k+1}=L_{G}\Delta y^{k+1}+\nabla G(y^{k})$. From (30) and the update
of $\omega$ in (10) we have

$\displaystyle\mathcal{B}^{*}((1-\frac{\tau_{1}}{\tau_{2}})\omega^{k}+\frac{1}{\tau_{2}}\omega^{k+1})=-t^{k+1},$

which implies that

$\frac{1}{\tau_{2}}\mathcal{B}^{*}\Delta\omega^{k+1}+(1-\frac{\tau_{1}}{\tau_{2}})\mathcal{B}^{*}\Delta\omega^{k}=-\Delta
t^{k+1}.$

Hence, we have

$\frac{1}{\tau_{1}}\mathcal{B}^{*}\Delta\omega^{k+1}=(1-\frac{\tau_{2}}{\tau_{1}})\mathcal{B}^{*}\Delta\omega^{k}-\frac{\tau_{2}}{\tau_{1}}\Delta
t^{k+1}.$ (31)

Note that if $\tau_{2}/\tau_{1}\in[1,2)$ then (31) can be rewritten as

$\frac{1}{\tau_{1}}\mathcal{B}^{*}\Delta\omega^{k+1}=-(\tau_{2}/\tau_{1}-1)\mathcal{B}^{*}\Delta\omega^{k}-(2-\tau_{2}/\tau_{1})\frac{\tau_{2}/\tau_{1}}{2-\tau_{2}/\tau_{1}}\Delta
t^{k+1}.$

Hence, from $\tau_{2}/\tau_{1}\in(0,2)$ and the convexity of the norm
$\|\cdot\|$, we can derive from (31) the following inequality

$\frac{1}{\tau_{1}}\|\mathcal{B}^{*}\Delta\omega^{k+1}\|^{2}\leq|1-\tau_{2}/\tau_{1}|\|\mathcal{B}^{*}\Delta\omega^{k}\|^{2}+\frac{(\tau_{2}/\tau_{1})^{2}}{1-|1-\tau_{2}/\tau_{1}|}\|\Delta
t^{k+1}\|^{2}.$ (32)

On the other hand, from the definition of $t^{k+1}$ and the Lipschitz
continuity of $\nabla G$, we have $\|\Delta t^{k+1}\|\leq L_{G}\|\Delta
y^{k+1}\|+L_{G}\|\Delta y^{k}\|+L_{G}\|\Delta y^{k}\|$. It implies that

$\begin{split}\|\Delta t^{k+1}\|^{2}\leq 2L_{G}^{2}\|\Delta
y^{k+1}\|^{2}+8L_{G}^{2}\|\Delta y^{k}\|^{2}.\end{split}$

Hence, from (32) we obtain

$\begin{split}\|\mathcal{B}^{*}\Delta\omega^{k+1}\|^{2}&\leq\frac{|\tau_{1}-\tau_{2}|}{(1-|\tau_{1}-\tau_{2}|)}\big{(}\|\mathcal{B}^{*}\Delta\omega^{k}\|^{2}-\|\mathcal{B}^{*}\Delta\omega^{k+1}\|^{2}\big{)}\\\
&\quad+\frac{\tau_{2}^{2}/\tau_{1}}{(1-|\tau_{1}-\tau_{2}|)(1-|1-\tau_{2}/\tau_{1}|)}\big{(}2L_{G}^{2}\|\Delta
y^{k+1}\|^{2}+8L_{G}^{2}\|\Delta y^{k}\|^{2}\big{)}.\end{split}$ (33)

Furthermore,

$\begin{split}&\mathcal{L}(x^{k+1},y^{k+1},\omega^{k+1})-\mathcal{L}(x^{k+1},y^{k+1},\omega^{k})\\\
&=\big{\langle}h(x^{k+1})+\mathcal{B}y^{k+1},\omega^{k+1}-\omega^{k}\big{\rangle}\\\
&=\big{\langle}\frac{1}{\tau_{2}\beta}(\omega^{k+1}-\tau_{1}\omega^{k}),\omega^{k+1}-\omega^{k}\big{\rangle}\\\
&=\frac{1}{\tau_{2}\beta}\big{\langle}\tau_{1}(\omega^{k+1}-\omega^{k})+(1-\tau_{1})\omega^{k+1},\omega^{k+1}-\omega^{k}\big{\rangle}\\\
&=\frac{\tau_{1}}{\tau_{2}\beta}\|\Delta\omega^{k+1}\|^{2}+\frac{1-\tau_{1}}{2\tau_{2}\beta}\big{(}\|\Delta\omega^{k+1}\|^{2}+\|\omega^{k+1}\|^{2}-\|\omega^{k}\|^{2}\big{)},\end{split}$

which leads to

$\begin{split}&\mathcal{L}(x^{k+1},y^{k+1},\omega^{k+1})-\mathcal{L}(x^{k+1},y^{k+1},\omega^{k})\\\
&=\frac{\tau_{1}+1}{2\tau_{2}\beta}\|\Delta\omega^{k+1}\|^{2}+\frac{1-\tau_{1}}{2\tau_{2}\beta}\big{(}\|\omega^{k+1}\|^{2}-\|\omega^{k}\|^{2}\big{)}.\end{split}$
(34)

Therefore, from (33), (34), and noting that
$\sigma_{\mathcal{B}}\|\Delta\omega^{k+1}\|^{2}\leq\|\mathcal{B}^{*}\Delta\omega^{k+1}\|^{2}$,
we have

$\begin{split}&\mathcal{L}(x^{k+1},y^{k+1},\omega^{k+1})-\mathcal{L}(x^{k+1},y^{k+1},\omega^{k})\\\
&\leq
C_{1}\big{(}\|\mathcal{B}^{*}\Delta\omega^{k}\|^{2}-\|\mathcal{B}^{*}\Delta\omega^{k+1}\|^{2}\big{)}+C_{2}2L_{G}^{2}\|\Delta
y^{k+1}\|^{2}+C_{2}8L_{G}^{2}\|\Delta y^{k}\|^{2}\\\
&\qquad\qquad+\frac{1-\tau_{1}}{2\tau_{2}\beta}\big{(}\|\omega^{k+1}\|^{2}-\|\omega^{k}\|^{2}\big{)}.\end{split}$

Together with (16) and (28) we obtain the result.

###### Lemma 1

Denote

$\hat{\mathcal{L}}^{k}=\mathcal{L}(x^{k},y^{k},\omega^{k})-\frac{1-\tau_{1}}{2\tau_{2}\beta}\|\omega^{k}\|^{2}+\sum_{i=1}^{s}B_{1}\eta_{i}^{k}\|\Delta
x_{i}^{k}\|^{2}+C_{1}\|\mathcal{B}^{*}\Delta\omega^{k}\|^{2}+B_{2}C_{3}\|\Delta
y^{k}\|^{2}.$ (35)

We have $\hat{\mathcal{L}}^{k}\geq\nu$ for all $k\geq 1$, where $\nu$ is the
lower bound of $\Theta(x,y)$ (see Assumption 2).

###### Proof

We use the technique of (melo2017, , Lemma 2.9). From Proposition 5 and the
conditions $\gamma_{i}^{k}\leq B_{1}\eta_{i}^{k}$ and $8C_{2}L_{G}^{2}\leq
B_{2}C_{3}$ for some constants $B_{1},B_{2}\in(0,1)$, we have
$\hat{\mathcal{L}}^{k+1}\leq\hat{\mathcal{L}}^{k}$, that is,
$\\{\hat{\mathcal{L}}^{k}\\}$ is a non-increasing sequence. Suppose there
exists $k_{0}$ such that $\hat{\mathcal{L}}^{k}<\nu$ for all $k\geq k_{0}$,
then we have

$\sum_{k=1}^{K}(\hat{\mathcal{L}}^{k}-\vartheta)\leq\sum_{k=1}^{k_{0}}(\hat{\mathcal{L}}^{k}-\vartheta)+(K-k_{0})(\hat{\mathcal{L}}^{k}-\vartheta).$

This implies that
$\sum_{k=1}^{\infty}(\hat{\mathcal{L}}^{k}-\vartheta)=-\infty$. However, since

$\begin{split}\hat{\mathcal{L}}^{k}&\geq\mathcal{L}(x^{k},y^{k},\omega^{k})-\frac{1-\tau_{1}}{2\tau_{2}\beta}\|\omega^{k}\|^{2}\\\
&\geq\nu+\langle
h(x^{k})+\mathcal{B}y^{k},\omega^{k}\rangle-\frac{1-\tau_{1}}{2\tau_{2}\beta}\|\omega^{k}\|^{2}\\\
&=\nu+\frac{1}{\tau_{2}\beta}\langle\omega^{k}-\tau_{1}\omega^{k-1},\omega^{k}\rangle-\frac{1-\tau_{1}}{2\tau_{2}\beta}\|\omega^{k}\|^{2}\\\
&=\nu+\frac{1+\tau_{1}}{2\tau_{2}\beta}\|\omega^{k}\|^{2}-\frac{\tau_{1}}{\tau_{2}\beta}\langle\omega^{k-1},\omega^{k}\rangle\\\
&=\nu+\frac{1+\tau_{1}}{2\tau_{2}\beta}\|\omega^{k}\|^{2}+\frac{\tau_{1}}{2\tau_{2}\beta}(\|\Delta\omega^{k}\|^{2}-\|\omega^{k}\|^{2}-\|\omega^{k-1}\|^{2})\\\
&\geq\nu+\frac{1}{2\tau_{2}\beta}\|\omega^{k}\|^{2}-\frac{\tau_{1}}{2\tau_{2}\beta}\|\omega^{k-1}\|^{2}\\\
&\geq\nu+\frac{\tau_{1}}{2\tau_{2}\beta}(\|\omega^{k}\|^{2}-|\omega^{k-1}\|^{2}),\end{split}$

we have

$\sum_{k=1}^{K}(\hat{\mathcal{L}}^{k}-\vartheta)\geq\frac{\tau_{1}}{2\tau_{2}\beta}\sum_{k=1}^{K}(\|\omega^{k}\|^{2}-\|\omega^{k-1}\|^{2})\geq\frac{\tau_{1}}{2\tau_{2}\beta}(-\|\omega^{0}\|^{2}),$

which gives a contradiction to the fact
$\sum_{k=1}^{\infty}(\hat{\mathcal{L}}^{k}-\vartheta)=-\infty$.

###### Proposition 6

Suppose there exist $\tilde{\eta}_{i}>0$ such that
$\eta_{i}^{k}\geq\tilde{\eta}_{i}$ for all $i\in[s]$ and $k\geq 1$. Then the
sequence $\\{\Delta y^{k}\\}_{k\geq 0}$, $\\{\Delta x_{i}^{k}\\}_{k\geq 0}$,
and $\\{\Delta\omega^{k}\\}_{k\geq 0}$ converge to 0.

###### Proof

From Proposition 5 we have

$\hat{\mathcal{L}}^{k+1}+(1-B_{1})\sum_{i=1}^{s}\eta_{i}^{k}\|\Delta
x_{i}^{k+1}\|^{2}+(1-B_{2})C_{3}\|\Delta
y^{k+1}\|^{2}\leq\hat{\mathcal{L}}^{k},$ (36)

where $\hat{\mathcal{L}}^{k}$ is defined in (35). For all $K\geq 1$, by
summing (36) from $k=0$ to $K$ we get

$\hat{\mathcal{L}}^{K+1}+(1-B_{1})\sum_{k=0}^{K}\sum_{i=1}^{s}\eta_{i}^{k}\|\Delta
x_{i}^{k+1}\|^{2}+(1-B_{2})C_{3}\sum_{k=0}^{K}\|\Delta
y^{k+1}\|^{2}\leq\hat{\mathcal{L}}^{0}.$ (37)

From (37), Lemma 1 and the assumption $\eta_{i}^{k}\geq\tilde{\eta}_{i}>0$ we
derive that the sequence $\\{\Delta y^{k}\\}_{k\geq 0}$ and $\\{\Delta
x_{i}^{k}\\}_{k\geq 0}$ converge to 0.

From (33), we have

$\begin{split}&\sum_{k=0}^{K}\|\mathcal{B}^{*}\Delta\omega^{k+1}\|^{2}\leq\frac{|\tau_{1}-\tau_{2}|}{(1-|\tau_{1}-\tau_{2}|)}\|\mathcal{B}^{*}\Delta\omega^{0}\|^{2}\\\
&\quad+\frac{\tau_{2}^{2}/\tau_{1}}{(1-|\tau_{1}-\tau_{2}|)(1-|1-\tau_{2}/\tau_{1}|)}\big{(}2L_{G}^{2}\sum_{k=0}^{K}\|\Delta
y^{k+1}\|^{2}+8L_{G}^{2}\sum_{k=0}^{K}\|\Delta y^{k}\|^{2}\big{)}.\end{split}$

It implies that
$\sum_{k=0}^{\infty}\|\Delta\omega^{k}\|^{2}\leq\sum_{k=0}^{\infty}1/\sigma_{\mathcal{B}}\|\mathcal{B}^{*}\Delta\omega^{k}\|^{2}<+\infty$.
Hence $\\{\Delta\omega^{k}\\}_{k\geq 0}$ converge to 0.

######  Theorem 3.1 (Subsequential convergence)

Suppose there exist $\tilde{\eta}_{i}>0$ such that
$\eta_{i}^{k}\geq\tilde{\eta}_{i}$ for all $i\in[s]$ and $k\geq 1$. And
suppose $a_{i}^{k}$ in (17) (if we use use the update rule in (15)), or
$a_{i}^{k}$ in (22) (if we use use the update rule in (21)), or $a_{i}^{k}$ in
(26) (if we use use the update rule in (25)) are upper bounded111In fact, if
the Lipschitz constants of $x_{i}\mapsto h_{2}(x)$ in Assumption 1, or that of
$x_{i}\mapsto\hat{h}(x)$ in Remark 1, or that of $x_{i}\mapsto F(x)$ in Remark
2 are upper bounded over the set that contains $\\{x^{k}\\}_{k\geq 0}$ then
this assumption is satisfied. Note that $x_{i}^{k}\in{\rm{dom}\,}(f_{i})$.. It
holds that if $(x^{*},y^{*},\omega^{*})$ is a limit point of the generated
sequence of iADMMn then $(x^{*},y^{*},\omega^{*})$ is an
$\frac{1-\tau_{1}}{\tau_{2}\beta}\|\omega^{*}\|$-approximate stationary point
of Problem (2). When $\tau_{1}=1$, we have $(x^{*},y^{*})$ is a stationary
point of Problem (2) (or equivalently, $(x^{*},y^{*},\omega^{*})$ is a
critical point of $\mathcal{L}_{\beta}$).

###### Proof

We remind that the update rule of $x_{i}^{k}$ in (15), (21) and (25) are just
special cases of (12). In (12), $u_{i}(x_{i},z)$ is a block $x_{i}$ surrogate
function of $F(x)+\frac{\beta}{2}\|h(x)\|^{2}$ and the formulas of
$\mathcal{G}$ in (14), (20), and (24), which respectively correspond to the
update in (15), (21), and (25), satisfy the condition
$\|\mathcal{G}_{i}^{k}(x_{i}^{k},x_{i}^{k-1})\|\leq a_{i}^{k}\|\Delta
x_{i}^{k}\|$.

From (12), we know that the following inequality holds for all $x_{i}$

$\begin{split}&f_{i}(x_{i}^{k+1})+\langle h(x^{k+1}_{i},x^{k,i-1}_{\neq
i}),\omega^{k}+\beta\mathcal{B}y^{k}\rangle+u_{i}(x_{i}^{k+1},x^{k,i-1})\\\
&\leq f_{i}(x_{i})+\langle h(x_{i},x^{k,i-1}_{\neq
i}),\omega^{k}+\beta\mathcal{B}y^{k}\rangle+u_{i}(x_{i},x^{k,i-1})\\\
&\qquad\qquad\qquad-\langle\mathcal{G}_{i}^{k}(x_{i}^{k},x_{i}^{k-1}),x_{i}-x_{i}^{k+1}\rangle.\end{split}$
(38)

Let $\\{(x^{k_{n}},y^{k_{n}},\omega^{k_{n}})\\}$ be a subsequence that
converges to $(x^{*},y^{*},\omega^{*})$. As $\Delta x^{k}_{i}\to 0$ (see
Proposition 6), we have $x_{i}^{k_{n}+1}$, $x_{i}^{k_{n}-1}$ and
$x_{i}^{k_{n}-2}$ also converge to $x_{i}^{*}$ for all $i\in[s]$.
Consequently, $\|\mathcal{G}_{i}^{k}(x_{i}^{k_{n}-1},x_{i}^{k_{n}-2})\|\to 0$.
Taking $x_{i}=x_{i}^{*}$ and $k=k_{n}-1$ in (38), we obtain
$\limsup_{n\to\infty}f_{i}(x_{i}^{k_{n}})\leq f_{i}(x_{i}^{*}),$ which implies
$f_{i}(x_{i}^{k_{n}})\to f_{i}(x_{i}^{*})$ as $f_{i}$ is lower semi-
continuous. Taking $k=k_{n}-1$ in (38) and let $k_{n}\to\infty$, we obtain the
following inequality for all $x_{i}$

$\begin{split}&f_{i}(x_{i}^{*})+\langle
h(x^{*}),\omega^{*}+\beta\mathcal{B}y^{*}\rangle+u_{i}(x_{i}^{*},x^{*})\\\
&\leq f_{i}(x_{i})+\langle h(x_{i},x^{*}_{\neq
i}),\omega^{*}+\beta\mathcal{B}y^{*}\rangle+u_{i}(x_{i},x^{*}).\end{split}$
(39)

Note that $u_{i}(x_{i}^{*},x^{*})=F(x^{*})+\frac{\beta}{2}\|h(x^{*})\|^{2}$.
Hence, (39) implies that for all $x_{i}$ we have

$\begin{split}\varphi(x^{*},y^{*},\omega^{*})+f_{i}(x_{i}^{*})&\leq\varphi(x_{i},x^{*}_{\neq
i},y^{*},\omega^{*})+f_{i}(x_{i})\\\
&\qquad+u_{i}(x_{i},x^{*})-F(x_{i},x^{*}_{\neq
i})-\frac{\beta}{2}\|h(x_{i},x^{*}_{\neq i})\|^{2}.\end{split}$ (40)

Let $e_{i}(x_{i},z)=u_{i}(x_{i},z)-F(x_{i},z_{\neq
i})-\frac{\beta}{2}\|h(x_{i},z_{\neq i})\|^{2}$. We have
$e_{i}(x_{i}^{*},x^{*})=0$. It follows from (40) that $x_{i}^{*}$ is a
solution of

$\min_{x_{i}}\varphi(x_{i},x^{*}_{\neq
i},y^{*},\omega^{*})+f_{i}(x_{i})+e_{i}(x_{i},x^{*}).$

Writing the optimality condition for this problem and noting that
$\nabla_{x_{i}}e_{i}(x^{*}_{i},x^{*})=0$ for all the three surrogates in (13),
(19), and (23) that respectively corresponds to the update rule of $x_{i}$ in
(15), (21) and (25), we obtain

$0\in\partial_{x_{i}}\mathcal{L}(x^{*},y^{*},\omega^{*})=\partial_{x_{i}}(\varphi(x^{*},y^{*},\omega^{*})+f_{i}(x_{i}^{*})).$

Since the essence of the update rule of $y$ is to minimize a block $y$
surrogate function of $\varphi$ (see (27)), we can use the same reasoning for
$x_{i}$ to prove that $0\in\nabla_{y}L(x^{*},y^{*},\omega^{*})$.

From (10) we have
$h(x^{k_{n}})+\mathcal{B}y^{k_{n}}=\frac{1}{\tau_{2}\beta}(\Delta\omega^{k_{n}}+(1-\tau_{1})\omega^{{k_{n}}-1})$.
Furthermore, $\Delta\omega^{k_{n}}\to 0$ (see Proposition 6), which implies
$\omega^{{k_{n}}-1}\to\omega^{*}$. Hence, we have

$h(x^{*})+\mathcal{B}y^{*}=\frac{1-\tau_{1}}{\tau_{2}\beta}\omega^{*}.$

The result follows then.

As discussed at the beginning of the section, we will use the same methodology
employed in Hien_iADMM ; Hien_ICML2020 ; Titan2020 ; Hien_BMME2022 ; Ochs2014
; Xu2013  to prove the global convergence of iADMMn. Although the general
methodology is the same in these papers as its essence is originally from the
general proof recipe of Attouch2013 ; Bolte2014 , how to choose a Lyapunov
function and how to establish the two properties (i) the sufficient decrease
and (ii) the boundedness of subgradient are different in these papers since
these steps highly rely on the structure of the problem.

######  Theorem 3.2 (Global convergence)

Suppose the conditions for Theorem 3.1 are satisfied and we also assume that
$\Theta(x,y)$ has the KŁ property, the generated sequence
$\\{(x^{k},y^{k},\omega^{k})\\}$ is bounded and let $\tau_{1}=\tau_{2}=1$.
Furthermore, together with the existence of the constants $\tilde{\eta}_{i}$
in Theorem 3.1 we assume there exist $\bar{\eta}_{i}>0$ such that
$\eta_{i}^{k}\leq\bar{\eta}_{i}$ for all $i\in[s]$ and $k\geq 0$, and the
constant $B_{1}$ in (29) satisfies
$B_{1}<\min\\{\tilde{\eta}_{i}/\bar{\eta}_{i}\\}$. If we use the surrogate in
(13), which corresponds to the update rule in (15), then we also assume the
following assumption 5

###### Assumption 5

(a) We have
$\partial_{x_{i}}(f_{i}(x_{i})+F(x))=\partial_{x_{i}}f_{i}(x_{i})+\partial_{x_{i}}F(x)$
and
$\partial_{x_{i}}(f_{i}(x_{i})+\hat{u}_{i}(x_{i},z))=\partial_{x_{i}}f_{i}(x_{i})+\partial_{x_{i}}\hat{u}_{i}(x_{i},z)$
for any $x_{i}\in{\rm{dom}\,}(f_{i})$ and $x,z\in\mathbb{R}^{n}$, and

(b) for any $S_{i}^{k}\in\partial_{x_{i}}\hat{u}_{i}(x_{i}^{k+1},x^{k,i-1})$
there exists $\bar{S}_{i}^{k}\in\partial_{x_{i}}F(x^{k+1})$ such that
$\|S_{i}^{k}-\bar{S}_{i}^{k}\|\leq\bar{L}_{i}\|x^{k+1}-x^{k,i-1}\|$ for some
constant $\bar{L}_{i}$.

Then the whole sequence $\\{(x^{k},y^{k},\omega^{k})\\}$ converges to a
critical point of $\bar{\mathcal{L}}$.

Before proving Theorem 3.2, let us have some discussion on the assumptions
used for Theorem 3.2.

  * •

We remark that $F$ is not necessary to be smooth if we use the surrogate in
(13), which corresponds to the update rule in (15). See (mADMM2022, , Appendix
D) for an example that still satisfies Assumption 5 (b) but $F$ is nonsmooth.
Furthermore, the constants $\bar{L}_{i}$ do not affect the way to choose
parameters in (29), their existence is just for the purpose of the upcoming
proof for the boundedness of subgradients.

  * •

If $h(x)=\mathcal{A}x=\sum_{i=1}^{s}\mathcal{A}_{i}x_{i}$, where
$\mathcal{A}_{i}$ are linear mappings, and we use the surrogate in (13) then
$l_{i}^{k}=\|\mathcal{A}_{i}^{*}\mathcal{A}_{i}\|$. In this case, we simply
choose $\kappa_{i}^{k}$ in the update (15) to be any constant
$\kappa_{i}\geq\|\mathcal{A}_{i}^{*}\mathcal{A}_{i}\|$, then $\eta_{i}^{k}$ in
Proposition 3 are constants. Hence, $\tilde{\eta}_{i}/\bar{\eta}_{i}=1$.

  * •

It is important to note that we require $\tau_{1}=\tau_{2}=1$ for the global
convergence, but this condition is not required for Theorem 3.1.

  * •

The boundedness assumption of $\\{(x^{k},y^{k},\omega^{k})\\}$ is necessary
for the upcoming proof. Totally similarly to (mADMM2022, , Proposition 7), we
can prove that if $\tau_{1}=1$, $ran\,h(x)\,\subseteq Im(\mathcal{B})$,
$\lambda_{\min}(\mathcal{B}^{*}\mathcal{B})>0$, and $\Theta(x,y)$ is coercive
over the feasible set $\\{(x,y):h(x)+\mathcal{B}y=0\\}$ then the generated
sequence $\\{(x^{k},y^{k},\omega^{k})\\}$ is bounded. We omit the proof of
this property.

  * •

Once the global convergence of iADMMn is guaranteed, by using the same
technique as in (Attouch2009, , Theorem 2) (see (Xu2013, , Theorem 2.9) and
(Hien_ICML2020, , Theorem 3) for some examples of using this technique to
establish the convergence rate), we can establish a convergence rate for
iADMMn. The type of convergence rate depends on the value of the KŁ exponent.
We omit the details, see (mADMM2022, , Theorem 2) or (Hien_iADMM, , Section
2.2) for a similar discussion.

We now prove Theorem 3.2.

###### Proof

Denote $z=(x,y,\omega)$. We use the following Lyapunov function

$\begin{split}\bar{\mathcal{L}}(z,\tilde{z})&=\mathcal{L}_{\beta}(x,y,\omega)+\sum_{i=1}^{s}\frac{\tilde{\eta}_{i}+B_{1}\bar{\eta}_{i}}{2}\|x_{i}-\tilde{x}_{i}\|^{2}\\\
&\qquad\qquad+\frac{(1+B_{2})C_{3}}{2}\|y-\tilde{y}\|^{2}+C_{1}\|\mathcal{B}^{*}(\omega-\tilde{\omega})\|^{2}.\end{split}$
(41)

We will prove the two properties (i) the sufficient decrease and (ii) the
boundedness of subgradient. The other remaining steps to obtain the global
convergence are totally similar to the proof of (Hien_iADMM, , Theorem 2).

#### Sufficient decrease property

From Proposition 5 and the assumption $\tilde{\eta}_{i}\leq\eta_{i}^{k}$,
$\gamma_{i}^{k}\leq B_{1}\eta_{i}^{k}\leq B_{1}\bar{\eta}_{i}$, and
$8C_{2}L_{G}^{2}\leq B_{2}C_{3}$, we obtain

$\begin{split}&\bar{\mathcal{L}}(z^{k+1},z^{k})+\sum_{i=1}^{s}\frac{\tilde{\eta}_{i}-B_{1}\bar{\eta}_{i}}{2}\big{(}\|\Delta
x_{i}^{k+1}\|^{2}+\|\Delta x_{i}^{k}\|^{2}\big{)}\\\
&\qquad\qquad+\frac{(1-B_{2})C_{3}}{2}\big{(}\|\Delta y^{k+1}\|^{2}+\|\Delta
y^{k}\|^{2}\big{)}\leq\bar{\mathcal{L}}(z^{k},z^{k-1}).\end{split}$

#### Boundedness of subgradient

In the following, we will work on the bounded set that contains the generated
sequence (as we assume the generated sequence is bounded). The optimality
condition of (12) gives us

$\begin{split}&\mathcal{G}_{i}^{k}(x_{i}^{k},x_{i}^{k-1})-\nabla_{x_{i}}h(x_{i}^{k+1},x^{k,i-1}_{\neq
i})^{\top}(\omega^{k}+\beta\mathcal{B}y^{k})\\\
&\qquad\in\partial_{x_{i}}(f_{i}(x_{i}^{k+1})+u_{i}(x^{k+1}_{i},x^{k,i-1})).\end{split}$
(42)

On the other hand, if we use the surrogate in (13) and assume Assumption 5 is
satisfied then we have
$\partial_{x_{i}}(f_{i}(x_{i})+u_{i}(x_{i},z))=\partial_{x_{i}}f_{i}(x_{i})+\partial_{x_{i}}u_{i}(x_{i},z)$
for any $x_{i}\in{\rm{dom}\,}(f_{i})$ and $x,z\in\mathbb{R}^{n}$; furthermore,
it is not difficult to show that there exist some constant $\hat{L}_{i}$ such
that for any $\xi_{i}^{k+1}\in\partial_{x_{i}}u_{i}(x_{i}^{k+1},x^{k,i-1})$
there exists $\bar{\xi}_{i}^{k+1}\in\partial_{x_{i}}\big{(}F(x^{k+1})+\beta
h_{2}(x^{k+1})\big{)}$ such that
$\|\xi_{i}^{k+1}-\bar{\xi}_{i}^{k+1}\|\leq\hat{L}_{i}\|x^{k+1}-x^{k,i-1}\|$.
Note that if the surrogate in (19) of (23) is used then these two properties
are always satisfied. Hence, from (42) we know that there exists
$d_{i}^{k+1}\in\partial_{x_{i}}f_{i}(x_{i}^{k+1})$ and
$\xi_{i}^{k+1}\in\partial_{x_{i}}u_{i}(x^{k+1}_{i},x^{k,i-1})$ such that

$\mathcal{G}_{i}^{k}(x_{i}^{k},x_{i}^{k-1})-\nabla_{x_{i}}h(x_{i}^{k+1},x^{k,i-1}_{\neq
i})^{\top}(\omega^{k}+\beta\mathcal{B}y^{k})=d_{i}^{k+1}+\xi_{i}^{k+1}.$ (43)

Since $h_{2}$ is continuously differentiable and

$\begin{split}\partial_{x_{i}}\mathcal{L}_{\beta}(z^{k+1})&=\partial_{x_{i}}(f_{i}(x_{i}^{k+1})+F(x^{k+1})+\beta
h_{2}(x^{k+1}))\\\
&\qquad\qquad+\nabla_{x_{i}}h(x^{k+1})^{\top}(\omega^{k+1}+\beta\mathcal{B}y^{k+1}),\end{split}$

we have
$D_{i}^{k+1}:=d_{i}^{k+1}+\bar{\xi}_{i}^{k+1}+\nabla_{x_{i}}h(x^{k+1})^{\top}(\omega^{k+1}+\beta\mathcal{B}y^{k+1})\in\partial_{x_{i}}\mathcal{L}_{\beta}(z^{k+1})$.
Together with (43) and
$\|\xi_{i}^{k+1}-\bar{\xi}_{i}^{k+1}\|\leq\hat{L}_{i}\|x^{k+1}-x^{k,i-1}\|$ we
get

$\begin{split}\|D_{i}^{k+1}\|&=\|d_{i}^{k+1}+\xi_{i}^{k+1}+\bar{\xi}_{i}^{k+1}-\xi_{i}^{k+1}+\nabla_{x_{i}}h(x^{k+1})^{\top}(\omega^{k+1}+\beta\mathcal{B}y^{k+1})\|\\\
&\leq\big{\|}\mathcal{G}_{i}^{k}(x_{i}^{k},x_{i}^{k-1})-\nabla_{x_{i}}h(x_{i}^{k+1},x^{k,i-1}_{\neq
i})^{\top}(\omega^{k}+\beta\mathcal{B}y^{k})\\\
&\qquad\qquad+\nabla_{x_{i}}h(x^{k+1})^{\top}(\omega^{k+1}+\beta\mathcal{B}y^{k+1})\big{\|}+\|\bar{\xi}_{i}^{k+1}-\xi_{i}^{k+1}\|\\\
&\leq\|\mathcal{G}_{i}^{k}(x_{i}^{k},x_{i}^{k-1})\|+\|\nabla_{x_{i}}h(x^{k+1})^{\top}\|(\|\Delta\omega^{k+1}\|+\beta\|\mathcal{B}\Delta
y^{k+1}\|)\\\ &\qquad\qquad+\hat{L}_{i}\|x^{k+1}-x^{k,i-1}\|.\end{split}$ (44)

On the other hand, as $h_{i}$, $i=1,\ldots,q$, are continuously
differentiable, hence $\|\nabla_{x_{i}}h(x^{k+1})^{\top}\|$ is bounded on the
bounded set that contains the generated sequence. Since $\tau_{1}=\tau_{2}=1$,
we get from (33) that

$\|\Delta\omega^{k+1}\|^{2}\leq\frac{1}{\sigma_{\mathcal{B}}}\|\mathcal{B}^{*}\Delta\omega^{k+1}\|^{2}\leq\frac{1}{\sigma_{\mathcal{B}}}(2L_{G}^{2}\|\Delta
y^{k+1}\|^{2}+8L_{G}^{2}\|\Delta y^{k}\|^{2}).$ (45)

Therefore, from (44) we derive that

$\|D_{i}^{k+1}\|\leq a_{1}(\|\Delta x^{k+1}\|+\|\Delta x^{k}\|+\|\Delta
y^{k+1}\|+\|\Delta y^{k}\|)$ (46)

for some positive constant $a_{1}$. Note that

$D_{y}^{k+1}:=\nabla
G(y^{k+1})+\mathcal{B}^{*}(\omega^{k+1}+\beta(h(x^{k+1})+\mathcal{B}y^{k+1}))\in\nabla_{y}\mathcal{L}_{\beta}(z^{k+1}),$

which together with(30) and (45) leads to

$\begin{split}\|D_{y}^{k+1}\|&=\|\nabla G(y^{k+1})-\nabla
G(y^{k})+\mathcal{B}^{*}\Delta\omega^{k+1}-L_{G}\Delta y^{k+1}\|\\\ &\leq
2L_{G}\|\Delta y^{k+1}\|+\|\mathcal{B}^{*}\Delta\omega^{k+1}\|\\\ &\leq
a_{2}(\|\Delta y^{k+1}\|+\|\Delta y^{k}\|)\end{split}$ (47)

for some positive constant $a_{2}$. On the other hand, we have

$\|\nabla_{\omega}\mathcal{L}_{\beta}(z^{k+1})\|=\|h(x^{k+1})+\mathcal{B}y^{k+1}\|=1/\beta\|\Delta\omega^{k+1}\|,$
(48)

and

$\begin{split}\partial\bar{\mathcal{L}}(z,\tilde{z})=\partial\mathcal{L}_{\beta}(z)+&\partial\Big{(}\sum_{i=1}^{s}\frac{\tilde{\eta}_{i}+B_{1}\bar{\eta}_{i}}{2}\|x_{i}-\tilde{x}_{i}\|^{2}+B_{2}C_{3}\|y-\tilde{y}\|^{2}\\\
&\qquad+C_{1}\|\mathcal{B}^{*}(\omega-\tilde{\omega})\|^{2}\Big{)}.\end{split}$
(49)

Therefore, from (49), (48), (45), (47), and (44), it is not difficult to
derive that

$\|D^{k+1}\|\leq a_{3}(\|\Delta x^{k+1}\|+\|\Delta x^{k}\|+\|\Delta
y^{k+1}\|+\|\Delta y^{k}\|),$

for some positive constant $a_{3}$ and
$D^{k+1}\in\partial\bar{\mathcal{L}}(z^{k+1},z^{k})$.

##  4 Conclusion

We have analyzed iADMMn, an inertial alternating direction method of
multipliers for solving Problem (2). In essence, iADMMn is an extension of
iADMM proposed in Hien_iADMM : iADMMn allows the nonlinearity of the coupling
constraint whereas iADMM only considers linear coupling constraint. iADMMn is
also an extension of (mADMM2022, , Algorithm 2): iADMMn allows inertial term
in the updating of $x_{i}$ whereas (mADMM2022, , Algorithm 2) does not. On the
other hand, iADMMn considers a non-standard update rule for the multiplier
$\omega$ when it allows a scaling factor $\tau_{1}\in(0,1]$ for $\omega^{k}$
whereas both iADMM and (mADMM2022, , Algorithm 2) use standard rule
$\tau_{1}=1$. Theorem 3.1 shows that every limit point of the generated
sequence is an $\varepsilon$-stationary point of Problem (2). In Theorem 3.2,
$\tau_{1}=1$ is required to guarantee a global convergence for the generated
sequence. However, when $\tau_{1}\neq 1$, if there exist many limit points of
the generated sequence or not is still an open question for us. We take it as
a future research direction.

## Data availability Statement

The manuscript has no associated data.

## References

  * (1) Ahookhosh, M., Hien, L.T.K., Gillis, N., Patrinos, P.: Multi-block Bregman proximal alternating linearized minimization and its application to sparse orthogonal nonnegative matrix factorization.  Computational Optimization and Applications 79 (2021)
  * (2) Attouch, H., Bolte, J.: On the convergence of the proximal algorithm for nonsmooth functions involving analytic features.  Mathematical Programming 116(1), 5–16 (2009)
  * (3) Attouch, H., Bolte, J., Redont, P., Soubeyran, A.: Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Łojasiewicz inequality.  Mathematics of Operations Research 35(2), 438–457 (2010)
  * (4) Attouch, H., Bolte, J., Svaiter, B.F.: Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward–backward splitting, and regularized gauss–seidel methods.  Mathematical Programming 137(1), 91–129 (2013)
  * (5) Bertsekas, D.P.: Nonlinear Programming.  Athena Scientific (2016)
  * (6) Bolte, J., Sabach, S., Teboulle, M.: Proximal alternating linearized minimization for nonconvex and nonsmooth problems.  Mathematical Programming 146(1), 459–494 (2014)
  * (7) Bolte, J., Sabach, S., Teboulle, M.: Nonconvex Lagrangian-based optimization: Monitoring schemes and global convergence.  Mathematics of Operations Research 43(4), 1210–1232 (2018).  DOI 10.1287/moor.2017.0900
  * (8) Bot, R.I., Nguyen, D.K.: The proximal alternating direction method of multipliers in the nonconvex setting: Convergence analysis and rates.  Mathematics of Operations Research 45(2), 682–712 (2020)
  * (9) Cohen, E., Hallak, N., Teboulle, M.: A dynamic alternating direction of multipliers for nonconvex minimization with nonlinear functional equality constraints.  Journal of Optimization Theory and Applications 193, 324–353 (2022)
  * (10) Duchi, J.C., Ruan, F.: Solving (most) of a set of quadratic equalities: composite optimization for robust phase retrieval.  Information and Inference: A Journal of the IMA 8(3), 471–529 (2018)
  * (11) Dutter, R., Huber, P.J.: Numerical methods for the nonlinear robust regression problem.  Journal of Statistical Computation and Simulation 13(2), 79–113 (1981).  DOI 10.1080/00949658108810482.  URL https://doi.org/10.1080/00949658108810482
  * (12) Goodfellow, I.J., Bengio, Y., Courville, A.: Deep Learning.  MIT Press, Cambridge, MA, USA (2016).  http://www.deeplearningbook.org
  * (13) Hien, L., Phan, D., Gillis, N.: Inertial alternating direction method of multipliers for non-convex non-smooth optimization.  Computational Optimization and Applications 83, 247–285 (2022)
  * (14) Hien, L.T.K., Gillis, N., Patrinos, P.: Inertial block proximal method for non-convex non-smooth optimization.  In: Thirty-seventh International Conference on Machine Learning ICML 2020 (2020)
  * (15) Hien, L.T.K., Papadimitriou, D.: Multiblock ADMM for nonsmooth nonconvex optimization with nonlinear coupling constraints (2022).  ArXiv:2201.07657
  * (16) Hien, L.T.K., Phan, D.N., Gillis, N.: Inertial block majorization minimization framework for nonconvex nonsmooth optimization (2020).  ArXiv:2010.12133, accepted in Journal of Machine Learning Research
  * (17) Hien, L.T.K., Phan, D.N., Gillis, N., Ahookhosh, M., Patrinos, P.: Block bregman majorization minimization with extrapolation.  SIAM Journal on Mathematics of Data Science 4(1), 1–25 (2022).  DOI 10.1137/21M1432661
  * (18) Hong, M., Chang, T.H., Wang, X., Razaviyayn, M., Ma, S., Luo, Z.Q.: A block successive upper-bound minimization method of multipliers for linearly constrained convex optimization.  Mathematics of Operations Research 45(3), 833–861 (2020)
  * (19) Liu, Y., Wu, M., Miao, C., Zhao, P., Li, X.L.: Neighborhood regularized logistic matrix factorization for drug-target interaction prediction.  PLOS Computational Biology 12(2), 1–26 (2016).  DOI 10.1371/journal.pcbi.1004760.  URL https://doi.org/10.1371/journal.pcbi.1004760
  * (20) Maillard, S., Roncalli, T., Teiletche, J.: The properties of equally weighted risk contribution portfolios.  The Journal of Portfolio Management 36(4), 60–70 (2010)
  * (21) Mairal, J.: Optimization with first-order surrogate functions.  In: Proceedings of the 30th International Conference on International Conference on Machine Learning - Volume 28, ICML’13, pp. 783–791. JMLR.org (2013)
  * (22) Mao, X., Li, Q., Xie, H., Lau, R.Y.K., Wang, Z., Smolley, S.P.: Least squares generative adversarial networks.  In: Proceedings of the IEEE International Conference on Computer Vision (2017)
  * (23) Melo, J.G., Monteiro, R.D.C.: Iteration-complexity of a Jacobi-type non-Euclidean ADMM for multi-block linearly constrained nonconvex programs (2017)
  * (24) Nesterov, Y.: A method of solving a convex programming problem with convergence rate O$(1/k^{2})$.  Soviet Mathematics Doklady 27(2) (1983)
  * (25) Nesterov, Y.: On an approach to the construction of optimal methods of minimization of smooth convex functions.  Ekonom. i. Mat. Metody 24, 509–517 (1998)
  * (26) Nesterov, Y.: Introductory lectures on convex optimization: A basic course.  Kluwer Academic Publ. (2004)
  * (27) Nesterov, Y.: Smooth minimization of non-smooth functions.  Math. Prog. 103(1), 127–152 (2005)
  * (28) Ochs, P.: Unifying abstract inexact convergence theorems and block coordinate variable metric ipiano.  SIAM Journal on Optimization 29(1), 541–570 (2019)
  * (29) Ochs, P., Chen, Y., Brox, T., Pock, T.: iPiano: Inertial proximal algorithm for nonconvex optimization.  SIAM Journal on Imaging Sciences 7(2), 1388–1419 (2014)
  * (30) Pock, T., Sabach, S.: Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems.  SIAM Journal on Imaging Sciences 9(4), 1756–1787 (2016)
  * (31) Polyak, B.: Some methods of speeding up the convergence of iteration methods.  USSR Computational Mathematics and Mathematical Physics 4(5), 1 – 17 (1964).  DOI https://doi.org/10.1016/0041-5553(64)90137-5.  URL http://www.sciencedirect.com/science/article/pii/0041555364901375
  * (32) Roosta-Khorasani, F., van den Doel, K., Ascher, U.: Stochastic algorithms for inverse problems involving pdes and many measurements.  SIAM Journal on Scientific Computing 36(5), S3–S22 (2014).  DOI 10.1137/130922756
  * (33) Sun, K., Sun, A.: Dual descent alm and admm (2021).  URL https://arxiv.org/abs/2109.13214
  * (34) Wang, Y., Yin, W., Zeng, J.: Global convergence of admm in nonconvex nonsmooth optimization.  Journal of Scientific Computing 78, 29–63 (2019).  DOI 10.1007/s10915-018-0757-z
  * (35) Xu, Y., Yin, W.: A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion.  SIAM Journal on Imaging Sciences 6(3), 1758–1789 (2013).  URL https://doi.org/10.1137/120887795
  * (36) Xu, Y., Yin, W.: A globally convergent algorithm for nonconvex optimization based on block coordinate update.  Journal of Scientific Computing 72(2), 700–734 (2017)
  * (37) Yang, Y., Jia, Q.S., Xu, Z., Guan, X., Spanos, C.J.: Proximal ADMM for nonconvex and nonsmooth optimization.  Automatica 146, 110,551 (2022)
  * (38) Yashtini, M.: Convergence and rate analysis of a proximal linearized ADMM for nonconvex nonsmooth optimization.  Journal of Global Optimization 84, 913–939 (2022)
  * (39) Zavriev, S., Kostyuk, F.: Heavy-ball method in nonconvex optimization problems.  Computational Mathematics and Modeling (1993)