Huber group LASSO
Given m datasets, \stackrel{\u0303}{X}\left(1\right),\dots ,\stackrel{\u0303}{X}\left(m\right), satisfying model (5), to ensure that all A(k)'s have the same structure, elements of A(k)'s on the same position are grouped together and can be inferred by the group LASSO,
{\displaystyle \underset{A\left(k\right)}{\text{min}}}{\displaystyle \sum _{i=1}^{p}}{\displaystyle \sum _{k=1}^{m}}{w}_{k}{\displaystyle \sum _{j=1}^{{n}_{k}}}{\left({y}_{ij}\left(k\right){A}_{i}{\left(k\right)}^{T}{x}_{j}\left(k\right)\right)}^{2}+\lambda {\displaystyle \sum _{i=1}^{p}}{\displaystyle \sum _{\ell =1}^{p}}\sqrt{{a}_{i\ell}{\left(1\right)}^{2}+\dots +{a}_{i\ell}{\left(m\right)}^{2}},
where A_{
i
}(k)^{T} is the i th row of the matrix A(k) and x_{
j
} (k) is the j th column of the matrix X(k). w_{
k
} is the weight for the k th dataset, which can be assigned by experience. In this study, we choose w_{
k
} = n_{
k
} / Σn_{
k
} i.e., the more observations the dataset has, the higher weight it is assigned with. The penalty term in (6) takes advantage of the sparse nature of GRNs and has the effect making the grouped parameters to be estimated either all zeros or all nonzeros [10], i.e., a_{
iℓ
}(k)'s, k = 1, . . . , m, become either all zeros or all nonzeros. Therefore, a consistent network topology can be obtained from the group LASSO method. λ is a tuning parameter which controls the degree of sparseness of the inferred network. The larger the value of λ, the more grouped parameters become zeros.
To introduce robustness, we consider using the Huber loss function instead of the squared error loss function and propose the following Huber group LASSO method
{\displaystyle \underset{A\left(k\right)}{\text{min}}}{\displaystyle \sum _{i=1}^{p}}{\displaystyle \sum _{k=1}^{m}}{w}_{k}{\displaystyle \sum _{j=1}^{{n}_{k}}}{H}_{\delta}\left({y}_{ij}\left(k\right){A}_{i}{\left(k\right)}^{T}{x}_{j}\left(k\right)\right)+\lambda {\displaystyle \sum _{i=1}^{p}}{\displaystyle \sum _{\ell =1}^{p}}\sqrt{{a}_{i\ell}{\left(1\right)}^{2}+\dots +{a}_{i\ell}{\left(m\right)}^{2}},
(7)
where the Huber loss function is defined as
{H}_{\delta}\left(\theta \right)=\left\{\begin{array}{cc}{\theta}^{2}\hfill & \text{if}\phantom{\rule{0.3em}{0ex}}\left\theta \right\le \delta \hfill \\ 2\delta \left\theta \right{\delta}^{2}\hfill & \text{otherwise}.\hfill \end{array}\right.
(8)
The squared error and Huber loss function are illustrated in Figure 9. It can be seen that for small errors, these two loss functions are exactly the same while for large errors, Huber loss which increases linearly is less than the squared error loss which increases quadratically. Because Huber loss penalizes much less than the squared error loss for large errors, the Huber group LASSO is more robust than group LASSO when there exists large noise or outliers in the data. It is also known that the Huber loss is nearly as efficient as squared error loss for Gaussian errors [24].
For convenience, we define some notations and rewrite the problems (6) and (7) in more compact forms. Let Y_{
i
} = [Y_{
i
}(1)^{T} , . . . , Y_{
i
}(m)^{T}]^{T}, the vector stacking observations of the i th target gene across all datasets, where Y_{
i
}(k)^{T} is the i th row of Y(k). Let b_{
iℓ
}= [a_{
iℓ
}(1), . . . , a_{
iℓ
}(m)]^{T}, the vector containing the grouped parameters. Denote by {b}_{i}={\left[{b}_{i1}^{T},\dots ,{b}_{ip}^{T}\right]}^{T} the vector containing all parameters related to the regulation of the i th target gene. According to the orders of the parameters in b_{
i
}, rearrange the rows of X(k) and piece them together to have X={\left[{X}_{1}^{T},\dots ,{X}_{p}^{T}\right]}^{T} where X_{
i
} = diag(X_{
i
}(1)^{T} , . . . , X_{
i
}(m)^{T}) with X_{
i
}(k)^{T} being the i th row of X(k). Then (7) can be rewritten as
{\displaystyle \sum _{i=1}^{p}}{\displaystyle \sum _{j=1}^{n}}{\omega}_{j}{H}_{\delta}\left({y}_{ij}{x}_{j}^{T}{b}_{i}\right)+\lambda {\displaystyle \sum _{i=1}^{p}}{\displaystyle \sum _{\ell =1}^{n}}{\u2225{b}_{i\ell}\u2225}_{2},
(9)
where x_{
j
} is the j th column of X, y_{
ij
} is the j th element of Y_{
i
}, n={\sum}_{k=1}^{m}{n}_{k} and {\omega}_{i}={w}_{1}I\left(i\le {n}_{1}\right)+{\sum}_{k=2}^{m}{w}_{k}I\left({\sum}_{l=1}^{k=1}{n}_{l}<i\le {\sum}_{l=1}^{k}{n}_{l}\right). (6) can be rewritten similarly.
Optimization algorithm
The minimization of problem (9) is not easy as the penalty term is not differentiable at zero and the Huber loss does not have the second order derivatives at the transition points, ±δ. Observed that fixing i, the problem (9) can be decomposed into p suboptimization problems. For each, we get b_{
i
} by minimizing
J\left(b\right)={\displaystyle \sum _{j=1}^{n}}{\omega}_{j}{H}_{\delta}\left({y}_{j}{x}_{j}^{T}{b}_{i}\right)+\lambda {\displaystyle \sum _{\ell =1}^{n}}{\u2225{b}_{\ell}\u2225}_{2},
(10)
where for notational convenience, we omit the subscript i here and b_{
ℓ
} is a block of parameters of b, i.e. b={\left[{b}_{1}^{T},\dots ,{b}_{p}^{T}\right]}^{T}.
To optimize (10), an iterative method is developed by constructing an auxiliary function, the optimization of which keeps J (b) decreasing. As in [13], given any current estimate b^{(k)}, a function Q(b  b^{(k)}) is an auxiliary function for J (b) if conditions
J\left({b}^{\left(k\right)}\right)=Q\left({b}^{\left(k\right)}{b}^{\left(k\right)}\right)\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}and\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}J\left(b\right)\le Q\left(b{b}^{\left(k\right)}\right)\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\text{for}\phantom{\rule{0.3em}{0ex}}\text{all}\phantom{\rule{0.3em}{0ex}}b,
(11)
are satisfied. In this study, we construct the auxiliary function as
Q\left(b{b}^{\left(k\right)}\right)={\displaystyle \sum _{j=1}^{n}}{\omega}_{j}{H}_{\delta}\left({y}_{j}{x}_{j}^{T}{b}^{\left(k\right)}\right){\displaystyle \sum _{j=1}^{n}}{\omega}_{j}{{H}^{\prime}}_{\delta}\left({y}_{j}{x}_{j}^{T}{b}^{\left(k\right)}\right){x}_{j}^{T}\left(b{b}^{\left(k\right)}\right)+2\gamma {\u2225b{b}^{\left(k\right)}\u2225}_{2}^{2}+\lambda {\displaystyle \sum _{\ell =1}^{p}}{\u2225{b}_{\ell}\u2225}_{2},
(12)
where γ is the largest eigenvalue of {\sum}_{j=1}^{n}{\omega}_{j}{x}_{j}{x}_{j}^{T}. It can be easily shown that this auxiliary function satisfies (11).
Considering the block structure of b, we apply a blockwise descent strategy [14], i.e., cyclically optimize one block of parameters, b_{
j
} , at a time. Denote by {b}^{\left(k\right)}\left(\ell \right)={\left[{b}_{1}^{{\left(k+1\right)}^{T}},\dots ,{b}_{\ell}^{{\left(k+1\right)}^{T}},{b}_{\ell +1}^{{\left(k\right)}^{T}},\dots ,{b}_{p}^{{\left(k\right)}^{T}}\right]}^{T} the vector after updating the R th block. Given b^{(k)}(ℓ− 1), update it to b^{(k)}(ℓ) by computing
\begin{array}{c}{b}_{\ell}^{\left(k+1\right)}=\underset{{b}_{\ell}}{\text{arg}\phantom{\rule{0.3em}{0ex}}\text{min}}\phantom{\rule{0.3em}{0ex}}Q\left({\left[{b}_{1}^{{\left(k+1\right)}^{T}},\dots ,{b}_{\ell 1}^{{\left(k+1\right)}^{T}},{b}_{\ell}^{T},{b}_{\ell +1}^{{\left(k\right)}^{T}},\dots ,{b}_{p}^{{\left(k\right)}^{T}}\right]}^{T}{b}^{\left(k\right)}\left(\ell 1\right)\right)\hfill \\ =\left(\frac{1}{4\gamma}\frac{\lambda}{4\gamma {\u2225{\sum}_{j=1}^{n}{\omega}_{j}{{H}^{\prime}}_{\delta}\left({y}_{j}{x}_{j}^{T}{b}^{\left(k\right)}\left(\ell 1\right)\right){x}_{j,\ell}+4\gamma {b}_{\ell}^{\left(k\right)}\u2225}_{2}}\right)\hfill \\ \times \left(\sum _{j=1}^{n}{\omega}_{j}{{H}^{\prime}}_{\delta}\left({y}_{j}{x}_{j}^{T}{b}^{\left(k\right)}\left(\ell 1\right)\right){x}_{j,\ell}+4\gamma {b}_{\ell}^{\left(k\right)}\right).\hfill \end{array}
(13)
where x_{
j,ℓ
} is the block of elements in x_{
j
} corresponding to b_{
ℓ
} and (·)_{+} = max(·, 0). We repeat to update every block using (13) until it converges. For a specific value of λ, the whole procedure is described as follows:

