######  Proposition 4.11 (Stability of Lipschitz representations).

If ${\boldsymbol{L}}_{n},{\boldsymbol{L}}:\mathcal{X}\to\mathcal{Y}$ are
$L$-Lipschitz invariant maps and $D\subset\mathcal{X}$ is an invariant closed
set, such that

${\boldsymbol{L}}_{n}X\to{\boldsymbol{L}}X\quad\text{for every }X\in D,$
(4.15)

then setting $\tilde{D}:=\big{\\{}X_{\sharp}\mathbb{P}:X\in D\big{\\}}$, we
have
${\boldsymbol{f}}_{{\boldsymbol{L}}_{n}}\to{\boldsymbol{f}}_{{\boldsymbol{L}}}$
pointwise in $\mathcal{S}(\mathsf{X},\tilde{D})$, where
${\boldsymbol{f}}_{{\boldsymbol{L}}}$ is as in Theorem 4.8.

###### Proof.

We fix $\mu\in\tilde{D}$ and we set
$\boldsymbol{f}_{n}:={\boldsymbol{f}}_{{\boldsymbol{L}}_{n}}(\cdot,\mu):\operatorname{supp}(\mu)\to\mathsf{Y}$.
We observe that $\boldsymbol{f}_{n}$ are $L$-Lipschitz and form a Cauchy
sequence in $L^{q}(\mathsf{X},\mu;\mathsf{Y})$: we denote by
$\tilde{\boldsymbol{f}}$ its limit. For every $x\in\operatorname{supp}(\mu)$
and $\rho>0$, we can set

$\boldsymbol{f}_{n,\rho}(x):=\frac{1}{\mu(B(x,\rho))}\int_{B(x,\rho)}\boldsymbol{f}_{n}(u)\,\mathrm{d}\mu(u),\quad\boldsymbol{f}_{\rho}(x):=\frac{1}{\mu(B(x,\rho))}\int_{B(x,\rho)}\tilde{\boldsymbol{f}}(u)\,\mathrm{d}\mu(u).$
(4.16)

We have $\boldsymbol{f}_{n,\rho}(x)\to\boldsymbol{f}_{\rho}(x)$ for every
$x\in\operatorname{supp}(\mu)$ and every $\rho>0$. On the other hand

$|\boldsymbol{f}_{n,\rho}(x)-\boldsymbol{f}_{n}(x)|_{\mathsf{Y}}\leq L\rho,$

so that a simple argument using the triangle inequality shows that the
sequence $(\boldsymbol{f}_{n}(x))_{n\in\mathbb{N}}$ is Cauchy in $\mathsf{Y}$
for every $x\in\operatorname{supp}(\mu)$ and its pointwise limit
$\boldsymbol{g}(\cdot,\mu)$ is $L$-Lipschitz and represents ${\boldsymbol{L}}$
as in (4.10) for every $X\in D$. Arguing as in claim (d) in the proof of
Theorem 4.8 and using the continuity of ${\boldsymbol{L}}$, we can eventually
deduce that $\boldsymbol{g}$ is continuous in
$\mathcal{S}(\mathsf{X},\tilde{D})$. The fact that
$\boldsymbol{g}={\boldsymbol{f}}_{{\boldsymbol{L}}}$ comes from Lemma 4.7. ∎

When the maps ${\boldsymbol{L}}_{n}$ are not uniformly Lipschitz, we can still
obtain a limiting representation.

######  Lemma 4.12.

Let $D\subset\mathcal{X}$ and
${\boldsymbol{L}}_{n},{\boldsymbol{L}}:D\to\mathcal{Y}$ be maps such that
${\boldsymbol{L}}_{n}$ are law invariant, ${\boldsymbol{L}}_{n}$ converge
pointwise to ${\boldsymbol{L}}$ on $D$ as $n\to\infty$ and there exist Borel
functions $\boldsymbol{f}_{n}:\mathcal{S}(\mathsf{X},\iota(D))\to\mathsf{Y}$
such that

for every $n\in\mathbb{N}$ and every $X\in D$,
${\boldsymbol{L}}_{n}X(\omega)=\boldsymbol{f}_{n}(X(\omega),X_{\sharp}\mathbb{P})$
for a.e. $\omega\in\Omega$.

Then ${\boldsymbol{L}}$ is law invariant and, for every $\mu\in\iota(D)$,
there exists a map $\boldsymbol{f}[\mu]\in L^{q}(\mathsf{X},\mu;\mathsf{Y})$
such that, for every $X\in D$ there exists an increasing subsequence $k\mapsto
n_{k}$ such that

$\lim_{k}\boldsymbol{f}_{n_{k}}(X(\omega),X_{\sharp}\mathbb{P})=\boldsymbol{f}[X_{\sharp}\mathbb{P}](X(\omega))={\boldsymbol{L}}X(\omega)\text{
for a.e.\leavevmode\nobreak\ }\omega\in\Omega.$

###### Proof.

First of all notice that ${\boldsymbol{L}}$ is law invariant, being pointwise
limit of law invariant maps. Let $\mu\in\iota(D)$ be fixed. If $X\in D$ is
such that $X_{\sharp}\mathbb{P}=\mu$, we can find an increasing subsequence
$k\mapsto n_{k}$ such that ${\boldsymbol{L}}_{n_{k}}X\to{\boldsymbol{L}}X$
$\mathbb{P}$-a.e. in $\Omega$ and define

$\boldsymbol{f}_{X}(x):=\begin{cases}\lim_{k\to+\infty}\boldsymbol{f}_{n_{k}}(x,\mu)\quad&\text{
if }x\in E,\\\ 0\quad&\text{ else},\end{cases}$

where $E\subset\mathsf{X}$ is the Borel set of points in $\mathsf{X}$ where
$\lim_{k\to+\infty}\boldsymbol{f}_{n_{k}}(x,\mu)$ exists. Let us now show that
$\boldsymbol{f}_{X}$ represents ${\boldsymbol{L}}$, i.e.

${\boldsymbol{L}}X(\omega)=\boldsymbol{f}_{X}(X(\omega))$ for a.e.
$\omega\in\Omega$.

We know that

for every $k\in\mathbb{N}$
${\boldsymbol{L}}_{n_{k}}X(\omega)=\boldsymbol{f}_{n_{k}}(X(\omega),X_{\sharp}\mathbb{P})$
for a.e. $\omega\in\Omega$,

so that, choosing representatives of $X$, ${\boldsymbol{L}}X$,
${\boldsymbol{L}}_{n_{k}}X$, we can find a full $\mathbb{P}$-measure set
$\Omega_{0}\subset\Omega$ such that

${\boldsymbol{L}}X(\omega)=\lim_{k\to+\infty}{\boldsymbol{L}}_{n_{k}}X(\omega)=\lim_{k\to+\infty}\boldsymbol{f}_{n_{k}}(X(\omega),X_{\sharp}\mathbb{P})\quad\text{
for every }\omega\in\Omega_{0}.$

Thus, for every $\omega\in\Omega_{0}$, $X(\omega)\in E$ and then
${\boldsymbol{L}}X(\omega)=f_{X}(X(\omega))$. If now $X^{\prime}\in D$ is such
that ${X^{\prime}}_{\sharp}\mathbb{P}=\mu$, arguing as in claim (c) of the
proof of Theorem 4.8, it is easy to see that
$\boldsymbol{f}_{X}=\boldsymbol{f}_{X^{\prime}}$ $\mu$-almost everywhere. ∎

We eventually apply the previous results to (possibly multivalued) maximal
$\lambda$-dissipative operators
$\boldsymbol{B}\subset\mathcal{X}\times\mathcal{X}$ (with $p=q=2$ and
$\mathsf{X}=\mathsf{Y}$) which are invariant by measure-preserving
isomorphisms. We show that the invariance property is inherited by the
associated resolvent operator, Yosida approximation, minimal selection and
semigroup which also enjoy a map-representation property.

We denote by
$\iota\big{(}D(\boldsymbol{B})\big{)}:=\Big{\\{}X_{\sharp}\mathbb{P}:X\in\mathrm{D}(\boldsymbol{B})\Big{\\}}$
the image in $\mathcal{P}_{2}(\mathsf{X})$ of the domain of $\boldsymbol{B}$.

######  Theorem 4.13 (Structure of resolvents, Yosida approximations and
semigroups).

Let $\mathsf{X}$ be a separable Hilbert space and let $p=2$. Let
$\boldsymbol{B}\subset\mathcal{X}\times\mathcal{X}$ be a maximal
$\lambda$-dissipative operator which is invariant by measure-preserving
isomorphisms. Then for every $0<\tau<1/\lambda^{+},\ t\geq 0$ the operators
$\boldsymbol{B},\boldsymbol{B}_{\tau},{\boldsymbol{J}_{\tau}},\boldsymbol{S}_{t},\boldsymbol{B}^{\circ}$
are law invariant. Moreover there exist (uniquely defined) continuous maps
$\boldsymbol{j}_{\tau}:\mathcal{S}(\mathsf{X})\to\mathsf{X}$,
$\boldsymbol{b}_{\tau}:\mathcal{S}(\mathsf{X})\to\mathsf{X}$, and
$\boldsymbol{s}_{t}:\mathcal{S}(\mathsf{X},\overline{\iota\big{(}D(\boldsymbol{B})\big{)}})\to\mathsf{X}$
such that

  1. (1)

for every $\mu\in\mathcal{P}_{2}(\mathsf{X})$, the map
$\boldsymbol{j}_{\tau}(\cdot,\mu):\operatorname{supp}(\mu)\to\mathsf{X}$ is
$(1-\lambda\tau)^{-1}$-Lipschitz continuous, for $0<\tau<1/\lambda^{+}$,

  2. (2)

for every $\mu\in\mathcal{P}_{2}(\mathsf{X})$, the map
$\boldsymbol{b}_{\tau}(\cdot,\mu):\operatorname{supp}(\mu)\to\mathsf{X}$ is
$\frac{2-\lambda\tau}{\tau(1-\lambda\tau)}$-Lipschitz continuous, for
$0<\tau<1/\lambda^{+}$,

  3. (3)

for every $\mu\in\overline{\iota\big{(}D(\boldsymbol{B})\big{)}}$, the map
$\boldsymbol{s}_{t}(\cdot,\mu):\operatorname{supp}(\mu)\to\mathsf{X}$ is
$e^{\lambda t}$-Lipschitz continuous,

and

$\displaystyle\text{ for every $X\in\mathcal{X}$,
}{\boldsymbol{J}_{\tau}}X(\omega)=\boldsymbol{j}_{\tau}(X(\omega),X_{\sharp}\mathbb{P})\text{
for $\mathbb{P}$-a.e.\leavevmode\nobreak\ $\omega\in\Omega$,}$ (4.17)
$\displaystyle\text{ for every $X\in\mathcal{X}$,
}\boldsymbol{B}_{\tau}X(\omega)=\boldsymbol{b}_{\tau}(X(\omega),X_{\sharp}\mathbb{P})\text{
for $\mathbb{P}$-a.e.\leavevmode\nobreak\ $\omega\in\Omega$,}$ (4.18)
$\displaystyle\text{ for every $X\in\overline{\mathrm{D}(\boldsymbol{B})}$,
}\boldsymbol{S}_{t}X(\omega)=\boldsymbol{s}_{t}(X(\omega),X_{\sharp}\mathbb{P})\text{
for $\mathbb{P}$-a.e.\leavevmode\nobreak\ $\omega\in\Omega$,}$ (4.19)

together with the invariance and semigroup properties

$\begin{gathered}\mu\in\overline{\iota\big{(}D(\boldsymbol{B})\big{)}}\quad\Rightarrow\quad\boldsymbol{s}_{t}(\cdot,\mu)_{\sharp}\mu\in\overline{\iota\big{(}D(\boldsymbol{B})\big{)}};\qquad\mu\in{\iota\big{(}D(\boldsymbol{B})\big{)}}\quad\Rightarrow\quad\boldsymbol{s}_{t}(\cdot,\mu)_{\sharp}\mu\in{\iota\big{(}D(\boldsymbol{B})\big{)}},\\\
\boldsymbol{s}_{t+h}(x,\mu)=\boldsymbol{s}_{h}(\boldsymbol{s}_{t}(x,\mu),\boldsymbol{s}_{t}(\cdot,\mu)_{\sharp}\mu)\quad\text{for
every
}(x,\mu)\in\mathcal{S}(\mathsf{X},\overline{\iota\big{(}D(\boldsymbol{B})\big{)}}),\quad
t,h\geq 0.\end{gathered}$ (4.20)

Finally, for every $\mu\in\iota\big{(}D(\boldsymbol{B})\big{)}$, there exists
a map $\boldsymbol{b}^{\circ}(\cdot,\mu)\in L^{2}(\mathsf{X},\mu;\mathsf{X})$
such that for every $X\in\mathcal{X}$

$\text{ if $X_{\sharp}\mathbb{P}=\mu$ then $X\in\mathrm{D}(\boldsymbol{B})$,
}\boldsymbol{B}^{\circ}X(\omega)=\boldsymbol{b}^{\circ}(X(\omega),\mu)\text{
for $\mathbb{P}$-a.e.\leavevmode\nobreak\ $\omega\in\Omega$.}$ (4.21)

For every $\mu\in\iota\big{(}D(\boldsymbol{B})\big{)}$, the map
$\boldsymbol{b}^{\circ}(\cdot,\mu)$ is $\lambda$-dissipative in a set
$\mathsf{X}_{0}\subset\mathsf{X}$ of full $\mu$-measure and satisfies

$\lim_{h\downarrow
0}\int\bigg{|}\frac{1}{h}(\boldsymbol{s}_{t+h}(x,\mu)-\boldsymbol{s}_{t}(x,\mu))-\boldsymbol{b}^{\circ}(\boldsymbol{s}_{t}(x,\mu),\boldsymbol{s}_{t}(\cdot,\mu)_{\sharp}\mu)\bigg{|}^{2}\,\mathrm{d}\mu(x)=0\quad
t\geq 0.$ (4.22)

######  Remark 4.14.

