# An ADMM approach for multi-response regression with overlapping groups and
interaction effects

Theophilus Quachie Asensoa*, Manuela Zucknicka CONTACT Theophilus Quachie
Asenso. Email<EMAIL_ADDRESS>a Oslo Center for Epidemiology and
Biostatistics, Institute of Basic Medical Sciences, University of Oslo

###### Abstract

In this paper, we consider the regularized multi-response regression problem
where there exists some structural relation within the responses and also
between the covariates and a set of modifying variables. To handle this
problem, we propose MADMMplasso, a novel regularized regression method. This
method is able to find covariates and their corresponding interactions, with
some joint association with multiple related responses. We allow the
interaction term between covariate and modifying variable to be included in a
(weak) asymmetrical hierarchical manner by first considering whether the
corresponding covariate main term is in the model. For parameter estimation,
we develop an ADMM algorithm that allows us to implement the overlapping
groups in a simple way. The results from the simulations and analysis of a
pharmacogenomic screen data set show that the proposed method has an advantage
in handling correlated responses and interaction effects, both with respect to
prediction and variable selection performance.

Key Words: pliable lasso, Lagrange multipliers, drug sensitivity,
interactions, tree lasso, overlapping groups.

##  1 Introduction

The pliable lasso (plasso) (Tibshirani and Friedman, 2020) allows one to solve
problems involving main effects and their corresponding interaction effects.
Let $\mathbf{y}\in\mathbb{R}^{N}$ represent the response from a linear
regression model with $N$ observations, $X\in\mathbb{R}^{N\times p}$ be the
matrix of covariates with entries $x_{ij}$ and $Z\in\mathbb{R}^{N\times K}$ be
the matrix for the modifying variables with entries $z_{ik}$ for
$i=1,\ldots,N$, $j=1,\ldots,p$, $k=1,\ldots,K$. Let $X_{j}$ be the $j^{th}$
column of $X$, $Z_{k}$ be the $k^{th}$ column of $Z$ and $\mathbf{1}$ be an
$N$-vector of ones. The pliable lasso model is given as

$\begin{split}\hat{\mathbf{y}}&=\beta_{0}\mathbf{1}+Z\theta_{0}+\sum_{j=1}^{p}X_{j}(\beta_{j}\mathbf{1}+Z\bm{\theta}_{j})\\\
&=\beta_{0}\mathbf{1}+Z\theta_{0}+X{\beta}+\sum_{j=1}^{p}(X_{j}\circ
Z)\bm{\theta}_{j},\end{split}$ (1)

where ${\beta}\in\mathbb{R}^{p}$, $\theta$ is $p\times K$ matrix of parameters
with $j^{th}$ row $\bm{\theta}_{j}$ and individual entries $\theta_{jk}$,
$(X_{j}\circ Z)\in\mathbb{R}^{N\times K}$ is a matrix formed by multiplying
each column of $Z$ component-wise with the column vector $X_{j}$ and
$\theta_{0}\in\mathbb{R}^{K\times 1}$. The objective function of model (1) is
given as

$M(\beta_{0},\theta_{0},\beta,\theta)=\frac{1}{2N}\sum_{i}(y_{i}-\hat{y}_{i})^{2}\\\
+(1-\alpha)\lambda\sum_{j=1}^{p}(\lVert(\beta_{j},\bm{\theta}_{j})\rVert_{2}+\lVert\bm{\theta}_{j}\rVert_{2})+\alpha\lambda\sum_{j,k}\lvert\theta_{j,k}\rvert,$
(2)

where $\lambda$ and $\alpha$ are tuning parameters, both $\geq 0$. The model
assumes an asymmetric weak hierarchy constraint which is introduced as an
overlapping group in the first term of the penalty. The constraints ensure
that the interaction term can be nonzero only if the corresponding main term
is nonzero. Even though the idea is still young, it has been applied in
different areas, for example to multinomial logistic regression (Asenso et
al., 2022b), Cox’s proportional hazards model (Du and Tibshirani, 2018) and
support vector machines (Asenso et al., 2022a). However, in all the above
studies, the block-wise coordinate descent procedure was used in solving the
problem which includes overlapping groups. The algorithm involves multiple
“if” statements and a generalized gradient at the final stage. This implies
that extending the model to a multi-response case would require rigorous
computations like the case of Li et al. 2015, which might be difficult to
handle.

In this paper, we introduce the alternating direction method of multipliers
(ADMM) to handle this problem and extend the results from the single response
model to a multi-response problem. We provide a publicly available software
package MADMMplasso (Asenso and Zucknick, 2022) implemented in R. We present a
brief review on the ADMM algorithm in what follows.

###  1.1 A gentle introduction to ADMM

Details of the results in this section can be obtained from the works of Boyd
et al. (2011), Ye et al. (2011b), Deng et al. (2013), Han and Yuan (2012),
Humayoo and Cheng (2019a) and Haris et al. (2016). Here we give a brief
introduction to the ADMM algorithm. Given a separable objective function

$\underset{\beta}{\min}f(\beta)+h(\beta),$ (3)

where $f$ and $h$ are both convex, closed and proper, by introducing an
auxiliary variable $\omega$, ADMM allows (3) to be re-written as

$\underset{\beta,\omega}{\min}f(\beta)+h(\omega)\quad\text{s.t}\quad\beta=\omega.$
(4)

The problem in (4) can have a corresponding augmented Lagrangian in the form

$L(\beta,\omega,\gamma)=f(\beta)+h(\omega)+\gamma(\beta-\omega)+(\rho/2)\lVert\beta-\omega\rVert_{2}^{2},$
(5)

where $\rho$ is the augmented Lagrangian parameter or the multiplier. It has
been established that the ADMM algorithm converges to the global optimum,
under few simple conditions. The algorithm begins with initial parameters
$\epsilon^{pri}>0,\epsilon^{dual}>0,\beta^{0},\omega^{0},\gamma^{0}$ and in
each step $t+1$, with $t=1,2,3,\ldots$ indexing the iterations, the ADMM
algorithm updates $\beta$ and $\omega$ in an alternating or sequential manner
in the following way until the convergence condition is met.

$\begin{split}\beta^{t+1}&=\underset{\beta}{\arg\min}\quad
L(\beta,\omega^{t},\gamma^{t})\\\
\omega^{t+1}&=\underset{\omega}{\arg\min}\quad
L(\beta^{t+1},\omega,\gamma^{t})\\\
\gamma^{t+1}&=\gamma^{t}+\rho(\beta^{t+1}-\omega^{t+1}).\end{split}$ (6)

The algorithm meets the convergence condition if $s^{t+1}<\epsilon^{dual}$ and
$r^{t+1}<\epsilon^{pri}$, where
$s^{t+1}=\rho\lVert\omega^{t+1}-\omega^{t}\rVert_{2}$ and
$r^{t+1}=\lVert\beta^{t+1}-\omega^{t+1}\rVert_{2}$, representing the dual and
primal residuals respectively. The augmented Lagrangian parameter $\rho$ can
either be fixed or allowed to change by using the rule:

$\rho^{t+1}=\begin{cases}2\rho^{t}&\text{if}\quad\rho^{t}>10s^{t}\\\
\rho^{t}/2&\text{if}\quad 10\rho^{t}<s^{t}\\\
\rho^{t}&\text{otherwise}.\end{cases}$ (7)

The method has received attention (for example in the works of Anbari et al.
(2014), Ye et al. (2011a), Goldstein et al. (2014), Han and Yuan (2012),
Humayoo and Cheng (2019b) and Qin and Goldfarb (2012)) and has been applied to
different optimisation problems including lasso and group lasso and a
possibility to include overlapping groups which was introduced in the work of
Deng et al. (2013). In this paper, we study the possible application of the
ADMM algorithm to handle the overlapping groups in the pliable lasso objective
function.

The rest of the paper is organized as follows. In section (2), we present two
ADMM algorithms. The first algorithm is an ADMM version of the original
pliable lasso linear model. We show that ADMM can handle the overlapping
groups in the objective function and will result in the same solution as using
the block-wise coordinate descent method in the original pliable lasso paper.
The second algorithm is to solve a multi-response problem that has the pliable
lasso penalty, and in addition a penalty to consider the correlation between
the responses. We assume that the correlation between the response follows a
hierarchical tree structure and the responses form groups at various levels
within the hierarchy. This implies that each response could belong to several
groups, hence the need to allow for overlapping groups. In section (3) and
(4), we show both the simulations and a real data application, respectively of
the proposed method. Section (5) is for the conclusion.

##  2 Method

###  2.1 Model for single response

Let $B\in\mathbb{R}^{p\times(K+1)}$ be a matrix of coefficients, with rows and
columns indexed from $1$ to $p$ and $1$ to $(K+1)$ respectively, representing
the combination of $X$ and $Z$ with the $j^{th}$ row of B defined as
$B_{j}=[\beta_{j},\bm{\theta}_{j}]\in\mathbb{R}^{K+1}$. We can reconstruct the
problem in (1) using array notations. Let $W$ be an $p\times(1+K)\times N$
array constructed as follows for $i=1,\ldots,N$, $j=1,\ldots,p$,
$k=1,\ldots,K+1$

$W_{jki}=\begin{cases}x_{ij}z_{ik}&\text{for}\quad k\neq 1\\\
x_{ij}&\text{for}\quad k=1,\end{cases}$ (8)

then equation (1) is equivalent to the model solution

$\hat{\bm{y}}=\beta_{0}\mathbf{1}+Z\theta_{0}+W*\mathbf{\textit{B}},$ (9)

where $B$ is the matrix of the coefficients and $W*B$ denotes the $N$ vector
whose $i^{th}$ element takes the form

$(W*B)_{i}=\sum_{j=1}^{p}\sum_{k=1}^{K+1}W_{jki}B_{jk}.$ (10)

In equation (10) above,
$W_{ji}=[x_{ij},x_{ij}z_{i1},x_{ij}z_{i2},\ldots,x_{ij}z_{iK}]\in\mathbb{R}^{1\times(K+1)}$.
By using (9), the objective function of the pliable lasso model in (2) can be
rewritten as

$\underset{\mathbf{\textit{B}}\in\mathbb{R}^{p\times(1+K)}}{\min}\frac{1}{2N}\lVert\bm{y-\hat{y}}\rVert_{2}^{2}+(1-\alpha)\lambda\sum_{j=1}^{p}(\lVert
B_{j}\rVert_{2}+\lVert B_{j(-1)}\rVert_{2})+\alpha\lambda\sum_{j=1}^{p}\lVert
B_{j(-1)}\rVert_{1}.$ (11)

Here, the largest group is of the size $K+1$ (to represent
$[\beta_{j},\bm{\theta}_{j}]$) and there are two groups (with elements
$[\beta_{j},\bm{\theta}_{j}],[\bm{\theta}_{j}]$). Also, $B_{j(-1)}$ represents
the $j^{\text{th}}$ row of $B$ without the first element i.e, without the
coefficient of the main effect $\beta_{j}$. In what follows in this section,
we abuse notation and write $K=K+1$ unless otherwise stated. Let $\\{\
B_{j}^{g_{s}}\in\mathbb{R}^{k_{s}};s=1,2\\}$ be the groupings of $B_{j}$,
where $g_{s}$ with $\underset{s}{\cup}g_{s}\subseteq\\{1,2,\ldots,K\\}$ is an
index set corresponding to the $s^{th}$ group, and $B_{j}^{g_{s}}$, which is a
row vector, denotes the subvector of $B_{j}$ indexed by $g_{s}$ and $k_{s}$
represents the size of $g_{s}$. This implies that
$B_{j}^{g_{1}}=B_{j}=[\beta_{j},\bm{\theta}_{j}]\quad\text{and}\quad
B_{j}^{g_{2}}=B_{j(-1)}=[\bm{\theta}_{j}]$. To solve the problem (11), we
introduce 2 auxiliary variables as follows; we introduce a matrix
$V\in\mathbb{R}^{p\times\tilde{K}}$, with each $j^{th}$ row vector
$V_{j}\in\mathbb{R}^{\tilde{K}}$, defined for the two groups
$(B_{j}^{g_{1}},B_{j}^{g_{2}})$, and let
$\tilde{B}_{j}=[B_{j}^{g_{1}},B_{j}^{g_{2}}]\in\mathbb{R}^{\tilde{K}}$, as a
row vector. This is constructed for the overlapping group part of the pliable
lasso penalty and let us note that $\tilde{K}=\sum_{s=1}^{2}k_{s}\geq K$ and
$V_{j}^{1}=B_{j}^{g_{1}}$ and $V_{j}^{2}=B_{j}^{g_{2}}$. The second auxiliary
variable, $Q\in\mathbb{R}^{p\times K}$ is introduced to handle the
$l_{1}$-penalty term. By introducing these auxiliary variables in problem
(11), we obtain the following optimization problem:

