# Correct order on some certain weighted representation functions

Shi–Qiang Chen & Yuchen Ding* & Xiaodong Lü & Yuhan Zhang  (Shi–Qiang Chen)
School of Mathematics and Statistics, Anhui Normal University, Wuhu 241002,
People’s Republic of China<EMAIL_ADDRESS>(Yuchen Ding) School of
Mathematical Sciences, Yangzhou University, Yangzhou 225002, People’s Republic
of China<EMAIL_ADDRESS>(Xiaodong Lü) School of Mathematical Sciences,
Yangzhou University, Yangzhou 225002, People’s Republic of China
<EMAIL_ADDRESS>(Yuhan Zhang) School of Mathematical Sciences, Yangzhou
University, Yangzhou 225002, People’s Republic of China
<EMAIL_ADDRESS>
###### Abstract.

Let $\mathbb{N}$ be the set of all nonnegative integers. For any positive
integer $k$ and any subset $A$ of nonnegative integers, let $r_{1,k}(A,n)$ be
the number of solutions $(a_{1},a_{2})$ to the equation $n=a_{1}+ka_{2}$. In
2016, Qu proved that

$\liminf_{n\rightarrow\infty}r_{1,k}(A,n)=\infty$

providing that $r_{1,k}(A,n)=r_{1,k}(\mathbb{N}\setminus A,n)$ for all
sufficiently large integers, which answered affirmatively a 2012 problem of
Yang and Chen. In a very recent article, another Chen (the first named author)
slightly improved Qu’s result and obtained that

$\liminf_{n\rightarrow\infty}\frac{r_{1,k}(A,n)}{\log n}>0.$

In this note, we further improve the lower bound on $r_{1,k}(A,n)$ by showing
that

$\liminf_{n\rightarrow\infty}\frac{r_{1,k}(A,n)}{n}>0.$

Our bound reflects the correct order of magnitude of the representation
function $r_{1,k}(A,n)$ under the above restrictions due to the trivial fact
that $r_{1,k}(A,n)\leq n/k.$

###### Key words and phrases:

representation functions; order of functions; partitions of integers

###### 2010 Mathematics Subject Classification:

Primary 11B34, Secondary 11A41

*Corresponding author

##  1\. Introduction

Let $\mathbb{N}$ be the set of all nonnegative integers and $A$ a subset of
$\mathbb{N}$. For any nonnegative integer $n$, let $R_{1}(A,n);R_{2}(A,n)$ and
$R_{3}(A,n)$ be the number of solutions $(a,a^{\prime})$ to the equations
$n=a+a^{\prime}$ with $a,a^{\prime}\in A$; $a,a^{\prime}\in A,~{}a<a^{\prime}$
and $a,a^{\prime}\in A,~{}a\leq a^{\prime}$, respectively. For backgrounds on
these representation functions $R_{i}(A,n)$, $i=1,2,3$, one can refer to an
early survey article of Sárközy and Sós [10]. Following Sárközy’s question,
Dombi [5], Chen and Wang [4], Lev [6], Sándor [9], Tang [11], Chen and Tang
[3] and Chen [2] investigated various properties on values of the
representation functions $R_{i}(A,n)$ and $R_{i}(\mathbb{N}\setminus
A,n),i=1,2,3$.

In an interesting paper, Yang and Chen [13] introduced the following weighted
representation function

$r_{k_{1},k_{2}}(A,n)=\\#\left\\{(a_{1},a_{2})\in
A^{2}:n=k_{1}a_{1}+k_{2}a_{2}\right\\},$

where $A$ is a subset of $\mathbb{N}$ and $k_{1},k_{2}$ are two positive
integers. They determined all pairs $(k_{1},k_{2})$ of positive integers for
which there exists a set $A\subseteq\mathbb{N}$ such that

$r_{k_{1},k_{2}}(A,n)=r_{k_{1},k_{2}}\left(\mathbb{N}\setminus A,n\right)$

for all sufficiently large integers, which would reduce to partial answers to
the original question of Sárközy mentioned above for $k_{1}=k_{2}=1$ on
$R_{1}(A,n)$. For $1\leq k_{1}<k_{2}$ with $(k_{1},k_{2})=1$, if there exists
a set $A\subseteq\mathbb{N}$ such that

$r_{k_{1},k_{2}}(A,n)=r_{k_{1},k_{2}}\left(\mathbb{N}\setminus A,n\right)$

for all sufficiently large integers, then Yang and Chen proved that $k_{1}=1$.
So the studies of the weighted representation function $r_{1,k}(A,n)$ would be
of particular interest.

For a positive integer $k\geq 2$, let $\Psi_{k}$ be the set of all
$A\subseteq\mathbb{N}$ such that

$r_{1,k}(A,n)=r_{1,k}\left(\mathbb{N}\setminus A,n\right)$

for sufficiently large integers $n$. A result of Yang [12] states that if
$k,\ell$ are multiplicatively independent (equivalently, $\log k/\log\ell$ is
irrational), then $\Psi_{k}\cap\Psi_{\ell}=\emptyset$. Qu [8] then gave a
complete criteria for which $\Psi_{k}\cap\Psi_{\ell}=\emptyset$. It turns out
to be that $\Psi_{k}\cap\Psi_{\ell}\neq\emptyset$ if and only if $\log
k/\log\ell=a/b$ for some odd positive integers $a$ and $b$, which disproved a
conjecture of Yang [12]. For related result, see also the article of Li and Ma
[7]. In [13], Yang and Chen posed the following problem:

######  Problem 1.1.

For any set $A\in\Psi_{k}$, is it true that $r_{1,k}(A,n)\geq 1$ for all
sufficiently large integers $n$? Is it true that
$r_{1,k}(A,n)\rightarrow\infty$ as $n\rightarrow\infty$?

Problem 1.1 was later answered affirmatively by Qu [8]. Very recently, Chen
[1] improved Qu’s result by showing that

$\liminf_{n\rightarrow\infty}\frac{r_{1,k}(A,n)}{\log n}>0$

for any $A\in\Psi_{k}$. In this note, we give the following very much stronger
bound.

######  Theorem 1.1.

Let $k\geq 2$ be a given integer. For any set $A\in\Psi_{k}$, we have

$\liminf_{n\rightarrow\infty}\frac{r_{1,k}(A,n)}{n}>0.$

######  Remark 1.1.

It can be seen that our new bound is sharp on the order of magnitude of
$r_{1,k}(A,n)$ in the sense that

$\limsup_{n\rightarrow\infty}\frac{r_{1,k}(A,n)}{n}\leq 1/k.$

Perhaps, it should be pointed out that the original argument of Qu [8] with
necessary adjustments can also lead to the bound given by Chen [1]. It is also
worth mentioning that our argument here leading to the sharp bound above in
Theorem 1.1 is different and simplified comparing with the somewhat
complicated ones taken by Qu and Chen.

##  2\. Proofs

Following Qu [8], we may write $A$ as the following union of the ‘blocks’

$A=\bigcup_{i=0}^{\infty}\left[t_{2i},t_{2i+1}\right),$

where $0\leq t_{0}<t_{1}<t_{2}<\cdot\cdot\cdot$ is an increasing sequence of
integers. The proof of our theorem is based on the following lemma of Qu [8,
Lemma 2.1].

######  Lemma 2.2.

Let $k\geq 2$ be a given integer. For any $A\in\Psi_{k}$ with

$A=\bigcup_{i=0}^{\infty}\left[t_{2i},t_{2i+1}\right),$

there exist an odd positive integer $a$ and a nonnegative integer $i_{0}$ such
that $t_{i+a}=kt_{i}$ for all $i\geq i_{0}$.