By Theorem 4.13, a maximal $\lambda$-dissipative operator
$\boldsymbol{B}\subset\mathcal{X}\times\mathcal{X}$, $\lambda\in\mathbb{R}$,
is law invariant if and only if it is invariant by measure-preserving
isomorphisms. Thanks to (4.21), also $D(\boldsymbol{B})$ is _law invariant_ ,
i.e. if $X\in D(\boldsymbol{B})$ and $Y\in\mathcal{X}$ is such that
$Y_{\sharp}\mathbb{P}=X_{\sharp}\mathbb{P}$, then $Y\in D(\boldsymbol{B})$.

###### Proof.

First, notice that $\boldsymbol{B}$ is closed being maximal
$\lambda$-dissipative, hence it is law invariant by Proposition 4.3. Recall
that ${\boldsymbol{J}_{\tau}}$ is everywhere defined,
$(1-\lambda\tau)^{-1}$-Lipschitz continuous for every $0<\tau<1/\lambda^{+}$
(see Section 2.3) and invariant by measure-preserving isomorphisms by
Proposition 2.10. Fixed $0<\tau<1/\lambda^{+}$, we can thus apply Proposition
4.3 and Theorem 4.8(1) and get that ${\boldsymbol{J}_{\tau}}$ is law invariant
together with property (1) in Theorem 4.13 and (4.17). Similarly,
$\boldsymbol{S}_{t}$ is defined on the closed set
$\overline{\mathrm{D}(\boldsymbol{B})}$, it is $e^{\lambda t}$-Lipschitz
continuous (cf. Section 2.3) and invariant by measure-preserving isomorphisms
by Proposition 2.10, so that we can apply Proposition 4.3 and Theorem 4.8(1)
and get that it is law invariant together with property (3) of Theorem 4.13
and (4.19). The content of (4.20) immediately follows by the semigroup and
invariance properties of $\boldsymbol{S}_{t}$ (cf. Section 2.3), also using
that
$\overline{\iota(\mathrm{D}(\boldsymbol{B}))}=\iota(\overline{\mathrm{D}(\boldsymbol{B})})$.

We now prove (4.21). If $X\in\mathrm{D}(\boldsymbol{B})$, we have that
$\boldsymbol{B}_{\tau}X\to\boldsymbol{B}^{\circ}X$ as $\tau\downarrow 0$,
moreover $\boldsymbol{B}_{\tau}$ is law invariant, everywhere defined and, for
every $0<\tau<1/\lambda^{+}$, $\boldsymbol{B}_{\tau}$ is
$\lambda/(1-\lambda\tau)$-dissipative and
$\frac{2-\lambda\tau}{\tau(1-\lambda\tau)}$-Lipschitz continuous (cf. Section
2.3). Hence, can apply Theorem 4.8(2) and get that $\boldsymbol{B}^{\circ}$ is
law invariant and that there exists, for every
$\mu\in\iota\big{(}D(\boldsymbol{B})\big{)}$, a map
$\boldsymbol{b}^{\circ}[\mu]\equiv\boldsymbol{b}^{\circ}(\cdot,\mu)\in
L^{2}(\mathsf{X},\mu;\mathsf{X})$ such that for every
$X\in\mathrm{D}(\boldsymbol{B})$

$\displaystyle\boldsymbol{B}^{\circ}X(\omega)=$
$\displaystyle\boldsymbol{b}^{\circ}[\mu](X(\omega))\text{ for
$\mathbb{P}$-a.e.\leavevmode\nobreak\ $\omega\in\Omega$},$ (4.23) there exists
$\tau_{k}\downarrow 0$ s.t.
$\displaystyle\boldsymbol{b}_{\tau_{k}}(X(\omega),\mu)\to\boldsymbol{b}^{\circ}[\mu](X(\omega))\text{
for $\mathbb{P}$-a.e.\leavevmode\nobreak\ $\omega\in\Omega$},$ (4.24)

where $\boldsymbol{b}_{\tau}$ is the (unique) continuous map that represents
$\boldsymbol{B}_{\tau}$ coming from Theorem 4.8(1), $0<\tau<1/\lambda^{+}$. To
complete the proof of (4.21), we have to check that, if
$\mu\in\iota\big{(}D(\boldsymbol{B})\big{)}$ and $X\in\mathcal{X}$ is such
that $X_{\sharp}\mathbb{P}=\mu$, then $X\in\mathrm{D}(\boldsymbol{B})$. Since
$\mu\in\iota\big{(}D(\boldsymbol{B})\big{)}$, there exists
$Y\in\mathrm{D}(\boldsymbol{B})$ such that
$Y_{\sharp}\mathbb{P}=X_{\sharp}\mathbb{P}=\mu$. By (4.18), we have

$|\boldsymbol{B}_{\tau}Y|^{2}_{\mathcal{X}}=\int|b_{\tau}(Y(\omega),\mu)|^{2}\mathrm{d}\mathbb{P}(\omega)=\int|b_{\tau}(x,\mu)|^{2}\mathrm{d}\mu(x)=\int|b_{\tau}(X(\omega),\mu)|^{2}\mathrm{d}\mathbb{P}(\omega)=|\boldsymbol{B}_{\tau}X|_{\mathcal{X}}^{2}.$

Hence, since $Y\in\mathrm{D}(\boldsymbol{B})$, by (2.24) we have

$(1-\lambda\tau)|\boldsymbol{B}_{\tau}X|_{\mathcal{X}}=(1-\lambda\tau)|\boldsymbol{B}_{\tau}Y|_{\mathcal{X}}\uparrow|\boldsymbol{B}^{\circ}Y|_{\mathcal{X}}<+\infty,\quad\text{as
}\tau\downarrow 0.$

Recalling (2.25), we get $X\in\mathrm{D}(\boldsymbol{B})$. Finally (4.22)
follows by (2.27) using (4.19) and (4.21).

It only remains to show that, for every
$\mu\in\iota\big{(}D(\boldsymbol{B})\big{)}$, the map
$\boldsymbol{b}^{\circ}[\mu]$ is $\lambda$-dissipative in a full $\mu$-measure
set. To this aim, observe that we can apply Theorem 4.8(1) to
$\boldsymbol{B}_{\tau}$, $0<\tau<1/\lambda^{+}$, so that for every
$\mu\in\mathcal{P}_{2}(\mathsf{X})$ we have

$\langle\boldsymbol{b}_{\tau}(y)-\boldsymbol{b}_{\tau}(y^{\prime}),y-y^{\prime}\rangle\leq\frac{\lambda}{1-\lambda\tau}|y-y^{\prime}|^{2}\quad\text{
for every }y,y^{\prime}\in\operatorname{supp}(\mu).$ (4.25)

Let $\mu\in\iota\big{(}D(\boldsymbol{B})\big{)}$ and let us consider a
representative of $X\in\mathrm{D}(\boldsymbol{B})$ such that
$X_{\sharp}\mathbb{P}=\mu$; let $\tau_{k}$ be a sequence as in (4.24) and let
$\Omega_{0}$ be a full $\mathbb{P}$-measure set where the convergence in
(4.24) takes place. If we take $y,y^{\prime}\in
X(\Omega_{0})\cap\operatorname{supp}(\mu)$ then we can find
$\omega,\omega^{\prime}\in\Omega_{0}$ such that $X(\omega)=y$ and
$X(\omega^{\prime})=y^{\prime}$ so that, passing to the limit as $k\to+\infty$
in (4.25) written for $\tau=\tau_{k}$ we get that

$\langle\boldsymbol{b}^{\circ}[\mu](y)-\boldsymbol{b}^{\circ}[\mu](y^{\prime}),y-y^{\prime}\rangle\leq\lambda|y-y^{\prime}|^{2}\quad\text{
for every }y,y^{\prime}\in X(\Omega_{0})\cap\operatorname{supp}(\mu).$

It is then enough to observe that $X(\Omega_{0})\cap\operatorname{supp}(\mu)$
contains a Borel set $\mathsf{X}_{0}$ of full $\mu$-measure to conclude: in
fact, being $X(\Omega_{0})$ a Suslin set, we can find two Borel sets
$E_{0},E_{1}$ such that $E_{0}\subset X(\Omega_{0})\subset E_{1}$ and
$\mu(E_{1}\setminus E_{0})=0$. Since
$\mu(E_{1})=X_{\sharp}\mathbb{P}(E_{1})\geq\mathbb{P}(\Omega_{0})=1$, we
conclude that we can take $\mathsf{X}_{0}:=E_{0}\cap\operatorname{supp}(\mu)$.
∎

##  Appendix A An alternative proof of the extension theorem for
$\mathsf{G}$-invariant Lipschitz maps

We provide an alternative proof to the extension result for invariant
Lipschitz maps stated in Theorem 2.13. Here we use a recent and beautiful
explicit construction of a Lipschitz extension provided by [ALM21] (see also
[ALM18]) and we show that it preserves the invariance. We consider a Hilbert
space $\mathcal{H}$ with norm $|\cdot|$ and scalar product
$\langle\cdot,\cdot\rangle$ as in Section 2.4.

Recall that if $g:\mathsf{W}\to\mathbb{R}$ is a function defined in a Hilbert
space $\mathsf{W}$, then its convex envelope is defined by

$\operatorname{co}(g)(w):=\inf\Big{\\{}\sum_{i=1}^{N}\alpha_{i}g(w_{i}):\alpha_{i}\geq
0,\ \sum_{i=1}^{N}\alpha_{i}=1,\ w_{i}\in\mathsf{W},\
\sum_{i=1}^{N}\alpha_{i}w_{i}=w,\ N\in\mathbb{N}\Big{\\}}.$ (A.1)

If $g$ is locally bounded from above, then also $\operatorname{co}(g)$ is
locally bounded from above and it is therefore locally Lipschitz.

######  Theorem A.1 ([ALM21]).

Let $f:D\to\mathcal{H}$ be a $L$-Lipschitz map defined in
$D\subset\mathcal{H}$. Setting for every $x,y\in\mathcal{H}$

$\displaystyle g(x,y):={}$ $\displaystyle\inf_{x^{\prime}\in
D}\Big{\\{}\langle
f(x^{\prime}),y\rangle+\frac{L}{2}|(x-x^{\prime},y)|_{\mathcal{H}\times\mathcal{H}}^{2}\Big{\\}}+\frac{L}{2}|(x,y)|_{\mathcal{H}\times\mathcal{H}}^{2},$
(A.2) $\displaystyle\tilde{g}:={}$ $\displaystyle\operatorname{co}(g),$

then $\tilde{g}$ is a convex function of class $\mathrm{C}^{1,1}$ in
$\mathcal{H}\times\mathcal{H}$ and its partial differential with respect to
the second variable in $\mathcal{H}$

$F(x):=\nabla_{y}\tilde{g}(x,0),\quad x\in\mathcal{H},$ (A.3)

is a $L$-Lipschitz extension of $f$.

We state our first result concerning the extension of
$\mathsf{G}_{\mathcal{H}}$-invariant Lipschitz maps.

######  Theorem A.2.

Under the same assumption of Theorem A.1, let us also suppose that $f$ is
$\mathsf{G}_{\mathcal{H}}$-invariant according to (2.29). Then $F$ is
$\mathsf{G}_{\mathcal{H}}$-invariant as well. In particular, any
$\mathsf{G}_{\mathcal{H}}$-invariant $L$-Lipschitz function
$f:D\to\mathcal{H}$ defined in a subset $D$ of $\mathcal{H}$ admits a
$\mathsf{G}_{\mathcal{H}}$-invariant $L$-Lipschitz extension
$F:\mathcal{H}\to\mathcal{H}.$

###### Proof.

We divide the proof in several claims.

_Claim 1: the map $g$ is $\mathsf{G}_{\mathcal{H}}$-invariant, i.e.
$g(Ux,Uy)=g(x,y)$ for every $(x,y)\in\mathcal{H}\times\mathcal{H}$ and
$U\in\mathsf{G}_{\mathcal{H}}.$_

Every element $x^{\prime}\in D$ can be written as
$x^{\prime}=Ux^{\prime\prime}$ with $x^{\prime\prime}=U^{-1}x^{\prime}\in D$,
so that for every $U\in\mathsf{G}_{\mathcal{H}}$ and
$(x,y)\in\mathcal{H}\times\mathcal{H}$

$\displaystyle g(Ux,Uy)$ $\displaystyle=\inf_{x^{\prime}\in D}\Big{\\{}\langle
f(x^{\prime}),Uy\rangle+\frac{L}{2}|(Ux-x^{\prime},Uy)|_{\mathcal{H}\times\mathcal{H}}^{2}\Big{\\}}+\frac{L}{2}|(Ux,Uy)|_{\mathcal{H}\times\mathcal{H}}^{2}$
$\displaystyle=\inf_{x^{\prime\prime}\in D}\Big{\\{}\langle
f(Ux^{\prime\prime}),Uy\rangle+\frac{L}{2}|(U(x-x^{\prime\prime}),Uy)|_{\mathcal{H}\times\mathcal{H}}^{2}\Big{\\}}+\frac{L}{2}|(Ux,Uy)|_{\mathcal{H}\times\mathcal{H}}^{2}$
$\displaystyle=\inf_{x^{\prime\prime}\in D}\Big{\\{}\langle
f(x^{\prime\prime}),y\rangle+\frac{L}{2}|(x-x^{\prime\prime},y)|_{\mathcal{H}\times\mathcal{H}}^{2}\Big{\\}}+\frac{L}{2}|(x,y)|_{\mathcal{H}\times\mathcal{H}}^{2}$
$\displaystyle=g(x,y)$

where we used (2.29) and the isometric character of $U$ to get $\langle
f(Ux^{\prime\prime}),Uy\rangle=\langle Uf(x^{\prime\prime}),Uy\rangle=\langle
f(x^{\prime\prime}),y\rangle$.

_Claim 2: the map $\tilde{g}:=\operatorname{co}(g)$ is
$\mathsf{G}_{\mathcal{H}}$-invariant as well._