$\underset{\underset{V,Q}{B\in\mathbb{R}^{p\times
K}},}{\min}\quad\frac{1}{2N}\lVert\bm{y-\hat{y}}\rVert_{2}^{2}+(1-\alpha)\lambda\sum_{j=1}^{p}\sum_{s=1}^{2}\lVert
V_{j}^{s}\rVert_{2}+\alpha\lambda\sum_{j=1}^{p}\lVert Q_{j}\rVert_{1}\\\
\text{s.t}\quad B_{j}^{g_{s}}=V_{j}^{s},\quad B_{j}=Q_{j}.$ (12)

The augmented Lagrangian can be written as

$L(B,V,Q,O,P)=\lVert\bm y-\hat{\bm
y}\rVert_{2}^{2}+(1-\alpha)\lambda\sum_{j=1}^{p}\sum_{s=1}^{2}\lVert
V_{j}^{s}\rVert_{2}+\alpha\lambda\sum_{j=1}^{p}\lVert Q_{j}\rVert_{1}\\\
+\sum_{j}O_{j}(\tilde{B}_{j}-V_{j})^{T}+\langle
P,B-Q\rangle+\frac{\rho}{2}\sum_{j}\sum_{s}\lVert
B_{j}^{g_{s}}-V_{j}^{s}\rVert_{2}^{2}+\frac{\rho}{2}\lVert B-Q\rVert_{2}^{2},$
(13)

where $O\in\mathbb{R}^{p\times\tilde{K}}$ and $P\in\mathbb{R}^{p\times K}$ are
the Lagrange multipliers associated with $V$ and $Q$ respectively and $\langle
A,B\rangle\equiv Tr(A^{T}B)$ denotes the matrix inner product.

####  2.1.1 Updates for the ADMM algorithm

In each iteration, the algorithm will update $B,V,Q,P$ and $O$ alternatively
until the stopping conditions are met. We show the various updates in what
follows:

  1. 1.

$B_{j}$ update: We begin by showing how
$\tilde{B}_{j}=[B_{j}^{g_{1}},B_{j}^{g_{2}}]\in\mathbb{R}^{\tilde{K}}$ is
constructed. Let $\tilde{B}_{j}=(GB_{j}^{T})^{T}$ with each row of
$G\in\mathbb{R}^{\tilde{K}\times K}$, defined as $G_{(s,l)}^{k}=1$ for $k\in
g_{s}$ and $0$ elsewhere, with rows indexed by $(s,l)$ such that
$s\in\\{1,2\\}$ and $l\in g_{s}$. For example given $B\in\mathbb{R}^{5\times
4}$ and knowing that in each $j$ are two groups (with elements
$[\beta_{j},\bm\theta_{j}],[\bm\theta_{j}]$), $G$ will be represented as a
$(2\times 4)\times 4$ matrix with elements,

$\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\ 0&0&0&1\\\ 0&0&0&0\\\
0&1&0&0\\\ 0&0&1&0\\\ 0&0&0&1\\\ \end{bmatrix}$ (14)

From (14) above, $\tilde{B}_{j}$ will have elements
$[B_{j1},B_{j2},B_{j3},B_{j4},0,B_{j2},B_{j3},B_{j4}]$, which is the same as
the original elements
$[\beta_{j},\theta_{j1},\theta_{j2},\theta_{j3},0,\theta_{j1},\theta_{j2},\theta_{j3}]$,
which is an expansion of $[\beta_{j},\bm\theta_{j}],[\bm\theta_{j}]$. The
update for $B_{j}$ is given as

$\underset{B_{j}}{\arg\min}\frac{1}{2N}\sum_{i=1}^{N}(y_{i}-\beta_{0}-Z_{i}\theta_{0}-W_{ji}\mathbf{\mathit{{B}}}_{j}^{T}-\sum_{t\neq
j}W_{ti}\mathbf{\mathit{{B}}}_{t}^{T})^{2}\\\
+O_{j}(\tilde{B}_{j})^{T}+P_{j}(B_{j})^{T}+\frac{\rho}{2}\lVert\tilde{B}_{j}-V_{j}\rVert_{2}^{2}+\frac{\rho}{2}\lVert
B_{j}-Q_{j}\rVert_{2}^{2},\quad t\neq j,$ (15)

where $Z_{i}$ is the $i^{th}$ row of $Z$ with individual elements $Z_{ik}$.
Let $\tilde{\bm{y}}=\bm{y}-\beta_{0}\mathbf{1}_{N}-Z\theta_{0}$ and define
$\bm{r}_{(-j)}=\tilde{\bm{y}}-\sum_{t\neq j}W_{t}\mathbf{\mathit{{B}}}_{t}$ as
the partial residual with the $j^{th}$ group removed, where we have
represented $W_{j}$ as
$W_{j}=[W_{j1}^{T}:W_{j2}^{T}:\ldots:W_{jN}^{T}]^{T}\in\mathbb{R}^{N\times
K}$. Further computations will show that,

$B_{j}=\left(\frac{W_{j}^{T}W_{j}}{N}+\rho(G^{T}G+I)\right)^{-1}\left(\frac{1}{N}W_{j}^{T}\bm{r}_{(-j)}+\rho(G^{T}V_{j}^{T}+Q_{j}^{T}-G^{T}O_{j}^{T}-P_{j}^{T})\right)$
(16)

  2. 2.

$V_{j}$ update: From (13), we can update $V$ by solving the problem

$\underset{V_{j}}{\arg\min}\quad(1-\alpha)\lambda\sum_{s=1}^{2}\lVert
V_{j}^{s}\rVert_{2}+O_{j}(-V_{j})^{T}+\frac{\rho}{2}\sum_{s}\lVert
B_{j}^{g_{s}}-V_{j}^{s}\rVert_{2}^{2}.$ (17)

This can be simplified as

$\underset{V_{j}}{\arg\min}\sum_{s=1}^{2}\left[\frac{\rho}{2}\Bigg{\|}B_{j}^{g_{s}}-V_{j}^{s}+\frac{O_{j}^{s}}{\rho}\Bigg{\|}_{2}^{2}+(1-\alpha)\lambda\lVert
V_{j}^{s}\rVert_{2}\right].$ (18)

Further computations similar to the group lasso problem (Li et al., 2015; Deng
et al., 2013) will show that (18) has a closed form solution by the soft
threshold formula

$\begin{split}V_{j}^{s}&=\max\Bigg{\\{}\lVert
r_{s}\rVert_{2}-\frac{(1-\alpha)\lambda}{\rho},0\Bigg{\\}}\frac{r_{s}}{\lVert
r_{s}\rVert_{2}},\\\ r_{s}&=B_{j}^{g_{s}}+O_{j}^{s}\end{split}$ (19)

  3. 3.

$Q_{j}$ update: To update $Q_{j}$, we should note that we added the main
effect $\beta_{j}$ to the $B_{j}$ vector in the previous updates. However,
since we wish to avoid the effect of the $l_{1}$-norm on it, we remove the
penalty effect here. The $Q$ component of the optimization problem (13) is
updated as

$\underset{Q_{j}}{\arg\min}\quad\alpha\lambda\lVert
Q_{j}\rVert_{1}+P_{j}(-Q_{j})^{T}+\frac{\rho}{2}\lVert
B_{j}-Q_{j}\rVert_{2}^{2},$ (20)

which can further be expressed as

$\underset{Q_{j}}{\arg\min}\quad\frac{\rho}{2}\Bigg{\|}B_{j}-Q_{j}+\frac{P_{j}}{\rho}\Bigg{\|}_{2}^{2}+\alpha\lambda\lVert
Q_{j}\rVert_{1}.$ (21)

Further computations similar to the lasso problem (Tibshirani, 1996; Hastie et
al., 2015) will show that

$\begin{split}Q_{jk}&=sign\left(B_{jk}+P_{jk}\right)\left(\Bigg{|}B_{jk}+P_{jk}\Bigg{|}-\frac{\alpha\lambda}{\rho}\right)_{+},\quad\text{for}\quad
k\neq 1\\\ Q_{j1}&=B_{j1}+P_{j1},\quad\text{for}\quad k=1,\end{split}$ (22)

where $(a)_{+}=\max(a,0)$.

  4. 4.

$P$ and $O$ updates: The two Lagrange multipliers ($O,P$), associated with the
penalties in (13) are updated as follows;

$\begin{split}P&\leftarrow P+B-Q\\\ O&\leftarrow O+(GB^{T})^{T}-V\end{split}$
(23)

###  2.2 Extension to multi-response model

With the motivation from the model for a single-response variable, we proceed
to propose an ADMM method for a multi-response model, where the response
variables are assumed to have some correlation structures. The section begins
with the model formulation in what follows. Let $Y\in\mathbb{R}^{N\times D}$
represent the response from a regression model with $N$ number of
observations, with each observation having $D$ responses,
$X\in\mathbb{R}^{N\times p}$ be the matrix of covariates with entries $x_{ij}$
and $Z\in\mathbb{R}^{N\times K}$ be the matrix for the modifying variables
with entries $z_{ik}$ for $i=1,\ldots,N$, $j=1,\ldots,p$, $k=1,\ldots,K$. Let
$X_{j}$ be the $j^{th}$ column of $X$, $Z_{k}$ be the $k^{th}$ column of $Z$
and $\mathbf{1}$ be an $N$-vector of ones. The pliable lasso problem in (1)
can be extended to a multi-response model with interaction effect, where each
response $Y_{d}$ can be modeled as

$\hat{Y}_{d}=\beta_{0d}\mathbf{1}+Z\theta_{0d}+\underset{j}{\sum}X_{j}\beta_{jd}\mathbf{1}+\underset{j}{\sum}(X_{j}\circ
Z)\bm\theta_{jd},$ (24)

with $d=1,\ldots,D$ and $(X_{j}\circ Z)\in\mathbb{R}^{N\times K}$ is a matrix
formed by multiplying each column of $Z$ component-wise with the column vector
$X_{j}$. The parameter
$\beta_{0}=(\beta_{01},\beta_{02},\ldots,\beta_{0D})^{T}$ is a $D$ vector of
intercepts for the main effect, $\theta_{0}\in\mathbb{R}^{K\times D}$ is the
matrix of coefficients for the modifying variable $Z$ with each column
represented by $\theta_{0d}\in\mathbb{R}^{K}$, $\beta$ is a $p\times D$ matrix
of coefficients for the main effect with each row represented by $\beta_{j}$,
having individual entries $\beta_{jd}$, $\theta$ is a $p\times K\times D$
array of size $D$ with each $d=1,\ldots,D$ element having $p\times K$ matrix
where each $j^{th}$ row of the $d^{th}$ matrix is represented by
$\bm\theta_{jd}$. Let $B\in\mathbb{R}^{p\times(K+1)\times D}$ be an array of
size $D$ with each $d=1,\ldots,D$ representing a matrix of coefficients for
the main and interaction effects for response $d$, with rows and columns
indexed from $1$ to $p$ and $1$ to $(K+1)$, respectively. This is to cater for
the combination of $X$ and $Z$ with the $j^{th}$ row of $B_{d}$ defined as
$\mathbf{\textit{B}}_{jd}=[\beta_{jd},\bm\theta_{jd}]\in\mathbb{R}^{K+1}$. We
can reconstruct the equation (24) using array notations. Let $W$ be an
$p\times(1+K)\times N$ array with same construction as (8), then the model for
all responses in equation (24) is equivalent to the model

$\hat{Y}=\mathbf{1}\beta_{0}^{T}+Z\theta_{0}+[W*B_{1}:W*B_{2}:\ldots:W*B_{D}].$
(25)

We will abuse notation to define $W*B=[W*B_{1}:W*B_{2}:\ldots:W*B_{D}]$ to
denote the $N\times D$ matrix whose $id$ element takes the form