###### Proof of theorem 1.1.

By Lemma 2.2 for $A\in\Psi_{k}$ with

$A=\bigcup_{i=0}^{\infty}\left[t_{2i},t_{2i+1}\right),$

there exists an odd positive integer $a$ such that $t_{i+a}=kt_{i}$ for all
$i\geq i_{0}$. Without loss of generality, we can assume that $i_{0}=0$,
otherwise one can consider $\widetilde{A}=A\setminus[0,t_{i_{0}})$ instead of
$A$ (this can be seen from the proofs below). Let $T=4(t_{a+2}-t_{0}).$ Then
there exists some odd integer $g\in\mathbb{N}$ such that $k^{g}>T$.

From now on, let $n$ be a sufficiently large number. It is clear that there
are nonnegative integers $m$ and $r$ with $0\leq r<(k^{g}+1)$ such that

$n=(k^{g}+1)m+r.$

We can assume that $m\in\left[k^{s}t_{\ell},k^{s}t_{\ell+1}\right)$ for two
nonnegative integers $s$ and $\ell$ with

$0\leq\ell\leq a-1.$

Recall that the integer $n$ is assumed to be sufficiently large, it follows
that both $m$ and $s$ are sufficiently large. We will prove that

$r_{1,k}(A,n)\geq\frac{n}{k^{5}t_{a}(k^{g}+2)}-(k^{g}+1),$

from which it follows clearly that

$\liminf_{n\rightarrow\infty}\frac{r_{1,k}(A,n)}{n}>0.$

Since $r_{1,k}(A,n)=r_{1,k}(\mathbb{N}\setminus A,n)$ for all sufficiently
large $n$, we can also assume, without loss of generality, that

$\left[k^{s}t_{\ell},k^{s}t_{\ell+1}\right)\subseteq A.$

Since $g$ is an odd integer, we make the observation that each of the
following three intervals

$\left[k^{s}t_{\ell-2},k^{s}t_{\ell-1}\right),\quad\left[k^{s}t_{\ell+2},k^{s}t_{\ell+3}\right)\quad\text{and}\quad\left[k^{s+g-1}t_{\ell},k^{s+g-1}t_{\ell+1}\right)$

contains in $A$ as well, which is crucial in the following arguments. Before
the continuation of the proof, we make the following notice that for brevity
we write, for example,

$k^{3}t_{2},\quad k^{4}t_{2-a}\quad\text{and}\quad k^{2}t_{2+a}$

as the same number at different occasions. The proofs are divided into three
cases:

Case I. $k^{s}t_{\ell}+k^{s-4}\leq m<k^{s}t_{\ell+1}-k^{s-4}.$ Noting that for
any $0\leq q<k^{s-5}-r$, we have

$n=(m+kq+r)+k\left(k^{g-1}m-q\right),$

where both $m+kq+r$ and $k^{g-1}m-q$ belong to $A$ since

$m+kq+r\in\left[k^{s}t_{\ell},k^{s}t_{\ell+1}\right)\quad\text{and}\quad
k^{g-1}m-q\in\left[k^{s+g-1}t_{\ell},k^{s+g-1}t_{\ell+1}\right).$

Thus, we deduce that

$\displaystyle r_{1,k}(A,n)$ $\displaystyle\geq k^{s-5}-r$
$\displaystyle>\frac{m}{k^{5}t_{\ell+1}}-(k^{g}+1)$
$\displaystyle\geq\frac{n-r}{k^{5}t_{a}(k^{g}+1)}-(k^{g}+1)$
$\displaystyle\geq\frac{n}{k^{5}t_{a}(k^{g}+2)}-(k^{g}+1).$

Case II. $k^{s}t_{\ell}\leq m<k^{s}t_{\ell}+k^{s-4}$. For any

$k^{s-1}(t_{\ell}-t_{\ell-1})+k^{s-5}+r<q\leq k^{s-1}(t_{\ell}-t_{\ell-2}),$

it can be seen that

$m-kq+r\in\left[k^{s}t_{\ell-2},k^{s}t_{\ell-1}\right)\subseteq A$

and

$k^{g-1}m+q\in\left[k^{s+g-1}t_{\ell},k^{s+g-1}t_{\ell+1}\right)\subseteq A,$

where the latter inclusion relation comes from the observation that
$\left[k^{s+g-1}t_{\ell},k^{s+g-1}t_{\ell+1}\right)$ contains in $A$ made
previously and the facts that

$k^{g-1}m+q<k^{g-1}\left(k^{s}t_{\ell}+k^{s-4}\right)+k^{s-1}(t_{\ell}-t_{\ell-2})\leq
k^{s+g-1}t_{\ell+1}$

since $k^{g}>T\geq 2(t_{\ell}-t_{\ell-2}).$ Note that

$n=(m-kq+r)+k\left(k^{g-1}m+q\right)$

for all these $q$, from which we conclude that

$\displaystyle r_{1,k}(A,n)$ $\displaystyle\geq
k^{s-1}(t_{\ell}-t_{\ell-2})-k^{s-1}(t_{\ell}-t_{\ell-1})-k^{s-5}-r$
$\displaystyle\geq\frac{1}{2}k^{s-1}-r$
$\displaystyle\geq\frac{n-r}{2kt_{a}(k^{g}+1)}-(k^{g}+1)$
$\displaystyle\geq\frac{n}{2kt_{a}(k^{g}+2)}-(k^{g}+1).$

Case III. $k^{s}t_{\ell+1}-k^{s-4}\leq m<k^{s}t_{\ell+1}$. It can be verified
directly that

$m+kq+r\in\left[k^{s}t_{\ell+2},k^{s}t_{\ell+3}\right)\subseteq A$

and

$k^{g-1}m-q\in\left[k^{s+g-1}t_{\ell},k^{s+g-1}t_{\ell+1}\right)\subseteq A,$

for any

$k^{s-1}(t_{\ell+2}-t_{\ell+1})+k^{s-5}\leq q\leq
k^{s-1}(t_{\ell+3}-t_{\ell+1})-r$

via similar arguments in Case II. In fact,

$k^{g-1}m-q\geq
k^{g-1}\left(k^{s}t_{\ell+1}-k^{s-4}\right)-k^{s-1}(t_{\ell+3}-t_{\ell+1})\geq
k^{s+g-1}t_{\ell}$

since $k^{g}>T\geq 2(t_{\ell+3}-t_{\ell+1})$. Note that

$n=(m+kq+r)+k\left(k^{g-1}m-q\right)$

for all these $q$, from which we conclude that

$\displaystyle r_{1,k}(A,n)$ $\displaystyle\geq
k^{s-1}(t_{\ell+3}-t_{\ell+1})-k^{s-1}(t_{\ell+2}-t_{\ell+1})-k^{s-5}-r$
$\displaystyle\geq\frac{1}{2}k^{s-1}-r$
$\displaystyle\geq\frac{n-r}{2kt_{a}(k^{g}+1)}-(k^{g}+1)$
$\displaystyle\geq\frac{n}{2kt_{a}(k^{g}+2)}-(k^{g}+1).$

This completes the proof of Theorem 1.1. ∎

######  Acknowledgments.

The authors would like to thank the anonymous referee for his helpful
comments.

The second author is supported by National Natural Science Foundation of China
under Grant No. 12201544, Natural Science Foundation of Jiangsu Province,
China, Grant No. BK20210784, China Postdoctoral Science Foundation, Grant No.
2022M710121, the foundations of the projects ”Jiangsu Provincial
Double–Innovation Doctor Program”, Grant No. JSSCBS20211023 and ”Golden
Phoenix of the Green City–Yang Zhou” to excellent PhD, Grant No.
YZLYJF2020PHD051.

