# Non-uniqueness for the ab-family of equations

John Holmes and Rajan Puri

###### Abstract.

We study the cubic ab-family of equations, which includes both the Fokas-
Olver-Rosenau-Qiao (FORQ) and the Novikov (NE) equations. For $a\neq 0$, it is
proved that there exist initial data in the Sobolev space $H^{s}$, $s<3/2$,
with non-unique solutions. Multiple solutions are constructed by studying the
collision of 2-peakon solutions. Furthermore, we prove the novel phenomenon
that for some members of the family, collision between 2-peakons can occur
even if the “faster” peakon is in front of the “slower” peakon.

###### Key words and phrases:

Well-posedness, initial value problem, Cauchy problem, Sobolev spaces,
Camassa-Holm equation, solitons, peakons.

###### 2010 Mathematics Subject Classification:

Primary: 35Q35

Department of Mathematics and Statistics, Wake Forest University, Winston
Salem, NC 27109, USA<EMAIL_ADDRESS>
Department of Mathematics and Statistics, Wake Forest University, Winston
Salem, NC 27109, USA, (purir@wfu.edu).

##  1\. Introduction

We consider the Cauchy problem for the ab-equation

$\displaystyle
u_{t}+u^{2}u_{x}-au_{x}^{3}+D^{-2}\partial_{x}\big{[}\frac{b}{3}u^{3}+\frac{6-6a-b}{2}uu_{x}^{2}\big{]}+D^{-2}\big{[}\frac{2a+b-2}{2}u_{x}^{3}\big{]}=0$
(1) $\displaystyle u(x,0)=u_{0}(x),$ (2)

where $u_{0}(x)\in H^{s}(\mathbb{R})$, for $s<3/2$. The two parameters $a$,
$b\in\mathbb{R}$ and we will assume $a\neq 0$. $u_{t}$ and $u_{x}$ denote the
derivatives of $u$ with respect to $t$ and $x$, $\partial_{x}$ denotes
differentiation with respect to $x$, and the non-local operator
$D^{-2}=(1-\partial_{x}^{2})^{-1}$ is the inverse Fourier transform of
$(1+\xi^{2})^{-1}$.

Peakons, or peaked traveling wave solutions, were discovered in 1978 by
Fornberg and Whitham [8] and then by Camassa and Holm [CH] in their quest for
a water wave model that could capture wave breaking. These are special
traveling wave solutions which take the form

$u(x,t)=ce^{-|x-q(t)|}.$

The Camassa-Holm (CH) equation

$\displaystyle(1-\partial_{x})^{2}u_{t}=uu_{xxx}+2u_{x}u_{xx}-3uu_{x},$ (3)

and later the Degasperis-Procesi (DP) equation

$\displaystyle(1-\partial_{x})^{2}u_{t}=uu_{xxx}+3u_{x}u_{xx}-4uu_{x},$ (4)

introduced in [6] are two such integral equations which admit peakon
solutions. In fact, these equations admit special 2-peakon solutions, called
peakon-antipeakon solutions which are of the form

$u(x,t)=p(t)e^{-|x+q(t)|}-p(t)e^{-|x-q(t)|},$

where $p(t)$ and $q(t)$, magnitude and position, satisfy a system of ordinary
differential equations (ODEs). For a discussion of peakon and other solutions,
we refer the reader to Holm and Ivanov [15]. Himonas, Holliman and Grayshan
[12] show that the peak and antipeak move towards each other, collide in
finite time and ill-posedness of the corresponding Cauchy problems in $H^{s}$,
$s<3/2$ is established due to this event. This symmetry does not exist in the
ab-family of equations.

Two members of ab-family (1) have been studied extensively in different
contexts by several researchers. In particular, the choice of parameters
$a=1/3$, $b=2$ corresponds to the Fokas-Olver-Rosenau-Qiao (FORQ) equation

$u_{t}+u^{2}u_{x}-u_{x}^{3}+D^{-2}\partial_{x}\big{[}\frac{2}{3}u^{3}+uu_{x}^{2}\big{]}+D^{-2}\big{[}\frac{1}{3}u_{x}^{3}\big{]}=0,$
(5)

derived in Fokas [7], Fuchssteiner [9], Olver and Rosenau [17], and Qiao [18];
while the choice a = 0, b = 3 gives the Novikov equation (NE)

$u_{t}+u^{2}u_{x}+D^{-2}\partial_{x}\big{[}\frac{1}{3}u^{3}+\frac{3}{2}uu_{x}^{2}\big{]}+D^{-2}\big{[}\frac{1}{2}u_{x}^{3}\big{]}=0,$
(6)

derived by Novikov [16]. Both of these equations are integrable equations with
a bi-Hamiltonian structure, Lax pair, an infinite hierarchy of hereditary
symmetries. and infinitely many conserved quantities, for a more detailed
discussion of integrability see [19] and the references contained therein. In
particular, classical solutions to both equations conserve the $H^{1}$ norm.

The well-posedness theory for the ab-family is not completely understood. Some
partial results are as follows. Himonas and Mantzavinos [1] showed that FORQ
is well-posed in $H^{s}$, with $s>5/2$ and this was extended to a four-
parameter family in [2] which includes the $ab$-family. A nonuniqueness result
by Himonas and Holliman [11] showed that the FORQ equation is ill-posed in
$H^{s}$ for any $s<3/2$. There is no theory concerning well-posedness in the
gap $3/2\leq s\leq 5/2$. In contrast, the NE is well-posed in $H^{s}$ for all
$s>3/2$, see [10] for details, and Himonas, Kenig and Holliman [13] showed
ill-posedness in $H^{s}$ for $s<3/2$. Both the ill-posedness results for the
NE equation and the FORQ equation study the behavior of the solution near the
time of collision of a 2-peakon solution. The method used to prove ill-
posedness for both the NE and FORQ equations is similar to that used by many
authors for other nonlinear evolution equations. For instance, Bourgain and
Pavlovic [4] studied the Navier-Stokes equations in Besov spaces and Christ,
Colloander and Tao [5] used a similar technique to study defocusing dispersive
PDEs. We too use this idea to prove our main theorem, which is stated as
follows.

######  Theorem 1.1.

For all $b\in\mathbb{R}$ and $a\neq 0$, solutions to the Cauchy problem for
the ab-equation are not unique in $H^{s}$ when $s<3/2$.

Perhaps the most interesting phenomenon discovered in our proof, is that
multipeakon solutions to the ab-equation interact unlike classical solitons.
The Korteweg-de Vries (KdV) equation is the canonical example of a nonlinear
PDE supporting multi-soliton solutions. When two solitons collide, although
they interact in a nonlinear way, they decouple; both before and after the
collision, they act similarly to a linear superposition of two solitons.
Moreover, since their speed is proportional to their height, if a larger
soliton begins in front of a smaller soliton, they will never collide. For
more information describing their interaction, see Benes, Kasman and Young
[3].

A similar phenomenon has been observed of the 2-peakon solutions of Camassa-
Holm type equations. Shown below is a numerical solutions to the FORQ
equation, with the amplitudes of the two peaks shown in dotted and dashed,
while the distance between the peaks is shown in solid. Numerically, the two
peakons pass through each other relatively undisturbed, as shown in Figure 1.
(However, we note that for all of the 2-peakon solutions we consider, since
solutions are not unique, there are other solutions which are not described by
the numerical solution.) In contrast, for some values of $a$ and $b$, our
2-peakon solutions are entangled and do not seem to separate; we prove that
smaller solitons can begin behind larger solitons, and yet, they collide. In
Section 2, we contrast the interaction of the 2-peakon solutions of the ab-
equation with the solution of the FORQ equation.

Figure 1. FORQ wave heights (dotted and dashed) and relative positions
(solid).

In the next section, we will begin by studying a system of ODEs which
describes the time evolution of a 2-peakon solution to the ab-equation. We
will show that the solution to the ODE system does not blow up as the two
peakons collide. Then, we will show that at the time of collision, the
2-peakon solution reduces to a single peakon solution, which will complete the
proof of nonuniqueness.

##  2\. The ODE System

We begin by recalling the relationship between the 2-peakon solutions of the
ab-equation and a system of ODEs. We will study the system of ODEs, prove some
preliminary results necessary in the proof of ill-posedness, then choose the
initial data which we can show converges to a single peakon solution.

######  Lemma 2.1.

[14] The 2-peakon function

$u(x,t)=p_{1}(t)e^{-|x-q_{1}(t)|}+p_{2}(t)e^{-|x-q_{2}(t)|},$ (7)

solves the ab-equation if and only if the positions $q_{1},q_{2}$ and the
momenta $p_{1},p_{2}$ satisfy the $4\times 4$ system of ordinary differential
equations:

$\begin{split}q_{1}^{{}^{\prime}}&=(1-a)p_{1}^{2}+2p_{1}p_{2}e^{-|q_{1}-q_{2}|}+(1-3a)p_{2}^{2}e^{-2|q_{1}-q_{2}|}\\\
q_{2}^{{}^{\prime}}&=(1-a)p_{2}^{2}+2p_{1}p_{2}e^{-|q_{1}-q_{2}|}+(1-3a)p_{1}^{2}e^{-2|q_{1}-q_{2}|}\\\
p_{1}^{{}^{\prime}}&=(2-b)sgn(q_{2}-q_{1})p_{1}p_{2}e^{-|q_{1}-q_{2}|}(p_{1}+p_{2}e^{-|q_{1}-q_{2}|})\\\
p_{2}^{{}^{\prime}}&=(2-b)sgn(q_{1}-q_{2})p_{1}p_{2}e^{-|q_{1}-q_{2}|}(p_{1}e^{-|q_{1}-q_{2}|}+p_{2}).\end{split}$
(8)

We first observe that if at time zero, $p_{1}(0),p_{2}(0)\neq 0$, and
$q_{1}(0)\neq q_{2}(0)$ (ie we have a multi-peakon initial profile), then, due
to the time reversibility of solutions, nonuniqueness will be shown if either
$p_{1}(t)$ or $p_{2}(t)$ equal zero at some time, or if $q_{1}(t)=q_{2}(t)$
while all the functions remain bounded.

We define $q(t)\ \dot{=}\ q_{2}-q_{1},$ $h\ \dot{=}\ p_{2}-p_{1}$,
$w\dot{=}p_{1}+p_{2}$ and $z\ \dot{=}\ p_{1}p_{2}$. We will assume that at
initially, $q_{1}(0)=0$ and $q_{2}(0)=\mu$, and therefore, by continuity,
$q(t)\geq 0$ for some $0\leq t\leq T^{c}$, where $T^{c}$ is the collision
time, i.e when $q(t)=0$ for the first time. From equation (8), we have

$\begin{split}q^{{}^{\prime}}=q_{2}^{{}^{\prime}}-q_{1}^{{}^{\prime}}&=\bigg{[}(1-a)p_{2}^{2}+2p_{1}p_{2}e^{-|q_{1}-q_{2}|}+(1-3a)p_{1}^{2}e^{-2|q_{1}-q_{2}|}\bigg{]}\\\
&-\bigg{[}(1-a)p_{1}^{2}+2p_{1}p_{2}e^{-|q_{1}-q_{2}|}+(1-3a)p_{2}^{2}e^{-2|q_{1}-q_{2}|}\bigg{]}\\\
&=(1-a)(p_{2}^{2}-p_{1}^{2})-(1-3a)(p_{2}^{2}-p_{1}^{2})e^{-2|q_{1}-q_{2}|}\\\
&=(p_{2}-p_{1})(p_{2}+p_{1})\bigg{(}(1-a)-(1-3a)e^{-2|q_{1}-q_{2}|}\bigg{)}.\end{split}$
(9)