It is sufficient to observe that for every $U\in\mathsf{G}_{\mathcal{H}},$
$N\in\mathbb{N}$, and $\alpha_{i}\geq 0$ with $\sum_{i=1}^{N}\alpha_{i}=1$, a
collection
$\\{(x_{i},y_{i})\\}_{i=1}^{N}\in(\mathcal{H}\times\mathcal{H})^{N}$ satisfies
$\sum_{i=1}^{N}\alpha_{i}(x_{i},y_{i})=(x,y)$ if and only if
$\sum_{i=1}^{N}\alpha_{i}(Ux_{i},Uy_{i})=(Ux,Uy)$. Using (A.1) and the
invariance of $g$ we thus obtain $\tilde{g}(Ux,Uy)=\tilde{g}(x,y)$ for every
$x,y\in\mathcal{H}$.

_Claim 3: the map $F:=\nabla_{y}\tilde{g}(\cdot,0)$ is
$\mathsf{G}_{\mathcal{H}}$-invariant._

Since we know that $\tilde{g}$ is Frechét differentiable, we observe that
$z=F(x)$ if and only if

$\tilde{g}(x,y)=\tilde{g}(x,0)+\langle z,y\rangle+o(|y|)\quad\text{as }y\to
0.$

Being $\tilde{g}$ invariant, for every $U\in\mathsf{G}_{\mathcal{H}}$ and
$x\in\mathcal{H}$ the above formula immediately yields

$\displaystyle\tilde{g}(Ux,y)$
$\displaystyle=\tilde{g}(Ux,UU^{-1}y)=\tilde{g}(x,U^{-1}y)$
$\displaystyle=\tilde{g}(x,0)+\langle z,U^{-1}y\rangle+o(|U^{-1}y|)$
$\displaystyle=\tilde{g}(Ux,0)+\langle Uz,y\rangle+o(|y|)\quad\text{as }y\to
0,$

so that $F(Ux)=Uz=UF(x)$. ∎

## References

  * [AA99] Giovanni Alberti and Luigi Ambrosio  “A geometrical approach to monotone functions in ${\bf R}^{n}$”  In _Math. Z._ 230.2, 1999, pp. 259–316  DOI: 10.1007/PL00004691
  * [AGS08] Luigi Ambrosio, Nicola Gigli and Giuseppe Savaré  “Gradient Flows In Metric Spaces and in the Space of Probability Measures”, Lectures in Mathematics. ETH Zürich  Birkhäuser Basel, 2008
  * [ALM18] D. Azagra, E. Le Gruyer and C. Mudarra  “Explicit formulas for $C^{1,1}$ and $C_{\rm conv}^{1,\omega}$ extensions of 1-jets in Hilbert and superreflexive spaces”  In _J. Funct. Anal._ 274.10, 2018, pp. 3003–3032  DOI: 10.1016/j.jfa.2017.12.007
  * [ALM21] Daniel Azagra, Erwan Le Gruyer and Carlos Mudarra  “Kirszbraun’s theorem via an explicit formula”  In _Canad. Math. Bull._ 64.1, 2021, pp. 142–153  DOI: 10.4153/S0008439520000314
  * [Amb03] Luigi Ambrosio  “Lecture notes on optimal transport problems”  In _Mathematical aspects of evolving interfaces (Funchal, 2000)_ 1812, Lecture Notes in Math.  Springer, Berlin, 2003, pp. 1–52  DOI: 10.1007/978-3-540-39189-0“˙1
  * [Bar10] Viorel Barbu  “Nonlinear differential equations of monotone types in Banach spaces”, Springer Monographs in Mathematics  Springer, New York, 2010, pp. x+272  DOI: 10.1007/978-1-4419-5542-5
  * [BG03] Yann Brenier and Wilfrid Gangbo  “$L^{p}$ approximation of maps by diffeomorphisms”  In _Calc. Var. Partial Differential Equations_ 16.2, 2003, pp. 147–164  DOI: 10.1007/s005260100144
  * [Bog06] V.I. Bogachev  “Measure Theory”  Springer Berlin Heidelberg, 2006
  * [Bor10] Jonathan M. Borwein  “Fifty years of maximal monotonicity”  In _Optim. Lett._ 4.4, 2010, pp. 473–490  DOI: 10.1007/s11590-010-0178-x
  * [Bré73] Haim Brézis  “Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert” North-Holland Mathematics Studies, No. 5. Notas de Matemática (50)  North-Holland Publishing Co., Amsterdam-London; American Elsevier Publishing Co., Inc., New York, 1973, pp. vi+183
  * [BS03] Regina Sandra Burachik and B.. Svaiter  “Maximal monotonicity, conjugation and the duality product”  In _Proc. Amer. Math. Soc._ 131.8, 2003, pp. 2379–2383  DOI: 10.1090/S0002-9939-03-07053-9
  * [BW09] Heinz H. Bauschke and Xianfu Wang  “The kernel average for two convex functions and its application to the extension and representation of monotone operators”  In _Trans. Amer. Math. Soc._ 361.11, 2009, pp. 5947–5965  DOI: 10.1090/S0002-9947-09-04698-4
  * [BW10] Heinz H. Bauschke and Xianfu Wang  “Firmly nonexpansive and Kirszbraun-Valentine extensions: a constructive approach via monotone operator theory”  In _Nonlinear analysis and optimization I. Nonlinear analysis_ 513, Contemp. Math.  Amer. Math. Soc., Providence, RI, 2010, pp. 55–64  DOI: 10.1090/conm/513/10075
  * [Car13] Pierre Cardaliaguet  “Notes on Mean Field Games” From P.-L. Lions’ lectures at College de France, link to the notes., 2013
  * [CD18] René Carmona and François Delarue  “Probabilistic theory of mean field games with applications. I” Mean field FBSDEs, control, and games 83, Probability Theory and Stochastic Modelling  Springer, Cham, 2018, pp. xxv+713
  * [CSS] Giulia Cavagnari, Giuseppe Savaré and Giacomo E. Sodini  “A Lagrangian approach to totally dissipative evolutions in Wasserstein spaces” arXiv preprint 2023
  * [DF64] Hans Debrunner and Peter Flor  “Ein Erweiterungssatz für monotone Mengen”  In _Arch. Math._ 15, 1964, pp. 445–447  DOI: 10.1007/BF01589229
  * [Fit88] Simon Fitzpatrick  “Representing monotone operators by convex functions”  In _Workshop/Miniconference on Functional Analysis and Optimization (Canberra, 1988)_ 20, Proc. Centre Math. Anal. Austral. Nat. Univ.  Austral. Nat. Univ., Canberra, 1988, pp. 59–65
  * [FSS22] Massimo Fornasier, Giacomo E. Sodini and Giuseppe Savaré  “Density of subalgebras of Lipschitz functions in metric Sobolev spaces and applications to Sobolev-Wasserstein spaces”  In _arXiv:2209.00974_ , 2022
  * [Gho08] Nassif Ghoussoub  “A variational theory for monotone vector fields”  In _J. Fixed Point Theory Appl._ 4.1, 2008, pp. 107–135  DOI: 10.1007/s11784-008-0083-4
  * [GT19] Wilfrid Gangbo and Adrian Tudorascu  “On differentiability in the Wasserstein space and well-posedness for Hamilton-Jacobi equations”  In _J. Math. Pures Appl. (9)_ 125, 2019, pp. 119–174  DOI: 10.1016/j.matpur.2018.09.003
  * [Kir34] M.. Kirszbraun  “Über die zusammenziehenden und Lipschitzschen Transformationen”  In _Fundam. Math._ 22, 1934, pp. 77–108  DOI: 10.4064/fm-22-1-77-108
  * [MT01] Juan-Enrique Martinez-Legaz and Michel Théra  “A convex representation of maximal monotone operators” Special issue for Professor Ky Fan  In _J. Nonlinear Convex Anal._ 2.2, 2001, pp. 243–247
  * [Pen04] Jean-Paul Penot  “A representation of maximal monotone operators by closed convex functions and its impact on calculus rules”  In _C. R. Math. Acad. Sci. Paris_ 338.11, 2004, pp. 853–858  DOI: 10.1016/j.crma.2004.03.017
  * [Pen04a] Jean-Paul Penot  “The relevance of convex analysis for the study of monotonicity”  In _Nonlinear Anal._ 58.7-8, 2004, pp. 855–871  DOI: 10.1016/j.na.2004.05.018
  * [Pra07] Aldo Pratelli  “On the equality between Monge’s infimum and Kantorovich’s minimum in optimal mass transportation”  In _Annales de l’Institut Henri Poincare (B) Probability and Statistics_ 43, 2007, pp. 1–13  DOI: 10.1016/j.anihpb.2005.12.001
  * [Roc97] R. Rockafellar  “Convex analysis” Reprint of the 1970 original, Princeton Paperbacks, Princeton Landmarks in Mathematics  Princeton University Press, Princeton, NJ, 1997, pp. xviii+451
  * [Roy88] H.. Royden  “Real analysis”  Macmillan Publishing Company, New York, 1988, pp. xx+444
  * [SS22] Giuseppe Savaré and Giacomo E. Sodini  “A simple relaxation approach to duality for optimal transport problems in completely regular spaces”  In _J. Convex Anal._ 29.1, 2022, pp. 1–12
  * [Str11] Daniel W. Stroock  “Probability theory” An analytic view  Cambridge University Press, Cambridge, 2011, pp. xxii+527
  * [Val43] F.. Valentine  “On the extension of a vector function so as to preserve a Lipschitz condition”  In _Bull. Amer. Math. Soc._ 49, 1943, pp. 100–108  DOI: 10.1090/S0002-9904-1943-07859-7
  * [Val45] F.. Valentine  “A Lipschitz condition preserving extension for a vector function”  In _Amer. J. Math._ 67, 1945, pp. 83–93  DOI: 10.2307/2371917
  * [Vis17] Augusto Visintin  “On the variational representation of monotone operators”  In _Discrete Contin. Dyn. Syst. Ser. S_ 10.4, 2017, pp. 909–918  DOI: 10.3934/dcdss.2017046


# Stochastic Solutions for Dense Subgraph Discovery
in Multilayer Networks

Yasushi Kawase  The University of TokyoHongo 7-3-1TokyoJapan113-8654
<EMAIL_ADDRESS>,  Atsushi Miyauchi  The University of TokyoHongo
7-3-1TokyoJapan113-8654<EMAIL_ADDRESS>and  Hanna Sumita  Tokyo
Institute of TechnologyOookayama 2-12-1TokyoJapan<EMAIL_ADDRESS>
(2023)

###### Abstract.

Network analysis has played a key role in knowledge discovery and data mining.
In many real-world applications in recent years, we are interested in mining
_multilayer networks_ , where we have a number of edge sets called _layers_ ,
which encode different types of connections and/or time-dependent connections
over the same set of vertices. Among many network analysis techniques, dense
subgraph discovery, aiming to find a dense component in a network, is an
essential primitive with a variety of applications in diverse domains. In this
paper, we introduce a novel optimization model for dense subgraph discovery in
multilayer networks. Our model aims to find a stochastic solution, i.e., a
probability distribution over the family of vertex subsets, rather than a
single vertex subset, whereas it can also be used for obtaining a single
vertex subset. For our model, we design an LP-based polynomial-time exact
algorithm. Moreover, to handle large-scale networks, we also devise a simple,
scalable preprocessing algorithm, which often reduces the size of the input
networks significantly and results in a substantial speed-up. Computational
experiments demonstrate the validity of our model and the effectiveness of our
algorithms.

network analysis, multilayer networks, dense subgraph discovery, stochastic
solutions

††journalyear: 2023††copyright: acmcopyright††conference: Proceedings of the
Sixteenth ACM International Conference on Web Search and Data Mining; February
27–March 3, 2023; Singapore, Singapore††booktitle: Proceedings of the
Sixteenth ACM International Conference on Web Search and Data Mining (WSDM
’23), February 27–March 3, 2023, Singapore, Singapore††price: 15.00††doi:
10.1145/3539597.3570444††isbn: 978-1-4503-9407-9/23/02††submissionid: ††ccs:
Theory of computation Graph algorithms analysis††ccs: Theory of computation
Mathematical optimization

##  1\. Introduction

Network analysis has played a key role in knowledge discovery and data mining.
In many real-world applications in recent years, we are interested in mining
_multilayer networks_ rather than ordinary (i.e., single-layer) networks,
where we have a number of edge sets called _layers_ , which encode different
types of connections and/or time-dependent connections over the same set of
vertices (Boccaletti et al., 2014; De Domenico et al., 2013; Kivelä et al.,
2014). For example, in the Twitter network, there are various layers
representing different types of connections, e.g., follower–followee
relations, retweets, and mentions, among users. Moreover, each of those
connections is time-dependent and therefore leads to multiple layers by
itself. As another example, consider brain networks arising in neuroscience,
where vertices correspond to small regions of a brain (Friston, 2011). In this
network, we can obtain at least two layers representing the structural
connectivity and the functional connectivity (e.g., co-activation) among the
small pieces of a brain.

Among many network analysis techniques, dense subgraph discovery, aiming to
find a dense component in a network, is an essential primitive with a variety
of applications in diverse domains (Gionis and Tsourakakis, 2015; Lee et al.,
2010). Examples include detecting communities and spam link farms in Web
graphs (Dourisboure et al., 2007; Gibson et al., 2005), experts extraction in
crowdsourcing systems (Kawase et al., 2019), real-time story identification in
microblogging streams (Angel et al., 2012), and extracting molecular complexes
in protein–protein interaction networks (Bader and Hogue, 2003).

Figure 1. An optimal solution to our model with the regret metric for the
real-world multilayer network called WILDBIRDS: selecting the black vertices
with probability 0.29 and the union of the gray and black vertices with
probability 0.71. There are six layers, and the color of each edge represents
the layer containing the edge.

Recently, dense subgraph discovery has been extended from single-layer
networks to multilayer networks. Jethava and Beerenwinkel (Jethava and
Beerenwinkel, 2015) introduced the _densest common subgraph problem_ , where
given a multilayer network, we are asked to find a vertex subset that
maximizes the minimum _degree density_ (defined later) over the layers. They
mentioned some concrete applications of this problem in biological networks.
Moreover, we can see that the densest common subgraph problem may appear in
the context of robust optimization. Indeed, multilayer networks can be seen as
a model of single-layer networks with uncertain edges with a number of
scenarios, where each layer corresponds to one scenario. If we can find a
subgraph that is reasonably dense for all layers, the solution is more robust
than that obtained by single-layer network analysis, which may be arbitrarily
bad for some layers.