The third author is supported by National Natural Science Foundation of China
under Grant No. 12101538.

## References

  * [1] S.–Q. Chen, The lower bound of weighted representation function, (preprint), http://arxiv.org//abs/2306.16025
  * [2] Y.–G. Chen, On the values of representation functions, Sci. China Math. 54 (2011), 1317–1331.
  * [3] Y.–G. Chen and M. Tang, Partitions of natural numbers with the same representation functions, J. Number Theory 129 (2009), 2689–2695.
  * [4] Y.–G. Chen and B. Wang, On additive properties of two special sequences, Acta Arith. 110 (2003), 299–303.
  * [5] G. Dombi, Additive properties of certain sets, Acta Arith. 103 (2002), 137–146.
  * [6] V.F. Lev, Reconstructing integer sets from their representation functions, Electron. J. Combin. 11 (2004), no. 1, R78.
  * [7] Y.–L. Li and W.–X. Ma, Partitions of natural numbers with the same weighted representation functions, Colloq. Math. 159 (2020), 1–5.
  * [8] Z. Qu, A note on representation functions with different weights, Colloq. Math. 143 (2016), 105–112.
  * [9] C. Sándor, Partitions of natural numbers and their representation functions, Integers 4 (2004), A18.
  * [10] A. Sárközy and V.T. Sós, On additive representation functions, in: The Mathematics of Paul Erdős I, R. Graham et al. (eds.), Springer, Berlin, 1997, 129–150.
  * [11] M. Tang, Partitions of the set of natural numbers and their representation functions, Discrete Math. 308 (2008), 2614–2616.
  * [12] Q.–H. Yang, Representation functions with different weights, Colloq. Math. 137 (2014), 1–6.
  * [13] Q.–H. Yang and Y.–G. Chen, Partitions of natural numbers with the same weighted representation functions, J. Number Theory 132 (2012), 3047–3055.


# BLOCK DISTRIBUTED 3MG ALGORITHM AND ITS APPLICATION TO 3D IMAGE RESTORATION

###### Abstract

Modern 3D image recovery problems require powerful optimization frameworks to
handle high dimensionality while providing reliable numerical solutions in a
reasonable time. In this perspective, asynchronous parallel optimization
algorithms have received an increasing attention by overcoming memory
limitation issues and communication bottlenecks. In this work, we propose a
block distributed Majorize-Minorize Memory Gradient (BD3MG) optimization
algorithm for solving large scale non-convex differentiable optimization
problems. Assuming a distributed memory environment, the algorithm casts the
efficient 3MG scheme [1] into smaller dimension subproblems where blocks of
variables are addressed in an asynchronous manner. Convergence of the sequence
built by the proposed BD3MG method is established under mild assumptions.
Application to the restoration of 3D images degraded by a depth-variant blur
shows that our method yields significant computational time reduction compared
to several synchronous and asynchronous competitors, while exhibiting great
scalability potential.

Index Terms—  Majorization-Minimization ; Block-alternating optimization ;
Distributed scheme ; Asynchronous communication ; Image deblurring ; Depth-
varying blur.

##  1 Introduction

Constantly improving image acquisition devices, from microscopes to medical
imaging machines, impose to work with increasingly large data. From the
mathematical perspective, many problems of image processing require to solve

$\operatorname{minimize}\limits_{x\in\mathbb{R}^{N}}f(x)$ (1)

with $f\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{N}\mapsto\mathbb{R}$ a
differentiable objective function. To limit the dependence of the optimization
process on the dimension $N$ of the problem, block alternating algorithms have
been developed. In these schemes, at each iteration only a subset of the
variables are updated, by minimizing $f$ with respect to only those variables,
the others being fixed. The blocks are selected iteratively following a cyclic
(or quasi-cyclic) order or a random rule. For a majority of problems
encountered in image restoration, the minimization of $f$ with respect to a
given block of variables is usually not possible in a closed form.
Furthermore, the application of such basic block coordinate descent update is
usually non desirable, as it may lead to convergence issues [2]. Better
performance and stability are obtained when combining the block alternating
approach with a so-called majoration-minimization (MM) scheme [3]. It consists
in building, at each iteration, a majorizing approximation for $f$ within the
active block of variables, whose minimizer has a more tractable form. Many
powerful algorithms fall within this framework, such as BSUM [4], PALM [5],
NMF [6], to name a few. By relying more closely on the structure of the
objective function, block alternating MM methods can reach fast convergence
rate [7, 8] while offering theoretical guarantees in non-convex cases [5, 9].

When the problem size becomes increasingly large, running such algorithm
becomes difficult, due to memory limitation issues. Parallel implementations
of MM schemes have been devised, where the block updates are performed
simultaneously, allowing to distribute computations on different nodes (or
machines) [10, 11, 12]. Implementation on parallel architecture requires to
pay attention to communication cost. The latter can be reduced by resorting to
an asynchronous parallel implementation. Each computation node has its own
iteration loop, so that it can keep updating its local variables without
having to wait for the update of the other, distant variables. This raises
challenging questions, in terms of convergence analysis, as the communication
delays may introduce instabilities. A plethora of recent works have focused on
proposing distributed optimization algorithms with assessed convergence, based
on stochastic proximal primal [13, 14] or primal-dual [15, 16, 17, 18, 19, 20]
techniques. However, as they rely on the formulation of dual instances of a
stochastic coordinate descent strategy, those algorithms are limited to convex
(sometimes even strongly convex) optimization and often require specific
probabilistic assumptions on the block update rule difficult to meet in
practice. In the context of MM algorithms, although the need for distributed
implementation strategies is crucial (see the discussion in [4]), few results
are available so far regarding theoretical convergence guarantees. For
example, we can mention the work of [21, 22], that proposes an asynchronous
version of PALM, with proven convergence in non-convex case, and good
practical behaviour [23].

In this paper, we focus on the MM Memory Gradient (3MG) method [1]. This
algorithm integrates a subspace acceleration strategy in the MM framework,
leading to one of the most efficient strategies for smooth optimization at
large scales [24]. Here, we propose to set forth an asynchronous distributed
version of this algorithm, called Block Distributed 3MG (BD3MG). The latter
allows to account efficiently for implementation on parallel architecture with
communication delays. By relying on the theoretical framework introduced in
[22], we derive convergence guarantees for BD3MG in the challenging non-convex
setting, under mild assumptions on the block update rule. We illustrate its
performance on the problem of 3D image restoration in the presence of depth-
variant blur. The rest of the paper is organized as follows. In Section 2, we
introduce our algorithm, analyse its principle features and state its
convergence properties. Section 3 focuses on its application to 3D image
deblurring. We conclude in Section 4.

##  2 Proposed method

###  2.1 Notations

We consider the Hilbert space $\mathbb{R}^{N}$ endowed with the usual scalar
product and norm $(\langle\cdot,\cdot\rangle,\|\cdot\|)$. Let a subset
$\mathcal{S}\subset[\\![1,N]\\!]$ of cardinality $|\mathcal{S}|$. We denote
$x_{(\mathcal{S})}=(x_{i})_{i\in\mathcal{S}}\in\mathbb{R}^{|\mathcal{S}|}$ the
restriction of $x$ to the coordinates in $\mathcal{S}$. In the same fashion,
$\nabla_{(\mathcal{S})}f(x)\in\mathbb{R}^{|\mathcal{S}|}$ denotes the gradient
of $f$, evaluated at $x$ along the components indexed in $\mathcal{S}$.
Finally, for a matrix $\smash{\mathcal{A}\in\mathbb{R}^{N\times N}}$, we
define the submatrix
$\mathcal{A}_{(\mathcal{S})}=\left([\mathcal{A}]_{p,p}\right)_{p\in\mathcal{S}}\in\mathbb{R}^{|\mathcal{S}|\times|\mathcal{S}|}$.