Additionally, the difference $p_{2}^{{}^{\prime}}-p_{1}^{{}^{\prime}}$ and
$p_{2}^{{}^{\prime}}+p_{1}^{{}^{\prime}}$ satisfy

$\begin{split}h^{{}^{\prime}}=p_{2}^{{}^{\prime}}-p_{1}^{{}^{\prime}}&=\bigg{[}(2-b)sgn(q_{1}-q_{2})p_{1}p_{2}e^{-|q_{1}-q_{2}|}(p_{1}e^{-|q_{1}-q_{2}|}+p_{2})\bigg{]}\\\
&-\bigg{[}(2-b)sgn(q_{2}-q_{1})p_{1}p_{2}e^{-|q_{1}-q_{2}|}(p_{1}+p_{2}e^{-|q_{1}-q_{2}|})\bigg{]}\\\
&=-(2-b)p_{1}p_{2}e^{-|q_{1}-q_{2}|}(p_{1}+p_{2}+(p_{1}+p_{2})e^{-|q_{1}-q_{2}|})\\\
&=-(2-b)p_{1}p_{2}(p_{1}+p_{2})(1+e^{-|q_{1}-q_{2}|})e^{-|q_{1}-q_{2}|},\end{split}$
(10)

and,

$w^{{}^{\prime}}=(2-b)p_{1}p_{2}(p_{1}-p_{2})(1-e^{-|q_{1}-q_{2}|})e^{-|q_{1}-q_{2}|}.$
(11)

Therefore, from the above calculation, we deduce

$\displaystyle\begin{aligned} q^{\prime}&=hw(1-a-(1-3a)e^{-2q}),\\\
h^{\prime}&=-(2-b)wz(1+e^{-q})e^{-q},\\\
w^{\prime}&=-(2-b)hz(1-e^{-q})e^{-q},\\\
z^{\prime}&=(2-b)hwze^{-2q},\end{aligned}\qquad\begin{aligned}
q(0)&=q_{2}(0)-q_{1}(0)=\mu,\\\ h(0)&=p_{2}(0)-p_{1}(0)=h_{0},\\\
w(0)&=p_{1}(0)+p_{2}(0)=w_{0},\\\ z(0)&=p_{1}(0)p_{2}(0)=z_{0}.\end{aligned}$
(12)

If $a=1/3$, then $(1-a-(1-3a)e^{-2q})=2/3$, and we set $\mu>0$ to be a small
positive number. Else, we choose a number $1<c<2$ (there is always a choice of
$c$) such that

$0<\mu=\mu(a)=-\frac{1}{2}\ln\left(\frac{(c+1)a-1}{3a-1}\right)<1,$

which implies

$L_{a}(\mu)\ \dot{=}\ 1-e^{-2\mu}+3ae^{-2\mu}-a=ca.$

We note that as $q$ tends to $0$, $(1-a-(1-3a)e^{-2q})$ tends to $2a$, and
therefore, the sign of this term remains constant. Also, note that
$L_{a}(\mu)$ can never take the value $0$ when $a\neq 0$ on the domain
$\mu\geq 0$. The main results will require different initial data for $h$, $w$
and $z$ depending upon both $a$ and $b$. Therefore, we hold off on the choices
of $p_{1}(0)$ and $p_{2}(0)$ until necessary. Our first proposition describes
the solution to the above system. In particular, we will see that so long as
$q\geq 0$, the functions $h,w,$ and $z$ are bounded and therefore
$p_{1}=\frac{1}{2}(w-h)$ and $p_{2}=\frac{1}{2}(h+w)$ remain bounded as well.

######  Proposition 2.2.

If $a\neq 0$, for any real numbers $h_{0},w_{0},z_{0}$, the initial value
problem (12) has a unique smooth solution on some positive time interval.
Furthermore, in terms of $q$, the functions $h(t)$ and $w(t)$ and $z(t)$
remain bounded for all $0\leq q$.

###### Proof.

The right hand side of the system (12) is smooth in the arguments $q,h,w,z$
and therefore, by the ODE theorem, has a solution on some time interval
$[0,T)$, $T>0$. We now derive the relations between $h,w,z$ and $q$. Beginning
with $z$, from the equations with $z^{\prime}$ and $q^{\prime}$ we find

$\displaystyle\frac{z^{\prime}}{z}=\frac{(2-b)e^{-2q}q^{\prime}}{1-a-(1-3a)e^{-2q}}.$

Therefore, we have

$\displaystyle\frac{d}{dt}\ln|z|=\frac{(2-b)}{2(1-3a)}\frac{d}{dt}\ln|1-a-(1-3a)e^{-2q}|.$

If $z_{0}>0$, then by continuity, we can assume $z>0$ for some time and
$|z|=z$. Likewise, if $z_{0}<0$, we will assume $|z|=-z$. In either case,
integrating from $0$ to $t$ yields

$\displaystyle\ln\left(\frac{z(t)}{z_{0}}\right)=\frac{(2-b)}{2(1-3a)}\ln\left(\frac{1-a-(1-3a)e^{-2q}}{1-a-(1-3a)e^{-2\mu}}\right).$
(13)

We exponentiate and rearrange terms to find

$\displaystyle
z(t)=z_{0}\left(\frac{1-a-(1-3a)e^{-2q}}{1-a-(1-3a)e^{-2\mu}}\right)^{\frac{(2-b)}{2(1-3a)}}.$
(14)

From here we can conclude that $z(t)$ remains bounded for all $t\in[0,T^{c}]$,
since neither the numerator nor the denominator take the value $0$.

Now, we will use the above formula for $z(t)$ to find $h(t)$. We have

$\displaystyle\frac{h^{\prime}}{q^{\prime}}=\frac{-(2-b)wz(1+e^{-q})e^{-q}}{hw(1-a-(1-3a)e^{-2q})},$

or rearranging we find

$\displaystyle
hh^{\prime}=z\cdot\frac{-(2-b)(1+e^{-q})e^{-q}q^{\prime}}{1-a-(1-3a)e^{-2q}}.$

We substitute the formula found for $z$ in equation (14) to get

$\displaystyle
hh^{\prime}=z_{0}\left|1-a-(1-3a)e^{-2\mu}\right|^{\frac{-(2-b)}{2(1-3a)}}\frac{-(2-b)(1+e^{-q})e^{-q}q^{\prime}}{|1-a-(1-3a)e^{-2q}|^{1-\frac{(2-b)}{2(1-3a)}}}.$

Set

$\displaystyle f(q)\ \dot{=}\
z_{0}\left|1-a-(1-3a)e^{-2\mu}\right|^{\frac{-(2-b)}{2(1-3a)}}\frac{-(2-b)e^{-q}}{|1-a-(1-3a)e^{-2q}|^{1-\frac{(2-b)}{2(1-3a)}}},$
(15)

and then define

$F_{1}(q)\ \dot{=}\ \int_{0}^{t}(1+e^{-q})f(q)dq.$

Since $f(q)$ is smooth and bounded for all $q\geq 0$ (since the denominator is
singular only when $a=0$), $F(q)$ remains smooth, bounded and differentiable
for all $0\leq t\leq T^{c}$. Therefore,

$h^{2}=h_{0}^{2}+2F_{1}(q),$

remains bounded for all $0\leq t\leq T^{c}$.

Next we solve for $w(t)$. We rearrange the equations for $w^{\prime}$ and
$q^{\prime}$ to find

$\displaystyle
ww^{\prime}=z\cdot\frac{-(2-b)(1-e^{-q})e^{-q}q^{\prime}}{(1-a-(1-3a)e^{-2q})}=(1-e^{-q})f(q)q^{\prime},$

where we used the definition of $f(q)$ found in equation (15) and the formula
for $z$ found in equation (14). Defining

$F_{2}(q)\ \dot{=}\ \int_{0}^{t}(1-e^{-q})f(q)dq,$

we find

$w^{2}=w_{0}^{2}+2F_{2}(q),$

and similarly to $h$, $w$ remains smooth and bounded for all $0\leq t\leq
T^{c}$. ∎

Since $p_{1}=\frac{1}{2}(w-h)$ and $p_{2}=\frac{1}{2}(h+w)$, the above
proposition shows that as long as $q\geq 0$, $p_{1},p_{2}<\infty$. We will
next choose the initial values for $p_{1}(0)$ and $p_{2}(0)$ to be consistent
with the above proposition and which will necessarily lead to a collision time
$T^{c}<\infty$. We break our analysis into cases described by the parameters
$a$ and $b$. For each case, we will consider the function
$p(t)=p_{2}^{2}(t)-p_{1}^{2}(t)$. From (9),

$\displaystyle q^{\prime}(t)$
$\displaystyle=p(t)\big{(}(1-e^{-2|q(t)})+a(3e^{-2|q(t)|}-1)\big{)}$ (16)
$\displaystyle p^{\prime}(t)$
$\displaystyle=2(b-2)p_{1}p_{2}(p_{1}^{2}+p_{2}^{2})e^{-|q(t)|}+4(b-2)p_{1}^{2}p_{2}^{2}e^{-2|q(t)|}.$
(17)

We will show that $T^{c}$ exists and find an upper bound by showing
$q^{\prime}(t)$ is bounded by a negative number so long as $q(t)$ remains non-
negative. More precisely, we will prove the following lemma.

######  Lemma 2.3.

For all $b\in\mathbb{R}$ and $a\neq 0$, there exists an initial multipeakon
profile such that for some $\epsilon>0$, $\frac{dq}{dt}<-\epsilon<0$ for
$t\in[0,T^{c})$ (and hence $T^{c}<\infty$) or there exists a time $T^{p}$ such
that at least one of $p_{1}(T^{p})=0$ or $p_{2}(T^{p})=0$.

###### Proof.

If $\min\\{T^{p},T^{c}\\}=T^{p}$ and $T^{p}<\infty$, we are done. Therefore,
we will assume $p_{1}(t)$ and $p_{2}(t)$ do not equal zero. By continuity,
whatever the sign of their initial data is, we may assume the solutions take
as well. We prove the theorem in four cases, based upon the values of $a$ and
$b$, omitting the trivial case when $b=2$.

Case 1: $a>0,b>2$. We take the initial data

$p_{1}(0)=\alpha+\delta,\ q_{1}(0)=0,\ p_{2}(0)=-\alpha,\ q_{2}(0)=\mu,$

shown in Figure 3 and we shall prove that either a collision time exists, or
one of the peakons ceases to exist. In Figure 3 we show the time evolution of
$p_{1}(t)$ and $p_{2}(t)$ in dotted and dashed respectively, and in solid we
show the time evolution of $q(t)$.

By the choice of our initial data:

$p(0)=p_{2}^{2}(0)-p_{1}^{2}(0)=-(2\alpha\delta+\delta^{2})<0,\text{ and
}q(0)=q_{2}(0)-q_{1}(0)=\mu>0.$

Figure 2. Case 1: Peakon-antipeakon initial data.

Figure 3. The heights of the two peaks are dotted and dashed; $q$ is the solid
line.

