*! version 1.1.0 20nov2025

Title

tmpinvl2 -- Tabular Matrix Problems via Pseudoinverse Estimation
(TMPinv) is a modular two-step estimator for solving underdetermined, ill-posed,
or structurally constrained allocation 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

tmpinvl2
string
[if] [in]
[, options]

Problem definition
string One of: full/general - full model or # # -
reduced model with given matrix dimensions (rows columns), so that the
problem is estimated as a set of reduced problems constructed from
contiguous submatrices of the original table. For example, reduced =
(6, 6) implies 5×5 data blocks with 1 slack row and 1 slack column each
(edge blocks may be smaller). The minimum is 3 3.

reduced(#) Position within the reduced model, from which the
respective reduced solution's results are loaded into e(), with
possible values ranging between -number of reduced problems
(equivalent to 1) and number of reduced problems (equivalent
to -1) (default -1). Requires reduced model.

Right-hand side
brow(strings:string) Right-hand side vector of row totals,
varname or name of an existing matrix. Please
note that both brow(string) and
bcol(string) must be provided. Required.

bcol(strings:string) Right-hand side vector of column totals,
varname or name of an existing matrix. Please
note that both brow(string) and
bcol(string) must be provided. Required.

bval(strings:string) Right-hand side vector of known cell
values, varname or
name of an existing matrix.

Left-hand side
slackvars(strings:string) Block S, of shape
(m + p)×(m + p), of the design matrix A = [C, S\ M,
Q], varlist or
name of an existing matrix. A diagonal sign slack (surplus) matrix
with entries in {0,±1}: 0 enforces equality (==
brow(string) or
bcol(string)), 1 enforces a
lower-than-or-equal (≤) condition, -1 enforces a
greater-than-or-equal (≥) condition. The first m diagonal entries
correspond to row constraints, and the remaining p to column
constraints. Please note that, in the reduced model, S is
ignored: slack behavior is derived implicitly from block-wise marginal
totals.

model(strings:string) Block M of the design matrix
A = [C, S\ M, Q], varlist
or name of an existing matrix. Each row defines a linear
restriction on the flattened solution matrix. The corresponding
right-hand side values must be provided in
bval(string). This block is used to encode
known cell values. Please note that, in the reduced model, M must
be a unique row subset of an identity matrix (i.e.,
diagonal-only). Arbitrary or non-diagonal model matrices cannot be
mapped to reduced blocks, making the model infeasible.

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.

symmetric Enforce symmetry of the estimated solution matrix
as: e(X) = 0.5 * (e(X) + e(X)'). Applies to
e(X) only. For proper model symmetry, please add explicit symmetry
constraints to M in a full-model solve instead of using this
option.

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 m×p cells.

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 m×p cells.

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.

saving(strings:string) path for the reduced-model .ster
file; by default, results are saved as ./tmpinvl2.ster.

varname and varlist in brow(string),
bcol(string), bval(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

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

tmpinvl2 focuses on a special case of the CLSP estimation, 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, - yielding a L2-norm-based Minimum-Norm
BLUE (MNBLUE) unique solution.

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:

. tmpinvl2 estat corr
[, options]

reduced(#) Position within the reduced model, from which the
respective reduced solution's results are loaded into e(), with
possible values ranging between -number of reduced problems
(equivalent to 1) and number of reduced problems (equivalent
to -1) (default e(model_i)).

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:

. tmpinvl2 estat ttest
[, options]

reduced(#) Position within the reduced model, from which the
respective reduced solution's results are loaded into e(), with
possible values ranging between -number of reduced problems
(equivalent to 1) and number of reduced problems (equivalent
to -1) (default e(model_i)).

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

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

For a Python-based Stata alternative, consult tmpinv.

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: P = order(runiform(`m' * `p', 1), 1)
. mata: idx = (1::(`m' * `p'))[P][1..floor(0.1 * (`m' * `p'))]

* model
. mata: M = I(`m' * `p')[idx,.]
. mata: st_matrix("M", M )
. mata: b_row = quadrowsum(X_true)
. mata: st_matrix("brow", b_row)
. mata: b_col = quadcolsum(X_true)'
. mata: st_matrix("bcol", b_col)
. mata: b_val = vec(X_true)[idx,.]
. mata: st_matrix("bval", b_val)
. local l = 0
. local u = .
. local LHS "brow(brow) bcol(bcol) bval(bval) model(M) sym"
. quiet tmpinvl2 full, `LHS' lower(`l') upper(`u') r(1)

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

* AP (TM), based on a trade matrix, with a zero diagonal and 20% of
* known values (computation is slow)
. clear
. local seed = 123456789
. set seed `seed'

* sample (dataset)
. local m = 40
. local p = 40
. mata: X_true = rnormal(`m', `p', 0, 1 )
. mata: _diag(X_true, 0)
. mata: st_matrix("X_true", X_true)
. mata: P = order(runiform(`m' * `p', 1), 1)
. mata: idx = (1::(`m' * `p'))[P][1..floor(0.2 * (`m' * `p'))]

* model
. mata: M = I(`m' * `p')[idx,.]
. mata: st_matrix("M", M )
. mata: b_row = quadrowsum(X_true)
. mata: st_matrix("brow", b_row)
. mata: b_col = quadcolsum(X_true)'
. mata: st_matrix("bcol", b_col)
. mata: b_val = vec(X_true)[idx,.]
. mata: st_matrix("bval", b_val)
. local l = 0
. local u = .
. local LHS "brow(brow) bcol(bcol) bval(bval) model(M) zerod sym"
. quiet tmpinvl2 20 20, `LHS' lower(`l') upper(`u') r(1)

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

Stored results

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

e(estat_cmd) program used to implement estat

e(title) estimator title

e(model_i) position within the reduced model

e(model_n) number of reduced problems

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(X) m×p matrix of the full solution

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 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). TMPINVL2: Stata module solving Tabular Matrix Problems via
Pseudoinverse Estimation (TMPinv) is a modular two-step estimator for
solving underdetermined, ill-posed, or structurally constrained allocation
problems using least-squares and convex optimization. Available from
https://ideas.repec.org/c/boc/bocode/s459697.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.