###  1.1. Our contribution

In this paper, we introduce a novel optimization model for dense subgraph
discovery in multilayer networks. Given an edge-weighted multilayer network,
our model aims to find a vertex subset that is dense for the layer selected by
an adversary. Specifically, we employ a stochastic solution, i.e., a
probability distribution over the family of vertex subsets, while existing
work focuses only on a deterministic solution, i.e., a single vertex subset.
More precisely, we stochastically choose a vertex subset according to a
stochastic solution and then, knowing only the stochastic solution but not a
realization, the adversary selects the worst layer for us. Note that the worst
layer with respect to each realization of the stochastic solution may be
different from that with respect to the stochastic solution. It is worth
mentioning that our model can be seen as a zero-sum two-player Stackelberg
game (Başar and Olsder, 1999), where the leader is our algorithm and the
follower is the adversary.

We measure the density of a vertex subset using the quality function called
the _degree density_ (or simply _density_) similarly to Jethava and
Beerenwinkel (Jethava and Beerenwinkel, 2015). The degree density of a vertex
subset (in a single-layer network) is defined as half the average degree of
the subgraph induced by the subset. This quality function is often used in the
literature of dense subgraph discovery; in fact, this is the objective
function of the well-known _densest subgraph problem_ (see Section 2).

We evaluate the performance of a stochastic solution (or an algorithm) using
the following three metrics for the layer selected by the adversary: (i) the
_density_ , i.e., the expected degree density, (ii) the _robust ratio_ , i.e.,
the ratio of the expected density to the optimal density, and (iii) the
_regret_ , i.e., the difference between the optimal density and the expected
density. Therefore, we address three optimization problems according to the
metrics.

Our main algorithmic contribution is to design a polynomial-time exact
algorithm for our optimization model, which is applicable to all of the above
three metrics. Specifically, the algorithm first solves an LP, which is a
generalization of the LP used for the densest common subgraph problem (Jethava
and Beerenwinkel, 2015) and then computes an optimal probability distribution
based on the LP’s optimal solution.

We observe that the output of our algorithm has a useful structure; the family
of the vertex subsets with positive probabilities has a hierarchical
structure. This leads to several practical benefits, e.g., the largest size
subset contains all the other subsets and the optimal solution obtained by our
algorithm has support size at most the number of vertices (although a solution
of our optimization model may have support size exponential in the number of
vertices). Moreover, we can demonstrate that the support size of the solution
obtained is upper bounded by the number of layers. It is desirable to have a
small support size for the understandability of stochastic solutions and for
purposes of simple verification and validation.

For practical use of our proposed algorithm, we wish to speed up our algorithm
so that it is applicable to larger-scale multilayer networks. The bottleneck
is the computation cost of the LP, which has a lot of variables and
constraints. To this end, we devise a simple, scalable preprocessing
algorithm, which often reduces the size of the input networks significantly
and results in a substantial speed-up. Specifically, the algorithm first
computes an approximate solution by solving a much smaller LP than the
aforementioned LP and then removes vertices from the original network using
the information of the approximate solution obtained. Note that our
preprocessing algorithm is a generalization of that proposed by Balalau et al.
(Balalau et al., 2015) for the densest subgraph problem. However, to verify
that the preprocessing algorithm does not harm any optimal solution, we need a
more sophisticated analysis, which is totally different from theirs.

Finally, we conduct thorough computational experiments using synthetic graphs
and real-world networks to verify the validity of our model and to evaluate
the performance of our algorithms. Figure 1 illustrates an optimal solution to
our model with the regret metric for the real-world multilayer network called
WILDBIRDS (see Section 7 for its detailed characteristics). The optimal
solution obtained is reasonably dense for all layers, as desired.

It should be remarked that stochastic solutions can also be used for obtaining
a single vertex subset. For example, the following three rules are reasonable
to select a vertex subset from a stochastic solution: (i) stochastically
select a subset following the probability distribution, (ii) select a subset
with the highest probability, and (iii) select a subset with non-zero
probability that optimizes some metric at hand (e.g., the minimum density,
minimum robust ratio, and maximum regret, over layers). In our experiments, we
compare the vertex subsets detected using the above rules with those obtained
by existing algorithms for the densest common subgraph problem, in terms of
various evaluation metrics.

##  2\. Related Work

The densest subgraph problem is one of the most popular optimization models
for dense subgraph discovery. Let $G=(V,E)$ be an undirected graph and
$w\colon E\to\mathbb{R}_{++}$ a positive edge weight. For a vertex subset
$S\subseteq V$, the subgraph induced by $S$ is denoted by
$G[S]\coloneqq(S,E[S])$, where $E[S]\coloneqq\bigl{\\{}\\{u,v\\}\in E\mid
u,v\in S\bigr{\\}}$. In addition, we denote by $w(S)$ the total weight of the
edges in $S$, i.e., $w(S)=\sum_{e\in E[S]}w(e)$. For a nonempty $S\subseteq
V$, the _degree density_ (or simply called _density_) of $S$ is defined as
$w(S)/|S|$ (where we define the density of the empty set (i.e., $0/0$) to be
$0$). In the densest subgraph problem, given a graph $G=(V,E)$ with an edge
weight $w$, we are asked to find $S\subseteq V$ that maximizes the density
$w(S)/|S|$.

It is well known that the densest subgraph problem can be solved exactly in
polynomial time using a maximum-flow-based algorithm (Goldberg, 1984) or an
LP-based algorithm (Charikar, 2000). Moreover, it was shown that a simple
greedy algorithm called the _greedy peeling_ admits $2$-approximation in
$O(m+n\log n)$ time (Charikar, 2000; Kortsarz and Peleg, 1994). Recently, Boob
et al. (Boob et al., 2020) designed an iterative greedy peeling algorithm, and
demonstrated empirically that the output tends to be nearly optimal. Later,
Chekuri, Quanrud, and Torres (Chekuri et al., 2022) proved the convergence to
optimality of this algorithm (in a more general context).

Balalau et al. (Balalau et al., 2015) introduced a simple preprocessing
algorithm to improve the scalability of (exact) algorithms for the densest
subgraph problem. Their preprocessing algorithm first computes an approximate
solution $S\subseteq V$ (using the greedy peeling), and then iteratively
removes a vertex with the (weighted) degree less than the objective value of
the approximate solution obtained. The validity of this preprocessing
algorithm is guaranteed by the fact that any vertex with the (weighted) degree
less than the optimal value is not contained in any optimal solution (Balalau
et al., 2015).

For the densest common subgraph problem, Jethava and Beerenwinkel (Jethava and
Beerenwinkel, 2015) devised an LP-based polynomial-time heuristic and a
$2k$-approximation algorithm based on the greedy peeling, where $k$ is the
number of layers. Later, Charikar, Naamad, and Wu (Charikar et al., 2018)
designed two polynomial-time algorithms with approximation ratios
$O(\sqrt{|V|\log k})$ and $O(|V|^{2/3})$ (irrespective of $k$), respectively.
Moreover, they showed some strong inapproximability results for the problem,
based on some reasonable computational complexity assumptions. Thus, it is
very unlikely that a well-approximate solution can be found in polynomial
time. In contrast to this, as mentioned above, we can compute an _optimal_
stochastic solution in terms of the aforementioned three metrics in polynomial
time.

Recently, Galimberti, Bonchi, and Gullo (Galimberti et al., 2017) introduced a
generalization of the densest common subgraph problem, which they refer to as
the _multilayer densest subgraph problem_. This problem exploits a trade-off
between the minimum density value over layers and the number of layers
considered. They proposed an approximation algorithm using a core
decomposition technique for multilayer networks. Very recently, Hashemi,
Behrouz, and Lakshmanan (Hashemi et al., 2022) designed a sophisticated core
decomposition algorithm, which they call the _FirmCore decomposition
algorithm_. Their algorithm finds the set of $(k,\lambda)$-FirmCores for all
possible $k$ and $\lambda$ in polynomial time, where $(k,\lambda)$-FirmCore is
a maximal subgraph in which every vertex has degree no less than $k$ in the
subgraph for at least $\lambda$ layers. They demonstrated that the
decomposition unfolds a better solution to the multilayer densest subgraph
problem than the algorithm by Galimberti, Bonchi, and Gullo (Galimberti et
al., 2017) for many instances.

Semertzidis et al. (Semertzidis et al., 2019) introduced another
generalization of the densest common subgraph problem, called the Best Friends
Forever (BFF) problem, in the context of evolving graphs with a number of
snapshots. The BFF problem is a series of optimization problems that maximize
an _aggregate density_ over snapshots, where the aggregate density is set to
be the average/minimum value of the average/minimum degree of vertices over
layers. They investigated the computational complexity of the problems and
designed some approximation or heuristic algorithms.

Recalling that multilayer networks can also be seen as a model of networks
with _uncertainty_ , we can find some other optimization models related. Zou
(Zou, 2013) studied the densest subgraph problem in _uncertain graphs_. An
uncertain graph is a pair of $G=(V,E)$ and $p\colon E\rightarrow[0,1]$, where
$e\in E$ is present with probability $p(e)$ whereas $e\in E$ is absent with
probability $1-p(e)$. In the problem, given an uncertain graph $G=(V,E)$ with
$p$, we seek $S\subseteq V$ that maximizes the expected density. Zou (Zou,
2013) showed that this problem can be reduced to the original densest subgraph
problem, and designed a polynomial-time exact algorithm based on the
reduction. Recently, Tsourakakis et al. (Tsourakakis et al., 2019) introduced
a more general optimization problem called the _risk-averse dense subgraph
discovery_.

As another example, Miyauchi and Takeda (Miyauchi and Takeda, 2018) introduced
an optimization problem called the _robust densest subgraph problem_. In this
problem, given an undirected graph $G=(V,E)$ and an edge-weight space
$I=\times_{e\in E}[l_{e},r_{e}]\subseteq\times_{e\in E}[0,\infty)$, we are
asked to find $S\subseteq V$ that maximizes $\min_{w\in
I}\frac{w(S)/|S|}{w(S^{*}_{w})/|S^{*}_{w}|}$, where $S^{*}_{w}$ is an optimal
solution to the densest subgraph problem for $G$ with $w$. The intuition of
this problem is the same as that of ours; this problem also seeks $S\subseteq
V$ that is reasonably dense for any $w\in I$. However, they considered only
deterministic solutions and gave a strong hardness result.

As well as dense subgraph discovery, many important primitives for single-
layer network analysis have recently been extended to multilayer networks.
Examples include community detection (Bazzi et al., 2016; De Bacco et al.,
2017; Interdonato et al., 2017; Tagarelli et al., 2017), link prediction (De
Bacco et al., 2017; Jalili et al., 2017), analyzing spreading processes (De
Domenico et al., 2016; Salehi et al., 2015), and identifying central vertices
(Basaras et al., 2019; De Domenico et al., 2015).

##  3\. Model

In this section, we formally define our optimization model. Let
$G=(V,(E_{i})_{i\in[k]})$ be a multilayer network consisting of $k$ layers,
and let $w_{i}\colon E_{i}\to\mathbb{R}_{++}$ be a positive edge weight for
layer $i$. We denote by $E$ the union of all edge sets, i.e.,
$E\coloneqq\bigcup_{i\in[k]}E_{i}$. Let $S_{i}^{*}$ be a densest subgraph for
layer $i$, i.e., $S_{i}^{*}\in\operatorname*{arg\,max}_{S\subseteq
V}w_{i}(S)/|S|$. Our task is to find a vertex subset that is dense for the
layer selected adversarially. As mentioned in the introduction, we consider a
stochastic solution to compete with the adversary. Let $\Delta(2^{V})$ be the
set of probability distributions over $2^{V}$. For each $p\in\Delta(2^{V})$,
we denote by $p_{S}$ the probability of choosing $S\subseteq V$.

We aim to compute $p\in\Delta(2^{V})$ that maximizes some metric (when the
adversary selects the worst layer to $p$). We employ the following three
metrics:

Density. The first metric is the degree density itself. Specifically, when we
select a vertex subset according to a probability distribution
$p\in\Delta(2^{V})$ and the adversary selects a layer $i\in[k]$, our metric is
defined as follows:

(1) $\displaystyle\mathbb{E}_{S\sim
p}\left[\frac{w_{i}(S)}{|S|}\right]\quad\Bigg{(}=\sum_{S\subseteq
V}p_{S}\frac{w_{i}(S)}{|S|}\Bigg{)}.$

As the adversary selects the worst layer, we aim to find $p\in\Delta(2^{V})$
that maximizes the minimum of the density (1) among $i\in[k]$. Our
optimization model with this metric can be seen as a stochastic version of the
densest common subgraph problem introduced by Jethava and Beerenwinkel
(Jethava and Beerenwinkel, 2015).

Robust ratio. The _robust ratio_ is a metric based on the ratio of the
expected degree to the optimal degree. For a nonempty $S\subseteq V$ and
$i\in[k]$, let us consider a normalized density defined as
$\frac{w_{i}(S)/|S|}{w_{i}(S_{i}^{*})/|S_{i}^{*}|}$. When we select a vertex
subset according to a probability distribution $p\in\Delta(2^{V})$ and the
adversary selects a layer $i\in[k]$, the metric is defined as follows:

(2) $\displaystyle\mathbb{E}_{S\sim
p}\left[\frac{w_{i}(S)/|S|}{w_{i}(S_{i}^{*})/|S_{i}^{*}|}\right]\quad\Bigg{(}=\sum_{S\subseteq
V}p_{S}\frac{w_{i}(S)/|S|}{w_{i}(S_{i}^{*})|S_{i}^{*}|}\Bigg{)}.$

