## Further implementation details
This repository contains the code used for the experiments in the paper Partial Correlation Network Estimation by Semismooth Newton Methods.

### Environment Setup
#### Computing infrastructure
We experimented with Intel® Xeon® Platinum 8568Y+ CPU @ 2.70GHz processors and six NVIDIA RTX 6000 Ada Generation GPUs with 48GB memory. All the implementations were based on Python 3.10, PyTorch 1.13.1, and CUDA 11.7.

#### Create a Conda environment 
To ensure reproducibility, create the Conda environment as follows:
```
conda env create -f environment.yml
conda activate semismooth
```

If the following runtime error occurs:
```
RuntimeError: Error in dlopen: libmkl_intel_lp64.so.1: cannot open shared object file: No such file or directory
```

Fix it by linking MKL libraries manually:
```
cd $CONDA_PREFIX/lib
ln -sf libmkl_intel_lp64.so.2 libmkl_intel_lp64.so.1
ln -sf libmkl_intel_thread.so.2 libmkl_intel_thread.so.1
ln -sf libmkl_core.so.2 libmkl_core.so.1
```

### Repository Structure
```bash
├── data_simul 
│     └── X_1000.npy # Experiment 4.1 dataset
├── algorithm.py
├── Bsemi_func.py
├── generate_data.py
├── lemke.py
├── Experiment1.ipynb
├── objective_func.py
├── update_obj.py
├── update_omega.py
├── util.py
├── environment.yml # Conda environment
├── quadratic_conv.pdf # Figure 1 in the paper
└── Readme.md
```

### Code Overview
Experiment1.ipynb
+ Runs the overall program and orchestrates the algorithm execution for Experiment of Section 4.1.

lemke.py
+ Implements the Lemke algorithm to solve linear complementarity problems (LCPs).

algorithm.py
+ Contains implementations of the B-semismooth Newton and the proximal gradient algorithm.

Bsemi_func.py
+ Defines each iteration step of the B-semismooth Newton method.

generate_data.py
+ Provides utility functions for generating and projecting precision matrix

objective_function.py
+ Includes the objective function used in the optimization problem.

update_omega.py
+ Updates primal and dual variables using parallel computation across multiple GPUs.

update_obj.py
+ Computes and updates the objective function value in parallel using multiple GPUs.

util.py
+ Provides utility functions for reproducibility and sparse matrix manipulation.

### Data Generation and Preprocessing
Given a graph structure (skeleton), we constructed the corresponding ground-truth precision matrix as follows. 

Edge weights were independently sampled from the uniform distribution on $[-1, 1]$ and assigned to the lower-triangular part of the matrix according to the graph’s sparsity pattern. 

The matrix was then symmetrized by adding its transpose, and all diagonal entries were set to one. 

To ensure positive definiteness and preserve signal strength, we iteratively projected the matrix so that 
+ all eigenvalues were bounded below by 0.2, and 
+ all nonzero entries were enforced to have an absolute magnitude of at least 0.1. 

This alternating projection procedure was repeated until convergence or a maximum of 1000 iterations was reached; for details, see generate_data.py 

This repository contains the simulation dataset presented in the paper.

### Step Size Schedule
We used a geometrically decaying step size with initial value `tau_init` and decay factor 0.6 in the proximal gradient method, as this setting had previously shown good empirical performance for the ACCORD problem \citep{lee2024learning}. In all experiments, we set `tau_init` to 0.9.

For the B-semismooth Newton algorithm, the step size was controlled by a linear schedule of the form

```math
\rho = \rho_\text{max} - \frac{\text{current iteration}}{\text{max iteration}} \cdot (\rho_\text{max} - \rho_\text{min})
```
To ensure consistency with the description of the step size behavior in the main paper, we set the corresponding parameters as follows:

+ Experiment 4.1: `rho_max` = 0.7, `rho_min` = 0.4, `num_iters` = 100