###  2.2 Block quadratic majorant function

Our approach relies on the use of block quadratic majorant functions which
constitute quadratic surrogates functions for the restriction of $f$ to any
set of coordinates $\mathcal{S}\subset[\\![1,N]\\!]$. Let
$\widetilde{x}\in\mathbb{R}^{N}$. Let us define :

$\forall\;v\in\mathbb{R}^{|\mathcal{S}|},\;\;Q_{(\mathcal{S})}(v,\widetilde{x})=f(\widetilde{x})+\langle\nabla_{(S)}f(\widetilde{x}),v-\widetilde{x}_{(S)}\rangle\\\
+\frac{1}{2}\langle
v-\widetilde{x}_{(S)},\mathcal{A}_{(S)}(\widetilde{x})(v-\widetilde{x}_{(S)})\rangle.$
(2)

Matrix
$\mathcal{A}_{(S)}(\widetilde{x})\in\mathbb{R}^{|\mathcal{S}|\times|\mathcal{S}|}$
hereabove is a symmetric definite positive matrix whose expression depends on
$\widetilde{x}$. Following the MM paradigm [3], it should be chosen so as to
fulfill:

$\forall v\in\mathbb{R}^{|\mathcal{S}|},\quad
Q_{(\mathcal{S})}(v,\widetilde{x})\geq f_{(\mathcal{S})}(v;\widetilde{x}),$
(3)

where $v\mapsto f_{(\mathcal{S})}(v;\widetilde{x})$ denotes the restriction of
function $f$ to coordinates in $\mathcal{S}$, the other coordinates being
fixed to those in vector $\widetilde{x}$. The existence and construction for
such majorant matrices, not discussed here due to the lack of space, is
addressed for instance in [11, 9, 24].

###  2.3 Block distributed MM memory gradient algorithm

We are now ready to present our BD3MG algorithm, assuming a memory-distributed
star cluster of C computing agents with a Master node connected to all other
agents. The Master loop starts at $x^{0}\in\mathbb{R}^{N}$ and generates the
sequence of iterates $(x^{k})_{k\in\mathbb{N}}$, that is incremented whenever
a worker $c\in[\\![1,C]\\!]$ updates a subset of coordinates of the global
variable. The main particularity of our algorithm is that the updates can
occur in an asynchronous fashion, thus reducing considerably idle time. We
denote $\mathcal{S}_{c}^{k}\subset[\\![1,N]\\!]$ the processing set of each
worker $c$ at times $k$, and $\mathbb{S}_{k}=\bigcup_{1\leq c\leq
C}(\mathcal{S}_{c}^{k})$ the total set of active blocks at that time. The
subset associated to any worker is allowed to change from one iteration to an
other. We only impose that there is no overlap in the coordinates updated by
the workers at a given time:

$(\forall
k\in\mathbb{N})\quad\bigcap\limits_{c\in\\{1,\ldots,C\\}}\mathcal{S}_{c}^{k}=\emptyset.$
(4)

At iteration $k\in\mathbb{N}$, the Master receives an updated increment
$d_{(\mathcal{S}_{c}^{k})}$ from a given worker $c\in[\\![1,C]\\!]$. The later
is used to increment the corresponding indexes $\mathcal{S}_{c}^{k}$ within
the global variable $x^{k-1}$, while the others remain untouched, which
defines $x^{k}$. The Master then decides for the new set of indexes
$\mathcal{S}_{c}^{k+1}$ to be treated by worker $c$. He informs worker $c$ of
his new task, and send him the triplet111Note that communication cost can be
reduced, when function $f$ reads as composition of terms with sparse
operators. See our discussion in [11] for an example.
$(x^{k},\mathcal{S}_{c}^{k+1},(x^{k}-x^{k-1}))_{(\mathcal{S}_{c}^{k+1})}$.

From the viewpoint of the workers, for each new triplet
$(\mathcal{S},x,d_{(\mathcal{S})})$ sent by the Master, a 3MG iterate is
performed. The later corresponds to the minimization of the block quadratic
majorant function (2), within the memory gradient subspace spanned by the two
columns of
$\mathcal{D}_{(\mathcal{S})}(x)=[-\nabla_{(\mathcal{S})}f(x)\;|\;d_{(\mathcal{S})}]\in\mathbb{R}^{|\mathcal{S}|\times
2}.$ This amounts to find $\widehat{u}$, a minimizer of $u\to
Q_{(\mathcal{S})}(x+\mathcal{D}_{(\mathcal{S})}(x)u,x)$, which can be obtained
through:

$\widehat{u}=(\mathcal{D}_{(\mathcal{S})}^{\top}(x)\mathcal{A}_{(\mathcal{S})}(x)\mathcal{D}_{(\mathcal{S})}(x))^{\dagger}\mathcal{D}_{(\mathcal{S})}(x)\nabla_{(\mathcal{S})}f(x),$
(5)

with $\cdot^{\dagger}$ the Moore-Penrose pseudo-inverse operator. The upcoming
tables summarize our BD3MG algorithm, that consists of two parts, one to be
executed by a ’Master’ computing node that receives updates and sends tasks, a
second one to be executed by all other computing nodes.

Block Distributed 3MG (Master)
$\begin{cases}\textbf{Initialization}\mathrel{\mathop{\mathchar 58\relax}}\\\
\text{Set }k=0,\;\;x^{0}\in\mathbb{R}^{N}.\;\;\\\ \text{For all
}c\in[\\![1,C]\\!],\text{ set }\mathcal{S}_{c}^{0}\subset[\\![1,N]\\!]\text{
s.t.}\bigcap\limits_{c\in[\\![1,C]\\!]}\mathcal{S}_{c}^{0}=\emptyset,\\\
\text{ and send }(x^{0},\mathcal{S}_{c}^{0},0_{|\mathcal{S}_{c}^{0}|})\text{
to worker }c.\\\ \text{Define
}\mathbb{S}_{0}=\bigcup_{c\in[\\![1,C]\\!]}\mathcal{S}_{c}^{0}.\\\
\textbf{While}\text{ a stopping criterion is not met:}\\\ \text{
}\small{(0)}\text{ Wait for any worker to send an update }\\\ \text{
}\small{(1)}\text{ Receive }(d_{(\mathcal{S}_{c}^{k})})\text{ from a worker
}c\\\ \\\ \text{ }\small{(2)}\text{ Update
}\begin{cases}x^{k+1}_{(\mathcal{S}_{c}^{k})}=x^{k}_{(\mathcal{S}_{c}^{k})}+d_{(\mathcal{S}_{c}^{k})}\\\
x^{k+1}_{(\overline{\mathcal{S}_{c}^{k}})}=x^{k}_{(\overline{\mathcal{S}_{c}^{k}})}\end{cases}\\\
\\\ \text{ }\small{(3)}\text{ Choose
}\mathcal{S}_{c}^{k+1}\subset[\\![1,N]\\!]/(\mathbb{S}_{k}/\mathcal{S}_{c}^{k})\\\
\qquad\text{For every }c^{\prime}\in[\\![1,C]\\!]/\\{c\\}\text{, set
}\mathcal{S}_{c^{\prime}}^{k+1}=\mathcal{S}_{c^{\prime}}^{k}.\\\
\qquad\text{Define
}\mathbb{S}_{k+1}=(\mathbb{S}_{k}/\mathcal{S}_{c}^{k})\cup\mathcal{S}_{c}^{k+1}\\\
\\\ \text{ }\small{(4)}\text{ Send
}(x^{k+1},\mathcal{S}_{c}^{k+1},(x^{k+1}-x^{k})_{(\mathcal{S}_{c}^{k+1})})\text{
to worker }c\\\ \text{ }\small{(5)}\text{ Increment }k=k+1\\\ \end{cases}$