Since $p_{2}(0)=-\alpha<0$ and $p_{1}(0)=\alpha+\delta>0,\ p_{1}(0)p_{2}(0)<0$
and by continuity $p_{1}(t)p_{2}(t)<0$.

We will now show that for all $0\leq t\leq T^{c}$,

$\frac{dp}{dt}=2(b-2)p_{1}p_{2}\bigg{(}(p_{1}^{2}+p_{2}^{2})e^{-|q|}+2p_{1}p_{2}e^{-2|q|}\bigg{)}<0.$

Indeed, the following calculation shows that
$\bigg{(}(p_{1}^{2}+p_{2}^{2})e^{-|q|}+2p_{1}p_{2}e^{-2|q|}\bigg{)}>0.$ Use
the fact that $e^{-2|x|}<e^{-|x|}$ to compute

$\bigg{(}(p_{1}^{2}+p_{2}^{2})e^{-|q|}+2p_{1}p_{2}e^{-2|q|}\bigg{)}\geq\bigg{(}(p_{1}^{2}+p_{2}^{2})e^{-2|q|}+2p_{1}p_{2}e^{-2|q|}\bigg{)}=(p_{1}+p_{2})^{2}e^{-2|q|}\geq
0.$

Therefore, we may now use $p(t)<p(0)<0$. Substituting this into equation (16),
we have

$\frac{dq}{dt}\leq p(0)\big{(}(1-e^{-2|q(t)|})+a(3e^{-2|q(t)|}-1)\big{)}.$

The right hand side is negative, since initially $0<a\leq L_{a}(\mu)=ca\leq
2a$, and as $q$ decreases,

$ca\leq(1-e^{-2|q(t)|})+a(3e^{-2|q(t)|}-1)\leq 2a.$

Therefore we compute $\epsilon$:

$\frac{dq}{dt}\leq ap(0)=-a(2\alpha\delta+\delta^{2})=-\epsilon<0$

Hence, either $p_{1}(t)=0$, $p_{2}(t)=0$, or $q(t)=0$ in finite time.

Case 2: $a>0,b<2$. In contrast to the peakon-antipeakon initial profile of
Case 1, we take the two peakon initial profile:

$p_{1}(0)=\alpha+\delta,\ q_{1}(0)=0,\ p_{2}(0)=\alpha,\ q_{2}(0)=\mu,$

shown in Figure 5 and we shall prove that either a collision time exists, or
one of the peakons ceases to exist. In contrast to the previous case,
numerical experiments show that there are many times in which $q(t)=0$;
$T^{c}$ is the first time this occurs. Moreover, the numerical solution shows
that the two peakons do not decouple, rather, they repeatedly leapfrog each
other.

Figure 4. Case 2: 2-Peakon initial data.

Figure 5. The heights of the two peaks are dotted and dashed; $q$ is the solid
line.

By the choice of our initial data

$p(0)=p_{2}^{2}(0)-p_{1}^{2}(0)=-(2\alpha\delta+\delta^{2})<0,\
q(0)=q_{2}(0)-q_{1}(0)=\mu>0.$

Since $p_{1}(0)p_{2}(0)>0$, by continuity $p_{1}(t)p_{2}(t)>0$ for
$t\in[0,T^{c})$. We have

$\frac{dp}{dt}=2(b-2)p_{1}p_{2}\bigg{(}(p_{1}^{2}+p_{2}^{2})e^{-|q|}+2p_{1}p_{2}e^{-2|q|}\bigg{)}<0.$

Thus, $p(t)<p(0)<0\ \forall t\in[0,T^{c}).$ Now from the equation (16),

$\frac{dq(t)}{dt}\leq p(0)\big{(}(1-e^{-2|q(t)})+a(3e^{-2|q(t)|}-1)\big{)}\
\text{for}\ t\in[0,T^{c}],$

and, similarly to the first case

$\frac{dq(t)}{dt}\leq ap(0)=-a(2\alpha\delta+\delta^{2})=-\epsilon<0.$

Thus, we have shown that there exists an initial profile such that
$p_{1}(t)=0$, $p_{2}(t)=0$, or $q(t)=0$ in finite time.

Case 3: $a<0,b>2$. Similar to Case 2, we take the two peakon initial profile:

$p_{1}(0)=\alpha,\ q_{1}(0)=0,\ p_{2}(0)=\alpha+\delta,\ q_{2}(0)=\mu,$

shown in Figure 7. As shown in Figure 7 the solution of this 2-peakon initial
data behaves similarly to Case 2.

Figure 6. Case 3: 2-Peakon initial data.

Figure 7. The heights of the two peaks are dotted and dashed; $q$ is the solid
line.

By the choice of our initial data

$p(0)=p_{2}^{2}(0)-p_{1}^{2}(0)=(2\alpha\delta+\delta^{2})>0,\
q(0)=q_{2}(0)-q_{1}(0)=\mu>0.$

By continuity $p_{1}(t)p_{2}(t)>0$ and we will assume this hold for
$t\in[0,T^{c})$, thus

$\frac{dp}{dt}=2(b-2)p_{1}p_{2}\bigg{(}(p_{1}^{2}+p_{2}^{2})e^{-|q|}+2p_{1}p_{2}e^{-2|q|}\bigg{)}>0.$

and therefore, $p(t)>p(0)>0$.

Recalling the choice of $\mu$, we have $L_{a}(\mu)=ca<0$ and
$\raisebox{2.15277pt}{\scalebox{0.8}{$\displaystyle\lim_{q\to\
0}\;$}}L_{a}(q)=2a<0.$ Hence,

$2a<L_{a}(q(t))<ca<0,t\in[0,T^{c}).$ (18)

Therefore, using equation (16) we have

$\frac{dq(t)}{dt}=p(t)L_{a}(q(t))\leq p(0)\cdot ca=-\epsilon<0.$

Again, this shows that there exists an initial profile such that $p_{1}(t)=0$,
$p_{2}(t)=0$, or $q(t)=0$ in finite time.

Case 4: $a<0,b<2$. Similar to Case 1, we take a peakon-antipeakon initial
profile:

$p_{1}(0)=-\alpha,\ q_{1}(0)=0,\ p_{2}(0)=\alpha+\delta,\ q_{2}(0)=\mu,$

and as before, we will assume $T^{c}\leq T^{p}$. As shown in Figure 9, this
profile is not the same as in Case 1; though the behavior of the solution
resembles closely Case 1 as shown in Figure 9.

Figure 8. Case 4: Antipeakon-peakon initial data.

Figure 9. The heights of the two peaks are dotted and dashed; $q$ is the solid
line.

By the choice of our initial data

$p(0)=p_{2}^{2}(0)-p_{1}^{2}(0)=(2\alpha\delta+\delta^{2})>0,\text{ and
}q(0)=q_{1}(0)-q_{2}(0)=\mu>0.$

Since $p_{1}(0)p_{2}(0)<0$, by continuity $p_{1}(t)p_{2}(t)<0$, and therefore
an argument similar to the argument presented in Case 1 shows

$\frac{dp}{dt}=2(b-2)p_{1}p_{2}\bigg{(}(p_{1}^{2}+p_{2}^{2})e^{-|q|}+2p_{1}p_{2}e^{-2|q|}\bigg{)}>0,$

thus $p(t)>p(0)>0\ \forall\ t\in[0,T^{c}).$ Therefore, using equation (16) and
the estimate in inequality (18) we have again

$\frac{dq(t)}{dt}=p(t)L_{a}(q(t))\leq acp(0)=-\epsilon<0.$

This completes the fourth case, and we have shown that for every choice of $a$
and $b$, there is an initial profile which leads to
$\min\\{T^{p},T^{c}\\}<\infty$. ∎

###  2.1. Proof of Theorem 1.1

We are now prepared to show nonuniqueness in $H^{s}$ for the ab-family of
equations stated in Theorem 1.1. To complete the proof, we shall show that at
time, $T=\min\\{T^{c},T^{p}\\}$, the solution to the Cauchy problem (1)-(2)
with the initial profiles (depending upon $a$ and $b$) given in the previous
subsection, is either a single peakon or the zero solution.

We define the collision function:

$C(x)=p^{*}e^{-|x-q^{*}|}.$ (19)

where if $T=T^{c}$,
$q^{*}=\raisebox{2.15277pt}{\scalebox{0.8}{$\displaystyle\lim_{t\to\
T^{-}}\;$}}q_{1}(t)$ is the location of the collision and
$p^{*}=\raisebox{2.15277pt}{\scalebox{0.8}{$\displaystyle\lim_{t\to\
T^{-}}\;$}}(p_{1}(t)+p_{2}(t))$ is the magnitude of the collision. If
$T=T^{p}$, then label $i$ and $j$ ($1$ and $2$) such that
$\raisebox{2.15277pt}{\scalebox{0.8}{$\displaystyle\lim_{t\to\
T^{-}}\;$}}p_{i}(t)=0$, and
$\raisebox{2.15277pt}{\scalebox{0.8}{$\displaystyle\lim_{t\to\
T^{-}}\;$}}p_{j}(t)\neq 0$. We define
$q^{*}=\raisebox{2.15277pt}{\scalebox{0.8}{$\displaystyle\lim_{t\to\
T^{-}}\;$}}q_{j}(t)$. If both $p_{1}$ and $p_{2}$ converge to zero at time
$T$, since $C(x)=0$, the choice in $q^{*}$ is irrelavent. We will now show
that the solution converges to the collision function.

######  Lemma 2.4.

The $H^{s}$ limit of $u$, as $t$ approaches $T$ from below is C:

$\lim_{t\to\ T^{-}}||u(t)-C||_{H^{s}}=0.$

###### Proof.

We take the Fourier transform of $u$, and we have

$\widehat{u}(\xi,t)=\frac{2p_{1}e^{-i\xi
q_{1}}}{1+\xi^{2}}+\frac{2p_{2}e^{-i\xi q_{2}}}{1+\xi^{2}}.$

Similarly, we can find the Fourier transform of $C$ as

$\widehat{C}(\xi)=\frac{2p^{*}e^{-i\xi q^{*}}}{1+\xi^{2}}.$

Calculating the $H^{s}$ norm of $u(t)-C$ gives us

$\lim_{t\to\ T^{-}}||u(t)-C||_{H^{s}}^{2}=4\lim_{t\to\ T^{-}}\int_{\mbox{\bbc
R}}({1+\xi^{2}})^{s-2}|p_{1}e^{-i\xi q_{1}}+p_{2}e^{-i\xi q_{2}}-p^{*}e^{-i\xi
q^{*}}|^{2}d\xi.$

We can bound the quantity inside the absolute value by
$(|p_{1}|+|p_{2}|+|p^{*}|)\leq M<\infty.$ Let $v(\xi)=(1+\xi^{2})^{s-2}\cdot
M^{2}$ then $v$ dominates our original integral and $v$ is itself integrable
when $s<3/2.$ Therefore, we may apply the Dominated Convergence Theorem and
bring the limit inside the integral.

$\lim_{t\to\ T^{-}}||u(t)-C||_{H^{s}}^{2}=4\int_{\mbox{\bbc
R}}({1+\xi^{2}})^{s-2}|p_{1}(T)e^{-i\xi q_{1}(T)}+p_{2}(T)e^{-i\xi
q_{2}(T)}-p^{*}e^{-i\xi q^{*}}|^{2}d\xi.$ (20)