In other words, the robust ratio is equivalent to the density for the
multilayer network with weights $w_{1}^{\prime},\dots,w_{k}^{\prime}$ given by
$w_{i}^{\prime}(e)\coloneqq\frac{w_{i}(e)}{w_{i}(S_{i}^{*})/|S_{i}^{*}|}$ for
each $i\in[k]$ and $e\in E_{i}$. As the adversary selects the worst layer, we
aim to find $p\in\Delta(2^{V})$ that maximizes the minimum of (2) among
$i\in[k]$. Note that the optimal robust ratio is contained in the interval
$[1/k,\,1]$ because $p\in\Delta(2^{V})$ such that $p_{S_{i}^{*}}=1/k$ for each
$i\in[k]$ has the objective value of $1/k$.

Regret. The _regret_ is a metric based on the difference between the optimal
density and the expected density. For $S\subseteq V$ and $i\in[k]$, the regret
is defined as $w_{i}(S^{*}_{i})/|S^{*}_{i}|-w_{i}(S)/|S|$. When we select a
vertex subset according to a probability distribution $p\in\Delta(2^{V})$ and
the adversary selects a layer $i\in[k]$, the metric is defined as follows:

(3) $\displaystyle\mathbb{E}_{S\sim
p}\left[\frac{w_{i}(S_{i}^{*})}{|S_{i}^{*}|}-\frac{w_{i}(S)}{|S|}\right]\quad\Bigg{(}=\frac{w_{i}(S_{i}^{*})}{|S_{i}^{*}|}-\sum_{S\subseteq
V}p_{S}\frac{w_{i}(S)}{|S|}\Bigg{)}.$

As the adversary selects the worst layer, we aim to find $p\in\Delta(2^{V})$
that minimizes the maximum of (3) among $i\in[k]$.

Here we explain how to select an appropriate metric. The density and regret
metrics are useful when we are concerned with multilayer networks with
homogeneous layers such as time-dependent follower-followee relations in the
Twitter network. Although the density metric can be the first choice, the
regret metric is more suitable for robust analysis. For example, consider the
case where there are a number of layers consistent with each other together
with some noisy (e.g., random) layers. The density metric would suffer from
the effect of the noisy layers, but the regret metric would avoid it and find
dense subgraphs in the other meaningful layers. On the other hand, the robust
ratio metric is useful when we analyze multilayer networks with heterogeneous
layers such as brain networks with structural and functional connectivity
layers. From its definition, the robust ratio metric would find subgraphs that
are reasonably dense for all layers. The density and regret metrics focus only
on the layers with small optimal densities and the layers with large optimal
densities, respectively.

Unified concept: $(\bm{\alpha},\bm{\beta})$-density. Here we introduce a
general metric, enabling us to deal with the above three metrics in a unified
manner. An important fact is that the robust ratio and regret metrics can be
obtained by affine transformations of the density. Specifically, when we
select a subgraph according to a probability distribution $p\in\Delta(2^{V})$
and the adversary selects a layer $i\in[k]$, we define the
_$(\bm{\alpha},\bm{\beta})$ -density_ using two vectors
$\bm{\alpha}\in\mathbb{R}^{k}_{+}$ and $\bm{\beta}\in\mathbb{R}^{k}$ as
follows:

(4) $\displaystyle\mathbb{E}_{S\sim
p}\left[\alpha_{i}\frac{w_{i}(S)}{|S|}+\beta_{i}\right]\quad\Bigg{(}=\alpha_{i}\sum_{S\subseteq
V}p_{S}\frac{w_{i}(S)}{|S|}+\beta_{i}\Bigg{)}.$

Note that the above three metrics, the density, robust ratio, and regret, are
equivalent to the $(\bm{1},\bm{0})$-density,
$((|S_{i}^{*}|/w_{i}(S_{i}^{*}))_{i\in[k]},\bm{0})$-density, and
$(\bm{1},(-w_{i}(S_{i}^{*})/|S_{i}^{*}|)_{i\in[k]})$-density111The actual
regret value is the negation of
$(\bm{1},(-w_{i}(S_{i}^{*})/|S_{i}^{*}|)_{i\in[k]})$-density., respectively.
Note that $w_{i}(S_{i}^{*})/|S_{i}^{*}|$ is polynomially computable for each
$i\in[k]$ (Charikar, 2000; Goldberg, 1984). We refer to
$(\bm{\alpha},\bm{\beta})$-DENSITY as the problem of finding
$p\in\Delta(2^{V})$ that maximizes the minimum of the
$(\bm{\alpha},\bm{\beta})$-density among $i\in[k]$. Therefore, in the
following, we aim to design an algorithm for
$(\bm{\alpha},\bm{\beta})$-DENSITY.

##  4\. Algorithm

In this section, we provide an LP-based polynomial-time exact algorithm for
$(\bm{\alpha},\bm{\beta})$-DENSITY. Let $((x_{e})_{e\in E},(y_{v})_{v\in
V},t)$ be continuous variables. We consider the following LP:

(10) $\displaystyle\begin{array}[]{rll}\text{max.}&t&\\\\[3.0pt]
\text{s.t.}&\displaystyle t\leq\alpha_{i}\cdot\sum_{e\in
E_{i}}w_{i}(e)x_{e}+\beta_{i}&(\forall i\in[k]),\\\\[3.0pt] &x_{e}\leq y_{u},\
x_{e}\leq y_{v}&(\forall e=\\{u,v\\}\in E),\\\\[3.0pt]
&\displaystyle\sum_{v\in V}y_{v}=1,&\\\\[3.0pt] &x_{e},y_{v}\geq 0&(\forall
e\in E,\forall v\in V).\end{array}$

When $(\bm{\alpha},\bm{\beta})=(\bm{1},\bm{0})$, this formulation coincides
with the LP introduced by Jethava and Beerenwinkel (Jethava and Beerenwinkel,
2015) for the densest common subgraph problem. However, they did not point out
the connection between a solution of the LP and a distribution in
$\Delta(2^{V})$.

Our algorithm first computes an optimal solution to LP (10), denoted by
$((\hat{x}_{e})_{e\in E},(\hat{y}_{v})_{v\in V},\hat{t})$. Then we set
$r_{0},r_{1},\dots,r_{\ell}$ to be the reals such that
$\\{r_{0},r_{1},\dots,r_{\ell}\\}=\\{\hat{y}_{v}\mid v\in V\\}\cup\\{0\\}$ and
$r_{0}~{}(=0)<r_{1}<r_{2}<\cdots<r_{\ell}$. In addition, let $S_{j}=\\{v\in
V\mid\hat{y}_{v}\geq r_{j}\\}$ $(j=1,\dots,\ell)$. The output of our algorithm
is the probability distribution $\hat{p}\in\Delta(2^{V})$ defined as