1
Initialize b(0). Set iteration number k = 0.

2
Cycle through (13) one at a time to update the ℓ th block, ℓ = 1, . . . , p

3
If {b^{(k)}} converges to b^{∗}, go to the next step. Otherwise, set k := k + 1 and go to Step 2.

4
Return the solution b^{∗}.
Note that the algorithm can be adapted to solve (6) with quite similar derivations. In the following section, we show that the sequence {b^{(k)}} generated from the algorithm guarantees the objective function J (b) keep decreasing. We also show that the limit point of the sequence generated is indeed the minimum point of J (b).
Convergence analysis
The convergence of the optimization algorithm for the minimization of (10) is analyzed in the way similar to [25]. We first show the descent property of the algorithm.
Lemma 1 The sequence {b(k)} generated from the optimization algorithm keeps the objective function J(b) decreasing, i.e., J(b^{(k)}) ≥ J (b^{(k+1)}).
Proof By (11) and (13), we have
\begin{array}{cc}J\left({b}^{\left(k\right)}\right)\hfill & =Q\left({b}^{\left(k\right)}{b}^{\left(k\right)}\right)\ge Q\left({b}^{\left(k\right)}\left(1\right){b}^{\left(k\right)}\right)\hfill \\ \ge Q\left({b}^{\left(k\right)}\left(2\right){b}^{\left(k\right)}\left(1\right)\right)\ge \cdots \ge Q\left({b}^{\left(k\right)}\left(p\right){b}^{\left(k\right)}\left(p1\right)\right)\ge J\left({b}^{\left(k\right)}\left(p\right)\right)\hfill \\ =J\left({b}^{\left(k+1\right)}\right).\hfill \end{array}
Next, we show that if the generated sequence satisfies some conditions, it converges to the optimal solution.
Lemma 2 Assume the data (y, X) lies on a compact set and the following conditions are satisfied:
1 The sequence {b^{(k)}} is bounded.
2 For every convergent subsequence \left\{{b}^{\left({n}_{k}\right)}\right\}\subset \left\{{b}^{\left(k\right)}\right\}, the successive differences converge to zeros, {b}^{\left({n}_{k}\right)}{b}^{\left({n}_{k}1\right)}\to 0.
Then, every limit point b^{∞} of the sequence {b^{(k)}} is a minimum for the function J(b), i.e., for any \delta ={\left({\delta}_{1}^{T},\dots ,{\delta}_{p}^{T}\right)}^{T}\in {\mathbb{R}}^{mp},
Proof For any b={\left({b}_{1}^{T},\dots ,{b}_{p}^{T}\right)}^{T}\in {\mathbb{R}}^{mp} and \delta \left(j\right)={\left({0}^{T},\dots ,{\delta}_{j}^{T},\dots ,{0}^{T}\right)}^{T}\in {\mathbb{R}}^{mp}
{\displaystyle \underset{\alpha \downarrow 0+}{\text{lim}\text{inf}}}\left\{\frac{J\left(b+\alpha \delta \left(j\right)\right)J\left(b\right)}{\alpha}\right\}={\nabla}_{j}f{\left(b\right)}^{T}{\delta}_{j}+\underset{\alpha \downarrow 0+}{\text{lim}\text{inf}}\left\{\frac{\lambda {\u2225{b}_{j}+\alpha {\delta}_{j}\u2225}_{2}{\u2225{b}_{j}\u2225}_{2}}{\alpha}\right\},
where f\left(b\right)={\sum}_{i=1}^{n}{\omega}_{i}{H}_{\delta}\left({y}_{i}{x}_{i}^{T}b\right) and ∇_{
j
} represents the partial derivatives with respect to the j th block of parameters. Denote the second term by ∂P (b_{
j
} ; δ_{j} ) and it has
\partial P\left({b}_{j};{\delta}_{j}\right)=\left\{\begin{array}{cc}\lambda \frac{{b}_{j}^{T}{\delta}_{j}}{{\u2225{b}_{j}\u2225}_{2}}\hfill & \text{if}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{b}_{j}\ne 0,\hfill \\ \lambda {\u2225{\delta}_{j}\u2225}_{2}\hfill & \text{otherwise}\text{.}\hfill \end{array}\right.
(14)
We assume the subsequence \left\{{b}^{\left({n}_{k}\right)}\right\} converges to {b}^{\infty}={\left({b}_{1}^{\infty T},\dots ,{b}_{p}^{\infty T}\right)}^{T}\in {\mathbb{R}}^{mp}. From condition 2. and (14), we have
{b}^{\left({n}_{k}1\right)}\left(j\right)={\left({b}_{1}^{{\left({n}_{k}\right)}^{T}},\dots ,{b}_{j}^{{\left({n}_{k}\right)}^{T}},{b}_{j+1}^{{\left({n}_{k}1\right)}^{T}},\dots ,{b}_{p}^{{\left({n}_{k}1\right)}^{T}}\right)}^{T}\to {b}^{\infty},\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\text{as}\phantom{\rule{0.3em}{0ex}}k\to \infty ,
and
\text{if}\phantom{\rule{0.3em}{0ex}}{b}_{j}^{\infty}\ne 0,\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\partial P\left({b}_{j}^{\left({n}_{k}\right)};{\delta}_{j}\right)\to \partial P\left({b}_{j}^{\infty};{\delta}_{j}\right);\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\text{if}\phantom{\rule{0.3em}{0ex}}{b}_{j}^{\infty}=0,\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\partial P\left({b}_{j}^{\infty};{\delta}_{j}\right)\ge \underset{k\to \infty}{\text{lim}\phantom{\rule{0.3em}{0ex}}\text{inf}}\partial P\left({b}_{j}^{\left({n}_{k}\right)};{\delta}_{j}\right),
(15)
since {b}_{j}^{T}{\delta}_{j}\le {\u2225{b}_{j}\u2225}_{2}{\u2225{\delta}_{j}\u2225}_{2}.
As {b}_{j}^{\left({n}_{k}\right)} minimizes Q\left({\left({b}_{1}^{{\left({n}_{k}\right)}^{T}},\dots ,{b}_{j1}^{{\left({n}_{k}\right)}^{T}},{b}_{j}^{T},{b}_{j+1}^{{\left({n}_{k}1\right)}^{T}},\dots ,{b}_{p}^{{\left({n}_{k}1\right)}^{T}}\right)}^{T}{b}^{\left({n}_{k}\right)}\left(j1\right)\right) with respect to the j th block of parameters, using (14), we have
{\nabla}_{j}q{\left({b}^{\left({n}_{k}\right)}\left(j\right){b}^{\left({n}_{k}\right)}\left(j1\right)\right)}^{T}{\delta}_{j}+\partial P\left({b}_{j}^{\left({n}_{k}\right)};{\delta}_{j}\right)\ge 0,\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\text{for}\phantom{\rule{0.3em}{0ex}}\text{all}\phantom{\rule{0.3em}{0ex}}k,
(16)
with
\begin{array}{cc}q\left({b}^{\left({n}_{k}\right)}\left(j\right){b}^{\left({n}_{k}\right)}\left(j1\right)\right)\hfill & ={\displaystyle \sum _{i=1}^{n}}{\omega}_{i}{H}_{\delta}\left({y}_{i}{x}_{i}^{T}{b}^{\left({n}_{k}\right)}\left(j1\right)\right)\hfill \\ \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{\displaystyle \sum _{i=1}^{n}}{\omega}_{i}{{H}^{\prime}}_{\delta}\left({y}_{i}{x}_{i}^{T}{b}^{\left({n}_{k}\right)}\left(j1\right)\right){x}_{i}^{T}\left({b}^{\left({n}_{k}\right)}\left(j\right){b}^{\left({n}_{k}\right)}\left(j1\right)\right)\hfill \\ \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}+2\gamma {\u2225{b}^{\left({n}_{k}\right)}\left(j\right){b}^{\left({n}_{k}\right)}\left(j1\right)\u2225}_{2}^{2}\hfill \end{array}
Due to condition 2.,
{\nabla}_{j}q\left({b}^{\left({n}_{k}\right)}\left(j\right){b}^{\left({n}_{k}\right)}\left(j1\right)\right)\to {\nabla}_{j}f\left({b}^{\infty}\right)\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}as\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}k\to \infty .
(17)
Therefore, (15), (16) and (17) yield
{\nabla}_{j}f{\left({b}^{\infty}\right)}^{T}{\delta}_{j}+\partial P\left({b}^{\infty};{\delta}_{j}\right)\ge {\displaystyle \underset{k\to \infty}{\text{lim}\phantom{\rule{0.3em}{0ex}}\text{inf}}}\left\{{\nabla}_{j}q({b}^{\left({n}_{k}\right)}\left(j\right){\left({b}^{\left({n}_{k}\right)}\left(j1\right)\right)}^{T}{\delta}_{j}+\partial P\left({b}_{j}^{\left({n}_{k}\right)};{\delta}_{j}\right)\right\}\ge 0,
(18)
for any 1 ≤ j ≤ p.
For \delta ={\left({\delta}_{1}^{T},\dots ,{\delta}_{p}^{T}\right)}^{T}\in {\mathbb{R}}^{mp}, due to the differentiability of f (b),
\begin{array}{cc}{\displaystyle \underset{\alpha \downarrow 0+}{\text{lim}\text{inf}}}\left\{\frac{J\left({b}^{\infty}+\alpha \delta \right)J\left({b}^{\infty}\right)}{\alpha}\right\}\hfill & ={\displaystyle \sum _{j=1}^{p}}{\nabla}_{j}f{\left({b}^{\infty}\right)}^{T}{\delta}_{j}+{\displaystyle \sum _{j=1}^{p}}\underset{\alpha \downarrow 0+}{\text{lim}\text{inf}}\left\{\frac{\lambda {\u2225{b}_{j}^{\infty}+\alpha {\delta}_{j}\u2225}_{2}{\u2225{b}_{j}^{\infty}\u2225}_{2}}{\alpha}\right\}\hfill \\ ={\displaystyle \sum _{j=1}^{p}}\left\{{\nabla}_{j}f{\left({b}^{\infty}\right)}^{T}{\delta}_{j}+\partial P\left({b}^{\infty};{\delta}_{j}\right)\right\}\ge 0.\hfill \end{array}
Finally, we show that the sequence generated from the proposed algorithm satisfies these two conditions.
Theorem 3 Assuming the data (y, X) lies on a compact set and no column of X is identically 0, the sequence {b^{(k)}} generated from the algorithm converges to the minimum point of the objective function J (b).
Proof We only need to show that the generated sequence meets the conditions in Lemma 2.
For the sake of notational convenience, for fixed j and {\left({b}_{1}^{T},\dots ,{b}_{j1}^{T},{b}_{j+1}^{T},\dots ,{b}_{p}^{T}\right)}^{T} define
\chi \left(\cdot \right):u\mapsto J\left({\left({b}_{1}^{T},\dots ,{b}_{j1}^{T},{u}^{T},{b}_{j+1}^{T},\dots ,{b}_{p}^{T}\right)}^{T}\right).
Let b(u) be the vector containing u as its j th block of parameters with other blocks being the fixed values.
Assume u + δ and u represent the values of the j th block of parameters before and after the block update, respectively. Hence, as defined in (12), u is obtained by minimizing the following function with respect to the j th block in the algorithm:
\begin{array}{c}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}Q\left(b\left(u\right)b\left(u+\delta \right)\right)\hfill \\ =f\left(b\left(u+\delta \right)\right)+{\nabla}_{{}_{j}}f{\left(b\left(u+\delta \right)\right)}^{T}\left(u\left(u+\delta \right)\right)\hfill \\ \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}+2\gamma {\u2225u\left(u+\delta \right)\u2225}_{2}^{2}+\lambda {\u2225u\u2225}_{2}+\lambda \sum _{\ell \ne j}{\u2225{b}_{\ell}\u2225}_{2}.\hfill \end{array}
(19)
where f\left(b\right)={\sum}_{i=1}^{n}{\omega}_{i}{H}_{\delta}\left({y}_{i}{x}_{j}^{T}b\right) and {\nabla}_{j}f\left(b\right)={\sum}_{i=1}^{n}{\omega}_{i}{{H}^{\prime}}_{\delta}\left({y}_{i}{x}_{j}^{T}b\right){x}_{i,j}. Thus, u should satisfy
{\nabla}_{j}f\left(b\left(u+\delta \right)\right)4\gamma \delta +\lambda s=0,
(20)
where s = u / u_{2} if u ≠ 0; s_{2} ≤ 1 if u = 0. Then, we have
\begin{array}{c}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\chi \left(u+\delta \right)\chi \left(u\right)\hfill \\ =f\left(b\left(u+\delta \right)\right)f\left(b\left(u\right)\right)+\lambda \left(\left\right\left(u+\delta \right){}_{2}\leftu\right{}_{2}\right)\hfill \\ ={\nabla}_{j}f{\left(b\left(u+\delta \right)\right)}^{T}\delta {\nabla}_{j}f{\left(b\left(u+\delta \right)\right)}^{T}\delta \left[{\nabla}_{j}f{\left(b\left(u+\delta \right)\right)}^{T}\delta 4\gamma {\delta}^{T}\delta +\lambda {s}^{T}\delta \right]\hfill \\ \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}+4\gamma {\delta}^{T}\delta +\lambda \left(\left\right\left(u+\delta \right){}_{2}\leftu\right{}_{2}{s}^{T}\delta \right)\hfill \\ =\left[{\displaystyle \sum _{i=1}^{n}}{\omega}_{i}{{H}^{\prime}}_{\delta}\left({y}_{i}{x}_{j}^{T}b\left(u+\tau \delta \right)\right){x}_{i,j}^{T}\delta {\displaystyle \sum _{i=1}^{n}}{\omega}_{i}{{H}^{\prime}}_{\delta}\left({y}_{i}{x}_{j}^{T}b\left(u+\delta \right)\right){x}_{i,j}^{T}\delta \right]\hfill \\ \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}+4\gamma {\delta}^{T}\delta +\lambda \left(\left\right\left(u+\delta \right){}_{2}\leftu\right{}_{2}{s}^{T}\delta \right)\hfill \\ \ge 2\gamma \left(1\tau \right)\left\right\delta {}_{2}^{2}+4\gamma \left\delta \right{}_{2}^{2}\ge 2\gamma \left\right\delta {}_{2}^{2}.\hfill \end{array}
(21)
The second and third equalities are obtained using mean value theorem with τ ∈ (0, 1) and (20). For the first inequality, the following property of the Huber loss function and the property of subgradient are used.
\left({H}_{\delta}^{\prime}\left({\theta}_{1}\right){H}_{\delta}^{\prime}\left({\theta}_{2}\right)\right)\left({\theta}_{1}{\theta}_{2}\right)\le 2{\left({\theta}_{1}{\theta}_{2}\right)}^{2}.
The result from (21) gives that
J\left({b}^{\left(k\right)}\left(j1\right)\right)J\left({b}^{\left(k\right)}\left(j\right)\right)\ge 2\gamma {\u2225{b}_{j}^{\left(k\right)}{b}_{j}^{\left(k+1\right)}\u2225}_{2}^{2}=2\gamma {\u2225{b}^{\left(k\right)}\left(j1\right){b}^{\left(k\right)}\left(j\right)\u2225}_{2}^{2},
(22)
where {b}^{\left(k\right)}\left(j\right)={\left[{b}_{1}^{{\left(k+1\right)}^{T}},\dots ,{b}_{j}^{{\left(k+1\right)}^{T}},{b}_{j+1}^{{\left(k\right)}^{T}},\dots ,{b}_{p}^{{\left(k\right)}^{T}}\right]}^{T}.
Using (22) repeatedly across every block, for any k, we have
J\left({b}^{\left(k\right)}\right)J\left({b}^{\left(k+1\right)}\right)\ge 2\gamma {\u2225{b}^{\left(k\right)}{b}^{\left(k+1\right)}\u2225}_{2}^{2}.
Note that by Lemma 1, {J (b^{(k)})} converges as it keeps decreasing and is bounded from below. The convergence of {J (b^{(k)})} yield the convergence of {b^{(k)}}. Hence, conditions of Lemma 2 hold which imply that the limit of {b^{(k)}} is the minimum point of J (b).
Implementation
The tuning parameter λ controls the sparseness of the resulted network. A network solution path can be obtained by computing networks on a grid of λ values from {\lambda}_{\text{max}}=\underset{i,\ell}{\text{max}}\u2225{\sum}_{j=1}^{n}{\omega}_{j}{{H}^{\prime}}_{\delta}\left({y}_{ij}\right){x}_{j,\ell}\u2225, which is the smallest value that gives the empty network, to a small value, e.g. λ_{min} = 0.01λ_{max}. In our previous work [12], BIC criterion is used to pick a specific λ value which corresponds to a determinant network topology. A method called "stability selection" recently proposed by Meinshausen and Buhlmann [15] finds a network with probabilities for edges. Stability selection performs the network inference method, e.g. group LASSO, many times, resampling the data in each run and computing the frequencies with which each edge is selected across these runs. It has been used with the linear regression method to infer GRNs from steadystate gene expression data in Haury et al. [26] and has shown perspective effectiveness. In this study, we adapt the stability selection method to finding GRN topology from multiple timecourse gene expression datasets. Given a family of m timecourse gene expression datasets \stackrel{\u0303}{X}\left(k\right)\in {\mathcal{R}}^{p\times {n}_{k}}, k = 1, . . . , m, for a specific λ ∈ Λ, the stability selection procedure is as follows