By definition of $p^{*}$ and $q^{*}$, the term inside the integral is zero. ∎

## References

  * [1] A. Alexandrou Himonas and Dionyssios Mantzavinos.  The Cauchy problem for the Fokas-Olver-Rosenau-Qiao equation.  Nonlinear Anal., 95:499–529, 2014.
  * [2] A. Alexandrou Himonas and Dionyssios Mantzavinos.  The Cauchy problem for a 4-parameter family of equations with peakon traveling waves.  Nonlinear Anal., 133:161–199, 2016.
  * [3] N. Benes, A. Kasman, and K. Young.  On decompositions of the KdV 2-soliton.  J. Nonlinear Sci., 16(2):179–200, 2006.
  * [4] Jean Bourgain and Nataša Pavlović.  Ill-posedness of the Navier-Stokes equations in a critical space in 3D.  J. Funct. Anal., 255(9):2233–2247, 2008.
  * [5] Michael Christ, James Colliander, and Terrence Tao.  Asymptotics, frequency modulation, and low regularity ill-posedness for canonical defocusing equations.  Amer. J. Math., 125(6):1235–1293, 2003.
  * [6] A. Degasperis, S. De Lillo, and F. Tinebra.  Solutions of $1+1$ relativistic field models.  J. Math. Phys., 39(1):497–531, 1998.
  * [7] A. S. Fokas.  On a class of physically important integrable equations.  volume 87, pages 145–150. 1995.  The nonlinear Schrödinger equation (Chernogolovka, 1994).
  * [8] B. Fornberg and G. B. Whitham.  A numerical and theoretical study of certain nonlinear wave phenomena.  Philos. Trans. Roy. Soc. London Ser. A, 289(1361):373–404, 1978\.
  * [9] Benno Fuchssteiner.  Some tricks from the symmetry-toolbox for nonlinear equations: generalizations of the Camassa-Holm equation.  Phys. D, 95(3-4):229–243, 1996.
  * [10] A. Alexandrou Himonas and Curtis Holliman.  The Cauchy problem for the Novikov equation.  Nonlinearity, 25(2):449–479, 2012.
  * [11] A. Alexandrou Himonas and Curtis Holliman.  Non-uniqueness for the Fokas-Olver-Rosenau-Qiao equation.  J. Math. Anal. Appl., 470(1):647–658, 2019.
  * [12] A. Alexandrou Himonas, Curtis Holliman, and Katelyn Grayshan.  Norm inflation and ill-posedness for the Degasperis-Procesi equation.  Comm. Partial Differential Equations, 39(12):2198–2215, 2014.
  * [13] A. Alexandrou Himonas, Curtis Holliman, and Carlos Kenig.  Construction of 2-peakon solutions and ill-posedness for the Novikov equation.  SIAM J. Math. Anal., 50(3):2968–3006, 2018.
  * [14] A. Alexandrou Himonas and Dionyssios Mantzavinos.  An $ab$-family of equations with peakon traveling waves.  Proc. Amer. Math. Soc., 144(9):3797–3811, 2016.
  * [15] Darryl D. Holm and Rossen I. Ivanov.  Smooth and peaked solitons of the CH equation.  J. Phys. A, 43(43):434003, 18, 2010.
  * [16] Vladimir Novikov.  Generalizations of the Camassa-Holm equation.  J. Phys. A, 42(34):342002, 14, 2009.
  * [17] Peter J. Olver and Philip Rosenau.  Tri-Hamiltonian duality between solitons and solitary-wave solutions having compact support.  Phys. Rev. E (3), 53(2):1900–1906, 1996.
  * [18] Zhijun Qiao.  A new integrable equation with cuspons and W/M-shape-peaks solitons.  J. Math. Phys., 47(11):112701, 9, 2006.
  * [19] Zhijun Qiao.  New integrable hierarchy, its parametric solutions, cuspons, one-peak solitons, and M/W-shape peak solitons.  J. Math. Phys., 48(8):082701, 20, 2007.


# Multi-task Learning for Gaussian Graphical Regressions with High Dimensional
Covariates

Jingfei Zhang$~{}^{1}$ and Yi Li$~{}^{2}$
1 Department of Management Science,
University of Miami, Coral Gables, FL 33146.
2 Department of Biostatistics,
University of Michigan, Ann Arbor, MI 48109.

###### Abstract

Gaussian graphical regression is a powerful means that regresses the precision
matrix of a Gaussian graphical model on covariates, permitting the numbers of
the response variables and covariates to far exceed the sample size. Model
fitting is typically carried out via separate node-wise lasso regressions,
ignoring the network-induced structure among these regressions. Consequently,
the error rate is high, especially when the number of nodes is large. We
propose a multi-task learning estimator for fitting Gaussian graphical
regression models; we design a cross-task group sparsity penalty and a within
task element-wise sparsity penalty, which govern the sparsity of active
covariates and their effects on the graph, respectively. For computation, we
consider an efficient augmented Lagrangian algorithm, which solves subproblems
with a semi-smooth Newton method. For theory, we show that the error rate of
the multi-task learning based estimates has much improvement over that of the
separate node-wise lasso estimates, because the cross-task penalty borrows
information across tasks. To address the main challenge that the tasks are
entangled in a complicated correlation structure, we establish a new tail
probability bound for correlated heavy-tailed (sub-exponential) variables with
an arbitrary correlation structure, a useful theoretical result in its own
right. Finally, the utility of our method is demonstrated through simulations
as well as an application to a gene co-expression network study with brain
cancer patients.

##  1 Introduction

Gaussian graphical models are an effective tool for inferring the dependence
among variables of interest, such as the co-expression patterns among genes
(e.g., Cai et al., 2012; Chen et al., 2016) and functional connectivity
between brain regions (e.g., Zhang et al., 2019), because precision matrices
for multivariate Gaussian variables can be interpreted as partial correlations
(Peng et al., 2009). What has been overlooked is that graphical structures may
depend on external covariates, which can be high dimensional. For example,
genetic variants, clinical and environmental factors, may affect both the
expression levels of individual genes and the co-expression relationships
among genes (Wang et al., 2012, 2013; Brynedal et al., 2017). Genetic variants
that alter co-expression relationships are referred to as co-expression
quantitative trait loci (QTLs), and identifying them is of keen scientific
interest (Wang et al., 2012, 2013; van der Wijst et al., 2018). Identifying
co-expression QTLs among many SNPs can be formulated as a problem that relates
graphical models to high dimensional external covariates, and a fundamental
interest is to ascertain how the covariates modulate the individual-level
graphical structures, and recover both the population- and subject-level
graphs.

Though much progress has been made for developing graphical models, less has
been done for covariate-dependent graphical models. Several works (e.g., Yin
and Li, 2011; Li et al., 2012; Cai et al., 2012; Chen et al., 2016) considered
covariate-dependent Gaussian graphical models, wherein the mean of the nodes
depends on covariates, while the network structure is the same across all of
the subjects. Guo et al. (2011) and Danaher et al. (2014) estimated stratified
graphical models by preserving the common structure among them. Liu et al.
(2010) partitioned the covariate space into several subspaces and fitted
separate Gaussian graphical models for each subspace. Kolar et al. (2010)
nonparametrically estimated the dependence of a covariance matrix on one
continuous covariate. Cheng et al. (2014) fitted a conditional Ising model for
binary data. Ni et al. (2019) proposed a DAG model that allows the graph
structure to vary with a small number of covariates, and assumed a
hierarchical ordering of the nodes. None of the cited works, however, can
handle high dimensional covariates. Notably, Zhang and Li (2022) proposed a
Gaussian graphical regression framework that regresses the precision matrix of
a graphical model on high dimensional covariates and estimates the parameters
using penalized node-wise regressions. However, Zhang and Li (2022) utilized
separate node-wise regressions, ignoring the network-induced common structure
among these regressions. Consequently, the obtained error rate is high,
especially with a large number of nodes.

This work addresses these issues by making several methodological and
theoretical contributions.

  * •

We propose to adapt multi-task learning to fit Gaussian graphical regression
models. Specifically, we design a cross-task group sparsity penalty which
enables us to borrow information across different tasks. As the number of
covariates, $q$, can be much larger than $n$, it is reasonable to assume the
active covariates (i.e., those with nonzero effects on the graph) are sparse,
which can also be accommodated by the specified cross-task group sparsity
penalty; see Figure 1. We combine it with an element-wise sparsity to address
the possible sparse effects of active covariates and further reduce the model
complexity.

  * •

To optimize the cross-task objective function with the combined sparsity
penalty, we adapt an efficient augmented Lagrangian algorithm, which solves
subproblems with a semi-smooth Newton method.

  * •

To address the main challenge that the regression tasks are entangled in a
general correlation structure and that the combined sparsity penalty is not
decomposable, we establish a new and sharp tail probability bound for
correlated heavy-tailed (sub-exponential) variables with an arbitrary
correlation structure, which is meritorious even on its own.

  * •

Finally, we prove that, compared to Zhang and Li (2022), the error rate of the
simultaneously estimated precision parameters improves by a factor of $p$, the
number of response variables. The improvement is remarkable for a large $p$,
as further corroborated in simulation studies. Code is available on GitHub,
and detailed proofs are given in the Supplementary Material.

##  2 Multi-task Learning for Gaussian Graphical Regressions

Denote by ${\boldsymbol{X}}=(X_{1},\ldots,X_{p})^{\top}$ the vector of
response variables (e.g., expression levels of $p$ genes) and
${\boldsymbol{U}}=(U_{1},\ldots,U_{q})^{\top}$ the vector of covariates (e.g.,
age, sex and genetic variants), and assume that

${\boldsymbol{X}}\,|\,{\boldsymbol{U}}={\boldsymbol{u}}\,\sim\,\mathcal{N}_{p}({\boldsymbol{\Gamma}}{\boldsymbol{u}},{\boldsymbol{\Sigma}}({\boldsymbol{u}})),$

where ${\boldsymbol{\Sigma}}({\boldsymbol{u}})$ is the conditional covariance,
and
${\boldsymbol{\Omega}}({\boldsymbol{u}})={\boldsymbol{\Sigma}}^{-1}({\boldsymbol{u}})$
is the conditional precision matrix linked to ${\boldsymbol{u}}$ via

$-{\boldsymbol{\Omega}}({\boldsymbol{u}})_{jk}=\begin{cases}-\sigma^{jj}&\quad
j=k,\\\ \beta^{\prime}_{jk0}+\sum_{h=1}^{q}\beta^{\prime}_{jkh}u_{h}&\quad
j\neq k,\\\ \end{cases}$ (1)

where $\beta^{\prime}_{jkh}=\beta^{\prime}_{kjh}$, $h\in\\{0\\}\cup[q]$,
$j,k\in[p]$. We assume
${\boldsymbol{\Omega}}({\boldsymbol{u}})_{jj}=\sigma^{jj}$ to be free of
${\boldsymbol{u}}$, and this is discussed after (2). Writing
${\boldsymbol{Z}}={\boldsymbol{X}}-{\boldsymbol{\Gamma}}{\boldsymbol{u}}=(Z_{1},\ldots,Z_{p})^{\top}$,
some straightforward algebra shows that (1) relates
${\boldsymbol{\Omega}}({\boldsymbol{u}})$ to the following regression, termed
Gaussian graphical regression (Zhang and Li, 2022):