Block Distributed 3MG (Worker) $\begin{cases}\textbf{While}\text{ the Master
stopping criterion is not met:}\\\ \text{ }\small{(1)}\text{ Receive
}(x,\mathcal{S},d_{(\mathcal{S})})\text{ from Master}\\\ \text{
}\small{(2)}\text{ Set
}\mathcal{D}_{(\mathcal{S})}(x)=[-\nabla_{(S)}f(x)\;|\;d_{(\mathcal{S})}]\\\
\text{ }\small{(3)}\text{ Compute }\mathcal{A}_{(\mathcal{S})}(x)\text{ and
}\nabla_{(S)}f(x)\\\ \text{ }\small{(4)}\text{
}\mathcal{B}_{(\mathcal{S})}(x)=(\mathcal{D}_{(\mathcal{S})}(x)\mathcal{A}_{(S)}(x)\mathcal{D}_{(\mathcal{S})}(x))^{\dagger}\\\
\text{ }\small{(6)}\text{
}d^{\prime}_{(\mathcal{S})}=-\mathcal{D}_{(\mathcal{S})}(x)\mathcal{B}_{(\mathcal{S})}(x)\mathcal{D}_{(\mathcal{S})}(x)\nabla_{(S)}f(x)\\\
\text{ }\small{(7)}\text{ Send }(d^{\prime}_{(\mathcal{S})})\text{ to the
Master}\end{cases}$

###  2.4 Analysis

In contrast with its parallel variant [11], BD3MG does not impose any locking
condition between workers to perform their computations. Therefore, latency
may appear in local variables. A local coordinate $x_{i}$ used by any worker
to perform its update at time $k$ belongs to a previous element
$x^{k^{\prime}_{i}}_{i}$ of the sequence $\\{x^{k}\\}_{k\in\mathbb{N}}$ with:
$k^{\prime}_{i}=\max\\{k^{\prime}\in[\\![0,k]\\!]\;|\;i\in\mathbb{S}_{k^{\prime}}\\}$.
We will model this latency at any iteration $k\in\mathbb{N}$ by introducing
$\delta_{k,n}=k-k^{\prime}_{n}\in[\\![0,k]\\!]$ the delay at a coordinate
$n\in[\\![1,N]\\!]$ and $\delta_{k}=(\delta_{k,n})_{n\in[\\![1,N]\\!]}$ the
complete vector of delays. For the sake of readability, we will denote
$x^{k,\delta_{k}}=\smash{(x_{n}^{k-\delta_{k,n}})_{n\in[\\![1,N]\\!]}}$.
Thanks to the majorizing condition (3), the sequence
$(x^{k})_{k\in\mathbb{N}}$ defined in BD3MG algorithm satisfies, for every
$k\in\mathbb{N}$, for every $c\in[\\![1,C]\\!]$,

$f(x^{k+1})\leq
Q_{(\mathcal{S}_{c}^{k})}(x^{k+1}_{(\mathcal{S}_{c}^{k})},x^{k,\delta_{k}})\leq
Q_{(\mathcal{S}_{c}^{k})}(x^{k,\delta_{k}}_{(\mathcal{S}_{c}^{k})},x^{k,\delta_{k}})=f(x^{k,\delta_{k}}).$

When there is no delay , i.e. $\delta_{k}\equiv 0$, we go back to a
synchronous algorithm, benefiting from the classical block descent result
inherent to MM schemes:

$f(x^{k+1})\leq
Q_{(\mathcal{S}_{c}^{k})}(x^{k+1}_{(\mathcal{S}_{c}^{k})},x^{k})\leq
Q_{(\mathcal{S}_{c}^{k})}(x^{k}_{(\mathcal{S}_{c}^{k})},x^{k})=f(x^{k})$.

###  2.5 Hypothesis and convergence

We now state our convergence result for the sequence
$(x^{k})_{k\in\mathbb{N}}$ resulting from the BD3MG algorithm presented in
Section 2.3. We first introduce the following assumptions.

###### Assumption 1

Function $f$ is differentiable, bounded from below and semi-algebraic222Real
semi-algebraic functions represent a wide class of functions that satisfy the
Kurdyka-Lojasiewicz inequality, which is at the core of our convergence study.
See [25, 5] for further details.. Moreover, $f$ has an
$\mathcal{L}$-Lipschitzian gradient on $\mathbb{R}^{N}$ with $\mathcal{L}>0$,
i.e.

$\forall(x,y)\in(\mathbb{R}^{N})^{2},\;\|\nabla f(x)-\nabla
f(y)\|\leq\mathcal{L}\|x-y\|.$ (6)

######  Assumption 2 (Boundedness of delay)

Under the convention that
$\forall\;l\in\mathbb{N}^{*},\;\mathbb{S}_{-l}=\emptyset$, there exists
$\tau\in\mathbb{N}$ such that

$\forall
k\in\mathbb{N},\quad[\\![1,N]\\!]\subset\bigcup\limits_{i=k-\tau}^{k}\mathbb{S}_{i}.$
(7)

######  Assumption 3 (Curvature of quadratic majorant)

Let us denote $(\Gamma^{k}_{c})_{k\in\mathbb{N},c\in[\\![1,C]\\!]}$ the
sequence of matrices defined as:

$\forall k\in\mathbb{N},\;\forall
c\in[\\![1,C]\\!],\quad\Gamma^{k}_{c}=\mathcal{A}_{(\mathcal{S}^{k}_{c})}(x^{k,\delta_{k}})-\frac{1}{2}\mathcal{A}_{(\mathcal{S}^{k}_{c})}(x^{k}).$

There exists $(\underline{\nu},\overline{\nu})>0$ such that, for every
$k\in\mathbb{N}$, for every $c\in[\\![1,C]\\!]$,

$(\mathcal{L}\sqrt{\tau}+\underline{\nu})\mathrm{Id}_{|\mathcal{S}^{k}_{c}|}\preceq\Gamma^{k}_{c}\preceq\overline{\nu}\mathrm{Id}_{|\mathcal{S}^{k}_{c}|}.$
(8)

Under those assumptions, we state our convergence result for BD3MG, whose
proof, relying on the analysis from [22, 1], is skipped by lack of space.

######  Theorem 1 (Convergence of BD3MG)

Let $f\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{N}\mapsto\mathbb{R}$,
and $x^{0}\in\mathbb{R}^{N}$. If Assumptions 1-2-3 are verified, then the
sequence $(x^{k})_{k\in\mathbb{N}}$ built by BD3MG algorithm converges
globally to a stationary point $x^{*}$ of $f$.

It is worthy to point out that, in contrast with most existing works in the
literature of distributed optimization, no convexity assumption is required in
our analysis. The MM framework allows to make use of the recent theory of non-
smooth analysis [25], as we also have shown in our previous works [9, 1].

##  3 Application to 3D image deblurring

###  3.1 Model and objective function

We consider the restoration of a 3D microscopic volume $\overline{x}$ of size
$N=N_{\mathsf{X}}\times N_{\mathsf{Y}}\times N_{\mathsf{Z}}$, given a degraded
observation $y$, altered by a depth-variant blur operator $H$ and an additive
white noise $b$:

$y=H\overline{x}+b.$ (9)

The associated inverse problem can be solved efficiently by minimizing a least
squares penalization penalized by a smooth 3D regularization function.
Specifically, $f$ takes the form:

$\forall x\in\mathbb{R}^{N},\quad f(x)=\sum_{s=1}^{4}f_{s}(L_{s}x)$ (10)

with $f_{1}\circ L_{1}=\frac{1}{2}\|H\cdot-y\|^{2}$, $f_{2}\circ
L_{2}=\eta\,d^{2}_{[x_{\min},x_{\max}]^{N}}$, $f_{3}\circ
L_{3}=\lambda\sum_{n=1}^{N}\sqrt{\left([V^{\mathsf{X}}\cdot]_{n}\right)^{2}+\left([V^{\mathsf{Y}}\cdot]_{n}\right)^{2}+\delta^{2}}$,
and $f_{4}\circ L_{4}=\kappa\|V^{\mathsf{Z}}\cdot\|^{2}$. Hereabove,
$V^{\mathsf{X}}\in\mathbb{R}^{N\times N}$,
$V^{\mathsf{Y}}\in\mathbb{R}^{N\times N}$,
$V^{\mathsf{Z}}\in\mathbb{R}^{N\times N}$ state for discrete gradient
operators along the three directions of the volume, $x_{\min}$ (resp
$x_{\max}$) $\in\mathbb{R}$ are minimal (resp. maximal) bounds on the sought
intensity values and $d_{E}$ is the Euclidian distance to set $E$. Function
$f$ satisfies Ass. 1 as a combination of squared distances and discretized
total-variation distances. If not specified otherwise, the majorizing matrices
for $f$ will be constructed following the strategy in [1], ensuring the
fulfillment of Ass. 3. The convolution operator $H$ simulates a depth-varying
3D Gaussian blur. For each depth $z\in\\{1,\ldots,N_{\mathsf{Z}}\\}$, the blur
kernel is characterized by different variance and rotation parameters whose
values are taken as random realizations of uniform distributions. The depth-
variant structure of the blur model motivates us to split the vector $x$ along
the dimension $\mathsf{Z}$, assigning computing agents (i.e. workers) to one
(or several) selected slice(s) of the 3D volume. For the sake of simplicity, a
sequential ordering is used, for the selection rule of the processed blocks,
hence Ass. 2 is fulfilled. Furthermore, the same strategy as in [11] is
employed to control the memory usage during communications master/worker. All
algorithms are initialized with $x^{0}\in\mathbb{R}^{N}$ whose entries are
uniformly sampled in $[0,\max(y)]$. Parameters $\lambda,\eta,\delta,\kappa>0$
will be tuned manually, for both presented examples, so as to maximize the
Signal-to-Noise Ratio (SNR) of the restored volume. Two microscopic images,
namely FlyBrain and Aneurysm will be considered. Example of restoration
results are presented in Fig. 2.

###  3.2 Ablation study and comparative analysis

In order to measure the performance improvement allowed by BD3MG, we conducted
an “ablation study”, that consists in removing some/all acceleration features
of our algorithm, namely asynchrony, subspace line-search and MM scaling. We
then obtain, in addition to BD3MG, four parallel/distributed optimization
algorithms that we detail hereafter. The _asynchronous gradient descent_
(Async-GD) corresponds to the algorithm from [13], obtained by limiting the
subspace to the sole gradient descent direction
$\mathcal{D}_{(\mathcal{S})}(x)=-\nabla_{(\mathcal{S})}f(x)$ and by setting
$\mathcal{A}_{(\mathcal{S})}=\mathcal{L}\,\mathrm{Id}_{|\mathcal{S}|}$, with
$\mathcal{L}$ the Lipschitz constant of $\nabla f$. The _asynchronous
conjugate gradient_ algorithm (Async-CG) identifies with our BD3MG method when
using the basic MM metric
$\mathcal{A}_{(\mathcal{S})}=\mathcal{L}\,\mathrm{Id}_{|\mathcal{S}|}$. The
resulting algorithm can be viewed as a distributed version of the nonlinear
conjugate gradient method, with closed form stepsize from [26]. The
_asynchronous MM algorithm_ (Async-MM) is obtained by removing the subspace
acceleration strategy in BD3MG, so that
$d^{\prime}_{S}=-\mathcal{A}_{(\mathcal{S})}(x)^{-1}\nabla_{(\mathcal{S})}f(x)$
in the Worker loop. The latter inversion is performed with linear biconjugate
gradient solver. Async-MM can be interpreted as a distributed implementation
of a half-quadratic algorithm [27]. Finally, _block parallel 3MG algorithm_
(BP3MG) is the method we originally proposed in [11], where blocks updates are
assumed to be computed at the same time without any communication delay. The
names and characteristics of the resulting five tested methods are summarized
in Tab. 1.

Name | Asynchrony | Memory | MM scaling
---|---|---|---
Async-GD | ✓ | ✗ | ✗
Async-CG | ✓ | ✓ | ✗
Async-MM | ✓ | ✗ | ✓
BP3MG | ✗ | ✓ | ✓
BD3MG | ✓ | ✓ | ✓
Table 1: Algorithmic features of the compared approaches. Fig. 1: Evolution of
SNR in dB (left) and relative distance to solution $\|x^{k}-x^{*}\|/\|x^{*}\|$
(right) along time (in seconds) for FlyBrain restoration. Results are averaged
over ten noise realizations.

We perform our comparative analysis, using the microscopic image FlyBrain with
size $N=256\times 256\times 24$. It is degraded by depth-variant blur kernels
of size $11\times 11\times 21$, and a zero-mean white Gaussian noise with
standard deviation $0.04$. The results are averaged over ten noise
realizations, and the initial signal to noise ratio is around $11.6$ dB. Fig.
1 illustrates the evolution of the SNR of the restored image along time in
seconds, for each method, for experiments performed using Python 3, running an
Intel® Xeon(R) W-2135 CPU with 12 cores clocked at 3.70GHz. We also display
the relative distance to the solution $x^{*}$, computed after a large number
of iterations (typically, $10^{4}$), characterized by an average SNR of
$14.33$ dB. One can see that the proposed BD3MG method clearly outperforms the
others in terms of time to reach close-to-optimal solution. Async-MM led to
the slowest convergence, as it requires linear system inversion at each
iteration. The superiority of Async-CG over Async-GD illustrates the benefits
for including the memory term within the subspace. Finally, BD3MG reaches
convergence faster than its synchronous counterpart, BP3MG, probably thanks to
the removal of the locking conditions in the Master loop.

---
Fig. 2: Comparison between original (left), degraded (middle) and restored
(right) slices ($z=10$) of FlyBrain and Aneurysm.

###  3.3 Linear Speedup

In order to assess the scalability of the BD3MG algorithm, we further analysed
the speed-up of the optimization process when a High Parallel Computing
computer is being used. Namely, we ran BD3MG and BP3MG on an Intel Xeon CPU
6148 with up to 80 physical cores at 2.4 GHz (Skylake) and 1.5 Tio of RAM.
Image Aneurysm with size $N=155\times 154\times 79$ is degraded by blur
kernels of size to $5\times 5\times 11$, and noise standard deviation of
$0.04$, so that the initial SNR is $6.44$ dB, while the restored SNR is
$11.92$ dB. Fig. 3 presents the acceleration ratio between the required
computation time for two cores (i.e. one Master and one Worker) versus the
computation time when activating from 10 to 80 cores, for reaching the
stopping criterion $\|x^{k+1}-x^{k}\|\leq 10^{-6}\|x^{k}\|$. This illustrates
the great potential of scalability of the proposed algorithm. Asynchronicity
of BD3MG allows to improve the speed-up, in comparison to the one exhibited by
BP3MG [11]. Finally, as the number of core increases, a mild saturation effect
is observed (in agreement with Amdahl’s law [28]).