$(W*B)_{id}=\sum_{j=1}^{p}\sum_{k=1}^{K+1}W_{jki}B_{jkd},\quad i=1,\ldots
N,\quad d=1,\ldots,D.$ (26)

The objective function for a general multi-response pliable lasso model can be
written as

$\underset{B\in\mathbb{R}^{p\times(1+K)\times D}}{\min}\quad\frac{1}{2N}\lVert
Y-\hat{Y}\rVert_{F}^{2}\\\
+\sum_{d=1}^{D}\left[(1-\alpha)\lambda\sum_{j=1}^{p}(\lVert
B_{jd}\rVert_{2}+\lVert
B_{j(-1)d}\rVert_{2})+\alpha\lambda\sum_{j=1}^{p}\lVert
B_{j(-1)d}\rVert_{1}\right].$ (27)

Here also, the largest group is of the size $K+1$ (to represent
$[\beta_{jd},\bm\theta_{jd}]$) and there are two groups (with elements
$[\beta_{jd},\bm\theta_{jd}],[\bm\theta_{jd}]$). Also, $B_{j(-1)d}$ represents
the $j^{\text{th}}$ row of $\mathbf{\textit{B}}_{d}$ without the first element
i.e, without the coefficient of the main effect $\beta_{jd}$.

####  2.2.1 Multi-response pliable lasso with tree-guided structure

The problem (27) can be solved as a single-response problem with interaction
effects where each of the responses could be treated as independent of the
others. However, in our case, we assume correlation among the response
variables. To handle this correlation, we allow the response variables to form
overlapping groups and having a structure, specifically a tree-like structure
in this paper. More details about the tree-like structure can be found in the
paper “Tree-guided group lasso for multi-response regression with structured
sparsity, with an application to EQTL mapping” (Kim and Xing, 2012) and
similar structures can also be found in Li et al. (2015) and in Zhao and
Zucknick (2020).

We represent the set of internal and leaf nodes of the tree as
$M_{\text{int}}$, $M_{\text{leaf}}$ of size $|M_{\text{int}}|$ and
$|M_{\text{leaf}}|$ respectively, the group of responses forming an internal
node $m\in M_{\text{int}}$ as $\mathcal{G}_{m}$, where
$\mathcal{G}_{m}\subseteq\\{1,\ldots,D\\}$ and let $B_{j}^{\mathcal{G}_{m}}$
denotes the $j^{th}$ sub-vector of $B$, indexed by $\mathcal{G}_{m}$ with a
group weight $w_{m}$. Each sub-vector $B_{j}^{\mathcal{G}_{m}}$ has elements
$\\{B_{jd};d\in\mathcal{G}_{m}\\}$. Adding the overlapping groups in the
responses to the problem (27) gives us

$\underset{B\in\mathbb{R}^{p\times(1+K)\times D}}{\min}\frac{1}{2N}\lVert
Y-\hat{Y}\rVert_{F}^{2}+\lambda_{1}\sum_{j=1}^{p}\sum_{m\in
M_{\text{int}}}w_{m}\lVert
B_{j}^{\mathcal{G}_{m}}\rVert_{2}+\lambda_{2}\sum_{j=1}^{p}\sum_{m\in
M_{\text{leaf}}}w_{m}\lVert B_{j}^{\mathcal{G}_{m}}\rVert_{2}\\\
+\sum_{d}^{D}\left[(1-\alpha)\lambda_{3}\sum_{j=1}^{p}(\lVert
B_{jd}\rVert_{2}+\lVert
B_{j(-1)d}\rVert_{2})+\alpha\lambda_{3}\sum_{j=1}^{p}\lVert
B_{j(-1)d}\rVert_{1}\right].$ (28)

The first penalty terms ensure that correlated responses are penalized within
groups while ensuring sparsity among the coordinate $j$ and the second penalty
term is the pliable lasso part to ensure that main and interaction effects are
considered in a hierarchical format.

To solve the problem (28) by the ADMM approach, we introduce four auxiliary
variables as follows. Let $E\in\mathbb{R}^{(\sum_{m\in
M}\lvert\mathcal{G}_{m}\rvert)\times(p+p\times K)}$ be the matrix with each
$E_{j}\in\mathbb{R}^{\sum_{m\in M}\lvert\mathcal{G}_{m}\rvert}\times(K+1)$ to
represent the groups formed by internal nodes $M$ on the $j^{th}$ coordinate ,
and
$E^{\mathcal{G}_{m}}\in\mathbb{R}^{\lvert\mathcal{G}_{m}\rvert\times(p+p\times
K)}$ being the block representation of group $\mathcal{G}_{m}$ and equal to
$\mathit{B}^{\mathcal{G}_{m}}$. This is defined for the overlapping groups in
the responses for
$(\mathit{B}_{j}^{\mathcal{G}_{1}},\mathit{B}_{j}^{\mathcal{G}_{2}},\ldots,\mathit{B}_{j}^{\mathcal{G}_{M}})$
with
$\mathit{\tilde{\tilde{B}}}_{j}\equiv[\mathit{B}_{j}^{\mathcal{G}_{1}T},\mathit{B}_{j}^{\mathcal{G}_{2}T},\ldots,\mathit{B}_{j}^{\mathcal{G}_{M}T}]^{T}\in\mathbb{R}^{\sum_{m\in
M}\lvert\mathcal{G}_{m}\rvert\times(K+1)}$. Note that we have defined $E$ to
contain both the main effect and the interaction effect of each $j$, hence
each row of $E$ can be seen as having $p(K+1)$ columns. We also introduce an
array $\tilde{E}\in\mathbb{R}^{p\times(1+K)\times D}$ to handle the
$l_{2}$-penalty of the leaf node ($\lambda_{2}\sum_{j=1}^{p}\sum_{m\in
M_{\text{leaf}}}w_{m}\lVert B_{j}^{\mathcal{G}_{m}}\rVert_{2}$) in (28).

In what follows in this section, we abuse notation and write $K=K+1$ unless
otherwise stated. Let $\\{\ B_{jd}^{g_{s}}\in\mathbb{R}^{k_{s}};s=1,2\\}$ be
the grouping of elements of $B_{jd}$, where $g_{s}$ with
$\underset{s}{\cup}g_{s}\subseteq\\{1,\ldots,K\\}$ is an index set
corresponding to the $s^{th}$ group, and $B_{jd}^{g_{s}}$ which is a row
vector, denotes the subvector of $B_{jd}$, indexed by $g_{s}$ and $k_{s}$
represents the size of $g_{s}$. We introduce the second auxiliary variable as
an array $V$ of size $D$ with each of the index $d=1,\ldots,D$ having a matrix
$V_{d}\in\mathbb{R}^{p\times\tilde{K}}$, with each $j^{th}$ row vector
$V_{jd}\in\mathbb{R}^{\tilde{K}}$, defined for the two groups
$(B_{jd}^{g_{1}},B_{jd}^{g_{2}})$. This is for the overlapping part of the
pliable lasso penalty and let us note that $\tilde{K}=\sum_{s=1}^{2}k_{s}\geq
K$ and $V_{jd}^{1}=B_{jd}^{g_{1}}$ and $V_{jd}^{2}=B_{jd}^{g_{2}}$. Also, let
$\tilde{B}_{jd}\equiv[B_{jd}^{g_{1}},B_{jd}^{g_{2}}]\in\mathbb{R}^{\tilde{K}}$
and $\tilde{B}_{d}\in\mathbb{R}^{p\times\tilde{K}}$ be formed with elements
from $\tilde{B}_{jd}$ such that each row is $\tilde{B}_{jd}$.

Finally, we introduce an array $Q\in\mathbb{R}^{D\times(p\times K})$ to handle
the $l_{1}$-penalty of the pliable lasso. Implementing these auxiliary
variables on problem (28) leads to the following;

$\underset{B,E,\tilde{E},V,Q}{\min}\quad\frac{1}{2N}\lVert
Y-\hat{Y}\rVert_{F}^{2}+\lambda_{1}\sum_{j=1}^{p}\sum_{m\in
M_{\text{int}}}w_{m}\lVert
E_{j}^{\mathcal{G}_{m}}\rVert_{2}+\lambda_{2}\sum_{j=1}^{p}\sum_{m\in
M_{\text{leaf}}}w_{m}\lVert\tilde{E}_{j}^{\mathcal{G}_{m}}\rVert_{2}\\\
+\sum_{d}(1-\alpha)\lambda_{3}\sum_{j=1}^{p}\sum_{s=1}^{2}\lVert
V_{jd}^{s}\rVert_{2}+\alpha\lambda_{3}\sum_{j=1}^{p}\lVert Q_{jd}\rVert_{1}\\\
\text{s.t}\quad B_{j}^{\mathcal{G}_{m}}=E_{j}^{\mathcal{G}_{m}},\quad
B_{d}=\tilde{E}_{d},\quad B_{jd}^{g_{s}}=V_{jd}^{s},\quad B_{d}=Q_{d}.$ (29)

In the above, $\tilde{E}_{j}^{\mathcal{G}_{m}}$ is the $j^{th}$ component of
each $\tilde{E}_{d}$ and since it is a leaf node, each will have a single
element. Therefore the penalty $\lambda_{2}\sum_{j=1}^{p}\sum_{m\in
M_{\text{leaf}}}w_{m}\lVert\tilde{E}_{j}^{\mathcal{G}_{m}}\rVert_{2}$ can be
rewritten as
$\lambda_{2}\sum_{d}\sum_{j=1}^{p}w_{d}\lVert\tilde{E}_{jd}\rVert_{2}$ which
is just a further shrinkage on each $j$. The augmented Lagrangian form of (29)
can be written as

${L}(B,E,\tilde{E},V,Q,H,\tilde{H},O,P)=\frac{1}{2N}\lVert
Y-\hat{Y}\rVert_{F}^{2}+\lambda_{1}\sum_{j=1}^{p}\sum_{m\in
M_{\text{int}}}w_{m}\lVert E_{j}^{\mathcal{G}_{m}}\rVert_{2}\\\
+\lambda_{2}\sum_{d}\sum_{j=1}^{p}w_{d}\lVert\tilde{E}_{jd}\rVert_{2}+\sum_{d}(1-\alpha)\lambda_{3}\sum_{j=1}^{p}\sum_{s}\lVert
V_{jd}^{s}\rVert_{2}+\alpha\lambda_{3}\sum_{j=1}^{p}\lVert Q_{jd}\rVert_{1}\\\
+\sum_{j}H_{j}(\tilde{\tilde{B}}_{j}-E_{j})^{T}+\sum_{d}\langle\tilde{H}_{d},B_{d}-\tilde{E}_{d}\rangle+\sum_{d}\sum_{j}O_{jd}(\tilde{B}_{jd}-V_{jd})^{T}+\sum_{d}\langle
P_{d},B_{d}-Q_{d}\rangle\\\
+\frac{\rho}{2}\sum_{j}\lVert\tilde{\tilde{B}}_{j}-E_{j}\rVert_{2}^{2}+\frac{\rho}{2}\sum_{d}\lVert
B_{d}-\tilde{E}_{d}\rVert_{2}^{2}+\frac{\rho}{2}\sum_{d}\sum_{j}\sum_{s}\lVert\tilde{B}_{jd}^{s}-V_{jd}^{s}\rVert_{2}^{2}+\frac{\rho}{2}\sum_{d}\lVert
B_{d}-Q_{d}\rVert_{2}^{2}.$ (30)

In the problem (30) above, $H,\tilde{H},O,P$ have the same dimensions as
$E,\tilde{E},V$ and $Q$ respectively, and are called the Lagrange multipliers
and $\rho$ is a penalty parameter.

####  2.2.2 Update for the ADMM algorithm

Just like for the single-response problem, the algorithm will update
$B,E,\tilde{E},V,Q,H,\tilde{H},P$ and $O$ alternatively in each iteration
until the stopping criteria are met. We show the various updates in what
follows:

  1. 1.

