# Limits and Signatures of Relativistic Spaceflight

Ulvi Yurtsever<EMAIL_ADDRESS>Steven Wilkinson
<EMAIL_ADDRESS>Raytheon Company, El Segundo, CA 92195

###### Abstract

While special relativity imposes an absolute speed limit at the speed of
light, our Universe is not empty Minkowski spacetime. The constituents that
fill the interstellar/intergalactic vacuum, including the cosmic microwave
background photons, impose a lower speed limit on any object travelling at
relativistic velocities. Scattering of cosmic microwave phtotons from an
ultra-relativistic object may create radiation with a characteristic signature
allowing the detection of such objects at large distances.

Introduction

At a fixed speed $v$, the distance travelled by a spacecraft within a given
(proper) travel time $\tau$ scales as $\gamma v\tau$, which for relativistic
speeds approximately equals $\gamma c\tau$, where
$\gamma=(1-v^{2}/c^{2})^{-\nicefrac{{1}}{{2}}}$. Hence range is directly
proportional to $\gamma$, and to cover intergalactic distances within a
limited “lifetime” $\tau$ requires $\gamma\gg 1$. The primary obstacles to
relativistic ($\gamma\gg 1$) space travel would be collisions with
interstellar dust particles (cosmic dust) and larger space objects, which will
impact with kinetic energies of $(\gamma-1)mc^{2}$ for rest mass $m$, similar
to collisions in particle accelerators but potentially at much higher
energies. For large enough $\gamma$, even molecular collisions could be a
significant source of drag and possibly damaging. For example, at $\gamma=2$ a
baseball size object of mass 150g has an impact energy equivalent to 36
Megatons of TNT; a single cosmic dust grain of mass $10^{-14}\rm g$ at
$\gamma=10^{8}$ has an impact energy of close to 24 kgs of TNT. tiplerbook ;
weinbergbook ; griffithsbook ; perkinsbook

The intergalactic medium has much less debris per unit volume as compared to
the interstellar medium so the chance of a catastrophic collision occurring is
significantly reduced beyond the confines of any galaxy. We define extreme
velocities to occur at $\gamma>9.1$ which is the threshold velocity for
proton-proton collisions to produce antiprotons. Now to travel at a safe speed
in the interstellar medium one would need a $\gamma<1.3$, the threshold for
the pion production in proton-proton collisions. The problem of colliding with
larger objects as described above remains, but we assume that sufficiently
advanced technology likely to accompany relativistic space travel
capabilitites will be able to circumvent these interactions. One can imagine,
though, intergalactic travel through mostly matter-free space, and under those
conditions the cosmic microwave background (CMB) photons are primarily the
particles a spacecraft will encounter in collisions.

CMB Interactions

A fast-moving spacecraft traveling in intergalactic space still has to contend
with collisions with cosmic microwave photons, which, at relativistic speeds,
will appear in the spacecraft frame as highly energetic gamma rays.
Interactions of CMB photons with the material of the spacecraft hull will have
effects ranging from ionization and Compton scattering to pair production with
increasing $\gamma$. We will assume that advanced technology can mitigate the
harmful effets of ionization and Compton scattering, and therefore concentrate
on pair production as the obstacle that is most likely to resist removal via
technology. tiplerbook

The first speed level where pair production will pose a challenge for
relativistic spaceflight engineers is photon-nucleon interactions in the hull
of the spacecraft. Let us consider a hull made of ordinary baryonic matter. In
this case, a single CMB photon will create an electron-positron pair via its
collision with a nucleus in the hull if its energy exceeds the rest mass of an
electron and positron. That threshold is at blueshifted photon energies of
about 1 MeV. In terms of the frequency $\omega$ in the rest frame of the
spacecraft, the condition for nucleus-mediated pair production is

$\omega\mathchar 13374\relax\frac{2m_{e}c^{2}}{\mathchar 1406\relax}\approx
2.47\times 10^{20}\mbox{Hz}$ (1)

The central frequency of the (Planck) distribution of cosmic microwave
background photons is 160GHz. The blueshift factor to transform a typical
counter-directed CMB photon to the rest frame of the spacecraft is then
$\omega=\gamma(1+v/c)\times 10^{12}\,\mbox{Hz}$, where $v$ is the velocity of
the spacecraft and $\gamma=(1-v^{2}/c^{2})^{-\nicefrac{{1}}{{2}}}$. So the
condition for pair creation in terms of spacecraft speed is:

$\gamma\left(1+\frac{v}{c}\right)\mathchar 13374\relax 2.47\times 10^{8}$ (2)

corresponding to a $\gamma$ of about $\gamma\approx 1.24\times 10^{8}$, or
$v/c\approx 1-3.3\times 10^{-17}$.

How much energy will be dissipated at this speed or above due to pair
production? There are 412 cosmic microwave photons in the Universe per cubic
centimeter (see Eq. (13) below). Assuming the spacecraft has a crosssectional
area of $A$ (measured in cm2), it will be colliding with approximately
$1.2\times 10^{13}A$ microwave background photons per second. Each created
electron-positron pair dissipates $1.6\times 10^{-13}$ Joules of energy, thus
the viscous-dissipative effect of the cosmic background on the spacecraft
would be energy dissipation of approximately $2A$ Joules per second. Assuming
an effective cross-sectional area of say 100 square meters, $A\approx
10^{6}\mbox{cm}^{2}$, the dissipative effect is about 2 million Joules per
second, but with the time rate of change measured in the rest frame of the
cosmic microwave background. In the spacecraft’s rest frame, 1 “cosmic
background second” corresponds to $1/\gamma$ proper time seconds; so energy
dissipation is a factor of $\gamma$ higher, or about $2.5\times 10^{14}$
Joules per proper time second, which is a large drag to be overcome by the
spacecraft engines just to keep it moving at constant velocity. Therefore
shutting down the propulsion system when the desired speed is reached is not
an option at large $\gamma$ values. In general, the drag on the spacecraft
scales as $\gamma^{2}$, so, for example, at $\gamma\mathchar 13366\relax 10$
the drag will be less than a few Joules per second and thus negligible.

At the next level, we can imagine technology so advanced that the spacecraft
“hull” does not contain any ordinary baryonic matter but some other kind of
physical matter field. However, regardless of how advanced the spacecraft hull
material can be, at sufficiently relativistic velocities pair production of
CMB photons via the Schwinger breakdown of vacuum will become significant. The
Schwinger limit corresponds to the electric field strength at which
spontaneous electron-positron pair production from the electromagnetic field
(photons) becomes energetically favorable:

$eE_{S}\frac{h}{m_{e}c}=2m_{e}c^{2}$ (3)

or

$E_{S}=\frac{{m_{e}}^{2}c^{3}}{\pi e\mathchar 1406\relax}\approx 1.4\times
10^{14}\;(\mbox{erg/cm${}^{3}$})^{\nicefrac{{1}}{{2}}}$ (4)

Here $h/m_{e}c$ is the electron Compton wavelength, and $e$ is the electronic
charge. Let $\delta$ be the effective “skin depth” for the
interaction/penetration of high energy photons into the “material” making up
the spacecraft (or its forward shield). Then the colliding gamma rays will
pair-create provided the electric field $E$ associated with their localization
at the length scale $\delta$ and with their frequency $\omega$ in the rest
frame of the spacecraft exceeds the Schwinger limit:

$E\approx\sqrt{\frac{8\pi\mathchar 1406\relax\omega}{\delta^{3}}}\mathchar
13374\relax E_{S}=\frac{{m_{e}}^{2}c^{3}}{\pi e\mathchar 1406\relax}$ (5)

In terms of the frequency $\omega$ in the rest frame of the spacecraft, the
condition for pair creation is

$\omega\mathchar
13374\relax\frac{{m_{e}\\!}^{4}\,c^{6}\delta^{3}}{8\pi^{3}e^{2}\mathchar
1406\relax^{3}}=\frac{{E_{S}}^{2}}{8\pi\mathchar
1406\relax}\,\delta^{3}\approx 7.8\times
10^{51}\mbox{Hz}\times[\delta(\mbox{cm})]^{3}$ (6)

Again with the central CMB frequency of 160GHz, the blueshift factor to
transform the CMB photon to the rest frame of the spacecraft is
$\omega=\gamma(1+v/c)\times 10^{12}\,\mbox{Hz}$ as before. So the condition
for nucleus-mediated pair creation in terms of spacecraft speed and $\delta$
is:

$\gamma\left(1+\frac{v}{c}\right)\mathchar 13374\relax 7.8\times
10^{39}[\delta(\mbox{cm})]^{3}$ (7)

Let us assume the skin-depth $\delta$ for high frequency radiation is on the
order of a micron: $\delta=10^{-4}\,\mbox{cm}$. This gives

$\gamma\left(1+\frac{v}{c}\right)\mathchar 13374\relax 7.8\times 10^{27}$ (8)

corresponding to a $\gamma$ of about $\gamma\approx 4\times 10^{27}$, or
$v/c\approx 1-3\times 10^{-56}$.

Even if the speed of spacecraft is below the limits required for pair
production, simply scattering the microwave photons creates drag. How big is
that drag? To calculate this in order of magnitude, assume each photon is
scattered by reflecting directly back as it hits the spacecraft. This
corresponds to a momentum transfer to each photon of about

$2\frac{\mathchar
1406\relax\omega}{c}=\gamma\left(1+\frac{v}{c}\right)\,6.6\times
10^{-26}\;\mbox{g cm/sec}=\gamma\left(1+\frac{v}{c}\right)\,6.6\times
10^{-31}\;\mbox{Nt sec}$ (9)

(as measured in the rest frame of the spacecraft). Taking into account all
photons scattering off the frontal cross section per unit (spacecraft) proper
time, we obtain a viscous drag force on the spacecraft of about

$A(\mbox{cm}^{2})\gamma^{2}\left(1+\frac{v}{c}\right)\;8\times
10^{-18}\mbox{Nt}$ (10)

This will become a significant obstacle at speeds for which $\gamma\mathchar
13374\relax 10^{8}$ or $v/c\mathchar 13351\relax 1-5\times 10^{-17}$, which is
just below the velocity threshold for nucleus-mediated pair production.

We summarize various effects caused by interactions with matter and CMB
photons in the table below:

$\gamma$ | Speed |  Effect
---|---|---
1.3 | .69$c$ |  pion production in nucleon-nucleon collisions
9.1 | .994$c$ |  anti-proton production in nucleon-nucleon collisions
$10^{8}$ | $(1-5\times 10^{-17})c$ |  nucleus-mediated pair production in collisions with CMB photons
$4\times 10^{27}$ | $(1-3\times 10^{-56})c$ |  pair production in collisions with CMB photons via Schwinger vacuum breakdown

Lorentz-transformed Appearance of CMB

We now turn to examine how the CMB photon distribution would appear to an
observer in the rest frame of the spacecraft. In the rest frame of the CMB,
the photon distribution is thermal with the Planck density

$n(\vec{k})\,d^{3}k=\frac{1}{4\pi^{3}}\frac{1}{e^{\mathchar
1406\relax\omega/k_{B}T}-1}\,d^{3}k$ (11)

where $\omega=c|\vec{k}|$, $k_{B}$ is Boltzmann’s constant, and the left hand
side denotes the number of photons per unit volume with wave vectors lying in
a cell of size $d^{3}k$ centered around $\vec{k}$. From Eq. (11) we can
deduce, among other things, that the number of photons per unit volume with
frequencies between $\omega$ and $\omega+d\omega$ is

$n(\omega)\,d\omega=\frac{1}{\pi^{2}c^{3}}\frac{1}{e^{\mathchar
1406\relax\omega/k_{B}T}-1}\,\omega^{2}\,d\omega$ (12)

Equation (12) is the more familiar form of the Planck (Bose-Einstein)
distribution giving rise to, for example, the total number of photons per unit
volume

$\int_{0}^{\infty}n(\omega)\,d\omega=\frac{2\,{k_{B}\\!}^{3}\,T^{3}}{\pi^{2}c^{3}\mathchar
1406\relax^{3}}\,\zeta(3)\approx 412\;\mbox{photons/cm}^{3}$ (13)

where $T\approx 2.73\,^{\circ}\\!K$ for the CMB, and $\zeta$ is the Riemann
zeta function.

Consider now a CMB photon with 4-momentum vector

$p=\mathchar
1406\relax\omega(\partial_{t}+n_{x}\partial_{x}+n_{y}\partial_{y}+n_{z}\partial_{z})$
(14)

in the rest frame of the CMB spanned by the orthonormal basis
$(\partial_{t},\partial_{x},\partial_{y},\partial_{z})$. We will use units in
which $c=1$ unless we explicitly switch back to standard units for physical
interpretation, so here $\vec{n}$ is a unit 3-vector in the direction of the
photon’s propagation. The rest frame of the spacecraft is spanned by the basis
vectors

