*! version 1.1.0 20aug2025

Title

lppinvl2 -- Linear Programming via Regularized Least Squares
(LPPinv) is a modular two-step estimator for solving underdetermined, ill-posed,
or structurally constrained linear and quadratic programming problems, based on
the L2 norm (i.e., implementing the alpha = 1 special case)

Requirements

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

lppinvl2
string
[if] [in]
[, options]

Problem definition
string One of: nonneg - enforce nonnegativity constraints,
the standard setting in real-world linear and quadratic programming
problems, or free - enforce no constraints.

Right-hand side
bub(strings:string) Right-hand side vector for inequality
constraints, varname or
name of an existing matrix. Please note that at least one of
bub(string) or
beq(string) must be provided, matched by
a corresponding aub(string) or
aeq(string).

beq(strings:string) Right-hand side vector for equality
constraints, varname or
name of an existing matrix. Please note that at least one of
bub(string) or
beq(string) must be provided, matched by
a corresponding aub(string) or
aeq(string).

Left-hand side
aub(strings:string) Matrix for inequality constraints
aub(string) ⋅ x ≤
bub(string), varlist
or name of an existing matrix. Please note that at least one of
aub(string) or
aeq(string) must be provided, matched by
a corresponding bub(string) or
beq(string).

aeq(strings:string) Matrix for equality constraints
aeq(string) ⋅ x =
beq(string), varlist
or name of an existing matrix. Please note that at least one of
aub(string) or
aeq(string) must be provided, matched by
a corresponding bub(string) or
beq(string).

Prior information
lowerbound(numlist) Lower bounds on cell values. If a single
value is given (please use missing values . or .a-z for
-∞), it is applied to all variables. Under nonneg, negative
values in lowerbound(numlist) will produce an
error.

upperbound(numlist) Upper bounds on cell values. If a single
value is given (please use missing values . or .a-z for
∞), it is applied to all variables. Under nonneg, negative
values in upperbound(numlist) will produce an
error.

replacevalue(numlist) Final replacement value (real or missing
. or .a-z) for any cell in the solution matrix that violates
the bounds specified in lowerbound(numlist)
and upperbound(numlist) by more than a
tolerance equal to tolerance(real).

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 bub(string),
beq(string),
aub(string), and
aeq(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

lppinvl2 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).

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:

. lppinvl2 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:

. lppinvl2 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/pylppinv/ with the same functionality.

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

For a Python-based Stata alternative, consult lppinv.

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

* LPRLS/QPRLS, based on an underdetermined and potentially infeasible
* problem
. clear
. local seed = 123456789
. set seed `seed'

* sample (dataset)
. mata: A_ub = rnormal(50, 500, 0, 1)
. mata: st_matrix("A_ub", A_ub)
. mata: A_eq = rnormal(25, 500, 0, 1)
. mata: st_matrix("A_eq", A_eq)
. mata: b_ub = rnormal(50, 1, 0, 1)
. mata: st_matrix("b_ub", b_ub)
. mata: b_eq = rnormal(25, 1, 0, 1)
. mata: st_matrix("b_eq", b_eq)
. quiet lppinvl2 free, aub(A_ub) aeq(A_eq) bub(b_ub) beq(b_eq) r(1)

* results
. display "x hat (x_M hat):"
. matlist e(x), format(%9.4f)
. estat ttest, samplesize(30) seed(`seed') distribution(normal)

Stored results

lppinvl2 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 lppinvl2 always 1 (added
for compatibility with the Python-based tmpinv)

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) lppinvl2

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) 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 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. (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.