$B$ update: We begin by showing the construction of $\tilde{B}_{d}$ and
$\tilde{\tilde{B}}$. Let $\tilde{B}_{d}=(GB_{d}^{T})^{T}$ with each row of
$G\in\mathbb{R}^{\tilde{K}\times K}$, defined as $G_{(s,l)}^{k}=1$ for $k\in
g_{s}$ and $0$ elsewhere, with rows indexed by $(s,l)$ such that
$s\in\\{1,2\\}$ and $l\in g_{s}$. Also, let $\tilde{\tilde{B}}=IB^{T}$ with
$I\in\mathbb{R}^{(\sum\lvert\mathcal{G}_{m}\rvert)\times D}$, and
$I_{m,u}^{d}=w_{m}$ for $d\in\mathcal{G}_{m}$ and $0$ otherwise, indexed by
$(m,u)$ such that $m\in M$ and $u\in\mathcal{G}_{m}$. An example of these two
constructions is shown in (14). Note that we have written the matrix
$B\in\mathbb{R}^{(p+p\times K)\times D}$ in the construction of
$\tilde{\tilde{B}}$ above as a transformation of the array $B$ such that, each
column $B_{d}$ (of matrix $B$) is a vectorized version of the matrix $B_{d}$
of the array $B$. The $B$ component of (30) is updated as follows;

$\underset{B}{\arg\min}\quad\frac{1}{2N}\lVert
Y-\mathbf{1}\beta_{0}^{T}-Z\theta_{0}-W*B\rVert_{F}^{2}\\\
+H(\tilde{\tilde{B}})^{T}+\sum_{d}\tilde{H}_{d}(B_{d})^{T}+\sum_{d}O_{d}(\tilde{B}_{d})^{T}+\sum_{d}P_{d}(B_{d})^{T}\\\
+\frac{\rho}{2}\lVert\tilde{\tilde{B}}-E\rVert_{2}^{2}+\frac{\rho}{2}\sum_{d}\lVert
B_{d}-\tilde{E}_{d}\rVert_{2}^{2}+\frac{\rho}{2}\sum_{d}\lVert\tilde{B}_{d}-V_{d}\rVert_{2}^{2}+\frac{\rho}{2}\sum_{d}\lVert
B_{d}-Q_{d}\rVert_{2}^{2}.$ (31)

Let $\tilde{W}$ denote the $N\times(p+p\times K)$-dimensional matrix of $W$
and let $R=Y-\mathbf{1}\beta_{0}^{T}-Z\theta_{0}$, with each column $R_{d},$
being an $N$ column vector. Taking the derivative of (31) with respect to
$B_{d}$ and equating to zero gives

$\begin{split}B_{d}&=\frac{\phi}{\frac{\tilde{W}^{T}\tilde{W}}{N}+C_{d}},\\\
\text{where}\quad\phi=&\frac{1}{N}\tilde{W}^{T}R_{d}+(\rho(I^{T}E-I^{T}H))_{d}\\\
&+vec(\rho(EE_{d}-HH_{d})(\rho(G^{T}(V_{d})^{T}-G^{T}(O_{d})^{T}))^{T}+\rho(Q_{d}-P_{d})),\end{split}$
(32)

and we have allowed $\frac{a}{b}$ to represent $(b)^{-1}a$ for the sake of
simplicity and $vec(a)$ is a vectorized form of $a$. The resulting $B_{d}$
which is a vector, is then transformed into a matrix of dimension
$p\times(K)$. To construct $C$ in (32), we first defined $\tilde{I}$ as the
diagonal of $I^{T}I$ and $\tilde{G}$ to be the diagonal of $G^{T}G$ which has
elements $\\{1,2,2,2,\ldots,2\\}$. Then we defined $C$ to be an array of size
$D$ with each $d=1,\ldots,D$ being a $(p+p\times K)\times(p+p\times K)$
diagonal matrix with the first $p$ elements equal to
$\rho(1+\tilde{G}_{1})+\rho(\tilde{I}_{d}+1)$ and the $p+1$ to $p+p\times K$
elements equal to $3\rho+\rho(\tilde{I}_{d}+1)$ since the elements of
$\tilde{G}_{j}=2,j=2,\ldots,K$. This construction is to cater for the
groupings within the responses and within each coordinate $j=1,\ldots,p$.

  2. 2.

$V_{jd}$ update: from (30), we update $V$ by solving the problem

$\underset{V_{jd}}{\arg\min}\quad(1-\alpha)\lambda_{3}\sum_{s}\lVert
V_{jd}^{s}\rVert_{2}+O_{jd}(-V_{jd})^{T}+\frac{\rho}{2}\sum_{s}\lVert\tilde{B}_{jd}^{s}-V_{jd}^{s}\rVert_{2}^{2}.$
(33)

With some few computations, this can be simplified as

$\underset{V_{jd}}{\arg\min}\sum_{s=1}^{2}\left[\frac{\rho}{2}\Bigg{\|}\tilde{B}_{jd}^{s}-V_{jd}^{s}+\frac{O_{jd}^{s}}{\rho}\Bigg{\|}_{2}^{2}+(1-\alpha)\lambda_{3}\lVert
V_{jd}^{s}\rVert_{2}\right].$ (34)

Further computations similar to the group lasso problem will show that,

$\begin{split}V_{jd}^{s}&=\max\Bigg{\\{}\lVert
r_{s}\rVert_{2}-\frac{(1-\alpha)\lambda_{3}}{\rho},0\Bigg{\\}}\frac{r_{s}}{\lVert
r_{s}\rVert_{2}},\\\ r_{s}&=\tilde{B}_{jd}^{s}+O_{jd}^{s}\end{split}$ (35)

  3. 3.

$Q_{jd}$ update: The $Q$ component of the optimization problem (30) is updated
as

$\underset{Q_{jd}}{\arg\min}\quad\alpha\lambda_{3}\lVert
Q_{jd}\rVert_{1}+P_{jd}(-Q_{jd})^{T}+\frac{\rho}{2}\lVert
B_{jd}-Q_{jd}\rVert_{2}^{2},$ (36)

which can further be expressed as

$\underset{Q_{jd}}{\arg\min}\quad\frac{\rho}{2}\Bigg{\|}B_{jd}-Q_{jd}+\frac{P_{jd}}{\rho}\Bigg{\|}_{2}^{2}+\alpha\lambda_{3}\lVert
Q_{jd}\rVert_{1}.$ (37)

Further computations, similar to the lasso problem will show that,

$\begin{split}Q_{jkd}&=sign\left(B_{jkd}+P_{jkd}\right)\left(\Bigg{|}B_{jkd}+P_{jkd}\Bigg{|}-\frac{\alpha\lambda_{3}}{\rho}\right)_{+},\quad\text{for}\quad
k\neq 1\\\ Q_{j1d}&=B_{j1d}+P_{j1d},\quad\text{for}\quad k=1.\end{split}$ (38)

  4. 4.

$E_{j}^{\mathcal{G}_{m}}$ update: The $E$ component of the optimization
problem (30) is updated as

$\underset{E_{j}^{\mathcal{G}_{m}}}{\arg\min}\quad\lambda_{1}w_{m}\lVert
E_{j}^{\mathcal{G}_{m}}\rVert_{2}+H_{j}^{\mathcal{G}_{m}}(-E_{j}^{\mathcal{G}_{m}})^{T}+\frac{\rho}{2}\lVert\tilde{\tilde{B}}_{j}^{\mathcal{G}_{m}}-E_{j}^{\mathcal{G}_{m}}\rVert_{2}^{2},$
(39)

which can be simplified as

$\underset{E_{j}^{\mathcal{G}_{m}}}{\arg\min}\left[\lambda_{1}w_{m}\lVert
E_{j}^{\mathcal{G}_{m}}\rVert_{2}+\frac{\rho}{2}\Bigg{\|}\tilde{\tilde{B}}_{j}^{\mathcal{G}_{m}}-E_{j}^{\mathcal{G}_{m}}+\frac{H_{j}^{\mathcal{G}_{m}}}{\rho}\Bigg{\|}_{2}^{2}\right].$
(40)

Further computations, by following the group lasso example will show that,

$\begin{split}E_{j}^{\mathcal{G}_{m}}&=\max\Bigg{\\{}\lVert
r_{m}\rVert_{2}-\frac{\lambda_{1}w_{m}}{\rho},0\Bigg{\\}}\frac{r_{m}}{\lVert
r_{m}\rVert_{2}},\\\
r_{m}&=\tilde{\tilde{B}}_{j}^{\mathcal{G}_{m}}+H_{j}^{\mathcal{G}_{m}}\end{split}$
(41)

  5. 5.

$\tilde{E}_{jd}$ update: The $\tilde{E}$ component of the optimization problem
(30) is updated as

$\underset{\tilde{E}_{jd}}{\arg\min}\quad\lambda_{2}w_{d}\lVert\tilde{E}_{jd}\rVert_{2}+\tilde{H}_{jd}(-\tilde{E}_{jd})^{T}+\frac{\rho}{2}\lVert
B_{jd}-\tilde{E}_{jd}\rVert_{2}^{2},$ (42)

which can be simplified as

$\underset{\tilde{E}_{jd}}{\arg\min}\left[\lambda_{2}w_{d}\lVert\tilde{E}_{jd}\rVert_{2}+\frac{\rho}{2}\Bigg{\|}B_{jd}-\tilde{E}_{jd}+\frac{\tilde{H}_{jd}}{\rho}\Bigg{\|}_{2}^{2}\right].$
(43)

Further computations, by following the group lasso example will show that,

$\begin{split}\tilde{E}_{jd}&=\max\Bigg{\\{}\lVert
r_{d}\rVert_{2}-\frac{\lambda_{2}w_{d}}{\rho},0\Bigg{\\}}\frac{r_{d}}{\lVert
r_{d}\rVert_{2}},\\\ r_{d}&=B_{jd}+\tilde{H}_{jd}\end{split}$ (44)

  6. 6.

$H,\tilde{H},P$ and $O$ updates: The three Lagrange multipliers, associated
with the penalties in (30) are updated as follows;

$\begin{split}P_{d}&=P_{d}+B_{d}-Q_{d}\\\ O_{d}&=O_{d}+\tilde{B}_{d}-V_{d}\\\
H&=H+\tilde{\tilde{B}}-E\\\
\tilde{H}_{d}&=\tilde{H}_{d}+B_{d}-\tilde{E}_{d}\end{split}$ (45)

##  3 Simulations

This section entails the simulation studies for our proposed model. In the
first part, we show the results for the single-response problem. We compare
the ADMM pliable lasso algorithm to the original pliable lasso implementation,
which uses the block-wise coordinate descent (BCD) method. The second part of
the simulation studies is dedicated to the multi-response problem. We compare
our proposed method with the tree lasso and the original pliable lasso which
assumes that the multiple responses are independent. For the original pliable
lasso, we allowed each response to be modeled one at a time, since there is no
multi-response version implemented.

###  3.1 Simulation of a single-response problem

We generated data with $N=100$, $p=10$ or $p=50$, $K=3$ and standard Gaussian
independent predictors. The generative model is given as

$y=X_{1}(\beta_{1}\mathbf{1}+2Z_{3})+X_{2}\beta_{2}+X_{3}(\beta_{3}\mathbf{1}+2Z_{1})+X_{4}(\beta_{4}\mathbf{1}-2Z_{2})+.5\varepsilon.$
(46)

Here, $\beta=(2,-2,2,2,0,0,\ldots)$ and $Z$ is $N\times K$ matrix such that
$Z_{1},\ldots,Z_{K}\overset{iid}{\sim}\mathbb{N}(0,1)$ and
$\varepsilon\sim\mathbb{N}(0,1)$. The training, validation and test error
plots are shown in figure (1). We find that both the proposed approach and the
original plasso, which uses the coordinate descent method, give similar
results. We should note here that our main purpose for the single response
ADMM approach is to show how the algorithm can be used to effectively handle
the overlapping group penalty in the problem.

(a)

(b)

(c)

(d)

Figure 1: Plots for training (a), testing (b) and validation error (c-d) for
the two methods. Plot (c) and (d) are for the ADMM and BCD respectively. The
error is the mean squared error and “nlambda” is the length of the $\lambda$
vector. The numbers across the top of the validation error plots (c) and (d)
represent the numbers of non-zero main effects at the given $\lambda$ and the
vertical broken lines show the selected $\lambda$ from cross-validation.

