Signal Lasso with Non-Convex Penalties for Efficient Network Reconstruction and Topology Inference

TL;DR

This paper introduces a novel non-convex penalty-based signal lasso method for efficient network reconstruction and topology inference, demonstrating improved speed and accuracy over existing methods through extensive simulations and real data application.

Contribution

The paper proposes a new signal lasso approach using non-convex penalties, enabling better shrinkage of signal parameters and faster computation for large networks.

Findings

Effective in estimating signal parameters in linear regression.

Faster and more accurate than existing convex penalty methods.

Validated on simulations and real behavioral data.

Abstract

Inferring network structures remains an interesting question for its importance on the understanding and controlling collective dynamics of complex systems. The existing shrinking methods such as Lasso-type estimation can not suitably reveal such property. A new method recently suggested, called by {\it signal lasso} (or its updating version: adaptive signal lasso) was proposed to solve the network reconstruction problem, where the signal parameter can be shrunk to either 0 or 1 in two different directions. The signal lasso or adaptive signal lasso employed the additive penalty of signal and non-signal terms which is a convex function and easily to complementation in computation. However their methods need tuning the one or two parameters to find an optimal solution, which is time cost for large size network. In this paper we propose new signal lasso method based on two penalty functions to estimate the signal parameter and uncovering network topology in complex network with a small amount of observations. The penalty functions we introduced are non-convex function, thus coordinate descent algorithms are suggested. We find in this method the tuning parameter can be set to a large enough values such that the signal parameter can be completely shrunk either 0 or 1. The extensive simulations are conducted in linear regression models with different assumptions, the evolutionary-game-based dynamic model and Kuramoto model of synchronization problem. The advantage and disadvantage of each method are fully discussed in various conditions. Finally a real example comes from behavioral experiment is used for illustration. Our results show that signal lasso with non-convex penalties is effective and fast in estimating signal parameters in linear regression model.

Figures

Tables

Table 1. Table 1: Measures for accuracy of network reconstruction

	Actual class	Predicted class
		Signal	Non signal	Unclassified
	signal class	True positive (TP)	False negative (FN)	Unclassified positive (UCP)
	Non signal class	False positive (FP)	True negative (TN)	Unclassified negative (UCN)

Table 2. Table 2: Simulation results in Case I for nine methods: Lasso, Adaptive Lasso, SCAD, MCP, Elastic Net, Signal Lasso, Adaptive Signal Lasso, Signal lasso with product penalty and minimum penalty in linear regression models. All of the results are averaged over 500 independent realizations, where n is the sample size, p p is the number of explanatory variables, p 1 p_{1} is the number of signals (number of β = 1 \beta=1 ). The first panel is for case of p < n p<n and signal is sparse with p 1 = 6 p_{1}=6 . The second panel is for p > n p>n and sparse signal. The third panel consider the cases of non-sparse signal with p 1 = 20 p_{1}=20 . The bottom panel consider the cases of dense signal with p 1 = 100 p_{1}=100 . The noise is introduced by σ = 0.4 , 1 , 2 \sigma=0.4,1,2 , respectively, in three columns.

Method	MSE/UCR/MCC/MCCa
$(n,p,p_1,\sigma )$	$(100,30,6,0.4)$	$(100,30,6,1)$	$(100,30,6,2)$
Lasso	0.0014/0.027/1.000/0.919	0.0083/0.221/0.980/0.347	0.0320/0.359/0.863/0.022
A-lasso	0.0015/0.343/1.000/0.900	0.0095/0.232/0.985/0.316	0.0459/0.349/0.757/-0.006
SCAD	0.0008/0.023/1.000/0.937	0.0043/0.108/1.000/0.632	0.0509/0.226/0.702/0.150
MCP	0.0008/0.019/1.000/0.943	0.0046/0.111/1.000/0.633	0.0506/0.220/0.736/0.184
ElasticNet	0.0012/0.333/1.000/0.892	0.0066/0.139/0.970/0.503	0.0239/0.201/0.721/0.225
SigL	0.0007/0.012/1.000/0.966	0.0046/0.133/1.000/0.609	0.0232/0.256/0.943/0.249
ASigL	0.0001/0.002/1.000/0.993	0.0006/0.012/1.000/0.965	0.0086/0.016/0.982/0.933
SL_prod	0.0000/0.000/0.999/1.000	0.0000/0.000/0.999/1.000	0.0058/0.000/0.982/0.982
SL_min	0.0000/0.000/0.999/1.000	0.0000/0.000/0.999/1.000	0.0058/0.000/0.982/0.982
$(n,p,p_1,\sigma )$	$(50,150,6,0.4)$	$(50,150,6,1)$	$(150,150,6,2)$
Lasso	0.0016/0.043/0.990/0.529	0.0085/0.162/0.970/0.107	0.034/0.268/0.559/-0.042
A-lasso	0.0026/0.076/0.980/0.375	0.0130/0.191/0.860/0.061	0.0415/0.264/0.577/0.035
SCAD	0.0002/0.007/1.000/0.896	0.0111/0.037/0.609/0.305	0.0315/0.043/0.141/0.038
MCP	0.0002/0.007/1.000/0.892	0.0130/0.032/0.588/0.309	0.0365/0.031/0.144/0.065
ElasticNet	0.0007/0.019/1.000/0.712	0.0037/0.039/0.850/0.359	0.0126/0.050/0.439/0.144
SigL	0.0004/0.013/1.000/0.839	0.0028/0.051/1.000/0.482	0.0127/0.085/0.883/0.255
ASigL	0.0009/0.001/0.989/0.976	0.0213/0.005/0.821/0.781	0.0111/0.003/0.873/0.841
SL_prod	0.0000/0.000/0.999/1.000	0.0030/0.000/0.967/0.967	0.0584/0.000/0.559/0.559
SL_min	0.0000/0.000/0.999/1.000	0.0030/0.000/0.967/0.967	0.0584/0.000/0.559/0.559
$(n,p,p_1,\sigma )$	$(50,150,20,0.4)$	$(50,150,20,1)$	$(50,150,20,2)$
Lasso	0.0279/0.202/0.905/0.091	0.0431/0.255/0.847/-0.01	0.0796/0.311/0.600/-0.108
A-lasso	0.0104/0.140/0.975/0.372	0.0305/0.237/0.918/0.057	0.0802/0.285/0.646/-0.066
SCAD	0.1556/0.085/0.095/-0.01	0.1565/0.082/0.124/0.057	0.1680/0.077/0.107 /0.004
MCP	0.1709/0.072/0.078/-0.01	0.1771/0.069/0.091/-0.00	0.1919/0.058/0.064/-0.008
ElasticNet	0.0142/0.119/0.961/0.362	0.0276/0.153/0.883/0.172	0.0533/0.157/0.609/0.068
SigL	0.0015/0.044/1.000/0.815	0.0092/0.118/0.997/0.549	0.0349/0.161/0.927/0.370
ASigL	0.0368/0.005/0.840/0.819	0.0400/0.004/0.827/0.809	0.0532/0.004/0.767/0.751
SL_prod	0.0028/0.000/0.988/0.988	0.0168/0.000/0.931/0.931	0.0835/0.000/ 0.689/0.689
SL_min	0.0028/0.000/0.988/0.988	0.0168/0.000/0.931/0.931	0.0835/0.000/0.689/0.689
$(n,p,p_1,\sigma )$	$(50,150,100,0.4)$	$(50,150,100,1)$	$(50,150,100,2)$
Lasso	0.7554/0.408/0.120/-0.382	0.7435/0.408/0.121/-0.376	0.7622/0.411/0.135/-0.376
A-lasso	0.7498/0.364/0.122/-0.268	0.7388/0.361/0.135/-0.257	0.7671/0.361/0.146/-0.257
SCAD	0.9004/0.082/0.045/-0.058	0.8662/0.076/0.042/0.050	0.8545/0.080/0.039 /-0.053
MCP	1.0681/0.046/0.015/-0.048	1.0106/0.045/0.007/-0.044	0.9901/0.037/0.010/-0.045
ElasticNet	0.5375/0.311/0.016/-0.185	0.5337/0.313/0.153/0.188	0.5314/0.321/0.177/-0.175
SigL	0.0208/0.145/0.981/0.657	0.0414/0.202/0.966/0.527	0.0837/0.213/0.886/0.448
ASigL	0.2982/0.008/0.437/0.427	0.3080/0.007/0.425/0.416	0.3049/0.006/0.432/0.424
SL_prod	0.1985/0.000/0.597/0.597	0.2082/0.000/0.582/0.582	0.1968/0.000/ 0.601/0.601
SL_min	0.1985/0.000/0.597/0.597	0.2082/0.000/0.582/0.582	0.1968/0.000/0.601/0.601