1
Use moving block bootstrap to draw N bootstrap samples for every dataset to form N bootstrap families of multiple timecourse datasets, i.e. \stackrel{\u0303}{X}{\left(k\right)}^{*\left(b\right)}{\}}_{k=1}^{m}, b = 1, . . . , N

2
Use the proposed Huber group LASSO to infer {\left\{A{\left(k\right)}_{\lambda}^{*\left(b\right)}\right\}}_{k=1}^{m} from the b th bootstrap family of datasets. Denote by {A}_{\lambda}^{*\left(b\right)} the network topology shared by {\left\{A{\left(k\right)}_{\lambda}^{*\left(b\right)}\right\}}_{k=1}^{m}.

3
Compute the frequencies for each edge (i, j), i.e., from the gene j to gene i, in the network
\prod _{\lambda}\left(i,j\right)=\frac{\#\left\{b:{A}_{\lambda}^{*\left(b\right)}\left(i,j\right)\ne 0\right\}}{N},
(24)
where {A}_{\lambda}^{*\left(b\right)}\left(i,j\right) is the (i, j)'s entry of {A}_{\lambda}^{*\left(b\right)} and #{·} is the number of elements in that set.
For a set of λ ∈ Λ, the probability of each edge in the inferred network is
\prod \left(i,j\right)=\underset{\lambda \in \Lambda}{\text{max}}\prod _{\lambda}\left(i,j\right)
(25)
The final network topology can be obtained by setting a threshold, edges with probabilities or scores less than which are considered nonexistent. This study only focus on giving a list of edges with scores. The selection of threshold is not discussed here. The stability selection procedure can also be applied with the group LASSO method (6).
Since the data used are time series data, the moving block bootstrap method is employed in the first step to draw bootstrap samples from each dataset. For a dataset with n observations, in the moving block bootstrap with block length l, the data is split into n − l + 1 blocks: block j consists of observations j to j + l − 1, j = 1, . . . , n − l + 1. [n/b] blocks are randomly drawn from n − l + 1 blocks with replacement and are aligned in the order they are picked to form a bootstrap sample.
Another tuning parameter δ controls the degree of robustness. Generally, it picks \delta =1.345\widehat{\sigma} where \widehat{\sigma} is the estimated standard deviation of the error and \widehat{\sigma}=MAD/0.6745, where MAD is the median absolute deviation of the residuals. In this study, we use the least absolute deviations (LAD) regression to obtain the residuals. To avoid the overfitting of LAD which leads to a very small δ, we choose by \delta =\text{max}\left(1.345\widehat{\sigma},1\right).