$\displaystyle e_{0}$ $\displaystyle=$
$\displaystyle\gamma(\partial_{t}+v\partial_{x})$ $\displaystyle e_{1}$
$\displaystyle=$ $\displaystyle\gamma(v\partial_{t}+\partial_{x})$
$\displaystyle e_{2}$ $\displaystyle=$ $\displaystyle\partial_{y}$
$\displaystyle e_{3}$ $\displaystyle=$ $\displaystyle\partial_{z}$ (15)

The same photon is seen in the rest frame of the spacecraft as having the
4-momentum

$p=\mathchar
1406\relax\omega\gamma(1-vn_{x})\left[e_{0}+\frac{n_{x}-v}{1-vn_{x}}\,e_{1}+\frac{n_{y}}{\gamma(1-vn_{x})}\,e_{2}+\frac{n_{z}}{\gamma(1-vn_{x})}\,e_{3}\right]$
(16)

which is the 4-vector Eq. (14) acted on by the Lorentz transformation matrix

$\Lambda=\begin{bmatrix}\gamma&-\gamma v&0&0\\\ -\gamma v&\gamma&0&0\\\
0&0&1&0\\\ 0&0&0&1\end{bmatrix}$ (17)

Similarly, a photon with 4-momentum vector

$p=\mathchar 1406\relax\omega(e_{0}+m_{x}e_{1}+m_{y}e_{2}+m_{z}e_{3})$ (18)

in the rest frame of the spacecraft is seen as having 4-momentum

$p=\mathchar
1406\relax\omega\gamma(1+vm_{x})\left[\partial_{t}+\frac{m_{x}+v}{1+vm_{x}}\,\partial_{x}+\frac{m_{y}}{\gamma(1+vm_{x})}\,\partial_{y}+\frac{m_{z}}{\gamma(1+vm_{x})}\,\partial_{z}\right]$
(19)

in the rest frame of the CMB, which is just the 4-vector Eq. (18) acted on by
the inverse of $\Lambda$

$\Lambda^{-1}=\begin{bmatrix}\gamma&\gamma v&0&0\\\ \gamma v&\gamma&0&0\\\
0&0&1&0\\\ 0&0&0&1\end{bmatrix}$ (20)

The distribution function for photon density transforms under a coordinate
change according to the general identity

$\hat{n}(\vec{k}^{\prime})=n[\vec{k}(\vec{k}^{\prime})]\left|\frac{\partial\vec{k}}{\partial\vec{k}^{\prime}}\right|$
(21)

where $\hat{n}$ is the density in the Lorentz transformed frame, and the right
hand side involves the coordinates $\\{k_{x},k_{y},k_{z}\\}$ expressed as
functions of the transformed coordinates
$\\{k_{x}^{\prime},k_{y}^{\prime},k_{z}^{\prime}\\}$ where the last term is
the Jacobian determinant. According to Eq. (19), we have

$\vec{k}(\vec{k}^{\prime})=\left[\gamma(k_{x}^{\prime}+v\omega^{\prime}),k_{y}^{\prime},k_{z}^{\prime}\right]$
(22)

and

$\omega^{2}=k_{x}^{2}+k_{y}^{2}+k_{z}^{2}=\gamma^{2}(\omega^{\prime}+vk_{x}^{\prime})^{2}$
(23)

The Jacobian determinant is easily computed:

$\left|\frac{\partial\vec{k}}{\partial\vec{k}^{\prime}}\right|=\gamma\left(1+\frac{vk_{x}^{\prime}}{\omega^{\prime}}\right)$
(24)

Figure 1:  How the cosmic microwave background would appear to a forward
facing observer located in a relativistic spacecraft. The plots are for
increasing spacecraft speeds (relative to the CMB) as measured by the Lorentz
factor $\gamma$. The appearance is approximately the same as that of a thermal
distribution with latitudinally-varying temperature. If we use Wien’s
displacement law, the apparent color of the CMB is the wavelength at which the
spectral radiance of the Planck distribution reaches its peak. The visible
colors form concentric rings confined to a narrow cone in the polar angle
$\theta$, with the zenith $\theta=0$ fixed along the forward $x$-direction.
The plots depict the apparent “celestial sky” parameterized with the $y,z$
coordinates of a unit sphere whose north pole lies on the $x$-axis. The first
light from blue-shifted CMB photons is visible in red at around $\gamma=700$.
For an explanation of the variation in shape and size of the color rings with
increasing $\gamma$ see the discussion following Eq. (27).

Combining Eqs. (21)–(24) with Eq. (11) and dropping the primes on transformed
quantities for notational convenience, we obtain the CMB photon distribution
as seen in the rest frame of the spacecraft

$\hat{n}(\vec{k})=\frac{1}{4\pi^{3}}\frac{\gamma\left(1+\frac{vk_{x}}{\omega}\right)}{e^{\mathchar
1406\relax\gamma(\omega+vk_{x})/k_{B}T}-1}$ (25)

where $\omega=c\sqrt{k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}$.

It is an interesting exercise to plot the appearance of the CMB to a forward
facing observer located in the spacecraft (Fig. 1). For qualitative
description, we can neglect the Jacobian factor in the numerator of Eq. (25),
and notice that for photons moving opposite to the “look” direction the wave
vector obeys $k_{x}=-(\omega/c)\cos\theta$, where $\theta$ measures the angle
between the look direction and the forward ($x$-) direction. Thus to the
forward-facing observer the CMB appears to have a thermal distribution with
latitudinally-varying temperature

$T(\theta)=\frac{T}{\gamma\left(1-\frac{v}{c}\cos\theta\right)}=T\frac{\sqrt{1-\frac{v^{2}}{c^{2}}}}{\left(1-\frac{v}{c}\cos\theta\right)}$
(26)

If we use Wien’s displacement law and approximate the apparent color of the
CMB with the wavelength at which the spectral radiance of the Planck
distribution reaches its peak, then the color visible at polar angle $\theta$
is given by the wavelength

$\lambda(\theta)=\frac{b}{T(\theta)}=\frac{b}{T}\frac{\left(1-\frac{v}{c}\cos\theta\right)}{\sqrt{1-\frac{v^{2}}{c^{2}}}}=\frac{b}{T}\left(\gamma-\sqrt{\gamma^{2}-1}\cos\theta\right)$
(27)

where $b$ is Wien’s displacement constant equal to $2.9\times
10^{-3}\,{}^{\circ}\rm{K}\,\rm{m}$. As shown in Fig.1, the visible frequencies
are confined to a narrow cone in $\theta$ around the forward direction, and
the constant color temperatures form rings around $\theta=0$ whose shape and
size vary with increasing $\gamma$. This is closely related to the phenomeneon
of (relativistic) aberration, which describes how the appearance of an optical
source varies in a geometric, frequency-independent way under a Lorentz
transformation aberr . To understand this variation better, it may help to
examine the angle $\theta$ at which the relativistic Doppler shift of the
incoming CMB photon is zero (i.e., no shift in frequency). The condition for
this to happen is

$\gamma\left(1-\frac{v}{c}\cos\theta\right)=\gamma-\sqrt{\gamma^{2}-1}\cos\theta=1$
(28)

The solution $\theta_{0}$ to Eq. (28) is

$\theta_{0}=\arccos\left(\frac{\gamma-1}{\sqrt{\gamma^{2}-1}}\right)=\arccos\left(\sqrt{\frac{{\gamma-1}}{{\gamma+1}}}\right)$
(29)

The CMB photons appear blue-shifted to our forward-facing observer for angles
$\theta$ less than $\theta_{0}$ and red-shifted for angles $\theta\mathchar
13374\relax\theta_{0}$. For small velocities ($\gamma$ close to 1),
$\theta_{0}$ is near $\pi/2$ as would be expected. However, just as in
relativistic aberration, in the limit of large velocities (large $\gamma$)
$\theta_{0}$ approaches zero:

$\theta_{0}=\sqrt{\frac{2}{\gamma}}-\frac{1}{3\sqrt{2}\;\gamma^{\nicefrac{{3}}{{2}}}}+O\left(\frac{1}{\gamma^{\nicefrac{{5}}{{2}}}}\right)$
(30)

which explains the progressive narrowing at very large $\gamma$ values of the
visible ring in Fig. 1.

Signature of CMB Scattered from a Relativistic Spacecraft

The possibility of detecting radiation associated with distant relativistic
spacecraft has been discussed in the literature before zubrin ; viewing ;
garcia . These discussions mostly focus on detecting radiation from spacecraft
engines or light from nearby stars reflecting off the spacecraft. Our approach
is different in that we do not speculate on possible propulsion technologies
but are interested in how a large relativistic object would interact with the
interstellar/intergalactic medium and mainly with the CMB radiation. As a
baryonic spacecraft travels at relativistic speeds it will interact with the
CMB through scattering to cause a frequency shift that could be detectable on
Earth with current technology.

Now we consider the characteristic signature of this reflected CMB radiation;
specifically, CMB photons reflected by a surface area element of area $A$
moving with velocity $v$ in the $x$ direction. Let the normal to the surface
element $A$ be described by a unit 3-vector $\hat{s}$ in the spacecraft
restframe. Then every photon impacting $A$ is reflected back according to

$\vec{k}\longrightarrow\vec{k}-2(\vec{k}\cdot\hat{s})\,\hat{s}$ (31)

From a far-field perspective the area element $A$ would appear as a point
source. We set up a coordinate system in the CMB rest frame where $A$ is at
$\vec{x}=\vec{v}t=vt\hat{x}$ and the observation point is at $\vec{r}$. Notice
that the time coordinate is chosen in such a way that the spacecraft’s
worldline passes through the coordinate origin at time $t=0$. The radiation
detected at $\vec{r}$ at time $t$ would have been emitted from the spacecraft
at the retarded time $\sigma$ which is the time at the intersection between
the past null cone of the event $(\vec{r},t)$ and the worldline of the
spacecraft. Hence the equation for $\sigma$ is

$c(t-\sigma)=\sqrt{(x-v\sigma)^{2}+y^{2}+z^{2}}$ (32)

The solution to Eq. (32) is

$\sigma^{{\rm ret}\choose{\rm
adv}}=\gamma^{2}\left(t-\frac{vx}{c^{2}}\mp\frac{1}{c\gamma}\sqrt{\gamma^{2}(x-vt)^{2}+y^{2}+z^{2}}\right)$
(33)

where the retarded (advanced) solution represents intersection with the past
(future) null cone of the observation point $(\vec{r},t)$. We will only use
the retarded solution

$\sigma=\sigma(\vec{r},t)\equiv\sigma^{({\rm
ret})}(\vec{r},t)=\gamma^{2}\left(t-\frac{vx}{c^{2}}-\frac{1}{c\gamma}\sqrt{\gamma^{2}(x-vt)^{2}+y^{2}+z^{2}}\right)$
(34)

Only reflected photons with wave vectors parallel to $(\vec{r}-\vec{v}\sigma)$
will be detected. There are two Lorentz transformations involved: First to
find the photon distribution in the spacecraft rest frame according to Eq.
(25), and second to transform the distribution reflected according to Eq. (31)
back into the CMB rest frame, which is accomplished by the same formulas as in
Eqs. (16)–(25). According to Eqs. (16)–(25) and Eq. (31), the spectral
intensity of the reflected radiation in the CMB rest frame at position
$\vec{r}$ and at time $t$ would be (measured in number of photons of frequency
$\omega$ per unit area per unit time per unit frequency)