$\displaystyle\hat{p}_{S_{j}}=(r_{j}-r_{j-1})\cdot|S_{j}|\ (j\in[\ell])\
\text{ and }$ (11) $\displaystyle\hat{p}_{S}=0\text{ for the other $S$'s}.$

We remark that $\sum_{S\in 2^{V}}\hat{p}_{S}=1$ holds because

$1=\sum_{v\in
V}\hat{y}_{v}=\sum_{j\in[\ell]}\sum_{v:\,\hat{y}_{v}=r_{j}}r_{j}=\sum_{j\in[\ell]}r_{j}(|S_{j}|-|S_{j+1}|)=\sum_{j\in[\ell]}\hat{p}_{S_{j}},$

where $S_{\ell+1}=\emptyset$ for ease of notation.

Our algorithm is formally described in Algorithm 1. Clearly, the algorithm
runs in polynomial time.

1 Solve LP (10) to obtain an optimal solution $((\hat{x}_{e})_{e\in
E},(\hat{y}_{v})_{v\in V},\hat{t})$;

2 Let $r_{0},r_{1},\dots,r_{\ell}$ be reals such that
$\\{r_{0},r_{1},\dots,r_{\ell}\\}=\\{\hat{y}_{v}\mid v\in V\\}\cup\\{0\\}$ and
$r_{0}~{}(=0)<r_{1}<r_{2}<\cdots<r_{\ell}$;

3 Let $S_{j}=\\{v\in V\mid\hat{y}_{v}\geq r_{j}\\}$ $(j=1,\dots,\ell)$;

return $\hat{p}\in\Delta(2^{V})$ with
$\hat{p}_{S_{j}}=(r_{j}-r_{j-1})\cdot|S_{j}|$ ($j=1,\dots,\ell$) (and
$\hat{p}_{S}=0$ for the other $S$’s).

Algorithm 1 LP-based algorithm

##  5\. Analysis

In this section, we first demonstrate that the output $\hat{p}$ of Algorithm 1
is optimal to $(\bm{\alpha},\bm{\beta})$-DENSITY. Then we analyze the
structure of the output $\hat{p}$ with a special attention to its support
size.

###  5.1. Optimality of the output of Algorithm 1

We first prove that the $(\bm{\alpha},\bm{\beta})$-density of $\hat{p}$ is
equal to the optimal value of LP (10):

###### Lemma 5.1.

It holds that

$\displaystyle\min_{i\in[k]}\mathbb{E}_{S\sim\hat{p}}\Bigl{[}\alpha_{i}\frac{w_{i}(S)}{|S|}+\beta_{i}\Bigr{]}=\min_{i\in[k]}\Bigl{[}\alpha_{i}\sum_{e\in
E_{i}}w_{i}(e)\hat{x}_{e}+\beta_{i}\Bigr{]}.$

###### Proof.

To prove the lemma, we show that for any $i\in[k]$,

$\displaystyle\mathbb{E}_{S\sim\hat{p}}\Bigl{[}\alpha_{i}\frac{w_{i}(S)}{|S|}+\beta_{i}\Bigr{]}=\alpha_{i}\sum_{e\in
E_{i}}w_{i}(e)\hat{x}_{e}+\beta_{i}.$

By the definition of $S_{j}$ $(j=1,\ldots,\ell)$, for each $e=\\{u,v\\}\in
E_{i}$, we observe that $e\in E_{i}[S_{j}]\iff u,v\in S_{j}\iff\hat{y}_{u}\geq
r_{j}\text{ and }\hat{y}_{v}\geq r_{j}\iff\hat{x}_{e}\geq r_{j}$, where the
last equivalence follows from $\hat{x}_{e}=\min\\{\hat{y}_{u},\hat{y}_{v}\\}$.

Recall that $\hat{p}$ is defined as (4). Then we have

$\displaystyle\mathbb{E}_{S\sim\hat{p}}\Bigl{[}\alpha_{i}\frac{w_{i}(S)}{|S|}+\beta_{i}\Bigr{]}$
$\displaystyle=\alpha_{i}\sum_{j\in[\ell]}\frac{w_{i}(S_{j})}{|S_{j}|}\cdot\hat{p}_{S_{j}}+\beta_{i}$
$\displaystyle=\alpha_{i}\sum_{j\in[\ell]}w_{i}(S_{j})\cdot(r_{j}-r_{j-1})+\beta_{i}$
$\displaystyle=\alpha_{i}\sum_{j\in[\ell]}\sum_{e\in
E_{i}[S_{j}]}w_{i}(e)\cdot(r_{j}-r_{j-1})+\beta_{i}$
$\displaystyle=\alpha_{i}\sum_{j\in[\ell]}\sum_{e\in E_{i}:\,\hat{x}_{e}\geq
r_{j}}w_{i}(e)\cdot(r_{j}-r_{j-1})+\beta_{i}$
$\displaystyle=\alpha_{i}\sum_{e\in E_{i}}w_{i}(e)\hat{x}_{e}+\beta_{i},$

where the second last equality follows from the above equivalence and the last
equality follows from the fact that for each $e\in E$,
$\hat{x}_{e}=r_{j^{\prime}}=\sum_{j=1}^{j^{\prime}}(r_{j}-r_{j-1})$ for some
$j^{\prime}$. ∎

Next, we prove that the optimal value of LP (10) gives an upper bound on the
optimal value of $(\bm{\alpha},\bm{\beta})$-DENSITY:

###### Lemma 5.2.

It holds that

$\displaystyle\min_{i\in[k]}\Bigl{[}\alpha_{i}\sum_{e\in
E_{i}}w_{i}(e)\hat{x}_{e}+\beta_{i}\Bigr{]}\geq\max_{p\in\Delta(2^{V})}\min_{i\in[k]}\mathbb{E}_{S\sim
p}\Bigl{[}\alpha_{i}\frac{w_{i}(S)}{|S|}+\beta_{i}\Bigr{]}.$

###### Proof.

Let us take an arbitrary $\tilde{p}\in\Delta(2^{V})$. We consider a solution
$((\tilde{x}_{e})_{e\in E},(\tilde{y}_{v})_{v\in V},\tilde{t})$ of LP (10)
such that

(12)
$\displaystyle\tilde{x}_{e}=\\!\\!\\!\\!\sum_{\begin{subarray}{c}S\subseteq
V:\\\ e\in E[S]\end{subarray}}\frac{\tilde{p}_{S}}{|S|},\
\tilde{y}_{v}=\\!\\!\\!\\!\sum_{\begin{subarray}{c}S\subseteq V:\\\ v\in
S\end{subarray}}\frac{\tilde{p}_{S}}{|S|},\
\tilde{t}=\\!\min_{i\in[k]}\Bigl{[}\alpha_{i}\\!\sum_{e\in
E_{i}}w_{i}(e)\tilde{x}_{e}+\beta_{i}\Bigr{]}.$

The solution $((\tilde{x}_{e})_{e\in E},(\tilde{y}_{v})_{v\in V},\tilde{t})$
is feasible for LP (10); in fact, the only concern is the third constraint but
we see that

(13) $\displaystyle\sum_{v\in V}\tilde{y}_{v}=\sum_{v\in
V}\sum_{\begin{subarray}{c}S\subseteq V:\\\ v\in
S\end{subarray}}\frac{\tilde{p}_{S}}{|S|}=\sum_{S\subseteq
V}\left(|S|\cdot\frac{\tilde{p}_{S}}{|S|}\right)=1.$

Hence, the optimal value of LP (10) is at least $\tilde{t}$. Moreover, we have

$\displaystyle\tilde{t}$
$\displaystyle=\min_{i\in[k]}\Bigl{[}\alpha_{i}\sum_{e\in
E_{i}}w_{i}(e)\tilde{x}_{e}+\beta_{i}\Bigr{]}$
$\displaystyle=\min_{i\in[k]}\Bigl{[}\alpha_{i}\sum_{S\subseteq
V}\tilde{p}_{S}\cdot\frac{w_{i}(S)}{|S|}+\beta_{i}\Bigr{]}=\min_{i\in[k]}\mathbb{E}_{S\sim\tilde{p}}\Bigl{[}\alpha_{i}\frac{w_{i}(S)}{|S|}+\beta_{i}\Bigr{]}.$

Recalling that $\tilde{p}$ is taken arbitrarily from $\Delta(2^{V})$, we see
that the optimal value of LP (10) is at least that of
$(\bm{\alpha},\bm{\beta})$-DENSITY. ∎

Combining Lemmas 5.1 and 5.2, we have the desired result:

###### Theorem 5.3.

Algorithm 1 outputs an optimal solution to $(\bm{\alpha},\bm{\beta})$-DENSITY.

###  5.2. Hierarchical structure and support size of the output of Algorithm 1

Here we observe some useful properties of the output of Algorithm 1. By the
design of Algorithm 1, we see that the support of the output $\hat{p}$ of the
algorithm (i.e., an optimal solution to $(\bm{\alpha},\bm{\beta})$-DENSITY)
has a hierarchical structure. We denote the support of $p\in\Delta(2^{V})$ by
$\operatorname{supp}(p)=\\{S\subseteq V\mid p_{S}>0\\}$.

###### Proposition 5.4.

Algorithm 1 outputs a hierarchical solution $\hat{p}$, i.e., $S\subseteq T$ or
$S\supseteq T$ for any $S,T\in\operatorname{supp}(\hat{p})$.

From this proposition, the output $\hat{p}$ has support size at most $|V|-1$
(since the empty set and the singletons are useless). In addition, if we pick
an optimal basic solution to LP (10) in Algorithm 1, the support size of
$\hat{p}$ becomes at most $k$. For the definition of a basic solution, see
e.g., Vanderbei’s book (Vanderbei, 2020).

###### Theorem 5.5.

If Algorithm 1 takes a basic optimal solution to LP (10), its output $\hat{p}$
has support size at most $k$.

###### Proof.

Let $((\hat{x}_{e})_{e\in E},(\hat{y}_{v})_{v\in V},\hat{t})$ be a basic
solution. Without loss of generality, we may assume that

(14) $\displaystyle\hat{x}_{e}=\min\\{\hat{y}_{u},\hat{y}_{v}\\}\quad(\forall
e=\\{u,v\\}\in E).$

Recall that $\ell$ denotes the number of different positive values in
$\\{\hat{y}_{v}\mid\hat{y}_{v}>0,\ v\in V\\}$. Let $V_{0}=\\{v\in
V\mid\hat{y}_{v}=0\\}$. We divide $V\setminus V_{0}$ into $\ell$ subsets of
vertices sharing the same value of $\hat{y}_{v}$, denoted by
$V_{1},\ldots,V_{\ell}$.

Let us focus on the constraints in LP (10) that are satisfied with equality.
For each $j=0,1,\ldots,\ell$, let $F_{j}\subseteq E[V_{j}]$ be a spanning
forest in $E[V_{j}]$. Let $\rho$ be the number of connected components in
$E[V_{1}]\cup\cdots\cup E[V_{\ell}]$, and let $\zeta$ be that of $E[V_{0}]$.
We arbitrarily take $\zeta$ vertices, denoted by $u_{1},\ldots,u_{\zeta}$, one
from each connected component in $E[V_{0}]$. We focus on the following
constraints (satisfied with equality):

$\displaystyle\textstyle t=\alpha_{i}\sum_{e\in
E_{i}}w_{i}(e)x_{e}+\beta_{i}\quad(\forall i\in K^{\prime}),$
$\displaystyle\textstyle x_{e}=y_{u},\ x_{e}=y_{v}\quad(\forall
e=\\{u,v\\}\in\bigcup_{j=0}^{\ell}E[V_{j}]),$ $\displaystyle\textstyle
x_{e}=y_{v}\quad(\forall e=\\{u,v\\}\in
E\setminus\bigcup_{j=0}^{\ell}E[V_{j}],\ \hat{y}_{v}<\hat{y}_{u}),$
$\displaystyle\textstyle\sum_{v\in V}y_{v}=1,$ $\displaystyle\textstyle
y_{v}=0\quad(\forall v\in V_{0}),$

where $K^{\prime}=\\{i\in[k]\mid t=\alpha_{i}\cdot\sum_{e\in
E_{i}}w_{i}(e)x_{e}+\beta_{i}\\}$.

We prove that the coefficient matrix of those constraints is not full-rank.
For ease of discussion, we remove constraints that are represented by a linear
combination of others. Specifically, it is enough to focus on the following
constraints:

(15) $\displaystyle\textstyle t=\alpha_{i}\cdot\sum_{e\in
E_{i}}w_{i}(e)x_{e}+\beta_{i}\quad(\forall i\in K^{\prime}),$ (16)
$\displaystyle\textstyle x_{e}=y_{u},\ x_{e}=y_{v}\quad(\forall
e=\\{u,v\\}\in\bigcup_{j=0}^{\ell}F_{j}),$ (17) $\displaystyle\textstyle
x_{e}=y_{v}\quad(\forall e=\\{u,v\\}\in\bigcup_{j=0}^{\ell}(E[V_{j}]\setminus
F_{j}),\ u<v),$ (18) $\displaystyle\textstyle x_{e}=y_{v}\quad(\forall
e=\\{u,v\\}\in E\setminus\bigcup_{j=0}^{\ell}E[V_{j}],\
\hat{y}_{v}<\hat{y}_{u}),$ (19) $\displaystyle\textstyle\sum_{v\in
V}y_{v}=1,\quad$ (20) $\displaystyle\textstyle y_{u_{i}}=0\quad(\forall
i\in[\zeta]).$

There are two types of missing constraints. First, let $e=\\{u,v\\}\in
E[V_{j}]\setminus F_{j}$ $(j\in[\ell])$ such that $x_{e}=y_{v}$ appears in
(17). There exists a cycle $C$ in $F_{j}\cup\\{e\\}$. By a telescoping sum of
the constraints (16) along $C$, i.e., $y_{v}=x_{e^{\prime}}$,
$x_{e^{\prime}}=y_{v^{\prime}}$, …, $x_{e^{\prime\prime}}=y_{u}$, we obtain
$y_{u}=y_{v}$. Then, constraint $x_{e}=y_{u}$ is obtained by summing
$x_{e}=y_{v}$ for $y_{u}=y_{v}$. Next, let $v\in V_{0}$ belong to a connected
component containing $u_{j}$ $(j\in[\zeta])$. By summing (16) along a path
from $v$ to $u_{j}$, i.e., $y_{v}=x_{e^{\prime}}$,
$x_{e^{\prime}}=y_{v^{\prime}}$, …, $x_{e^{\prime\prime}}=y_{u_{j}}$, we
obtain $y_{v}=y_{u_{j}}$. Then, constraint $y_{v}=0$ is obtained by summing
$y_{v}=y_{u_{j}}$ and $y_{u_{j}}=0$ from (20).

The rank of the coefficient matrix is equal to the number of constraints
(15)–(20). For (15), we have at most $k$ constraints; the number of
constraints (16) is $2(|V|-\rho-\zeta)$; that for (17) and (18) is
$|E|-|V|+\rho+\zeta$; that for (19) and (20) is $1+\zeta$. Therefore, we have
at most $|V|+|E|-\rho+1+k$ constraints.

We have $|V|+|E|+1$ variables in LP (10). If $\rho>k$, then we have at most
$|V|+|E|$ constraints, and hence the coefficient matrix of (15)–(20) cannot
have rank $|V|+|E|+1$. As the solution is basic, $\rho\leq k$ must hold.
Recall that each $V_{j}$ $(j=1,\ldots,\ell)$ has at least one connected
component. Therefore, we have at most $k$ different positive values of
$y_{v}$’s, which implies that the output $\hat{p}$ has support size at most
$k$. ∎

##  6\. Preprocessing

In this section, we present a simple, scalable preprocessing algorithm, which
often reduces the size of the input networks significantly and results in a
substantial speed-up. Specifically, the algorithm first computes an
approximate solution by solving an LP, which is much smaller than LP (10) in
practice, and then removes vertices from the original network using the
information of the approximate solution obtained. We assume that $S^{*}_{i}$
for all $i\in[k]$ are known in advance because they can be computed
efficiently using Charikar’s LP-based algorithm (Charikar, 2000) together with
the preprocessing algorithm introduced by Balalau et al. (Balalau et al.,
2015). To describe our algorithm, we introduce some notations. For $S\subseteq
V$, $v\in S$, and $i\in[k]$, let $d_{i}(S,v)$ denote the weighted degree of
$v$ in the subgraph induced by $S$ in layer $i$, i.e.,
$d_{i}(S,v)\coloneqq\sum_{e\in E_{i}[S]:\,v\in e}w_{i}(e)$. When $S=V$, we
simply write $d_{i}(v)$.

We first describe a fast algorithm for finding an approximate solution for
$(\bm{\alpha},\bm{\beta})$-DENSITY. Specifically, we compute a probability
distribution $q\in\Delta(2^{V})$ that maximizes the
$(\bm{\alpha},\bm{\beta})$-density under the constraint that $q_{S}=0$ for all
$S\in 2^{V}\setminus\\{S_{1}^{*},\dots,S_{k}^{*}\\}$. The distribution can be
found by solving the following LP:

(25) $\displaystyle\begin{array}[]{rll}\text{max.}&t&\\\
\text{s.t.}&\displaystyle
t\leq\alpha_{i}\sum_{j\in[k]}\frac{w_{i}(S_{j}^{*})}{|S_{j}^{*}|}q_{j}+\beta_{i}&(\forall
i\in[k]),\\\ &\displaystyle\sum_{j\in[k]}q_{j}=1,&\\\ &q_{j}\geq 0&(\forall
j\in[k]).\end{array}$

Note that this LP has $k+1$ variables and $2k+1$ constraints. As $k$ is
usually much smaller than $|V|$ and $|E|$, this LP is much smaller than LP
(10) in practice.

Next we describe an algorithm for removing vertices using the information of
the above approximate solution $q$. Let $\ell^{*}$ be the
$(\bm{\alpha},\bm{\beta})$-density of $q$, i.e., the optimal value of LP (25).
Note that this is a lower bound on the optimal value of
$(\bm{\alpha},\bm{\beta})$-DENSITY. Our algorithm iteratively removes any
vertex $v^{*}$ that satisfies $\max_{i\in[k]}[\alpha_{i}\cdot
d_{i}(V^{\prime},v^{*})+\beta_{i}]<\ell^{*}$, where $V^{\prime}$ is a
remaining vertex set (initially $V^{\prime}=V$), as long as there exists such
a vertex. For reference, we describe the procedure in Algorithm 2. This
algorithm can be implemented to run in $O(k|E|+|V|\log|V|)$ time.

Input : $(V,(E_{i})_{i\in[k]})$ with $w_{1},\dots,w_{k}$, and
$\ell^{*}\in\mathbb{R}$

Output : $(V^{\prime},(E_{i}[V^{\prime}])_{i\in[k]})$

1 $V^{\prime}\leftarrow V$;

2 while  _True_ do

3 Let $v^{*}\in\operatorname*{arg\,min}_{v\in
V^{\prime}}\max_{i\in[k]}[\alpha_{i}\cdot d_{i}(V^{\prime},v)+\beta_{i}]$;

4 if  _$\max_{i\in[k]}[\alpha_{i}\cdot
d_{i}(V^{\prime},v^{*})+\beta_{i}]\geq\ell^{*}$_ then

5 return $(V^{\prime},(E_{i}[V^{\prime}])_{i\in[k]})$.

6 else $V^{\prime}\leftarrow V^{\prime}\setminus\\{v^{*}\\}$;

7

Algorithm 2 Remove useless vertices

From now on, we demonstrate that the algorithm does not remove any vertex that
is contained in a subset in $\operatorname{supp}(p)$, where $p$ is an
arbitrary optimal solution to $(\bm{\alpha},\bm{\beta})$-DENSITY. The
following is a key lemma in our analysis.

###### Lemma 6.1.

Let $\ell^{*}$ be a lower bound on the optimal value of
$(\bm{\alpha},\bm{\beta})$-DENSITY. If $\max_{i\in[k]}[\alpha_{i}\cdot
d_{i}(v^{*})+\beta_{i}]<\ell^{*}$, then $y_{v^{*}}=0$ for any optimal solution
to LP (10).

###### Proof.

We prove the lemma by contradiction. We denote by $((\hat{x}_{e})_{e\in
E},(\hat{y}_{v})_{v\in V},\hat{t})$ an optimal solution to LP (10) and let
$v^{*}\in V$ be a vertex that satisfies $\alpha_{i}\cdot
d_{i}(v^{*})+\beta_{i}<\ell^{*}$ for all $i\in[k]$. Suppose for contradiction
that $\hat{y}_{v^{*}}>0$.

We construct a solution $((x_{e})_{e\in E},(y_{v})_{v\in V},t)$ of LP (10) as
follows:

$\displaystyle x_{e}$
$\displaystyle=\begin{cases}\frac{1}{1-\hat{y}_{v^{*}}}\cdot\hat{x}_{e}&(e\not\ni
v^{*}),\\\ 0&(e\ni v^{*}),\end{cases}\quad
y_{v}=\begin{cases}\frac{1}{1-\hat{y}_{v^{*}}}\cdot\hat{y}_{v}&(v\neq
v^{*}),\\\ 0&(v=v^{*}),\end{cases}$ $\displaystyle t$
$\displaystyle=\min_{i\in[k]}\Bigl{[}\alpha_{i}\sum_{e\in
E_{i}}w_{i}(e)x_{e}+\beta_{i}\Bigr{]}.$

It is easy to see that $((x_{e})_{e\in E},(y_{v})_{v\in V},t)$ is a feasible
solution of LP (10). Moreover, we have

$\displaystyle t$ $\displaystyle=\min_{i\in[k]}\Bigl{[}\alpha_{i}\sum_{e\in
E_{i}}w_{i}(e)x_{e}+\beta_{i}\Bigr{]}$
$\displaystyle=\min_{i\in[k]}\Bigl{[}\alpha_{i}\frac{\sum_{e\in
E_{i}:\,v^{*}\not\in
e}w_{i}(e)\hat{x}_{e}}{1-\hat{y}_{v^{*}}}+\beta_{i}\Bigr{]}$
$\displaystyle=\min_{i\in[k]}\Bigl{[}\alpha_{i}\frac{\sum_{e\in
E_{i}}w_{i}(e)\hat{x}_{e}-\sum_{e\in E_{i}:\,v^{*}\in
e}w_{i}(e)\hat{x}_{e}}{1-\hat{y}_{v^{*}}}+\beta_{i}\Bigr{]}$
$\displaystyle\geq\min_{i\in[k]}\Bigl{[}\alpha_{i}\frac{\sum_{e\in
E_{i}}w_{i}(e)\hat{x}_{e}-\sum_{e\in E_{i}:\,v^{*}\in
e}w_{i}(e)\hat{y}_{v^{*}}}{1-\hat{y}_{v^{*}}}+\beta_{i}\Bigr{]}$
$\displaystyle=\min_{i\in[k]}\Bigl{[}\alpha_{i}\frac{\sum_{e\in
E_{i}}w_{i}(e)\hat{x}_{e}-d_{i}(v^{*})\hat{y}_{v^{*}}}{1-\hat{y}_{v^{*}}}+\beta_{i}\Bigr{]}$
$\displaystyle=\min_{i\in[k]}\frac{(\alpha_{i}\sum_{e\in
E_{i}}w_{i}(e)\hat{x}_{e}+\beta_{i})-(\alpha_{i}\cdot
d_{i}(v^{*})+\beta_{i})\cdot\hat{y}_{v^{*}}}{1-\hat{y}_{v^{*}}}$
$\displaystyle>\min_{i\in[k]}\frac{(\alpha_{i}\sum_{e\in
E_{i}}w_{i}(e)\hat{x}_{e}+\beta_{i})-\ell^{*}\cdot\hat{y}_{v^{*}}}{1-\hat{y}_{v^{*}}}$
$\displaystyle=\frac{\hat{t}-\ell^{*}\cdot\hat{y}_{v^{*}}}{1-\hat{y}_{v^{*}}}\geq\frac{\hat{t}-\hat{t}\cdot\hat{y}_{v^{*}}}{1-\hat{y}_{v^{*}}}=\hat{t},$

where the first inequality follows from $\hat{x}_{e}\leq\hat{y}_{v^{*}}$ for
each $e\ni v^{*}$, the second inequality follows from the assumptions
$\alpha_{i}\cdot d_{i}(v^{*})+\beta_{i}<\ell^{*}$ and $\hat{y}_{v^{*}}>0$, and
the third inequality follows from Lemma 5.2. This contradicts the optimality
of $((\hat{x}_{e})_{e\in E},(\hat{y}_{v})_{v\in V},\hat{t})$. ∎

###### Theorem 6.2.

Let $\ell^{*}$ be a lower bound on the optimal value of
$(\bm{\alpha},\bm{\beta})$-DENSITY. Then, any vertex $v^{*}$ that satisfies
$\max_{i\in[k]}[\alpha_{i}\cdot d_{i}(v^{*})+\beta_{i}]<\ell^{*}$ is not
contained in any subset in the support of any optimal solution to
$(\bm{\alpha},\bm{\beta})$-DENSITY.

###### Proof.

Let $v^{*}$ be any vertex with $\max_{i\in[k]}[\alpha_{i}\cdot
d_{i}(v^{*})+\beta_{i}]<\ell^{*}$. Let $\hat{p}$ be any optimal solution to
$(\bm{\alpha},\bm{\beta})$-DENSITY. Construct $((\hat{x}_{e})_{e\in
E},(\hat{y}_{v})_{v\in V},\hat{t})$ from $\hat{p}$ as in (12). From Lemma 5.1
and the proof of Lemma 5.2, we see that $((\hat{x}_{e})_{e\in
E},(\hat{y}_{v})_{v\in V},\hat{t})$ is an optimal solution to LP (10). Thus,
by Lemma 6.1, we have $\hat{y}_{v^{*}}=0$. By the construction of
$\hat{y}_{v^{*}}$ in (12), $\hat{p}_{S}=0$ for all $S\subseteq V$ containing
$v^{*}$. ∎

This theorem indicates that Algorithm 2 does not remove any vertex that is
contained in a subset in the support of any optimal solution to
$(\bm{\alpha},\bm{\beta})$-DENSITY.

##  7\. Experimental Evaluation

In this section, we conduct computational experiments using synthetic graphs
and real-world networks to verify the validity of our proposed model and to
evaluate the performance of our proposed algorithms. All experiments were
conducted on a machine equipped with Intel Xeon W 10-core processor and 64GB
RAM. Algorithms were implemented in Python using Gurobi Optimizer 9.0.2.

###  7.1. Validity of our model

Here we aim to verify the validity of our model using synthetic graphs. To
this end, we use a randomly generated multilayer network with a planted
clique, and examine whether an optimal solution detects vertex subsets close
to the clique.

We first explain our random procedure for generating multilayer networks. We
produce an (unweighted) random power-law graph as a layer using the Chung–Lu
model (Chung and Lu, 2002), where we first specify an expected degree $d_{v}$
for each $v\in V$ according to the power-law distribution with exponent
$\beta$, and then connect each pair of vertices $\\{u,v\\}$ with probability
$\frac{d_{u}\cdot d_{v}}{\sum_{r\in V}d_{r}}$. Note that the graph becomes
sparser as the exponent $\beta$ increases. In this multilayer network, we
randomly select a vertex subset $V_{c}$ with some size, and plant a clique on
$V_{c}$ (in some specified layers).

To evaluate the performance of optimal solution $p\in\Delta(2^{V})$ to
$(\bm{\alpha},\bm{\beta})$-DENSITY in the above multilayer network, we
introduce the following measure, which we refer to as the (expected) F
measure:

$\displaystyle\text{F
measure}\coloneqq\frac{2\cdot\text{precision}\cdot\text{recall}}{\text{precision}+\text{recall}},$

where $\text{precision}\coloneqq\mathbb{E}_{S\sim p}\left[\frac{|S\cap
V_{c}|}{|S|}\right]=\sum_{S\subseteq V}p_{S}\frac{|S\cap V_{c}|}{|S|}$ and
$\text{recall}\coloneqq\mathbb{E}_{S\sim p}\left[\frac{|S\cap
V_{c}|}{|V_{c}|}\right]=\sum_{S\subseteq V}p_{S}\frac{|S\cap
V_{c}|}{|V_{c}|}$. The F measure approaches to 1 if $p$ tends to be close to
$V_{c}$.

We first investigate the case where a clique is planted in all layers. We
generate $k~{}(=1,2,3,4,5)$ power-law graphs (i.e., layers) with $\beta=2.3$
on $V$ with $|V|=\text{1,000}$. Then we randomly select a subset
$V_{c}\subseteq V$ consisting of $10$ vertices, and plant a clique on $V_{c}$
in all layers. The performance of optimal solutions to
$(\bm{\alpha},\bm{\beta})$-DENSITY is shown in Figure 2. As can be seen, for
any metric, the F measure is reasonably large for $k\geq 2$, meaning that our
algorithm tends to detect $V_{c}$ using the information of multiple layers.

(a) Clique in all layers.

(b) Clique in only one layer.

Figure 2. Performance of optimal solutions. Each point corresponds to the
average value over 100 network realizations.

(a) Clique in all layers.

(b) Clique in only one layer.

Figure 3. Performance of single vertex subsets obtained by optimal solutions.
The same averaging procedure is applied. Table 1. Performance of our algorithm
for real-world datasets. OPT and $|\text{supp}|$ denote the optimal value of
$(\bm{\alpha},\bm{\beta})$-DENSITY and the support size of the optimal
solution, respectively. Note that the objective values for the regret are
negated. LB indicates the optimal value of LP (25). $|V^{\prime}|$ and
$|E[V^{\prime}]|$, respectively, denote the number of vertices and edges after
running Algorithm 2. Preprocess time represents the running time of Algorithm
2. Total time represents the running time of our proposed algorithm, i.e.,
Algorithm 1 together with Algorithm 2. For reference, the running time of
Algorithm 1 without Algorithm 2 is also presented in the next column. The last
three columns report the performance of single vertex subsets obtained by
optimal solutions. For each instance and metric, the best value among the
algorithms is in bold.

Dataset | $|V|$ | $|E|$ | $k$ | Metric | OPT | $|\text{supp}|$ | LB | $|V^{\prime}|$ | $|E[V^{\prime}]|$ |  Preprocess time (s) |  Total time (s) |  w/o preprocess time (s) | Density |  Robust ratio | Regret
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  | Density | 1.1950 | 3 | 1.1313 | 123 | 3,126 | 0.04 | 0.66 | 0.73 | 1.1875 | 0.6579 | -0.6282
WILDBIRDS | 202 | 4,574 | 6 | Robust ratio | 0.7707 | 4 | 0.7364 | 121 | 2,901 | 0.05 | 0.77 | 0.79 | 1.0058 | 0.6970 | -0.5065
(Firth and Sheldon, 2015)$*$ |  |  |  | Regret | -0.4122 | 2 | -0.4473 | 124 | 2,976 | 0.04 | 0.71 | 0.77 | 0.9129 | 0.6518 | -0.5821
|  |  |  | Density | 4.7023 | 2 | 4.4257 | 730 | 6,403 | 1208.79 | 1304.28 | 1834.88 | 4.6316 | 0.6399 | -3.1581
AS-733 | 7,716 | 24,179 | 733 | Robust ratio | 0.7721 | 4 | 0.7295 | 278 | 3,452 | 1165.87 | 1211.82 | 1906.21 | 4.5789 | 0.7057 | -1.9621
(Leskovec et al., 2005)$\dagger$ |  |  |  | Regret | -6.9861 | 4 | -1.8023 | 228 | 3,020 | 1289.15 | 1329.56 | 1950.33 | 3.6154 | 0.7317 | -1.6316
|  |  |  | Density | 12.0263 | 1 | 12.0263 | 118 | 1,815 | 0.29 | 0.92 | 8.31 | 12.0263 | 0.9666 | -0.4382
Oregon1 | 11,492 | 26,461 | 9 | Robust ratio | 0.9808 | 3 | 0.9738 | 110 | 1,702 | 0.48 | 1.02 | 7.55 | 11.9016 | 0.9728 | -0.3578
(Leskovec et al., 2005)$\dagger$ |  |  |  | Regret | -0.2544 | 3 | -0.3441 | 110 | 1,702 | 0.50 | 0.95 | 9.72 | 11.8889 | 0.9728 | -0.3578
|  |  |  | Density | 22.2503 | 2 | 21.8228 | 1,130 | 14,166 | 1.90 | 4.31 | 93.30 | 18.4000 | 0.1183 | -247.6000
MoscowAthletics2013 | 88,804 | 186,846 | 3 | Robust ratio | 0.4709 | 3 | 0.3566 | 105 | 787 | 2.02 | 2.12 | 112.50 | 13.2000 | 0.4197 | -181.1667
(Omodei et al., 2015)$\ddagger$ |  |  |  | Regret | -125.4477 | 2 | -130.1752 | 88,804 | 186,846 | 1.40 | 336.56 | 334.43 | 0.4286 | 0.0186 | -200.8333
|  |  |  | Density | 3.5649 | 2 | 3.5077 | 17,551 | 184,659 | 1.98 | 332.02 | 5466.73 | 3.5625 | 0.1768 | -17.0770
NYClimateMarch2014 | 102,439 | 329,474 | 3 | Robust ratio | 0.6661 | 3 | 0.5503 | 3,901 | 78,126 | 2.36 | 147.40 | 6058.12 | 2.3839 | 0.6652 | -6.8047
(Omodei et al., 2015)$\ddagger$ |  |  |  | Regret | -2.3052 | 3 | -3.5835 | 102,439 | 329,473 | 1.58 | 18575.97 | 18909.39 | 1.2785 | 0.3576 | -2.3221
|  |  |  | Density | 60.7462 | 2 | 60.7462 | 375 | 4,536 | 11.07 | 11.72 | 4268.72 | 52.7143 | 0.0840 | -837.8095
Cannes2013 | 438,537 | 848,017 | 3 | Robust ratio | 0.3633 | 3 | 0.3633 | 246 | 1,617 | 11.53 | 11.72 | 3981.18 | 42.6000 | 0.2441 | -365.0667
(Omodei et al., 2015)$\ddagger$ |  |  |  | Regret | -132.8628 | 2 | -132.8628 | 438,537 | 848,017 | 7.31 | 4021.19 | 4022.75 | 0.6667 | 0.0073 | -155.5000
|  |  |  | Density | 1.1891 | 10 | 1.1236 | 191,074 | 559,628 | 18.52 | 7826.41 | 13787.37 | 0.4146 | 0.0387 | -31.2295
DBLP | 513,627 | 888,353 | 10 | Robust ratio | 0.1178 | 10 | 0.1125 | 317,231 | 687,335 | 17.01 | 15323.66 | 22828.17 | 0.6067 | 0.0607 | -20.2759
(Galimberti et al., 2017)$\mathsection$ |  |  |  | Regret | -11.6944 | 3 | -11.6963 | 513,627 | 888,353 | 14.78 | 37266.95 | 37085.30 | 0.1045 | 0.0114 | -13.4481

  * $*$

http://networkrepository.com

  * $\dagger$

http://snap.stanford.edu

  * $\ddagger$

https://comunelab.fbk.eu/data.php

  * $\mathsection$

https://goo.gl/8741Gs

Next we investigate the case where a clique is planted in only one layer. We
generate $k~{}(=1,2,3,4,5)$ power-law graphs with $\beta=3.0$ on $V$ with
$|V|=\text{1,000}$. Then we randomly select $V_{c}\subseteq V$ consisting of
20 vertices, but plant a clique on $V_{c}$ only in one randomly selected
layer. The performance of optimal solutions are described in Figure 2. As can
be seen, the F measure becomes smaller as the number of layers increases.
Among the three metrics, the regret performs particularly well because it
concentrates on the layer containing the clique from its definition;
therefore, the regret metric seems most suitable for robust analysis with
noisy layers. The robust ratio performs second best; it also cares about the
layer containing the clique.

Finally we conclude this subsection by evaluating vertex subsets that are
obtained from optimal solutions $p\in\Delta(2^{V})$. To verify the validity of
our model, we select a vertex subset attaining the highest probability in
$p\in\Delta(2^{V})$. As our algorithm does not know about $V_{c}$, it cannot
select a vertex subset using the F measure. For reference, we also run two
baseline algorithms, DCS-LP and DCS-Greedy, designed by Jethava and
Beerenwinkel (Jethava and Beerenwinkel, 2015). Note that DCS-LP can be seen as
the algorithm that selects a vertex subset with the largest minimum density
value over layers from the support of $p\in\Delta(2^{V})$, where
$p\in\Delta(2^{V})$ is an optimal solution to our algorithm with the density
metric. The results are depicted in Figure 3, where we employed the same
experimental settings as above and used the usual (deterministic) F measure
for evaluation. As can be seen, the trend of our algorithm is similar to that
observed in the above experiments; all metrics almost detect $V_{c}$ for the
all-layers setting, but only the regret metric is successful for the only-one-
layer setting. As for the baseline methods, DCS-LP is comparable to our
algorithm with the density metric, while DCS-Greedy performs quite poorly.

###  7.2. Performance of our algorithm

Here we examine the performance of our proposed algorithm using publicly-
available real-world multilayer networks. Table 1 shows the results together
with the summary of the characteristics of the datasets. Even for large
networks, our algorithm obtains an optimal solution in reasonable time. The
preprocessing algorithm often reduces the size of the networks significantly
using a reasonably large lower bound computed by LP (25) and results in a
substantial speed-up. As an extreme example, for Cannes2013 with the density
and robust ratio metrics, the lower bound attains the optimal value, and the
number of vertices is reduced by more than 99.9%, which makes the computation
more than $300$ times faster. Consistent with our theoretical analysis,
$|\text{supp}|$ is at most $k$. For AS-733, $k$ is quite large but
$|\text{supp}|$ is still small.

Finally we evaluate the single vertex subsets obtained from optimal solutions.
For $p\in\Delta(2^{V})$, we select a vertex subset from
$\operatorname{supp}(p)$ that optimizes the metric employed in the algorithm.
Note that in this setting, the output of DCS-LP coincides with that obtained
by our algorithm with the density metric. As DCS-Greedy performed quite
poorly, it is omitted. The results are shown in the last three columns of
Table 1. Although there are a few exceptions, the algorithm with a metric
performs best in terms of the metric employed. A critical fact is that
depending only on the density metric, we may fail to obtain meaningful
structure from networks. For example, the algorithm with the robust ratio
admits a particularly large robust ratio value of $0.6652$ for
NYClimateMarch2014, meaning that the vertex subset obtained achieves an
approximation ratio of $0.6652$ for all layers. Moreover, the algorithm with
the regret metric admits a particularly small regret value of 155.5000 for
Cannes2013. For those instances, the algorithm with the density metric (i.e.,
DCS-LP) performs poorly in terms of those metrics, respectively. From the
above, it seems quite important to select an appropriate metric depending on
the practical purpose at hand.

##  8\. Conclusion

In this paper, we have introduced a novel optimization model and algorithms
for dense subgraph discovery in multilayer networks. There are several
possible directions for future research. One direction is to improve the
scalability of our algorithm, particularly for the regret metric, for which
our preprocessing algorithm does not necessarily perform well. Another
direction is to apply our model and algorithms to some real-world applications
of multilayer-network analysis. Investigating multilayer-network counterparts
of some existing generalizations of the densest subgraph problem (see e.g.,
(Kawase and Miyauchi, 2018; Miyauchi and Kakimura, 2018; Tsourakakis, 2015;
Veldt et al., 2021)) is also interesting future work.

###### Acknowledgements.

This work was partially supported by JST PRESTO Grant Number JPMJPR2122 and
JSPS KAKENHI Grant Numbers JP17K12646, JP19K20218, JP20K19739, JP21K17708, and
JP21H03397.

## References

  * (1)
  * Angel et al. (2012) A. Angel, N. Sarkas, N. Koudas, and D. Srivastava. 2012\.  Dense subgraph maintenance under streaming edge weight updates for real-time story identification. In _Proceedings of VLDB_. 574–585.
  * Başar and Olsder (1999) T. Başar and G. J. Olsder. 1999.  _Dynamic Noncooperative Game Theory_. Classics in Applied Mathematics, Vol. 23.  SIAM.
  * Bader and Hogue (2003) G. D. Bader and C. W. V. Hogue. 2003.  An automated method for finding molecular complexes in large protein interaction networks.  _BMC Bioinformatics_ 4, 1 (2003), 1–27.
  * Balalau et al. (2015) O. D. Balalau, F. Bonchi, T-H. H. Chan, F. Gullo, and M. Sozio. 2015.  Finding subgraphs with maximum total density and limited overlap. In _Proceedings of WSDM_. 379–388.
  * Basaras et al. (2019) P. Basaras, G. Iosifidis, D. Katsaros, and L. Tassiulas. 2019.  Identifying influential spreaders in complex multilayer networks: A centrality perspective.  _IEEE Transactions on Network Science and Engineering_ 6, 1 (2019), 31–45.
  * Bazzi et al. (2016) M. Bazzi, M. A. Porter, S. Williams, M. McDonald, D. J. Fenn, and S. D. Howison. 2016\.  Community detection in temporal multilayer networks, with an application to correlation networks.  _Multiscale Modeling & Simulation_ 14, 1 (2016), 1–41.
  * Boccaletti et al. (2014) S. Boccaletti, G. Bianconi, R. Criado, C.I. del Genio, J. Gómez-Gardeñes, M. Romance, I. Sendiña-Nadal, Z. Wang, and M. Zanin. 2014.  The structure and dynamics of multilayer networks.  _Physics Reports_ 544, 1 (2014), 1–122.
  * Boob et al. (2020) D. Boob, Y. Gao, R. Peng, S. Sawlani, C. E. Tsourakakis, D. Wang, and J. Wang. 2020.  Flowless: Extracting densest subgraphs without flow computations. In _Proceedings of The Web Conference 2020_. 573–583.
  * Charikar (2000) M. Charikar. 2000\.  Greedy approximation algorithms for finding dense components in a graph. In _Proceedings of APPROX_. 84–95.
  * Charikar et al. (2018) M. Charikar, Y. Naamad, and J. Wu. 2018.  On finding dense common subgraphs.  _arXiv preprint arXiv:1802.06361_ (2018).
  * Chekuri et al. (2022) C. Chekuri, K. Quanrud, and M. R. Torres. 2022.  Densest subgraph: Supermodularity, iterative peeling, and flow. In _Proceedings of SODA_. 1531–1555.
  * Chung and Lu (2002) F. Chung and L. Lu. 2002\.  The average distances in random graphs with given expected degrees.  _Internet Mathematics_ 1 (2002), 15879–15882.
  * De Bacco et al. (2017) C. De Bacco, E. A. Power, D. B. Larremore, and C. Moore. 2017\.  Community detection, link prediction, and layer interdependence in multilayer networks.  _Physical Review E_ 95 (2017), 042317.
  * De Domenico et al. (2016) M. De Domenico, C. Granell, M. A. Porter, and A. Arenas. 2016\.  The physics of spreading processes in multilayer networks.  _Nature Physics_ 12, 10 (2016), 901–906.
  * De Domenico et al. (2013) M. De Domenico, A. Solé-Ribalta, E. Cozzo, M. Kivelä, Y. Moreno, M. A. Porter, S. Gómez, and A. Arenas. 2013\.  Mathematical formulation of multilayer networks.  _Physical Review X_ 3 (2013), 041022.
  * De Domenico et al. (2015) M. De Domenico, A. Solé-Ribalta, E. Omodei, S. Gómez, and A. Arenas. 2015\.  Ranking in interconnected multilayer networks reveals versatile nodes.  _Nature Communications_ 6 (2015), 6868.
  * Dourisboure et al. (2007) Y. Dourisboure, F. Geraci, and M. Pellegrini. 2007\.  Extraction and classification of dense communities in the web. In _Proceedings of WWW_. 461–470.
  * Firth and Sheldon (2015) J. A. Firth and B. C. Sheldon. 2015.  Experimental manipulation of avian social structure reveals segregation is carried over across contexts.  _Proceedings of the Royal Society B: Biological Sciences_ 282, 1802 (2015), 20142350.
  * Friston (2011) K. J. Friston. 2011\.  Functional and effective connectivity: A review.  _Brain Connectivity_ 1, 1 (2011), 13–36.
  * Galimberti et al. (2017) E. Galimberti, F. Bonchi, and F. Gullo. 2017.  Core decomposition and densest subgraph in multilayer networks. In _Proceedings of CIKM_. 1807–1816.
  * Gibson et al. (2005) D. Gibson, R. Kumar, and A. Tomkins. 2005.  Discovering large dense subgraphs in massive graphs. In _Proceedings of VLDB_. 721–732.
  * Gionis and Tsourakakis (2015) A. Gionis and C. E. Tsourakakis. 2015.  Dense subgraph discovery: KDD 2015 Tutorial. In _Proceedings of KDD_. 2313–2314.
  * Goldberg (1984) A. V. Goldberg. 1984\.  _Finding a maximum density subgraph_.  Technical Report. University of California Berkeley.
  * Hashemi et al. (2022) F. Hashemi, A. Behrouz, and L. V.S. Lakshmanan. 2022.  FirmCore decomposition of multilayer networks. In _Proceedings of The Web Conference_. 1589–1600.
  * Interdonato et al. (2017) R. Interdonato, A. Tagarelli, D. Ienco, A. Sallaberry, and P. Poncelet. 2017.  Local community detection in multilayer networks.  _Data Mining and Knowledge Discovery_ 31, 5 (2017), 1444–1479.
  * Jalili et al. (2017) M. Jalili, Y. Orouskhani, M. Asgari, N. Alipourfard, and M. Perc. 2017.  Link prediction in multiplex online social networks.  _Royal Society Open Science_ 4, 2 (2017), 160863\.
  * Jethava and Beerenwinkel (2015) V. Jethava and N. Beerenwinkel. 2015.  Finding dense subgraphs in relational graphs. In _Proceedings of ECML PKDD_. 641–654.
  * Kawase et al. (2019) Y. Kawase, Y. Kuroki, and A. Miyauchi. 2019.  Graph mining meets crowdsourcing: Extracting experts for answer aggregation. In _Proceedings of IJCAI_. 1272–1279.
  * Kawase and Miyauchi (2018) Y. Kawase and A. Miyauchi. 2018.  The densest subgraph problem with a convex/concave size function.  _Algorithmica_ 80, 12 (2018), 3461–3480.
  * Kivelä et al. (2014) M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A. Porter. 2014\.  Multilayer networks.  _Journal of Complex Networks_ 2, 3 (2014), 203–271.
  * Kortsarz and Peleg (1994) G. Kortsarz and D. Peleg. 1994.  Generating sparse 2-spanners.  _Journal of Algorithms_ 17, 2 (1994), 222–236.
  * Lee et al. (2010) V. E. Lee, N. Ruan, R. Jin, and C. Aggarwal. 2010\.  _A survey of algorithms for dense subgraph discovery_.  303–336.
  * Leskovec et al. (2005) J. Leskovec, J. Kleinberg, and C. Faloutsos. 2005\.  Graphs over time: Densification laws, shrinking diameters and possible explanations. In _Proceedings of KDD_. 177–187.
  * Miyauchi and Kakimura (2018) A. Miyauchi and N. Kakimura. 2018.  Finding a dense subgraph with sparse cut. In _Proceedings of CIKM_. 547–556.
  * Miyauchi and Takeda (2018) A. Miyauchi and A. Takeda. 2018.  Robust densest subgraph discovery. In _Proceedings of ICDM_. 1188–1193.
  * Omodei et al. (2015) E. Omodei, M. De Domenico, and A. Arenas. 2015\.  Characterizing interactions in online social networks during exceptional events.  _Frontiers in Physics_ 3 (2015), 59.
  * Salehi et al. (2015) M. Salehi, R. Sharma, M. Marzolla, M. Magnani, P. Siyari, and D. Montesi. 2015\.  Spreading processes in multilayer networks.  _IEEE Transactions on Network Science and Engineering_ 2, 2 (2015), 65–83.
  * Semertzidis et al. (2019) K. Semertzidis, E. Pitoura, E. Terzi, and P. Tsaparas. 2019\.  Finding lasting dense subgraphs.  _Data Mining and Knowledge Discovery_ 33, 5 (2019), 1417–1445.
  * Tagarelli et al. (2017) A. Tagarelli, A. Amelio, and F. Gullo. 2017.  Ensemble-based community detection in multilayer networks.  _Data Mining and Knowledge Discovery_ 31, 5 (2017), 1506–1543.
  * Tsourakakis (2015) C. E. Tsourakakis. 2015\.  The k-clique densest subgraph problem. In _Proceedings of WWW_. 1122–1132.
  * Tsourakakis et al. (2019) C. E. Tsourakakis, T. Chen, N. Kakimura, and J. Pachocki. 2019\.  Novel dense subgraph discovery primitives: Risk aversion and exclusion queries. In _Proceedings of ECML PKDD_. 378–394.
  * Vanderbei (2020) R. J. Vanderbei. 2020\.  _Linear Programming: Foundations and Extensions_. International Series in Operations Research & Management Science, Vol. 196.  Springer.
  * Veldt et al. (2021) N. Veldt, A. R. Benson, and J. Kleinberg. 2021.  The generalized mean densest subgraph problem. In _Proceedings of KDD_. 1604–1614.
  * Zou (2013) Z. Zou. 2013.  Polynomial-time algorithm for finding densest subgraphs in uncertain graphs. In _Proceedings of MLG_.  No page numbers.