$Z_{j}={\sum}_{k\neq j}^{p}\beta_{jk0}Z_{k}+{\sum}_{k\neq
j}^{p}{\sum}_{h=1}^{q}\underbrace{\beta_{jkh}\times
u_{h}Z_{k}}_{\text{{interaction term}}}+\epsilon_{j},$ (2)

where $\beta_{jkh}=\beta^{\prime}_{jkh}/\sigma^{jj}$, $\epsilon_{j}$ is
independent of ${\boldsymbol{Z}}_{-j}=\\{Z_{k}:k\neq j\in[p]\\}$ and
$\text{Var}(\epsilon_{j})=1/\sigma^{jj}$, for all $j,k$ and $h$. It is easy to
see that, when $\beta_{jkh}=0$ for all $j,k,h\neq 0$, (2) reduces to the usual
Gaussian graphical model with
${\boldsymbol{Z}}\,|\,{\boldsymbol{U}}={\boldsymbol{u}}\,\sim\,\mathcal{N}_{p}({\boldsymbol{0}},{\boldsymbol{\Sigma}})$.

Figure 1: An illustration of multi-task learning in Gaussian graphical
regression.

Model (2) provides a regression framework for estimating the precision
parameters in (1), by including the interactions between
${\boldsymbol{X}}_{-j}$ and ${\boldsymbol{u}}$. Correspondingly, the partial
correlation between $X_{j}$ and $X_{k}$, conditional on all other $X$
variables, is modeled as a function of ${\boldsymbol{u}}$; see an illustration
in Figure 1. The diagonal elements of
${\boldsymbol{\Omega}}({\boldsymbol{u}})$ (i.e., $\sigma^{jj}$’s) are
connected to the residual variances in (2), that is,
$\text{Var}(\epsilon_{j})=1/\sigma^{jj}$. From this perspective, assuming
$\sigma^{jj}$ to be free of ${\boldsymbol{u}}$ may be viewed as assuming the
residual variance of $Z_{j}$, after removing the effects of
${\boldsymbol{u}}$, ${\boldsymbol{Z}}_{-j}$ and their interactions, not to
depend on ${\boldsymbol{u}}$, which is plausible in the context of regression.

Let
${\boldsymbol{\beta}}_{j}=(\beta_{j10},\ldots,\beta_{j,j-1,0},\beta_{j,j+1,0},\ldots,\beta_{jp0},\ldots,\beta_{j1q},\ldots,\beta_{j,j-1,q},\beta_{j,j+1,q},\ldots,\beta_{jpq})^{\top}\in\mathbb{R}^{(p-1)(q+1)}$
be the vector of coefficients in (2) and denote
${\boldsymbol{\beta}}=({\boldsymbol{\beta}}_{1},\ldots,{\boldsymbol{\beta}}_{p})$.
To expose the key ideas, we assume a known ${\boldsymbol{\Gamma}}$ in the
ensuing development, and focus on the estimation of ${\boldsymbol{\beta}}$;
extensions with unknown ${\boldsymbol{\Gamma}}$ are possible, but with more
involved notation. With $n$ independent data
$\mathcal{D}=\\{({\boldsymbol{u}}^{(i)},{\boldsymbol{x}}^{(i)}),i\in[n]\\}$,
let
${\boldsymbol{z}}^{(i)}={\boldsymbol{x}}^{(i)}-{\boldsymbol{\Gamma}}{\boldsymbol{u}}^{(i)}$
and
${\boldsymbol{w}}^{(i)}_{-j}={\boldsymbol{z}}_{-j}^{(i)}\otimes{\boldsymbol{u}}^{(i)}$,
where $\otimes$ denotes the Kronecker product. To estimate
${\boldsymbol{\beta}}$, we consider

$\arg{\min}_{{\boldsymbol{\beta}}}\frac{1}{2n}\sum_{j=1}^{p}\sum_{i=1}^{n}(z^{(i)}_{j}-{\boldsymbol{w}}^{(i)}_{-j}{}^{\top}{\boldsymbol{\beta}}_{j})^{2}+\lambda_{1}\sum_{h=1}^{q}\|{\boldsymbol{b}}_{h}\|_{2}+\lambda_{2}\|{\boldsymbol{\beta}}\|_{1},$
(3)

where
${\boldsymbol{b}}_{h}=(({\boldsymbol{\beta}}_{1})_{(h)},\ldots,({\boldsymbol{\beta}}_{p})_{(h)})$
with
$({\boldsymbol{\beta}}_{j})_{(h)}=(\beta_{j1h},\ldots,\beta_{j,j-1,h},\beta_{j,j+1,h},\ldots,\beta_{jph})$
and $\lambda_{1},\lambda_{2}\geq 0$ are tuning parameters. The convex
regularizing terms, $\sum_{h=1}^{q}\|{\boldsymbol{b}}_{h}\|_{2}$ and
$\|{\boldsymbol{\beta}}\|_{1}$, encourage group- and element-wise sparsity,
respectively, though the group sparse penalty is not applied to
$({\boldsymbol{\beta}}_{j})_{(0)}$ (i.e., intercept terms). The combined
sparsity penalty in (3) is termed the sparse group lasso penalty (Simon et
al., 2013). Distinguishing from the node-wise estimation (Zhang and Li, 2022),
(3) considers $p$ graphical regressions simultaneously, with the group lasso
penalty regulating ${\boldsymbol{b}}_{h}$ (the effects from the $h$th
covariate) cross tasks; see Figure 1. Thus, (3) is termed multi-task learning
for a Gaussian graphical regression model.

##  3 An Augmented Lagrangian Algorithm with a Semi-smooth Newton Method

By specifying a convex objective function, (3) facilitates use of existing
algorithms for solving sparse group lasso problems, including proximal
gradient descent (Chen et al., 2012), block coordinate descent (Simon et al.,
2013; Wang et al., 2019) and alternating direction method of multipliers
(ADMM) (Boyd et al., 2011). We propose to use an efficient augmented
Lagrangian algorithm which uses a semi-smooth Newton method to solve
subproblems. We opt for the semi-smooth Newton method because it exploits the
second order sparsity structure in our setting (Zhang et al., 2020) and enjoys
a linear rate of convergence.

More specifically, write
$\mathcal{W}=\begin{pmatrix}{\boldsymbol{W}}_{-1}&\cdots&{\boldsymbol{0}}_{n\times(p-1)(q+1)}\\\
\vdots&\ddots&\vdots\\\
{\boldsymbol{0}}_{n\times(p-1)(q+1)}&\cdots&{\boldsymbol{W}}_{-p}\end{pmatrix}$,
where
${\boldsymbol{W}}_{-j}=[{\boldsymbol{w}}^{(1)}_{-j};\ldots;{\boldsymbol{w}}^{(n)}_{-j}]$,
and
${\boldsymbol{a}}=(z^{(1)}_{1},\ldots,z^{(n)}_{1},z^{(2)}_{1},\ldots,z^{(n)}_{p})\in\mathbb{R}^{np}$.
Correspondingly, (3) can be written as

$\arg{\min}_{{\boldsymbol{\beta}}}\frac{1}{2n}\|\mathcal{W}{\boldsymbol{\beta}}-{\boldsymbol{a}}\|_{2}^{2}+\lambda_{1}\sum_{h=1}^{q}\|{\boldsymbol{b}}_{h}\|_{2}+\lambda_{2}\|{\boldsymbol{\beta}}\|_{1},$

and its dual problem can be written as:

$\displaystyle\min_{{\boldsymbol{y}},{\boldsymbol{g}}}\,\,\langle{\boldsymbol{a}},{\boldsymbol{y}}\rangle+\frac{1}{2}\|{\boldsymbol{y}}\|^{2}+h^{*}({\boldsymbol{g}})$
$\displaystyle\text{s.t.}\,\,\mathcal{W}^{\top}{\boldsymbol{y}}+{\boldsymbol{g}}=0,$

where $h^{*}(\cdot)$ is the conjugate function of
$h({\boldsymbol{\beta}}):=\lambda_{1}\sum_{h=1}^{q}\|{\boldsymbol{b}}_{h}\|_{2}+\lambda_{2}\|{\boldsymbol{\beta}}\|_{1}$,
i.e.,
$h^{*}({\boldsymbol{\beta}})=\sup_{\boldsymbol{e}}\\{\langle{\boldsymbol{\beta}},\boldsymbol{e}\rangle-h(\boldsymbol{e})\\}$.
Define the proximal mapping of $h(\cdot)$ to be

$\text{Prox}_{h}({\boldsymbol{u}})=\arg\min_{\boldsymbol{x}}\left\\{h({\boldsymbol{x}})+\frac{1}{2}\|{\boldsymbol{x}}-{\boldsymbol{u}}\|^{2}\right\\}.$

Similar to ADMM, the augmented Lagrangian of the dual problem is given by

$\langle{\boldsymbol{a}},{\boldsymbol{y}}\rangle+\frac{1}{2}\|{\boldsymbol{y}}\|^{2}+h^{*}({\boldsymbol{g}})+\frac{\tau}{2}\|\mathcal{W}^{\top}{\boldsymbol{y}}+{\boldsymbol{g}}-\sigma^{-1}{\boldsymbol{\beta}}\|^{2}-\frac{1}{2\tau}\|{\boldsymbol{\beta}}\|^{2},$

where $\tau>0$, and the subproblem relates to iteratively minimizing

$\psi_{k}({\boldsymbol{y}})=\langle{\boldsymbol{a}},{\boldsymbol{y}}\rangle+\frac{1}{2}\|{\boldsymbol{y}}\|_{2}^{2}+h^{*}(\text{Prox}_{h^{*}/\tau_{k}}(\tau_{k}^{-1}{\boldsymbol{\beta}}^{k}-\mathcal{W}^{\top}{\boldsymbol{y}})+\frac{\tau_{k}}{2}\|\text{Prox}_{h}(\tau_{k}^{-1}{\boldsymbol{\beta}}^{k}-\mathcal{W}^{\top}{\boldsymbol{y}})\|_{2}^{2}-\frac{1}{2\tau_{k}}\|{\boldsymbol{\beta}}^{k}\|^{2},$

where ${\boldsymbol{\beta}}^{k}$ is the updated estimate and $\tau_{k}>0$ is
given at the $k$th iteration; see the detailed steps below.

Algorithm 1 An augmented Lagrangian method for solving (3)

Let $\tau_{0}>0$ be given, and choose
$({\boldsymbol{y}}_{0},{\boldsymbol{g}}_{0},{\boldsymbol{\beta}}_{0})$.
Iterate the following steps for $k=0,1,\ldots,$ until convergence.

Step 1: compute
${\boldsymbol{y}}^{k+1}=\min_{\boldsymbol{y}}\psi_{k}({\boldsymbol{y}})$ and
${\boldsymbol{g}}^{k+1}=\text{Prox}_{p^{*}/\tau_{k}}(\tau_{k}^{-1}{\boldsymbol{\beta}}^{k}-\mathcal{W}^{\top}{\boldsymbol{y}}^{k+1})$;

Step 2: compute
${\boldsymbol{\beta}}^{k+1}={\boldsymbol{\beta}}^{k}-\tau_{k}(\mathcal{W}^{\top}{\boldsymbol{y}}^{k+1}+{\boldsymbol{g}}^{k+1})$;

Step 3: update $\tau_{k+1}$ so that $\tau_{k+1}>\tau_{k}$.

Step 1 solves the subproblem $\min_{\boldsymbol{y}}\psi_{k}({\boldsymbol{y}})$
with a semi-smooth Newton method, utilizing a generalized Jacobian of
$\text{Prox}_{h}(\cdot)$ given in Zhang et al. (2020). In general, the
algorithm outperforms the first-ordered ADMM algorithm in computational
efficiency and estimation accuracy, because it accounts for the second order
sparsity structure when solving
$\min_{\boldsymbol{y}}\psi_{k}({\boldsymbol{y}})$ (Li et al., 2017; Zhang et
al., 2020).

Positive-definiteness. A natural sufficient condition to ensure positive
definiteness of ${\boldsymbol{\Omega}}({\boldsymbol{u}})$ is diagonal
dominance, leading to
$\max(1,\|{\boldsymbol{u}}\|_{\infty})\|{\boldsymbol{\beta}}_{j}\|_{1}<1$.
With Assumption 1 stipulating $|u^{(i)}_{h}|\leq M$, positive definiteness
holds when $\|{\boldsymbol{\beta}}_{j}\|_{1}<1/\max(1,M)$, $j\in[p]$. Assuming
$u_{h}\in[-1,1]$ (if not, rescale first), this sufficient condition can be
simplified to $\|{\boldsymbol{\beta}}_{j}\|_{1}<1$ for all $j$ (note that
${\boldsymbol{\beta}}_{j}$ is sparse), suggesting that, to satisfy diagonal
dominance, the “effect sizes" of ${\boldsymbol{u}}$ (i.e.,
$\|{\boldsymbol{\beta}}_{j}\|_{1}$) on partial correlations cannot be too
large. If the true covariance/precision matrix is positive definite, Theorem 2
implies that the estimated precision matrix is asymptotically positive
definite. For finite sample cases, a post-hoc rescaling procedure can be
utilized for securing positive definiteness.