Figures (2) and (4) show the results for the selected coefficients along the
$\lambda$ path for ADMM, and the BCD for $p=10$ and $p=50$, respectively,
after implementing them on the model (46). The figures (3) and (5) show the
matrix plot of the observed and the predicted coefficients of the interaction
effects. We can see that our proposed algorithm is able to identify the true
non-zero main effects as well as the non-zero interaction effects for model
(46) starting from a position where all $\hat{\beta}_{j}$s are zero.

(a)

(b)

Figure 2: Coefficient plots from the simulated experiment for $p=10$ for ADMM
(a) and BCD (b). The numbers across the top of the validation plots represent
the number of non-zero main effects at the given $\lambda$ and the vertical
broken lines show the selected $\lambda$ from cross-validation. The “X” symbol
indicates that a modifying term $Z\bm\theta_{j}$ has entered the model i.e the
presence of an interaction effect.

(a)

(b)

(c)

Figure 3: Matrix plot for the interaction coefficients $\theta_{jk}$ for
$X_{j}$ and $Z_{k}$ from the simulated experiment for $p=10$ for the simulated
values (a), ADMM estimates (b), and BCD estimates (c). The x-axis represents
the modifying variables $Z_{1},\ldots,Z_{K}$ and the y-axis represents the
$X_{1},\ldots,X_{p}$ to show the interactions $X_{j}\circ Z$.

(a)

(b)

Figure 4: Coefficient plots from the simulated experiment for $p=50$ for ADMM
(a) and BCD (b). The numbers across the top of the validation plots represent
the number of non-zero main effects at the given $\lambda$ and the vertical
broken lines show the selected $\lambda$ from cross-validation. The “X” symbol
indicates that a modifying term $Z\bm\theta_{j}$ has entered the model, i.e.
the presence of an interaction effect.

(a)

(b)

(c)

Figure 5: Matrix plot for main and interaction coefficients from the simulated
experiment for $p=50$ for the simulated values (a), ADMM estimates (b), and
BCD estimates (c). The x-axis represents the modifying variables
$Z_{1},\ldots,Z_{K}$ and the y-axis represents the $X_{1},\ldots,X_{p}$ to
show the interactions $X_{j}\circ Z$.

###  3.2 Simulation of a multi-response problem

To see the performance of the proposed method for the multi-response problem,
we simulated two different datasets; one with a weak hierarchical structure
within the response variables and the other with a strong hierarchical
structure. Both datasets were allowed to have both main effects $X$ and
interactions with a modifying variable $Z$.

####  3.2.1 Simulation 1

For the first simulation, we generated a multi-response problem with $D=6$
responses, $N=100$ observations, $p=500$ main effects and $K=4$ modifying
variables. Both $X$ and $Z$ were generated as standard Gaussian independent
predictors. We allowed each response to have five non-zero main coefficients
and also allowed four $X_{j}$s in each response to have interactions with four
$Z_{k}$s. The error term, $\mathbb{E}_{d}\sim\mathbf{N}(0,1)$ was added to
each response $d=1,\ldots,D$ and we generated the responses from

$Y=XB+\sum_{d}\underset{j}{\sum}(X_{j}\circ Z)\bm\theta_{jd}+\mathbb{E}.$ (47)