Fig. 3: Speed-up ratio for BD3MG (blue) and BP3MG (red), with respect to the
number of active cores for the restoration of Aneurysm.

##  4 Conclusion

In this paper, we have presented a new block distributed Majorize-Minimize
Memory Gradient algorithm to tackle a wide class of large scale optimization
problems. The main feature of our method lies in its distributed asynchronous
formulation that allows for delays among workers, while theoretically
maintaining the convergence guarantees and the practical performance of the
powerful 3MG scheme. The new algorithm has been tested in the context of 3D
image restoration under depth-variant blur. Experimental results underlined
its scalability and efficiency333The code of BD3MG for depth-variant image
deblurring has been made available at https://github.com/mathieuchal/BD3MG.
Future works will be focused on the extension to more general graph
topologies.

## References

  * [1] E. Chouzenoux, A. Jezierska, J.-C. Pesquet, and H. Talbot,  “A majorize-minimize subspace approach for $\ell_{2}-\ell_{0}$ image regularization,”  SIAM Journal on Imaging Sciences (SIIMS), vol. 6, no. 1, pp. 563–591, January 2013.
  * [2] P. Tseng,  “Convergence of a block coordinate descent method for nondifferentiable minimization,”  Journal of Optimization, Theory and Applications, vol. 109, no. 3, pp. 475–494, 2001.
  * [3] M. W. Jacobson and J. A. Fessler,  “An expanded theoretical treatment of iteration-dependent majorize-minimize algorithms,”  IEEE Transactions on Image Processing, vol. 16, no. 10, pp. 2411–2422, Oct 2007.
  * [4] M. Hong, M. Razaviyayn, Z. Luo, and J. Pang,  “A unified algorithmic framework for block-structured optimization involving big data: With applications in machine learning and signal processing,”  IEEE Signal Processing Magazine, vol. 33, no. 1, pp. 57–77, Jan 2016.
  * [5] J. Bolte, S. Sabach, and M. Teboulle,  “Proximal alternating linearized minimization for nonconvex and nonsmooth problems,”  Mathematical Programming, vol. 146, no. 1, pp. 459–494, Aug 2014\.
  * [6] D. D. Lee and H. S. Seung,  “Algorithms for non-negative matrix factorization,”  in Proceedings of the 13th International Conference on Neural Information Processing Systems (NIPS 2000), Denver, Colorado, 2000, p. 535–541.
  * [7] J. A. Fessler,  “Grouped coordinate descent algorithms for robust edge-preserving image restoration,”  in Image Reconstruction and Restoration II, Timothy J. Schulz, Ed. International Society for Optics and Photonics, 1997, vol. 3170, pp. 184 – 194, SPIE.
  * [8] L. Duval E. Chouzenoux A. Repetti, M. Q. Pham and J.-C. Pesquet,  “Euclid in a taxicab: Sparse blind deconvolution with smoothed $\ell_{1}$/$\ell_{2}$ regularization,”  IEEE Signal Processing Letters, vol. 22, no. 5, pp. 539–543, May 2015.
  * [9] E. Chouzenoux, J.-C. Pesquet, and A. Repetti,  “A block coordinate variable metric forward-backward algorithm,”  Journal of Global Optimization, vol. 66, no. 3, pp. 457–485, 2015\.
  * [10] S. Sotthivirat and J. A. Fessler,  “Image recovery using partitioned-separable paraboloidal surrogate coordinate ascent algorithms,”  IEEE Transactions on Signal Processing, vol. 11, no. 3, pp. 306–317, 2002.
  * [11] S. Cadoni, E. Chouzenoux, J. Pesquet, and C. Chaux,  “A block parallel majorize-minimize memory gradient algorithm,”  in Proceedings of the IEEE International Conference on Image Processing (ICIP 2016), Phoenix, Arizona, 25-28 Sep. 2016, pp. 3194–3198.
  * [12] G. E. Moon, A. Sukumaran-Rajam, S. Parthasarathy, and P. Sadayappan,  “PL-NMF: parallel locality-optimized non-negative matrix factorization,” 2019,  http://arxiv.org/abs/1904.07935.
  * [13] F. Niu, B. Recht, C. Re, and S. J. Wright,  “Hogwild: A lock-free approach to parallelizing stochastic gradient descent,”  in Proceedings of the 25th Conference on Advances in Neural Information Processing Systems (NIPS 2011), pp. 693–701. Granada, Spain, 12-17 Dec. 2011.
  * [14] D. Grishchenko, F. Iutzeler, J. Malick, and M.R. Amini,  “Asynchronous distributed learning with sparse communications and identification,” 2018,  https://arxiv.org/abs/1812.03871.
  * [15] R. Zhang and J. T. Kwok,  “Asynchronous distributed ADMM for consensus optimization,”  in Proceedings of the 31st International Conference on Machine Learning (ICML 2014), 21-26 June 2014, pp. 1701–1709.
  * [16] J.-C. Pesquet and A. Repetti,  “A class of randomized primal-dual algorithms for distributed optimization,”  Journal on Nonlinear Convex Analysis, vol. 16, no. 12, pp. 2353–2490, Dec. 2015.
  * [17] R. Hannah and W. Yin,  “On unbounded delays in asynchronous parallel fixed-point algorithms,”  Journal of Scientific Computing, vol. 76, no. 1, pp. 299–326, Dec 2017.
  * [18] F. Abboud, E. Chouzenoux, J.-C. Pesquet, and H. Talbot,  “Distributed algorithms for proximity operator computation with applications to video processing,” 2019,  https://hal.archives-ouvertes.fr/hal-01942710.
  * [19] A. Onose, R. E. Carrillo, A. Repetti, J. D. McEwen, J.-T. Thiran, J.-C. Pesquet, and Y. Wiaux,  “Scalable splitting algorithms for big-data interferometric imaging in the SKA era,”  Monthly Notices of the Royal Astronomical Society, vol. 462, no. 4, pp. 4314–4335, 2016.
  * [20] J. Xu, Y. Sun, Y. Tian, and G. Scutari,  “A unified contraction analysis of a class of distributed algorithms for composite optimization,” 2019,  https://arxiv.org/abs/1910.09817.
  * [21] D. Davis, M. Udell, and B. Edmunds,  “The sound of APALM clapping: Faster nonsmooth nonconvex optimization with stochastic asynchronous palm,”  in Proceedings of the 30th International Conference on Neural Information Processing Systems (NIPS 2016), 5-10 Dec. 2016, p. 226–234.
  * [22] D. Davis,  “The asynchronous palm algorithm for nonsmooth nonconvex problems,” 2016,  https://arxiv.org/abs/1604.00526.
  * [23] P. Thouvenin, N. Dobigeon, and J. Tourneret,  “Partially asynchronous distributed unmixing of hyperspectral images,”  IEEE Transactions on Geoscience and Remote Sensing, vol. 57, no. 4, Apr. 2019.
  * [24] Y. Sun, P. Babu, and D. P. Palomar,  “Majorization-minimization algorithms in signal processing, communications, and machine learning,”  IEEE Transactions on Signal Processing, vol. 65, no. 3, pp. 794–816, 2016.
  * [25] H. Attouch, J. Bolte, and B.F. Svaiter,  “Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward–backward splitting, and regularized gauss–seidel methods,”  Mathematical Programming, vol. 137, no. 1, pp. 91–129, Feb 2013\.
  * [26] C. Labat and J. Idier,  “Convergence of conjugate gradient methods with a closed-form stepsize formula,”  Journal of Optimization Theory and Applications, vol. 136, no. 1, pp. 43–60, Jan 2008.
  * [27] M. Allain, J. Idier, and Y. Goussard,  “On global and local convergence of half-quadratic algorithms,”  IEEE Transactions on Image Processing, vol. 15, no. 5, pp. 1130–1142, May 2006.
  * [28] G. M. Amdahl,  “Validity of the single processor approach to achieving large scale computing capabilities,”  in Proceedings of the American Federation of Information processing Societies conference (AFIPS 1967), Atlantic City, 18-20 Apr. 1967, pp. 483–485.


