*! version 1.1.0 20nov2025

Title

clspl2 -- Convex Least Squares Programming (CLSP) is a modular
two-step estimator for solving underdetermined, ill-posed, or structurally
constrained least-squares problems, based on the L2 norm (i.e., implementing the
alpha = 1 special case)

Syntax

clspl2
string
[if] [in]
[, options]

Problem definition
string Problem type, a structural template for matrix
construction. One of: ap/tm - allocation (tabular) matrix
problem (AP), cmls/rp - constrained-model least squares
(regression) problem, or anything else - general CLSP problem
(user-defined constraints(string) and/or
model(string).

Right-hand side
b(strings:string) Right-hand side vector,
varname or name of an existing matrix. Must
have as many rows as the design matrix A = [C, S\
M, Q] in Az = b (please note that for AP, it
should start with row sums). Required.

Left-hand side
constraints(strings:string) Block C of the design matrix
A = [C, S\ M, Q], varlist
or name of an existing matrix. If C and/or M are
provided, the matrix A is constructed accordingly (please note
that for APs/TMs, C is constructed automatically and known values
are specified in M).

slackvars(strings:string) Block S of the design matrix
A = [C, S\ M, Q], varlist
or name of an existing matrix.

model(strings:string) Block M of the design matrix
A = [C, S\ M, Q], varlist
or name of an existing matrix. If C and/or M are
provided, the matrix A is constructed accordingly (please note
that for APs/TMs, C is constructed automatically and known values
are specified in M).

m(#) Number of rows in the problem X ∈ ℝ^ in
APs/TMs.

p(#) Number of columns in the problem X ∈ ℝ^ in
APs/TMs.

i(#) Grouping size for row sum constraints in
APs/TMs.

j(#) Grouping size for column sum constraints in
APs/TMs.

zerodiagonal Enforce structural zero diagonals in
APs/TMs.

Estimation (first step)
r(#) Number of refinement iterations for the
pseudoinverse-based estimator z-hat (default 1).

z(name) A symmetric idempotent matrix (projector) defining
the subspace for Bott-Duffin pseudoinversion. By default, the identity
matrix is used, reducing the Bott-Duffin inverse to the Moore-Penrose
case.

rcond(real) Regularization parameter for the Moore-Penrose
and Bott-Duffin inverses, providing numerically stable inversion and
ensuring convergence of singular values. If -1, an automatic
tolerance equal to tolerance(real) is
applied. If set to a positive number in (0, 1),
it specifies the relative cutoff below which small singular values are
treated as zero (default 0, i.e. no regularization).

tolerance(real) Convergence tolerance for NRMSE change
between refinement iterations. If -1, an automatic tolerance
equal to the square root of machine epsilon is applied. If set to
a non-negative number in [0, 1), it specifies a
custom tolerance for the estimator (default -1).

iterationlimit(#) Maximum number of iterations allowed in the
refinement loop (default 50).

Estimation (second step)
nofinal Suppress the second step, a second pseudoinverse
estimation to refine z-hat.

Other
condtolerance(real) Singular-value cutoff for the custom
condition number function. If set to a missing value, the
implementation uses an internal relative cutoff of 1e-14.

varname and varlist in b(string),
constraints(string),
slackvars(string), and
model(string) may not contain time-series operators
or factor variables and are filtered by [if] [in] with missing values being
excluded by default; see tsvarlist and fvvarlist.

by, collect, fp, rolling, statsby, and xi
are allowed; see prefix.

weights are not allowed; see weight.

See Postestimation for features available after estimation.

Description

clspl2 implements the
Convex Least Squares Programming (CLSP) estimator, a modular two-step
convex optimization framework, capable of addressing ill-posed and
underdetermined problems. After reformulating a problem in its canonical form,
Az = b -- where A^(r) = [C, S\
M^(r), Q^(r)], C - constraints matrix, S - sign slack
matrix, M^(r) - model matrix, and Q^(r) - zero matrix or reverse-sign
slack matrix for residuals -- Step 1 yields an iterated (if r > 1) minimum-norm
least-squares estimate zhat = (A_Z^†)b on a constrained subspace
defined by a symmetric idempotent Z (reducing to the Moore-Penrose
pseudoinverse when Z = I) -- where (A_Z^†)^(r) =
((Z(A^(r))^+A^(r)Z)^†)Z(A^(r))^+ is the
Bott-Duffin inverse with a left projector P_Z^L (reduced to
(A^†)^(r) = (((A^(r))^+A^(r))^†)(A^(r))^+ when
Z = I). The optional Step 2 corrects zhat via a second
pseudoinverse estimation z^* = zhat + (A^(r))^+
(A^(r)(A^(r))^+)^† (b - A^(r)zhat), guaranteeing
an L2-norm-based Minimum-Norm BLUE (MNBLUE) unique solution (corresponding to
alpha = 1 in the estimator).

Special cases:

1. allocation problems (APs) (or, for flow variables, tabular
matrix problems, TMs) - in most cases, underdetermined problems involving
matrices X ∈ ℝ^ to be estimated, subject to known row and column
sums, manifested as a block submatrix [((I_ ⊗ 1_i) ⊗
1_p)^+, (1_m ⊗ (I_

⊗ 1_j)^+)]^+ in
C and a [row sums\ column sums^+] subvector in b;

2. constrained-model least squares (CMLS) (or, more generally, regression
problems, RPs) - (over)determined problems to be estimated, subject to,
among others, linear(ized) data matrix D-related (in)equality constraints,
such as the i-th difference or shift (lag/lead) of specific signs, manifested
as submatrices of transformed D in C and b.

Numerical stability of A:

For general cases, the conditioning of A -- with respect to the problem
constraints (i.e., rows) in [C, S] -- is analyzed from the point of
view of its sensitivity to (a) iterations in Step 1 (using a decomposition
derived from the inequality kappa(A) ≤ kappa(B) kappa([C,
S]), where B is (A^(r))^†[C, S], the variable part of
A in terms of iterations), and (b) dropping a single row in [C,
S]. To assess the angular alignment of rows in [C, S] as
vectors, a cosine-based mean root square alignment (RMSA) is constructed, with
values close to ±1 indicating near-collinearity and potential ill-conditioning
and values near 0 implying near-orthogonality and improved numerical stability
(RMSA can be calculated for both [C, S] and [C, S] reduced
by one row accompanied by changes in conditioning and solutions).

Goodness of fit:

The CLSP estimator employs goodness-of-fit statistics robust to the (potential)
underdeterminedness of solved problems: partial R^2 (i.e., the adaptation of
R^2 to the CLSP structure, where the term "partial" stands for statistics
related to the subvector x^* in z^*, corresponding to the columns of
C and M, as opposed to y^* in z^*, corresponding to the
columns of S and Q), normalized RMSE (NRMSE), t-test for the mean of
the NRMSE -- with the sample generated via bootstrap resampling of residuals
(quicker) or from repeated reestimations based on randomly generated b
(more precise, Monte Carlo simulation) -- and a diagnostic interval derived
from the condition number of A and the norms of residuals and b
(constructed from the condition number inequality, where the residuals are
assumed to be the perturbation, and perturbation for each value of z^* is
considered to be uniform).

Postestimation

For total RMSA and sensitivity of [C, S] to dropping one row:

. clspl2 estat corr
[, options]

reset Force recomputation of all diagnostic values instead of
replaying previously stored results.

threshold(real) If positive, limits the output to constraints
with e(rmsa_i) ≥ threshold(real).

bar Display a vertical bar chart for e(rmsa_i) or for the
stored result named in matrix(string). If the
stored result is a matrix, a panel of bar charts is produced, with the
number of columns controlled by ncols(#). Mutually
exclusive with hbar.

hbar Display a horizontal bar chart for e(rmsa_i) or for
the stored result named in matrix(string). If
the stored result is a matrix, a panel of horizontal bar charts is
produced, with the number of columns controlled by
ncols(#). Mutually exclusive with bar.

matrix(string) Stored result in e() to graph
(default rmsa_i). Requires bar or hbar.

ncols(#) Number of columns in the graph panel when plotting
a matrix (default 1). Requires bar or hbar.

For the bootstrap/Monte Carlo t-test for the mean of NRMSE:

. clspl2 estat ttest
[, options]

reset Force a new generation of the t-test sample instead of
replaying previously stored results.

samplesize(#) Size of the Monte Carlo simulated sample under
H0 (default 50).

seed(#) Optional random seed to override the default
(default 123456789).

distribution(string) Distribution for generating synthetic
b vectors. One of: normal, uniform, or laplace
(default normal). Requires simulate.

partial Run the t-test on the partial NRMSE: during
simulation, the entries in C are preserved and the entries in
M are simulated.

simulate Perform a parametric Monte Carlo simulation by
generating synthetic right-hand side vectors b. Otherwise, execute
a nonparametric bootstrap procedure on residuals without
re-estimation.

level(#) Set confidence level (default 95)

Remarks

For Python users, there is a standalone Python implementation https://pypi.org/project/pyclsp/ with the same functionality.

Likewise, for R users, there is an R equivalent https://cran.r-project.org/web/packages/rclsp/.

For a Python-based Stata alternative, consult clsp.

To ensure cross-platform reproducibility, all CLSP implementations use a
modified condition number function based on singular values, with a
relative cutoff equal to condtolerance * the largest singular value.

Examples

* AP (TM), based on a symmetric input-output table, with 10% of
* known values
. clear
. local seed = 123456789
. set seed `seed'

* sample (dataset)
. local m = 20
. local p = 20
. mata: X_true = rnormal(`m', `p', 0, 1 )
. mata: X_true = abs( X_true + X_true') / 2.0
. mata: st_matrix("X_true", X_true)
. mata: idx = runiformint(trunc(0.1 * (`m' * `p')), 1, 1, `m' * `p')

* model
. mata: M = I(`m' * `p')[idx,.]
. mata: st_matrix("M", M)
. mata: b_row = quadrowsum(X_true)
. mata: b_col = quadcolsum(X_true)'
. mata: b_val = vec(X_true)[idx,.]
. mata: st_matrix("rhs", b_row\b_col\b_val)
. quiet clspl2 ap, b(rhs) model(M) m(`m') p(`p')
. quiet estat corr, bar

* results
. display "true X:"
. matlist X_true
. display "X_hat: "
. matlist e(x)
. mata: ss_res = sum((st_matrix("X_true") - st_matrix("e(x)")) :^2)
. mata: ss_tot = sum((st_matrix("X_true") :- mean(st_matrix("X_true"))) :^2)
. mata: st_numscalar("r2_user", 1 - ss_res/ss_tot)
. di as res %-32s "R2_user_defined:" %5.4f r2_user
. estat ttest, samplesize(30) seed(`seed') distribution(normal)

* CMLS (RP), based on known stationary points for y = Dx + e,
* x to be estimated (x is labeled beta in regression analysis)
. clear
. local seed = 123456789
. set seed `seed'

* sample (dataset)
. local k = 10
. local p = 6
. local c = 1
. set obs `k'
. mata: D = J(`k', 1, 1)
. mata: D = D, rnormal(`k', `p'-1, 0, 1)
. mata: b_true = rnormal(`p', 1, 0, 1)
. mata: b_true = ( b_true :/ colsum( b_true)) :* `c'
. mata: st_matrix("b_true", b_true)
. mata: e = rnormal(`k', 1, 0, 1)
. mata: y = D * b_true + e
. mata: st_addvar("float", (("D" :+ strofreal(1..`p')), "y"))
. mata: st_store(1::`k', 1..`p'+1, (D, y))

* model
. mata: st_matrix("b", `c' \ J(`k'-2, 1, 0) \ J(`k'-1, 1, 0) \ y)
. mata: d1D = D[2::rows( D), .] :- D[1::rows( D)-1, .]
. mata: d2D = d1D[2::rows(d1D), .] :- d1D[1::rows(d1D)-1, .]
. mata: st_matrix("C", J(1, `p', 1) \ d2D \ d1D )
. mata: d1y = y[2::rows( y), .] :- y[1::rows( y)-1, .]
. mata: d2y = d1y[2::rows(d1y), .] :- d1y[1::rows(d1y)-1, .]
. mata: st_matrix("S", J(1, `k'-2, 0) \ diag(sign(d2y)) \ J(`k'-1, `k'-2, 0))
. quiet clspl2 cmls, b(b) c(C) s(S) mod(D*) r(1)

* results
. display "true beta (x_M):"
. matlist b_true'
. display "beta hat (x_M hat):"
. matlist e(x)'
. estat ttest, samplesize(30) seed(`seed') distribution(normal) partial

Stored results

clspl2 stores the following in e():

(developers may be interested in . ereturn list, all)

Scalars

e(N) number of observations

e(r) number of refinement iterations performed in the first step

e(rcond) regularization parameter for the Moore-Penrose and
Bott-Duffin inverses

e(tolerance) convergence tolerance for NRMSE change between
refinement iterations

e(alpha) regularization parameter (weight) in the final convex
program of the CLSP estimator; for clspl2 always 1 (added for
compatibility with the Python-based clsp)

e(kappaC) spectral kappa() for [C, S]

e(kappaB) spectral kappa() for B = A[C,
S]^†

e(kappaA) spectral κ() for A

e(r2_partial) R^2 for M in A

e(nrmse) root mean square error calculated from A and
normalized by standard deviation (NRMSE)

e(nrmse_partial) root mean square error calculated from M in
A and normalized by standard deviation (NRMSE)

Macros

e(cmdline) command as typed

e(cmd) clspl2

e(estat_cmd) program used to implement estat

e(title) estimator title

e(final) presence of the second step in estimation, True or
False

e(seed) random seed in the Monte Carlo-based t-test

e(distribution) distribution for generating simulated b
vectors in the Monte Carlo-based t-test

Matrices

e(zhat) vector of the first-step estimate

e(z) vector of the final solution. If the second step is
disabled, it equals zhat

e(x) m×p matrix or vector containing the variable component of
z

e(y) vector containing the slack component of z

e(z_lower) lower bound of the diagnostic interval (confidence
band) based on kappa(A) for z

e(z_upper) upper bound of the diagnostic interval (confidence
band) based on kappa(A) for z

e(x_lower) lower bound of the diagnostic interval (confidence
band) based on kappa(A) for x

e(x_upper) upper bound of the diagnostic interval (confidence
band) based on kappa(A) for x

e(y_lower) lower bound of the diagnostic interval (confidence
band) based on kappa(A) for y

e(y_upper) upper bound of the diagnostic interval (confidence
band) based on kappa(A) for y

Functions

e(sample) marks estimation sample

In addition to the above, the following is stored in r():

Matrices

r(rmsa_i) summary statistics for RMSA after dropping one row in
[C, S] (i indicates each row)

r(rmsa_dkappaC) summary statistics for changes in e(kappaC)
after dropping one row in [C, S]

r(rmsa_dkappaB) summary statistics for changes in e(kappaB)
after dropping one row in [C, S]

r(rmsa_dkappaA) summary statistics for changes in e(kappaA)
after dropping one row in [C, S]

r(rmsa_dnrmse) summary statistics for changes in e(nrmse)
after dropping one row in [C, S]

r(rmsa_dzhat) summary statistics for changes in e(zhat)
(L2 norm) after dropping one row in [C, S]

r(rmsa_dz) summary statistics for changes in e(z)
(L2 norm) after dropping one row in [C, S]

r(rmsa_dx) summary statistics for changes in e(x)
(L2 or Frobenius norm) after dropping one row in [C, S]

r(z_lower) summary statistics for the lower bound of the
diagnostic interval (confidence band) based on kappa(A) for
z

r(z_upper) summary statistics for the upper bound of the
diagnostic interval (confidence band) based on kappa(A) for
z

r(x_lower) summary statistics for the lower bound of the
diagnostic interval (confidence band) based on kappa(A) for
x

r(x_upper) summary statistics for the upper bound of the
diagnostic interval (confidence band) based on kappa(A) for
x

r(y_lower) summary statistics for the lower bound of the
diagnostic interval (confidence band) based on kappa(A) for
y

r(y_upper) summary statistics for the upper bound of the
diagnostic interval (confidence band) based on kappa(A) for
y

Note that results stored in r() are updated when the command is replayed and
will be replaced when any r-class command is run after the estimation command,
including estat corr and estat ttest.

Author

Ilya Bolotov

Prague University of Economics and Business

Prague, Czech Republic

ilya.bolotov@vse.cz

Thanks for citing this software and my works on the topic:

Bolotov, I. (2026). CLSPL2: Stata module providing Convex Least Squares
Programming (CLSP) is a modular two-step estimator for solving
underdetermined, ill-posed, or structurally constrained least-squares
problems. Available from
https://ideas.repec.org/c/boc/bocode/s459695.html.

Bolotov, I. (2025). CLSP: Linear Algebra Foundations of a Modular Two-Step
Convex Optimization-Based Estimator for Ill-Posed Problems. Mathematics,
13, 3476. Available from https://doi.org/10.3390/math13213476.