growl-reg 0.1.2
GrOWL regression estimator with OWL, OSCAR, and Lasso variants
Requires Python: >=3.8

[![PyPI version](https://img.shields.io/pypi/v/growl-reg)](https://pypi.org/project/growl-reg/)
[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](LICENSE)
[![Python >=3.8](https://img.shields.io/badge/python-3.8+-blue.svg)](https://www.python.org/)

# 🧮 Group Ordered Weighted $\ell_1$ (GrOWL) Norm

This repository provides a Python implementation of the **Group Ordered Weighted
 $\ell_1$ (GrOWL) Norm** regularization using the **Proximal Operator 
 algorithm** and the **Fast Iterative Shrinkage-Thresholding Algorithm (FISTA)**. 
 It solves the following general optimization problem:

$$
\min_{B} \frac{1}{2n} \lVert Y - XB \rVert_F^2 + 
\sum_i w_i \space \lVert \beta\_{[i], \cdot}\rVert_2 \space, \quad \quad \quad (1)
$$

where  
- $X \in \mathbb{R}^{n \times r}$ is the design matrix,  
- $Y \in \mathbb{R}^{n \times p}$ is the matrix of response variables,  
- $B \in \mathbb{R}^{r \times p}$ is the coefficient matrix to be estimated,  
- $\beta\_{[i], \cdot}$ denotes the $i$-th largest row of $B$ in terms of 
  its $\ell_2$-norm, and  
- $w \in \mathbb{R}^r$ is a vector of non-negative, non-increasing weights.

This regularizaton problem was introduced by Oswal et al. (2016) and it is a 
multi-task ($p > 1$) version of the standard ($p=1$) Ordered Weighted 
$\ell_1$ (OWL) Norm introduced independently by Zeng and Figueiredo (2014a) 
and Bogdan et al. (2013). 

Due to the non-smothness of the GrOWL penalty, a closed-form solution to this 
problem is not available. However, the objective function remains convex, 
allowing the use of efficient proximal optimization algorithms to reliably compute 
the solution. Specifically, it is used the Proximal Gradient Method with 
**Fast Iterative Shrinkage-Thresholding Algorithm (FISTA)** from Beck and 
Teboulle (2009). Readers who are not familiar with proximal algorithms are referred 
to Parikh and Boyd (2013).

---

---

## 📐 Mathematical Background

<!-- The standard ($p = 1$) Ordered Weighted $\ell_1$ (OWL) regularization
problem can be written as

$$
\min_{\beta} \frac{1}{2n} \lVerty - X\beta
Vert_2^2 + \sum_i w_i |\beta\_{[i]}|,
$$

where $w$ is as before but now we have $p=1$, and then $y := Y \in 
\mathbb{R}^{n \times 1}$ and $\beta := B \in \mathbb{R}^{r \times 1}$,
with $\beta\_{[i]}$ being the $i$-th largest component of $\beta$. -->

Due to non-smothness of the penalty term in (1), this optimization problem 
has no closed-form solution. Proximal operator algorithms is employed to solve
it. The proximal operator of the GrOWL norm is given by

$$
\mathrm{prox}_G(V) = \mathrm{arg min}_B \space \frac{1}{2} \lVert B - V \rVert_F^2 + 
\sum_i w_i \space \lVert \beta\_{[i], \cdot} \rVert_2.
$$

The proximal operator of GrOWL is solved in terms of the proximal operator of 
the standard OWL (when $p=1$) norm, denoted by $\mathrm{prox}\_{\Omega_w}$. We thus 
have the following result:

---

**Theorem 4 from Oswal et al. (2016).**  
Let $\tilde{v}_i = \lVert v\_{i,\cdot}\rVert$ for $i = 1, ..., p$. Then 
$\mathrm{prox}_G(V) = \hat{V}$, where the $i$-th row of $\hat{V}$ is given by

$$
\hat{\mathbf{v}}\_{i,\cdot} = 
\left(\mathrm{prox}\_{\Omega_w}(\tilde{\mathbf{v}}) \right)_i \times
\frac{\mathbf{v}\_{i,\cdot}}{\lVert \mathbf{v}\_{i,\cdot} \rVert}.
$$

---

The formulation of $\mathrm{prox}\_{\Omega_w}$ is given in equation (24) of
Zeng and Figueiredo (2014b):

$$
\mathrm{prox}\_{\Omega_w}(\mathbf{\tilde{v}}) = 
\mathrm{sign}(\mathbf{\tilde{v}}) \odot \left( \mathbf{P}(|\mathbf{\tilde{v}}|)^T 
\mathrm{proj}\_{\mathbb{R}\_+^n} \left( \mathrm{proj}\_{\mathcal{K}\_m} 
(|\mathbf{\tilde{v}}|\_{\downarrow} - \mathbf{w}) \right) \right),
$$

where
- $\mathrm{sign}(\mathbf{\tilde{v}})$ denotes the elementwise sign 
  of vector $\mathbf{\tilde{v}}$.
- $\odot$ is the Hadamard (elementwise) product.
- $\mathbf{P}(|\mathbf{\tilde{v}}|)$ is the permutation matrix that 
  sorts the absolute values $|\mathbf{\tilde{v}}|$ in non-increasing order,
  i.e., $|\mathbf{v}|\_{\downarrow} = \mathbf{P}(|\mathbf{\tilde{v}}|) 
  |\mathbf{\tilde{v}}|$.
- $\mathrm{proj}\_{\mathcal{K}_m}$ is the Euclidean projection 
  onto the monotone cone $\mathcal{K}\_m =$ {$\mathbf{x} 
  \in \mathbb{R}^n : x_1 \geq x_2 \geq \cdots \geq x_n$}, 
  implemented using the Pool Adjacent Violators (PAV) algorithm.
- $\mathrm{proj}\_{\mathbb{R}_+^n}$ is the Euclidean projection
  onto the nonnegative orthant, i.e., it replaces negative values by zero 
  (clipping).
- $\mathbf{w}$ is a weight vector satisfying $w_1 \geq w_2 \geq
  \cdots \geq w_n \geq 0$.
- $|\mathbf{\tilde{v}}|\_{\downarrow}$ denotes the absolute values of
  $\mathbf{\tilde{v}}$ sorted in non-increasing order.

We use **FISTA** (Beck and Teboulle, 2009), which is an accelerated first-order 
method designed for problems of the form:

$$
\min_{B} f(B) + g(B),
$$

where $f$ is convex and differentiable with Lipschitz continuous gradient, 
and $g$ is convex (possibly non-smooth) with a proximal operator that can be 
computed efficiently.

In our case:
- $f(B) := \frac{1}{2n} \lVert Y - XB \rVert_2^2$ is the smooth loss,
- $g(B) := \sum_i w_i \space \lVert \beta\_{[i], \cdot}\rVert_2$ is the GrOWL
  penalty.

FISTA proceeds by alternating between gradient descent steps on $f$ and proximal 
steps on $g$, with a Nesterov-type momentum update to accelerate convergence. Each
iteration consists of:

1. **Gradient step:**

$$
V^{(k)} = Z^{(k)} - \frac{1}{L} \nabla f(Z^{(k)}),
$$

where $L$ is the Lipschitz constant of $\nabla f$, computed as 
$L = \lVert X \rVert_2^2/n$, where $\lVert X \rVert_2$ denotes the spectral norm of the matrix $X$.

2. **Proximal step using the GrOWL operator:**

$$
B^{(k+1)} = \mathrm{prox}_G (V^{(k)}),
$$

which is implemented as described earlier, using the Pool Adjacent Violators (PAV) 
algorithm for isotonic regression and restoring the original signs and order.

3. **Nesterov momentum step:**

$$
t_{k+1} = \frac{1}{2} \left(1 + \sqrt{1 + 4t_k^2} \right), \quad
Z^{(k+1)} = B^{(k+1)} + \left( \frac{t_k - 1}{t_{k+1}} \right) (B^{(k+1)} - B^{(k)}).
$$

The algorithm continues until convergence is detected, based on one of three 
user-defined stopping criteria:
- Absolute change in objective value,
- Relative change in objective value,
- Frobenius norm of the difference between successive iterates.

This FISTA procedure is implemented in the function `growl_fista()` inside the file 
'fista_solver.py' in the codebase. The function handles flexible weight vector definitions 
(manual or parameterized via `lambda_1`,  `lambda_2`, and `ramp_size`) and returns the 
estimated coefficient matrix along with the cost history. 

The proximal operators evaluations are implemented in the functions 'prox_owl()' and
'prox_growl()' inside the file 'prox_operator.py'.

---

---

## 🗂 Repository Structure

Below are the important modules in this project and their functionalities:

1. **`__init__.py`**  
   This file is part of the `growl/` module and exposes the `GrowlRegressor` class
   as the main interface for the package. It enables clean imports such as:
   
   ```python
   from growl import GrowlRegressor
   ```

3. **`base.py`**  
   Contains the main class `GrowlRegressor`, a `scikit-learn` compatible estimator that
   implements GrOWL regression. This class provides:
   
   - `.fit(X, Y)` to estimate coefficients using the GrOWL penalty
   - `.predict(X)` for in-sample or out-of-sample predictions
   - Integration with `GridSearchCV`
   - Optional centering of `X` and `Y` when `fit_intercept=True`
   - Storage of the coefficient matrix `coef_` and optimization history `cost_history_`

5. **`prox_operator.py`**  
   Implements proximal operators required for optimization:
   - `prox_owl(v, w)`: Evaluate the proximal operator for the OWL penalty.
   - `prox_growl(V, w)`: Evaluate the proximal operator for the GrOWL penalty.

6. **`fista_solver.py`**  
   Implements the FISTA-based optimization routine used to solve the GrOWL regularized
   least-squares problem. This module includes:
   
   - `growl_fista(...)`: A solver using Nesterov’s acceleration
   - Weight vector construction based on `lambda_1`, `lambda_2`, and `ramp_size`
   - Convergence monitoring based on cost, relative cost, or solution change
   - Optional scaling of the objective function to improve numerical stability

8. **`growl_example.py`**  
   Located in the `examples/` folder, this script demonstrates the usage of
   the `GrowlRegressor`:
   
   - Grid search over hyperparameters (`lambda_1`, `lambda_2`, `ramp_size`)
   - Visual comparisons between:
     - True vs estimated coefficients
     - GrOWL vs MultiTaskLasso (for pooled regression)
     - GrOWL (OWL style) vs Lasso (for standard regression)
   - Plots showing grouping behavior and coefficient shrinkage

   To run the example, use:
   ```bash
   python examples/growl_example.py
   ```

---

---

## ⚙️ Setup

**Install the repository:**

   ```bash
   pip install growl_reg
   ```

---

---

## 📚 References

Beck, A. and Teboulle, M. "A fast iterative shrinkage-thresholding algorithm
for linear inverse problems", _SIAM Journal on Imaging Sciences, vol. 2, no. 1,
pp. 183–202_, 2009.

Bogdan, J., Berg, E., Su, W. and Candes, E. "Statistical 
estimation and testing via the ordered $\ell_1$ norm", arXiv preprint 
[arxiv:1310.1969v2](https://arxiv.org/abs/1310.1969) 2013.

Oswal, U., Cox, C., Ralph, M. A. L., and Rogers, T., Nowak, R., 2016. 
"Representational Similarity Learning with Application to Brain Networks". 
_Proceedings of the 33 rd International Conference on Machine Learning, 
New York, NY, USA, 2016. JMLR: W\&CP volume 48_.

Parikh, Neal and Boyd, Stephen. "Proximal algorithms". _Foundations and Trends
in optimization_, 1(3):123–231, 2013.

Zeng, X. and Figueiredo, M, 2014a. "Decreasing Weighted Sorted $\ell_1$ 
Regularization". arXiv preprint 
[arXiv:1404.3184v1](https://arxiv.org/abs/1404.3184), 2014.

Zeng, X. and Figueiredo, M, 2014b. "The ordered weighted $\ell_1$ norm - atomic 
formulation, projections, and Algorithms". arXiv preprint 
[arXiv:1409.4271v5](https://arxiv.org/abs/1409.4271), 2014.

---

## 📑 Citation

If you use `growl_reg` in your work, please cite it as:

Matheus Lopes Carrijo. "GrOWL Regression Estimator (Python package)." 2025. 
Available at: https://github.com/matheuscarrijo/growl_reg


Or use the following BibTeX entry:

```bibtex
@misc{carrijo2025growl,
  author       = {Carrijo, M. L.},
  title        = {GrOWL Regression Estimator (Python Package)},
  year         = {2025},
  howpublished = {https://github.com/matheuscarrijo/growl_reg},
  note         = {Version 0.1.2}
}
```