Table 3. Table 3: Payoff matrices used in the experiments. The payoff matrix of prisoner’s dilemma game with punishment option (treatment I) and the standard prisoner’s dilemma (treatment II and treatment III) are shown in left panel and right panel, respectively.

	C	D	P
C	2	-2	-5
D	4	0	-3
P	2	-2	-5

Table 4. Table 4: Simulation results in linear regression model with design matrix X X have 2 zero columns for nine methods

Method	MSE/UCR/MCC/MCCa
$(n,p,p_1,\sigma )$	$(50,150,6,0.4)$	$(50,150,6,1)$	$(50,150,6,2)$
Lasso	0.0020/0.049/0.960/0.435	0.011/0.173/0.716/0.036	0.0423/0.265/0.208/-0.090
A-lasso	0.0024/0.071/0.960/0.346	0.0146/0.193/0.664/-0.001	0.0464/0.269/0.291/-0.081
SCAD	0.0002/0.004/1.000/0.934	0.0059/0.036/0.799/0.394	0.0312/0.036/0.067/-0.003
MCP	0.0008/0.005/0.980/0.911	0.0065/0.030/0.811/0.426	0.0332/0.025/0.110/0.031
ElasticNet	0.0019/0.034/0.780/0.365	0.0113/0.062/0.316/0.061	0.0316/0.038/0.048/-0.002
SigL	0.0009/0.026/1.000/0.640	0.0069/0.098/0.950/0.188	0.0281/0.128/0.475/0.021
ASigL	0.0103/0.005/0.862/0.805	0.0156/0.004/0.779/0.736	0.0314/0.006/0.525/0.479
SLprod	0.0000/0.000/0.999/1.000	0.0058/0.000/0.936/0.936	0.0891/0.000/0.386/0.386
SLmin	0.0133/0.000/0.859/0.859	0.0191/0.000/0.806/0.806	0.1024/0.000/0.355/0.355

Table 5. Table 5: Comparison of computational complexity measures (CPU) for four methods

Linear Dynamic ($\sigma =2$, Correlated-Gaussian)
$(n,p,p_1)$	(150,30,6)	(150,150,6)	(150,150,60)	(200,30,6)	(200,150,6)	(200,150,60)
SigL	0.389935	16.330999	7.049003	0.466000	11.325395	6.643295
ASigL	0.130004	0.457998	0.503997	0.157996	0.556997	0.644602
SL_prod	0.000000	0.000000	0.000999	0.000000	0.001004	0.000995
SL_min	0.001000	0.001003	0.001001	0.000000	0.001000	0.001004
PDG Dynamic with $N=100,\Delta =0.5,b=1.02$
	k=6			k=60
	BA	ER	SW	BA	ER	SW
SigL	1.138426	0.848676	0.569877	2.514897	3.014704	2.714757
ASigL	0.044129	0.039120	0.039770	0.041260	0.039300	0.039610
SL_prod	0.000441	0.000502	0.000363	0.000374	0.000452	0.000372
SL_min	0.000438	0.000503	0.000361	0.000444	0.000371	0.000453
Kuramoto Dynamic with $N=100,\Delta =0.5,\lambda =15$
	k=6			k=60
	BA	ER	SW	BA	ER	SW
SigL	2.433161	4.215170	4.553449	3.071393	4.965631	3.801474
ASigL	0.118870	0.260713	0.306684	0.202545	0.304234	0.191870
SL_prod	0.000580	0.001184	0.001225	0.001444	0.003402	0.001280
SL_min	0.000423	0.000365	0.000367	0.000375	0.000380	0.000378

Equations

$${\dot{\mathbf{y}}}_i(t)={\psi }_i({\mathbf{y}}_i(t),{\nu }_i)+\sum _{i=1}^N{a}_{ij}{\varphi }_{ij}({\mathbf{y}}_i(t),{\mathbf{y}}_j(t))+{ϵ}_i(t),$$

$$\mathbf{M_{PDG}}=\left(\begin{array}[]{cc}R\ \ S\\ T\ \ P\end{array}\right),$$

$$F_i=\sum _{j\in {\Omega }_i}{s}_i{M}_{PDG}{s}_j=\sum _{j=1,j\ne i}^N{a}_{ij}{P}_{ij},$$

$$\mathbf{Y}=\mathbf{X}\beta +ϵ$$

$$\ell (\beta ,{\lambda }_1,{\lambda }_2|Y,X)=\frac{1}{2}||Y-X\beta |{|}_2^2+{\lambda }_1\sum _{j=1}^p{\omega }_{1j}|{\beta }_j|+{\lambda }_2\sum _{j=1}^p{\omega }_{2j}|{\beta }_j-1|,$$

$$PF(\beta )={\lambda }_1\sum _{j=1}^p{\omega }_{1j}|{\beta }_j|+{\lambda }_2\sum _{j=1}^p{\omega }_{2j}|{\beta }_j-1|$$

