cirls.fit: Constrained Iteratively Reweighted Least-Squares
In cirls: Constrained Iteratively Reweighted Least Squares

cirls.fit

R Documentation

Constrained Iteratively Reweighted Least-Squares

Description

Fits a generalized linear model with linear constraints on the coefficients through a Constrained Iteratively Reweighted Least-Squares (CIRLS) algorithm. This function is the constrained counterpart to glm.fit and is meant to be called by glm through its method argument. See details for the main differences.

Usage

cirls.fit(x, y, weights = rep.int(1, nobs), start = NULL,
  etastart = NULL, mustart = NULL, offset = rep.int(0, nobs),
  family = stats::gaussian(), control = list(), intercept = TRUE,
  singular.ok = TRUE)

Arguments

`x`, `y`	`x` is a design matrix and `y` is a vector of response observations. Usually internally computed by glm.
`weights`	An optional vector of observation weights.
`start`	Starting values for the parameters in the linear predictor.
`etastart`	Starting values for the linear predictor.
`mustart`	Starting values for the vector or means.
`offset`	An optional vector specifying a known component in the model. See model.offset.
`family`	The result of a call to a family function, describing the error distribution and link function of the model. See family for details of available family functions.
`control`	A list of parameters controlling the fitting process. See details and cirls.control.
`intercept`	Logical. Should an intercept be included in the null model?
`singular.ok`	Logical. If `FALSE`, the function returns an error for singular fits.

Details

This function is a plug-in for glm and works similarly to glm.fit. In addition to the parameters already available in glm.fit, cirls.fit allows the specification of a constraint matrix Cmat with bound vectors lb and ub on the regression coefficients. These additional parameters can be passed through the control list or through ... in glm but not both. If any parameter is passed through control, then ... will be ignored.

The CIRLS algorithm is a modification of the classical IRLS algorithm in which each update of the regression coefficients is performed by a quadratic program (QP), ensuring the update stays within the feasible region defined by Cmat, lb and ub. More specifically, this feasible region is defined as

⁠lb <= Cmat %*% coefficients <= ub⁠

where coefficients is the coefficient vector returned by the model. This specification allows for any linear constraint, including equality ones.

Specifying constraints

The package includes several mechanisms to specify constraints. The most straightforward is to pass a full matrix to Cmat with associated bound vectors in lb and ub. In this case, the number of columns in Cmat must match the number of coefficients estimated by glm. This includes all variables that are not involved in any constraint, potential expansion such as factors or splines for instance, as well as the intercept. By default lb and ub are set to 0 and Inf, respectively, but any bounds are possible. When some elements of lb and ub are identical, they define equality constraints. Setting lb = -Inf and ub = Inf disable the constraints.

To avoid pre-constructing potentially large and complex Cmat objects, the arguments Cmat and constr can be combined to conveniently specify constraints for the coefficients. More specifically, Cmat can alternatively take a named list of matrices to constrain only specific terms in the model. The argument constr provides a formula interface to specify built-in common constraints. The documentation of buildCmat provides full details on how to specify constraints along with examples.

Quadratic programming solvers

The function cirls.fit relies on a quadratic programming solver. Several solver are currently available.

"quadprog" (the default) performs a dual algorithm to solve the quadratic program. It relies on the function solve.QP.
"osqp" solves the quadratic program via the Alternating Direction Method of Multipliers (ADMM). Internally it calls the function solve_osqp.
"coneproj" solves the quadratic program by a cone projection method. It relies on the function qprog.

Each solver has specific parameters that can be controlled through the argument qp_pars. Sensible defaults are set within cirls.control and the user typically doesn't need to provide custom parameters. "quadprog" is set as the default being generally more reliable than the other solvers. "osqp" is faster but can be less accurate, in which case it is recommended to increase convergence tolerance at the cost of speed.

Value

A cirls object inheriting from the class glm. At the moment, two non-standard methods specific to cirls objects are available: vcov.cirls to obtain the coefficients variance-covariance matrix and confint.cirls to obtain confidence intervals. These custom methods account for the reduced degrees of freedom resulting from the constraints, see vcov.cirls and confint.cirls.