Tuning. The $\lambda_{1}$ and $\lambda_{2}$ in (3) can be jointly selected via
$L$-fold (e.g., $L=5$) cross validation. To facilitate parallelization of the
algorithm, we consider a hierarchical selection procedure. Rewrite
$\lambda_{1}=(1-\alpha)\lambda_{0}$ and $\lambda_{2}=\alpha\lambda_{0}$, where
$\alpha$ reflects the weight of the lasso penalty relative to the group lasso
penalty (e.g., $\alpha=0$ and $1$ correspond to group lasso and lasso,
respectively) and $\lambda_{0}$ reflects the total amount of regularization.
We assess a set of values for $\alpha\in[0,1]$; for each $\alpha$, a grid of
$\lambda_{0}$ values are considered in cross validation. The search can be
parallelized for different $\alpha$’s.

##  4 Theory: Concentration Inequality and Error Rate for Multi-task Learning

We derive the non-asymptotic $\ell_{2}$ error rate of the sparse group lasso
estimator from (3). The main challenge is that the $p$ tasks in (3) are
correlated, and this differs from the usual multi-task learning with group
sparsity (Lounici et al., 2011). To tackle this challenge, we have made a key
advance in Theorem 1 that gives a new tail bound for the sum of correlated
heavy-tailed (sub-exponential) variables with an arbitrary correlation
structure.

Our theoretical investigation faces other challenges. First, because the
design matrix in (3) includes high-dimensional interactions between
${\boldsymbol{z}}^{(i)}$ and ${\boldsymbol{u}}^{(i)}$, and the variance of
${\boldsymbol{z}}^{(i)}$ is a function of ${\boldsymbol{u}}^{(i)}$,
characterizing the joint distribution of each row in ${\boldsymbol{W}}_{-j}$
is difficult and requires a delicate treatment. Second, as the combined
penalty term
$\lambda_{1}\sum_{h=1}^{q}\|{\boldsymbol{b}}_{h}\|_{2}+\lambda_{2}\|{\boldsymbol{\beta}}\|_{1}$
is not decomposable, the classic techniques for decomposable regularizers and
null space properties (Negahban et al., 2012) are not applicable. By utilizing
a novel concentration inequality for the sum of correlated variables in (3),
we derive two interrelated bounds for the stochastic term, whose combination
yields a sharp upper bound of the stochastic term. We, therefore, show that
our proposed estimator possesses an improved $\ell_{2}$ error bounds compared
to the lasso and the group lasso when the true coefficients are simultaneously
sparse, and, more importantly, the error rate of the multi-task estimates
improves by a factor of $p$, compared to the separate node-wise estimates
obtained by Zhang and Li (2022).

###  4.1 A new concentration inequality for the sum of correlated variables

The concentration results for correlated random variables are often derived
under specific correlation structures, such as weak dependence (Merlevède et
al., 2011) and asymptotic independence (Ko and Tang, 2008). These structures
are unlikely to be applicable to our setting because the error terms in (3),
$\epsilon_{1},\ldots,\epsilon_{p}$ (from the $p$ tasks), depend on, for
example, the expressions of $p$ genes, and are correlated via a complicated
co-expression network. To bound $\sum_{j}\epsilon_{j}$, we employ a novel idea
that partitions the index set of $p$ response variables into mutually
exclusive subsets such that the variables within each subset are independent.
This can be done by exploring the topology of a graph and solving a vertex
coloring problem (Lewis, 2015); see the proof of Theorem 1 in the
Supplementary Material. With that, we present a concentration inequality
result under a general correlation structure.

######  Theorem 1.

Consider $N$ correlated mean zero sub-exponential random variables $Y_{j}$,
$j\in[N]$ and an induced network $G(V,E)$ with a node set $V=\\{1,\ldots,N\\}$
and an edge set $E=\\{(j,k):\text{Cov}(Y_{j},Y_{k})\neq 0\\}$. Denote the
maximum node degree of $G(V,E)$ by $d_{\max}$ and let
$c_{G}=\min\left(d_{\max}+1,\frac{1+\sqrt{8|E|+1}}{2}\right)$. For any $t\geq
0$ and a constant $c>0$, it holds that

$\mathbb{P}\left(\sum_{j=1}^{N}Y_{j}\geq t\right)\leq
c_{G}\exp\left[-c\min\left\\{\frac{t^{2}}{c_{G}^{2}\sum_{j=1}^{N}\|Y_{i}\|_{\psi_{1}}^{2}},\frac{t}{c_{G}\max_{j}\|Y_{i}\|_{\psi_{1}}}\right\\}\right],$

where $\|\cdot\|_{\psi_{1}}$ is the sub-exponential norm defined in the
Supplementary Material. The results are inclusive. For example, when
$Y_{1},\ldots,Y_{N}$ are mutually independent, we have $c_{G}=1$, and the
inequality reduces to the usual form for independent sub-exponential random
variables; if the variables are correlated with $c_{G}=O(1)$, we obtain a tail
probability in the same order as when the variables are independent. The
theorem leads to a corollary that gives a sharp bound on the sum of correlated
chi-squared random variables, which is critical for our proof.

######  Corollary 1.

Consider $N$ correlated $\chi^{2}$ variables $Y_{j}\sim\chi^{2}_{d_{j}}$,
$j\in[N]$ and an induced network $G(V,E)$ defined as in Theorem 1. Denote the
maximum node degree of $G(V,E)$ by $d_{\max}$ and let
$c_{G}=\min\left(d_{\max}+1,\frac{1+\sqrt{8|E|+1}}{2}\right)$. For any $t\geq
0$, it holds that

$\mathbb{P}\left(\sum_{j=1}^{N}Y_{j}-\sum_{j=1}^{N}d_{j}\geq t\right)\leq
c_{G}\exp\left[-\frac{\left\\{t-(c_{G}-1)\sum_{j=1}^{N}d_{j}\right\\}^{2}}{4c_{G}(t+\sum_{j=1}^{N}d_{j})}\right].$

###  4.2 Error rate analysis

Let $\mathcal{S}$ be the element-wise support set of ${\boldsymbol{\beta}}$
and $\mathcal{G}$ be the group-wise support set of ${\boldsymbol{\beta}}$, and
denote by $s_{e}=|\mathcal{S}|$ and $s_{g}=|\mathcal{G}|$, i.e., $s_{e}$ and
$s_{g}$ are the numbers of nonzero entries and nonzero groups, respectively.
Without loss of generality, we assume $\sigma^{jj}=1$. We state the needed
regularity conditions.

######  Assumption 1.

Suppose ${\boldsymbol{u}}^{(i)}$ are i.i.d. mean zero random vectors with a
covariance matrix satisfying
$\lambda_{\min}(\text{Cov}({\boldsymbol{u}}^{(i)}))\geq 1/\phi_{0}$ for some
constant $\phi_{0}>0$. Moreover, there exists a constant $M>0$ such that
$|u^{(i)}_{h}|\leq M$ for all $i$ and $h$.

######  Assumption 2.

Suppose
$\phi_{1}\leq\lambda_{min}(\text{Cov}({\boldsymbol{z}}^{(i)}))\leq\lambda_{max}(\text{Cov}({\boldsymbol{z}}^{(i)}))\leq\phi_{2}$
for some constants $\phi_{1},\phi_{2}>0$.

######  Assumption 3.

The dimensions $p,q$ and sparsity $s_{e}$ satisfy
$\hbox{log}\,p+\hbox{log}\,q=\mathcal{O}(n^{\delta})$ and
$s_{e}=o(n^{\delta})$ for $\delta\in[0,1/6]$. The maximum column $\ell_{0}$
norm of ${\boldsymbol{\Omega}}({\boldsymbol{u}})$ is bounded above by a
positive constant $d_{0}$.

Assumption 1 stipulates that the covariates are element-wise bounded, which is
needed in characterizing the joint distribution of each row in $\mathcal{W}$.
This condition is not restrictive as genetic variants are often coded to be
$\\{0,1\\}$ or $\\{0,1,2\\}$ (Chen et al., 2016). Assumptions 1 and 2 impose
bounded eigenvalues on $\text{Cov}({\boldsymbol{u}}^{(i)})$ and
$\text{Cov}({\boldsymbol{z}}^{(i)})$ as commonly done in the high-dimensional
regression literature (Chen et al., 2016; Hao et al., 2018; Cai et al., 2019).
Assumption 3 is a sparsity condition, allowing both $\hbox{log}\,p$ and
$\hbox{log}\,q$ to grow at a polynomial order of $n$. Moreover, the number of
nonzero entries $s_{e}$ can grow with $n$. This condition and
$\delta\in[0,1/6]$ are useful when establishing a restricted eigenvalue
condition (Bickel et al., 2009) for $\mathcal{W}^{\top}\mathcal{W}/n$ and when
bounding the stochastic term
$\langle{\boldsymbol{\epsilon}},\mathcal{W}\Delta\rangle$. The condition on
the column $\ell_{0}$ norm of ${\boldsymbol{\Omega}}({\boldsymbol{u}})$
implies that, for example, each gene is only connected to a finite number of
genes.

######  Theorem 2.

Let $s_{\lambda}$ denote the number of nonzero entries in a candidate model
such that $s_{e}<s_{\lambda}\leq n$. Suppose Assumptions 1-3 hold,
$s_{\lambda}(\hbox{log}\,p+\hbox{log}\,q)=\mathcal{O}(\sqrt{n})$ and $n\geq
A_{1}\\{s_{g}\hbox{log}(eq/s_{g})+s_{e}\hbox{log}(ep)\\}$ for some constant
$A_{1}>0$. Then the sparse group lasso estimator $\hat{\boldsymbol{\beta}}$ in
(3) with