$I(\vec{r},t,\omega)=\frac{cA}{4\pi|\vec{r}-\vec{v}\sigma|^{2}}\frac{\omega^{2}}{\pi^{2}c^{3}}\frac{\gamma^{2}\left[1+\left(\frac{\vec{v}}{c|\vec{r}-\vec{v}\sigma|}\cdot(\vec{r}-\vec{v}\sigma-2[(\vec{r}-\vec{v}\sigma)\cdot\hat{s}]\hat{s}\right)\right]^{2}}{\exp\left[\frac{\gamma^{2}\mathchar
1406\relax\omega}{k_{B}T}\left[1+\left(\frac{\vec{v}}{c|\vec{r}-\vec{v}\sigma|}\cdot(\vec{r}-\vec{v}\sigma-2[(\vec{r}-\vec{v}\sigma)\cdot\hat{s}]\hat{s}\right)\right]^{2}\right]-1}$
(35)

where $\sigma$ is given by Eq. (34) and $\vec{v}=v\hat{x}$. For a general
reflecting surface $S$ made up of many small surface elements, the total
intensity can be computed as a surface integral of Eq. (35)

$I_{T}(\vec{r},t,\omega)=\frac{\omega^{2}}{4\pi^{3}c^{2}}\oint_{S}\frac{1}{|\vec{r}-\vec{v}\sigma|^{2}}\frac{\gamma^{2}\left[1+\left(\frac{\vec{v}}{c|\vec{r}-\vec{v}\sigma|}\cdot(\vec{r}-\vec{v}\sigma-2[(\vec{r}-\vec{v}\sigma)\cdot\hat{s}]\hat{s}\right)\right]^{2}}{\exp\left[\frac{\gamma^{2}\mathchar
1406\relax\omega}{k_{B}T}\left[1+\left(\frac{\vec{v}}{c|\vec{r}-\vec{v}\sigma|}\cdot(\vec{r}-\vec{v}\sigma-2[(\vec{r}-\vec{v}\sigma)\cdot\hat{s}]\hat{s}\right)\right]^{2}\right]-1}\;dA$
(36)

We will limit our attention in this paper to reflection from a single surface
element $A$ as given in Eq. (35). Using Wien’s displacement law, the peak
frequency of the reflected radiation Eq. (35) is given by

$\nu_{\rm
max}(\vec{r},t)=\frac{c}{b\gamma^{2}}\frac{T}{\left[1+\left(\frac{\vec{v}}{c|\vec{r}-\vec{v}\sigma|}\cdot(\vec{r}-\vec{v}\sigma-2[(\vec{r}-\vec{v}\sigma)\cdot\hat{s}]\hat{s}\right)\right]^{2}}$
(37)

corresponding to a Planck distribution with effective temperature

$T_{\rm
eff}(\vec{r},t)=\frac{T}{\gamma^{2}\left[1+\left(\frac{\vec{v}}{c|\vec{r}-\vec{v}\sigma|}\cdot(\vec{r}-\vec{v}\sigma-2[(\vec{r}-\vec{v}\sigma)\cdot\hat{s}]\hat{s}\right)\right]^{2}}$
(38)

On the other hand, the relative intensity profile is determined by the
dimensionless multiplicative factor

$i(\vec{r},t,\omega)=\frac{A}{|\vec{r}-\vec{v}\sigma|^{2}}\gamma^{2}\left[1+\left(\frac{\vec{v}}{c|\vec{r}-\vec{v}\sigma|}\cdot(\vec{r}-\vec{v}\sigma-2[(\vec{r}-\vec{v}\sigma)\cdot\hat{s}]\hat{s}\right)\right]^{2}$
(39)

Figure 2:  Some examples of signal signatures of CMB scattering from
relativistic spacecraft. In all plots the observation point is located at
$x=10$ lightyears, $y=z=1$ lightyear. The time axis origin is at $t=8$ years
(spacecraft is moving along the $x$-direction and is at the coordinate origin
at time $t=0$). In the upper row, we plot the apparent (effective) temperature
of the signal from a surface element orthogonal to the $x$-direction (the
direction of motion) for spacecraft velocities $v=0.95\,c$ and $v=0.99\,c$,
respectively. In the lower row, the plots are for a surface element tilted
$10^{\circ}$ from orthogonal to the $x$-direction but still parallel to the
$z$-direction, with the same two velocities. The red curves are plots of the
multiplicative factor Eq. (39) in arbitrary units, showing the relative
intensity profile of the received radiation as a function of time. The
decrease in temperature with increasing $\gamma$ in the second row of plots
might seem counter-intuitive at first; but it is due to the same phenomenon of
relativistic aberration discussed above [cf. Eq. (30)]; namely, the
$10^{\circ}$ tilt in the reflecting surface is enough to put the reflected
beam outside the cone of blue-shift for the $\gamma$ value corresponding to
$v=0.99\,c$.

Some examples of signal signatures resulting from Eqs. (38)–(39) are shown in
Figs. 2 and 3. The salient features of the signal are a rapid drop in
temperature accompanied by a rapid rise in intensity, along with the motion of
the source with respect to a reference frame fixed to distant quasars, which
should be observable. No natural source of THz to infrared radiation in the
known Universe is likely to have time variability of the kind depicted in
Figs. 2–3. Note that increasing velocity (increasing $\gamma$) manifests
itself in a steepening of the cooling curve, along with an overall increase in
temperature for a surface element oriented orthogonal to the direction of
motion.

Figure 3:  A further example of signal signature for CMB scattering from
relativistic spacecraft. The spacecraft speed in this case is a relatively
small $v=0.64\,c$ ($\gamma=1.3$). As in Fig. 2, the observation point is
located at $x=10$ lightyears, $y=z=1$ lightyear. The time axis origin is at
$t=8$ years (spacecraft is moving along the $x$-direction and is at the
coordinate origin at time $t=0$). This time we plot the peak frequency of the
signal; on the left from a surface element orthogonal to the $x$-direction,
and on the right from a surface element tilted $10^{\circ}$ from that
orthogonal but still parallel to the $z$-direction. The red curve is again a
plot of the multiplicative factor Eq. (39) in arbitrary units, showing the
relative intensity profile of the received radiation as a function of time.
The salient features recognizable in Fig. 2 are much less prominent, owing to
the comparatively small value of $\gamma$ in this case.

Finally, we note the similarity of CMB scattering from macroscopic
relativistic objects to the Sunyaev-Zeldovich (S-Z) effect in cosmology. The
S-Z effect is the spectral distortion of the CMB caused by inverse Compton
scattering of microwave background photons from energetic intracluster thermal
electrons (intracluster referring to clusters of galaxies) SZ1 ; SZ2 . S-Z
scattering has been observed birkinson ; carlstrom . The intracluster
electrons have temperatures of order 3 to 15 keV, which correspond to gamma
factors of $\gamma\approx 1.01-1.03$, which correspond to thermal velocities
of $v/c\approx 0.15-0.25$. Therefore, while these electrons are relativistic,
they are far from the ultra-relativistic speeds we are interested in. However,
the observed S-Z effect in cosmology is qualitatively similar to what we
discuss above.

On the other hand, merger shocks in the collisions of clusters of galaxies can
produce non-thermal electrons at energies from 300MeV ($\gamma\approx 300$) up
to $10^{6}$MeV ($\gamma\approx 2\times 10^{6}$), so these electrons are more
relevant to our discussion mergershocks . In fact, these ultra-relativistic
electrons, protons, and heavier ions accelerated by merger shocks loose energy
via inverse Compton, pair production, and other effects mediated by collisions
with the cosmic microwave photons, much as we discussed for the
ultrarelativistic spaceship example. The S-Z effect in this case, which is
exactly the same as the spaceship example except the relativistic object is a
single particle, should be similar provided it is observable.

Summary

We have explored the physical interaction of a relativistic spacecraft with
the interstellar and intergalactic mediums. Our main discussion focused on the
interaction with the CMB where very little information currently exists and
seems to have been overlooked. Central to our discussion was not how to obtain
relativistic speeds but the consequences of traveling that fast. We discussed
the interaction with baryonic matter in terms of high speed collisions to
present a complete picture. In general one can imagine the same interactions
that occur in a particle accelerator to occur between relativistic spacecraft
and interstellar matter. Our assumption that matter-matter interactions can be
dealt with when civilization can build relativistic spacecraft may prove false
and may be a barrier that will prevent space travel with a large $\gamma$.

We looked at two special reference frames, the spacecraft frame, and the rest
frame of the CMB to understand how the CMB is distorted or acquires an
aberration from those view points. The scattering of the CMB from a
relativistic spacecraft is very similar to the Sunyaev-Zedovich effect where
the inverse Compton scattering instead occurs at a macroscopic level. Our
calculation for what an observer on earth could detect predicts a very unusual
signature that is unlikely to be caused by any naturally occurring object in
the known universe. This result is independent of propulsion technology, but
the ability to detect the signal from Earth depends on available detector
technologies. We are currently working to predict how far can we see this
signature given our current capability.

## References

  * (1) P. A. Tipler & R. A. Llewellyn, Modern Physics. W. H. Freeman and Co., New York, Fifth Edition (2008).
  * (2) S. Weinberg, Gravitation and Cosmology. Wiley, New York (1972).
  * (3) D. J. Griffiths, Introduction to Elementary Particles. John Wiley and Sons (1987).
  * (4) D. H. Perkins, Introduction to High Energy Physics. Addison-Wesley Press (1989).
  * (5) G. D. Scott and H. J. van Driel, Am. J. Phys. 38, 971–977 (1970).
  * (6) R. Zubrin, Detection of Extraterrestrial Civilizations via the Spectral Signature of Advanced Interstellar Spacecraft, in Progress in the Search for Extraterrestrial Life, ASP Conference Series, 74, 487 (1995).
  * (7) D. R. J. Viewing, C. Horswell, E. W. Palmer, JBIS 30, 99–104 (1977).
  * (8) J. C. Garcia-Escarfin & P. Chamorro-Posada, arXiv.1203:3980 (2012).
  * (9) R. A. Sunyaev & Y. B. Zeldovich, Comments Astrophys. Space Phys. 2, 66 (1970).
  * (10) R. A. Sunyaev & Y. B. Zeldovich, Comments Astrophys. Space Phys. 4, 173 (1972).
  * (11) M. Birkinshaw, Phys. Reports 310, 97–196 (1999).
  * (12) J. E. Carlstrom, M. Joy, L. Grego, Astrophys. Journal 456, L75–L78 (1996).
  * (13) M. Markevitch & A. Vikhlinin, Phys. Reports 443, 1 (2007).


∎

11institutetext: Chengjing Wang 22institutetext: School of Mathematics,
Southwest Jiaotong University 22email<EMAIL_ADDRESS>33institutetext: Xile Zhao 44institutetext: School of Mathematical
Sciences/Research Center for Image and Vision Computing, University of
Electronic Science and Technology of China 44email<EMAIL_ADDRESS>55institutetext: Qingsong Wang 66institutetext: Department of Mathematics,
National University of Defense Technology 66email<EMAIL_ADDRESS>77institutetext: Zepei Ma 88institutetext: School of Mathematics, Southwest
Jiaotong University 88email<EMAIL_ADDRESS>99institutetext: Peipei Tang
(corresponding author) 1010institutetext: Hangzhou City University, 1010email:
<EMAIL_ADDRESS>
# An inexact proximal majorization-minimization Algorithm for remote sensing
image stripe noise removal

Chengjing Wang  Xile Zhao  Qingsong Wang  Zepei Ma  Peipei Tang

(Received: date / Accepted: date)

###### Abstract

The stripe noise existing in remote sensing images badly degrades the visual
quality and restricts the precision of data analysis. Therefore, many
destriping models have been proposed in recent years. In contrast to these
existing models, in this paper, we propose a nonconvex model with a DC
function (i.e., the difference of convex functions) structure to remove the
strip noise. To solve this model, we make use of the DC structure and apply an
inexact proximal majorization-minimization algorithm with each inner
subproblem solved by the alternating direction method of multipliers. It
deserves mentioning that we design an implementable stopping criterion for the
inner subproblem, while the convergence can still be guaranteed. Numerical
experiments demonstrate the superiority of the proposed model and algorithm.

###### Keywords:

Remote sensing group sparsity proximal majorization-minimization algorithm
alternating direction method of multipliers

###### MSC:

65K05 90C26 90C90

##  1 Introduction

Nowadays, remote sensing images are playing an important role in the
environment, agriculture, biology, mineralogy, and so on. One may see
Richards22 , ZhouCGLW22 , HuangTTZHCR22  for more details. Due to the
limitation of sensors, remote sensing images are inevitably corrupted by
noise, which deteriorates the quality of remote sensing images and hinders
their subsequent applications (e.g., classification HuangH14  and detection
DuZTZ13 ). Remote sensing image denoising (RSID), which removes the noise
while preserving the details and edges, is one of the fundamental problems in
remote sensing image processing and analysis Richards22 .

In this work, we mainly focus on the stripe noise scenario. The degradation
process BoualiL11 , ShenZ09 , can be mathematically formulated as

$\displaystyle\begin{aligned} f(x,y)=u(x,y)+s(x,y),\end{aligned}$ (1)

where $f(x,y),u(x,y)$ and $s(x,y)$ represent the degraded image from the
detector, the potential fringeless image, and the fringe component at $(x,y)$,
respectively. In order to discuss the numerical algorithm, the matrix vector
form of (1) can be rewritten as follows

$\displaystyle\begin{aligned} f=u+s,\end{aligned}$ (2)

where $f,u$ and $s$ represent the vectorized discrete versions of
$f(x,y),u(x,y)$ and $s(x,y)$, respectively.

  * (1)

Directionality: consider the gradient of fringe components in horizontal and
vertical directions. Here, we use the convex function 1-norm as the sparse
regularization and give

$R_{1}(s)=\|\nabla_{y}s\|_{1,1},$

where $\nabla_{y}$ denotes a linear first order difference operator in the
vertical direction.

The horizontal gradient of the desired image $u$ should be smooth. Therefore,
the one-way total change BoualiL11  of the desired image is considered to be

$R_{2}(s)=\|\nabla_{x}u\|_{1}=\|\nabla_{x}f-\nabla_{x}s\|_{1,1},$

where $\nabla_{x}$ denotes the first order difference operator in the
horizontal direction.

  * (2)

Structural properties: moreover, different from the random noise, the fringe
component shows a special columnar structure. Therefore, the regular term
ChenHDZW  is designed as

$R_{3}(s)=\|s\|_{2,1},$

where
$\|s\|_{2,1}=\sum\limits_{j=1}^{n}\big{(}\sum\limits_{i=1}^{m}s_{i,j}^{2}\big{)}^{\frac{1}{2}},\forall
s\in\mathbb{R}^{m\times n}$.

Based on the above analysis, the stripe component has significant
directionality and structure. By combining $R_{1},R_{2}$ and $R_{3}$, Chen et
al. ChenHDZW  proposed the fringe noise removal model as follows

$\displaystyle\begin{aligned}
\min_{s}\;&\\{\lambda_{1}\|\nabla_{y}s\|_{1,1}+\lambda_{2}\|\nabla_{x}f-\nabla_{x}s\|_{1,1}+\lambda_{3}\|s\|_{2,1}\\},\end{aligned}$
(3)

where $\lambda_{1},\lambda_{2}$ and $\lambda_{3}$ are three positive
regularization parameters that balance the three terms.

In order to simplify the sign, we abbreviate
$p_{1}(s)=\lambda_{1}\|\nabla_{y}s\|_{1,1}$,
$p_{2}(u)=\lambda_{2}\|\nabla_{x}f-\nabla_{x}u\|_{1,1}$ and
$p_{3}(v)=\lambda_{3}\|v\|_{2,1}=\lambda_{3}\|h(v)\|_{1}$ with
$(h(v))_{i}:=\|v(:,i)\|$. So (3) is rewritten as follows

$\displaystyle\begin{aligned}
\min_{s}\;&\\{p_{1}(s)+p_{2}(s)+p_{3}(s)\\}.\end{aligned}$ (4)

For the convenience of discussion, each stripe line is regarded as a column.
If the stripes are horizontal, rotate them so that the stripe lines are
vertical. When the fringe component is extracted from the fringe image, the
final denoising image can be estimated by the following formula

$u=f-s.$

Fan and Li FanL2017  proposed a nonconvex sparsity function called SCAD as a
surrogate of the $\ell_{0}$ function, which has been used as the regularizer
to avoid the model overfitting in high-dimensional statistical learning. It
can achieve a sparse estimation with fewer measurements, faster convergence
and furthermore is more robust against noises. In this paper, we apply the
SCAD function to substitute the $\ell_{1}$ function part in the model (4) in
order to remove the noise more well. Hence the considered model becomes a
nonconvex model. In the computational aspect, we adopt an inexact proximal
majorization-minimization (PMM) algorithm to solve the nonconvex model with
each inner subproblem solved by the dual alternating direction method (dADMM).
Although Tang et al. TangWST  proposed the PMM algorithm to solve the
nonconvex square-root-loss regression problems and the inner subproblem is
also solved inexactly by the semismooth Newton method, the stopping criterion
for the inner subproblem is not implementable and remains on the theoretical
level. Tang, Wang and Jiang TangWJ  proposed the proximal-proximal
majorization-minimization (PPMM) algorithm to solve the nonconvex rank
regression problems with each inner subproblem solved by the semismooth Newton
method. They proposed a stopping criterion for the inner subproblem, however,
the convergence analysis is relatively easier since the considered problem is
polyhedral and the inner subproblem is an unconstrained one. Numerical
experiments demonstrate that the modified model and the proposed algorithm is
more competitive.

The remaining parts of this paper are organized as follows. In Section 2, we
introduce some basic concepts and preliminary results. In Section 3, we set up
the model problem and present the details of the algorithm. In Section 4, we
analyze the convergence of the proposed algorithm. In Section 5, we implement
the numerical experiments to compare our model with the original convex model
and also compare our algorithm with the existing algorithms. We conclude our
paper in Section 6.

##  2 Preliminaries

In this section, we introduce some preliminaries that will be used in this
paper.

Let $\mathbb{X}$ and $\mathbb{Y}$ be finite dimensional Hilbert spaces. For
any convex function $f:\mathbb{X}\to(-\infty,+\infty]$, the conjugate function
of $f$ is defined as

$f^{*}(y):=\sup_{x\in\textrm{dom}(f)}\;\Big{\\{}\langle
y,x\rangle-f(x)\Big{\\}}.$

For a closed proper convex function $f:\mathbb{X}\to(-\infty,+\infty]$ and a
parameter $\sigma>0$, the proximal mapping $\textrm{Prox}_{\sigma f}$ is
defined as

$\textrm{Prox}_{\sigma
f}(y):=\mathop{\arg\min}_{x}\Big{\\{}f(x)+\frac{1}{2\sigma}\|y-x\|^{2}\Big{\\}},\forall\,y\in\mathbb{X}.$

The proximal mapping $\mbox{Prox}_{\sigma f}$ is single-valued and continuous
with the following Moreau identity (see e.g.,(Rockafellar70, , Theorem 31.5))
holds

$\displaystyle\textrm{Prox}_{\sigma
f}(x)+\sigma\textrm{Prox}_{\frac{1}{\sigma}f^{*}}(\frac{1}{\sigma}x)=x.$

A multifunction $\mathcal{F}:\mathbb{X}\rightrightarrows\mathbb{Y}$ is locally
upper Lipschitz continuous at $x\in\mathbb{X}$ if there exist a parameter
$\kappa$ which is independent of $x$ and a neighbourhood $\mathcal{U}$ of $x$
such that
$\mathcal{F}(y)\subseteq\mathcal{F}(x)+\kappa\|y-x\|\mathbb{B}_{\mathbb{X}}$
holds for any $y\in\mathcal{U}$, where $\mathbb{B}_{\mathbb{X}}$ is a unit
ball of the space $\mathbb{X}$. The multifunction $\mathcal{F}$ is said to be
piecewise polyhedral if its graph
$\mbox{gph}\mathcal{F}:=\\{(x,y)\,|\,y\in\mathcal{F}(x)\\}$ is the union of
finitely many polyhedral convex sets.

The Kurdyka-Łojasiewicz (KŁ) property (see e.g., bolte-pauwels2016 ) plays a
central role in the convergence analysis. For further discussion, we introduce
this concept below.

###### Definition 1

A proper lower semicontinuous function
$r:\mathbb{X}\rightarrow(-\infty,+\infty]$ is said to have the KŁ property at
$x\in\mbox{dom}(\partial r)$ if there exist $\eta\in(0,+\infty]$, a neighbour
$\mathcal{U}$ of $x$ and a continuous concave function
$\varphi:[0,\eta)\rightarrow[0,+\infty)$ satisfying
(1) $\varphi(0)=0$;
(2) $\varphi$ is continuous at 0 and continuously differentiable on
$(0,\eta)$;
(3) $\varphi^{\prime}(s)>0$, for all $0<s<\eta$
such that the KŁ inequality
$\varphi^{\prime}(r(x^{\prime})-r(x))\mbox{dist}(0,\partial r(x^{\prime}))\geq
1$ holds for any $x^{\prime}\in\mathcal{U}$ and
$r(x)<r(x^{\prime})<r(x)+\eta$. If $r$ satisfies the KŁ property at each point
of $\mbox{dom}(\partial r)$, then $r$ is called a KŁ function.

##  3 Problem setup and the algorithm

In this paper, we consider the following nonconvex model problem by employing
the SCAD function instead of (4)

$\displaystyle\min_{s}\;\Big{\\{}g(s):=p_{1}(s)-q_{1}(s)+p_{2}(s)-q_{2}(s)+p_{3}(s)-q_{3}(s)\Big{\\}},$
(5)

where
$q_{1}(s)=q_{\alpha,\lambda_{1}}(\nabla_{y}s),q_{2}(u)=q_{\alpha,\lambda_{2}}(\nabla_{x}u)$
and $q_{3}(v)=q_{\alpha,\lambda_{3}}(h(v))$, here

$\displaystyle q_{\alpha,\lambda}(x)=\sum_{i=1}^{n}q^{\rm
scad}(x_{i};\alpha,\lambda),\ q^{\rm
scad}(t;\alpha,\lambda)=\left\\{\begin{array}[]{ll}0,&\mbox{if}\quad|t|\leq\lambda,\\\
\frac{(|t|-\lambda)^{2}}{2(\alpha-1)},&\mbox{if}\quad\lambda\leq|t|\leq\alpha\lambda,\\\
\lambda|t|-\frac{\alpha+1}{2}\lambda^{2},&\mbox{if}\quad|t|>\alpha\lambda.\end{array}\right.$

In the above formula, we take $\alpha=3.7$ in both the theoretical analysis
and the numerical implementation. Note that the function
$q_{\alpha,\lambda}(x)$ is continuously differentiable with

$\displaystyle\frac{\partial q_{\alpha,\lambda}(x)}{\partial
x_{i}}=\left\\{\begin{array}[]{ll}0,&\mbox{if}\ |x_{i}|\leq\lambda,\\\
\frac{\mbox{sign}(x_{i})(|x_{i}|-\lambda)}{\alpha-1},&\mbox{if}\
\lambda<|x_{i}|\leq\alpha\lambda,\\\ \lambda\mbox{sign}(x_{i}),&\mbox{if}\
|x_{i}|>\alpha\lambda.\end{array}\right.$

###  3.1 The inexact PMM for the problem (5)

In this subsection, we present the framework of the inexact PMM. For
simplicity, we write the problem (5) as

$\displaystyle\min_{s}\;\Big{\\{}P(s)-Q(s)\Big{\\}},$ (8)

where $P(s)=p_{1}(s)+p_{2}(s)+p_{3}(s)$, $Q(s)=q_{1}(s)+q_{2}(s)+q_{3}(s)$.

Although problem (8) is nonconvex, it is in a DC form. Thus, we naturally
linearize the latter convex term $Q(s)$ to transform a nonconvex problem to a
convex one, add a proximal term to guarantee the strong convexity, then adopt
the sequential convexification approach to solve the problem. Now we present
the algorithm as follows.

Algorithm 1 The inexact PMM for the problem (8):

Choose $s^{0}$. For $k=0,1,2,\ldots$, iterate:

1.

Compute $s^{k+1}$

$\displaystyle\begin{aligned}
s^{k+1}&\approx\mathop{\arg\min}_{s}\;\left\\{P(s)-Q(s^{k})-\langle\nabla
Q(s^{k}),s-s^{k}\rangle+\frac{1}{2\widetilde{\sigma}_{k}}\|s-s^{k}\|^{2}\right\\}.\end{aligned}$
(9)

2.

If a stopping criterion is not satisfied, go to Step 1.

###  3.2 The dADMM for the subproblem (9)

In this subsection, we focus on how to solve the subproblem (9). In the
problem (9), for $q_{i}(x)\,(i=1,2,3)$, at the $k$th iteration we have
$q_{i}(x)\geq q_{i}(x^{k})+\langle g_{i}^{k},x-x^{k}\rangle$, where
$g_{i}^{k}:=\nabla q_{i}(x^{k})$. We denote
$\widehat{p}_{1}(s):=p_{1}(s)+\frac{1}{2\widetilde{\sigma}_{k}}\|s-s^{k}\|^{2}$
and $\widehat{p}_{3}(v):=\lambda_{3}\|h(v)\|-\langle g_{h}^{k},h(v)\rangle$,
where $g_{h}^{k}:=\nabla q_{3}(h(v^{k}))$. By introducing two variables $u$
and $v$, the problem (9) can be equivalently rewritten as

$\displaystyle\begin{aligned} \min_{s,u,v}\,\;&\widehat{p}_{1}(s)-\langle
g_{1}^{k},s\rangle+p_{2}(u)-\langle g_{2}^{k},u\rangle+\widehat{p}_{3}(v)\\\
\mbox{s.t.}\quad&f-s-u=0,\\\ &s-v=0.\end{aligned}$ (10)

The Lagrangian function associated with the problem (10) is given by

$\displaystyle l(s,u,v;x,y)$ $\displaystyle=$
$\displaystyle\widehat{p}_{1}(s)-\langle g_{1}^{k},s\rangle+p_{2}(u)-\langle
g_{2}^{k},u\rangle+\widehat{p}_{3}(v)+\langle x,u+s-f\rangle+\langle
y,v-s\rangle.$

Then the Karush-Kuhn-Tucker (KKT) condition for the problem (10) is

$\displaystyle f-s-u=0,\ s-v=0,\
s-\textrm{Prox}_{\widehat{p}_{1}}(s+g^{k}_{1}-x+y)=0,$ $\displaystyle
u-\textrm{Prox}_{p_{2}}(u+g^{k}_{2}-x)=0,\
v-\textrm{Prox}_{\widehat{p}_{3}}(v-y)=0.$ (11)

And the dual problem is

$\displaystyle\begin{aligned}
\min_{x,y,z,\hat{x},\hat{y}}\;&\widehat{p}_{1}^{*}(z)+p_{2}^{*}(\widehat{x})+\widehat{p}_{3}^{*}(\widehat{y})+\langle
x,f\rangle\\\ \mbox{s.t.}\quad&-x+y+g_{1}^{k}=z,\\\
&-x+g_{2}^{k}=\widehat{x},\\\ &-y=\widehat{y}.\end{aligned}$ (12)

Next, the augmented Lagrangian function associated with the problem (12) is

$\displaystyle L_{\sigma}(x,y,z,\widehat{x},\widehat{y};s,u,v)$
$\displaystyle=$
$\displaystyle\widehat{p}_{1}^{*}(z)+p_{2}^{*}(\widehat{x})+\widehat{p}_{3}^{*}(\widehat{y})+\langle
x,f\rangle+\frac{\sigma}{2}\|-x+y+g_{1}^{k}-z+\frac{1}{\sigma}s\|^{2}$
$\displaystyle+\frac{\sigma}{2}\|-x+g_{2}^{k}-\widehat{x}+\frac{1}{\sigma}u\|^{2}+\frac{\sigma}{2}\|-y-\widehat{y}+\frac{1}{\sigma}v\|^{2}$
$\displaystyle-\frac{1}{2\sigma}(\|s\|^{2}+\|u\|^{2}+\|v\|^{2}).$

Then we describe the details of how to solve the subproblems.

  * (1)

$(x,y)$-subproblem: Let $(\bar{x},\bar{y})$ denote the optimal solution of the
$(x,y)$-subproblem, i.e.,

$\displaystyle(\bar{x},\bar{y})$
$\displaystyle=\operatorname*{\mathop{\arg\min}}_{x,y}\Big{\\{}L_{\sigma}(x,y,z,\widehat{x},\widehat{y};s,u,v)\Big{\\}}.$

It is equivalent to solving the following linear system of equations

$\displaystyle\begin{bmatrix}2I&-I\\\ -I&2I\end{bmatrix}\begin{bmatrix}x\\\
y\end{bmatrix}$
$\displaystyle=\frac{1}{\sigma}\begin{bmatrix}s+u-f-\sigma(z+\widehat{x}-g^{k}_{1}-g^{k}_{2})\\\
v-s+\sigma(z-\widehat{y}-g^{k}_{1})\end{bmatrix}.$

  * (2)

$(z,\widehat{x},\widehat{y})$-subproblem: Let
$(\bar{z},\bar{\widehat{x}},\bar{\widehat{y}})$ be the optimal solution of the
$(z,\widehat{x},\widehat{y})$-subproblem. We have

$\displaystyle\bar{z}$ $\displaystyle=$
$\displaystyle\operatorname*{\mathop{\arg\min}}_{z}\;\Big{\\{}\widehat{p}_{1}^{*}(z)+\frac{\sigma}{2}\|-x+y+g_{1}^{k}-z{\color[rgb]{0,0,1}+}\frac{1}{\sigma}s\|^{2}\Big{\\}}$
$\displaystyle=$ $\displaystyle{\rm
Prox}_{\frac{1}{\sigma}\widehat{p}_{1}^{*}}(\frac{s}{\sigma}-x+y+g_{1}^{k}),$
$\displaystyle\bar{\widehat{x}}$ $\displaystyle=$
$\displaystyle\operatorname*{\mathop{\arg\min}}_{\widehat{x}}\;\Big{\\{}p_{2}^{*}(\widehat{x})+\frac{\sigma}{2}\|-x+g_{2}^{k}-\widehat{x}+\frac{1}{\sigma}u\|^{2}\Big{\\}}$
$\displaystyle=$ $\displaystyle{\rm
Prox}_{\frac{1}{\sigma}p_{2}^{*}}(\frac{u}{\sigma}-x+g_{2}^{k}),$
$\displaystyle\bar{\widehat{y}}$ $\displaystyle=$
$\displaystyle\operatorname*{\mathop{\arg\min}}_{\widehat{y}}\;\Big{\\{}\widehat{p}_{3}^{*}(\widehat{y})+\frac{\sigma}{2}\|-y-\widehat{y}+\frac{1}{\sigma}v\|^{2}\Big{\\}}$
$\displaystyle=$ $\displaystyle{\rm
Prox}_{\frac{1}{\sigma}\widehat{p}_{3}^{*}}(\frac{v}{\sigma}-y).$

Based on the Moreau identity, we have

$\displaystyle{\rm
Prox}_{\frac{1}{\sigma}\widehat{p}_{1}^{*}}(\frac{s}{\sigma}-x+y+g_{1}^{k})$
$\displaystyle=$ $\displaystyle\frac{1}{\sigma}(s-\sigma x+\sigma y+\sigma
g_{1}^{k}-{\rm Prox}_{\sigma\widehat{p}_{1}}(s-\sigma x+\sigma y+\sigma
g_{1}^{k})),$ $\displaystyle{\rm
Prox}_{\frac{1}{\sigma}p_{2}^{*}}(\frac{u}{\sigma}-x+g_{2}^{k})$
$\displaystyle=$ $\displaystyle\frac{1}{\sigma}(u-\sigma x+\sigma
g_{2}^{k}-{\rm Prox}_{\sigma p_{2}}(u-\sigma x+\sigma g_{2}^{k})),$
$\displaystyle{\rm
Prox}_{\frac{1}{\sigma}\widehat{p}_{3}^{*}}(\frac{v}{\sigma}-y)$
$\displaystyle=$ $\displaystyle\frac{1}{\sigma}(v-\sigma y-{\rm
Prox}_{\sigma\widehat{p}_{3}}(v-\sigma y)).$

Note that

$\displaystyle{\rm Prox}_{\sigma\widehat{p}_{1}}(s)$ $\displaystyle=$
$\displaystyle\operatorname*{\mathop{\arg\min}}_{x}\Big{\\{}\|\nabla_{y}x\|_{1}+\frac{1}{2\widetilde{\sigma}_{k}}\|x-s^{k}\|^{2}+\frac{1}{2\sigma}\|x-s\|^{2}\Big{\\}}$
(13) $\displaystyle=$
$\displaystyle\operatorname*{\mathop{\arg\min}}_{x}\Big{\\{}\|\nabla_{y}x\|_{1}+\frac{1}{2\widehat{\sigma}}\|x-\widehat{s}\|^{2}\Big{\\}}$
$\displaystyle=$ $\displaystyle{\rm
Prox}_{\widehat{\sigma}p_{1}}(\widehat{s}),$

where
$\widehat{\sigma}=\sigma\widetilde{\sigma}_{k}/(\sigma+\widetilde{\sigma}_{k})<\sigma$,
$\widehat{s}=\widehat{\sigma}(s/\sigma+s^{k}/\widetilde{\sigma}_{k})$.

As we know, for $s_{v}\in{\mathbb{R}}^{n},\,\lambda>0$, let
$\|\nabla\cdot\|_{1}$ denote the first order difference operator of a vector,
then the proximal mapping of $\|\nabla\cdot\|_{1}$

$\displaystyle{\rm
Prox_{\lambda\,\|\nabla\cdot\|_{1}}}(s_{v})=\arg\min_{x\in\mathbb{R}^{n}}\Big{\\{}\|\nabla
x_{v}\|_{1}+\frac{1}{2\lambda}\|x_{v}-s_{v}\|^{2}\Big{\\}}$ (14)

can be obtained by the fast algorithm proposed by Condat Condat . Hence, ${\rm
Prox}_{\widehat{\sigma}p_{1}}(\widehat{s})$ in (13) can be computed column by
column via (14), i.e.,

$\displaystyle{\rm Prox}_{\sigma\widehat{p}_{3}}(s)$ $\displaystyle=$
$\displaystyle\operatorname*{\mathop{\arg\min}}_{v}\Big{\\{}\widehat{p}_{3}(v)+\frac{1}{2\sigma}\|v-s\|^{2}\Big{\\}}$
(15) $\displaystyle=$
$\displaystyle\operatorname*{\mathop{\arg\min}}_{v}\Big{\\{}\lambda_{1}\|h(v)\|_{1}-\langle
g^{k}_{h},h(v)\rangle+\frac{1}{2\sigma}\|v-s\|^{2}\Big{\\}},$

where $h(v)$ is a vector defined in Section 1. Due to the separability of
$\|\cdot\|_{1}$ and $\|\cdot\|^{2}$, we can solve the problem (15) in the
following separate form

$\displaystyle
v(:,i)=\arg\min_{v(:,i)}\Big{\\{}(\lambda_{1}-(g_{h}^{k})_{i})\|v(:,i)\|+\frac{1}{2\sigma}\|v(:,i)-s(:,i)\|^{2}\Big{\\}},\;i=1,\cdots,n,$
(16)

where $(g_{h}^{k})_{i}$ is the $i$th element of the vector $g_{h}^{k}$. Since
we take $\alpha=3.7$ in the SCAD function, we have
$\lambda_{1}-(g_{h}^{k})_{i}>0$, which guarantees that the problem (16) is
convex.

Now, we describe the dADMM for the problem (12) as follows.

Algorithm 2 The dADMM for the problem (12):

Let $\tau\in(0,(1+\sqrt{5})/2)$ be a scalar parameter, $\sigma>0$ be a given
parameter, and choose
$\\{x^{0},y^{0},z^{0},\widehat{x}^{0},\widehat{y}^{0},s^{0},u^{0},v^{0}\\}$.
For $j=0,1,2,\ldots$, iterate:

1.

Solve the following linear system to compute $(x^{j+1},y^{j+1})$:

$\displaystyle\begin{bmatrix}2I&-I\\\
-I&2I\end{bmatrix}\begin{bmatrix}x^{j+1}\\\ y^{j+1}\end{bmatrix}$
$\displaystyle=$
$\displaystyle\frac{1}{\sigma}\begin{bmatrix}s^{j}+u^{j}-f-\sigma(z^{j}+\widehat{x}^{j}-g_{1}^{k}-g_{2}^{k})\\\
v^{j}-s^{j}+\sigma(z^{j}-\widehat{y}^{j}-g_{1}^{k})\end{bmatrix}.$

2.

Compute $(\widehat{x}^{j+1},\widehat{y}^{j+1},z^{j+1})$:

$\displaystyle\begin{bmatrix}\widehat{x}^{j+1}\\\ \widehat{y}^{j+1}\\\
z^{j+1}\end{bmatrix}$ $\displaystyle=$ $\displaystyle\begin{bmatrix}{\rm
Prox}_{\frac{1}{\sigma}p_{2}^{*}}(\frac{1}{\sigma}u^{j}-x^{j+1}+g_{2}^{k})\\\
{\rm
Prox}_{\frac{1}{\sigma}\widehat{p}_{3}^{*}}(\frac{1}{\sigma}v^{j}-y^{j+1})\\\
{\rm
Prox}_{\frac{1}{\sigma}\widehat{p}_{1}^{*}}(\frac{1}{\sigma}s^{j}-x^{j+1}+y^{j+1}+g_{1}^{k})\end{bmatrix}.$

3.

Update $(s^{j+1},u^{j+1},v^{j+1})$:

$\displaystyle s^{j+1}$ $\displaystyle=$ $\displaystyle
s^{j}+\tau\sigma(-x^{j+1}+y^{j+1}+g_{1}^{k}-z^{j+1}),$ $\displaystyle u^{j+1}$
$\displaystyle=$ $\displaystyle
u^{j}+\tau\sigma(-x^{j+1}+g_{2}^{k}-\widehat{x}^{j+1}),$ $\displaystyle
v^{k+1}$ $\displaystyle=$ $\displaystyle
v^{j}+\tau\sigma(-y^{j+1}-\widehat{y}^{j+1}).$

4.

If a stopping criterion is not satisfied, go to Step 1.

##  4 Convergence Analysis

In this section, we analyze the convergence of the algorithm.

###  4.1 Convergence analysis for the inexact dADMM

Since the inexact dADMM for the inner subproblem can be regarded as a special
case of the inexact symmetric Gauss-Seidel method ChenST , the convergence
analysis can also be borrowed to here.

###### Theorem 4.1

Suppose that the solution set to the KKT system (11) is nonempty and the
sequence
$\\{(s^{j},u^{j},v^{j},x^{j},y^{j},z^{j},\widehat{x}^{j},\widehat{y}^{j})\\}$
is generated by the dADMM, then $\\{(s^{j},u^{j},v^{j})\\}$ converges to the
solution of the primal problem (10) and
$\\{(x^{j},y^{j},z^{j},\widehat{x}^{j},\widehat{y}^{j})\\}$ converges to the
solution of the dual problem (12).

###  4.2 Convergence analysis for the PMM algorithm

In this subsection, we focus on the convergence analysis of the PMM algorithm.
Since the inner subproblem is solved inexactly, the error must be considered
if analyzing the convergence of the PMM algorithm. When solving the subproblem
(12), we actually solve the following problem

$\displaystyle\begin{aligned}
\min_{x,y,z,\widehat{x},\widehat{y}}&\Big{\\{}\widehat{p}_{1}^{*}(z+r_{z}^{k+1})+p_{2}^{*}(\widehat{x}+r_{\widehat{x}}^{k+1})+\widehat{p}_{3}^{*}(\widehat{y}+r_{\widehat{y}}^{k+1})+\langle
x,f\rangle+\\\ &\langle x,r_{x}^{k+1}\rangle+\langle
y,r_{y}^{k+1}\rangle+\langle
z,r_{z}^{k+1}\rangle+\langle\widehat{x},r_{\widehat{x}}^{k+1}\rangle+\langle\widehat{y},r_{\widehat{y}}^{k+1}\rangle\Big{\\}}\\\
\mbox{s.t.}\quad&-x+y+g_{1}^{k}-z-r_{1}^{k+1}=0,\\\
&-x+g_{2}^{k}-\widehat{x}-r_{2}^{k+1}=0,\\\
&-y-\widehat{y}-r_{3}^{k+1}=0,\end{aligned}$ (17)

where
$r_{1}^{k+1},r_{2}^{k+1},r_{3}^{k+1},r_{x}^{k+1},r_{y}^{k+1},r_{z}^{k+1},r_{\widehat{x}}^{k+1},r_{\widehat{y}}^{k+1}\in\mathbb{R}^{m\times
n}$ are error matrices. Therefore, the original primal problem we actually
solve is

$\displaystyle\begin{aligned}
\min_{s}&\Big{\\{}\widehat{p}_{1}(s-r_{z}^{k+1})+p_{2}(f-s+r_{x}-r_{\widehat{x}}^{k+1})+\widehat{p}_{3}(s+r_{y}^{k+1}-r_{\widehat{y}}^{k+1})-\\\
&\langle s,g^{k}_{1}-r_{1}^{k+1}\rangle-\langle
f-s+r_{x}^{k+1},g^{k}_{2}-r_{2}^{k+1}\rangle-\langle
s+r_{y}^{k+1},-r_{3}^{k+1}\rangle-\\\ &\langle
r_{z}^{k+1},\,s-r_{z}^{k+1}\rangle-\langle
r_{\widehat{x}}^{k+1},\,f-s+r_{x}^{k+1}-r_{\widehat{x}}^{k+1}\rangle-\langle
r_{\widehat{y}}^{k+1},\,s+r_{y}^{k+1}-r_{\widehat{y}}^{k+1}\rangle\Big{\\}}.\\\
\end{aligned}$ (18)

Let
$r^{k+1}:=(r_{1}^{k+1},r_{2}^{k+1},r_{3}^{k+1},r_{x}^{k+1},r_{y}^{k+1},r_{z}^{k+1},r_{\widehat{x}}^{k+1},r_{\widehat{y}}^{k+1})\in\mathbb{R}^{m\times
n}$. Note that $r^{k+1}$ converges to zero, it follows that if the inner
subproblem is adequately iterated, then $r^{k+1}$ satisfies

$\displaystyle
2p_{1}(r_{z}^{k+1})+2p_{2}(r_{x}^{k+1}-r_{\widehat{x}}^{k+1})+2p_{3}(r_{y}^{k+1}-r_{\widehat{y}}^{k+1})$
$\displaystyle+\frac{1}{2\widetilde{\sigma}_{k}}\|r_{z}^{k+1}\|^{2}+\textrm{pert}(r^{k+1})\leq\frac{1}{4\widetilde{\sigma}_{k}}\|s^{k+1}-s^{k}\|^{2},$
(19)

where

$\displaystyle\textrm{pert}(r^{k+1}):=$ $\displaystyle|\langle
s,\,r_{1}^{k+1}\rangle+\langle f-s+r_{x}^{k+1},\,r_{2}^{k+1}\rangle+\langle
s+r_{y}^{k+1},\,r_{3}^{k+1}\rangle+\langle r_{z}^{k+1},\,s-r_{z}^{k+1}\rangle$
$\displaystyle+\langle
r_{\widehat{x}}^{k+1},\,f-s+r_{x}^{k+1}-r_{\widehat{x}}^{k+1}\rangle+\langle
r_{\widehat{y}}^{k+1},\,s+r_{y}^{k+1}-r_{\widehat{y}}^{k+1}\rangle|.$

In the following, for the convenience of statement, we denote

$\displaystyle f_{k}(s)$ $\displaystyle:=$
$\displaystyle\widehat{p}_{1}(s-r_{z}^{k+1})+p_{2}(f-s+r_{x}^{k+1}-r_{\widehat{x}}^{k+1})+\widehat{p}_{3}(s+r_{y}^{k+1}-r_{\widehat{y}}^{k+1})-q_{1}(\nabla_{y}s^{k})$
$\displaystyle-\langle g_{1}^{k},\,s-s^{k}\rangle-
q_{2}(\nabla_{x}f-\nabla_{x}s^{k})-\langle
g_{2}^{k},\,(f-s)-(f-s^{k})\rangle+f^{P}_{k}(s).$

where $f^{P}_{k}(s):=\langle s,\,r_{1}^{k+1}\rangle+\langle
f-s+r_{x}^{k+1},\,r_{2}^{k+1}\rangle+\langle
s+r_{y}^{k+1},\,r_{3}^{k+1}\rangle+\langle
r_{z}^{k+1},\,s-r_{z}^{k+1}\rangle+\langle
r_{\widehat{x}}^{k+1},\,f-s+r_{x}^{k+1}-r_{\widehat{x}}^{k+1}\rangle+\langle
r_{\widehat{y}}^{k+1},\,s+r_{y}^{k+1}-r_{\widehat{y}}^{k+1}\rangle$. Then at
the $k$th iteration of the PMM algorithm, the actual problem we solve is

$\displaystyle
s^{k+1}=\operatorname*{\mathop{\arg\min}}_{s}\Big{\\{}f_{k}(s)\Big{\\}}.$

For the proof of the convergence theory, we must first prove the sufficiently
decent property of the sequence $\\{g(s^{k})\\}$.

###### Lemma 1

Let $\\{s^{k}\\}$ be a sequence generated by the PMM algorithm. We have the
following descent property

$\displaystyle g(s^{k})\geq
g(s^{k+1})+\frac{1}{4\widetilde{\sigma}}_{k}\|s^{k+1}-s^{k}\|^{2}.$

###### Proof

Firstly, one can easily prove that the function $p_{i}\,(i=1,2\textrm{ or }3)$
satisfies

$\displaystyle|p_{i}(s_{1})-p_{i}(s_{2})|\leq p_{i}(s_{1}-s_{2}),\
\forall\,s_{1},s_{2}\in\mathbb{R}^{m\times n}.$ (20)

Then, by (20) we have

$\displaystyle f_{k}(s^{k})$ $\displaystyle=$ $\displaystyle
g(s^{k})+p_{1}(s^{k}-r_{z}^{k+1})-p_{1}(s^{k})+p_{2}(f-s^{k}+r_{x}^{k+1}-r_{\widehat{x}}^{k+1})-p_{2}(f-s^{k})$
(21)
$\displaystyle+p_{3}(s^{k}+r_{y}^{k+1}-r_{\widehat{y}}^{k+1})-p_{3}(s^{k})+\frac{1}{2\widetilde{\sigma}}_{k}\|r_{z}^{k+1}\|^{2}+f^{P}_{k}(s^{k}),$
$\displaystyle\leq$ $\displaystyle
g(s^{k})+p_{1}(r_{z}^{k+1})+p_{2}(r_{x}^{k+1}-r_{\widehat{x}}^{k+1})+p_{3}(r_{y}^{k+1}-r_{\widehat{y}}^{k+1})+\frac{1}{2\widetilde{\sigma}}_{k}\|r_{z}^{k+1}\|^{2}$
$\displaystyle+\textrm{pert}(r^{k+1}).$

In addition, since $q_{i}\,(i=1,2,3)$ is a convex function and $s^{k+1}$ is a
minimizer of the function $f_{k}$, we have

$\displaystyle f_{k}(s^{k})$ $\displaystyle\geq$ $\displaystyle
f_{k}(s^{k+1})$ $\displaystyle\geq$ $\displaystyle
p_{1}(s^{k}-r_{z}^{k+1})+\frac{1}{2\widetilde{\sigma}_{k}}\|s^{k+1}-s^{k}\|^{2}-q_{1}(\nabla_{y}s^{k+1})+$
$\displaystyle
p_{2}(f-s^{k+1}+r_{x}^{k+1}-r_{\widehat{x}}^{k+1})-q_{2}(\nabla_{x}f-\nabla_{x}s^{k+1})+p_{3}(s^{k+1}+r_{y}^{k+1}-r_{\widehat{y}}^{k+1})$
$\displaystyle-q_{3}(h(s^{k+1})),$ $\displaystyle\geq$ $\displaystyle
g(s^{k+1})-p_{1}(r_{z}^{k+1})-p_{2}(r_{x}^{k+1}-r_{\widehat{x}}^{k+1})-p_{3}(r_{y}^{k+1}-r_{\widehat{y}}^{k+1})+\frac{1}{2\widetilde{\sigma}_{k}}\|s^{k+1}-s^{k}\|^{2}.$

The last inequality in (LABEL:ineq:fk-lower-estimate) is according to (20). In
combination of (21) and (LABEL:ineq:fk-lower-estimate), we obtain

$\displaystyle g(s^{k})$ $\displaystyle\geq$ $\displaystyle
g(s^{k+1})-2p_{1}(r_{z}^{k+1})-2p_{2}(r_{x}^{k+1}-r_{\widehat{x}}^{k+1})-2p_{3}(r_{y}^{k+1}-r_{\widehat{y}}^{k+1})-\frac{1}{2\widetilde{\sigma}}_{k}\|r_{z}^{k+1}\|^{2}$
(23)
$\displaystyle-\textrm{pert}(r^{k+1})+\frac{1}{2\widetilde{\sigma}_{k}}\|s^{k+1}-s^{k}\|^{2}.$

According to the condition (19), we get the result of this lemma.

In the following, we present the result about the criterion of being a
d-stationary point. But before that, we give a definition of a function. For
$\widetilde{s}\in\mathbb{R}^{m\times n}$, $\widetilde{\sigma}>0$,

$\displaystyle\widehat{f}(s;\widetilde{s},\widetilde{\sigma})$
$\displaystyle:=$ $\displaystyle
p_{1}(s)+p_{2}(f-s)+p_{3}(s)-q_{1}(\nabla_{y}\widetilde{s})-\langle\nabla_{y}^{*}\nabla
q_{\alpha,\lambda_{1}}(\nabla_{y}\widetilde{s}),\,s-\widetilde{s}\rangle-$
$\displaystyle
q_{\alpha,\lambda_{2}}(\nabla_{x}(f-\widetilde{s}))-\langle\nabla_{x}^{*}\nabla
q_{\alpha,\lambda_{2}}(\nabla_{x}f-\nabla_{x}\widetilde{s}),\,s-\widetilde{s}\rangle$
$\displaystyle-q_{\alpha,\lambda_{3}}(h(\widetilde{s}))-\langle\nabla
q_{\alpha,\lambda_{3}}(h(\widetilde{s})),\,h(s)-h(\widetilde{s})\rangle+\frac{1}{2\widetilde{\sigma}}\|s-\widetilde{s}\|^{2}.$

Then in order to prove the result, we need a lemma. Since it is very similar
to (CuiPS, , Lemma 5, Lemma6), we present it without proof.

###### Lemma 2

Let $\bar{s}\in\mathbb{R}^{m\times n}$, then $\bar{s}$ is a d-stationary point
of (5) if and only if there exist $\widetilde{\sigma}\geq 0$ such that
$\bar{s}$ is a solution of the following minimization problem

$\displaystyle\min_{s\in\mathbb{R}^{m\times
n}}\Big{\\{}\widehat{f}(s;\widetilde{s},\widetilde{\sigma})\Big{\\}}.$

Now we present the convergence result of the PMM algorithm.

###### Theorem 4.2

Assume that the objective function in (5) is bounded below and
$\\{\widetilde{\sigma}_{k}\\}$ are positive convergent sequences. Then every
cluster point of the sequence $\\{s^{k}\\}$ generated by the PMM algorithm is
a d-stationary point of (5).

###### Proof

It is known from Lemma 1 that $\\{g(s^{k})\\}$ is a non-increasing sequence.
By the assumption that the value of $g(s)$ is bounded below, we obtain that
$\\{g(s^{k})\\}$ converges. Meanwhile,
$\lim_{k\rightarrow\infty}\|s^{k+1}-s^{k}\|=0$. Let
$\\{s^{k}\\}_{k\in\mathcal{K}}$ be a subsequence of $\\{s^{k}\\}$, and
$\lim_{k(\in\mathcal{K})\rightarrow\infty}s^{k}=s^{\infty}$. In the following,
we prove that $s^{\infty}$ is a d-stationary point of (5).

Since $s^{k+1}$ is a minimizer of the function $f_{k}$, we have

$\displaystyle\widehat{f}(s;s^{k},\widetilde{\sigma}_{k})+p_{1}(r_{z}^{k+1})+p_{2}(r_{x}^{k+1}-r_{\widehat{x}}^{k+1})+p_{3}(r_{y}^{k+1}-r_{\widehat{y}}^{k+1})+\frac{1}{2\widetilde{\sigma}_{k}}\|r_{z}^{k+1}\|^{2}$
$\displaystyle-\langle g_{1}^{k},\,r_{z}^{k+1}\rangle+\langle
g_{2}^{k},\,r_{x}^{k+1}-r_{\widehat{x}}^{k+1}\rangle-\langle|\nabla
q_{3}(h(s^{k}))|,\,h(r_{y}^{k+1}-r_{\widehat{y}}^{k+1})\rangle$
$\displaystyle\geq$ $\displaystyle f_{k}(s)\geq f_{k}(s^{k+1})$
$\displaystyle\geq$
$\displaystyle\widehat{f}(s^{k+1};s^{k},\widetilde{\sigma}_{k})-p_{1}(r_{z}^{k+1})-p_{2}(r_{x}^{k+1}-r_{\widehat{x}}^{k+1})-p_{3}(r_{y}^{k+1}-r_{\widehat{y}}^{k+1})-\frac{1}{2\widetilde{\sigma}_{k}}\|r_{z}^{k+1}\|^{2}$
$\displaystyle+\langle g_{1}^{k},\,r_{z}^{k+1}\rangle-\langle
g_{2}^{k},\,r_{x}^{k+1}-r_{\widehat{x}}^{k+1}\rangle+\langle|\nabla
q_{3}(h(s^{k}))|,\,h(r_{y}^{k+1}-r_{\widehat{y}}^{k+1})\rangle,\
\forall\,s\in\mathbb{R}^{m\times n}.$

We assume
$\lim_{k\rightarrow\infty}\widetilde{\sigma}_{k}=\widetilde{\sigma}_{\infty}$.
Note that $\|r^{k+1}\|\rightarrow 0$ due to (19) and $h$ is a continuous
function. Then letting $k(\in\mathcal{K})\rightarrow\infty$, we derive

$\displaystyle\widehat{f}(s;s^{\infty},\widetilde{\sigma}_{\infty})\geq\widehat{f}(s^{\infty};s^{\infty},\widetilde{\sigma}_{\infty}),\
\forall\,s^{\infty}\in\mathbb{R}^{m\times n}.$

Hence, by Lemma 2 we can get the desired result of this theorem.

Finally, we establish the local convergence rate of the sequence $\\{s^{k}\\}$
in the following. But before that, we present a lemma.

###### Lemma 3

In $\mathbb{R}^{m}$, define a function
$w:\,\mathbb{R}^{m}\rightarrow\mathbb{R}$, $w(x_{v}):=\|x_{v}\|$, then there
is a constant $\kappa\geq 0$ such that

$\displaystyle\partial w(x_{v})\subseteq\partial
w(\bar{x}_{v})+\kappa\|x_{v}-\bar{x}_{v}\|\mathbb{B}_{m},\ \
\forall\,x_{v}\in\mathbb{R}^{m}.$ (24)

###### Proof

We discuss this problem in three cases:
(i). If $x_{v}\neq 0$ and $\bar{x}_{v}\neq 0$, then $\partial
w(x_{v})=\Big{\\{}\frac{x_{v}}{\|x_{v}\|}\Big{\\}}$, $\partial
w(\bar{x}_{v})=\Big{\\{}\frac{\bar{x}_{v}}{\|\bar{x}_{v}\|}\Big{\\}}$, we can
see that the Lipschitz continuity holds, which implies that the result (24)
holds.
(ii). If $x_{v}\neq 0$ and $\bar{x}_{v}=0$, then $\partial
w(x_{v})=\Big{\\{}\frac{x_{v}}{\|x_{v}\|}\Big{\\}}$, $\partial
w(\bar{x}_{v})=\mathbb{B}_{m}$, the result (24) obviously holds.
(iii). If $x_{v}=0$ and $\bar{x}_{v}\neq 0$, then $\partial
w(x_{v})=\mathbb{B}_{m}$, $\partial
w(\bar{x}_{v})=\Big{\\{}\frac{\bar{x}_{v}}{\|\bar{x}_{v}\|}\Big{\\}}$,
therefore (24) holds with $\kappa\geq\frac{1}{\|x_{v}-\bar{x}_{v}\|}$.

###### Theorem 4.3

Assume that the objective function in (5) is bounded below. Let $\\{s^{k}\\}$
be the sequence generated by the PMM algorithm with $\mathcal{C}^{\infty}$ as
the set of all its cluster points. The whole sequence $\\{s^{k}\\}$ converges
to an element of $\mathcal{C}^{\infty}$, if one of the following conditions
holds:
(I) The set $\mathcal{C}^{\infty}$ contains an isolated element;
(II) $\\{s^{k}\\}$ is a bounded sequence.
Furthermore, based on the condition (II), let
$\lim\limits_{k\rightarrow\infty}s^{k}=s^{\infty}\in\mathcal{C}^{\infty}$ and
the function $g$ satisfies the KŁ property at $s^{\infty}$ with an exponent
$\gamma\in[0,1)$, we have
(a) The sequence $\\{s^{k}\\}$ converges in a finite number of steps, if
$\gamma=0$;
(b) The sequence $\\{s^{k}\\}$ converges R-linearly, i.e., there exist $\mu>0$
and $\theta\in[0,1)$ such that $\|s^{k}-s^{\infty}\|\leq\mu\theta^{{}^{k}}$,
if $0<\gamma\leq\frac{1}{2}$ and $k$ is sufficiently large;
(c) The sequence $\\{s^{k}\\}$ converges R-sublinearly, i.e., there exists
$\mu>0$ such that $\|s^{k}-s^{\infty}\|\leq\mu
k^{-\frac{1-\gamma}{2\gamma-1}}$, if $\frac{1}{2}<\gamma<1$ and $k$ is
sufficiently large.

###### Proof

Since the objective function in (5) is bounded below, it follows from the
proof of Theorem 4.2 that
$\lim\limits_{k\rightarrow\infty}\|s^{k+1}-s^{k}\|=0$. Under condition (I),
the sequence converges to the isolated element of $\mathcal{C}^{\infty}$ due
to (Facchinei2003, , Proposition 8.3.10). To prove the sequential convergence
of $\\{s^{k}\\}$ under the condition (II), it suffices to prove the following
three properties of $\\{s^{k}\\}$ and then apply the convergence results in
(AttouchBS, , Theorem 2.9) to the function $g$:
(i) The sequence $\\{g(s^{k})\\}$ is descent, i.e., $g(s^{k+1})\geq
g(s^{k})+\frac{c}{2}\|s^{k+1}-s^{k}\|^{2},\,\forall\,k\geq 1$, where $c>0$ is
a constant;
(ii) There exists a subsequence $\\{s^{k}\\}_{k\in\mathcal{K}}$ of
$\\{s^{k}\\}$ such that
$\lim\limits_{k(\in\mathcal{K})\rightarrow\infty}s^{k}=s^{\infty}$ and
$\lim\limits_{k(\in\mathcal{K})\rightarrow\infty}g(s^{k})=g(s^{\infty})$;
(iii) There exist a constant $K>0$ and $\varepsilon^{k+1}\in\partial
g(s^{k+1})$ such that $\|\varepsilon^{k+1}\|\leq K\|s^{k+1}-s^{k}\|$, for $k$
sufficiently large.
Obviously, the first two properties have been proven. Now we focus on the
proof of the third property. Let

$\displaystyle\bar{\varepsilon}^{k+1}:=g_{1}^{k}-g_{1}^{k+1}-g_{2}^{k}+g_{2}^{k+1}-\frac{1}{\widetilde{\sigma}_{k}}(s^{k+1}-s^{k}-r_{z}^{k+1})-r_{1}^{k+1}-r_{z}^{k+1}.$

Then $\bar{\varepsilon}^{k+1}\in\partial g_{k}(s^{k+1})$, where

$\displaystyle g_{k}(s)$ $\displaystyle=$ $\displaystyle
p_{1}(s-r_{z}^{k+1})+p_{2}(f-s+r_{x}^{k+1}-r_{\hat{x}}^{k+1})+p_{3}(s+r_{y}^{k+1}-r_{\widehat{y}}^{k+1})$
$\displaystyle-
q_{\alpha,\lambda_{1}}(\nabla_{y}s)-q_{\alpha,\lambda_{2}}(\nabla_{x}f-\nabla_{x}s)-q_{\alpha,\lambda_{3}}(h(s+r_{y}^{k+1}-r_{\widehat{y}}^{k+1})).$

According to the special structure of the function $h$, we may consider every
element $h_{i}$ separately. Since the functions $\partial p_{1}$ and $\partial
p_{2}$ are piecewise polyhedral, $\|r_{z}^{k+1}\|$ and
$\|r_{x}^{k+1}-r_{\widehat{x}}^{k+1}\|$ are sufficiently small for all $k$
sufficiently large, we can deduce from Robinson1981  that $\partial p_{1}$ and
$\partial p_{2}$ are locally upper Lipschitz continuous. Hence, it is known
from Lemma 3 and (Rockafellar70, , Theorem 23.8) that there exist
$\kappa_{1},\kappa_{2},\kappa_{3}\geq 0$, such that

$\displaystyle\bar{\varepsilon}^{k+1}$ $\displaystyle\in$
$\displaystyle\partial p_{1}(s^{k+1}-r_{z}^{k+1})+\partial
p_{2}(f-s^{k+1}+r_{x}^{k+1}-r_{\widehat{x}}^{k+1})+$ $\displaystyle\partial
p_{3}(s^{k+1}+r_{y}^{k+1}-r_{\widehat{y}}^{k+1})-\partial
q_{\alpha,\lambda_{3}}(h(s^{k+1}+r_{y}^{k+1}-r_{\widehat{y}}^{k+1}))-g_{1}^{k+1}+g_{2}^{k+1}$
$\displaystyle\subset$ $\displaystyle\partial
p_{1}(s^{k+1})+\kappa_{1}\|r_{z}^{k+1}\|\mathbb{B}_{m\times n}+\partial
p_{2}(f-s^{k+1})+\kappa_{2}\|r_{x}^{k+1}-r_{\hat{x}}^{k+1}\|\mathbb{B}_{m\times
n}$ $\displaystyle+\partial p_{3}(s^{k+1})-\partial
q_{\alpha,\lambda_{3}}(h(s^{k+1}))+\kappa_{3}\|r_{y}^{k+1}-r_{\widehat{y}}^{k+1}\|\mathbb{B}_{m\times
n}-g_{1}^{k+1}+g_{2}^{k+1}$ $\displaystyle=$ $\displaystyle\partial
p_{1}(s^{k+1})+\partial p_{2}(f-s^{k+1})+\partial p_{3}(s^{k+1})-\partial
q_{\alpha,\lambda_{3}}(h(s^{k+1}))+(\kappa_{1}\|r_{z}^{k+1}\|$
$\displaystyle+\kappa_{2}\|r_{x}^{k+1}-r_{\widehat{x}}^{k+1}\|+\kappa_{3}\|r_{y}^{k+1}-r_{\widehat{y}}^{k+1}\|)\mathbb{B}_{m\times
n}-g_{1}^{k+1}+g_{2}^{k+1}.$

Hence, we can find $d\bar{\varepsilon}^{k+1}\in\mathbb{R}^{m\times n}$ with
$\|d\bar{\varepsilon}^{k+1}\|\leq\kappa_{1}\|r_{z}^{k+1}\|+\kappa_{2}\|r_{x}^{k+1}-r_{\widehat{x}}^{k+1}\|+\kappa_{3}\|r_{y}^{k+1}-r_{\widehat{y}}^{k+1}\|$,
such that
$\varepsilon^{k+1}:=\bar{\varepsilon}^{k+1}+d\bar{\varepsilon}^{k+1}\in\partial
g(s^{k+1})$. From (19), we know that there exists a constant $\kappa_{4}>0$,
such that $\|d\bar{\varepsilon}^{k+1}\|\leq\kappa_{4}\|s^{k+1}-s^{k}\|$. In
combination of (19) and the Lipschitz continuity of $\nabla_{y}^{*}\nabla
q_{\alpha,\lambda_{1}}$ and $\nabla_{x}^{*}\nabla q_{\alpha,\lambda_{2}}$, we
can find some $\kappa_{5}>0$, such that
$\|\bar{\varepsilon}^{k+1}\|\leq\kappa_{5}\|s^{k+1}-s^{k}\|$. Therefore, we
can get the result of (iii). Finally, since we have proved the properties
(i)-(iii), the results for the convergence rate can be established similar to
(AttouchB, , Theroem 2) or (bolte-pauwels2016, , Proposition 4).

##  5 Numerical Experiments

In this section, we implement some numerical experiments to demonstrate the
efficiency of our algorithm. All the experiments are conducted on a Desktop
with Intel(R) Core(TM) i5-8250U<EMAIL_ADDRESS>of 8G memory running 64 bit
Windows operation system. All the codes are written in Matlab (R2018a) with
some subroutines written in C.

In order to give an overall evaluation, we select the full-reference
evaluation indices: the peak signal-to-noise ratio (PSNR) and the structural
similarity (SSIM) wang2004image , which are defined as follows.

######  Definition 2 (Peak Signal-to-noise Ratio (PSNR))

Given a clean image $u$ and a degraded image $v$ of $m\times n$, the PSNR is
defined as

$\textrm{PSNR}:=10\log_{10}(\frac{u_{\max}^{2}}{\textrm{MSE}}),$

where MSE denotes the mean square error and is defined as

$\textrm{MSE}:=\frac{1}{mn}\sum_{i=1}^{m}\sum_{j=1}^{n}(u_{ij}-v_{ij})^{2}.$

######  Definition 3 (Structual Similarity (SSIM) wang2004image )

The mean values of the images $u$ and $v$ are denoted as

$\displaystyle\mu_{u}$
$\displaystyle=\frac{1}{mn}\sum_{i=1}^{m}\sum_{j=1}^{n}u_{ij}$

and

$\displaystyle\mu_{v}$
$\displaystyle=\frac{1}{mn}\sum_{i=1}^{m}\sum_{j=1}^{n}v_{ij},$

respectively. The standard deviations are denoted as

$\displaystyle\sigma_{u}$
$\displaystyle=\Big{(}\frac{1}{mn-1}\sum_{i=1}^{m}\sum_{j=1}^{n}(u_{ij}-\mu_{u})^{2}\Big{)}^{\frac{1}{2}}$

and

$\displaystyle\sigma_{v}$
$\displaystyle=\Big{(}\frac{1}{mn-1}\sum_{i=1}^{m}\sum_{j=1}^{n}(v_{ij}-\mu_{v})^{2}\Big{)}^{\frac{1}{2}},$

respectively. The covariance of the images $u$ and $v$ is denoted as

$\displaystyle\sigma_{uv}$
$\displaystyle=\frac{1}{mn-1}\sum_{i=1}^{m}\sum_{j=1}^{n}(u_{ij}-\mu_{u})(v_{ij}-\mu_{v}).$

And the factors of the luminance, contrast, and structure of the images $u$
and $v$ are defined as

$l(u,v):=\frac{2\mu_{u}\mu_{v}+c_{1}}{\mu_{u}^{2}+\mu_{v}^{2}+c_{1}},\;c(u,v):=\frac{2\sigma_{u}\sigma_{v}+c_{2}}{\sigma_{u}^{2}+\sigma_{v}^{2}+c_{2}},\;s(u,v):=\frac{\sigma_{uv}+c_{3}}{\sigma_{u}\sigma_{v}+c_{3}},$

where $c_{1}=(k_{1}L_{u})^{2},c_{2}=(k_{2}L_{u})^{2},c_{3}=c_{2}/2$ are three
constants. $k_{1}=0.01$ and $k_{2}=0.03$ are the common default values. Then
we can define

$\textrm{SSIM}:=[l(u,v)]^{\alpha}[c(u,v)]^{\beta}[s(u,v)]^{\gamma}.$

In order to simplify the expression, we set $\alpha,\beta,\gamma=1$ in this
paper, then we have

$\textrm{SSIM}=\frac{(2\mu_{u}\mu_{v}+c_{1})(2\sigma_{uv}+c_{2})}{(\mu_{u}^{2}+\mu_{v}^{2}+c_{1})(\sigma_{u}^{2}+\sigma_{v}^{2}+c_{2})}.$

For the nonconvex model, the iteration process of the outer problem is
terminated if

$\max\\{R_{p}^{k},R_{d}^{k}\\}<\texttt{Tol}=2\times 10^{-4}$

is satisfied, where

$\displaystyle R_{p}^{k}=\frac{\|f-s^{k}-u^{k}\|+\|s^{k}-v^{k}\|}{1+\|f\|},$
$\displaystyle
R_{d}^{k}=\frac{\|-x^{k}-\widehat{x}^{k}\|+\|-x^{k}+y^{k}-z^{k}\|+\|-y^{k}-\widehat{y}^{k}\|}{1+\|f\|},$

or the iteration number reaches $5$. The iteration process of the subproblem
is terminated if

$\displaystyle
2p_{1}(r_{z}^{k+1})+2p_{2}(r_{x}^{k+1}-r_{\widehat{x}}^{k+1})+2p_{3}(r_{y}^{k+1}-r_{\widehat{y}}^{k+1})$
$\displaystyle+\frac{1}{2\widetilde{\sigma}_{k}}\|r_{z}^{k+1}\|^{2}+\textrm{pert}(r^{k+1})\leq\frac{1}{4\widetilde{\sigma}_{k}}\|s^{k+1}-s^{k}\|^{2},$

or the iteration number reaches $100$.

For the convex model, the iteration process is terminated if

$\max\\{R_{p}^{k},R_{d}^{k},R_{c}^{k}\\}<\texttt{Tol}=2\times 10^{-4}$

is satisfied, where

$\displaystyle R_{p}^{k}=\frac{\|f-s^{k}-u^{k}\|+\|s^{k}-v^{k}\|}{1+\|f\|},$
$\displaystyle
R_{d}^{k}=\frac{\|-x^{k}-\widehat{x}^{k}\|+\|-x^{k}+y^{k}-z^{k}\|+\|-y^{k}-\widehat{y}^{k}\|}{1+\|f\|},$
$\displaystyle R_{c}^{k}=\frac{\|s^{k}-{\rm
Prox}_{\widehat{p}_{1}}(s^{k}+z^{k})\|+\|u^{k}-{\rm
Prox}_{p_{2}}(u^{k}+\widehat{x}^{k})\|+\|v-{\rm
Prox}_{\widehat{p}_{3}}(v^{k}+\widehat{y}^{k})\|}{1+\|f\|},$

or the iteration number reaches $500$.

Three real remote sensing stripe images with different stripe noise
distributions are downloaded from the website
https://ladsweb.nascom.nasa.gov/. Although there are three parameters in the
model (3) and (4), we only need to do the cross validation for two parameters
by dividing one parameter.

In Table 1, we present the experimental results of different algorithms on the
convex model and the nonconvex model, respectively. The image is degraded by
nonperiodical and periodical stripes, respectively. From the table, we can see
that the nonconvex model performs better than the convex counterpart,
especially for the PSNR. In addition, for the convex model, the dADMM performs
better than the pADMM. This also demonstrates that we adopt the dADMM to solve
the subproblems of the PMM is reasonable.

Furthermore, we also present the stripe and destripe images in Figures 1-3 to
illustrate the effectiveness of the model. We can see that the proposed model
can remove the stripes, and preserve the stripe-free information and the image
details very well.

Table 1: Experimental results Stripe noise | Case | Index | Degraded | Convex model | Nonconvex model
---|---|---|---|---|---
pADMM | dADMM
Nonperiodical | Case 1 | KKT |  | 1.99e-4 | 1.99e-4 | 1.99e-4
PSNR | 23.05 | 55.47 | 59.13 | 63.36
SSIM | 0.6486 | 0.9995 | 0.9998 | 0.9999
time |  | 26.89 | 21.48 | 39.24
iter |  | 478 | 233 | 434
Case 2 | KKT |  | 1.99e-4 | 1.98e-4 | 1.99e-4
PSNR | 18.27 | 53.92 | 58.78 | 62.43
SSIM | 0.2753 | 0.9995 | 0.9997 | 0.9998
time |  | 24.38 | 17.56 | 33.21
iter |  | 401 | 218 | 382
Case 3 | KKT |  | 1.99e-4 | 1.99e-4 | 2.00e-4
PSNR | 24.33 | 43.67 | 48.82 | 55.38
SSIM | 0.3626 | 0.9943 | 0.9962 | 0.9976
time |  | 18.46 | 13.39 | 30.25
iter |  | 378 | 179 | 352
Periodical | Case 1 | KKT |  | 1.99e-4 | 2.00e-4 | 2.00e-4
PSNR | 20.68 | 52.00 | 57.13 | 62.00
SSIM | 0.5868 | 0.9987 | 0.9995 | 0.9999
time |  | 37.64 | 29.12 | 55.83
iter |  | 492 | 281 | 500
Case 2 | KKT |  | 1.99e-4 | 1.99e-4 | 1.79e-4
PSNR | 17.67 | 40.75 | 46.69 | 54.42
SSIM | 0.4577 | 0.9946 | 0.9960 | 0.9992
time |  | 16.57 | 11.28 | 28.26
iter |  | 455 | 180 | 440
Case 3 | KKT |  | 1.98e-4 | 1.99e-4 | 1.96e-4
PSNR | 18.32 | 41.05 | 48.36 | 55.98
SSIM | 0.3283 | 0.9949 | 0.9970 | 0.9995
time |  | 14.25 | 9.98 | 25.76
iter |  | 423 | 167 | 409

Figure 1: The destripe results. (a) Original image; (b) stripe component; (c)
degraded image; (d) estimated image by the convex model (pADMM); (e) estimated
image by the convex model (dADMM); (f) estimated image by the nonconvex model

Figure 2: The destripe results. (a) Original image; (b) stripe component; (c)
degraded image; (d) estimated image by the convex model (pADMM); (e) estimated
image by the convex model (dADMM); (f) estimated image by the nonconvex model

Figure 3: The destripe results. (a) Original image; (b) stripe component; (c)
degraded image; (d) estimated image by the convex model (pADMM); (e) estimated
image by the convex model (dADMM); (f) estimated image by the nonconvex model

##  6 Conclusion

In this paper, we have developed a nonconvex model for remote sensing stripe
images. For the computational issues, we have applied the inexact PMM
algorithm to solve the model with the inner subproblem solved by the dADMM.
Meanwhile, we have analyzed the convergence of the algorithm based on the KŁ
property. Finally, numerical experiments demonstrate the superiority of the
proposed model and algorithm.

## References

  * (1) Attouch, H., Bolte, J.: On the convergence of the proximal algorithm for nonsmooth functions involving analytic features.  Mathematical Programming 116, 5–16 (2009)
  * (2) Attouch, H., Bolte, J., Svaiter, B.: Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods.  Mathematical Programming 137, 91–129 (2013)
  * (3) Bolte, J., Pauwels, E.: Majorization-minimization procedures and convergence of sqp methods for semi-algebraic and tame programs.  Mathematics of Operations Research 41(2), 442–465 (2016)
  * (4) Bouali, M., Ladjal, S.: Toward optimal destriping of MODIS data using a unidirectional variational model.  IEEE Transactions on Geoscience and Remote Sensing 49(8), 2924–2935 (2011)
  * (5) Chen, L., Sun, D.F., Toh, K.C.: An efficient inexact symmetric gauss-seidel based majorized admm for high-dimensional convex composite conic programming.  Mathematical Programming 161, 237–270 (2017)
  * (6) Chen, Y., Huang, T.Z., Deng, L.J., Zhao, X.L., Wang, M.: Group sparsity based regularization model for remote sensing image stripe noise removal.  Neurocomputing 267, 95–106 (2017)
  * (7) Condat, L.: A direct algorithm for 1d total variation denoising.  IEEE Signal Processing Letters 20(11), 1054–1057 (2013)
  * (8) Cui, Y., Pang, J., Sen, B.: Composite difference-max programs for modern statistical estimation problems.  SIAM Journal on Optimization 28(4), 3344–3374 (2018)
  * (9) Du, B., Zhang, L., Tao, D., Zhang, D.: Unsupervised transfer learning for target detection from hyperspectral images.  Neurocomputing 120, 72–82 (2013)
  * (10) Facchinei, F., Pang, J.S.: Finite-Dimensional Variational Inequalities and Complementarity Problems.  Springer New York (2003)
  * (11) Fan, J., Li, R.: Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American Statistical Association 96(456), 1348–1360 (2017)
  * (12) Huang, H., Huang, Y.: Improved discriminant sparsity neighborhood preserving embedding for hyperspectral image classification.  Neurocomputing 136, 224–234 (2014)
  * (13) Huang, H., Tang, Y., Tan, Z., Zhuang, J., Hou, C., Chen, W., Ren, J.: Object-based attention mechanism for color calibration of UAV remote sensing images in precision agriculture.  IEEE Transactions on Geoscience and Remote Sensing 60, 1–13 (2022)
  * (14) Richards, J.A.: Remote Sensing Digital Image Analysis, 6th edn.  Springer (2022)
  * (15) Robinson, S.: Some continuity properties of polyhedral multifunctions.  Mathematical Programming Study 14, 206–214 (1981)
  * (16) Rockafellar, R.: Convex Analysis.  Princeton University Press (1970)
  * (17) Shen, H., Zhang, L.: A map-based algorithm for destriping and inpainting of remotely sensed images.  IEEE Transactions on Geoscience and Remote Sensing 47(5), 1492–1502 (2009)
  * (18) Tang, P., Wang, C., Jiang, B.: A proximal-proximal majorization-minimization algorithm for nonconvex rank regression problems.  to appear in Journal of Scientific Computing (2021)
  * (19) Tang, P., Wang, C., Sun, D., Toh, K.C.: A sparse semismooth newton based proximal majorization-minimization algorithm for nonconvex square-root-loss regression problems.  Journal of Machine Learning Research 21, 1–38 (2020)
  * (20) Wang, Z., Bovik, A.C., Sheikh, H.R., Simoncelli, E.P.: Image quality assessment: From error visibility to structural similarity.  IEEE Transactions on Image Processing 13(4), 600–612 (2004)
  * (21) Zhou, G., Chen, W., Gui, Q., Li, X., Wang, L.: Split depth-wise separable graph-convolution network for road extraction in complex environments from high-resolution remote-sensing images.  IEEE Transactions on Geoscience and Remote Sensing 60, 1–15 (2022)