clsp: Convex Least Squares Programming (CLSP) estimator.

Description

The Convex Least Squares Programming (CLSP) estimator solves underdetermined, ill-posed, or structurally constrained least-squares problems using a modular two-step approach. The first step computes a pseudoinverse-based estimate, and the second step applies a convex correction (Lasso, Ridge, or Elastic Net) to ensure numerical stability, constraint enforcement, and interpretability.

Usage

clsp(
  problem = "",
  C = NULL,
  S = NULL,
  M = NULL,
  b = NULL,
  m = NULL,
  p = NULL,
  i = 1L,
  j = 1L,
  zero_diagonal = FALSE,
  r = 1L,
  Z = NULL,
  rcond = FALSE,
  tolerance = NULL,
  iteration_limit = NULL,
  final = TRUE,
  alpha = NULL,
  ...
)

Value

An object of class "clsp" representing the fitted Convex Least Squares Programming (CLSP) model. The object is a named list containing all initialized fields and solver results. Class-specific methods such as summary.clsp(), corr.clsp(), and ttest.clsp() can be used to extract, analyze, and summarize the results.

Arguments

problem

character scalar, optional Structural template for matrix construction. One of:

'ap' or 'tm': allocation or tabular matrix problem.
'cmls' or 'rp': constrained modular least squares or RP-type.
'' or other: general CLSP problems (user-defined $\boldsymbol{C}$ and/or $\boldsymbol{M}$).

C, S, M

numeric matrix or NULL Blocks of the constraint matrix $\mathbf{A} = \begin{bmatrix} \mathbf{C} & \mathbf{S} \\ \mathbf{M} & \mathbf{Q} \end{bmatrix}$. If $\boldsymbol{C}$ and/or $\boldsymbol{M}$ are provided, the matrix $\boldsymbol{A}$ is constructed accordingly. If both are NULL and $\boldsymbol{A}$ is not yet defined, an error is raised.

b

numeric vector or NULL Right-hand-side vector. Must have as many rows as $\boldsymbol{A}$. Required.

m, p

integer scalar or NULL Dimensions of $\mathbf{X} \in \mathbb{R}^{m \times p}$, relevant for allocation problems ('ap').

i, j

integer scalar, default = 1 Grouping sizes for row and column-sum constraints in AP problems.

zero_diagonal

logical scalar, default = FALSE If TRUE, enforces structural zero diagonals via identity truncation.

r

integer scalar, default = 1 Number of refinement iterations for the pseudoinverse-based estimator. When $r > 1$, the slack block $\boldsymbol{Q}$ is updated iteratively to improve feasibility in underdetermined or ill-posed systems.

Z

numeric matrix or NULL A symmetric idempotent matrix (projector) defining the subspace for Bott–Duffin pseudoinversion. If NULL, the identity matrix is used, reducing to the Moore–Penrose case.

rcond

numeric scalar or logical scalar, default = FALSE Regularization parameter for the Moore–Penrose and Bott–Duffin inverses, providing numerically stable inversion and ensuring convergence of singular values. If TRUE, an automatic tolerance equal to tolerance is applied. If set to a numeric value, it specifies the relative cutoff below which small singular values are treated as zero.

tolerance

numeric scalar or NULL, default = NULL Convergence tolerance for NRMSE change between iterations.

iteration_limit

integer scalar or NULL, default = NULL Maximum number of iterations allowed in the refinement loop.

final

logical scalar, default = TRUE If TRUE, a convex programming problem is solved to refine zhat. The resulting solution $\boldsymbol{z}$ minimizes a weighted $\ell_1/\ell_2$ norm around $\widehat{\mathbf{z}}$ subject to $\mathbf{A}\mathbf{z} = \mathbf{b}$.

alpha

numeric scalar, numeric vector, or NULL, default = NULL Regularization parameter:

$\alpha = 0$: Lasso ($\ell_1$ norm)
$\alpha = 1$: Ridge ($\ell_2$ norm)
$0 < \alpha < 1$: Elastic Net. If a numeric scalar is provided, that value is used after clipping to $[0,1]$. If a numeric vector is provided, each candidate is evaluated via a full solve, and the $\alpha$ with the smallest NRMSE is selected. If NULL, $\alpha$ is chosen automatically according to $$\alpha = \min\left(1,\, \frac{\mathrm{NRMSE}_{\alpha=0}} {\mathrm{NRMSE}_{\alpha=0} + \mathrm{NRMSE}_{\alpha=1} + \mathrm{tolerance}}\right)$$.

...

Optional. Additional arguments passed to the CVXR solver backend.

Details

This estimator unifies pseudoinverse-based least squares with convex programming correction. The pseudoinverse step computes an initial solution $\mathbf{z}^{(r)}$ iteratively via the Moore–Penrose or Bott–Duffin inverse. The convex step then refines $\boldsymbol{z}$ by minimizing a mixed $\ell_1/\ell_2$ norm under equality constraints $\mathbf{A}\mathbf{z} = \mathbf{b}$. The method supports allocation problems (AP), constrained modular least squares (CMLS), and general CLSP formulations.

Examples

Run this code

if (FALSE) {
  ## Example: CMLS (RP) estimation with stationary-point constraints

  set.seed(123456789)

  # sample (dataset)
  k  <- 500L                                        # number of observations
  p  <- 6L                                          # number of regressors
  c0 <- 1                                           # sum of coefficients

  D        <- matrix(NA_real_, nrow = k, ncol = p)
  D[, 1]   <- 1.0                                   # constant
  D[, 2:p] <- matrix(rnorm(k * (p - 1)), k, p - 1)

  b_true   <- rnorm(p)
  b_true   <- (b_true / sum(b_true)) * c0           # normalize to sum = c

  e        <- matrix(rnorm(k), ncol = 1)
  y        <- D %*% b_true + e

  # build blocks for CLSP (CMLS)
  b <- rbind(
      matrix(c0, ncol = 1),                         # sum of coefficients
      matrix(0,  nrow = k - 2, ncol = 1),
      matrix(0,  nrow = k - 1, ncol = 1),
      matrix(y,  ncol = 1)
  )

  C <- rbind(
      matrix(1, nrow = 1, ncol = p),                # row of ones
      diff(D, differences = 2),                     # 2nd differences
      diff(D, differences = 1)                      # 1st differences
  )

  # diagonal sign-matrix for 2nd differences
  S <- rbind(
      matrix(0, nrow = 1,      ncol = k - 2),
      diag(sign(diff(as.numeric(y), differences = 2))),
      matrix(0, nrow = k - 1,  ncol = k - 2)
  )

  # model
  model <- rclsp::clsp(
      problem = "cmls",
      b       = b,
      C       = C,
      S       = S,
      M       = D,
      r       = 1L,                                 # no refinement
      alpha   = 1.0                                 # MNBLUE solution
  )

  # results
  print("true beta (x_M):")
  print(round(b_true, 4))

  print("beta hat (x_M hat):")
  print(round(model$x, 4))

  print(model)

  # bootstrap t-test
  tt <- rclsp::ttest(
      model,
      sample_size  = 30L,
      seed         = 123456789L,
      distribution = rnorm,
      partial      = TRUE
  )

  print("Bootstrap t-test:")
  print(tt)
}