$\lambda_{1}=C\sqrt{\hbox{log}(eq/s_{g})/n+2s_{e}\hbox{log}(ep)/(ns_{g})},\quad\lambda_{2}=\sqrt{s_{g}/s_{e}}\lambda_{1}$
(4)

satisfies, with probability at least
$1-C_{1}\exp[-C_{2}\\{s_{g}\hbox{log}(eq/s_{g})+s_{e}\hbox{log}(ep)\\}]$,

$\|\hat{\boldsymbol{\beta}}-{\boldsymbol{\beta}}\|_{2}^{2}\precsim\frac{1}{n}\\{s_{g}\hbox{log}(eq/s_{g})+s_{e}\hbox{log}(ep)\\},$
(5)

where $C$, $C_{1}$, and $C_{2}$ are positive constants.

The separate node-wise regressions considered in Zhang and Li (2022) yield an
error rate of
$\frac{1}{n}\\{ps_{g}\hbox{log}(eq/s_{g})+s_{e}\hbox{log}(ep)\\}$, which is
slower than that in (5) by a factor $p$ if
$s_{e}\hbox{log}p=O(s_{g}\hbox{log}q)$; as $p$ often far exceeds $n$, the
improvement with multi-task learning is considerable.

Theorem 2 also shows that our proposed estimator enjoys an improved $\ell_{2}$
error bound over estimators with only a lasso or a group lasso penalty on
${\boldsymbol{\beta}}$. Specifically, given that the dimension of
${\boldsymbol{\beta}}$ is $p(p-1)(q+1)$ and $s_{g}\leq s_{e}$, applying the
regular lasso regularizer $\lambda\|{\boldsymbol{\beta}}\|_{1}$ alone would
yield an error bound of $(s_{e}/n)\hbox{log}(pq)$ (Negahban et al., 2012),
which is slower than that in (5) when $\hbox{log}p/\hbox{log}q=o(1)$ and
$s_{g}/s_{e}=o(1)$, corresponding to group sparsity. Moreover, when $p>n+1$,
estimating with the group lasso regularizer
$\lambda_{1}\sum_{h}\|{\boldsymbol{b}}_{h}\|_{2}$ alone, which excludes
$({\boldsymbol{\beta}}_{j})_{(0)}$, is not feasible, because the dimension of
the latter (i.e., $p(p-1)$) exceeds $n$. If we utilize a group lasso
regularizer $\lambda_{g}\|{\boldsymbol{\beta}}_{j}\|_{1,2}$ that includes
${\boldsymbol{b}}_{0}$, the estimator would have an $\ell_{2}$ error bound of
$(s_{g}/n)\hbox{log}\,q+(s_{g}/n)p(p-1)$ (Lounici et al., 2011), which is
slower than that in (5) when $\hbox{log}q/\\{p(p-1)\\}=o(1)$ and
$s_{j}/s_{j,g}=o(p(p-1)/\hbox{log}p)$, corresponding to within-group sparsity.
In Theorem 2, the condition
$s_{\lambda}(\hbox{log}\,p+\hbox{log}\,q)=\mathcal{O}(\sqrt{n})$ upper bounds
the size of candidate models, which helps to bound the stochastic errors from
$p$ tasks.

##  5 Numerical Experiments