We averaged the results over 10 simulations and compared our work with that of
the tree lasso, which has been made available in R in the package ‘mixlasso’
version 0.1 available on Github (https://github.com/zhizuio/mixlasso) and the
pliable lasso (Tibshirani and Friedman, 2020) package (‘pliable’ version
1.1.1) also available in CRAN archive
(https://cran.r-project.org/src/contrib/Archive/pliable/). To allow the tree
lasso model to include the modifying variable, we allowed the modifier $Z$ to
be included without any penalty. This option has been made available in the
mixlasso R package by selecting method="tree-lasso" and num.nonpen=K.

(a)

(b)

Figure 6: The general structure of multi-response simulation 1 and the box
plot from the mean squared errors (MSE) on the test data. The dendrogram (a)
shows the hierarchical structure of the multivariate responses generated by
$\beta$, as estimated with complete-linkage hierarchical clustering using the
Euclidean distance. Below each node of the dendrogram is a number indicating
the response variable associated with the node. The MSE on the test data (b)
was averaged across the 10 simulations.

The structure of the response variables $Y$ is shown in figure (6(a)). The
summary of the results is shown in table (1), where we summarize the
sensitivity (the proportion of non-zero coefficients correctly estimated as
non-zeros), specificity (the proportion of zero-coefficients correctly
estimated as zeros), the total number of the estimated non-zero coefficients
(out of $p=500\times D=6$) and the test error measured as the mean squared
error on a separate test data with $N=500$, $p=500$ and $K=4$. It can be seen
that our proposed model performs best in this case in terms of the test error
(from figure (6(b))). In this setting with weak hierarchical structure in the
response, the pliable lasso model performed better than tree lasso, indicating
that the correctional structure in the responses is not so important. By
looking at the averaged absolute errors of estimated coefficients
($1/Dp\lVert\hat{\beta}-\beta\rVert_{1}$) in table (1), it can be seen that
the proposed model also performs best in terms of the accuracy of coefficients
matrix estimation and variable selection for this particular case.

Table 1: Results from the multi-response simulation 1 with weak hierarchical
structure in the response.

Model | $(1/Dp)\lVert\hat{\beta}-\beta\rVert_{1}$ |  Sensitivity1 |  Specificity2 |  Non-zero3 |  Test error (SD)4
---|---|---|---|---|---
Plasso | 0.021 | 1 | 0.763 | 733 | 19.693 (2.408)
Tree lasso | 0.066 | 1 | 0.142 | 2577 | 34.045 (1.802)
MADMMplasso | 0.006 | 1 | 0.991 | 237 | 5.050 (0.681)

  * 1

Sensitivity is the proportion of non-zero coefficients estimated as non-zeros.

  * 2

Specificity is the proportion of zero-coefficients estimated as zeros.

  * 3

The total number of non-zero coefficients in the model. We counted the
coefficients with at least two non-zero values across the 10 simulations.

Number of non-zero coefficients
$=\sum_{j=1}^{p}\sum_{d=1}^{D}\\{(\sum_{r=1}^{10}\bm 1_{\\{\beta_{jd}^{r}\neq
0\\}})\geq 2\\}$. Note that the selection is out of $p\times D=3000$ features
in total.

  * 4

The MSE on an independent test dataset. We include the standard deviation (SD)
across the 10 simulations.

Figure 7: The general structure of the multi-response simulation 1. Starting
from left, the simulated (“true”) coefficients, MADMMplasso estimates, plasso
estimates, and tree lasso estimates. For the true coefficients, we allowed
each response to have five non-zero coefficients as indicated by the red and
blue blocks. Below each block is a number indicating the response variable
associated with the block. Some responses were allowed to have non-zero
coefficients for the same variables so that we could obtain some correlations
among the responses.

Figure (7) shows the structure of the main effect coefficients ($\beta$) and
how they correlate across the responses. It can be seen that our proposed
model is able to recover the true relevant covariates for correlated responses
significantly better than other methods. The plasso model seem not to handle
the groups well since it not able to borrow statistical strength across
different responses. Even though the tree lasso does well in recovering the
groups, it seem to have a lot of false positives than our proposed model.

####  3.2.2 Simulation 2

For the second simulation in this section, we generated a multi-response
problem with $D=24$ responses (imagined to be responses of cancer cells to
drugs), $N=100$ observations (e.g. different cancer cell lines), $p=150$ or
$p=500$ (e.g. gene expression variables characterizing the different cancer
cell lines) and $K=4$ (e.g. cancer tissue type). The idea is to simulate a
drug sensitivity response model that includes interaction effects to allow the
cancer tissue types to modify the effects of gene expression features. $X$ was
generated as standard Gaussian independent predictors to represent the gene
expression features and the $Z$ was was generated using the binomial
distribution to represent the presence or absence of a certain tissue. The
interaction in this scenario was to allow for inclusion of tissue effect,
where we have assumed that certain drugs work best on specific tissues. The
error term, $\mathbb{E}_{d}\sim\mathbf{N}(0,1)$ was added to each response
$d=1,\ldots,D$ and we generated the drug sensitivity using the model

$Y=XB+\sum_{d}\sum_{j}(X_{j}\circ Z)\bm\theta_{jd}+\mathbb{E}.$ (48)

We generated $X$ by following the example in Zhao and Zucknick (2020).
However, since we have only one source of data in $X$, we decided to only use
the first part of their matrix. To do this, we generated $X$ from the
multivariate normal distribution with zero mean and non-diagonal $p\times p$
covariance matrix $\sum$. The covariance matrix $\sum$ has $p/b(\sigma)$ as
the covariance of features and the remaining part being set to zero. We set
$\sigma=0.4$, $b=10$, and the variance of each feature is one. We averaged the
results over 10 simulations and compared our work with that of the tree lasso
and the pliable lasso. To allow the tree lasso model to include the modifying
variable, we allowed the modifier $Z$ to be included without any penalty. This
option has been made available in the mixlasso R package selecting
method="tree-lasso" and num.nonpen=K.

(a)

(b)

(c)

Figure 8: The general structure of the multi-response simulation 2 and the box
plots from the mean squared error (MSE) on the test data. The dendrogram (a)
shows the hierarchical structure of the multivariate responses generated by
$\beta$, as estimated with complete-linkage hierarchical clustering using the
Euclidean distance. Below each node of the dendrogram is a number indicating
the response variable associated with the node. The MSE on the test data for
$p=150$ (b), and $p=500$ (c) were averaged across the 10 simulations.

The structure of the response variables $Y$ is shown in figure (8(a)). The
results are shown in table (2), where we summarize the sensitivity,
specificity, the total number of non-zero coefficients (out of $p\times D$)
and the test error measured as the mean squared error on a separate test
dataset. Our proposed model performs best in this case in terms of the test
error (in figures (8(b)) and (8(c))). In this setting with strong hierarchical
structure in the responses, the tree lasso performs better than the pliable
lasso, indicating the importance of the groups within the responses.

By looking at the the averaged absolute errors of estimated coefficients
($1/Dp\lVert\hat{\beta}-\beta\rVert_{1}$) in table (2), we see that the
proposed model also performs best in terms of the accuracy of coefficients
matrix estimation and variable selection for this particular case.

Table 2: Results from the multi-response simulation 2 with strong hierarchical
structure in the responses.

Model | $(1/Dp)\lVert\hat{\beta}-\beta\rVert_{1}$ |  Sensitivity1 |  Specificity2 |  Non-zero3 |  Test error (SD)4
---|---|---|---|---|---
$p=150$ |  |  |  |  |
Plasso | 0.034 | 1 | 0.446 | 2155 | 2.512 (0.181)
Tree lasso | 0.036 | 1 | 0.345 | 2483 | 2.072 (0.095)
MADMMplasso | 0.0299 | 1 | 0.727 | 2014 | 1.972 (0.112)
$p=500$ |  |  |  |  |
Plasso | 0.014 | 0.994 | 0.814 | 2514 | 4.57 (1.038)
Tree lasso | 0.023 | 1 | 0.360 | 7826 | 2.927 (0.163)
MADMMplasso | 0.010 | 1 | 0.912 | 1891 | 2.230 (0.116)

  * 1

Sensitivity is the proportion of non-zero coefficients estimated as non-zeros.

  * 2

Specificity is the proportion of zero-coefficients estimated as zeros.

  * 3

The total number of non-zero coefficients in the model. We counted the
coefficients with at least two non-zero values across the 10 simulations.

Number of non-zero coefficients
$=\sum_{j=1}^{p}\sum_{d=1}^{D}\\{(\sum_{r=1}^{10}\bm 1_{\\{\beta_{jd}^{r}\neq
0\\}})\geq 2\\}$. Note that the selection is out of $p\times D=3600$ (for
$p=150$) or $12000$ (for $p=500$) features in total.

  * 4

The MSE on an independent test dataset. We included the standard deviation
(SD) across the 10 simulations.

Figure 9: The general structure of the multi-response simulation 2. Starting
from left, the simulated (“true”) coefficients, MADMMplasso estimates, plasso
estimates, and tree lasso estimates. The non-zero coefficients in the true
data are represented by red blocks. Below each block is a number indicating
the response variable associated with the block. Some responses were allowed
to have non-zero coefficients for the same variables so that we could obtain
some correlations among the responses.

Figure (9) shows the structure of the main effect coefficients ($\beta$) and
how they correlate across the responses. It is clear that the tree lasso and
our proposed method are able to recover the true relevant covariates for
correlated responses significantly better than the pliable lasso even though
the tree lasso seems to have rather more false positives (also see table 2).

##  4 Genomics of drug sensitivity in cancer (GDSC) data analysis

We obtained this data set from the archives of the GDSC database
(ftp://ftp.sanger.ac.uk/pub4/cancerrxgene/releases/release-5.0/) (Yang et al.,
2012). The data set contains 97 cancer drugs tested on 498 cell lines. The
cell lines are representations of 13 cancer tissues. The drug responses were
summarized by the logarithm of the half maximal inhibitory concentration, i.e.
$\log(\text{IC}_{50})$, estimated from the dose-response curves that measured
drug response at specific drug concentrations in terms of normalized relative
cell viability. The data set also contains gene expression features and other
genomic data such as DNA copy numbers and mutations. In our experiment, we
preselected 2602 gene expression features that explain $50\%$ of the variation
as $X$ and represented the cancer types in a 12-column matrix of dummy
variables as $Z$. We further selected seven drugs for illustrative purpose.
The structure of the drugs is shown in figure (10(a)).

(a)

(b)

(c)

Figure 10: The general structure and the scatter plots for the seven drugs, as
estimated with complete-linkage hierarchical clustering using the Euclidean
distance (a), box plots of the overall mean squared error (MSE) on the test
data (b) and the MSE for the individual drugs on the test data (c). The
dendrogram (a) shows the hierarchical structure of the multivariate responses
generated by the drugs $Y$. Below each node of the dendrogram is the name of
the drug associated with the node. The MSE on the test data ((b) and (c)) were
averaged across the 10 repeated data splits.

We randomly split the data into $80\%$ training and $20\%$ testing data and a
high-dimensional multi-response regression model was fit using the training
data, where $Y$ was the $\log(\text{IC}_{50})$ values. We performed this
procedure repeatedly on 10 random splits and present the average of the
results. We show the efficiency our model by comparing with the tree lasso
(using the mixlasso R package) and the pliable lasso. To allow the tree lasso
to include the modifying variable, we allowed the cancer types $Z$ to be
included without any penalization. To use the pliable lasso, we modeled each
response $Y_{d},\quad d=1,\ldots,D$ separately. The summary of the results for
the 10 random splits is shown in table (3) and in figures (10(b)) and (10(c)).
We find that the proposed model performs best in terms of prediction error on
the test data. The proposed model is able to combine the information from the
various cancer types with corresponding gene expression features to determine
which drug will be suitable for each cancer type. It also borrows statistical
strength from the correlations that exist among the responses. This is very
essential in the area of drug discoveries where it important to identify
individual drugs or groups of drugs that could be helpful in treating specific
cancers and subtypes of cancers.

Table 3: Results from the GDSC data.

Model | Non-zero coefficients1 | Test error (SD)2
---|---|---
Plasso | 724 | 3.648 (0.270)
Tree lasso | 1016 | 3.404 (0.268)
MADMMplasso | 1424 | 3.227 (0.267)

  * 1

The number of non-zero coefficients in the model. We counted the coefficients
with at least two non-zero values across the 10 repeated data splits. Number
of non-zero coefficients $=\sum_{j=1}^{p}\sum_{d=1}^{D}\\{(\sum_{r=1}^{10}\bm
1_{\\{\beta_{jd}^{r}\neq 0\\}})\geq 2\\}$

Note that the selection is out of $p\times D=18844$ features in total.

  * 2

The MSE on an independent test data. We included the standard deviation (SD)
across the 10 repeated data splits.

In figure (11), we show the prediction performance of the three models. We can
see from the plots that all the models perform averagely well in terms of
predicting the individual response variables $Y$. However, by looking at
figure (10(c)) which gives the averages of the test errors for the individual
drugs we see that our proposed approach performs best in terms of prediction
error.

Further analyses on the proposed model are shown in figures (12) (for
individual drugs) and (13)111This was generated by using the R package
“EnrichIntersect” created by Zhao et al.,2022. (for all drugs together). These
plots show the selected interactions between gene expression features and
cancer types. We show the relationship between the identified genes (those
considered to interact with cancer types), the cancer types and the associated
drugs, identified by the model to be linked to specific cancer types. We note
here that the gene expression features shown in the plots are only those that
interact with at least one cancer type. It can be seen that not all the seven
drugs were identified to be linked to specific cancer types. Specifically,
methotrexate was not found by the model to be associated with a specific
cancer type as no interactions between cancer type and gene expression
features were selected in the MADMMplasso model for this drug. This is as
expected as methotrexate is a generic chemotherapy agent that does not target
specific molecular features or pathways in the cell. One of the interesting
drugs identified is nilotinib which is associated with blood cancers,
specifically it is a is a Bcr-Abl tyrosine kinase inhibitor which targets
chronic myelogenous leukemias (CML) with the Philadelphia chromosome. A common
gene associated with blood cancers as identified by all the three models is
the suppressor of cytokine signaling 2 (SOCS2). In figure (12), we can see
that this gene is seen to interact with the blood cancer type. In fact,SOCS2
is highly expressed in many of the blood cancer cell lines in this data set.
Studies have shown that SOCS2 is involved in the signal transduction cascades
in CML cells (Schultheis et al., 2002). From the GDSC data set, we also
realized that most of the cell lines with high expression of SOCS2 and having
low score of $\log(\text{IC}_{50})$ for nilotinib are CML tumors. This
provides evidence that the proposed model can help in identifying genes that
are associated with diseases and their corresponding drugs in drug
discoveries.

We see from figure (13) that the total number of identified interactions is
quite small. Most interactions are observed for the blood cancer type. Blood
cancers are biologically quite different from solid tumours, and the
interaction effects allow the MADMMplasso model to modify the gene expression
effect estimates to allow different effects for blood cancer compared to the
solid tumours. Most interaction effects with the blood cancer type are
identified for the nilotinib and axitinib drug response variables. This makes
sense as both these drugs are Bcr-Abl tyrosine kinase inhibitor which are only
applied to chronic myelogenous leukemias with Bcr-Abl mutations.

Figure 11: Scatter plots on the observed drugs responses versus the predicted
drug responses for MADMMplasso (a), plasso (b) and tree lasso (c). Note that
these results were obtained for one of the random data splits. The x-axis
represents the observed responses and and y-axis represents the predicted
responses.

(a)

(b)

(c)

(d)

(e)

Figure 12: Chord diagrams for the interaction effects between gene expression
features and cancer types for Nilotinib (a), RDEA119 (b), CI-1040 (c),
Axitinib (d) and PD-0325901 (e). The upper chord represents the cancer types
and the lower chord represents the gene expression features. The width of the
line connecting the cancer type to a gene indicates the strength of the
interaction. The selection was done for the coefficients with at least two
non-zero values across the 10 random data splits.  Figure 13: The interaction
plot (Sankey diagram) for the GDSC dataset. Each gene expression feature in
the figure has at least one corresponding cancer type and at least one
corresponding drug and each line connecting the three variables indicates a
selected interaction effect between the gene expression feature and the cancer
type for the given drug response variable. The lines have been been grouped
according to the drugs. Note that the selection into the figure was done for
the coefficients with at least two non-zero values across the 10 random data
splits.

##  5 Conclusion

In this paper, we proposed a novel regularized regression method, which we
named the MADMMplasso, to represent a multi-response pliable lasso regression
model in combination with an ADMM algorithm for model fitting. This method is
able to find covariates $X$ and their corresponding interactions with
modifying variables $Z$, with some joint association with multiple related
responses. The joint association takes into consideration the correlational
structure in the response variables, assumed to form some overlapping groups,
and specifically in this paper a hierarchical tree structure. We allowed the
interaction term to be included in an asymmetrical weak hierarchical manner by
first considering whether the corresponding main term is in the model. For
optimization, we implemented an ADMM algorithm that allowed us to model the
overlapping groups in a simple way. The results from the simulations and from
the application to a pharmacogenomic screening dataset show that the ability
of the proposed method to account for correlated responses and to include
interaction effects can result in a clear improvement in prediction and
variable selection performance. An implementation of the package is available
on Github through the link https://github.com/ocbe-uio/MADMMplasso and the
version used throughout this paper is 1.0.0.

## Acknowledgements

This work received funding from the European Union’s Horizon 2020 Research and
Innovation program, under the Marie Skłodowska-Curie Actions Grant, agreement
No. 801133 (Scientia fellowship), and under grant agreement No. 847912
(“RESCUER”).

## References

  * Anbari et al. [2014] M. E. Anbari, S. Alam, and H. Bensmail.  Cofadmm: A computational features selection with alternating direction method of multipliers.  _Procedia Computer Science_ , 29:821–830, 2014.  ISSN 1877-0509.  https://doi.org/10.1016/j.procs.2014.05.074.  URL https://www.sciencedirect.com/science/article/pii/S1877050914002518.  2014 International Conference on Computational Science.
  * Asenso and Zucknick [2022] T. Q. Asenso and M. Zucknick.  MADMMplasso, 12 2022.  URL https://github.com/ocbe-uio/MADMMplasso.
  * Asenso et al. [2022a] T. Q. Asenso, P. Wang, and H. Zhang.  Pliable lasso for the support vector machine.  _Communications in Statistics - Simulation and Computation_ , 0(0):1–13, 2022a.  10.1080/03610918.2022.2032160.  URL https://doi.org/10.1080/03610918.2022.2032160.
  * Asenso et al. [2022b] T. Q. Asenso, H. Zhang, and Y. Liang.  Pliable lasso for the multinomial logistic regression.  _Communications in Statistics - Theory and Methods_ , 51(11):3596–3611, 2022b.  10.1080/03610926.2020.1800041.  URL https://doi.org/10.1080/03610926.2020.1800041.
  * Boyd et al. [2011] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, et al.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  _Foundations and Trends® in Machine learning_ , 3(1):1–122, 2011.
  * Deng et al. [2013] W. Deng, W. Yin, and Y. Zhang.  Group sparse optimization by alternating direction method.  In _Wavelets and Sparsity XV_ , volume 8858, pages 242–256. SPIE, 2013.
  * Du and Tibshirani [2018] W. Du and R. Tibshirani.  A pliable lasso for the cox model, 2018.
  * Goldstein et al. [2014] T. Goldstein, B. O’Donoghue, S. Setzer, and R. Baraniuk.  Fast alternating direction optimization methods.  _SIAM Journal on Imaging Sciences_ , 7(3):1588–1623, 2014.  10.1137/120896219.  URL https://doi.org/10.1137/120896219.
  * Han and Yuan [2012] D. Han and X. Yuan.  A note on the alternating direction method of multipliers.  _Journal of Optimization Theory and Applications_ , 155(1):227–238, 2012.
  * Haris et al. [2016] A. Haris, D. Witten, and N. Simon.  Convex modeling of interactions with strong heredity.  _Journal of Computational and Graphical Statistics_ , 25(4):981–1004, 2016.  10.1080/10618600.2015.1067217.  URL https://doi.org/10.1080/10618600.2015.1067217.  PMID: 28316461.
  * Hastie et al. [2015] T. Hastie, R. Tibshirani, and M. Wainwright.  _Statistical Learning with Sparsity: The Lasso and Generalizations_.  Chapman & Hall/CRC, 2015.  ISBN 1498712169, 9781498712163.
  * Humayoo and Cheng [2019a] M. Humayoo and X. Cheng.  Parameter estimation with the ordered $\ell$2 regularization via an alternating direction method of multipliers.  _Applied Sciences_ , 9(20):4291, oct 2019a.  10.3390/app9204291.  URL https://doi.org/10.3390%2Fapp9204291.
  * Humayoo and Cheng [2019b] M. Humayoo and X. Cheng.  Parameter estimation with the ordered $l_{2}$ regularization via an alternating direction method of multipliers.  _Applied Sciences_ , 9(20), 2019b.  ISSN 2076-3417.  10.3390/app9204291.  URL https://www.mdpi.com/2076-3417/9/20/4291.
  * Kim and Xing [2012] S. Kim and E. P. Xing.  Tree-guided group lasso for multi-response regression with structured sparsity, with an application to eqtl mapping.  _The Annals of Applied Statistics_ , 6(3):1095–1117, 2012.  ISSN 19326157, 19417330.  URL http://www.jstor.org/stable/41713516.
  * Li et al. [2015] Y. Li, B. Nan, and J. Zhu.  Multivariate sparse group lasso for the multivariate multiple linear regression with an arbitrary group structure.  _Biometrics_ , 71(2):354–363, 2015.
  * Qin and Goldfarb [2012] Z. T. Qin and D. Goldfarb.  Structured sparsity via alternating direction methods.  _Journal of Machine Learning Research_ , 13(5), 2012.
  * Schultheis et al. [2002] B. Schultheis, M. Carapeti-Marootian, A. Hochhaus, A. Weisser, J. M. Goldman, and J. V. Melo.  Overexpression of SOCS-2 in advanced stages of chronic myeloid leukemia: possible inadequacy of a negative feedback mechanism.  _Blood_ , 99(5):1766–1775, 03 2002.  ISSN 0006-4971.  10.1182/blood.V99.5.1766.
  * Tibshirani [1996] R. Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society. Series B_ , 58(1), 1996.
  * Tibshirani and Friedman [2020] R. Tibshirani and J. Friedman.  A pliable lasso.  _Journal of Computational and Graphical Statistics_ , 29(1):215–225, 2020.  10.1080/10618600.2019.1648271.  URL https://doi.org/10.1080/10618600.2019.1648271.
  * Yang et al. [2012] W. Yang, J. Soares, P. Greninger, E. J. Edelman, H. Lightfoot, S. Forbes, N. Bindal, D. Beare, J. A. Smith, I. R. Thompson, S. Ramaswamy, P. A. Futreal, D. A. Haber, M. R. Stratton, C. Benes, U. McDermott, and M. J. Garnett.  Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells.  _Nucleic Acids Research_ , 41(D1):D955–D961, 11 2012.  ISSN 0305-1048.  10.1093/nar/gks1111.  URL https://doi.org/10.1093/nar/gks1111.
  * Ye et al. [2011a] G. Ye, Y. Chen, and X. Xie.  Efficient variable selection in support vector machines via the alternating direction method of multipliers.  In G. Gordon, D. Dunson, and M. Dudík, editors, _Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics_ , volume 15 of _Proceedings of Machine Learning Research_ , pages 832–840, Fort Lauderdale, FL, USA, 11–13 Apr 2011a. PMLR.  URL https://proceedings.mlr.press/v15/ye11a.html.
  * Ye et al. [2011b] G.-B. Ye, Y. Chen, and X. Xie.  Efficient variable selection in support vector machines via the alternating direction method of multipliers.  In _Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics_ , pages 832–840. JMLR Workshop and Conference Proceedings, 2011b.
  * Zhao and Zucknick [2020] Z. Zhao and M. Zucknick.  Structured penalized regression for drug sensitivity prediction.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 69(3):525–545, 2020.
  * Zhao et al. [2022] Z. Zhao, M. Zucknick, and T. Aittokallio.  EnrichIntersect: an R package for custom set enrichment analysis and interactive visualization of intersecting sets.  _Bioinformatics Advances_ , 2(1), 09 2022.  ISSN 2635-0041.  10.1093/bioadv/vbac073.  URL https://doi.org/10.1093/bioadv/vbac073.


#  Effect of electron-phonon interactions on Raman line at ferromagnetic
ordering

L.A. Falkovsky  Landau Institute for Theoretical Physics, Chernogolovka 142432
Verechagin Institute of the High Pressure Physics, Troitsk 142190

###### Abstract

The theory of Raman scattering in half-metals by optical phonons interacting
with conduction electrons is developed. We evaluate the effect of electron-
phonon interactions at ferromagnetic ordering in terms of the Boltzmann
equation for carriers. The chemical potential is found to decrease with
temperature decreasing. Both the linewidth and frequency shift exhibit a
dependence on temperature.

###### pacs:

42.50.Nn 63.20.-e 75.30.Ds 78.30.-j

##  I Introduction

Recently, the Raman scattering in the half-metallic CoS2 was studied La  in
the wide temperature region. The $\omega=400$ cm-1 Raman line, observed
previously at room temperature in Ref. AP ; ZST , demonstrates a particular
behavior nearby the ferromagnetic transition at $T_{c}=122$ K. The unusual
large Raman linewidth and shift of the order of 10 cm-1 were observed. The
reflectivity singularities of CoS2 were explained in Ref. YMM  by the
temperature variation of the electronic structure. Another example of the
electron-phonon interactions is given in Ref. LM  in order to explain the
phonon singularity at the $\Gamma$ point in graphene. The electron-phonon
interactions should be considered as well in the interpretation of the
observed Raman scattering around the Curie temperature.

Thermal broadening of phonon lines in the Raman scattering is usually
described in terms of three-phonon anharmonicity, i.e. by the decay of an
optical phonon with a frequency $\omega$ in two phonons. The simplest case
when the final state has two acoustic phonon from one branch (the Klemens
channel) was theoretically studied by Klemens Kl , who obtained the
temperature dependence of the Raman linewidth. The corresponding lineshift was
considered in Refs. BWH ; MC . This theory was compared in works BWH ; MC ;
DBM  with experimental data for Si, Ge, C, $\alpha-$Sn. A model was also
considered with the phonons in the final state from different branches. It was
found that anharmonic interactions of the forth order should be disregarded at
high temperatures $T>300$ K.

The situation is more complicated in substances with magnetic ordering. The
interaction of phonons with magnons in antiferromagnets was discussed in the
review article GZ  and more recently in the analysis of the thermal
conductivity MOZ , the spin Seebeck effect UTH ; JYM , high-temperature
superconductivity NMH , and optical spectra KKP . The magnon-phonon
interaction results in the magnon damping Wo , however, no effect for phonons
was observed. The influence of antiferromagnetic ordering is considered in
Ref. DPV , where the line shift was only calculated. Damping of the optical
phonons was found MPB  to become large in the rare-earth Gd and Tb below the
Curie temperature achieving a value of 15 cm-1, which is much greater than the
three-phonon interaction effect.

Despite attracting considerable interest for half a century since the
pioneering work by Fröhlich, the problem of electron-phonon interaction is
still far from being solved. Migdal Mi  developed a consistent many-body
approach based on the Fröhlich Hamiltonian for interaction of electrons with
acoustic (sound) phonons. As Migdal showed (”the Migdal theorem”), the vertex
corrections for acoustic phonons are small by the adiabatic parameter
$\sqrt{m/M}$, where $m$ and $M$ are the electron and ion masses, respectively.
The theory described correctly the electronic lifetime, renormalization of the
Fermi velocity $v_{F}$ and acoustic phonon attenuation but resulted in a
strong renormalization of the sound velocity $\tilde{s}=s(1-2\lambda)^{1/2}$,
where $\lambda$ is the dimensionless coupling constant. For sufficiently
strong electron-phonon coupling $\lambda\to 1/2$, the phonon frequency
approached to zero marking an instability point of the system. Instead, one
would intuitively expect the phonon renormalization to be weak along with the
adiabatic parameter.

This discrepancy was resolved by Brovman and Kagan BK  almost a decade later
(see also G ). They demonstrated the shortcomings of the Fröhlich model that
gave an anomalously large phonon renormalization. Employing the
Born–Oppenheimer (adiabatic) approximation (see, e.g., BHK ), they found that
there are two terms in the second order perturbation theory, which compensate
each other making a result small by the adiabatic parameter. Namely, when
calculating the phonon self-energy function $\Pi(\omega,k)$ with help of the
diagram technique, one should eliminate an adiabatic contribution of the
Fröhlich model by subtracting $\Pi(\omega,k)-\Pi(0,k)$.

The interaction of electrons with optical phonons was first considered by
Engelsberg and Schrieffer ES  within Migdal’s many-body approach for
dispersionless phonons. They predicted a splitting of the optical phonon at
finite wavenumbers $k$ into two branches. Ipatova and Subashiev IS  calculated
later on the optical phonon attenuation in the collisionless limit and pointed
out that the Brovman-Kagan renormalization should be carried out for optical
phonons in order to obtain correct phonon renormalization. In the paper AS ,
Alexandrov and Schrieffer corrected the calculational error of Ref. ES  and
argued that no splitting was found in fact. Instead, they predicted an
extremely strong dispersion of optical phonons, $\omega_{k}=\omega_{0}+\lambda
v_{F}^{2}k^{2}/3\omega_{0}$, due to the coupling to electrons. The large
phonon dispersion is a typical result of Migdal’s theory AGD  using the
Frölich Hamiltonian. No such a dispersion has ever been observed
experimentally. The usual dispersion of optical phonons in metals has the
order of the sound velocity. Reizer R  stressed the importance of screening
effects which should be taken into account. The works AS ; R  are limited to
the case of collisionless both electron and phonon systems. Moreover, only the
phonon renormalization was considered with no results available for the
attenuation of optical phonons.

A different from many-body technique semiclassical approach based on the
Boltzmann equation and the equations of the theory of elasticity was developed
in the papers by Akhiezer, Silin, Gurevich, Kontorovich, and many others (we
refer the reader to the review Kon ). This approach was compared with various
experiments, such as attenuation of sound waves, effects of strong magnetic
fields, crystal anisotropy, and sample surfaces on the sound attenuation, and
so on. It can be applied to the problem of the electron–optical-phonon
interaction MF  as well.

In the previous paper Fal , we have developed a quantum theory for the optical
phonon attenuation and shift induced by the interband electron transitions and
tuned with a temperature variation. Now we consider the optical phonon
renormalization as a result of the electron-phonon interaction taking into
account ferro-magnetic ordering. We argue that the reasonable phonon damping
and shift can be obtained using the semiclassical Boltzmann equation for
electrons and the motion equation of phonons coupled by the deformation
potential.

##  II Electron-phonon interactions at ferromagnetic ordering

Figure 1: (Color online) Proposed band scheme for two electron spin
projections.

We assume that the electron bands in CoS2 have a form shown in Fig. 1. The
ferromagnetic ordering results in the spin splitting $\mu H_{e}$ of the
unfilled half-metallic band

$\begin{array}[]{c}\displaystyle{\varepsilon_{\uparrow}({\bf
p})=\frac{p^{2}}{2m^{\vee}}}-\mu
H_{e}\quad\text{and}\quad\displaystyle{\varepsilon_{\downarrow}({\bf
p})=\frac{p^{2}}{2m^{\vee}}}+\mu H_{e}\end{array}$ (1)

in the effective Weiss field $H_{e}$. While the temperature decreases, the
magnetization, determined in the mean field approximation as

$M=M_{0}\sqrt{1-(T/T_{c})^{2}}\,,$ (2)

appears according to experimental data in CoS2 at $T_{c}=122$ K approximately,
and the spin splitting is proportional to the magnetization.

We write the interaction of electrons with the optical phonon $u_{i}$ as the
deformation potential

$H_{int}=\frac{u_{i}}{N}\sum_{s}\int\frac{d^{3}{\bf
p}}{(2\pi\hbar)^{3}}\zeta_{i}({\bf p})f({\bf p})\,,$ (3)

where $N\sim 1/a^{3}$ is a number of cells in the volume unit and $a$ is the
interatomic distance. For the acoustic phonon – electron interaction, we
should substitute the strain tensor $u_{ij}$ instead of the displacement
$u_{i}$ in order to satisfy the translation symmetry of the lattice.

Figure 2: (Color online) Calculated temperature dependence of the carrier
concentration for the spin up and spin down (relative to the total
concentration at temperatures above the temperature of ferromagnetic
ordering), and the dependence of the chemical potential (dash-dotted line).

The Boltzmann equation for the nonequilibrium part of the distribution
function $f({\bf p})$ has the form

$[-i(\omega-{\bf k\cdot v})+\tau^{-1}]f({\bf p})=-\frac{\partial
f_{0}}{\partial\varepsilon}[e{\bf v\cdot E}-i\omega u_{i}\zeta_{i}({\bf
p})]\,,$ (4)

where $f_{0}$ is the equilibrium distribution function. We omit in the
Boltzmann equation (4) the spin index $s$, which determines all the electron
parameters. The electron collision frequency $\tau^{-1}$ takes into account
the collisions with impurities and phonons. The collision frequency is
calculated for CoS2 in the Debye model with the temperature $T_{D}=500$ K. One
can see from Eq.(4) that the condition

$<\zeta_{i}>=0\,,$

have to be satisfied in order to obey the current continuity equation, where
the brackets mean averaging over the Fermi surface for temperatures
$T\ll\varepsilon_{F}$.

In the ferromagnetic phase while the temperature changes, the carriers
overflow from one spin state in another, but the total number of carriers

$N=\sum_{s}\int\frac{d^{3}{\bf p}}{(2\pi\hbar)^{3}}f_{0}(\varepsilon_{s})$ (5)

remains to be constant. This condition determines the chemical potential and
the carrier concentration with the spin up and spin down, shown in Fig.2. All
figures correspond here and what follows to the carrier concentration
$N=10^{21}$ cm-3 in the considered band with the Fermi energy $\mu=0.36$ eV
above the Curie temperature.

Let us write the motion equation for the phonon mode in a form

$(\omega_{0}^{2}-\omega^{2})u_{i}=\frac{QE_{i}}{M}-\frac{1}{M}\frac{\partial
H_{int}}{\partial u_{i}}\,,$ (6)

where $M$ is the reduced ion mass of the cell, $Q$ is the charge corresponding
to the optical vibration, and $\omega_{0}$ is the frequency of the considered
mode. Here, the last term represents the electron-phonon interaction. Using
the Boltzmann equation (4), we rewrite this term as follows

$-\frac{1}{M}\frac{\partial H_{int}}{\partial
u_{i}}=-\frac{u_{i}}{MN}\sum_{s}\int\frac{\omega\tau\zeta^{2}_{i}({\bf
p})}{\omega\tau+i}\left(-\frac{\partial
f_{0}}{\partial\varepsilon}\right)\frac{d^{3}{\bf p}}{(2\pi\hbar)^{3}}\,.$ (7)

The term with electric field in the Boltzmann equation disappears in the
integration over ${\bf p}$ due to the velocity inversion ${\bf
v}\rightarrow-{\bf v}$. The term with the wave vector ${\bf k}$ has to be
omitted for the Raman phonon, as the vector ${\bf k}$ is determined in this
case by the laser frequency $\omega_{i}$ and the optical phonon frequency
satisfies the condition $\omega\gg\omega_{i}v/c$.

Figure 3: (Color online) Calculated shift of the Raman line $\omega=400$ cm-1
due to the electron-phonon interaction, and the lineshift in the Klemens
channel (dashed line).

The electric field does not excited in the TO vibrations. Therefore, supposing
$E=0$ and integrating over the energy $\varepsilon$ instead of $p$, we find
from Eqs. (6) and (7) the lineshift $\delta\omega$ and linewidth
$\delta\Gamma$ determined by the electron-phonon interaction as

$\delta\omega_{TO}-i\delta\Gamma_{TO}=\frac{1}{2MN}\sum_{s}\int\frac{\tau(\omega\tau-i)\zeta^{2}({\bf
p})dS}{(\omega^{2}\tau^{2}+1)v(2\pi\hbar)^{3}}|_{\varepsilon=\varepsilon_{F}}\,,$
(8)

where $dS$ is an element of the Fermi surface, $v$ is the Fermi velocity.
Estimating $S=4\pi p_{F}^{2}$, $\zeta({\bf p})\sim\varepsilon_{0}/a$ and
$\varepsilon_{0}^{2}\sim\omega^{2}M/m$, where $\varepsilon_{0}\sim 3eV$ is the
typical electron energy in metals, we obtain

$\delta\omega_{TO}-i\delta\Gamma_{TO}\sim\frac{ap_{F}\tau\omega^{2}_{TO}}{2\pi^{2}\hbar(\tau\omega_{TO}+i)}\,.$

The equations (6) and (7) allow to express the phonon displacement $u$ in
terms the electric field $E$ and to calculate the phonon contribution $uNQ$
into the polarization. We find the total dielectric permittivity, adding the
contributions $\varepsilon_{\infty}$ of the filled bands

$\begin{array}[]{c}\varepsilon(\omega)={\displaystyle\varepsilon_{\infty}-\frac{4\pi
e^{2}}{3\omega}\sum_{s}\int\frac{\tau vdS}{(\omega\tau+i)(2\pi\hbar)^{3}}}+\\\
{\displaystyle\frac{4\pi
NQ^{2}}{M}\left[\omega_{0}^{2}-\omega^{2}+\frac{\omega\tau}{MN}\sum_{s}\int\frac{\zeta^{2}({\bf
p})dS}{(\omega\tau+i)v(2\pi\hbar)^{3}}\right]^{-1}}.\end{array}$ (9)

Figure 4: (Color online) Calculated temperature dependence of the width for
the Raman $\omega=400$ cm-1 line at the ferromagnetic ordering, and the
linewidth in the Klemens channel (dashed line).

The frequency of the longitudinal phonon mode is determined by the condition
$\varepsilon(\omega)=0$. In the absence of free carriers, one finds the
frequency of the LO mode as follows

$\omega_{LO}^{2}=\omega_{0}^{2}+\omega_{pi}^{2},$

where $\omega_{pi}^{2}=4\pi NQ^{2}/M\varepsilon_{\infty}$ is the ion plasma
frequency squared.

Using Eq. (9), we find the LO frequency in the presence of carriers as

$\omega^{2}_{LO}-\omega_{0}^{2}=\frac{\omega}{(2\pi\hbar)^{3}MN}\sum_{s}\int\frac{\tau\zeta^{2}({\bf
p})dS}{(\omega\tau+i)v}-\frac{\omega_{pi}^{2}}{\omega_{pe}^{2}}\omega(\omega+i\tau^{-1})\,,$
(10)

where the electron plasma frequencies squared

$\omega_{pe}^{2}=\frac{4\pi
e^{2}}{3\varepsilon_{\infty}}\sum_{s}\int\frac{vdS}{(2\pi\hbar)^{3}}\,$ (11)

is supposed to be large in comparison with $\omega_{pi}^{2}$. We can put also
$\omega=\omega_{LO}$ in the right-hand side of Eq. (10). Here the last term
describes the electric field screened by the free carriers. The main role
plays the first term, which coincides with the result for the TO mode, Eq.
(8), shown in Figs. 3 and 4, the results for the Klemens channel are taken
from Ref. Fal .

We should emphasize that the temperature dependence of the linewidth and
shift, Eq. (8), is determined mainly by the electron collision rate
$\tau^{-1}$ involving also, for instance, in the dc conductivity. Thus, for a
cubic crystal, the dc conductivity, i.e. at $\omega=k=0$, writes

$\sigma=\sum_{s}\frac{e^{2}}{3(2\pi\hbar)^{3}}\int\tau vdS\,.$

The details of the electron density of states and of the deformation potential
are responsible for peculiarities of the Raman line temperature dependence.

##  III summary

The Klemens formula describes the optical phonon width due to three-phonon
anharmonic interactions. The corresponding lineshift matches with the
linewidth. In such ferromagnets as CoS2 with the low Curie temperature, these
interactions are found to be too weak to describe quantitatively the
experimental data and to explain the very large Raman linewidth and shift.
Therefore, we propose the mechanism of the electron-phonon interaction
attended with the effect of the ferromagnetic ordering on the electron bands.
The deformation potential couples together the Boltzmann equation for
electrons and the motion equation for phonons producing the renormalization of
the phonon frequency. The corresponding Raman line width and shift are in
agreement with experiments in Ref. La .

##  IV acknowledgments

The author thank S. Lyapin and S. Stishov for information on their experiments
prior the publication, and A. Varlamov for useful discussions. This work was
supported by the Russian Foundation for Basic Research (grant No.
13-02-00244A) and the SIMTECH Program, New Centure of Superconductivity:
Ideas, Materials and Technologies (grant No. 246937).

## References

  * (1) S.G. Lyapin, A.N. Utyuzh, A.E. Petrova, A.P. Novikov, T.A. Lograsso, and S.M. Stishov, arXiv:1402.5785.
  * (2) E. Anastassakis and C. Perry, J. Chem. Phys. 64, 3604 (1976).
  * (3) L. Zhu, D. Susac, M. Teo, K.C. Wong, P.C. Wong, R.R. Parsons, D. Bizzotto, K.A.R. Mitchel, and S.A. Campbell, J. Catal. 258, 235 (2008).
  * (4) R. Yamamoto, A. Machida, Y. Moritomo, and A. Nakamura, Phys. Rev. B 59, R7793 (1999).
  * (5) S. Piscanec, M. Lazzeri, Francesco Mauri, A.C. Ferrari, and J. Robertson, Phys. Rev. Lett. 93, 185503 (2004); M. Lazzeri, S. Piscanec, Francesco Mauri, A.C. Ferrari, and J. Robertson, Phys. Rev. B 73, 155426 (2006); M. Lazzeri and Francesco Mauri, Phys. Rev. Lett. 97, 266407 (2006).
  * (6) P.G. Klemens, Phys. Rev. 148, 845 (1966).
  * (7) M. Balkanski, R.F. Wallis, and E. Haro, Phys. Rev. B 28, 1928 (1983).
  * (8) J. Menéndez and M. Cardona, Phys. Rev. B 29, 2051 (1984).
  * (9) A. Debernardi, S. Baroni, and E. Molinari, Phys. Rev. Lett. 75, 1819 (1995).
  * (10) G. Güntherodt and R. Zeyer, Light Scattering in Solids, vol. 4 (1984).
  * (11) M. Montagnese, M. Otter, X. Zotos, D.A. Fishman, N. Hlubek, O. Mityashkin, C. Hess, R. Saint-Martin, S. Singh, A. Revcolevschi, and P.H.M. van Loosdrecht, Phys. Rev. Lett. 110, 147206 (2013).
  * (12) S. Uchida, S. Takahashi, K. Harri, J. Ieda, W. Koshibae, K. Ando, S. Maekawa, and E. Saitoh, Nature (London) 455, 778 (2008).
  * (13) C.M. Jaworski, J. Yang, S. Mack, D.D. Awschalom, R.C. Myers, and J.R. Heremans, Phys. Rev. Lett. 106, 186601 (2011).
  * (14) F. Nori, R. Merlin, S. Haas, A. Sandvik, and E. Dagotto, Phys. Rev. Lett. 75, 553 (1995).
  * (15) S.A. Klimin, A.B. Kuzmenko, M.N. Popova, B.Z. Malkin, and I.V. Telegina, Phys. Rev. B 82, 174425 (2010).
  * (16) L.M. Woods, Phys. Rev. B 65, 014409 (2001).
  * (17) D.M. Djokic, Z.V. Popovic, F.R. Vukajlovic, Phys. Rev. B 77, 014305 (2008).
  * (18) A. Melnikov, A. Povolotskiy, and U. Bovensiepen, Phys. Rev. Lett. 100, 247401 (2008).
  * (19) A.B. Migdal, Sov. Phys. JETP 7, 996 (1958).
  * (20) E.G. Brovman and Yu. Kagan, Sov. Phys. JETP 25, 365 (1967).
  * (21) B.T. Geilikman, J. Low Temp. Phys. 4, 189, (1971).
  * (22) M. Born and Kung Huang, Dynamical Theory of Crystal Lattices (Oxford University Press, New York, 1954).
  * (23) S. Engelsberg and J.R. Schrieffer, Phys. Rev. 131, 993 (1963).
  * (24) I.P. Ipatova and A.V. Subashiev, Sov. Phys. JETP 39, 349 (1974).
  * (25) A.S. Alexandrov and J.R. Schrieffer, Phys. Rev. B 56, 13731 (1997).
  * (26) A.A. Abrikosov, L.P. Gor’kov, and I.Ye. Dzyaloshinskii, Methods of Quantum Field Theory in Statistical Physics, (Prentice-Hall, Englewood Cliffs, NJ, 1963).
  * (27) M. Reizer, Phys. Rev. B 61, 40 (2000).
  * (28) V.M. Kontorovich, Sov. Phys. Uspekhi 27 (2), 134 (1984).
  * (29) L.A. Falkovsky and E.G. Mishchenko, Phys. Rev. B 51, 7239 (1995); L.A. Falkovsky, Phys. Rev. B 66, 020302(R) (2002).
  * (30) L.A. Falkovsky, Phys. Rev. B 88, 155135 (2013).