A unified, computation-friendly framework for penalized principal
machines (P2M), a class of sparse sufficient dimension reduction (SDR)
estimators for regression and binary classification. Principal machines
(PM) estimate the central subspace by solving a family of convex-loss
problems over several cutoffs; their penalized counterparts (P2M) add a
row-group sparsity penalty so that dimension reduction and variable
selection are performed simultaneously. All estimators are fitted by a
single group coordinate descent (GCD) algorithm that accommodates least
squares, logistic, asymmetric least squares, L2-hinge, hinge (support
vector machine, SVM) and quantile losses, together with the least absolute
shrinkage and selection operator (LASSO), the smoothly clipped absolute
deviation (SCAD) penalty and the minimax concave penalty (MCP). Methods are
described in Li, Artemiou and Li (2011)
A unified and computationally efficient R package for sparse sufficient dimension reduction (SDR) using Penalized Principal Machines (P²M) and Group Coordinate Descent (GCD).
The ppmSDR package provides a unified interface and efficient algorithms for sparse sufficient dimension reduction (SDR) in regression and classification. It implements the Penalized Principal Machine (P²M) family, which generalizes the principal support vector machine (PSVM) by allowing a wide range of convex loss functions together with modern sparsity-inducing penalties. Efficient computation is achieved via the Group Coordinate Descent (GCD) and MM-GCD algorithms, making the package scalable to large and high-dimensional data.
ppm() that fits any of ten penalized principal
machine estimators, selected through the loss argument, for both regression
and binary classification.line.search) that guarantees a monotonically non-increasing
penalized objective for the iterative GCD solvers.ppm_tune() for cross-validation of the sparsity parameter.boston, wdbc) and summary()/print()
methods for inspecting the estimated basis and selected variables.# Install from GitHub (requires devtools)
devtools::install_github("c16267/ppmSDR")
ppmSDR/
├── R/
│ ├── ppm.R # ppm(): unified front-end + input validation + dispatch
│ ├── ppm_tune.R # ppm_tune(): dCor-based K-fold cross-validation
│ ├── estimators.R # ten internal P2M solvers (selected via loss=)
│ ├── methods.R # S3 methods: print.ppm, summary.ppm, print.ppm_tune
│ ├── utils-internal.R # internal GCD helpers (thresholding, block-diagonal algebra,
│ │ # penalized objective Q(theta) and step-halving line search)
│ ├── data.R # documentation of the bundled datasets
│ └── ppmSDR-package.R # package-level documentation and imports
├── data/ # boston.rda, wdbc.rda
├── man/ # generated by roxygen2
├── vignettes/ # ppmSDR.Rmd
└── tests/ # unit tests (testthat)
The loss-specific solvers are internal; users access them through the
loss argument of ppm() rather than calling them directly.
| Function | Description |
|---|---|
ppm() |
Unified wrapper to fit any penalized PM via the loss argument |
ppm_tune() |
K-fold cross-validation that selects the sparsity parameter lambda |
summary() |
Estimated basis of the central subspace and the selected variables |
print() |
Compact summary of a fitted ppm / ppm_tune object |
loss argument of ppm())loss |
Method | Response | Algorithm | line.search default |
|---|---|---|---|---|
"lssvm" |
P²LSM | continuous | GCD | not applicable |
"wlssvm" |
P²WLSM | binary | GCD | not applicable |
"logit" |
P²LR | continuous | Iterative GCD | TRUE |
"wlogit" |
P²WLR | binary | Iterative GCD | TRUE |
"asls" |
P²AR | both | Iterative GCD | FALSE |
"l2svm" |
P²L2M | continuous | Iterative GCD | TRUE |
"wl2svm" |
P²WL2M | binary | Iterative GCD | TRUE |
"svm" |
P²SVM | continuous | MM-GCD | FALSE |
"wsvm" |
P²WSVM | binary | MM-GCD | FALSE |
"qr" |
P²QR | both | MM-GCD | FALSE |
Each method supports the "grSCAD", "grMCP" and "grLasso" penalties.
The iterative GCD solvers replace the loss at each iteration by a local
quadratic approximation, which does not by itself guarantee descent of the
penalized P²M objective. With line.search = TRUE, ppm() accepts the
candidate update θ̂ produced by a GCD step only after halving it,
θ ← θ + 2^(−m)(θ̂ − θ), with the smallest m ≥ 0 (up to max.halving = 10)
for which the original penalized objective does not increase, and stops at the
current iterate if no such m exists. The default line.search = NULL turns the
safeguard on for "logit", "wlogit", "l2svm" and "wl2svm" and off
otherwise; line.search = FALSE gives the plain updates. The fitted object
records iter, status, the objective path objective and the number of
halvings n.halving at each iteration.
library(ppmSDR)
## Generate data: the first two predictors form the central subspace
set.seed(1)
n <- 1000; p <- 10
B <- matrix(0, p, 2); B[1, 1] <- B[2, 2] <- 1
x <- MASS::mvrnorm(n, rep(0, p), diag(1, p))
y <- (x %*% B[, 1] / (0.5 + (x %*% B[, 2] + 1)^2)) + 0.2 * rnorm(n)
## Penalized principal least-squares SVM (P2LSM) via the unified interface
fit <- ppm(x, y, loss = "lssvm", penalty = "grSCAD", lambda = 0.01)
round(fit$evectors[, 1:2], 3)
print(fit)
summary(fit, d = 2)
## Penalized principal asymmetric least squares (P2AR) on the same regression data
fit_ar <- ppm(x, y, loss = "asls", penalty = "grSCAD", lambda = 0.05)
round(fit_ar$evectors[, 1:2], 3)
## Binary response with a two-dimensional central subspace: P2AR and the
## penalized principal weighted logistic regression (P2WLR)
y.binary <- sign(x[, 1] + x[, 2]^3 / 3 + 0.2 * rnorm(n))
fit_arb <- ppm(x, y.binary, loss = "asls", penalty = "grSCAD", lambda = 0.2)
summary(fit_arb, d = 2)
fit_wlr <- ppm(x, y.binary, loss = "wlogit", penalty = "grSCAD", lambda = 0.005)
summary(fit_wlr, d = 2)
## Penalized principal logistic regression (P2LR): the step-halving line search
## is on by default; set line.search = FALSE for the plain iterative GCD updates
fit_lr <- ppm(x, y, loss = "logit", penalty = "grSCAD", lambda = 0.01)
fit_lr$status # "converged"
fit_lr$n.halving # number of halvings m_t at each iteration
all(diff(fit_lr$objective) <= 0) # the penalized objective never increases
## Select lambda by cross-validation (set the seed for reproducible folds)
set.seed(1)
cv <- ppm_tune(x, y, loss = "lssvm", d = 2, n.fold = 5,
nlambda = 10, lambda.max = 0.02)
cv$opt.lambda
summary(cv$fit, d = 2)
data(boston) # Boston housing (continuous response: medv)
data(wdbc) # Wisconsin Diagnostic Breast Cancer (binary: diagnosis B/M)
See the package vignette for a worked walk-through:
browseVignettes("ppmSDR")