An object of class cirls includes all components from glm objects, with the addition of:

`Cmat`, `lb`, `ub`	the constraint matrix, and lower and upper bound vectors. If provided as lists, the full expanded matrix and vectors are returned.
`active.cons`	vector of indices of the active constraints in the fitted model.
`inner.iter`	number of iterations performed by the last call to the QP solver.
`etastart`	the initialisation of the linear predictor `eta`. The same as `etastart` when provided.
`singular.ok`	the value of the `singular.ok` argument.

Any method for glm objects can be used on cirls objects. Several methods specific to cirls are available: vcov.cirls to obtain the coefficients variance-covariance matrix, confint.cirls to obtain confidence intervals, and logLik.cirls to extract the log-likelihood with appropriate degrees of freedom.

References

Goldfarb, D., Idnani, A., 1983. A numerically stable dual method for solving strictly convex quadratic programs. Mathematical Programming 27, 1–33. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/BF02591962")}

Meyer, M.C., 2013. A Simple New Algorithm for Quadratic Programming with Applications in Statistics. Communications in Statistics - Simulation and Computation 42, 1126–1139. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/03610918.2012.659820")}

Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S., 2020. OSQP: an operator splitting solver for quadratic programs. Math. Prog. Comp. 12, 637–672. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/s12532-020-00179-2")}

Examples

####################################################
# Simple non-negative least squares

# Simulate predictors and response with some negative coefficients
set.seed(111)
n <- 100
p <- 10
betas <- rep_len(c(1, -1), p)
x <- matrix(rnorm(n * p), nrow = n)
y <- x %*% betas + rnorm(n)

# Define constraint matrix (includes intercept)
# By default, bounds are 0 and +Inf
Cmat <- cbind(0, diag(p))

# Fit GLM by CIRLS
res1 <- glm(y ~ x, method = cirls.fit, Cmat = Cmat)
coef(res1)

# Same as passing Cmat through the control argument
res2 <- glm(y ~ x, method = cirls.fit, control = list(Cmat = Cmat))
identical(coef(res1), coef(res2))

####################################################
# Increasing coefficients

# Generate two group of variables: an isotonic one and an unconstrained one
set.seed(222)
p1 <- 5; p2 <- 3
x1 <- matrix(rnorm(100 * p1), 100, p1)
x2 <- matrix(rnorm(100 * p2), 100, p2)

# Generate coefficients: those in b1 should be increasing
b1 <- runif(p1) |> sort()
b2 <- runif(p2)

# Generate full data
y <- x1 %*% b1 + x2 %*% b2 + rnorm(100, sd = 2)

#----- Fit model

# Create constraint matrix and expand for intercept and unconstrained variables
Ciso <- diff(diag(p1))
Cmat <- cbind(0, Ciso, matrix(0, nrow(Ciso), p2))

# Fit model
resiso <- glm(y ~ x1 + x2, method = cirls.fit, Cmat = Cmat)
coef(resiso)

# Compare with unconstrained
plot(c(0, b1, b2), pch = 16)
points(coef(resiso), pch = 16, col = 3)
points(coef(glm(y ~ x1 + x2)), col = 2)

#----- More convenient specification

# Cmat can be provided as a list
resiso2 <- glm(y ~ x1 + x2, method = cirls.fit, Cmat = list(x1 = Ciso))

# Internally Cmat is expanded and we obtain the same result
identical(resiso$Cmat, resiso2$Cmat)
identical(coef(resiso), coef(resiso2))

#----- Adding bounds to the constraints
# Difference between coefficients must be above a lower bound and below 1
lb <- 1 / (p1 * 2)
ub <- 1

# Re-fit the model
resiso3 <- glm(y ~ x1 + x2, method = cirls.fit, Cmat = list(x1 = Ciso),
  lb = lb, ub = ub)

# Compare the fit
plot(c(0, b1, b2), pch = 16)
points(coef(resiso), pch = 16, col = 3)
points(coef(glm(y ~ x1 + x2)), col = 2)
points(coef(resiso3), pch = 16, col = 4)

cirls documentation built on Sept. 13, 2025, 1:09 a.m.