# A small note on variation in segmentation annotations

Silas Nyboe Ørting
Center for Quantification of Imaging data from MAX IV (QIM)
Department of Computer Science
University of Copenhagen
<EMAIL_ADDRESS>

###### Abstract

We report on the results of a small crowdsourcing experiment conducted at a
workshop on machine learning for segmentation held at the Danish Bio Imaging
network meeting 2020. During the workshop we asked participants to manually
segment mitochondria in three 2D patches. The aim of the experiment was to
illustrate that manual annotations should not be seen as the ground truth, but
as a reference standard that is subject to substantial variation. In this note
we show how the large variation we observed in the segmentations can be
reduced by removing the annotators with worst pair-wise agreement. Having
removed the annotators with worst performance, we illustrate that the
remaining variance is semantically meaningful and can be exploited to obtain
segmentations of cell boundary and cell interior.

_K_ eywords Annotation $\cdot$ Crowdsourcing $\cdot$ Mitochondria $\cdot$ EM

##  1 Introduction

Manual annotation is the gold standard for many image analysis tasks and a
crucial part of developing machine learning methods. Despite well defined
annotation protocols and training most tasks will be subject to substantial
inter- and intra-annotator variation[3, 1]. Although domain knowledge and
experience are important for understanding and interpreting the images, it is
not a magic tool that eliminates variation. On the contrary, [3] did not find
that more experienced radiologists had less variation in annotations. A
compelling illustration that lack of familiarity with a task can be helpful is
provided by [2] where simply flipping images improved annotation performance
of crowd workers.

We can split variation into different parts: variation due to errors, for
example, unintended clicks or annotating the wrong structure; variation due to
interpretation of the task, for example, including or excluding edema when
asked to segment a brain tumor; and variation due to the inherent ambiguity of
the task, for example, partial volume effects and poor signal. The first kind
of variation is noise we would like to eliminate completely; the second kind
is semantically meaningful information we would like to exploit; and the third
kind is a measure of uncertainty regarding the actual true segmentation which
we would like to explore.

In this report we illustrate simple approaches for investigating variation,
removing noisy annotations and exploiting variation in task interpretation. We
hope this note can serve as a reminder that variation in annotations should
not be ignored.

##  2 Materials & Methods

The crowdsourcing experiment was conducted during an online workshop on
machine learning for segmentation held at the Danish Bio Imaging network
meeting 2020111https://www.conferencemanager.dk/danish-bioimaging-
meeting-2020/. Approximately 45 persons participated in the workshop, which
was held online using Zoom222https://zoom.us/. We did not record information
regarding the participants scientific background, but assume all have some
familiarity with bioimaging. In the experiment we asked participants to
segment mitochondria in EM images using an online annotation tool.

###  2.1 Data

We used three 2D patches extracted from the Electron Microscopy Dataset
released by the CVLAB at EPFL333https://www.epfl.ch/labs/cvlab/data/data-em/.
The three patches, extracted from the slice 1, 34, and 67 of the training sub-
volume, are shown in Figure 1 alongside reference segmentations from the
published dataset.

Figure 1: Patches used for crowdsourcing experiment. Top: EM images. Bottom:
Reference segmentations

###  2.2 Annotation tools

We used two online drawing tools for the annotation task,
Pixilart444https://www.pixilart.com/ and Pixlr555https://pixlr.com/, both
illustrated in Figure 2

Figure 2: The two annotation tools used for the experiments. Top: Pixilart.
Bottom: Pixlr

###  2.3 Post-processing

All annotations where binarized by setting a pixel to 1 if the sum of channels
where larger than 0. Six annotations included the background image which was
removed prior to binarization. Segmentations from one annotator had speckle
noise, which was removed by finding connected components and discarding all
components with size of one pixel.

##  3 Experiments

The annotation experiment lasted for about 30 minutes and was initiated after
a talk on training a U-Net to segment mitochondria. Participants where given a
quick overview of the task and a link to written instructions and data. After
the overview, participants where randomly assigned to one of two groups, one
using Pixilart and one using Pixlr, using the "breakout rooms" feature of
Zoom. Each group was then given a live demonstration of how to perform the
annotation task, before working on their own.

##  4 Results & Discussion

Of the approximately 45 participants, 21 uploaded segmentations of one image,
16 of two images and 13 of all three images. The mean segmentations for each
image is shown alongside the reference segmentations in Figure 3. Although the
mean segmentations are very close to the references, it is clear that there is
substantial variation at the boundary. The extent of this variation is
illustrated in Figure 4 where the union and (union - intersection) of all
segmentations are shown. In the top and middle rows we can see that some
annotators missed mitochondria and some added mitochondria. From Figure 3 and
Figure 4 it is clear that at least some variation is due to low quality
annotations.

Figure 3: Top: Mean segmentations. Bottom: Reference segmentations. From left
to right: Slice 1, 34, 67

Figure 4: Top: Union of all segmentations. Bottom: Difference between union
and intersection of all segmentations. From left to right: Slice 1, 34, 67

One approach to removing poor annotations is to assume that good annotations
will agree with other good annotations, whereas poor annotations will disagree
with both good and poor annotations. Since we are dealing with binary
segmentations we use the Dice score as measure of agreement and calculate
pairwise agreement for all annotators. Figure 5 shows the distribution of
pairwise dice scores for each annotator on each of the three images. The left
column shows scores when all annotators are included, the right column shows
scores when annotators with median pairwise dice less than 0.9 are excluded.
Five annotations where excluded from slice 1, six from slice 34, and four from
slice 67.

Figure 5: Pairwise dice scores. Dashed red line indicates 0.9 Dice. Top:
Before removing annotators. Bottom: After removing annotators. From left to
right: Slice 1, 34, 67.

The result of excluding annotations with poor agreement is illustrated in
Figure 6 where the union and intersection are recomputed for the reduced set
of annotations. Results are clearly better with most variation at a comparably
small part of the cell boundaries. Figure 7 illustrates that the variation is
most likely due to different interpretations of the task, as the part with
variation corresponds quite well to the boundary of the cell. Although in
general it would be nice to have annotators perform the exact same task, we
should be careful treating variation as inherently bad once we have ensured
sufficient quality. In this case we get segmentations of the boundary and
interior for "free" and can easily get a full cell segmentation as the union.

Figure 6: Top: Union of all segmentations. Bottom: Difference between union
and intersection of segmentations after removing annotators. From left to
right: Slice 1, 34, 67

Figure 7: Annotations for slice 67 showing variation can be informative. In
this case the area of disagreement corresponds quite well to the boundary of
the cell. Left: Intersections. Right: (Union - intersection)

## References

  * [1] Anton S Becker et al.  “Variability of manual segmentation of the prostate in axial T2-weighted MRI: A multi-reader study”  In _European journal of radiology_ 121 Elsevier, 2019, pp. 108716
  * [2] Danna Gurari, Mehrnoosh Sameki and Margrit Betke  “Investigating the influence of data familiarity to improve the design of a crowdsourcing image annotation system”  In _Human Computation (HCOMP)_ , 2016
  * [3] Leo Joskowicz, D Cohen, N Caplan and J Sosna  “Inter-observer variability of manual contour delineation of structures in CT”  In _European radiology_ 29.3 Springer, 2019, pp. 1391–1399