$$\tilde{\beta}_{k}\to\left\{\begin{array}[]{l}1,\ \ \ \ \ \ \ \hat{\beta}_{k0}>\alpha,\\ 0,\ \ \ \ \ \ \ \hat{\beta}_{k0}<\alpha,\end{array}\right.$$

$$L^{\ast }(\lambda ,\beta )=\frac{1}{2}\sum _{i=1}^n(y_i-\sum _j{x}_{ij}{\beta }_j{)}^2+P{F}^{\ast }(\lambda ,\beta ).$$

$$P{F}^{\ast }(\lambda ,\beta )=\lambda \sum _{i=1}^p|{\beta }_i({\beta }_i-1)|$$

$$P{F}^{\ast }(\lambda ,\beta )=\lambda \sum _{i=1}^{p}min\{|{\beta }_i|,|{\beta }_i-1|\}$$

$${\beta }_k^{(l+1)}=\arg \underset{{\beta }_k}{\min }\frac{1}{2}\Vert \mathbf{Y}-\sum _{j<k}{\mathbf{x}}_j{\beta }_j^{(l+1)}-\sum _{j>k}{\mathbf{x}}_j{\beta }_j^{(l)}-{\mathbf{x}}_k{\beta }_k{\Vert }_2^2+\lambda |{\beta }_k||{\beta }_k-1|,$$

$$\beta_{k}^{(l+1)}=\left\{\begin{array}[]{ll}\min\left(0,\frac{{\bf x}_{k}^{\prime}{\bf\epsilon}_{k}^{(l)}+\lambda/2}{{\bf x}_{k}^{\prime}{\bf x}_{k}+\lambda}\right)&\text{if}\quad r_{k}^{(l)}\leq 0\\ \max\left(0,\frac{{\bf x}_{k}^{\prime}{\bf\epsilon}_{k}^{(l)}-\lambda/2}{{\bf x}_{k}^{\prime}{\bf x}_{k}-\lambda}\right)\cdot\text{I}({\bf x}_{k}^{\prime}{\bf x}_{k}>\lambda)&\text{if}\quad 0<r_{k}^{(l)}\leq 1/2\\ \min\left(1,\frac{{\bf x}_{k}^{\prime}{\bf\epsilon}_{k}^{(l)}-\lambda/2}{{\bf x}_{k}^{\prime}{\bf x}_{k}-\lambda}\right)\cdot\text{I}({\bf x}_{k}^{\prime}{\bf x}_{k}>\lambda)+\text{I}({\bf x}_{k}^{\prime}{\bf x}_{k}\leq\lambda)&\text{if}\quad 1/2<r_{k}^{(l)}\leq 1\\ \max\left(1,\frac{{\bf x}_{k}^{\prime}{\bf\epsilon}_{k}^{(l)}+\lambda/2}{{\bf x}_{k}^{\prime}{\bf x}_{k}+\lambda}\right)&\text{if}\quad r_{k}^{(l)}>1\end{array}\right.,$$

$${\beta }_k^{(l+1)}=\arg \underset{{\beta }_k}{\min }\frac{1}{2}\Vert \mathbf{Y}-\sum _{j<k}{\mathbf{x}}_j{\beta }_j^{(l+1)}-\sum _{j>k}{\mathbf{x}}_j{\beta }_j^{(l)}-{\mathbf{x}}_k{\beta }_k{\Vert }_2^2+\lambda \min \left(|{\beta }_k|,|{\beta }_k-1|\right),$$

$$\beta_{k}^{(l+1)}=\left\{\begin{array}[]{ll}\text{sign}(r_{k}^{(l)})\cdot\max\left(0,|r_{k}^{(l)}|-\frac{\lambda}{{\bf x}_{k}^{\prime}{\bf x}_{k}}\right)&\text{if}\quad|r_{k}^{(l)}|\leq|r_{k}^{(l)}-1|\\ 1+\text{sign}(r_{k}^{(l)}-1)\cdot\max\left(0,|r_{k}^{(l)}-1|-\frac{\lambda}{{\bf x}_{k}^{\prime}{\bf x}_{k}}\right)&\text{if}\quad|r_{k}^{(l)}|>|r_{k}^{(l)}-1|\end{array}\right.,$$

$$\frac{{\mathbf{x}}_k^{\prime }{ϵ}_k^{(l)}+\lambda /2}{{\mathbf{x}}_k^{\prime }{\mathbf{x}}_k+\lambda }\to 1/2,\frac{{\mathbf{x}}_k^{\prime }{ϵ}_k^{(l)}-\lambda /2}{{\mathbf{x}}_k^{\prime }{\mathbf{x}}_k-\lambda }\to 1/2,\text{I}({\mathbf{x}}_k^{\prime }{\mathbf{x}}_k>\lambda )\to 0,$$

$$TPR=\frac{TP}{TP+FN},TNR=\frac{TN}{TN+FP},$$

$$MCCa=\frac{TP\times TN-FPa\times FNa}{\sqrt{(TP+FPa)(TP+FNa)(TN+FPa)(TN+FNa)}}$$

$$Y=1_{\mathbf{n}}{\beta }_0+{\mathbf{X}}_1{\beta }_1+{\mathbf{X}}_2{\beta }_2+ϵ$$

$$Y_i={\mathbf{X}}_i{\tilde{\beta }}_i+{ϵ}_i,$$

$$\left(\begin{array}[]{cccccc}P_{i1}(t_{1})&\cdots&P_{i,i-1}(t_{1})&P_{i,i+1}(t_{1})&\cdots&P_{i_{N}}(t_{1})\\ P_{i1}(t_{2})&\cdots&P_{i,i-1}(t_{2})&P_{i,i+1}(t_{2})&\cdots&P_{i_{N}}(t_{2})\\ \vdots&\vdots&\vdots&\vdots\\ P_{i1}(t_{L})&\cdots&P_{i,i-1}(t_{L})&P_{i,i+1}(t_{L})&\cdots&P_{i_{N}}(t_{L})\end{array}\right),$$

$$\frac{d{\theta }_i}{dt}={\omega }_i+c\sum _{j=1}^N{a}_{ij}sin({\theta }_j-{\theta }_i),$$

$$Y_i={\omega }_i{1}_L+{\mathbf{X}}_i{\tilde{\beta }}_i+{ϵ}_i,$$

Full text

**Keywords: ** Linear regression, signal parameters, signal lasso, additive penalty, non-convex penalty, network reconstruction

††preprint: APS/123-QED

I  Introduction

Complex network has wide application in many fields and has been made great progress recently (Barabási (2012); Boccaletti et al. (2006)). In a complex network, the pattern of the node-to-node interaction or the network topology is unknown, and uncovering the network topology based on a series of observable quantities obtained from experiments or observations is very important and plays a potential role for understanding and controlling collective dynamics of complex systems (Han et al. (2015); Wang et al. (2011); Peixoto (2018)). Network reconstruction is an inverse problem in network science and has received great attention recently. Some typical examples include the reconstruction of gene networks using expression data in biology (Gardner et al. (2003); Geier, Timmer, and Fleck (2007)), extraction of various functional networks in the human brain from activation data in neuroscience (Grün, Diesmann, and Aertsen (2002); Gútig, Aertsen, and Rotter (2002); Supekar et al. (2008)), and detection of organizational networks in social science and trade networks in economics. Recently evolutionary-game based dynamics has also been used to study the network reconstruction, in which it is possible to observe a series of discrete but a small amount of quantities (Wang et al. (2011); Han et al. (2015); Shi et al. (2020, 2021)). In this case, the problem can be transformed into a statistical linear model with sparse and high-dimensional properties.

To understand how such signal parameter appeared in real practice, we use two typical examples to illustrate the problem. The first is a dynamic equation governed the evolution state in general complex systems, which can be written as following differential equations (Wang, Lai, and Grebogi (2016); Raimondo and De Domenico (2021))

[TABLE]

where ${\bf y}_{i}(t)$ denote a $m$-dimensional internal stats variable of a system consisting of $N$ dynamic units at time $t$, the function $\psi_{i}\in R^{m}$ and $\phi_{ij}$ $\in R^{m}$ respectively define the intrinsic and interaction dynamics of the units. $\epsilon_{i}(t)$ is a dynamic noise term and $\nu_{i}$ is a set of dynamic parameters. Finally, the term $a_{ij}$ defines the interaction topology and is called by adjacency matrix such that $a_{ij}$ = 1 if there is a direct physical interaction from unit $j$ to $i$ and $a_{ij}$ = 0 otherwise. The matrix $A=[a_{ij}]$ completely defines a network with size $N$, that is, an abstraction used to model a system that contains discrete interconnected elements. The elements are represented by nodes (also called vertices) and connections are represented by edges. In general ${\bf y}_{i}(t)$ can be observed as time series data but $a_{ij}$ for $i=1,\cdots,N$ are unknown and need to be estimated. It is clear Eq. (1) can be rewritten as linear regression model if functional form of $\phi_{ij}$ and $y_{i}(t)$ are known. This model include well known synchronization model, oscillator networks, spreading network, and so on (see Wang, Lai, and Grebogi (2016), Shi et al. (2021, 2023)).

Another example is the network reconstruction based on evolutionary game on structured populations, where node represent a player, link means the corresponding two players have game relationship. In this model, the payoff and strategy of each player taken from interaction with other players can be observed but the relationship between players are unknown and need to be estimated (Wang et al. (2011); Han et al. (2015); Peixoto (2018); Shi et al. (2021)). We use a prisoner’s dilemma game (PDG) as an example with the payoff matrix defined by

[TABLE]

This game is characterized with the temptation to defect $T$, reward for mutual cooperation $R$, punishment for mutual defection $P$, as well as the sucker’s payoff $S$. These quantities satisfies that $T>R>P>S$, and the mutual defection are the equilibrium solutions of PDG (Santos and Pacheco (2005); Nowak and May (1992); Hauert and Doebeli (2004); Szabó and Fath (2007)). In some experiment driving studies, each player can interact with other players by choosing either a cooperator (C) or defector (D) to obtain their payoff and procedure is continued until to some predetermined number (Li et al. (2018); Wang et al. (2018, 2017)). Let $s_{i}=(1,0)$ indicates that player $i$ takes action C ( or $s_{i}=(0,1)$ means taking action D). In spatial PDG game, the $i$ player, say focal player, acquires its fitness (total payoff) $F_{i}$ by playing the game with all its connected neighbors, which is defined by:

[TABLE]

where $\Omega_{i}$ represents the set of all connected neighbors of player $i$, $P_{ij}=s_{i}M_{PDG}s_{j}$. Eq. (3) can be converted into a linear model, where $a_{ij}$ is the elements of adjacency matrix for a given network. $a_{ij}=1$ means player $i$ and $j$ are connected, otherwise ($a_{ij}=0$) they are not connected. When the process continue forward to produce a time series data, we obtained a linear model with signal parameters (0 or 1) as regression coefficients. The detailed description of this framework will be presented by a real example of human behavioral experiments in Section VI of this paper.

In the fields of complex systems and applied physics as illustrated by above two examples, traditional estimation methods such as compressed sensing (CS) and the lasso method have been widely used for network reconstruction (Wang et al. (2011); Han et al. (2015); Shi et al. (2020)). These studies found that the lasso method is particularly robust against noise in the data, making it effective for reconstructing sparse signals. Given that the true value of $\hat{a}_{ij}$ should be either 0 or 1, they determined that player $i$ has a game relationship (connection) with player $j$ if $|\hat{a}_{ij}-1|\leq 0.1$, and no relationship if $|\hat{a}_{ij}|\leq 0.1$. Here, $\hat{a}_{ij}$ is the estimator of $a_{ij}$. In cases where the estimator falls outside these ranges, the relationship between the players is not identifiable. Although CS and the lasso method are effective at shrinking parameter estimates toward zero, especially in sparse complex networks, they often fail to shrink existing links to their true value of 1, which can reduce estimation accuracy. To address this limitation, Shi et al. (2021) introduced the signal lasso (SigL) method, which adds a penalty for the signal parameter. They found that SigL outperformed both lasso and CS in network reconstruction. Malizia et al. (2024) applied the SigL method to identify high-order interactions in coupled dynamical systems. As an advancement, Shi et al. (2023) further developed the adaptive signal lasso (ASigL) method, which has the desirable property of shrinking parameters fully to either 0 or 1, thus improving accuracy in network reconstruction.

In this paper we propose a new version of signal lasso based on two kinds of penalty function to estimate the signal parameter and uncovering network topology in complex network with a small amount of observations. We find the tuning parameter can be set to a large enough values such that the signal parameter can be completely shrunk either 0 or 1, therefore there is no need to tuning parameter. The theoretical properties of this method has been studied, and some simulations results are conducted to make fully comparison for all methods, which including in linear regression models with different assumption, evolutionary game dynamic model and synchronization model. Finally three dataset from the human experiment are deeply analyzed. All simulation results in terms of our proposed method is effective to uncover the signals presented in the models and greatly decreasing computational complexity of the method.

II  Signal Lasso with additive penalty

Consider the general linear regression model

[TABLE]

where $\epsilon$ is a noise or random error with mean zero and finite variance, ${\bf X}=[x_{ij}]$ is an $n\times p$ matrix, ${\bf Y}=[y_{i}]$ is an $n\times 1$ vector, and $\beta=[\beta_{i}]$ is a $p\times 1$ unknown vector. $\epsilon$ is a random error. To eliminate the intercept from (4), throughout this paper, we centre the response variable and predictor variable so that the mean of response is zero. We assume the parameter $\beta$ here have a signal property, e.g. the true values of $\beta_{j}$ for $j=1,\cdots,p$ either 0 or 1. This kind of problem is very common in reconstruction of complex network to identifying the signal of connection or not connection (Han et al. (2015); Wang et al. (2011)).

In many practical situations, the signal always sparse, therefore lasso method, as a useful shrinkage estimation, have been used to find the estimator of parameter $\beta$ (Han et al. (2015); Wang et al. (2011)). Although lasso method, have ability of shrink the parameter estimates toward to zero under the natural sparsity in complex network, the existent links between nodes can not be shrunk to its true value of 1, which will inevitably decrease the estimation accuracy in most cases. For this reason, Shi et al. (2021) proposed signal lasso (SigL) method to solve the network reconstruction problem and they found SigL performed better than lasso and CS methods. However signal lasso still have some drawback in which the elements of $\hat{\beta}$ that fall in interval (0.1,0.9) can not be successfully selected to the correct class. Furthermore Shi et al. (2023) proposed adaptive signal lasso, where appropriate weights are imposed to the penalty to promoting the estimation accuracy. The objective function with weighted penalty $\ell$ is defined by

[TABLE]

where $\omega_{1j}$ and $\omega_{2j}$ are two weights. The estimator defined by $\hat{\beta}^{*}=\arg\min_{\beta}\ell(\beta,\lambda_{1},\lambda_{2}|Y,X)$ is called by adaptive signal Lasso(ASigL). The penalty function here is defined as

[TABLE]

They employed $\omega_{1k}=1$ and $\omega_{2k}=|\hat{\beta}_{k0}|$ for $k=1,\cdots,p$ in their study, where $\hat{\beta}_{k0}$ is an initial estimator of $\beta_{k}$, for example $\hat{\beta}_{k0}$ can be ordinary least square estimator for $p<n$ or ridge estimator for $p>n$.

If Eq. (5) is re-parametrized by $\lambda=\lambda_{2}$ and $\alpha=\lambda_{1}/\lambda_{2}$, they found following fact when $\lambda\to+\infty$ for fixed $\alpha$

[TABLE]

since $\alpha_{2}\to\alpha$ and $\alpha_{1}\to\alpha$ when $\lambda\to+\infty$. This result indicates that if $\lambda$ is selected large enough, the estimators from adaptive signal Lasso can be completely shrunk to either 0 or 1 and thus remove the unidentified set which will be presented in signal lasso method Shi et al. (2021). Signal lasso (SigL) corresponds the result using penalty function Eq. 5 with two weight functions equal to 1, in which the property presented in Eq. 6 is not true. The choices of tuning parameters $(\lambda_{1},\lambda_{2})$ in SigL can use well-known cross-validation (CV) technique. In adaptive signal lasso method, They specify a large value for $\lambda$ and only tuning the parameter $\alpha$, which only involved one parameter and will reduce the computation greatly. Furthermore, since $\alpha$ represents the proportion of data compressed to 0 in the interval (0, 1), it should be less than 1 and greater than 0.

III  Signal Lasso with non-convex penalty

III.1 Method

We first consider the linear regression mode (4) with $cov(\epsilon)=\sigma^{2}I_{n}$, the regression parameter $\beta_{j}$ have signal property that being either 0 or 1, therefore the accurate method to estimate $\beta$ is to minimize the following penalized least square

[TABLE]

Different from the previous penalties, here we consider two candidate of penalties

[TABLE]

[TABLE]

Comparing with the additive penalty in signal lasso or adaptive signal lasso, these two penalties are not convex function, however they have appealing property that just one tuning parameter be involved and the $\beta_{j}$ can be shrunk to either 0 or 1 completely as shown later. We will prove that there is no need to tuning parameter in our proposal, as long as the regulation parameter selection is large enough, all the signals $\beta_{j}$ can be shrunk to either 0 or 1, which not only greatly reduce the computational complexity but also fully conforms to the network reconstruction problem in Eq.7.

The product penalty function in Eq.8 has been used to identifying and estimating latent group structures in panel data (Su, Shi, and Phillips (2016)), where they proposed classifier-Lasso (C-Lasso) to shrink individual coefficients to the unknown group-specific coefficients. For minimum penalty function in Eq.9, Tang, Xue, and Qu (2021) used similar functional form to study individualized model selection for different individuals and identify subgroups based on heterogeneous covariates’ effects simultaneously. We use SL_prod to denote signal lasso with product penalty and SL_min for signal lasso with minimum penalty

Some patterns of penalty functions for SigL (panel A), ASigL (panel B), SL_prod (panel C), SL_min (panel D) in two dimensional case (p=2) are given in Fig. 1, respectively. The first column represents the contour for penalty function of four methods, the second column represents the penalty function versus $\beta$, and last column represents the solution of $\beta$ under the orthogonal design. Constraint regions are determined by $PF(x)=c$ for some constant $c=1,0.8,0.6,0.5$. In second column, the red dashed line is the case of lasso method, and in Fig. 1 (e), $PF(x)=|x_{1}|+|\hat{\beta}_{20}||x_{2}-1|$, where green line for case of $\hat{\beta}_{20}=0.2$ and blue line for case of $\hat{\beta}_{20}=0.8$. In Fig. 1 (f) is the solution of $\beta$ under $\lambda_{1}=1,\lambda_{2}=2$. For SL_prod, $\lambda=0.9$, while for SL_min, $\lambda=1.2$.

It is clear that SigL method (panel A) fails to shrink the values of $\hat{\beta}$ that fall in the interval (0.1,0.9) to the correct class and leave an unclassified portion in network reconstruction. ASigL method improve this weakness such that the middle part between 0 and 1 can be shrunk toward to two directions. When one of tuning parameter is set to a large enough value, the slope in middle part of Fig. 1 (f) becomes a perpendicular line so that the parameters can be shrunk to either 0 or 1 (Shi et al. (2023)). The signal lasso method with non-convex penalty in panel C and D of Fig. 1 obviously have this favorable property if tuning parameter $\lambda$ set a large value. The advantage of signal lasso with additive penalty (SigL and ASigL) is that their penalty functions are convex and the optimization problem does not suffer from multiple local minima, and its global minimizer can be efficiently solved. However they need to select one or two tuning parameters which is computationally cost, especially for large size of network. Signal lasso method with non-convex penalty proposed in this paper do not need to tuning parameter as we show later, however the corresponding algorithm to solve the solution have to be developed.

III.2 Updating formula for signal lasso with product penalty

In this section, we propose a computationally very efficient algorithm for solving the two-directional signal lasso problem with product penalty (SL_prod). Based on the coordinate descent method, the key step is to choose a single coordinate say $\beta_{k}$ at the $l$-th iteration and update $\beta_{k}^{(l+1)}$ by applying a univariate minimization over this coordinate, that is for $k=1,\ldots,p$

[TABLE]

where $\|\cdot\|_{2}^{2}$ is the $L^{2}$-norm of a vector. Let ${\bf\epsilon}_{k}^{(l)}=\mathbf{Y}-\sum_{j<k}{\bf x}_{j}\beta_{j}^{(l+1)}-\sum_{j>k}{\bf x}_{j}\beta_{j}^{(l)}$, then $r_{k}^{(l)}=\frac{{\bf x}_{k}^{\prime}{\bf\epsilon}_{k}^{(l)}}{{\bf x}_{k}^{\prime}{\bf x}_{k}}$ is the axis of symmetry of quadratic function part in Eq.(10) that related to variable ${\bf x}_{k}$. Through some simple calculations, the update of $\beta_{k}^{(l+1)}$ is obtained as

[TABLE]

where $\text{I}(\cdot)$ is the indicative function.

III.3 Updating formula for signal lasso with minimum penalty

The computational algorithm for solving the SL_min can also be derived similarly. Based on the coordinate descent method, the key step is to choose a single coordinate say $X_{k}$ at the $l$-th iteration and update $\beta_{k}^{(l+1)}$ by applying a univariate minimization over this coordinate, that is for $k=1,\ldots,p$

[TABLE]

where $\|\cdot\|_{2}^{2}$ is the $L^{2}$-norm of a vector. Let ${\bf\epsilon}_{k}^{(l)}=\mathbf{Y}-\sum_{j<k}{\bf x}_{j}\beta_{j}^{(l+1)}-\sum_{j>k}{\bf x}_{j}\beta_{j}^{(l)}$, then $r_{k}^{(l)}=\frac{{\bf x}_{k}^{\prime}{\bf\epsilon}_{k}^{(l)}}{{\bf x}_{k}^{\prime}{\bf x}_{k}}$ is the axis of symmetry of quadratic function part in Eq.(12) that related to variable ${\bf x}_{k}$. Through some simple calculations, the update of $\beta_{k}^{(l+1)}$ is obtained as

[TABLE]

where $\text{sign}(\cdot)$ is the signum function.

III.4 Selection of tuning parameter

Lemma: *For coordinate descent estimator of signal parameters with non-convex penalty given in (8) or (9), the estimator iteratively obtained from equation (11) or (13) can converge to either 0 or 1 when $\lambda\to\infty$. *

*Proof: * The results can be proved from the formulae given in equation (11) or (13). For estimator using SL_prod when $\lambda\to\infty$ ,

[TABLE]

thus $\hat{\beta}_{k}^{(l+1)}$ tends to 0 for $r_{k}^{(l)}\leq 1/2$ and 1 otherwise. Finally the estimation of $\beta_{k}^{(l+1)}$ will be either 0 or 1. For estimator using SL_min, from equation (13), it is clear that $|r_{k}^{(l)}|-\frac{\lambda}{{\bf x}_{k}^{\prime}{\bf x}_{k}}$ and $|r_{k}^{(l)}-1|-\frac{\lambda}{{\bf x}_{k}^{\prime}{\bf x}_{k}}$ will eventually be less than 0, and finally the estimation of $\beta_{k}^{(l+1)}$ will be either 0 or 1.

This result indicates that we can set $\lambda$ by a large enough values and without tuning $\lambda$ using cross validation, and the signal parameter can be shrunk to 0 or 1 which achieve the purpose of complete classification. The overall iterative coordinate descent algorithm to find the solution can be obtained using these formulae.

IV  The metrics of reconstruction accuracy

To measure the accuracy of network reconstruction, we have to define some metrics for assessing the efficiency of proposed method. Follows the common strategy in network reconstruction in literature (Wang et al. (2011); Han et al. (2015)), the signal parameter $\beta$ can be classified as signal ($\beta=1$) if $\hat{\beta}\in 1\pm 0.1$ and non-signal ($\beta=0$) if $\hat{\beta}\in 0\pm 0.1$, and remains are unclassified. Therefore there exist a portion of unclassified set if estimated $\hat{\beta}$ does not belong to either of two intervals. All quantities involved in classification can be written in Table 1 (see Shi et al. (2023)).

When the predicted class can be completely classified to two classes, the most common used indexes for measuring accuracies include true positive rate (TPR, sensitivity or recall), true negative rate (TNR, or specificity) and precision (Positive prediction value, PPV ), as well as AUROC (the area under the receiver operating characteristic curve) and AUPR (the area under the precision recall curve) (Marbach et al. (2012); Han et al. (2015); Squartini, Fagiolo, and Garlaschelli (2011)), where TPR and TNR are defined by

[TABLE]

where $TP(TN)$ means the number that signal (non-signal ) is correctly identified, $FP$ is the number that non-signal is falsely predicted as signal and $FN$ is the number that signal is falsely predicted as non-signal. Among these measures, Matthews correlation coefficient ($MCC$) is a measure that accounting for the effects of unbalanced number of signal and non-signal parameters (Chicco, 2017; Chicco and Jurman, 2020), and widely used in machine learning research.

To considered the effects of un-classification in Table 1, the success rates for the detection of existing links (SREL) and non-existing links (SRNL) can be defined to study the performance of network reconstruction (Wang et al. (2011); Han et al. (2015); Shi et al. (2021)). Another more efficient measure in this case to study the estimated accuracy when unclassified set remain is the adjusted $MCC(MCCa)$ defined by Shi et al. (2023), which has form of

[TABLE]

where $FNa=FN+UCP$ and $FPa=FP+UCN$. Its values range from -1 to 1 and the large value indicates the good performance. It is clear when un-classified class disappears,$MCCa$ reduces to $MCC$. It is easy to see that MCCa play the similar role as $MCC$ when un-classification appear. Shi et al. (2023) found $MCCa$ performs well in network reconstruction, thus we will use this metric to measure the accuracy of the method in this paper.

V  Numerical Studies

In this section, we use two kinds of model to conduct the simulation studies. One based on standard linear regression model with different assumptions, which including different error distribution and correlated or uncorrelated design matrix. Second one use the network reconstruction model based on game-based evolutionary dynamics and Kuramoto model as used in Shi et al. (2021, 2023).

V.1 Linear regression models

The model generation is given by

[TABLE]

where $\beta_{0}=0$ denote the intercept, $\bf{1}_{n}$ is a $n\times 1$ vector with all elements equal to 1, $\beta_{1}\in R^{p_{1}}$ denote the signal parameter with elements 1, and $\beta_{2}\in R^{p_{2}}$ denote the non-signal parameter with elements 0, and $\epsilon$ is a noise term. Smaller $p_{1}$ is called by sparse signal and larger $p_{1}$ (comparing with n and $p=p_{1}+p_{2}$) is called dense signal. We consider several case in simulations.

• •

Case I: Each design matrix come from standard normal score with mean zeros and variance 1, but columns in in $X$ are dependent with the correlation between $x_{i}$ and $x_{j}$ is $r^{|i-j|}$ with $r=0.5$ (see Tibshirani (1996) and Antoniadis and Fan (2001) for similar settings). The error variable $\epsilon$ is generated from normal distribution with mean zeros and variance $\sigma=0.4,1,2$ respectively. Two combination of sample size and number of independent variable are considered by $n=100,p=30$ and $n=50,p=150$, respectively. The sparse signal with $p1=6$, non-sparse signal with $p1=20$ and dense signal with $p1=100$ are considered for $n<p$.

• •

Case II: The design matrix are the same as Case I, but the error variable are generated from exponential distribution $Exp(\sigma)$, and Gamma distribution $Ga(4,\sigma/2)$, $(n,p)$ are the same as Case I. The exponential and gamma distributions are independently generated. Other settings are same as Case I.

• •

Case III: The design matrix are the same as Case I, but the error variable $\epsilon$ has a AR(1) with correlation $\rho$ ($\epsilon_{i}=\rho\epsilon_{i-1}+u_{i}$, where $u_{i}$ is a Gaussian random variable with mean zero and variance $\sigma^{2}$ ), and $\rho$ takes values of 0.4 or 0.7 respectively.

These three cases corresponds to different sources of the data come from. Case I is for situation of standard linear model with error variable has a independent normal distribution, case II is for case when data has non-normal distribution and case III represents the scenario of correlated observations. After getting the simulated data set, the signal parameter $\beta$ can be classified as $\beta=1$ if $\hat{\beta}\in 1\pm 0.1$ and $\beta=0$ if $\hat{\beta}\in 0\pm 0.1$, and remains are unclassified (Wang et al. (2011); Han et al. (2015)). Some indices in classification problem, including TPR, TNR and Precision rate(PCR), SREL and SRNL can be calculated to quantitatively measure the performance (Marbach et al. (2012); Han et al. (2015); Shi et al. (2021, 2023)). However for saving the space of the paper, we only list several important measures such as MSE, UCR, MCC and MCCa in our simulations.

We first consider the Case I and the results are listed in Table 2, where upper panel is for sparse signal with $n>p$, middle panel is for sparse signal with $n<p$, while bottom panel is about non-sparse signal with $n<p$. The three column in each pane give the results of different variance with $\sigma=0.4,1,2$, respectively. All measures as listed in Shi et al. (2023) are calculated but here we only list the results of MSE, UCR, MCC and MCCa as we see these four measures can reveal the overall accuracy of network reconstruction (refer to Shi et al. (2023)). In Table 2, we calculate the results based on nine methods in which the firs five methods are well known as shrinkage estimators for comparisons, and the last four methods are called as signal-lasso-type estimator because they shrink the parameters towards either 0 or 1. The lasso estimator is used as the initial estimator in last four methods. It is clear all results show that signal-lasso-type estimators superior to the first five well known methods in estimating signal parameter. Therefore we focus on comparison among the last four methods, and have following conclusion from Table 2:

(1) In cases of $n>p$ and the signal is sparse (first panel), we find that the SL_prod and SL_m outperform SigL and ASigL for different variance. In such cases, the SigL is not performs well, especially in case of large variance. When signal is non-sparse or dense, the results are the same (no show here).

(2) For the scenario of $n<p$ and sparse signal ($p_{1}=6$) (second panel), when variance of error variable is not large such as $\sigma=0.4$ and $1$, SL_prod and SL_m clearly outperforms the ASigL and SigL based four measures. However when variance is large such as $\sigma=2$, ASigL performs best with large MCCa and small MSE, although its UCR is not exactly equal to 0. In this case, the SigL performs poorly, especially for large variance.

(3) For the scenario of $n<p$ and non-sparse signal ($p_{1}=20$) (third panel), SL_prod and SL_min performs best and followed by ASigL method. However for large variance with $\sigma=2$, it is observed that ASigL outperforms other three methods, where SigL is not good. The results from second and third panel indicate that AsigL is robust against variance of distribution, which again evidence the conclusion of Shi et al. (2023).

(4) Last panel of Table 2 list the results for case of $n<p$ and dense signal ($p_{1}=100$). It is of interest to find that SigL become better for small variance with smallest MSE and largest MCCa and MCC even its UCR is not closed to zero. When variance tends to large, the advantage of SL_prod and SL_min standout, where they have largest MCCa and small UCR. In this case the SigL have smallest MSE and largest MCC but non-zero UCR. Therefor the performance of SigL is comparable with SL_prod and SL_min for case of dense signal. The ASigL dose not perform well in the circumstance of dense signal.

Next we look at the Case II with Exponential or Gamma distribution are assumed to error variable. We find the results are similar to those of Gussian distribution. We only list the results for exponential distribution. Figure (2) is the plots for MSE, UCR, MCC and MCCa for nine methods, shown that SL_prod and SL_min have larger MCCa and MCC and smaller MSE and UCR than ASigL. However if variance $\sigma$ equal to 2, ASigL is better comparing with SL_prod and SL_min as shown in Figure (3). The box-plot indicates that the SL_prod and SL_min have more outliers in MCCa, MCC and MSE than that of ASigL. The results based on Gamma distribution have the very similar patterns and the results are omitted here.

Finally we consider Case III with linear regression with $n=50,p=150,p_{1}=6$, but error variable have correlated structure of AR(1) distribution with $\rho=0.7$. When variance is smaller with $\sigma=1$, Figure (4) show that SL_prod and SL_min performs better than ASigL with larger MCCa and smaller MSE and UCR. However when $\sigma=2$, Figure (5) clearly indicates that ASigL outperforms SL_prod and SL_m. For non-sparse signal and dense signal, the results are very similar to those based on Gaussian distribution and the results are not shown here.

In summary, we find that signal-lasso-type method performs well than the first five methods, where signal lasso with non-convex penalties seems overall outperform other methods in most cases. In the scenario of sparse and non-sparse signal, for small variance of noise, SL_prod and SL_min perform better than adaptive signal lasso, and they have exact zero UCR for Gaussian and other distribution, even for AR(1) correlated structure. However for large noise variance, ASigL outperform to SL_prod and SL_m, shown a better robustness for the noisy data. However for dense signal, our finding reveal that SigL performs better for small variance, while SL_prod and SL_min performs better for large variance. The results also show that the signal-lasso-type methods are robust against distribution of the data.

V.2 Evolutionary-game-based dynamical model

In an evolutionary game on structured populations, there are multiple players participate the game with their neighbors, where a node represents a player, and a link indicates that two players have a game relationship. The prisoner’s dilemma game (PDG), snowdrift game (SDG), or spatial ultimatum game (SUG) can be used for network reconstruction (Wang et al. (2011); Han et al. (2015); Shi et al. (2021)). We illustrate our method by iterative game dynamics through Monte Carlo simulation. The payoff matrix is given in Eq. (2), here we take a simple structure with $R=1,T=b=1.15$, and $P=S=0$ in our simulation. The game iterates forward in a Monte Carlo manner and player $i$ (the focal player) acquires its fitness (total payoff) $F_{i}$ by playing the game with all its direct neighbors, i.e., $F_{i}=\sum_{j=1}^{N}a_{ij}P_{ij}$. The focal player then randomly picks a neighbor $j$, which similarly acquires its fitness. For simulation cases, player $i$ tries to imitate the strategy of player $j$ with Fermi updating probability $W=1/(1+\exp[(F_{i}-F_{j})/K])$, where $K=0.1$ (Szabó and Fath (2007); Szolnoki, Perc, and Szabó (2012)). By this updating rule of strategies, a series dataset can be generated for each player. To make the model more realistic, we account for a mutation of small rate for the data.

Now $F_{i}=\sum_{j=1,j\neq i}^{N}a_{ij}P_{ij}$ can be written as a linear regression model

[TABLE]

where $Y_{i}=(F_{i}(t_{1}),F_{i}(t_{2}),\cdots,F_{i}(t_{L}))^{\prime}$, $\tilde{\beta}_{i}=(a_{i1},\cdots,a_{iN})^{\prime}$, and $/X_{i}$ has the form of

[TABLE]

Let $Y=(Y^{\prime}_{1},\cdots,Y^{\prime}_{N})^{\prime}$, $\beta=(\tilde{\beta}^{\prime}_{1},\cdots,\tilde{\beta}^{\prime}_{N})^{\prime}$, ${\bf X}=diag({\bf X}_{1},{\bf X}_{2},\cdots,{\bf X}_{N})$, then Eq. (17) can be converted into the general form of Eq. (4).

Here, $L$ represents the length of the data and indicates that there are $L$ interactions between players in the game dynamics. In many real-world scenarios, one needs to uncover a graph’s topology using a limited amount of data. The amount of data is typically scaled by $\Delta=L/N$. In Eq. (17), $\beta$ contains the elements of the connectivity matrix $A=[a_{ij}]$, representing a signal parameter where a value of 1 indicates a connection between two nodes, and 0 indicates no connection. To evaluate the accuracy of the estimator for $\beta$, we plot various metrics as a function of $\Delta=L/N$, which helps to determine how much data is needed to achieve high estimation accuracy. This approach has been widely used in network reconstruction (Wang et al. (2011); Han et al. (2015); Shi et al. (2021, 2023)).

Fig. 6 list the plot of accuracy measures against the scaled data length $\Delta=L/N$ in the PDG model with three networks under different conditions, using the methods of signal lasso, adaptive signal lasso and signal lasso with non-convex penalty. We only list these four methods for comparisons and other methods in Table 2 are ignored since SigL and ASigL outperform other methods as shown in Shi et al. (2023). We consider four cases with different assumptions for noise variances and connectivity density as shown in four panels, and three columns correspond to the results based on Erdös-Rényi (ER) random network, Barabási-Albert (BA) scale-free network and small world (WS) network, respectively.

Panel A in Fig. 6 represents the case with $N=100$, an average degree of $k=6$, and a noise variance of $\sigma^{2}=0.1$. It is evident that SigL performs the worst across the three networks, showing lower reconstruction accuracy (MCCa) and higher UCR at small $\Delta$. Both ASigL and SL_prod perform well, with ASigL slightly outperforming SL_prod at small $\Delta=L/N$. However, this differs from the results in Table 2, where the MCCa values obtained from SL_prod are higher than those of SL_min, particularly in the WS network. Upon examination, we find that this discrepancy arises because several columns in the design matrix $X$ are zero-columns (i.e., all elements in those columns are zero) due to the network’s sparse connectivity. In such cases, the SL_prod method is more robust than SL_min, as demonstrated in the linear regression model in Table 4 of Appendix. Panel B in Fig. 6 uses the same parameters as Panel A but with a higher noise variance of $\sigma^{2}=0.3$. The results are similar to those in Panel A in most cases, where SigL performs worse and ASigL outperforms the other three methods. However, in the WS network, SL_prod and SL_min have higher MCCa values than ASigL when $\Delta>0.6$. Additionally, we also observe that SL_prod outperforms SL_min at small $\Delta=L/N$.

Panel C in Fig. 6 presents the results for the case with $N=100$ and an average degree of $k=20$, representing a non-sparse network. The performances of ASigL, SL_prod, and SL_min are similar, while SigL still performs poorly. Panel D shows the results for a dense signal scenario with an average degree of 70 and a noise variance of $\sigma^{2}=0.1$. In this case, SigL outperforms the other three methods in the BA network. In both the ER and WS networks, the performance of SigL is comparable to the other three methods, showing smaller MCCa when $\Delta<0.5$ but larger MCCa when $\Delta>0.5$. However, it is notable that the UCR for SigL exhibits large non-zero values in some intervals. An interesting observation from Panel D is that SL_prod and SL_min outperform ASigL in both the BA and WS networks.

In conclusion, we find that SigL underperforms compared to the other three methods in both sparse and non-sparse signal scenarios, exhibiting the lowest MCCa and non-zero UCR values, particularly at shorter lengths of $\Delta$ and in situations with high noise variance (as highlighted in Panel B). ASigL slightly outperforms SL_prod, and SL_min demonstrates robustness against noise variance in both sparse and non-sparse signal cases. For dense signals, the performance of SigL varies depending on the network type. While it performs best on BA networks, it is generally on par with SL_prod and SL_min, each excelling under specific conditions. Notably, SigL consistently shows larger non-zero UCR values in some intervals. On the other hand, the UCR values for ASigL, SL_prod, and SL_min are very close to zero across all cases, which is a desirable property. Between SL_prod and SL_min, we recommend using SL_prod due to its greater robustness under specific conditions.

V.3 Kuramoto model in synchronization problem

For problems introduced in Eq. (1), we use the Kuramoto model (Boccaletti et al. (2006); Geier, Timmer, and Fleck (2007); Wu, et al., 2012) to illustrate the reconstruction of the network in a complex system. This model has the following governing equation:

[TABLE]

$i=1,\cdots,N$, where the system is composed of $N$ oscillators with phase $\theta_{i}$ and coupling strength $c$, each of the oscillators has its own intrinsic natural frequency $\omega_{i}$, $a_{ij}$ is the adjacency matrix of a give network and is need to be estimated in network reconstruction. Using the same framework of reference (Shi et al. (2021)), the Euler method can be employed to generate time series with an equal time step $h$. Let $Y_{i}=(y_{i1},\cdots,y_{iL})^{\prime}$, $y_{it}=[\theta_{i}(t+h)-\theta_{i}(t)]/h$, ${\bf X}_{i}=(\phi_{ij}(t))$ is a $L\times N$ matrix, with elements $\phi_{ij}(t)=c\times sin(\theta_{j}(t)-\theta_{i}(t))$ for $t=1,\cdots,L$ and $j=1,\cdots,N$, $\tilde{\beta}_{i}=(a_{i1},\cdots,a_{iN})^{\prime}$, then reconstruction model can be rewritten as

[TABLE]

where ${\bf 1}_{L}$ denote a $L\times 1$ vector with all element 1.

Now we study the performance of the four methods in terms of Kuramoto model in ER, WS and BA networks with $N=100$ and coupling strength $c=10$. The results are listed in Fig. 7, where four panels A-D have the same network structures and noise perturbations as that in Section V.2, but employing Kuramoto model as dynamics. The overall results are similar to the evolutionary-game-based dynamical model. In panel A and B, ASigL is slightly better than SL_prod and SL_min. However the situation is reversed in panel C and D with methods of SL_prod and SL_min being superior to ASigL especially in panel D. Different with case of PD game in Section V.2, the ASigL have non-zero UCR in panel C and D, while UCR from SL_prod and SL_min remain the very small values (close to zeros). In Panel A, B and C, the performance of SL_prod and SL_min almost coincide, but in Panel D, SL_prod outperforms SL_min for larger $\Delta$. Performance of four methods in panel D are very similar to Fig. 6. The signal lasso have largest values of MCCa for BA network, while in ER and WS network, SigL have smaller MCCa when $\Delta<0.5$, but larger MCCa when $\Delta>0.5$. However we find that UCR for SigL have large non-zero values in some intervals.

Summing up the simulation results from linear regression, the evolutionary-game-based dynamical model and the Kuramoto dynamic model, the SL_prod emerges as being on par with the ASigL for both sparse and non-sparse signals. The ASigL displays enhanced robustness against noise variance and can be used for noisy data. Yet, in scenarios involving dense signals, the SigL has some merit for larger $\Delta$, however its UCR will present some non-zeros values in some intervals. For dense signals, ASigL performs worse, while the performance of both SL_prod and SL_min remains stable and exhibits good overall performance. Moreover, the SL_prod slightly outperforms its counterpart with the minimum penalty (SL_min) when the design matrix reveals a zeros column. A significant advantage of both SL_prod and SL_min is the absence of a requirement to tune parameters, translating into substantial computational savings compared to the ASigL. We have calculate the time cost for four methods in different scenarios(see Table 5 in Appendix), which show SL_prod and SL_min are faster and can be recommended for application in some complex situations.

VI  Real Example: A human behavioral data

In this section, we present the results of social network reconstruction using a real data from human behavioral experiment (Li et al. (2018)), where the purpose is to study the impact of the punishment on network reciprocity. A total of 135 participants from Yunnan University of Finance and Economics and Tianjin University of Finance Economics took part in the experiments and three trials are separately designed and carried out. In treatment I, 35 participants from Tianjin university of finance and economics played iterative prisoner’s dilemma game with punishment on the static ring network with four neighbors, where there are 35 nodes and 140 links. For the other two treatments, 100 participants from Yunnan university of finance and economics were invited to participant in the iterated prisoner’s dilemma experiment. The treatment II was implemented on the homogeneous random network with degree of 4, where there are 50 nodes and 200 links. While in treatment III, each player was placed on the heterogeneous random network, in which the degree of half nodes is 3 and the degree of other half nodes is 5, there are 50 nodes and 200 links. The network structure for three treatments are shown in Shi et al. (2021) and the payoff matrices are listed in Table 3. In each round of the treatment, each player played with its direct neighbors to gain their payoff and updated its strategy to optimize its future payoff. The number of interactions in each session was set to 50, and the number of interactions was undisclosed until the sessions ended. In order to solve the network reconstruction problem in these examples, we thus recorded all the strategies and payoffs generated in the experiment and these information is the available data base for the reconstruction. This data set has been used by Shi et al. (2021) and Shi et al. (2023) to assess the network reconstruction of signal lasso and adaptive signal lasso.

We again use this data set to compare the accuracy of reconstruction of network using different methods. The results are summarized in Fig.8, where upper panel A is about results based on treatment I, middle panel B is the results based on treatment II, and bottom panel C is for treatment III. Left column in the Fig.8 give the results using measures of MCC and MSE. We find that SL_prod performs best in Fig.8A(a), while ASigL performs better slightly than other three methods in Fig.8B(a) and Fig.8C(a) based on MCC. MSE of SigL have smaller values than other three methods for small $\Delta$ because other three method can shrink parameters to either 0 or 1 thus will present large deviation once it make a wrong selection.

The right column in Fig.8 consider the influence of unclassified portion and thus use new measures of MCCa and UCR. The results in Fig.8A(b) show that SL_prod performs best, and both SL_prod and SL_min have larger values of MCCa than that of ASigL and SigL. For treatment II and III as shown in Fig.8B(b) and Fig.8C(b), three methods, except for SigL, almost have the same values of MCCa. These three methods have approximately zeros values of UCR, where both SL_prod and SL_min almost exactly equal to zeros even for small $\Delta$. The signal lasso performs worse in these two treatments.

VII  Conclusions

In this paper we propose signal lasso estimation with two non-convex penalties, and make an extensive comparison with other shrinking methods. The results show that the performance of signal-lasso-type methods depends on the variance of noise and the density of network connectivity, but are robust against distribution of the data. We also discuss the performance of signal-lasso-type method in cases of non-sparse or dense signal, which usually appear in world trade web (WTW) (see Squartini, Fagiolo, and Garlaschelli (2011); Shi et al. (2023)).

In the simulations using linear regression models for both sparse and non-sparse signals, we observe that SL_prod and SL_min demonstrate superior or at least comparable performance to other methods in accurately recovering the network structure under conditions of low noise variance. However, with a large noise variance, the adaptive signal lasso performs best, especially when $n<p$, suggesting its robustness against noise perturbations in discerning complex network topologies. In scenarios with dense signals where $n<p$, the signal lasso fares well for minor variances. Yet, as variance increase, SL_prod and SL_min outperform the rest.

Data sourced from evolutionary-game-based dynamical models and the Kuramoto model exhibit greater complexity compared to linear regression models. This data, often sparse within the design matrix, sometimes shows conditional illness or zero-columns. Our analysis under these conditions reveals that the performance of the four methods is critically dependent on network connectivity density. For sparse signals, the adaptive signal lasso (ASigL) is the most proficient, followed closely by SL_prod and SL_min. The signal lasso (SigL) lags, particularly under high variance conditions. In this scenario, ASigL slightly surpasses SL_prod, while SL_prod demonstrates superiority over SL_min in the WS network. As network connectivity density increases, ASigL’s advantages diminish, whereas SigL’s effectiveness improves but is unstable. The results indicate that the performance of SigL and ASigL is unstable in dense networks, while the robust performance of both SL_prod and SL_min remains consistent. A notable advantage of both SL_prod and SL_min is their ability to bypass the need for selecting the tuning parameter $\lambda$, thereby speeding up computational times relative to ASigL (see Table 5 in Appendix for detailed calculations). Furthermore, SL_prod and SL_min consistently shrink the parameter directly to 0 or 1, effectively eliminating any undetermined segments in classification challenges.

In conclusion, our results show that signal lasso with non-convex penalties is effective and fast in estimating signal parameters in linear regression model. We recommend using SL_prod or ASigL in sparse networks, with ASigL being the preferred choice when the data exhibit high variance or significant noise. For non-sparse or dense networks, SL_prod and SL_min are the preferred methods. Between SL_prod and SL_min, we suggest using SL_prod, as it is more robust against certain special conditions that may arise in the design matrix.

Acknowledgement

This work was supported by National Natural Science Foundation of China (NSFC, Grant No. 12271471, 11931015), National Social Science Foundation of China (Grants Nos. 22&ZD158 and 22VRCO49) and Project of ”XingDianYingCai” Plan of Yunnan Province to L. S. and JSPS Postdoctoral Fellowship Program for Foreign Researchers (No. P21374) and an accompanying Grant-in-Aid for Scientific Research to C. S.

Appendix

VII.1 The results for zero column in design matrix of linear regression model

In Table 4, we list the simulation results in linear regression model with design matrix $X$ have 2 zero columns for nine methods: Lasso, Adaptive Lasso, SCAD, MCP, Elastic Net, SigL, ASigL, SL_prod and SLmin in linear regression models. All of the results are averaged over 500 independent realizations, where n is the sample size, $p$ is the number of explanatory variables, $p_{1}$ is the number of signals (number of $\beta=1$). In this special case, it is clear that SL_prod has better performance than that of SL_min.

VII.2 Computational complexity measures (CPU)

In Table 5 , we calculated the CPU time (seconds) for computing estimations of parameters based on SigL, ASigL, SL_prod and SL_min methods, where the initial value for computation is lasso estimators for all methods. It is clear that SL_prod and SL_min methods have obvious advantage in computing speed.

VII.3 Program of algorithm and real dataset

The computation of methods such as lasso, adaptive lasso, SCAD, MCP, and ElasticNet can be found in R-software and has been widely used in scientific research. The R-code using the coordinate descent method for the SigL, ASigL, SL_prod and SL_min, as well as the real dataset are given on GitHub https://github.com/shilei65/adaptive-signal-lasso-code.git.