We compare the finite sample performance of our proposed method, defined in
(3) (referred to as “MtRegGMM"), with those of three competing solutions,
namely, a benchmarking i.i.d Gaussian graphical model estimated by the
neighborhood selection method (Meinshausen and Bühlmann, 2006) (“IID"), a
lasso estimator
$\arg{\min}_{{\boldsymbol{\beta}}}\frac{1}{2n}\sum_{j=1}^{p}\sum_{i=1}^{n}(z^{(i)}_{j}-{\boldsymbol{W}}^{(i)}_{-j}{}^{\top}{\boldsymbol{\beta}}_{j})^{2}+\lambda\|{\boldsymbol{\beta}}\|_{1}$
(Joint${}_{\text{lasso}}$), and the separate regressions considered in Zhang
and Li (2022) (“RegGMM").

We simulate $n$ samples
$\\{({\boldsymbol{u}}^{(i)},{\boldsymbol{x}}^{(i)}),i\in[n]\\}$ from (1), with
${\boldsymbol{x}}^{(i)}\in\mathbb{R}^{p}$ (e.g., genes) and external covariate
${\boldsymbol{u}}^{(i)}\in\mathbb{R}^{q}$ (e.g., SNPs). The
${\boldsymbol{u}}^{(i)}_{j}$’s are generated i.i.d. from Bernoulli$(0.5)$. We
randomly select 3 covariates to have nonzero effects, and in the graphs for
these covariates and the population-level graph, five edges are randomly
selected to be nonzero. We set $\sigma^{jj}=1$ for $j\in[p]$. The initial
nonzero coefficients $\beta_{jkh}$ are set to 0.3. For each $j$, we rescale
$\\{\beta_{jkh}\\}_{k\neq j\in[p],h\in\\{0\\}\cup[q]}$ by dividing each entry
by $\sum_{k\neq j\in[p],h\in\\{0\\}\cup[q]}|\beta_{jkh}|$. After rescaling,
for each $j,k$ and $h$, we use the average of $\beta_{jkh}$ and $\beta_{kjh}$
to fill the entries at $jkh$ and $kjh$. This process results in symmetry with
diagonal dominance and, thus, ensures positive definiteness of the precision
matrices. For each simulation configuration, we generate 50 independent data
sets; given ${\boldsymbol{u}}^{(i)}$, we determine
${\boldsymbol{\Omega}}({\boldsymbol{u}}^{(i)})$ and
${\boldsymbol{\Sigma}}({\boldsymbol{u}}^{(i)})$, and generate the $i$th sample
${\boldsymbol{x}}^{(i)}$ from
$\mathcal{N}({\boldsymbol{\Gamma}}{\boldsymbol{u}}^{(i)},{\boldsymbol{\Sigma}}({\boldsymbol{u}}^{(i)}))$,
$i\in[n]$. For a fair comparison, tuning parameters in all of the methods are
selected via 5-fold cross validation.

Table 1: Estimation accuracy of ${\boldsymbol{\beta}}_{j}$ and ${\boldsymbol{\Omega}}_{j}$ in simulations with varying sample size network size $p$ and covariate dimension $q$.  $n$ |  $p$, $q$ | Method | TPRβ | FPRβ | Error of ${\boldsymbol{\beta}}$ | Error of ${\boldsymbol{\Omega}}$
---|---|---|---|---|---|---
100 |  $p=20$ $q=100$ | MtRegGMM |  0.985 (0.013) |  0.0001 (0.0000) |  0.534 (0.032) |  0.405 (0.046)
Joint${}_{\text{lasso}}$ | 0.942 (0.015) | 0.0003 (0.0000) | 0.762 (0.033) | 0.712 (0.067)
RegGMM | 0.922 (0.020) | 0.0013 (0.0001) | 1.146 (0.045) | 1.984 (0.153)
IID | - | - | - | 1.764 (0.034)
$p=20$ $q=200$ | MtRegGMM |  0.983 (0.021) |  0.0001 (0.0000) |  0.549 (0.037) |  0.436 (0.060)
Joint${}_{\text{lasso}}$ | 0.918 (0.030) | 0.0002 (0.0000) | 0.848 (0.040) | 0.912 (0.088)
RegGMM | 0.923 (0.016) | 0.0009 (0.0001) | 1.207 (0.036) | 2.223 (0.199)
IID | - | - | - | 1.777 (0.026)
$p=100$ $q=200$ | MtRegGMM |  0.964 (0.013) |  0.0000 (0.0000) | 0.612 (0.034) |  0.489 (0.055)
Joint${}_{\text{lasso}}$ | 0.896 (0.020) |  0.0000 (0.0000) | 0.846 (0.038) | 0.786 (0.069)
RegGMM | 0.918 (0.020) | 0.0003 (0.0000) | 2.620 (0.060) | 5.845 (0.329)
IID | - | - | - | 5.119 (0.890)
Figure 2: Estimation errors and computing time.

To evaluate the estimation accuracy, we report the estimation errors
$\|{\hat{\boldsymbol{\beta}}-{\boldsymbol{\beta}}}\|_{2}$, where
$\hat{\boldsymbol{\beta}}$’s, with a slight overuse of notation, denote the
estimates of ${\boldsymbol{\beta}}$’s obtained by various methods. For the
selection accuracy, we report the true positive rate (TPR) and false positive
rate (FPR). Also reported is the average estimation error of the precision
matrix defined to be
$\sum_{i=1}^{n}\|\hat{{\boldsymbol{\Omega}}}_{i}-{\boldsymbol{\Omega}}_{i}\|_{F,\text{off}}^{2}/n$,
where ${\boldsymbol{\Omega}}_{i}$ is the true subject-specific precision
matrix and $\hat{\boldsymbol{\Omega}}_{i}$ is estimated by a given method.

The proposed MtRegGMM outperforms the competing methods in estimation and
selection accuracy for various $p$ and $q$. The estimation errors of MtRegGMM
increases with $p$ and $q$, confirming the results of Theorem 2. In the left
plot of Figure 2, we compare the estimation errors from RegGMM and MtRegGMM as
$p$ increases. The error from RegGMM increases roughly linearly with $p$ while
that from MtRegGMM remains relatively stable as $p$ increases, again
confirming the results of Theorem 2.

Finally, we assess the computation cost including tuning. The middle and right
panels of Figure 2 show the computation time of MtRegGMM for a given
$(\lambda_{1},\lambda_{2})$. The simulations were run on an iMac with a 3.6
GHz Intel Core i9 processor. Because the number of parameters is
$\mathcal{O}(p^{2}q)$, the total computing cost is expected to be in the order
of $p^{2}q$, as seen in Figure 2.

##  6 Co-expression QTL Analysis

We analyze the REMBRANDT trial (GSE108476) with a subcohort of $n=178$
glioblastoma multiforme (GBM) patients, who had undergone microarray and
single nucleotide polymorphism (SNP) chip profiling, with both gene expression
and SNP data available for analysis. The raw data were pre-processed and
normalized using standard pipelines (Gusev et al., 2018).

The response variables are the expression levels of $p=73$ genes belonging to
the human glioma pathway in the Kyoto Encyclopedia of Genes and Genomes (KEGG)
database (Kanehisa and Goto, 2000). The covariates include local SNPs (i.e.,
SNPs that fall within 2kb upstream and 0.5kb downstream of the gene) residing
near these 73 genes, resulting in a total of 118 SNPs. SNPs are coded with “0"
indicating homozygous in the major allele and “1" otherwise. Age and gender
are also included in analysis. Consequently, there are $q=120$ covariates,
bringing a total of $73\times 36\times 121=317,988$ parameters (including
intercepts). We construct the population-level gene co-expression network, and
examine if and how age, gender and SNPs modulate the network.

Figure 3: The population-level gene co-expression network, where the size of
each node is proportional to the mean expression level, and edges with
positive (negative) partial correlations are shown in red dashed (black solid)
lines. Figure 4: Effects of each covariate (e.g., SNPs) on edges or partial
correlations, with positive (negative) effects shown in red dashed (black
solid) lines.

We have used the proposed method to construct the population level network as
shown in Figure 3. Most of the connected genes in Figure 3 are known to be
oncogenes. For example, PIK3CA is a protein coding gene and is one of the most
highly mutated oncogenes identified in human cancers (Samuels and Velculescu,
2004). The PIK3CA gene is a part of the PI3K/AKT/MTOR signaling pathway, which
is one of the core pathways in human GBM and other types of cancer (Network et
al., 2008). In Figure 3, we can identify several core pathways in human GBM
including the PI3K/ AKT/MTOR, Ras-Raf-MEK-ERK and calcium signaling pathways.

We next examine the covariate effects on the network. Identified by our method
are nine co-expression QTLs, namely, rs6701524, rs10519201, rs1347069,
rs9303511, rs503314, rs7286558, rs759950, rs306098, rs25919. The network
effects of rs6701524 are shown in Figure 4 (the left panel). This SNP,
residing in MTOR, is found to affect co-expressions of genes in the PI3K/
AKT/MTOR pathway. This is an interesting finding as PI3K/MTOR is a key pathway
in GBM development and progression, and inhibition of PI3K/MTOR signaling was
found effective in increasing survival with GBM tumor (Batsios et al., 2019).
This co-expression QTL can potentially play an important role in activating
the PI3K/MTOR pathway. Shown in Figure 4 (the right panel) are the network
effects of rs1347069, a variant of MAP2K1. The figure indicates that this SNP
mostly affects the co-expressions of genes in the Ras-Raf-MEK-ERK pathway.

Moreover, we have found some novel co-expression QTLs. For example, rs10519201
regulates the co-expressions of PDGFB and GADD45A; rs9303511 is associated
with the co-expressions of GADD45A and CALML4; rs503314 influences the co-
expressions of CCND1 and PLCG2, CALML3; rs7286558 may modify the co-
expressions of EGFR and PRKCA; rs759950 modulates the co-expressions of MAP2K1
and HRAS; rs306098 influences the co-expressions of PIK3CD and CAMK2D; rs25919
may alter the co-expressions of CAMK2A and GADD45G. Co-expression QTL
identification has recently sparked much interest, and these findings warrant
in-depth investigations.

##  7 Conclusion

We have proposed a computationally efficient multi-task learning estimator for
Gaussian graphical regression models, and proved that the error rate of our
estimator improves upon that of the separate node-wise lasso estimator by a
factor of $p$ (i.e., the dimension of a graph). The remarkable improvement is
possible because multi-task learning borrows information across tasks.
Moreover, unlike the usual multi-task problems, our tasks are entangled in a
complicated correlation structure, defying existing theoretical results; we
address it by establishing a new tail probability bound for correlated heavy-
tailed variables with an arbitrary correlation structure. Simulations and a
gene network application have suggested feasibility and utility of the
proposal in high dimensional settings, where both the response variables and
covariates exceed the sample size. Future work includes investigations of the
optimal error rates and the conditions needed to ensure variable selection
consistency as well as extensions to accommodate more general types of
responses, such as binary data.

## References

  * Batsios et al. (2019) Batsios, G., Viswanath, P., Subramani, E., Najac, C., Gillespie, A. M., Santos, R. D., Molloy, A. R., Pieper, R. O., and Ronen, S. M. (2019), “PI3K/mTOR inhibition of IDH1 mutant glioma leads to reduced 2HG production that is associated with increased survival,” Scientific Reports, 9, 1–15.
  * Bickel et al. (2009) Bickel, P. J., Ritov, Y., and Tsybakov, A. B. (2009), “Simultaneous analysis of Lasso and Dantzig selector,” The Annals of Statistics, 37, 1705–1732.
  * Boyd et al. (2011) Boyd, S., Parikh, N., Chu, E., Peleato, B., Eckstein, J., et al. (2011), “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends® in Machine learning, 3, 1–122.
  * Brynedal et al. (2017) Brynedal, B., Choi, J., Raj, T., Bjornson, R., Stranger, B. E., Neale, B. M., Voight, B. F., and Cotsapas, C. (2017), “Large-scale trans-eQTLs affect hundreds of transcripts and mediate patterns of transcriptional co-regulation,” The American Journal of Human Genetics, 100, 581–591.
  * Cai et al. (2012) Cai, T. T., Li, H., Liu, W., and Xie, J. (2012), “Covariate-adjusted precision matrix estimation with an application in genetical genomics,” Biometrika, 100, 139–156.
  * Cai et al. (2019) Cai, T. T., Zhang, A., and Zhou, Y. (2019), “Sparse Group Lasso: Optimal Sample Complexity, Convergence Rate, and Statistical Inference,” in arXiv preprint arXiv:1909.09851.
  * Chen et al. (2016) Chen, M., Ren, Z., Zhao, H., and Zhou, H. (2016), “Asymptotically normal and efficient estimation of covariate-adjusted Gaussian graphical model,” Journal of the American Statistical Association, 111, 394–406.
  * Chen et al. (2012) Chen, X., Lin, Q., Kim, S., Carbonell, J. G., and Xing, E. P. (2012), “Smoothing proximal gradient method for general structured sparse regression,” The Annals of Applied Statistics, 6, 719–752.
  * Cheng et al. (2014) Cheng, J., Levina, E., Wang, P., and Zhu, J. (2014), “A sparse Ising model with covariates,” Biometrics, 70, 943–953.
  * Danaher et al. (2014) Danaher, P., Wang, P., and Witten, D. M. (2014), “The joint graphical lasso for inverse covariance estimation across multiple classes,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76, 373–397.
  * Guo et al. (2011) Guo, J., Levina, E., Michailidis, G., and Zhu, J. (2011), “Joint estimation of multiple graphical models,” Biometrika, 98, 1–15.
  * Gusev et al. (2018) Gusev, Y., Bhuvaneshwar, K., Song, L., Zenklusen, J.-C., Fine, H., and Madhavan, S. (2018), “The REMBRANDT study, a large collection of genomic data from brain cancer patients,” Scientific Data, 5, 180158\.
  * Hao et al. (2018) Hao, N., Feng, Y., and Zhang, H. H. (2018), “Model selection for high-dimensional quadratic regression via regularization,” Journal of the American Statistical Association, 113, 615–625.
  * Kanehisa and Goto (2000) Kanehisa, M. and Goto, S. (2000), “KEGG: kyoto encyclopedia of genes and genomes,” Nucleic Acids Research, 28, 27–30.
  * Ko and Tang (2008) Ko, B. and Tang, Q. (2008), “Sums of dependent nonnegative random variables with subexponential tails,” Journal of Applied Probability, 45, 85–94.
  * Kolar et al. (2010) Kolar, M., Parikh, A. P., and Xing, E. P. (2010), “On sparse nonparametric conditional covariance selection,” in Proceedings of the 27th International Conference on International Conference on Machine Learning, pp. 559–566.
  * Lewis (2015) Lewis, R. (2015), A Guide to Graph Colouring, vol. 7, Springer.
  * Li et al. (2012) Li, B., Chun, H., and Zhao, H. (2012), “Sparse estimation of conditional graphical models with application to gene networks,” Journal of the American Statistical Association, 107, 152–167.
  * Li et al. (2017) Li, Y., Wen, Z., Yang, C., and Yuan, Y. (2017), “A semi-smooth Newton method for solving semidefinite programs in electronic structure calculations,” arXiv preprint arXiv:1708.08048.
  * Liu et al. (2010) Liu, H., Chen, X., Wasserman, L., and Lafferty, J. D. (2010), “Graph-valued regression,” in Advances in Neural Information Processing Systems, pp. 1423–1431.
  * Lounici et al. (2011) Lounici, K., Pontil, M., Van De Geer, S., and Tsybakov, A. B. (2011), “Oracle inequalities and optimal inference under group sparsity,” The Annals of Statistics, 39, 2164–2204.
  * Meinshausen and Bühlmann (2006) Meinshausen, N. and Bühlmann, P. (2006), “High-dimensional graphs and variable selection with the lasso,” The Annals of Statistics, 34, 1436–1462.
  * Merlevède et al. (2011) Merlevède, F., Peligrad, M., and Rio, E. (2011), “A Bernstein type inequality and moderate deviations for weakly dependent sequences,” Probability Theory and Related Fields, 151, 435–474.
  * Negahban et al. (2012) Negahban, S. N., Ravikumar, P., Wainwright, M. J., and Yu, B. (2012), “A unified framework for high-dimensional analysis of $M$-estimators with decomposable regularizers,” Statistical Science, 27, 538–557.
  * Network et al. (2008) Network, C. G. A. R. et al. (2008), “Comprehensive genomic characterization defines human glioblastoma genes and core pathways,” Nature, 455, 1061.
  * Ni et al. (2019) Ni, Y., Stingo, F. C., and Baladandayuthapani, V. (2019), “Bayesian graphical regression,” Journal of the American Statistical Association, 114, 184–197.
  * Peng et al. (2009) Peng, J., Wang, P., Zhou, N., and Zhu, J. (2009), “Partial correlation estimation by joint sparse regression models,” Journal of the American Statistical Association, 104, 735–746.
  * Samuels and Velculescu (2004) Samuels, Y. and Velculescu, V. E. (2004), “Oncogenic mutations of PIK3CA in human cancers,” Cell Cycle, 3, 1221–1224.
  * Simon et al. (2013) Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2013), “A sparse-group lasso,” Journal of Computational and Graphical Statistics, 22, 231–245.
  * van der Wijst et al. (2018) van der Wijst, M. G., Brugge, H., de Vries, D. H., Deelen, P., Swertz, M. A., and Franke, L. (2018), “Single-cell RNA sequencing identifies celltype-specific cis-eQTLs and co-expression QTLs,” Nature Genetics, 50, 493–497.
  * Wang et al. (2019) Wang, J., Zhang, Z., and Ye, J. (2019), “Two-Layer Feature Reduction for Sparse-Group Lasso via Decomposition of Convex Sets.” J. Mach. Learn. Res., 20, 163–1.
  * Wang et al. (2013) Wang, L., Zheng, W., Zhao, H., and Deng, M. (2013), “Statistical analysis reveals co-expression patterns of many pairs of genes in yeast are jointly regulated by interacting loci,” PLoS Genet, 9, e1003414.
  * Wang et al. (2012) Wang, Y., Joseph, S. J., Liu, X., Kelley, M., and Rekaya, R. (2012), “SNPxGE2: a database for human SNP–coexpression associations,” Bioinformatics, 28, 403–410.
  * Yin and Li (2011) Yin, J. and Li, H. (2011), “A sparse conditional gaussian graphical model for analysis of genetical genomics data,” The Annals of Applied Statistics, 5, 2630.
  * Zhang and Li (2022) Zhang, J. and Li, Y. (2022), “High-Dimensional Gaussian Graphical Regression Models with Covariates,” Journal of the American Statistical Association, published online March 14, 2022. https://doi.org/10.1080/01621459.2022.2034632.
  * Zhang et al. (2019) Zhang, J., Sun, W. W., and Li, L. (2019), “Mixed-Effect Time-Varying Network Model and Application in Brain Connectivity Analysis,” Journal of the American Statistical Association, 1–15.
  * Zhang et al. (2020) Zhang, Y., Zhang, N., Sun, D., and Toh, K.-C. (2020), “An efficient Hessian based algorithm for solving large-scale sparse group Lasso problems,” Mathematical Programming, 179, 223–263.