Method Description¶

This method is developed by our team. We propose a doubly-penalized fusion framework that simultaneously estimates smooth functional curves and discovers row (sample) and column (covariate) clusters in a single optimization. The method is distribution-free and automatically determines the number of clusters. The key steps are listed below.

1. Functional Representation

   Each observed curve is expanded with B-spline basis functions

   \[g_{i,j}(t) \approx \mathbf{U}_p(t)^\top \boldsymbol{\beta}_{i,j},\]

   yielding coefficient vectors \(\boldsymbol{\beta}_{i,j} \in \mathbb{R}^p\).

2. Objective Function

   Minimize

   \[L(\boldsymbol{\beta}) = \frac{1}{2}\| \mathbf{Y} - \mathbf{U} \boldsymbol{\beta} \|_2^2 + \frac{\gamma_1}{2} \boldsymbol{\beta}^\top \mathbf{M} \boldsymbol{\beta} + \gamma_2 \sum_{i_1 < i_2} \rho_\tau \left( \| \boldsymbol{\beta}^{(r)}_{i_1} - \boldsymbol{\beta}^{(r)}_{i_2} \|_2 \right) + \gamma_2 \sqrt{\frac{N}{q}} \sum_{j_1 < j_2} \rho_\tau \left( \| \boldsymbol{\beta}^{(c)}_{j_1} - \boldsymbol{\beta}^{(c)}_{j_2} \|_2 \right),\]

where

• \(\gamma_1\) smooths each curve via a second-order difference penalty (matrix \(\mathbf{M}\));

• \(\gamma_2\) induces fusion of sample / covariate coefficients via a concave penalty \(\rho_\tau\) (MCP or SCAD);

• \(\boldsymbol{\beta}^{(r)}_i\) and \(\boldsymbol{\beta}^{(c)}_j\) denote the stacked coefficients for sample \(i\) and covariate \(j\) respectively.

3. Tuning

   • Step 1 (smoothness): fix \(\gamma_2 = 0\) and choose optimal \(\gamma_1\) via BIC.

   • Step 2 (fusion): fix \(\gamma_1\) and choose \(\gamma_2\) via a second BIC.

4. ADMM Algorithm

   • Primal updates: closed-form ridge-type solution for \(\boldsymbol{\beta}\).

   • Dual updates: closed-form soft-thresholding for sample- and covariate-level fusion variables.

   • Convergence: monitored through primal and dual residuals.

5. Statistical Guarantees

   • Consistency: under mild regularity, the oracle estimator (with known clusters) achieves

     \[\| \hat{g}_{k_r,k_c} - g^*_{k_r,k_c} \| = O_P \left( \left( \frac{p \log(Nq)}{|G_{\min}|} \right)^{1/2} \right).\]

   • Clustering consistency: a local minimizer equals the oracle estimator with probability → 1 as \(N, q \to \infty\).

Function¶

This method provides three core functions: pf_sim_data, pf_bifunc and FDPlot.pf_fdplot. In this section, we detail their respective usage, as well as parameters, output values and usage examples for each function.

pf_sim_data(n, T, nknots, order, seed = 123)

from BiFuncLib.simulation_data import pf_sim_data
pf_simdata = pf_sim_data(n = 60, T = 10, nknots = 3, order = 3, seed = 123)['data']

pf_bifunc(data, nknots, order, gamma1, gamma2, opt = False, theta = 1, tau = 3, max_iter = 500, eps_abs = 1e-3, eps_rel = 1e-3)

from BiFuncLib.simulation_data import pf_sim_data
pf_simdata = pf_sim_data(n = 60, T = 10, nknots = 3, order = 3, seed = 123)['data']
pf_result = pf_bifunc(pf_simdata, nknots = 3, order = 3, gamma1 = 0.023, gamma2 = 3,
                      theta = 1, tau = 3, max_iter = 500, eps_abs = 1e-3, eps_rel = 1e-3)
# Performs opt
pf_result_opt = pf_bifunc(pf_simdata, nknots = 3, order = 3, gamma1 = [0.023, 0.025], gamma2 = [2,3],
                          opt = True, theta = 1, tau = 3, max_iter = 500, eps_abs = 1e-3, eps_rel = 1e-3)

FDPlot(result).pf_fdplot()

from BiFuncLib.pf_bifunc import pf_bifunc
from BiFuncLib.simulation_data import pf_sim_data
from BiFuncLib.FDPlot import FDPlot
pf_simdata = pf_sim_data(n = 60, T = 10, nknots = 3, order = 3, seed = 123)['data']
pf_result = pf_bifunc(pf_simdata, nknots = 3, order = 3, gamma1 = 0.023, gamma2 = 3,
                     theta = 1, tau = 3, max_iter = 500, eps_abs = 1e-3, eps_rel = 1e-3)
FDPlot(pf_result).pf_fdplot()