CLSP: Linear Algebra Foundations of a Modular Two-Step Convex Optimization-Based Estimator for Ill-Posed Problems

1. Introduction

As existing literature on optimization — for example, Nocedal and Wright [1], pp. 1–9 and Boyd and Vandenberghe [2], pp. 1–2 — points out, constrained optimization plays an important role not only in natural sciences, where it was first developed by Euler, Fermat, Lagrange, Laplace, Gauss, Cauchy, Weierstrass, and others, but also in economics and its numerous applications, such as econometrics and operations research, including managerial planning at different levels, logistics optimization, and resource allocation. The history of modern constrained optimization began in the late 1940s to 1950s, when economics, already heavily reliant on calculus since the marginalist revolution of 1871–1874 [3], pp. 43–168, and the emerging discipline of econometrics, standing on the shoulders of model fitting [4] (pp. 438–457), internalized a new array of quantitative methods referred to as programming after the U.S. Army’s logistic jargon. Programming — including classical, linear, quadratic, and dynamic variants — originated with Kantorovich in economic planning in the USSR and was independently discovered and developed into a mathematical discipline by Dantzig in the United States (who introduced the bulk of terminology and the simplex algorithm) [5] [6], pp. 1–51. Since the pivotal proceedings published by Koopmans [7], programming became an integral part of mathematical economics [8,9,10], in both microeconomic analysis (an overview is provided in Intriligator and Arrow [11], pp. 76–91) and macroeconomic modeling, such as Leontief’s input-output framework outlined in Koopmans [7], pp. 132–173 and, in an extended form, in Dorfman et al. [12], pp. 204–264.

Calculus, model fitting, and programming are, however, subject to strong assumptions, including function smoothness, full matrix ranks, well-posedness, well-conditionedness, feasibility, non-degeneracy, and non-cycling (in the case of the simplex algorithm), as summarized by Nocedal and Wright [1], pp. 304–354, 381–382 and Intriligator and Arrow [11], pp. 15–76. These limitations have been partially circumvented by developing more strategy-based than analysis-based, computationally intensive constraint programming algorithms; consult Frühwirth and Abdennadher [13] and Rossi et al. [14] for an overview. Still, specific cases of constrained optimization problems in econometric applications, such as incorporating the textbook two-quarter business-cycle criterion into a general linear model [15] or solving underdetermined linear programs to estimate an investment matrix [16], remain unresolved within the existing methodologies and algorithmic frameworks. In response, this text develops a modular two-step convex programming methodology to formalize such problems as a linear system of the form $\mathbf{A}\mathbf{z}=\mathbf{b}$, including linear(ized) constraints $\mathbf{C}\mathbf{x}+\mathbf{S}{\mathbf{y}}^{*}={\mathbf{b}}_{1:k}$, with $k\le m$, where $\mathbf{A}\in {\mathbb{R}}^{m\times n}$ is a partitioned matrix, $\mathbf{C}\in {\mathbb{R}}^{k\times p}$ and $\mathbf{S}\in {\mathbb{R}}^{k\times (n-p)}$ are its submatrices, $\mathbf{b}\in {\mathbb{R}}^m$ is a known vector, and $\mathbf{z}\in {\mathbb{R}}^n$ is the vector of unknowns, which includes both the target solution subvector $\mathbf{x}\in {\mathbb{R}}^p$ and an auxiliary slack-surplus subvector ${\mathbf{y}}^{*}\in {\mathbb{R}}^{n-p}$, introduced to accommodate inequalities in constraints. The focus is on its linear algebra foundations with corresponding theorems and lemmas, followed by a discussion of numerical stability, goodness of fit, and potential applications.

The text is organized as follows. Section 2 summarizes the historical development and recent advances in convex optimization, motivating the formulation of a new methodology. Section 3 presents the formalization of the estimator, while Section 4 develops a sensitivity analysis for $\mathbf{C}$. Section 5 introduces goodness-of-fit statistics, including the normalized root mean square error (NRMSE) and its sample-mean test. Section 6 presents special cases, and Section 7 reports the results of a Monte Carlo experiment and solves simulated problems via CLSP-based Python 3 modules. Section 8 concludes with a general discussion. Throughout, the following notation is used: bold uppercase letters (e.g., $\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{G},\mathbf{M},\mathbf{S}$) denote matrices; bold lowercase letters (e.g., $\mathbf{b},\mathbf{x},\mathbf{y},\mathbf{z}$) denote column vectors; italic letters (e.g., $i,j,k,m,n,p,\lambda ,\omega$) denote scalars; the transpose of a matrix is denoted by the superscript ⊤ (e.g., ${\mathbf{X}}^{\top }$); the inverse by $-1$ (e.g., ${\mathbf{Q}}^{-1}$); generalized inverses are indexed using curly braces (e.g., ${\mathbf{G}}^{\{1,2\}}$); the Moore-Penrose pseudoinverse is denoted by dagger (e.g., ${\mathbf{A}}^{†}$); ${\mathsf{\ell }}_p$ norms, where $1\le p\le \infty$, are denoted by double bars (e.g., ${\parallel ·\parallel }_2$); and condition numbers by $\kappa (·)$. All functions, scalars, vectors, and matrices are defined over the real numbers ($\mathbb{R}$).

2. Historical and Conceptual Background of the CLSP Framework

Convex optimization — formally, ${\mathbf{x}}^{*}=arg{min}_{\mathbf{x}\in {\mathbb{R}}^n}f\left(\mathbf{x}\right)$, subject to convex inequality constraints $g_i\left(\mathbf{x}\right)\le 0$ and affine equality constraints $h_j\left(\mathbf{x}\right)=0$ [2], pp. 7, 127–129, encompassing, among others, least squares, linear programming, and quadratic programming — has its roots in differential calculus, a concept developed independently by Newton — who was influenced by Wallis and Barrow — in 1665–1666 (published posthumously in 1736), and by Leibniz, who drew among others on Pascal, in 1673–1675 (published in 1684 and 1686) [17], and, in its constrained form, by Lagrange in the 1760s. Calculus-based optimization dominated analytical methods in mathematics and physics until the 19th century, during which period it was superseded in specific applications by model-fitting techniques, driven by empirical inconsistencies in classical mechanics, particularly in celestial motion [18], pp. 11–61. Subsequent contributions by Euler, Laplace, Cauchy, and Weierstrass helped shape the modern algorithmic treatment of calculus-based optimization, as presented in, for example, Nocedal and Wright [1], pp. 10–29 and, for economic applications, in Intriligator and Arrow [11], pp. 53–72 and Ok [19], pp. 670–745. In (micro)economic theory, Lagrangian optimization was introduced by Edgeworth, who applied it in 1877 to utility maximization, in 1879 to utilitarian welfare distribution, and in 1881 to general equilibrium via the contract curve [20]. Since the 19th century, optimizing differentiable scalar-valued functions has become the foundation for the fitting techniques in mathematical statistics [4], pp. 494–510.

Stigler [18], pp. 11–61 traced the pre-19th century history leading to the invention of the least squares method (named after the original French term la méthode des moindres carrés, henceforth abbreviated as LS) by Legendre in 1805, Adrain in 1808-1809, and independently by Gauss in 1809 (although the latter claimed to have been using the method since 1795; a rigorous discussion of the priority dispute between Gauss and Legendre, a well-known case study in the history of scientific priority, is provided in Plackett [21] and Stigler [22] [23], pp. 320–331). The method’s formalism evolved from Gauss’s, Yule’s, and Laplace’s pre-matrix notations to Pearson’s, Fisher’s, Frisch’s, Aitken’s (who is believed to anticipate the Gauss-Markov theorem [24]), and Bartlett’s use of Cayley’s 1850s and subsequent methodologies, with the modern $\mathbf{y}=\mathbf{Xb}$ notation, where $\mathbf{X}$ is a matrix and $\mathbf{y},\mathbf{b}$ are vectors, introduced around 1950 and already in use by Durbin and Watson, among others [25]. LS became an essential tool at the dawn of modern econometrics and its application in linear regression (ordinary least squares, OLS) intensified with the creation of the Cowles Commission in the 1930s–1940s U.S. and the works of Fisher, Koopmans, Haavelmo, Stone, Prest, and others [26], pp. 117–145, [27], pp. 7–12. Contrary to the stricter assumptions of pure calculus-based optimization, LS (a) enabled modeling over discrete or irregular data, and (b) incorporated residuals, moving beyond exact curve fitting, while still ensuring minimum variance among unbiased linear estimators under the Gauss–Markov theorem.

Early efforts to estimate the normal equations ${\mathbf{X}}^{\top }\mathbf{y}={\mathbf{X}}^{\top }\mathbf{X}\mathbf{b}$, where ${\mathbf{X}}^{\top }\mathbf{X}$ is the sum-of-squares and cross-products matrix, spurred the development of generalized inverses, numerically stable matrix decompositions, and regularization techniques to tackle rank-deficient and ill-posed problems; consult Gentle [4], pp. 126–154, 217–301 for a comprehensive modern overview. Their interconnection — namely, for the Moore-Penrose pseudoinverse ${\mathbf{A}}^{†}$, which has been computed predominantly through singular value decomposition (SVD) since the 1960s, as introduced by Golub and Kahan [28] (a numerical algorithm for the weighted pseudoinverse was provided in the 1980s by Eldén [29]) and was demonstrated, among others, in the 1970s by Ward [30], p. 34 and, more profoundly, in the 2010s by Barata and Hussein [31], pp. 152–154 to be the limit of a Tikhonov regularization formula ${\mathbf{A}}^{†}={lim}_{\mu \to 0}{\left({\mathbf{A}}^{\top }\mathbf{A}+\mu \mathbf{I}\right)}^{-1}{\mathbf{A}}^{\top }$, where $\mathbf{I}\in {\mathbb{R}}^{n\times n}$ is an identity matrix for $\mathbf{A}\in {\mathbb{R}}^{m\times n}$ and $\mu$ is a parameter (for a detailed presentation of the regularization in its original notation, consult Tikhonov et al. [32], pp. 10–67) — has become a standard perspective in modern numerical analysis and forms the foundation of the CLSP framework, the first step of which is $\widehat{\mathbf{z}}=\mathbf{G}\mathbf{b}$, where $\mathbf{G}$ is either ${\mathbf{A}}^{†}\in {\mathbb{R}}^{n\times m}$ or a constrained reflexive inverse ${\mathbf{G}}^{\{1,2\}}\in {\mathbb{R}}^{n\times m}$.

As documented in the seminal works of Rao and Mitra [33], pp. vii–viii and Ben-Israel and Greville [34], pp. 4–5, 370–374, the first pseudoinverses (original term, now typically reserved for ${\mathbf{A}}^{†}$) or generalized inverses (current general term) emerged in the theory of integral operators — introduced by Fredholm in 1903 and further developed by Hurwitz in 1912 — and, implicitly, in the theory of differential operators by Hilbert in 1904, whose work was subsequently extended by Myller in 1906, Westfall in 1909, Bounitzky in 1909, Elliott in 1928, and Reid in 1931. The cited authors attribute their first application to matrices to Moore in 1920 under the term general reciprocals[35] (though some sources suggest he may have formulated the idea as early as 1906), with independent formulations by Siegel in 1937 and Rao in 1955, and generalizations to singular operators by Tseng in 1933, 1949, and 1956, Murray and von Neumann in 1936, and Atkinson in 1952–1953. A theoretical consolidation occurred with Bjerhammar’s [36] rediscovery of Moore’s formulation, followed by Penrose’s [37,38] introduction of four conditions defining the unique least-squares minimum-norm generalized inverse $\mathbf{G}\equiv {\mathbf{A}}^{†}$, which, for $\mathbf{G}\in {\mathbb{R}}^{n\times m}$, can be expressed in the following way:

$$\left(\mathrm{i}\right)\mathbf{A}\mathbf{G}\mathbf{A}=\mathbf{A}\left(\mathrm{ii}\right)\mathbf{G}\mathbf{A}\mathbf{G}=\mathbf{G}\left(\mathrm{iii}\right){\left(\mathbf{A}\mathbf{G}\right)}^{\top }=\mathbf{A}\mathbf{G}\left(\mathrm{iv}\right){\left(\mathbf{G}\mathbf{A}\right)}^{\top }=\mathbf{G}\mathbf{A}$$

(1)

In econometrics, it led to conditionally unbiased (minimum bias) estimators [39] and a redefinition of the Gauss-Markov theorem [40] in the 1960s, superseding earlier efforts [41].

Further contributions to the theory, calculation, and diverse applications of $\mathbf{G}$ were made in the 1950s–early 1970s by Rao and Mitra [33], Rao and Mitra [42] (a synthesis of the authors’ previous works from 1955–1969), Golub and Kahan [28], Ben-Israel and Greville [34], Greville [43], Greville [44], Greville [45], Greville [46], Cline [47], Cline and Greville [48], Ben-Israel and Cohen [49], and Lewis and Newman [50] (the most cited sources on the topic in the Web of Science database, although the list is not exhaustive; consult Rao and Mitra [33], pp. vii–viii for an extended one), which, among others, introduced the concept of a $\{·\}$-generalized inverse $\mathbf{G}\in {\mathbf{A}}^{\{·\}}$, where ${\mathbf{A}}^{\{·\}}\subseteq {\mathbb{R}}^{n\times m}$ is the space of the $\{·\}$-generalized inverses of $\mathbf{A}\in {\mathbb{R}}^{m\times n}$, — hereafter, $\{i,j\}$-inverses for $i>1$, $\{i,j,k\}$-inverses, Drazin inverses, and higher-order $\mathbf{G}$ are disregarded because of their inapplicability to the CLSP framework — by satisfying a subset or all of the (i)–(iv) conditions in Equation (1), as described in Table 1. Formally, for ${\mathbf{G}}^{\left\{1\right\}}\in {\mathbf{A}}^{\left\{1\right\}}$, ${\mathbf{I}}_n\in {\mathbb{R}}^{n\times n}$, ${\mathbf{I}}_m\in {\mathbb{R}}^{m\times m}$, and arbitrary matrices $\mathbf{U},\mathbf{V}\in {\mathbb{R}}^{n\times m}$, any $\left\{1\right\}$-inverse, including all $\{1,j\}$-inverses and ${\mathbf{A}}^{†}$, can be expressed as $\mathbf{G}={\mathbf{G}}^{\left\{1\right\}}+({\mathbf{I}}_n-{\mathbf{G}}^{\left\{1\right\}}\mathbf{A})\mathbf{U}+\mathbf{V}({\mathbf{I}}_m-\mathbf{A}{\mathbf{G}}^{\left\{1\right\}})$, out of which only $\mathbf{G}\equiv {\mathbf{A}}^{†}$ and a constrained $\mathbf{G}\equiv {\mathbf{G}}^{\{1,2\}}\in {\mathbb{R}}^{n\times m}$ — here, the Bott-Duffin inverse [51], expressed in modern notation as ${\mathbf{A}}_{\mathbf{S}}^{†}={\mathbf{P}}_{\mathbf{S}}^L{\left({\mathbf{P}}_{\mathbf{S}}^L\mathbf{A}{\mathbf{P}}_{\mathbf{S}}^R\right)}^{†}{\mathbf{P}}_{\mathbf{S}}^R$, where ${\mathbf{P}}_{\mathbf{S}}^L\in {\mathbb{R}}^{m\times m}$ and ${\mathbf{P}}_{\mathbf{S}}^R\in {\mathbb{R}}^{n\times n}$ are orthogonal (perpendicular) projectors onto the rows and columns defined by an orthonormal matrix $\mathbf{S}$ — can be uniquely defined (${\mathbf{A}}_{\mathbf{S}}^{†}$ under certain conditions) and qualify as minimum-norm unbiased estimators in the sense of Chipman [39] and Price [40].

Starting from the early 1960s, Cline [47], Rao and Mitra [33], pp. 64–71, Meyer [54], Ben-Israel and Greville [34], pp. 175–200, Hartwig [55], Campbell and Meyer [52], pp. 53–61, Rao and Yanai [56], Tian [57], Rakha [58], Wang et al. [53], pp. 193–210, and Baksalary and Trenkler [59], among others, extended the formulas for ${\mathbf{A}}^{†}$ to partitioned (block) matrices, including the general column-wise, ${\mathbf{A}}_1\in {\mathbb{R}}^{m\times n_1}$, ${\mathbf{A}}_2\in {\mathbb{R}}^{m\times n_2}$, and $n_1+n_2=n$ — and, equivalently, row-wise, ${\mathbf{A}}_1\in {\mathbb{R}}^{m_1\times n}$, ${\mathbf{A}}_2\in {\mathbb{R}}^{m_2\times n}$, $m_1+m_2=m$, — cases (employed in the numerical stability analysis of the CLSP estimator for ${\mathbf{C}}_{\mathrm{canon}}=\left[\mathbf{C}\mathbf{S}\right]$ in the decomposition $\mathbf{A}=\mathbf{B}{\mathbf{C}}_{\mathrm{canon}}+\mathbf{R}$ given $\mathbf{A}\in {\mathbb{R}}^{m\times n}$, ${\mathbf{C}}_{\mathrm{canon}}\in {\mathbb{R}}^{k\times n}$, and $\mathbf{B}=\mathbf{A}{\mathbf{C}}_{\mathrm{canon}}^{†}=\mathbf{A}{\left[\mathbf{C}\mathbf{S}\right]}^{†}\in {\mathbb{R}}^{m\times k}$ ):

$${\left[\begin{array}{cc}{\mathbf{A}}_1& {\mathbf{A}}_2\end{array}\right]}^{†}={\left({\left[\begin{array}{cc}{\mathbf{A}}_1& {\mathbf{A}}_2\end{array}\right]}^{\top }\left[\begin{array}{cc}{\mathbf{A}}_1& {\mathbf{A}}_2\end{array}\right]\right)}^{†}{\left[\begin{array}{cc}{\mathbf{A}}_1& {\mathbf{A}}_2\end{array}\right]}^{\top }={\left[\begin{array}{cc}{\mathbf{A}}_1^{\top }{\mathbf{A}}_1& {\mathbf{A}}_1^{\top }{\mathbf{A}}_2\\ {\mathbf{A}}_2^{\top }{\mathbf{A}}_1& {\mathbf{A}}_2^{\top }{\mathbf{A}}_2\end{array}\right]}^{†}\left[\begin{array}{c}{\mathbf{A}}_1^{\top }\\ {\mathbf{A}}_2^{\top }\end{array}\right]$$

(2)

where $rank\left(\left[{\mathbf{A}}_1{\mathbf{A}}_2\right]\right)\le rank\left({\mathbf{A}}_1\right)+rank\left({\mathbf{A}}_2\right)$, with equality iff $\mathcal{R}\left({\mathbf{A}}_1\right)\cap \mathcal{R}\left({\mathbf{A}}_2\right)=\left\{\mathbf{0}\right\}$, and strict inequality iff $\mathcal{R}\left({\mathbf{A}}_1\right)\cap \mathcal{R}\left({\mathbf{A}}_2\right)\ne \left\{\mathbf{0}\right\}$. For the original definitions and formulations, consult Rao and Mitra [33], pp. 64–66 and Ben-Israel and Greville [34], pp. 175–200.

To sum up, $\mathbf{G}$ (especially ${\mathbf{A}}^{†}$), SVD, and regularization techniques — namely, Lasso (based on the ${\mathsf{\ell }}_1$-norm), Tikhonov regularization (hereafter, referred to as Ridge regression, based on the ${\mathsf{\ell }}_2$-norm), and Elastic Net (a convex combination of ${\mathsf{\ell }}_1$ and ${\mathsf{\ell }}_2$ norms) [4], pp. 477–479 — have become an integral part of estimation (model-fitting), where ${\mathbf{A}}^{†}$ can be defined as a minimum-norm best linear unbiased estimator (MNBLUE), reducing to the classical BLUE in (over)determined cases. These methods have been applied in (a) natural sciences since the 1940s [32] (pp. 68–156) [36,51] [60], pp. 85–172 [61,62,63], and (b) econometrics since the 1960s [39,40] [64], pp. 36–157, 174–232 [65], pp. 61–106 [66], pp. 37–338 in both unconstrained and equality-constrained forms, thereby relaxing previous rank restrictions.

The incorporation of inequality constraints into optimization problems during the 1930s–1950s (driven by the foundational contributions of Kantorovich [67], Koopmans [7], and Dantzig [6], along with the later authors [5]) marked the next milestone in the development of convex optimization under the emerging framework of programming, namely for linear cases (LP) — formally, ${\mathbf{x}}^{*}=arg{max}_{\mathbf{x}\in {\mathbb{R}}^n}{\mathbf{c}}^{\top }\mathbf{x}$ subject to $\mathbf{A}\mathbf{x}\le \mathbf{b}$, $\mathbf{x}\ge \mathbf{0}$ (with the dual ${\lambda }^{*}=arg{min}_{\lambda \in {\mathbb{R}}^m}{\mathbf{b}}^{\top }\lambda$ subject to ${\mathbf{A}}^{\top }\lambda \ge \mathbf{c})$, where $\mathbf{A}\in {\mathbb{R}}^{m\times n}$, $\mathbf{b}\in {\mathbb{R}}^m$, and $\mathbf{c}\in {\mathbb{R}}^n$ — and convex quadratic ones (QP) — formally, ${\mathbf{x}}^{*}=arg{min}_{\mathbf{x}\in {\mathbb{R}}^n}\left\{\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Q}\mathbf{x}+{\mathbf{c}}^{\top }\mathbf{x}\right\}$ subject to $\mathbf{A}\mathbf{x}\le \mathbf{b}$ (with the dual ${\lambda }^{*}=arg{max}_{\lambda \in {\mathbb{R}}_{+}^m}\left\{-\frac{1}{2}{({\mathbf{A}}^{\top }\lambda +\mathbf{c})}^{\top }{\mathbf{Q}}^{-1}({\mathbf{A}}^{\top }\lambda +\mathbf{c})-{\mathbf{b}}^{\top }\lambda \right\}$), where $\mathbf{Q}\in {\mathbb{R}}^{n\times n}$ is symmetric positive definite, $\mathbf{c}\in {\mathbb{R}}^n$, $\mathbf{A}\in {\mathbb{R}}^{m\times n}$, and $\mathbf{b}\in {\mathbb{R}}^m$ — [2] (pp. 146–213) (for a classification of programming, consult Figure 1) [6] (pp. 1–31) [1], pp. 355–597, with Dantzig’s simplex and Karmarkar’s interior-point algorithms being efficient solvers [68,69].

In the late 1950s–1970s, the application of LP and QP extended to LS problems (which can be referred to as a major generalization of mainstream convex optimization methods; consult Boyd and Vandenberghe [2], pp. 291–349 for a comprehensive modern overview) with the pioneering works of Wagner [70] and Sielken and Hartley [71], expanded by Kiountouzis [72] and Sposito [73], who first substituted LS with LP-based (${\mathsf{\ell }}_1$-norm) least absolute (LAD) and (${\mathsf{\ell }}_{\infty }$-norm) least maximum (LDP) deviations and derived unbiased (${\mathsf{\ell }}_p$-norm, $1<p<\infty$) estimators with non-unique solutions; whilst, in the early 1990s, Stone and Tovey [69] demonstrated algorithmic equivalence between LP algorithms and iteratively reweighted LS. Judge and Takayama [74] and Mantel [75] reformulated (multiple) LS with inequality constraints as QP, and Lawson and Hanson [64], pp. 144–147 introduced, among others, the non-negative least squares (NNLS) method. These and further developments are reflected in the second step of the CLSP framework, where $\widehat{\mathbf{z}}=\mathbf{G}\mathbf{b}$ is corrected under a regularization scheme inspired by Lasso ($\alpha =0$), Ridge ($\alpha =1$), and Elastic Net ($0<\alpha <1$) regularization, where ${\widehat{\mathbf{z}}}^{*}=arg{min}_{\mathbf{z}\in {\mathbb{R}}^n}\left\{(1-\alpha ){∥\mathbf{z}-\widehat{\mathbf{z}}∥}_1+\alpha {∥\mathbf{z}-\widehat{\mathbf{z}}∥}_2^2\right\}$, subject to $\mathbf{A}\mathbf{z}=\mathbf{b}$.

The final frontier that motivates the formulation of a new framework is the class of ill-posed problems in Hadamard’s sense [4] (pp. 143, 241), i.e., systems with no (i) feasibility, $\mathbf{x}\notin \mathcal{F}:=\left\{\mathbf{x}\right|{\forall }_{i,j}{g}_i\left(\mathbf{x}\right)\le 0,h_j\left(\mathbf{x}\right)=0\}$, (ii) uniqueness, $\exists {\mathbf{x}}_1\ne {\mathbf{x}}_2:\mathbf{A}{\mathbf{x}}_1=\mathbf{A}{\mathbf{x}}_2$, (iii) stability, $\forall \epsilon >0,\exists \delta >0:∥\mathbf{b}-{\mathbf{b}}_{\delta }∥<\delta \Rightarrow ∥{\mathbf{A}}^{†}\mathbf{b}-{\mathbf{A}}^{†}{\mathbf{b}}_{\delta }∥\ge \epsilon$, or formulation — as a well-posed LS, LP, or QP problem — of solution under mainstream assumptions, where $\mathbf{A}\in {\mathbb{R}}^{m\times n}$, $\mathbf{x}\in {\mathbb{R}}^n$, and $\mathbf{b}\in {\mathbb{R}}^m$, such as the ones described in Bolotov [15], Bolotov [16]: namely, (a) underdetermined, $rank\left(\mathbf{A}\right)<min(m,n)$, and/or ill-conditioned linear systems, $\kappa \left(\mathbf{A}\right)={\parallel \mathbf{A}\parallel }_2{\parallel {\mathbf{A}}^{†}\parallel }_2\gg 1$, $\mathbf{A}\in {\mathbb{R}}^{m\times n}$, in cases of LS and LP, where solutions are either non-unique or highly sensitive to perturbations [76,77]; and (b) degenerate, $\#\left\{\mathrm{active}\mathrm{constraints}\mathrm{at}\mathbf{x}\right\}>dim\left(\mathbf{x}\right)=n$, with cycling behavior, $c^{\top }{\mathbf{x}}^{(k+1)}=c^{\top }{\mathbf{x}}^{\left(k\right)}$, and ${\mathbf{x}}^{(k+1)}\ne {\mathbf{x}}^{\left(k\right)}$, $\mathbf{x}\in {\mathbb{R}}^n$, in all LP and QP problems, where efficient methods, such as the simplex algorithm, may fail to converge [78,79]. Such problems have been addressed since the 1940s–1970s by (a) regularization (i.e., a method of perturbations) and (b) problem-specific algorithms [2] (pp. 455–630). For LS cases, the list includes, among others, Lasso, Ridge regularization, and Elastic Net [4] (pp. 477–479), constrained LS methods [80], restricted LS [81], LP-based sparse recovery approaches [77], Gröbner bases [82], as well as derivatives of eigenvectors and Nelson-type, BD, QR- and SVD-based algorithms [83]. For LS and QP — among others, Wolfe’s regularization-based technique for simplex [84] and its modifications [78,79], Dax’s LS-based steepest-descent algorithm [85], a primal-dual NNLS-based algorithm (LSPD) [86], as well as modifications of the ellipsoid algorithm [87]. A major early attempt to unify (NN)LS, LP, and QP within a single framework to solve both well- and ill-posed convex optimization problems was undertaken by Übi [88], Übi [89], Übi [90], Übi [91], Übi [92], who proposed a series of problem-specific algorithms (VD, INE, QPLS, and LS1). Figure 2 compares all the mentioned methods by functionality.

Therefore, to accommodate all considered problem classes, a modular two-step CLSP consists of (first, compulsory, step) minimum-norm LS estimation, based on ${\mathbf{A}}^{†}$ or ${\mathbf{A}}_{\mathbf{S}}^{†}$, followed by a (second, optional, step) regularization-based correction, contributing to the generalization of convex optimization methods by eliminating the need for problem-specific algorithms.

3. Construction and Formalization of the CLSP Estimator

Based on the citation and reference counts, as well as in-group relevance (i.e., the number of citations within the examined collection) of the selected sources (consult Figure 3), this text incorporates the seminal studies on the topic, with the exception of historical seminar proceedings (e.g., Nashed) and a few secondary contributions by Ben-Israel, Sergienko, and Grafarend (Lawson and Hanson, erroneously indicated as missing due to incomplete monograph indexing, are already included in the bibliography as a 1995 reprint of the 1976 edition [64]). Hence, it can be claimed that the CLSP framework is grounded in the general algorithmic logic of Wolfe [84] and Dax [85], as further formalized by Osborne [78], who were among the first to introduce an additional estimation step into simplex-based solvers.

The proposed CLSP estimator is modular — i.e., it comprises selectively omittable and repeatable actions for enhanced flexibility — and consists of two steps: (a) the first, grounded in the theory of (constrained) pseudoinverses as summarized in Rao and Mitra [33], Ben-Israel and Greville [34], Lawson and Hanson [64], and Wang et al. [53], is iterable to refine the solution and serves as a mandatory baseline in cases where the second step becomes infeasible (due to mutually inconsistent constraints) — as exemplified by Whiteside et al. [93] for regression-derived constraints and by Blair [76], Blair [94] in systems with either too few or too many effective constraints; — and (b) the second, grounded in the theory of regularization and convex programming as summarized in Tikhonov et al. [32], Gentle [4], Nocedal and Wright [1], and Boyd and Vandenberghe [2], provides an optional correction of the first step, drawing on the logic of Lasso, Ridge, and Elastic Net to address ill-posedness and compensate for constraint violations or residual approximation errors resulting from the first-step minimum-norm LS estimate. This structure ensures that CLSP is able to yield a unique solution under a strictly convex second-step correction and extends the classical LS framework in both scope and precision. The modular algorithmic flow is illustrated in Figure 4.

The first process block in the CLSP algorithmic flow denotes the transformation of the initial problem — ${\widehat{\mathbf{x}}}_M^{*}={Proj}_{{\mathbf{x}}_M}\left\{{\mathbf{Q}}_X\left[\begin{array}{c}{\mathbf{x}}_M\\ {\mathbf{x}}_L\end{array}\right]\in {\mathbb{R}}^p:{\forall }_{1\le i\le s}{g}_i({\mathbf{x}}_M,{\mathbf{x}}_L)⋛{\gamma }_i,{\forall }_j{h}_j({\mathbf{x}}_M,{\mathbf{x}}_L)={\eta }_j\right\}$, where ${Proj}_{{\mathbf{x}}_M}$ is the projection onto the ${\mathbf{x}}_M$ coordinates, ${\mathbf{x}}_M\in {\mathbb{R}}^{p-l}$ is the vector of model (target) variables to be estimated (${\widehat{\mathbf{x}}}_M^{*}$), ${\mathbf{x}}_L\in {\mathbb{R}}^l$, $l\ge 0$, is a vector of latent variables that appear only in the constraint functions $g_i(·)$ and $h_j(·)$, such that ${\mathbf{Q}}_X\left[\begin{array}{c}{\mathbf{x}}_M\\ {\mathbf{x}}_L\end{array}\right]=\mathbf{x}\in {\mathbb{R}}^p$, ${\mathbf{Q}}_X$ is the permutation matrix, and $g_i(·)⋛{\gamma }_i$ and $h_j(·)={\eta }_j$ are the linear(ized) inequality and equality constraints, ${\forall }_{i,j}{\gamma }_i,{\eta }_j\in \mathbf{b}$ — to the canonical form ${\mathbf{A}}^{\left(r\right)}{\mathbf{z}}^{\left(r\right)}=\mathbf{b}$ (a term first introduced in LP [6] (pp. 75–81)), where (a) ${\mathbf{A}}^{\left(r\right)}\in {\mathbb{R}}^{m\times n}$ is the block design matrix consisting of a ${\forall }_i{g}_i(·)$-and-${\forall }_j{h}_j(·)$ constraints matrix $\mathbf{C}\in {\mathbb{R}}^{k\times p}$, a model matrix$\mathbf{M}={\mathbf{Q}}_M^L\left[\begin{array}{c}{\mathbf{M}}_{\mathrm{partial}}\\ \mathbf{0}\end{array}\right]{\mathbf{Q}}_M^R\in {\mathbb{R}}^{(m-k)\times p}$, a sign slack matrix $\mathbf{S}={\mathbf{Q}}_S^L\left[\begin{array}{c}{\mathbf{S}}_s\subseteq {\forall }_{{\sigma }_1,\dots ,{\sigma }_s\in \{-1,+1\}}diag\left({\sigma }_1,\dots ,{\sigma }_s\right)\\ \mathbf{0}\end{array}\right]{\mathbf{Q}}_S^R\in {\mathbb{R}}^{k\times (n-p)}$, and either a zero matrix $\mathbf{0}$ (if $r=1$) or a reverse-sign — i.e., based on $-sign(·)$ — slack matrix ${\mathbf{S}}_{\mathrm{residual}}=\left[\mathbf{0}diag\left(-sign\left({(\mathbf{b}-{\mathbf{A}}^{(r-1)}{\widehat{\mathbf{z}}}^{(r-1)})}_{k+1:m}\right)\right)\right]\in {\mathbb{R}}^{(m-k)\times (n-p)}$ (provided $r>1$):

$${\mathbf{A}}^{\left(r\right)}=\left[\begin{array}{cc}\mathbf{C}& \mathbf{S}\\ \mathbf{M}& \left\{\begin{array}{cc}\mathbf{0},\hfill & \mathrm{initial}\mathrm{iteration},r=1\hfill \\ {\mathbf{S}}_{\mathrm{residual}},\hfill & \mathrm{subsequent}\mathrm{iterations},r>1\hfill \end{array}\right.\end{array}\right]$$

(3)

and (b) ${\mathbf{z}}^{\left(r\right)}=\left[\begin{array}{c}{\mathbf{x}}^{\left(r\right)}\\ {\left({\mathbf{y}}^{\left(r\right)}\right)}^{*}\end{array}\right]\in {\mathbb{R}}^n$ is the full vector of unknowns, comprising (model) ${\mathbf{x}}^{\left(r\right)}\in {\mathbb{R}}^p$ and (slack) ${\left({\mathbf{y}}^{\left(r\right)}\right)}^{*}=\left[\begin{array}{c}{\mathbf{y}}^{\left(r\right)}\\ {\mathbf{y}}_{\mathrm{residual}}^{\left(r\right)}\end{array}\right]\in {\mathbb{R}}^{n-p}$, $y_i^{\left(r\right)}\ge 0$ and $y_j^{\left(r\right)}=0$ — an equivalent problem [2] (pp. 130–132), where $g_i{(·)}^a=g_i(·)+{\mathbf{S}}_{ii}{y}_i^{\left(r\right)}=0$ and $h_j{(·)}^a=h_j(·)$. The constraint functions $g_i{(·)}^a$ and $h_j{(·)}^a$ must be linear(ized), so that $\mathbf{C},\mathbf{S}\subseteq {\mathbf{A}}^{\left(r\right)}\iff \mathbf{C}{\mathbf{x}}^{\left(r\right)}+\mathbf{S}{\left({\mathbf{y}}^{\left(r\right)}\right)}^{*}={\mathbf{b}}_{1:k}$, $k\le m$.

In practical implementations, for problems with $(m>m_E)\vee (n>n_E)$, where E is an estimability limit, the partitioning ${\mathbf{A}}^{\left(r\right)}=\left\{{\mathbf{A}}_{subset,i_E,j_E}^{\left(r\right)}\right\}$, $m_{\mathrm{subset}}=m_{\mathrm{reduced}}-1$, $n_{\mathrm{subset}}=n_{\mathrm{reduced}}-1$, $i_E\in \{1,\dots ,⌈{\displaystyle \frac{m}{m_{\mathrm{subset}}}}⌉\}$, and $j_E\in \{1,\dots ,⌈{\displaystyle \frac{n}{n_{\mathrm{subset}}}}⌉\}$, can be used, leading to a total of $\left(⌈{\displaystyle \frac{m}{m_{\mathrm{subset}}}}⌉·⌈{\displaystyle \frac{n}{n_{\mathrm{subset}}}}⌉\right)$${\mathbf{A}}_{\mathrm{reduced}}^{\left(r\right)}\in {\mathbb{R}}^{m_{\mathrm{reduced}}\times n_{\mathrm{reduced}}}$, composed of ${\mathbf{A}}_{\mathrm{subset}}^{\left(r\right)}$ and $\sum a_{i,j}^{\left(r\right)}$, $i\notin {\mathcal{I}}_{\mathrm{subset}}\subset \{1,\dots ,m\}$, $j\notin {\mathcal{J}}_{\mathrm{subset}}\subset \{1,\dots ,n\}$, with a total of $\left(⌈{\displaystyle \frac{m}{m_{\mathrm{subset}}}}⌉\right)$${\mathbf{b}}_{\mathrm{r}.}\in {\mathbb{R}}^{m_{\mathrm{reduced}}}$:

$${\mathbf{A}}_{\mathrm{reduced}}^{\left(r\right)}={\mathbf{Q}}_A^{\mathcal{I}}\left[\begin{array}{cc}{\mathbf{A}}_{{\mathcal{I}}_{\mathrm{subset}},{\mathcal{J}}_{\mathrm{subset}}}^{\left(r\right)}& \sum _{j\notin {\mathcal{J}}_{\mathrm{subset}}}{a}_{i,j}^{\left(r\right)}\\ \sum _{i\notin {\mathcal{I}}_{\mathrm{subset}}}{a}_{i,j}^{\left(r\right)}& \sum _{i\notin {\mathcal{I}}_{\mathrm{subset}}}\sum _{j\notin {\mathcal{J}}_{\mathrm{subset}}}{a}_{i,j}^{\left(r\right)}\end{array}\right]{\mathbf{Q}}_A^{\mathcal{J}},{\mathbf{b}}_{\mathrm{r}.}={\mathbf{Q}}_b^{\mathcal{I}}\left[\begin{array}{c}{\mathbf{b}}_{{\mathcal{I}}_{\mathrm{subset}}}\\ \sum _{i\notin {\mathcal{I}}_{\mathrm{subset}}}{b}_i\end{array}\right]$$

(4)

$\sum a_{i,j}^{\left(r\right)}$ are then treated as slack rows/columns, i.e., incorporated into $\mathbf{S}$, to act as balancing adjustments since the reduced model, by construction, tends to inflate the estimates ${\widehat{a}}_{i_E,j_E}$.

The second process block, i.e., the first step of the CLSP estimator, denotes obtaining an (iterated if $r>1$) estimate ${\widehat{\mathbf{z}}}^{\left(r\right)}={\mathbf{G}}^{\left(r\right)}\mathbf{b}$ (alternatively, $\left(⌈{\displaystyle \frac{m}{m_{\mathrm{subset}}}}⌉·⌈{\displaystyle \frac{n}{n_{\mathrm{subset}}}}⌉\right)$-times ${\widehat{\mathbf{z}}}_{\mathrm{reduced}}^{\left(r\right)}={\mathbf{G}}_{\mathrm{reduced}}^{\left(r\right)}{\mathbf{b}}_{\mathrm{r}.}$ with ${\widehat{\mathbf{z}}}_{i_E,j_E}^{\left(r\right)}=\left\{{\widehat{\mathbf{z}}}_{\mathrm{reduced},1:m_{\mathrm{subset}},1:n_{\mathrm{subset}}}^{\left(r\right),(i_E,j_E)}\right\}$) through the (constrained) Bott-Duffin inverse ${\mathbf{G}}^{\left(r\right)}\equiv {\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}\in {\mathbb{R}}^{n\times m}$[51][53](pp. 49–64), defined on a subspace specified by a symmetric idempotent matrix $\mathbf{Z}\in {\mathbb{R}}^{n\times n}$ (regulating solution-space exclusion through binary entries 0/1, such that ${\mathbf{P}}_{\mathbf{Z}}^L=\mathbf{Z}$ restricts the domain of estimated variables to the selected subspace, while ${\mathbf{P}}_{\mathbf{Z}}^R={\mathbf{I}}_m$ leaves the data vector, i.e., the input $\mathbf{b}$, unprojected). ${\mathbf{G}}^{\left(r\right)}\equiv {\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}$ reduces to the Moore-Penrose pseudoinverse ${\mathbf{G}}^{\left(r\right)}\equiv {\left({\mathbf{A}}^{\left(r\right)}\right)}^{†}\in {\mathbb{R}}^{n\times m}$ when $\mathbf{Z}={\mathbf{I}}_n$, and equals the standard inverse ${\mathbf{G}}^{\left(r\right)}\equiv {\left({\mathbf{A}}^{\left(r\right)}\right)}^{-1}$ iff $\mathbf{Z}={\mathbf{I}}_n\wedge rank\left({\mathbf{A}}^{\left(r\right)}\right)=m=n$:

$${\widehat{\mathbf{z}}}^{\left(r\right)}={\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}\mathbf{b}=\left[{\left(\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}\mathbf{Z}\right)}^{†}\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }\right]\mathbf{b},{\mathbf{Z}}^2=\mathbf{Z}={\mathbf{Z}}^{\top }$$

(5)

The number of iterations (r) increases until a termination condition is met, such as error $\left|{\displaystyle \frac{1}{\sqrt{m}}}{∥(\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{\left(r\right)})∥}_2-{\displaystyle \frac{1}{\sqrt{m}}}{∥(\mathbf{b}-{\mathbf{A}}^{(r-1)}{\widehat{\mathbf{z}}}^{(r-1)})∥}_2\right|<\epsilon$ or $r>r_{\mathrm{limit}}$, whichever occurs first. In practical implementations, ${\widehat{\mathbf{z}}}^{\left(r\right)}$ is efficiently computed using the SVD (see Appendix A.1).

The third process block, i.e., the second step of the CLSP estimator, denotes obtaining a regularized solution ${\widehat{\mathbf{z}}}^{*}$ (alternatively, $\left(⌈{\displaystyle \frac{m}{m_{\mathrm{subset}}}}⌉·⌈{\displaystyle \frac{n}{n_{\mathrm{subset}}}}⌉\right)$-times ${\widehat{\mathbf{z}}}_{\mathrm{reduced}}^{*}$ with ${\widehat{\mathbf{z}}}_{i_E,j_E}^{*}=\left\{{\widehat{\mathbf{z}}}_{\mathrm{reduced},1:m_{\mathrm{subset}},1:n_{\mathrm{subset}}}^{*,(i_E,j_E)}\right\}$) — unless the algorithmic flow terminates after Step 1, in which case ${\widehat{\mathbf{z}}}^{*}={\widehat{\mathbf{z}}}^{\left(r\right)}$ — by penalizing deviations from the (iterated if $r>1$) estimate ${\widehat{\mathbf{z}}}^{\left(r\right)}$ reflecting Lasso ($\alpha =0\iff {\widehat{\mathbf{z}}}^{*}={\widehat{\mathbf{z}}}^{\left(r\right)}+\widehat{\mathbf{u}}$, where $\widehat{\mathbf{u}}=arg{min}_{\mathbf{u}}{\parallel \mathbf{u}\parallel }_1,\mathrm{s}.\mathrm{t}.{\mathbf{A}}^{\left(r\right)}({\widehat{\mathbf{z}}}^{\left(r\right)}+\mathbf{u})=\mathbf{b}$), Ridge ($\alpha =1\iff {\widehat{\mathbf{z}}}^{*}={\widehat{\mathbf{z}}}^{\left(r\right)}+{\left({\mathbf{A}}^{\left(r\right)}\right)}^{†}{\mathbf{A}}^{\left(r\right)}\widehat{\mathbf{u}}$, where ${\left({\mathbf{A}}^{\left(r\right)}\right)}^{†}{\mathbf{A}}^{\left(r\right)}={\mathbf{P}}_{\mathcal{R}\left({\mathbf{A}}^{\top }\right)}$, a projector equal to $\mathbf{I}$ iff ${\mathbf{A}}^{\left(r\right)}$ is of full column rank, and $\widehat{\mathbf{u}}=arg{min}_{\mathbf{u}}{\parallel \mathbf{u}\parallel }_2^2,\mathrm{s}.\mathrm{t}.{\mathbf{A}}^{\left(r\right)}({\widehat{\mathbf{z}}}^{\left(r\right)}+\mathbf{u})=\mathbf{b}$), and Elastic Net ($0<\alpha <1\iff {\widehat{\mathbf{z}}}^{*}={\widehat{\mathbf{z}}}^{\left(r\right)}+(1-\alpha )\widehat{\mathbf{u}}+\alpha {\left({\mathbf{A}}^{\left(r\right)}\right)}^{†}{\mathbf{A}}^{\left(r\right)}\widehat{\mathbf{u}}$, where $\widehat{\mathbf{u}}=arg{min}_{\mathbf{u}}\left\{{(1-\alpha )\parallel \mathbf{u}\parallel }_1+\alpha {\parallel \mathbf{u}\parallel }_2^2\right\},\mathrm{s}.\mathrm{t}.{\mathbf{A}}^{\left(r\right)}({\widehat{\mathbf{z}}}^{\left(r\right)}+\mathbf{u})=\mathbf{b}$) [4] (pp. 477–479), [2] (pp. 306–310):

$${\widehat{\mathbf{z}}}^{*}=arg\underset{\mathbf{z}\in {\mathbb{R}}^n}{min}\left\{(1-\alpha ){∥\mathbf{z}-{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_1+\alpha {∥\mathbf{z}-{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_2^2\right\},\alpha \in [0,1],\mathrm{s}.\mathrm{t}.{\mathbf{A}}^{\left(r\right)}\mathbf{z}=\mathbf{b}$$

(6)

where the parameter $\alpha$ may be selected arbitrarily or determined from the input $\mathbf{b}$ using an appropriate criterion, such as minimizing prediction error via cross-validation, satisfying model-selection heuristics (e.g., AIC or BIC), or optimizing residual metrics under known structural constraints. Alternatively, $\alpha$ may be fixed at benchmark values, e.g., $\alpha \in \{0,{\displaystyle \frac{1}{2}},1\}$, or based on an error rule, e.g., $\alpha =min\left({\displaystyle \frac{{\left(\frac{1}{\sqrt{m}}{∥(\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{*})∥}_2\right)}_{\alpha =0}}{{\left(\frac{1}{\sqrt{m}}{∥(\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{*})∥}_2\right)}_{\alpha =0}+{\left(\frac{1}{\sqrt{m}}{∥(\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{*})∥}_2\right)}_{\alpha =1}+\epsilon }},1\right)$, $\epsilon \to 0^{+}$. The final value ${\widehat{\mathbf{z}}}^{*}$ is obtained by solving a convex programming problem, which, in practical implementations, ${\widehat{\mathbf{z}}}^{*}$ is efficiently computed with the help of numerical solvers.

This formalization of the process blocks ensures that the CLSP estimator yields a unique final solution ${\widehat{\mathbf{z}}}^{*}$ under regularization weight $\alpha \in (0,1]$ owing to the strict convexity of the objective function (consult Theorem 1 for a proof). For $\alpha =0$, the solution remains convex but may not be unique unless additional structural assumptions are imposed. In the specific case of $\alpha =1$, the estimator reduces to the Minimum-Norm Best Linear Unbiased Estimator (MNBLUE), or, alternatively, to a block-wise $m_{\mathrm{subset}}\times n_{\mathrm{subset}}$-MNBLUE across $\left(⌈{\displaystyle \frac{m}{m_{\mathrm{subset}}}}⌉·⌈{\displaystyle \frac{n}{n_{\mathrm{subset}}}}⌉\right)$ reduced problems (for a proof, consult Theorem 2 and its corollaries).

Theorem 1.

Let ${\widehat{\mathbf{z}}}^{\left(r\right)}={\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}\mathbf{b}$ be the r-iterated estimate obtained in the first step of the CLSP algorithm, and let ${\widehat{\mathbf{z}}}^{*}$ be the second-step correction obtained through convex programming by solving ${\widehat{\mathbf{z}}}^{*}=arg{min}_{\mathbf{z}\in {\mathbb{R}}^n}\left\{(1-\alpha ){∥\mathbf{z}-{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_1+\alpha {∥\mathbf{z}-{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_2^2\right\},\alpha \in (0,1],\mathrm{s}.\mathrm{t}.{\mathbf{A}}^{\left(r\right)}\mathbf{z}=\mathbf{b}$ (i.e., the regularization parameter excludes a pure Lasso-type correction), then the final solution ${\widehat{\mathbf{z}}}^{*}$ is unique.

Proof. 

Let ${\widehat{\mathbf{z}}}^{\left(r\right)}={\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}\mathbf{b}$ denote the linear estimator obtained via the Bott-Duffin inverse, defined on a subspace determined by a symmetric idempotent matrix $\mathbf{Z}$, and producing a conditionally unbiased estimate over that subspace. The Bott-Duffin inverse ${\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}$ is given by ${\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}={\left(\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}\mathbf{Z}\right)}^{†}\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }$, and the estimate ${\widehat{\mathbf{z}}}^{\left(r\right)}$ is unique if $rank\left({\mathbf{A}}^{\left(r\right)}\mathbf{Z}\right)=dim\left(\mathcal{R}\left(\mathbf{Z}\right)\right)$. In this case, ${\widehat{\mathbf{z}}}^{\left(r\right)}$ is the unique minimum-norm solution in $\mathcal{R}\left(\mathbf{Z}\right)$. In all other cases, the solution set is affine and given by ${\widehat{\mathbf{z}}}^{\left(r\right)}={\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}\mathbf{b}+\left[ker\left({\mathbf{A}}^{\left(r\right)}\mathbf{Z}\right)\cap \mathcal{R}\left(\mathbf{Z}\right)\right]$, where the null-space component represents degrees of freedom not fixed by the constraint ${\mathbf{A}}^{\left(r\right)}\mathbf{Z}$, hence ${\widehat{\mathbf{z}}}^{\left(r\right)}$ is not unique. At the same time, the second-step estimate ${\widehat{\mathbf{z}}}^{*}$ is obtained by minimizing the function $(1-\alpha ){∥\mathbf{z}-{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_1+\alpha {∥\mathbf{z}-{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_2^2$ over the affine (hence convex and closed) constraint set ${\mathcal{F}}^{\left(r\right)}=\left\{\mathbf{z}\in {\mathbb{R}}^n|{\mathbf{A}}^{\left(r\right)}\mathbf{z}=\mathbf{b}\right\}$. Under $\alpha \in (0,1]$, the quadratic term $\alpha {∥\mathbf{z}-{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_2^2$ contributes a positive-definite Hessian $2\alpha \mathbf{I}\succ 0$ to the objective function, making it a strictly convex parabola-like curvature, and, given ${\mathcal{F}}^{\left(r\right)}\ne \varnothing$, the minimizer ${\mathbf{z}}^{*}$ exists and is unique. Therefore, the CLSP estimator with $\alpha \in (0,1]$ yields a unique ${\widehat{\mathbf{z}}}^{*}$.    □

Corollary 1.

Let the final solution be ${\widehat{\mathbf{z}}}^{*}={\widehat{\mathbf{z}}}^{\left(r\right)}={\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}\mathbf{b}$, where ${\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}$ denotes the Bott-Duffin inverse on a subspace defined by a symmetric idempotent matrix $\mathbf{Z}$. Then (one-step) ${\widehat{\mathbf{z}}}^{*}$ is unique iff $rank\left({\mathbf{A}}^{\left(r\right)}\mathbf{Z}\right)=dim\left(\mathcal{R}\left(\mathbf{Z}\right)\right)$. Else, the solution set is affine and given by ${\widehat{\mathbf{z}}}^{*}\in {\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}\mathbf{b}+\left[ker\left({\mathbf{A}}^{\left(r\right)}\mathbf{Z}\right)\cap \mathcal{R}\left(\mathbf{Z}\right)\right]$, which implies that the minimizer ${\mathbf{z}}^{*}$ is not unique, i.e., $\#argmin(·)>1$.

Corollary 2.

Let ${\widehat{\mathbf{z}}}^{*}$ be the final solution obtained in the second step of the CLSP estimator by solving ${\widehat{\mathbf{z}}}^{*}=arg{min}_{\mathbf{z}\in {\mathbb{R}}^n}\left\{{∥\mathbf{z}-{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_1\right\},\mathrm{s}.\mathrm{t}.{\mathbf{A}}^{\left(r\right)}\mathbf{z}=\mathbf{b}$ (i.e., the regularization parameter is set to $\alpha =0$, corresponding to a Lasso correction). Then the solution ${\widehat{\mathbf{z}}}^{*}$ exists and the problem is convex. The solution is unique iff the following two conditions are simultaneously satisfied: (1) the affine constraint set ${\mathcal{F}}^{\left(r\right)}=\left\{\mathbf{z}\in {\mathbb{R}}^n|{\mathbf{A}}^{\left(r\right)}\mathbf{z}=\mathbf{b}\right\}$ intersects the subdifferential of $\parallel \mathbf{z}-{\widehat{\mathbf{z}}}^{\left(r\right)}{\parallel }_1$ at exactly one point, and (2) the objective function is strictly convex on ${\mathcal{F}}^{\left(r\right)}$, which holds when ${\mathcal{F}}^{\left(r\right)}$ intersects the interior of at most one orthant in ${\mathbb{R}}^n$. In all other cases, the minimizer ${\mathbf{z}}^{*}$ is not unique, and the set of optimal solutions forms a convex subset of the feasible region ${\mathcal{F}}^{\left(r\right)}$; that is, the final solution ${\widehat{\mathbf{z}}}^{*}\in arg{min}_{\mathbf{z}\in {\mathcal{F}}^{\left(r\right)}}{∥\mathbf{z}-{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_1$, where $argmin(·)$ is convex and $\#argmin(·)>1$.

Theorem 2.

Let ${\widehat{\mathbf{z}}}^{\left(r\right)}={\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}\mathbf{b}$ be the r-iterated estimate obtained in the first step of the CLSP algorithm, and let ${\widehat{\mathbf{z}}}^{*}$ be the second-step correction obtained through convex programming by solving ${\widehat{\mathbf{z}}}^{*}=arg{min}_{\mathbf{z}\in {\mathbb{R}}^n}\left\{{∥\mathbf{z}-{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_2^2\right\},\mathrm{s}.\mathrm{t}.{\mathbf{A}}^{\left(r\right)}\mathbf{z}=\mathbf{b}$ (i.e., the regularization parameter is $\alpha =1$, corresponding to a Ridge correction), then ${\widehat{\mathbf{z}}}^{*}$ is the Minimum-Norm Best Linear Unbiased Estimator (MNBLUE) of $\mathbf{z}$ under the linear(ized) constraints set by the canonical form of the CLSP.

Proof. 

Let ${\widehat{\mathbf{z}}}^{\left(r\right)}={\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}\mathbf{b}$ denote the linear estimator obtained via the Bott-Duffin inverse, defined on a subspace determined by a symmetric idempotent matrix $\mathbf{Z}$, and producing a conditionally unbiased estimate over that subspace. The Bott-Duffin inverse is ${\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}={\left(\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}\mathbf{Z}\right)}^{†}\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }$ and the estimate ${\widehat{\mathbf{z}}}^{\left(r\right)}$ is unique if $rank\left({\mathbf{A}}^{\left(r\right)}\mathbf{Z}\right)=dim\left(\mathcal{R}\left(\mathbf{Z}\right)\right)$. Given the linear model derived from the canonical form $\mathbf{b}={\mathbf{A}}^{\left(r\right)}{\mathbf{z}}^{\left(r\right)}+\epsilon$, where $\epsilon$ are residuals with $\mathbb{E}\left[\epsilon \right]=\mathbf{0}$, it follows that $\mathbb{E}\left[{\widehat{\mathbf{z}}}^{\left(r\right)}\right]=\mathbb{E}\left[{\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}\mathbf{b}\right]={\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}\mathbb{E}\left[\mathbf{b}\right]={\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}{\mathbf{A}}^{\left(r\right)}{\mathbf{z}}^{\left(r\right)}$. Substituting ${\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}$ yields $\mathbb{E}\left[{\widehat{\mathbf{z}}}^{\left(r\right)}\right]={\left(\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}\mathbf{Z}\right)}^{†}\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}{\mathbf{z}}^{\left(r\right)}$. By the generalized projection identity ${\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}{\mathbf{A}}^{\left(r\right)}\mathbf{Z}=\mathbf{Z}$, one obtains $\mathbb{E}\left[{\widehat{\mathbf{z}}}^{\left(r\right)}\right]=\mathbf{Z}{\mathbf{z}}^{\left(r\right)}$, which proves conditional unbiasedness on $\mathcal{R}\left(\mathbf{Z}\right)$ (and full unbiasedness if $\mathbf{Z}=\mathbf{I}$). Subsequently, the second-step estimate ${\widehat{\mathbf{z}}}^{*}$ is obtained by minimizing the squared Euclidean distance between ${\widehat{\mathbf{z}}}^{\left(r\right)}$ and $\mathbf{z}$, subject to the affine constraint ${\mathbf{A}}^{\left(r\right)}\mathbf{z}=\mathbf{b}$, which corresponds, in geometric and algebraic terms, to an orthogonal projection of ${\widehat{\mathbf{z}}}^{\left(r\right)}$ onto the affine (hence convex and closed) subspace ${\mathcal{F}}^{\left(r\right)}=\left\{\mathbf{z}\in {\mathbb{R}}^n|{\mathbf{A}}^{\left(r\right)}\mathbf{z}=\mathbf{b}\right\}$. By Theorem 1, the unique minimizer ${\mathbf{z}}^{*}$ exists and is given explicitly by ${\mathbf{z}}^{*}={\widehat{\mathbf{z}}}^{\left(r\right)}+{\left({\mathbf{A}}^{\left(r\right)}\right)}^{†}\left(\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{\left(r\right)}\right)$. This expression satisfies the first-order optimality condition of the convex program and ensures that ${\mathbf{z}}^{*}$ is both feasible and closest (in the examined ${\mathsf{\ell }}_2$-norm) to ${\widehat{\mathbf{z}}}^{\left(r\right)}$ among all solutions in ${\mathcal{F}}^{\left(r\right)}$. Therefore, ${\mathbf{z}}^{*}$ is (1) linear (being an affine transformation of a linear estimator), (2) unbiased (restricted to the affine feasible space), (3) minimum-norm (by construction), and (4) unique (by strict convexity). By the generalized Gauss-Markov theorem [33] (Chapter 7, Definition 4, p. 139), [34] (Chapter 8, Section 3.2, Theorem 2, p. 287), it follows that ${\widehat{\mathbf{z}}}^{*}$ is the Minimum-Norm Best Linear Unbiased Estimator (MNBLUE) under the set of linear(ized) constraints ${\mathbf{A}}^{\left(r\right)}\mathbf{z}=\mathbf{b}$.    □

Corollary 3.

Let the canonical system ${\mathbf{A}}^{\left(r\right)}{\mathbf{z}}^{\left(r\right)}=\mathbf{b}$ be consistent and of full column rank, ${\mathbf{A}}^{\left(r\right)}\in {\mathbb{R}}^{m\times n}$, $m\ge n$, and $rank\left({\mathbf{A}}^{\left(r\right)}\right)=n$. Suppose that the CLSP algorithm terminates after Step 1 and that $\mathbf{Z}={\mathbf{I}}_n$, such that ${\widehat{\mathbf{z}}}^{*}={\widehat{\mathbf{z}}}^{\left(r\right)}={\left({\mathbf{A}}^{\left(r\right)}\right)}^{†}\mathbf{b}={\left({\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}\right)}^{-1}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }\mathbf{b}$. Then, provided the classical linear model $\mathbf{b}={\mathbf{A}}^{\left(r\right)}{\mathbf{z}}^{\left(r\right)}+\epsilon$, where $\mathbb{E}\left[\epsilon \right]=\mathbf{0}$ and $Cov\left(\epsilon \right)={\sigma }^2{\mathbf{I}}_m$, holds, the CLSP estimator is equivalent to OLS while ${\widehat{\mathbf{z}}}^{*}$ is the Best Linear Unbiased Estimator (BLUE).

Corollary 4.

Let ${\widehat{\mathbf{z}}}_{\mathit{reduced}}^{\left(r\right)}={\mathbf{A}}_{{\mathbf{Z}}_{\mathit{reduced}},\mathit{reduced}}^{†}{\mathbf{b}}_{\mathrm{r}.}$ be the r-iterated estimate obtained via the Bott-Duffin inverse, defined on a subspace determined by a symmetric idempotent matrix ${\mathbf{Z}}_{\mathit{reduced}}$, and producing a conditionally unbiased estimate over that subspace. Subsequently, the estimate ${\widehat{\mathbf{z}}}_{\mathit{reduced}}^{*}$ is obtained by minimizing the squared Euclidean distance between ${\widehat{\mathbf{z}}}_{\mathit{reduced}}^{\left(r\right)}$ and ${\mathbf{z}}_{\mathit{reduced}}$, subject to the affine constraint ${\mathbf{A}}_{\mathit{reduced}}^{\left(r\right)}{\mathbf{z}}_{\mathit{reduced}}={\mathbf{b}}_{\mathrm{r}.}$, such that ${\widehat{\mathbf{z}}}_{i_E,j_E}^{*}=\left\{{\widehat{\mathbf{z}}}_{\mathit{reduced},1:m_{\mathit{subset}},1:n_{\mathit{subset}}}^{*,(i_E,j_E)}\right\}$. Then each ${\widehat{\mathbf{z}}}_{\mathit{reduced}}^{*}$ is the Minimum-Norm Best Linear Unbiased Estimator (MNBLUE) under the linear(ized) constraints defined by the canonical form of each reduced system of the CLSP, and ${\widehat{\mathbf{z}}}^{*}$ is a block-wise$m_{\mathrm{subset}}\times n_{\mathrm{subset}}$-MNBLUEcorresponding to each of the $\left(⌈{\displaystyle \frac{m}{m_{\mathit{subset}}}}⌉·⌈{\displaystyle \frac{n}{n_{\mathit{subset}}}}⌉\right)$ reduced problems.

4. Numerical Stability of the Solutions ${\widehat{\mathbf{z}}}^{\left(r\right)}$ and ${\widehat{\mathbf{z}}}^{*}$

The numerical stability of solutions in each step of the CLSP estimator — both the (iterated if $r>1$) first-step estimate ${\widehat{\mathbf{z}}}^{\left(r\right)}$ (alternatively, ${\widehat{\mathbf{z}}}_{i_E,j_E}^{\left(r\right)}=\left\{{\widehat{\mathbf{z}}}_{\mathrm{reduced},1:m_{\mathrm{subset}},1:n_{\mathrm{subset}}}^{\left(r\right),(i_E,j_E)}\right\}$) and the final solution ${\widehat{\mathbf{z}}}^{*}$ (alternatively, ${\widehat{\mathbf{z}}}_{i_E,j_E}^{*}=\left\{{\widehat{\mathbf{z}}}_{\mathrm{reduced},1:m_{\mathrm{subset}},1:n_{\mathrm{subset}}}^{*,(i_E,j_E)}\right\}$) — depends on the condition number of the design matrix, $\kappa \left({\mathbf{A}}^{\left(r\right)}\right)$, which, given its block-wise formulation in Equation (3), can be analyzed as a function of the “canonical” constraints matrix${\mathbf{C}}_{\mathrm{canon}}=\left[\mathbf{C}\mathbf{S}\right]\in {\mathbb{R}}^{k\times n}$ as a fixed, r-invariant, part of ${\mathbf{A}}^{\left(r\right)}$. For both full and reduced problems, sensitivity analysis of $\kappa \left({\mathbf{A}}^{\left(r\right)}\right)$ with respect to constraints (i.e., rows) in ${\mathbf{C}}_{\mathrm{canon}}$, especially focused on their reduction (i.e., dropping rows), is performed based on the decomposition:

$${\mathbf{A}}^{\left(r\right)}={\mathbf{B}}^{\left(r\right)}{\mathbf{C}}_{\mathrm{canon}}+{\mathbf{R}}^{\left(r\right)}={\mathbf{B}}^{\left(r\right)}\left[\begin{array}{cc}\mathbf{C}& \mathbf{S}\end{array}\right]+{\mathbf{R}}^{\left(r\right)},{\mathbf{R}}^{\left(r\right)}={\mathbf{A}}^{\left(r\right)}(\mathbf{I}-{\mathbf{P}}_{\mathcal{R}\left({\left[\mathbf{C}\mathbf{S}\right]}^{\top }\right)})$$

(7)

where ${\mathbf{B}}^{\left(r\right)}={\mathbf{A}}^{\left(r\right)}{\mathbf{C}}_{\mathrm{canon}}^{†}\in {\mathbb{R}}^{m\times k}$, denotes the constraint-induced, r-variant, part of $\mathbf{A}$, and, when the rows of ${\mathbf{A}}^{\left(r\right)}$ lie entirely within $\mathcal{R}\left({\mathbf{C}}_{\mathrm{canon}}\right)$, i.e., ${\mathbf{A}}^{\left(r\right)}={\mathbf{A}}^{\left(r\right)}{\mathbf{P}}_{\mathcal{R}\left({\left[\mathbf{C}\mathbf{S}\right]}^{\top }\right)}\iff {\mathbf{R}}^{\left(r\right)}=\mathbf{0}$,

$$\begin{array}{cc}\hfill \kappa \left({\mathbf{A}}^{\left(r\right)}\right)& \le \kappa \left({\mathbf{B}}^{\left(r\right)}\right)·\kappa \left({\mathbf{C}}_{\mathrm{canon}}\right)\hfill \\ \hfill \kappa \left({\mathbf{A}}^{\left(r\right)}\right)+\Delta {\kappa }_A& \le \left(\kappa \left({\mathbf{B}}^{\left(r\right)}\right)+\Delta {\kappa }_B\right)·\left(\kappa \left({\mathbf{C}}_{\mathrm{canon}}\right)+\Delta {\kappa }_C\right)\hfill \\ \hfill & =\kappa \left({\mathbf{B}}^{\left(r\right)}\right)·\kappa \left({\mathbf{C}}_{\mathrm{canon}}\right)+\Delta {\kappa }_B·\kappa \left({\mathbf{C}}_{\mathrm{canon}}\right)+\kappa \left({\mathbf{B}}^{\left(r\right)}\right)·\Delta {\kappa }_C+\Delta {\kappa }_B·\Delta {\kappa }_C\hfill \end{array}$$

(8)

i.e., a change in $\kappa \left({\mathbf{A}}^{\left(r\right)}\right)$ consists of a (lateral) constraint-side effect, $\kappa \left({\mathbf{B}}^{\left(r\right)}\right)·\Delta {\kappa }_C$, induced by structural modifications in ${\mathbf{C}}_{\mathrm{canon}}$, an (axial) model-side effect, $\Delta {\kappa }_B·\kappa \left({\mathbf{C}}_{\mathrm{canon}}\right)$, reflecting propagation into ${\mathbf{B}}^{\left(r\right)}$ — where ${\mathbf{B}}^{\left(r\right)}={\mathbf{A}}^{\left(r\right)}{\mathbf{C}}_{\mathrm{canon}}^{†}$ and ${\mathbf{A}}^{\left(r\right)}=f\left({\mathbf{C}}_{\mathrm{canon}},{\mathbf{M}}_{\mathrm{canon}}^{\left(r\right)}\right)$, ${\mathbf{M}}_{\mathrm{canon}}^{\left(r\right)}=g\left(\mathbf{M},{\epsilon }^{(r-1)}\right)$, and ${\epsilon }^{(r-1)}=\mathbf{b}-{\mathbf{A}}^{(r-1)}{\widehat{\mathbf{z}}}^{(r-1)}$ (for a “partitioned” decomposition of ${\mathbf{B}}^{\left(r\right)}$ given the structure of ${\mathbf{C}}_{\mathrm{canon}}^{†}$, see Appendix A.2), — and the (non-linear) term, $\Delta {\kappa }_B·\Delta {\kappa }_C$.

The condition number of ${\mathbf{C}}_{\mathrm{canon}}\in {\mathbb{R}}^{k\times n}$, $\kappa \left({\mathbf{C}}_{\mathrm{canon}}\right)$, the change of which determines the constraint-side effect, $\kappa \left({\mathbf{B}}^{\left(r\right)}\right)·\Delta {\kappa }_C$, can be itself analyzed with the help of pairwise angular measure between row vectors ${\mathbf{c}}_i^{\top },{\mathbf{c}}_j^{\top }\in {\mathbf{C}}_{\mathrm{canon}}$, where $i,j\in \{1,\dots ,k\}$. For any two ${\mathbf{c}}_i^{\top },{\mathbf{c}}_j^{\top }\in {\mathbb{R}}^{1\times n}$ (alternatively, in a centered form, ${\tilde{\mathbf{c}}}_i^{\top }={\mathbf{c}}_i^{\top }-{\overline{c}}_i,{\tilde{\mathbf{c}}}_j^{\top }={\mathbf{c}}_j^{\top }-{\overline{c}}_j\in {\mathbb{R}}^{1\times n}$, with ${\overline{c}}_i={\displaystyle \frac{1}{n}}{\sum }_{k=1}^n{c}_{i,k},{\overline{c}}_j={\displaystyle \frac{1}{n}}{\sum }_{k=1}^n{c}_{j,k}$), $i\ne j$, the cosine of the angle between them $cos\left({\theta }_{i,j}\right)$ (alternatively, a Pearson correlation coefficient $\rho =cos\left({\tilde{\theta }}_{i,j}\right)$), captures their angular alignment — values close to $\pm 1$ indicate near-collinearity and potential ill-conditioning, whereas values near 0 imply closeness to orthogonality and improved numerical stability:

$${\forall }_{{\mathbf{c}}_i^{\top },{\mathbf{c}}_j^{\top }\in {\mathbf{C}}_{\mathrm{canon}},i\ne j}cos\left({\theta }_{i,j}\right)={\displaystyle \frac{{\mathbf{c}}_i^{\top }·{\mathbf{c}}_j}{\parallel {\mathbf{c}}_i^{\top }{\parallel }_2·{\parallel {\mathbf{c}}_j^{\top }\parallel }_2}}\in [-1,1]$$

(9)

The pairwise angular measure can be aggregated across j (i.e., for each row ${\mathbf{c}}_i^{\top }$) and jointly across $i,j$ (i.e., for the full set of $\left({\displaystyle \genfrac{}{}{0pt}{}{k}{2}}\right)={\displaystyle \frac{k(k-1)}{2}}$ unique row pairs of ${\mathbf{C}}_{\mathrm{canon}}$), yielding a scalar summary statistic of angular alignment, e.g., the Root Mean Square Alignment (RMSA):

$${\mathrm{RMSA}}_i=\sqrt{{\displaystyle \frac{1}{k-1}}{\displaystyle \sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^k}{cos}^2\left({\theta }_{i,j}\right)},\mathrm{RMSA}=\sqrt{{\displaystyle \frac{2}{k(k-1)}}{\displaystyle \sum _{i=1}^{k-1}}{\displaystyle \sum _{j=i+1}^k}{cos}^2\left({\theta }_{i,j}\right)}$$

(10)

where ${\mathrm{RMSA}}_i$ captures the average angular alignment of constraint i with the rest, while RMSA evaluates the overall constraint anisotropy of ${\mathbf{C}}_{\mathrm{canon}}$. As with individual cosines $cos\left({\theta }_{i,j}\right)$ (alternatively, $\rho =cos\left({\tilde{\theta }}_{i,j}\right)$), values near 1 suggest collinearity and potential numerical instability, while values near 0 reflect angular dispersion and reduced condition number $\kappa \left({\mathbf{C}}_{\mathrm{canon}}\right)$. The use of squared cosines ${cos}^2\left({\theta }_{i,j}\right)$ (alternatively, ${\rho }^2={cos}^2\left({\tilde{\theta }}_{i,j}\right)$) in contrast with other metrics (e.g., absolute values) maintains consistency with the Frobenius norms in $\kappa \left({\mathbf{C}}_{\mathrm{canon}}\right)$ and assigns greater penalization to near-collinear constraints in ${\mathbf{C}}_{\mathrm{canon}}$. Furthermore, a combination of ${\mathrm{RMSA}}_i$, $i\in \{1,\dots ,k\}$, and the above-described marginal effects, obtained by excluding row i from ${\mathbf{C}}_{\mathrm{canon}}$ — namely, $\kappa \left({\mathbf{B}}^{\left(r\right)}\right)·\Delta {\kappa }_C^{(-i)}$, $\Delta {\kappa }_B^{(-i)}·\kappa \left({\mathbf{C}}_{\mathrm{canon}}\right)$, and $\Delta {\kappa }_B^{(-i)}·\Delta {\kappa }_C^{(-i)}$ if ${\mathbf{R}}^{\left(r\right)}=\mathbf{0}$ — together with the corresponding changes in the goodness-of-fit statistics from the solutions ${\widehat{\mathbf{z}}}^{\left(r\right)(-i)}$ and ${\widehat{\mathbf{z}}}^{*(-i)}$ (consult Section 5), can be visually presented over either all i, or a filtered subset, such that ${\mathrm{RMSA}}_i\ge \mathrm{threshold}$, as depicted in Figure 5.

5. Goodness of Fit of the Solutions ${\widehat{\mathbf{z}}}^{\left(r\right)}$ and ${\widehat{\mathbf{z}}}^{*}$

Given the hybrid structure of the CLSP estimator — comprising a least-squares-based (iterated if $r>1$) initial estimate ${\widehat{\mathbf{z}}}^{\left(r\right)}$ (or ${\widehat{\mathbf{z}}}_{i_E,j_E}^{\left(r\right)}=\left\{{\widehat{\mathbf{z}}}_{\mathrm{reduced},1:m_{\mathrm{subset}},1:n_{\mathrm{subset}}}^{\left(r\right),(i_E,j_E)}\right\}$), followed by a regularization-based final solution ${\widehat{\mathbf{z}}}^{*}$ (or ${\widehat{\mathbf{z}}}_{i_E,j_E}^{*}=\left\{{\widehat{\mathbf{z}}}_{\mathrm{reduced},1:m_{\mathrm{subset}},1:n_{\mathrm{subset}}}^{*,(i_E,j_E)}\right\}$) — and the potentially ill-posed or underdetermined nature of the problems it addresses, standard goodness-of-fit methodology — such as explained variance measures (e.g., the coefficient of determination, $R^2$), error-based metrics (e.g., Root Mean Square Error, RMSE, for comparability with the above-described RMSA), hypothesis tests (e.g., F-tests in the analysis of variance and t-tests at the coefficient level), or confidence intervals (e.g., based on Student’s t- and the normal distribution) — requires appropriate adaptation to be uniformly applied to both full and reduced problems, with some statistics applicable only to Step 1 and/or under specific conditions (namely, the ANOVA F- tests, t-tests, and confidence intervals). In response, several more robust metrics are developed below: partial $R^2$, normalized RMSE, t-test for the mean of the NRMSE, and a diagnostic interval based on the condition number.

For the explained variance measures, the block-wise formulation of the full vector of unknowns ${\mathbf{z}}^{\left(r\right)}=\left[\begin{array}{c}{\mathbf{x}}^{\left(r\right)}\\ {\left({\mathbf{y}}^{\left(r\right)}\right)}^{*}\end{array}\right]\in {\mathbb{R}}^n$, — obtained from ${\mathbf{x}}^{\left(r\right)}={\mathbf{Q}}_X\left[\begin{array}{c}{\mathbf{x}}_M^{\left(r\right)}\\ {\mathbf{x}}_L^{\left(r\right)}\end{array}\right]\in {\mathbb{R}}^p$ (where ${\mathbf{x}}_M^{\left(r\right)}\in {\mathbb{R}}^{p-l}$ is the vector of model (target) variables and ${\mathbf{x}}_L^{\left(r\right)}\in {\mathbb{R}}^l$, $l\ge 0$, is a vector of latent variables present only in the constraints) and a vector of slack variables ${\left({\mathbf{y}}^{\left(r\right)}\right)}^{*}\in {\mathbb{R}}^{n-p}$, — the design matrix ${\mathbf{A}}^{\left(r\right)}=\left[\begin{array}{c}{\mathbf{C}}_{\mathrm{canon}}\\ {\mathbf{M}}_{\mathrm{canon}}^{\left(r\right)}\end{array}\right]\in {\mathbb{R}}^{m\times n}$, where ${\mathbf{M}}_{\mathrm{canon}}^{\left(r\right)}=\left[\mathbf{M}\{\mathbf{0},r=1;{\mathbf{S}}_{\mathrm{residual}},r>1\}\right]\in {\mathbb{R}}^{(m-k)\times n}$ is the “canonical” model matrix, and the input $\mathbf{b}=\left[\begin{array}{c}{\mathbf{b}}_C\\ {\mathbf{b}}_M\end{array}\right]\in {\mathbb{R}}^m$, all necessitate a “partial” ${\mathbf{R}}^{\mathbf{2}}$ to isolate the estimated variables ${\widehat{\mathbf{x}}}^{*}$, where, given the vector $\mathbf{1}={\{1,\dots ,1\}}^{\top }$,

$$R_{\mathrm{partial}}^2=1-{\displaystyle \frac{{∥{\mathbf{b}}_M-\mathbf{M}{\widehat{\mathbf{x}}}^{*}∥}_2^2}{{∥{\mathbf{b}}_M-{\overline{b}}_M\mathbf{1}∥}_2^2}},{\overline{b}}_M={\displaystyle \frac{1}{m-k}}\sum _{i=1}^{m-k}{b}_{M,i}$$

(11)

provided that $rank\left({\mathbf{A}}^{\left(r\right)}\right)=n$, $m\ge n$, i.e., applicable strictly to (over)determined problems. If $\mathbf{C}=\varnothing \iff k=0\wedge \mathbf{y}=\varnothing \iff n=p$, then ${\mathbf{A}}^{\left(r\right)}\equiv \mathbf{M}$ and $R_{\mathrm{partial}}^2$ reduces to the classical $R^2$.

Next, for the error-based metrics, independent of $rank\left({\mathbf{A}}^{\left(r\right)}\right)$ and thus more robust across underdetermined cases (in contrast to $R_{\mathrm{partial}}^2$), a Root Mean Square Error (RMSE) can be defined for both (a) the full solution ${\widehat{\mathbf{z}}}^{*}$, via $\mathrm{RMSE}={\displaystyle \frac{1}{\sqrt{m}}}{∥\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{*}∥}_2$, serving as an overall measure of fit, and (b) its variable component ${\widehat{\mathbf{x}}}^{*}$, via ${\mathrm{RMSE}}_{\mathrm{partial}}^{\left(r\right)}={\displaystyle \frac{1}{\sqrt{m-k}}}{∥{\mathbf{b}}_{k+1:m}-\mathbf{M}{\widehat{\mathbf{x}}}^{*}∥}_2$, quantifying the fit quality with respect to the estimated target variables only (in alignment with $R_{\mathrm{partial}}^2$). Then, to enable cross-model comparisons and especially its use in hypothesis testing (see below), RMSE must be normalized — typically by the standard deviation of the reference input, ${\sigma }_{\mathbf{b}}$ or ${\sigma }_{{\mathbf{b}}_M}$, — yielding the Normalized Root Mean Square Error (NRMSE) comparable across datasets and/or models (e.g., in the t-test for the mean of the NRMSE):

$$\mathrm{NRMSE}={\displaystyle \frac{\frac{1}{\sqrt{m}}{∥\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{*}∥}_2}{{\sigma }_{\mathbf{b}}}},{\mathrm{NRMSE}}_{\mathrm{partial}}={\displaystyle \frac{\frac{1}{\sqrt{m-k}}{∥{\mathbf{b}}_{k+1:m}-\mathbf{M}{\widehat{\mathbf{x}}}^{*}∥}_2}{{\sigma }_{{\mathbf{b}}_M}}}$$

(12)

The standard deviation is preferred in statistical methodology over the alternatives (e.g., max-min scaling or range-based metrics) because it expresses residuals in units of their variability, producing a scale-free measure analogous to the standard normal distribution. It is also more common in practical implementations (such as Python’s sklearn or MATLAB).

Also, for the hypothesis tests, classical ANOVA-based F-tests and coefficient-level t-tests — relying on variance decomposition and residual degrees of freedom — are applicable exclusively to the least-squares-based (iterated if $r>1$) first-step estimate ${\widehat{\mathbf{z}}}^{\left(r\right)}$ and to the final solution ${\widehat{\mathbf{z}}}^{*}={\widehat{\mathbf{z}}}^{\left(r\right)}$, given a strictly overdetermined ${\mathbf{A}}^{\left(r\right)}{\mathbf{z}}^{\left(r\right)}=\mathbf{b}$, i.e., $rank\left({\mathbf{A}}^{\left(r\right)}\right)=n$, $m>n$ (an OLS case), and under the assumption of homoscedastic and normally distributed residuals. Then, in the case of the F-tests, the test statistics follow the distributions (a) $F\sim \mathcal{F}(q,m-n)$ and (b) $F_{\mathrm{partial}}\sim \mathcal{F}(q_{\mathrm{partial}},m-n)$, a Wald-type test on a linear restriction — with the degrees of freedom $q=n-1$ and $q_{\mathrm{partial}}=p-1$ if ${\mathbf{x}}^{\left(r\right)}$ is a vector of coefficients from a linear(ized) model with an intercept and $q=n$ and $q_{\mathrm{partial}}=p$ otherwise, — yielding (a) $H_0:{\widehat{\mathbf{z}}}^{\left(r\right)}=\mathbf{0}$ and (b) $H_0^{\mathrm{partial}}:{\widehat{\mathbf{x}}}^{\left(r\right)}=\mathbf{R}{\widehat{\mathbf{z}}}^{\left(r\right)}=\mathbf{0}$, where $\mathbf{R}\in {\mathbb{R}}^{p\times n}$ is a restriction matrix:

$$F={\displaystyle \frac{\frac{1}{q}{∥{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_2^2}{\frac{1}{m-n}{∥\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_2^2}},F_{\mathrm{partial}}={\displaystyle \frac{\frac{1}{q_{\mathrm{partial}}}{\left(\mathbf{R}{\widehat{\mathbf{z}}}^{\left(r\right)}\right)}^{\top }{\left[\mathbf{R}{\left({\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}\right)}^{-1}{\mathbf{R}}^{\top }\right]}^{-1}\left(\mathbf{R}{\widehat{\mathbf{z}}}^{\left(r\right)}\right)}{\frac{1}{m-n}{∥\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_2^2}}$$

(13)

In the case of the t-tests, (a) ${\forall }_{i\in \{1,\dots ,n\}}{t}_i\sim t_{m-n}$ (n test statistics) and (b) ${\forall }_{i\in \{1,\dots ,p\}}$$t_{\mathrm{partial},i}\equiv t_i$ (p test statistics), given that the partial result ${\widehat{\mathbf{x}}}^{\left(r\right)}$ or ${\widehat{\mathbf{x}}}^{*}={\widehat{\mathbf{x}}}^{\left(r\right)}$ is a subvector of length p of ${\widehat{\mathbf{z}}}^{\left(r\right)}$ (also ${\widehat{\mathbf{z}}}^{*}={\widehat{\mathbf{z}}}^{\left(r\right)}$), yielding (a) ${\forall }_{i\in \{1,\dots ,n\}}{H}_0:{\widehat{z}}_i^{\left(r\right)}=0$ and (b) ${\forall }_{i\in \{1,\dots ,p\}}{H}_0^{\mathrm{partial}}:{\widehat{x}}_i^{\left(r\right)}=0$:

$${\forall }_{i\in \{1,\dots ,n\}}{t}_i={\displaystyle \frac{{\widehat{z}}_i^{\left(r\right)}}{\left(\frac{1}{\sqrt{m-n}}{∥\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_2\right)·\sqrt{{\left[{\left({\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}\right)}^{-1}\right]}_{ii}}}},{\forall }_{i\in \{1,\dots ,p\}}{t}_{\mathrm{partial},i}\equiv t_i$$

(14)

An alternative hypothesis test — robust to both $rank\left({\mathbf{A}}^{\left(r\right)}\right)$ and the step of the CLSP algorithm, i.e., valid for the final solution ${\widehat{\mathbf{z}}}^{*}$ (in contrast to the classical ANOVA-based F-tests or coefficient-level t-tests) — can be a $\mathit{t}$-test for the mean of the NRMSE, comparing the observed ${\mathrm{NRMSE}}_{\mathrm{obs}}$ (as ${\mu }_0$), both full and partial, to the mean of a simulated sample ${\mathrm{NRMSE}}_{\mathrm{sim}}^{\left(\tau \right)}$, generated via Monte Carlo simulation, typically from a uniformly (a structureless baseline) or normally (a “canonical” choice) distributed random input ${\mathbf{b}}^{\left(\tau \right)}\sim \mathcal{U}(\mathbf{0},{\mathbf{I}}_m)$ or ${\mathbf{b}}^{\left(\tau \right)}\sim \mathcal{N}(\mathbf{0},{\mathbf{I}}_m)$ — both the distributions being standard, — for $\tau =1,\dots ,T$, where T is the sample size, yielding $t\sim t_{T-1}$ with $H_0:{\mathrm{NRMSE}}_{\mathrm{obs}}={\overline{\mathrm{NRMSE}}}_{\mathrm{sim}}$ and $H_1:{\mathrm{NRMSE}}_{\mathrm{obs}}\ne {\overline{\mathrm{NRMSE}}}_{\mathrm{sim}}$ for the two-sided test and $H_1:{\mathrm{NRMSE}}_{\mathrm{obs}}\gtrless {\overline{\mathrm{NRMSE}}}_{\mathrm{sim}}$ for the one-sided one:

$$t={\displaystyle \frac{{\overline{\mathrm{NRMSE}}}_{\mathrm{sim}}-{\mathrm{NRMSE}}_{\mathrm{obs}}}{\sqrt{\frac{1}{T(T-1)}{\sum }_{\tau =1}^T{\left({\mathrm{NRMSE}}_{\mathrm{sim}}^{\left(\tau \right)}-{\overline{\mathrm{NRMSE}}}_{\mathrm{sim}}\right)}^2}}},{\overline{\mathrm{NRMSE}}}_{\mathrm{sim}}={\displaystyle \frac{1}{T}}\sum _{\tau =1}^T{\mathrm{NRMSE}}_{\mathrm{sim}}^{\left(\tau \right)}$$

(15)

justified when ${\mathbf{b}}^{\left(\tau \right)}$ is normally distributed or approximately normally distributed based on the Lindeberg-Lévy Central Limit Theorem (CLT), i.e., when $T>30$ (as a rule of thumb) for ${\mathbf{b}}^{\left(\tau \right)}\sim \mathcal{U}(\mathbf{0},{\mathbf{I}}_m)$. Then, $H_0$ denotes a good fit for ${\widehat{\mathbf{z}}}^{*}$ in the sense that ${\mathrm{NRMSE}}_{\mathrm{obs}}$ does not significantly deviate from the simulated distribution, i.e., $H_0$ should not be rejected for the CLSP model to be considered statistically consistent. In practical implementations (such as Python, R, Stata, SAS, Julia, and MATLAB/Octave), T typically ranges from 50 to 50,000.

Finally, for the confidence intervals, classical formulations for (a) $z_i^{\left(r\right)}$, $i\in \{1,\dots ,n\}$, and (b) $x_i^{\left(r\right)}$, $i\in \{1,\dots ,p\}$, — based on Student’s t- and, asymptotically, the standard normal distribution, $t\sim t_{m-n}$ and $z\sim \mathcal{N}(0,1)$, — are equivalently exclusively applicable to the least-squares-based (iterated if $r>1$) first-step estimate ${\widehat{\mathbf{z}}}^{\left(r\right)}$ and to the final solution ${\widehat{\mathbf{z}}}^{*}={\widehat{\mathbf{z}}}^{\left(r\right)}$ under a strictly overdetermined ${\mathbf{A}}^{\left(r\right)}{\mathbf{z}}^{\left(r\right)}=\mathbf{b}$, i.e., $rank\left({\mathbf{A}}^{\left(r\right)}\right)=n$, $m>n$ (an OLS case), and assuming homoscedastic and normally distributed residuals. Then, provided $\alpha \in (0,1)$, such as $\alpha \in \{0.01,0.05,0.1\}$, the confidence intervals for (a) $z_i^{\left(r\right)}$ and (b) $x_i^{\left(r\right)}$ are:

$$\begin{array}{cc}\hfill & {\forall }_{i\in \{1,\dots ,n\}}{z}_i^{\left(r\right)}\in {\widehat{z}}_i^{\left(r\right)}\pm t_{m-n,1-\alpha /2}·{\displaystyle \frac{1}{\sqrt{m-n}}}{∥\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_2·\sqrt{{\left[{\left({\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}\right)}^{-1}\right]}_{ii}},\hfill \\ \hfill & {\forall }_{i\in \{1,\dots ,p\}}{x}_i^{\left(r\right)}\equiv z_i^{\left(r\right)}\hfill \end{array}$$

(16)

given that the partial result ${\widehat{\mathbf{x}}}^{\left(r\right)}$ or ${\widehat{\mathbf{x}}}^{*}={\widehat{\mathbf{x}}}^{\left(r\right)}$ is a subvector of length p of ${\widehat{\mathbf{z}}}^{\left(r\right)}$ or ${\widehat{\mathbf{z}}}^{*}={\widehat{\mathbf{z}}}^{\left(r\right)}$. For $n>30$, the distribution of the test statistic, $t_{m-n,1-\alpha /2}$, approaches a standard normal one, $Z_{1-\alpha /2}$, based on the CLT, yielding an alternative notation for the confidence intervals:

$$\begin{array}{cc}\hfill & {\forall }_{i\in \{1,\dots ,n\}}{z}_i^{\left(r\right)}\in {\widehat{z}}_i^{\left(r\right)}\pm Z_{1-\alpha /2}·{\displaystyle \frac{1}{\sqrt{m-n}}}{∥\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{\left(r\right)}∥}_2·\sqrt{{\left[{\left({\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}\right)}^{-1}\right]}_{ii}},\hfill \\ \hfill & {\forall }_{i\in \{1,\dots ,p\}}{x}_i^{\left(r\right)}\equiv z_i^{\left(r\right)}\hfill \end{array}$$

(17)

An alternative “confidence interval” — equivalently robust to both $rank\left({\mathbf{A}}^{\left(r\right)}\right)$ and the step of the CLSP algorithm, i.e., valid for the final solution ${\widehat{\mathbf{z}}}^{*}$ (in contrast to the classical Student’s t- and, asymptotically, the standard normal distribution-based formulations) — can be constructed deterministically via condition-weighted bands, relying on componentwise numerical sensitivity. Let the residual vector $\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{*}$ be a (backward) perturbation in $\mathbf{b}$, i.e., $\Delta \mathbf{b}=\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{*}$. Then, squaring both sides of the classical first-order inequality ${\displaystyle \frac{\parallel \Delta {\widehat{\mathbf{z}}}^{*}{\parallel }_2}{\parallel {\widehat{\mathbf{z}}}^{*}{\parallel }_2}}\le \kappa \left({\mathbf{A}}^{\left(r\right)}\right)·{\displaystyle \frac{{\parallel \Delta \mathbf{b}\parallel }_2}{{\parallel \mathbf{b}\parallel }_2}}$ yields ${\displaystyle \frac{{\sum }_{i=1}^n{(\Delta {\widehat{z}}_i^{*})}^2}{{\sum }_{i=1}^n{\left({\widehat{z}}_i^{*}\right)}^2}}\le \kappa {\left({\mathbf{A}}^{\left(r\right)}\right)}^2·{\displaystyle \frac{{∥\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{*}∥}_2^2}{{\parallel \mathbf{b}\parallel }_2^2}}$, where $\kappa \left({\mathbf{A}}^{\left(r\right)}\right)=\parallel {\mathbf{A}}^{\left(r\right)}{\parallel }_2{\parallel {\left({\mathbf{A}}^{\left(r\right)}\right)}^{†}\parallel }_2={\displaystyle \frac{{\sigma }_{\max }\left({\mathbf{A}}^{\left(r\right)}\right)}{{\sigma }_{rank\left({\mathbf{A}}^{\left(r\right)}\right)}\left({\mathbf{A}}^{\left(r\right)}\right)}}$ — ${\sigma }_{\max }$ and ${\sigma }_{rank\left({\mathbf{A}}^{\left(r\right)}\right)}$ being the biggest and smallest singular values of ${\mathbf{A}}^{\left(r\right)}$ — iff ${\sigma }_{rank\left({\mathbf{A}}^{\left(r\right)}\right)}>0\wedge {\parallel \mathbf{b}\parallel }_2>0$. Under a uniform relative squared perturbation (a heuristic allocation of error across components as a simplification, not an assumption) of ${\mathbf{z}}^{*}$, rearranging terms and taking square roots of both sides gives ${\forall }_{i\in \{1,\dots ,n\}}|\Delta {\widehat{z}}_i^{*}|\le |{\widehat{z}}_i^{*}|·\kappa \left({\mathbf{A}}^{\left(r\right)}\right)·{\displaystyle \frac{{∥\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{*}∥}_2}{{\parallel \mathbf{b}\parallel }_2}}$ and a condition-weighted “confidence” band:

$${\forall }_{i\in \{1,\dots ,n\}}{z}_i^{*}\in {\widehat{z}}_i^{*}·\left(1\pm \kappa \left({\mathbf{A}}^{\left(r\right)}\right)·{\displaystyle \frac{{∥\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{*}∥}_2}{{\parallel \mathbf{b}\parallel }_2}}\right),{\forall }_{i\in \{1,\dots ,p\}}{x}_i^{*}\equiv z_i^{*}$$

(18)

which is a diagnostic interval based on the condition number for ${\mathbf{z}}^{*}$, consisting of (1) a canonical-form system conditioning, $\kappa \left({\mathbf{A}}^{\left(r\right)}\right)$, (2) a normalized model misfit, ${\displaystyle \frac{{∥\mathbf{b}-{\mathbf{A}}^{\left(r\right)}{\widehat{\mathbf{z}}}^{*}∥}_2}{{\parallel \mathbf{b}\parallel }_2}}$, and (3) the “scale” of the final solution, ${\widehat{\mathbf{z}}}^{*}$, without probabilistic assumptions. A perturbation in one or more of these three components, e.g., caused by a change in RMSA resulting from dropping selected constraints (consult Section 4), will affect the limits of ${\mathbf{z}}^{*}$. Under a perfect fit, the interval collapses to $z_i^{*}\in {\widehat{z}}_i^{*}\pm 0$ and, for $\kappa \left({\mathbf{A}}^{\left(r\right)}\right)\gg 1$, tends to be very wide.

To summarize, the CLSP estimator requires a goodness-of-fit framework that reflects both its algorithmic structure and numerical stability. Classical methods remain informative but are valid only under strict overdetermination ($m>n$), full-rank design matrix ${\mathbf{A}}^{\left(r\right)}$, and distributional assumptions (primarily, homoscedasticity and normality of residuals). In contrast, the proposed alternatives — (a) partial $R_{\mathrm{partial}}^2$, (b) NRMSE and ${\mathrm{NRMSE}}_{\mathrm{partial}}$, (c) t-tests for the mean of NRMSE with the help of Monte Carlo sampling, and (d) condition-weighted confidence bands — are robust to rank deficiency and ill-posedness, making them applicable in practical implementations (namely, CLSP packages in Python, R, and Stata).

6. Special Cases of CLSP Problems: APs, CMLS, and LPRLS/QPRLS

The structure of the design matrix ${\mathbf{A}}^{\left(r\right)}$ in the CLSP estimator, as defined in Equation (3) (consult Section 3), allows for accommodating a selection of special-case problems, out of which three cases are covered, but the list may not be exhaustive. Among others, the CLSP estimator is, given its ability to address ill-posed or underdetermined problems under linear(ized) constraints, efficient in addressing what can be referred to as allocation problems (APs) (or, for flow variables, tabular matrix problems, TMs) — i.e., in most cases, underdetermined problems involving matrices of dimensions $m\times p$ to be estimated, subject to known row and column sums, with the degrees of freedom (i.e., nullity) equal to $n-rank\left({\mathbf{A}}^{\left(r\right)}\right)\ge mp+s^{*}-({\displaystyle \frac{m}{i}}+{\displaystyle \frac{p}{j}})-q$, where $0\le s^{*}\le ({\displaystyle \frac{m}{i}}+{\displaystyle \frac{p}{j}})+q$ is the number of active (non-zero) slack variables and $0\le q\le rank\left({\mathbf{M}}_{\mathrm{partial}}\right)\le min(k,mp)$ is the number of known model (target) variables (e.g., a zero diagonal), — whose design matrix, ${\mathbf{A}}^{\left(r\right)}\in {\mathbb{R}}^{({\displaystyle \frac{m}{i}}+{\displaystyle \frac{p}{j}}+k)\times n}$, comprises a (by convention) row-sum-column-sum constraints matrix $\mathbf{C}={\left[{\left(\left({\mathbf{I}}_{{\displaystyle \frac{m}{i}}}\otimes {\mathbf{1}}_i\right)\otimes {\mathbf{1}}_p\right)}^{\top }|{\left({\mathbf{1}}_m\otimes \left({\mathbf{I}}_{{\displaystyle \frac{p}{j}}}\otimes {\mathbf{1}}_j\right)\right)}^{\top }\right]}^{\top }\in {\mathbb{R}}^{({\displaystyle \frac{m}{i}}+{\displaystyle \frac{p}{j}})\times mp}$ (where $\mathbf{1}=\{1,\dots ,1\}$), a model matrix$\mathbf{M}={\mathbf{Q}}_M^L\left[\begin{array}{c}{\mathbf{M}}_{\mathrm{partial}}\\ \mathbf{0}\end{array}\right]{\mathbf{Q}}_M^R\in {\mathbb{R}}^{k\times mp}$ — in a trivial case, ${\mathbf{M}}_{\mathrm{partial}}\subseteq {\mathbf{I}}_{mp}$ — a sign slack matrix$\mathbf{S}={\mathbf{Q}}_S^L\left[\begin{array}{c}{\mathbf{S}}_s\subseteq {\forall }_{{\sigma }_1,\dots ,{\sigma }_s\in \{-1,+1\}}diag\left({\sigma }_1,\dots ,{\sigma }_s\right)\\ \mathbf{0}\end{array}\right]{\mathbf{Q}}_S^R\in {\mathbb{R}}^{({\displaystyle \frac{m}{i}}+{\displaystyle \frac{p}{j}})\times (n-mp)}$, and either a zero matrix $\mathbf{0}$ (if $r=1$) or a (standard) reverse-sign slack matrix ${\mathbf{S}}_{\mathrm{residual}}\in {\mathbb{R}}^{k\times (n-mp)}$ (provided $r>1$):

$${\mathbf{A}}^{\left(r\right)}=\left[\begin{array}{cc}\left[\begin{array}{c}\left({\mathbf{I}}_{{\displaystyle \frac{m}{i}}}\otimes {\mathbf{1}}_i\right)\otimes {\mathbf{1}}_p\\ {\mathbf{1}}_m\otimes \left({\mathbf{I}}_{{\displaystyle \frac{p}{j}}}\otimes {\mathbf{1}}_j\right)\end{array}\right]& \mathbf{S}\\ \mathbf{M}& \left\{\begin{array}{cc}\mathbf{0},\hfill & \mathrm{initial}\mathrm{iteration},r=1\hfill \\ {\mathbf{S}}_{\mathrm{residual}},\hfill & \mathrm{subsequent}\mathrm{iterations},r>1\hfill \end{array}\right.\end{array}\right]$$

(19)

where ${\displaystyle \frac{m}{i}}$ and ${\displaystyle \frac{p}{j}}$ denote groupings of the m rows and p columns, respectively, into i and j homogeneous blocks (when no grouped row or column sums are available, $i=j=1$); with real-world examples including: (a) input-output tables (national, inter-country, or satellite); (b) structural matrices (e.g., trade, country-product, investment, or migration); (c) financial clearing and budget-balancing; and (d) data interpolations (e.g., quarterly data from annual totals). Given the available literature in 2025, the first pseudoinverse-based method of estimating input-output tables was proposed in Pereira-López et al. [95] and the first APs (TMs)-based study was conducted in Bolotov [16], attempting to interpolate a world-level (in total, 232 countries and dependent territories) “bilateral flow” matrix of foreign direct investment (FDI) for the year 2013, based on known row and column totals from UNCTAD data (aggregate outward and inward FDI). The estimation employed a minimum-norm least squares solution under APs (TMs)-style equality constraints, based on the Moore-Penrose pseudoinverse of a “generalized Leontief structural matrix”, rendering it equivalent to the initial (non-regularized) solution in the CLSP framework, ${\widehat{\mathbf{z}}}^{\left(1\right)}={\mathbf{A}}_{\mathbf{Z}}^{\left(1\right)†}\mathbf{b}$. As a rule of thumb, AP (TM)-type problems are rank deficient for $m>2\wedge p>2$ with the CLSP being a unique (if $\alpha \in (0,1]$) and an MNBLUE (if $\alpha =1$) estimator (consult Section 3).

In addition to the APs (TMs), the CLSP estimator is also efficient in addressing what can be referred to as constrained-model least squares (CMLS) (or, more generally, regression problems, RPs) — i.e., (over)determined problems involving vectors of dimension p to be estimated with (standard) OLS degrees of freedom, subject to, among others, linear(ized) data matrix $\mathbf{D}\in {\mathbb{R}}^{k\times p}$-related (in)equality constraints, ${\forall }_{i\in \{1,\dots ,t\}}{\mathbf{C}}_i\equiv {\mathbf{T}}_i\mathbf{D}\in {\mathbb{R}}^{m_i\times p}$, where ${\mathbf{T}}_i\in {\mathbb{R}}^{m_i\times k}$ is a transformation matrix of $\mathbf{D}$, such as the i-th difference or shift (lag/lead), and ${\mathbf{C}}_i\mathbf{x}⋛{\mathbf{\gamma }}_i$, — whose design matrix, ${\mathbf{A}}^{\left(r\right)}\in {\mathbb{R}}^{({\sum }_{i=1}^u{m}_i+k)\times n}$, consists of a u-block constraints matrix $\mathbf{C}=\left[\begin{array}{c}{\mathbf{C}}_1\\ \dots \\ {\mathbf{C}}_u\end{array}\right]\in {\mathbb{R}}^{\left({\sum }_{i=1}^u{m}_i\right)\times p}$ (${\forall }_{i\in \{1,\dots ,u\}}{\mathbf{C}}_i\in {\mathbb{R}}^{m_i\times p}$), $u\ge t$, a data matrix as the model matrix $\mathbf{M}\equiv \mathbf{D}\in {\mathbb{R}}^{k\times p}$, a (standard) sign slack matrix $\mathbf{S}\in {\mathbb{R}}^{\left({\sum }_{i=1}^u{m}_i\right)\times (n-p)}$, and either a zero matrix $\mathbf{0}$ (if $r=1$) or a (standard) reverse-sign slack matrix ${\mathbf{S}}_{\mathrm{residual}}\in {\mathbb{R}}^{k\times (n-p)}$ (when $r>1$):

$${\mathbf{A}}^{\left(r\right)}=\left[\begin{array}{cc}\left[\begin{array}{c}{\mathbf{C}}_1\\ \dots \\ {\mathbf{C}}_u\end{array}\right]& \mathbf{S}\\ \mathbf{D}& \left\{\begin{array}{cc}\mathbf{0},\hfill & \mathrm{initial}\mathrm{iteration},r=1\hfill \\ {\mathbf{S}}_{\mathrm{residual}},\hfill & \mathrm{subsequent}\mathrm{iterations},r>1\hfill \end{array}\right.\end{array}\right]$$

(20)

with real-world examples including: (a) previously infeasible textbook-definition econometric models of economic (both micro and macro) variables (e.g., business-cycle models), (b) additional constraints applied to classical econometric models (e.g., demand analysis). Given the available literature in 2025, the first studies structurally resembling CMLS (RPs) were conducted in Bolotov [15], Bolotov [96], focusing on the decomposition of the long U.S. GDP time series into trend ($y_t$) and cyclical ($c_t$) components — using exponential trend, moving average, Hodrick-Prescott filter, and Baxter-King filter — under constraints on its first difference, based on the U.S. National Bureau of Economic Research (NBER) delimitation of the U.S. business cycle. $c_t$ was then modeled with the help of linear regression (OLS), based on an n-th order polynomial with extreme values smoothed via a factor of ${\displaystyle \frac{1}{{10}^{n-i}}}$ and simultaneously penalized via an n-th or ($n+1$)-th root, $c_t={\beta }_0+{\beta }_1\int \left(\sqrt[\{n,n+1\}]{{\prod }_{i=1}^{i=n}(t-t_i)·{\displaystyle \frac{1}{{10}^{n-i}}}}\right)dt+{\epsilon }_t$, where ${\beta }_0,{\beta }_1$ are the model parameters, ${\forall }_{1\le i\le n}{t}_i$ are the values of the (externally given) stationary points, and ${\epsilon }_t$ is the error. The order of the polynomial and the ad-hoc selection of smoothing and penalizing factors, however, render such a method inferior to the CLSP. Namely, the unique (if $\alpha \in (0,1]$) and an MNBLUE (if $\alpha =1$) CLSP estimator (consult Section 3) allows (a) presence of inequality constraints and (b) their ill-posed formulations.

Finally, the CLSP estimator can be used to address (often unsolvable using classical solvers) underdetermined and/or ill-posed — caused by too few and/or mutually inconsistent or infeasible constraints, as in the sense of Whiteside et al. [93] and Blair [76], Blair [94], — LP and QP problems, an approach hereafter referred to as linear/quadratic programming via regularized least squares (LPRLS/QPRLS): CLSP substitutes the original objective function of the LP/QP problem with the canonical form ${\mathbf{A}}^{\left(r\right)}{\mathbf{z}}^{\left(r\right)}=\mathbf{b}$ (i.e., focusing solely on the problem’s constraints, without distinguishing between the LP and QP cases) with ${\widehat{\mathbf{x}}}_M^{*}={Proj}_{{\mathbf{x}}_M^{\left(r\right)}}\left\{{\mathbf{Q}}_X\left[\begin{array}{c}{\mathbf{x}}_M^{\left(r\right)}\\ {\mathbf{x}}_L^{\left(r\right)}\end{array}\right]\in {\mathbb{R}}^p:{\forall }_{1\le i\le s}{g}_i({\mathbf{x}}_M^{\left(r\right)},{\mathbf{x}}_L^{\left(r\right)})⋛{\gamma }_i,{\forall }_j{h}_j({\mathbf{x}}_M^{\left(r\right)},{\mathbf{x}}_L^{\left(r\right)})={\eta }_j\right\}\subseteq {\mathbf{z}}^{*}$ being the solution of the original problem — where ${\mathbf{Q}}_X$ is a permutation matrix, and $g_i(·)⋛{\gamma }_i$ and $h_j(·)={\eta }_j$ represent the linear(ized) inequality and equality constraints,${\forall }_{i,j}{\gamma }_i,{\eta }_j\in \mathbf{b}$ — and the degrees of freedom being equal — under the classical formalization of a primal LP/QP problem (consult Section 2) — to $n-rank\left({\mathbf{A}}^{\left(r\right)}\right)\ge (p+s)-(m_{ub}+m_{eq}+m_{nn})=p-m_{eq}$, where p is the length of ${\mathbf{x}}_M^{\left(r\right)}$, $s=m_{ub}+m_{nn}$ is the number of introduced slack variables, while $m_{ub}\ge 0$, $m_{eq}\ge 0$, and $m_{nn}\equiv p$ denote the numbers of upper-bound, equality, and non-negativity constraints, respectively (under the standard assumption that all model variables are constrained to be non-negative). In the LPRLS/QPRLS, the design matrix, ${\mathbf{A}}^{\left(r\right)}\in {\mathbb{R}}^{(m_{ub}+m_{eq}+m_{nn})\times n}$, by the definition of LP/QP, is r-invariant and consists of a block constraints matrix $\mathbf{C}=\left[\begin{array}{c}{\mathbf{C}}_{ub}\\ {\mathbf{C}}_{eq}\\ {\mathbf{C}}_{nn}\end{array}\right]\in {\mathbb{R}}^{(m_{ub}+m_{eq}+m_{nn})\times p}$, where ${\mathbf{C}}_{ub}\in {\mathbb{R}}^{m_{ub}\times p}$, ${\mathbf{C}}_{eq}\in {\mathbb{R}}^{m_{eq}\times p}$, and ${\mathbf{C}}_{nn}\equiv {\mathbf{I}}_p\in {\mathbb{R}}^{p\times p}$, a (standard) sign slack matrix $\mathbf{S}\in {\mathbb{R}}^{(m_{ub}+m_{eq}+m_{nn})\times (n-p)}$, and ${\mathbf{M}}_{\mathrm{canon}}=\varnothing$:

$${\mathbf{A}}^{\left(r\right)}={\mathbf{A}}^{\left(1\right)}=\left[\begin{array}{cc}\left[\begin{array}{c}{\mathbf{C}}_{ub}\\ {\mathbf{C}}_{eq}\\ {\mathbf{C}}_{nn}\end{array}\right]& \mathbf{S}\end{array}\right]$$

(21)

Given the available literature in 2025, the first documented attempts to incorporate LS into LP and QP included, among others, Dax’s LS-based steepest-descent algorithm [85] and the primal-dual NNLS-based algorithm (LSPD) [86], whose structural restrictions — in terms of well-posedness and admissible constraints — can be relaxed through LPRLS/QPRLS, still guaranteeing a unique (if $\alpha \in (0,1]$) and an MNBLUE (if $\alpha =1$) solution (consult Section 3).

7. Monte Carlo Experiment and Numerical Examples

The Monte Carlo experiment and the simulated examples presented below illustrate the applicability of the CLSP to the AP (TM), CMLS (RP), and LPRLS/QPRLS problems. The standardized form of $\mathbf{C}$ in the APs (TMs), i.e., $\mathbf{C}=f(m,i,p,j)$, with given $\mathbf{S}$ and r, enables large-scale simulation of the distribution of NRMSE — a robust goodness-of-fit metric for the CLSP — through repeated random trials under varying matrix dimensions $m\times p$, row and column group sizes i and j, and composition of $\mathbf{M}$, inferring asymptotic percentile means and standard deviations for practical implementations (e.g., the t-test for the mean). This text presents the results of a simulation study on the distribution of the NRMSE from the first-step estimate ${\widehat{\mathbf{z}}}^{\left(1\right)}$ (implementable in standard econometric software without CLSP modules) for $i=1$, $j=1$, $\mathbf{S}=\varnothing$, $\mathbf{Z}={\mathbf{I}}_{m+p+k}\iff {\left({\mathbf{A}}_{\mathbf{Z}}^{\left(1\right)}\right)}^{†}={\left({\mathbf{A}}^{\left(1\right)}\right)}^{†}$, and $r=1$, conducted in Stata/MP (set to version 14.0) with the help of Mata’s 128-bit floating-point cross-product-based quadcross() for greater precision and SVD-based svsolve() with a strict tolerance of c("epsdouble") for increased numerical stability, and a random-variate seed 123456789 (see Appendix B.1 containing a list of dependencies [97] and the Stata .do-file). In this 50,000-iterations Monte Carlo simulation study, spanning matrix dimensions from 1×2 and 2×1 up to 50×50, random normal input vectors $\mathbf{b}\sim \mathcal{N}(\mathbf{0},{\mathbf{I}}_{m+p+k})$ were generated for each run, applied with and without zero diagonals (i.e., if, under $l=0\iff {\mathbf{x}}_L=\varnothing$, ${\mathbf{x}}_M=\mathbf{x}={\mathbf{z}}_{1:mp}\in {\mathbb{R}}^{mp}$ is reshaped into a ($m\times p$)-matrix ${\mathbf{X}}_M=\mathbf{X}\in {\mathbb{R}}^{m\times p}$, then ${\forall }_{i\in \{1,\dots ,min(m,p\left)\right\}}{\mathbf{X}}_{ii}=0$). Thus, a total of 249.9 million observations ($=(50·50-1)·2·50,000$) was obtained via the formula $\mathrm{NRMSE}={\displaystyle \frac{1}{\sqrt{m}}}{\parallel \mathbf{b}-{\mathbf{A}}^{\left(1\right)}\left({\left({\mathbf{A}}^{\left(1\right)}\right)}^{†}\mathbf{b}\right)\parallel }_2/{\sigma }_{\mathbf{b}}$, resulting in 4,998 aggregated statistics of the asymptotic distribution, assuming convergence: mean, sd, skewness, kurtosis, min, max, as well as p1–p99, as presented in Figure 6 (depicting the intensity of its first two moments, mean and sd) and Table 2 (reporting a subset of all thirteen statistics for $m,pmod10=0$).

From the results of the Monte Carlo experiment, it is observable that: (1) the NRMSE from the first-step estimate ${\widehat{\mathbf{z}}}^{\left(1\right)}$ and its distributional statistics exhibit increasing stability and boundedness as matrix dimensions $m\times p$ increase — specifically, for $m,p\gg 1$, both the mean and sd of NRMSE tend to converge, indicating improved estimator performance (e.g., under no diagonal restrictions, at $m=p=5$, $\mathtt{mean}=1.84$ and $\mathtt{sd}=1.29$, while at $m=p=45$, $\mathtt{mean}=0.21$ and $\mathtt{sd}=0.16$); (2) the zero-diagonal constraints (${\forall }_{i\in \{1,\dots ,min(m,p\left)\right\}}{\mathbf{X}}_{ii}=0$) reduce the degrees of freedom and lead to a much more uniform and dimensionally stable distribution of NRMSE across different matrix sizes (e.g., at $m=p=5$, mean drops to $1.17$ and sd to $0.85$, while at $m=p=45$, mean rises slightly to $0.26$ and sd to $0.19$); and (3) the distribution of NRMSE exhibits mild right skew and leptokurtosis — less so for the zero-diagonal case: under no diagonal restriction, the average skewness is $0.97$ with a coefficient of variation of $38.16\%$ (i.e., $>0.00$), and the average kurtosis is $4.15$ with a variation of $304.61\%$; whereas under the zero-diagonal constraint, skewness averages $0.95$ with only $8.57\%$ variation, and kurtosis averages $3.76$ with $25.94\%$ variation (i.e., $>3.00$). To sum up, the larger and less underdetermined (i.e., in an AP (TM) problem, the larger the ${\mathbf{M}}^{\left(1\right)}$ block of ${\mathbf{A}}^{\left(1\right)}$) the canonical system ${\mathbf{A}}^{\left(1\right)}{\mathbf{z}}^{\left(1\right)}=\mathbf{b}$, the lower the estimation error (i.e., the mean NRMSE) and its variance, and, provided a stable limiting law for NRMSE exists, the skewness and kurtosis from the simulations drift toward 0 and 3, respectively, consistent with — but not proving — a convergence toward a standard-normal distribution.

For an illustration of the applicability of the CLSP to the APs (TMs), CMLS (RPs), and LPRLS/QPRLS, the author-developed cross-platform Python 3 module, pyclsp (version ≥ 1.3.0, available on PyPI) [98], — together with its interface co-modules for the APs (TMs), pytmpinv (version ≥ 1.2.0, available on PyPI) [99], and the LPRLS/QPRLS, pylppinv (version ≥ 1.3.0, available on PyPI) [100] (a co-module for the CMLS (RPs) is not provided, as the variability of ${\forall }_{i\in \{1,\dots ,t\}}{\mathbf{C}}_i\equiv {\mathbf{T}}_i\mathbf{D}$ within the $\mathbf{C}=\left[\begin{array}{c}{\mathbf{C}}_1\\ \dots \\ {\mathbf{C}}_u\end{array}\right]$, $u\ge t$, block of ${\mathbf{A}}^{\left(r\right)}$ prevents efficient generalization), — is employed for a sample (dataset) of i random elements ${\mathcal{X}}_i$ (scalars, vectors, or matrices depending on context) and their transformations, drawn independently from the standard normal distribution, i.e., $\{{\mathcal{X}}_1,\dots ,{\mathcal{X}}_i\}\stackrel{\mathrm{iid}}{\sim }\mathcal{N}(0,1)$, under the mentioned random-variate seed 123456789 (configured for consistency with the above-described Monte Carlo experiment). Thus, using the cases of (a) a (non-negative symmetric) input-output table and (b) a (non-negative) zero-diagonal trade matrix as two AP (TM) numerical examples, this text simulates problems similar to the ones addressed in Pereira-López et al. [95] and Bolotov [16]: an underdetermined $\mathbf{X}\in {\mathbb{R}}^{m\times p}$ is estimated — from row sums, column sums, and k known values of two randomly generated matrices (a) ${\mathbf{X}}_{\mathrm{true}}=0.5·\left(|{\mathcal{X}}_1|+|{\mathcal{X}}_1{|}^{\top }\right)$ and (b) ${\mathbf{X}}_{\mathrm{true},i,j}=\left(|{\mathcal{X}}_1\right|{)}_{i,j}\left(1-{\mathbb{1}}_{\{i=j\le min(m,p\left)\right\}}\right)$, ${\left({\mathcal{X}}_1\right)}_{i,j}\stackrel{\mathrm{iid}}{\sim }\mathcal{N}(0,1)$, subject to (a) ${\mathbf{X}}_{i,j}\ge 0$ and (b) ${\mathbf{X}}_{i,j}\ge 0\wedge {\forall }_{i\in \{1,\dots ,min(m,p\left)\right\}}{\mathbf{X}}_{ii}=0$ — via CLSP assuming $\mathbf{x}=vec\left(\mathbf{X}\right)\subseteq \mathbf{z}$, $\mathbf{M}\equiv {\left({\mathbf{I}}_{mp}\right)}_{\mathcal{I}}$, $\mathbf{b}=\left[\begin{array}{c}{\mathbf{b}}_{\mathrm{row}\mathrm{sums}}\\ {\mathbf{b}}_{\mathrm{column}\mathrm{sums}}\\ {\mathbf{b}}_{\mathrm{known}\mathrm{values}}\end{array}\right]=\left[\begin{array}{c}{\mathbf{X}}_{\mathrm{true}}\mathbf{1}\\ {\left({\mathbf{1}}^{\top }{\mathbf{X}}_{\mathrm{true}}\right)}^{\top }\\ {\left(vec\left({\mathbf{X}}_{\mathrm{true}}\right)\right)}_{\mathcal{I}}\end{array}\right]$, and $\mathcal{I}\sim \mathcal{U}\left\{1,\dots ,mp\right\}\wedge \left|\mathcal{I}\right|=k$. The Python 3 code, based on pytmpinv and two other modules, installable by executing pip install matplotlib numpy pytmpinv==1.2.0, for (a) $m=p=20$ and $k=0.1·mp=40$ and (b) $m=p=40$ — estimated with the help of $m_{\mathrm{subset}}\times n_{\mathrm{subset}}$-MNBLUE across $\left(⌈{\displaystyle \frac{m}{m_{\mathrm{subset}}}}⌉·⌈{\displaystyle \frac{n}{n_{\mathrm{subset}}}}⌉\right)=\left(⌈{\displaystyle \frac{40}{20-1}}⌉·⌈{\displaystyle \frac{40}{20-1}}⌉\right)=9$ reduced models $m_{\mathrm{reduced}}=p_{\mathrm{reduced}}\le 20$ assuming an estimability constraint of $mp\le 400$ — and $k=0.2·mp=320$, with a non-iterated ($r=1$), unique ($\alpha >0$), and MNBLUE ($\alpha =1$) (two-step) CLSP solution, is implemented in (a) Listing  and (b) Listing  (e.g., to be executed in a Jupyter Notebook).

Listing 1. Simulated numerical example for a symmetric input-output-table-based AP (TM) problem.

In case (a) (see Listing  for code), the number of model (target) variables $\#{\widehat{\mathbf{x}}}_M^{*}$ is $mp=400$ with a nullity of $dim(ker\left({\mathbf{A}}^{\left(1\right)}\right))=n-rank\left({\mathbf{A}}^{\left(1\right)}\right)=321$, corresponding to 80.25% of the total unknowns. Given the simulation of ${\mathbf{X}}_{\mathrm{true}}\in {\mathbb{R}}^{20\times 20}$, a matrix unknown in real-world applications, — i.e., CLSP is used to estimate the elements of an existing matrix from its row sums, ${\mathbf{X}}_{\mathrm{true}}\mathbf{1}\in {\mathbb{R}}^{20}$, column sums, ${\left({\mathbf{1}}^{\top }{\mathbf{X}}_{\mathrm{true}}\right)}^{\top }\in {\mathbb{R}}^{20}$, and a randomly selected 10% of its entries, ${\left(vec\left({\mathbf{X}}_{\mathrm{true}}\right)\right)}_{\mathcal{I}}\in {\mathbb{R}}^{40}$, — the model’s goodness of fit can be measured by a user-defined $R^2\approx max\left(0,1-{\displaystyle \frac{{∥vec\left({\mathbf{X}}_{\mathrm{true}}\right)-{\widehat{\mathbf{x}}}_M^{*}∥}_2^2}{{∥vec\left({\mathbf{X}}_{\mathrm{true}}\right)-\mathbf{1}\frac{1}{400}{\sum }_{i=1}^{400}\left(vec{\left({\mathbf{X}}_{\mathrm{true}}\right)}_i\right)∥}_2^2}}\right)\approx 0.278256$ (i.e., CLSP achieves an improvement of $\Delta R^2=0.178256$ over a hypothetical naïve predictor reproducing the known 10% entries of ${\mathbf{X}}_{\mathrm{true},\mathcal{I}}$ and yielding a $R^2\approx 0.1$, but still a modest value of $R^2$, i.e. a relatively large error, $\mathrm{NRMSE}=12.649111$) with ${\widehat{\mathbf{x}}}_M^{*}$ lying within wide condition-weighted diagnostic intervals ${\widehat{\mathbf{x}}}_M^{*}\in \left[min\left({\widehat{\mathbf{x}}}_{M,\mathrm{lower}}^{*}\right),max\left({\widehat{\mathbf{x}}}_{M,\mathrm{upper}}^{*}\right)\right]=\left[-8.0847\times {10}^{15},8.0847\times {10}^{15}\right]$ reflecting the ill-conditioning of the strongly rectangular ${\mathbf{A}}^{\left(1\right)}\in {\mathbb{R}}^{80\times 400}$, $\kappa \left({\mathbf{A}}^{\left(1\right)}\right)=4.9968\times {10}^{14}$, and only the left-sided Monte Carlo-based t-test for the mean of the NRMSE (on a sample of 30 NRMSE values obtained by substituting $\mathbf{b}$ with ${\left({\mathbf{b}}^{\left(\tau \right)}\right)}_i\stackrel{\mathrm{iid}}{\sim }\mathcal{N}(0,1)$) suggesting consistency with expectations (i.e., with the $\mathrm{p}-\mathrm{value}\approx 1$). In terms of numerical stability, $\kappa \left({\mathbf{C}}_{\mathrm{canon}}\right)=4.8023\times {10}^{14}\approx \kappa \left({\mathbf{A}}^{\left(1\right)}\right)$ with a low ${\mathrm{RMSA}}_{i\in \{1,\dots ,40\}}=\mathrm{RMSA}=0.0400$, as confirmed by the correlogram produced in matplotlib, indicating a well-chosen constraints matrix, even given the underdeterminedness of the model (which also prevents a better fit).

Listing 2. Simulated numerical example for a (non-negative) trade-matrix-based AP (TM) problem.

In case (b) (see Listing  for code), the corresponding number of model (target) variables in each of the $\left(⌈{\displaystyle \frac{40}{19}}⌉·⌈{\displaystyle \frac{40}{19}}⌉\right)=9$ reduced design submatrices ${\mathbf{A}}_{\mathrm{subset}}^{\left(1\right)}\subset {\mathbf{A}}_{\mathrm{reduced}}^{\left(1\right)}$, $\#{\widehat{\mathbf{x}}}_{M,\mathrm{reduced}}^{*}$, is $m_{\mathrm{subset}}{p}_{\mathrm{subset}}\le 361$ — where $m_{\mathrm{subset}}\le 20-1=19$ and $p_{\mathrm{subset}}\le 20-1=19$, one row and one column being reserved for ${\sum }_{j\notin \mathcal{J}}{a}_{i,j}^{\left(1\right)}$ and ${\sum }_{i\notin \mathcal{I}}{a}_{i,j}^{\left(1\right)}$, which enter ${\mathbf{A}}_{\mathrm{reduced}}^{\left(1\right)}$ as $\mathbf{S}=\left[\begin{array}{c}{\mathbf{I}}_{m-m_{\mathrm{subset}}}\mathbf{0}\\ \mathbf{0}{\mathbf{I}}_{p-p_{\mathrm{subset}}}\end{array}\right]$ and the vector of slack variables ${\left({\widehat{\mathbf{y}}}_{\mathrm{reduced}}^{*}\right)}^{*}=\left[\begin{array}{c}{\widehat{\mathbf{y}}}_{\mathrm{reduced}}^{*}\\ {\widehat{\mathbf{y}}}_{\mathrm{residual},\mathrm{reduced}}^{*}\end{array}\right]\subset {\widehat{\mathbf{z}}}_{\mathrm{reduced}}^{*}$, so that ${\mathbf{C}}_{\mathrm{reduced}}{\mathbf{x}}_{\mathrm{reduced}}+{\mathbf{S}}_{\mathrm{reduced}}{\mathbf{y}}_{\mathrm{reduced}}^{*}={\mathbf{b}}_{\mathrm{r}.,1:k}$ (i.e., the slack matrix compensates for the unaccounted row and column sums in the reduced models as opposed to the full one, such as case (a)) — with a nullity $dim(ker\left({\mathbf{A}}_{\mathrm{reduced}}^{\left(1\right)}\right))=n_{\mathrm{reduced}}-rank\left({\mathbf{A}}_{\mathrm{reduced}}^{\left(1\right)}\right)\in [2,290]$ (per reduced model), corresponding to 25.00%–72.68% of the total unknowns (per reduced model) — to compare, a full model, under the same inputs and no computational constraints (i.e., $mp\nleqq m_E{p}_E$), would have a nullity $dim(ker\left({\mathbf{A}}^{\left(1\right)}\right))=n-rank\left({\mathbf{A}}^{\left(1\right)}\right)=1168$, corresponding to 73.00% of the total unknowns. In the examined example, — based on ${\mathbf{X}}_{\mathrm{true}}\mathbf{1}\in {\mathbb{R}}^{40}$, ${\left({\mathbf{1}}^{\top }{\mathbf{X}}_{\mathrm{true}}\right)}^{\top }\in {\mathbb{R}}^{40}$, and a randomly selected 20% of entries of the true matrix ${\mathbf{X}}_{\mathrm{true}}\in {\mathbb{R}}^{40\times 40}$, — the reduced-model block solution’s goodness of fit could not be efficiently measured by a user-defined $R^2\approx max\left(0,1-{\displaystyle \frac{{∥vec\left({\mathbf{X}}_{\mathrm{true}}\right)-vec\left(\left\{{\widehat{\mathbf{x}}}_{M,\mathrm{reduced},1:m_{\mathrm{subset}},1:n_{\mathrm{subset}}}^{*,(i_E,j_E)}\right\}\right)∥}_2^2}{{∥vec\left({\mathbf{X}}_{\mathrm{true}}\right)-\mathbf{1}\frac{1}{1600}{\sum }_{i=1}^{1600}\left(vec{\left({\mathbf{X}}_{\mathrm{true}}\right)}_i\right)∥}_2^2}}\right)\approx 0$ (i.e., the block matrix constructed from $\left(⌈{\displaystyle \frac{40}{19}}⌉·⌈{\displaystyle \frac{40}{19}}⌉\right)=9$ reduced-model estimates led to $R^2\approx max(0,-5.435377)$, but to an error proportionate to the one in case (a), ${\mathrm{NRMSE}}_{\mathrm{reduced}}\in [3.464102,15.748016]$ (per reduced model)) — in contrast, a full model would achieve $R^2\approx max\left(0,1-{\displaystyle \frac{{∥vec\left({\mathbf{X}}_{\mathrm{true}}\right)-{\widehat{\mathbf{x}}}_M^{*}∥}_2^2}{{∥vec\left({\mathbf{X}}_{\mathrm{true}}\right)-\mathbf{1}\frac{1}{1600}{\sum }_{i=1}^{1600}\left(vec{\left({\mathbf{X}}_{\mathrm{true}}\right)}_i\right)∥}_2^2}}\right)\approx 0.278427$, but at a cost of a greater error, $\mathrm{NRMSE}=29.427878$ — with ${\widehat{\mathbf{x}}}_{M,\mathrm{reduced}}^{*}$ lying within strongly varying condition-weighted diagnostic intervals ${\widehat{\mathbf{x}}}_{M,\mathrm{reduced}}^{*}\in \left[min\left({\widehat{\mathbf{x}}}_{M,\mathrm{lower},\mathrm{reduced}}^{*}\right),max\left({\widehat{\mathbf{x}}}_{M,\mathrm{upper},\mathrm{reduced}}^{*}\right)\right]$, where $min\left({\widehat{\mathbf{x}}}_{M,\mathrm{lower},\mathrm{reduced}}^{*}\right)\in [-330.638753,-2.5009\times {10}^{-19}]$ (per reduced model) and $max\left({\widehat{\mathbf{x}}}_{M,\mathrm{upper},\mathrm{reduced}}^{*}\right)\in [-2.5456\times {10}^{-19},336.547225]$ (per reduced model), and varying results of Monte Carlo-based t-tests for the mean of the ${\mathrm{NRMSE}}_{\mathrm{reduced}}$ (on a sample of 30 ${\mathrm{NRMSE}}_{\mathrm{reduced}}$ values obtained by substituting ${\mathbf{b}}_{\mathrm{r}.}$ with ${\left({\mathbf{b}}_{\mathrm{r}.}^{\left(\tau \right)}\right)}_i\stackrel{\mathrm{iid}}{\sim }\mathcal{N}(0,1)$), where p-values range $4.8022\times {10}^{-09}$–1.000000 (left-sided), 0.000000–1.000000 (two-sided), and $3.7021\times {10}^{-20}$–1.000000 (right-sided) (per reduced model) — alternatively, a full model would lead to wider condition-weighted diagnostic intervals ${\widehat{\mathbf{x}}}_M^{*}\in \left[min\left({\widehat{\mathbf{x}}}_{M,\mathrm{lower}}^{*}\right),max\left({\widehat{\mathbf{x}}}_{M,\mathrm{upper}}^{*}\right)\right]=\left[-1.5357\times {10}^{16},1.5357\times {10}^{16}\right]$ (i.e., reflecting the ill-conditioning of the strongly rectangular ${\mathbf{A}}^{\left(1\right)}\in {\mathbb{R}}^{400\times 1600}$, $\kappa \left({\mathbf{A}}^{\left(1\right)}\right)=2.2196\times {10}^{14}$) and only the left-sided Monte Carlo-based t-test for the mean of the NRMSE (on a sample of 30 NRMSE values obtained by substituting $\mathbf{b}$ with ${\left({\mathbf{b}}^{\left(\tau \right)}\right)}_i\stackrel{\mathrm{iid}}{\sim }\mathcal{N}(0,1)$) suggesting consistency with expectations (i.e., with the $\mathrm{p}-\mathrm{value}\approx 1$). In terms of numerical stability, $\kappa \left({\mathbf{C}}_{\mathrm{canon},\mathrm{reduced}}\right)\in [2.236068,6.244998]<\kappa \left({\mathbf{A}}_{\mathrm{reduced}}^{\left(1\right)}\right)\in [4.509742,9.325735]$ (per reduced model), which indicates well-conditioning of all the reduced models — conversely, in a full model, $\kappa \left({\mathbf{C}}_{\mathrm{canon}}\right)=1.0907\times {10}^{14}<\kappa \left({\mathbf{A}}^{\left(1\right)}\right)$ (therefore, a full model ensures an overall better fit but a lower fit quality, i.e., a trade-off).

As the CMLS (RP) numerical example, this text addresses a problem similar to the one solved in Bolotov [15], Bolotov [96]: a coefficient vector $\mathbf{\beta }$ from a (time-series) linear regression model in its classical (statistical) notation $y_t={\mathbf{x}}_t\mathbf{\beta }+{ϵ}_t,\mathrm{s}.\mathrm{t}.{\sum }_{j=1}^{i+1}{\beta }_j=c\wedge {\forall }_{t\in \mathbf{P}}\Delta y_t=0\wedge {\forall }_{t\in \mathbf{P}}{\Delta }^2{y}_t⪌0$ with ${\Delta }^n{y}_t$ denoting n-th order differences (i.e., the discrete analogue of ${\displaystyle \frac{d^n{y}_t}{d{t}^n}}$) — where $y_t$ is the dependent variable, ${\mathbf{x}}_t=\left[1,x_{t,1},\dots ,x_{t,i}\right]$ is the vector of regressors with a constant, ${ϵ}_t$ is the model’s error (residual), c is a constant, and $\mathbf{P}⫋\{t_j:j=1,\dots ,k\}$ is the set of stationary points — is estimated on a simulated sample (dataset) with the help of CLSP assuming ${\mathbf{x}}_M\equiv \mathbf{\beta }\in {\mathbb{R}}^p$ and $\mathbf{D}\equiv \mathbf{X}=\left[\mathbf{1},{\mathbf{x}}_1,\dots ,{\mathbf{x}}_i\right]\in {\mathbb{R}}^{k\times p}$, $p=i+1$. The Python 3 code, based on numpy and pyclsp modules, installable by executing pip install numpy pyclsp==1.3.0, for $k=500$, $p=6$, $c=1$, $\mathbf{y}=\mathbf{D}{\mathbf{\beta }}_{\mathrm{true}}+\mathbf{ϵ}$, where ${\forall }_{j\in \{2,\dots ,6\}}{\mathbf{D}}_{·,j}\equiv {\mathcal{X}}_j$, ${\left({\mathcal{X}}_j\right)}_i\stackrel{\mathrm{iid}}{\sim }\mathcal{N}(0,1)$, ${\mathcal{X}}_j\in {\mathbb{R}}^{500}$, ${\mathbf{\beta }}_{\mathrm{true}}={\displaystyle \frac{{\mathcal{X}}_{\mathbf{\beta }}}{{\mathbf{1}}^{\top }{\mathcal{X}}_{\mathbf{\beta }}}}$, ${\left({\mathcal{X}}_{\mathbf{\beta }}\right)}_i\stackrel{\mathrm{iid}}{\sim }\mathcal{N}(0,1)$, ${\mathcal{X}}_{\mathbf{\beta }}\in {\mathbb{R}}^6$, $\mathbf{ϵ}\equiv {\mathcal{X}}_{\mathbf{ϵ}}$, ${\left({\mathcal{X}}_{\mathbf{ϵ}}\right)}_i\stackrel{\mathrm{iid}}{\sim }\mathcal{N}(0,1)$, ${\mathcal{X}}_{\mathbf{ϵ}}\in {\mathbb{R}}^{500}$, and $\mathbf{P}=\{t:\Delta \mathbf{y}=0\}$, with a non-iterated ($r=1$), unique ($\alpha >0$), and MNBLUE ($\alpha =1$) two-step CLSP solution (consistent with the ones from cases (a) and (b) for APs (TMs)), is implemented in Listing  (e.g., for a Jupyter Notebook).

Listing 3. Simulated numerical example for a (time-series) stationary-points-based CMLS (RP) problem.

Compared to the true values ${\mathbf{\beta }}_{\mathrm{true}}=[0.1687,-0.3593,0.5969,-0.0628,-0.1742,0.8307]$, the CMLS (RP) estimate is $\widehat{\mathbf{\beta }}={\widehat{\mathbf{x}}}_M^{*}=[0.2024,-0.1121,0.2095,-0.0396,-0.0296,0.2611]$ with a modest $R_{\mathrm{partial}}^2=0.293328$ (i.e., a greater error, ${\mathrm{NRMSE}}_{\mathrm{partial}}=0.840638$), moderate condition-weighted diagnostic intervals ${\widehat{\mathbf{x}}}_M^{*}\in [-29.7780,30.3003]$, and only the right-sided Monte Carlo-based t-test for the mean of the ${\mathrm{NRMSE}}_{\mathrm{partial}}$ (on a sample of 30 ${\mathrm{NRMSE}}_{\mathrm{partial}}$ values obtained by substituting $\mathbf{b}$ with ${\left({\mathbf{b}}^{\left(\tau \right)}\right)}_i\stackrel{\mathrm{iid}}{\sim }\mathcal{N}(0,1)$) — in the example under scrutiny, ${\mathrm{NRMSE}}_{\mathrm{partial}}$ is preferable to NRMSE due to ${\widehat{\mathbf{y}}}^{*}\gg {\widehat{\mathbf{x}}}^{*}$ — suggesting consistency with expectations (i.e., with the $\mathrm{p}-\mathrm{value}\approx 1$). Similarly, in terms of numerical stability, $\kappa \left({\mathbf{C}}_{\mathrm{canon}}\right)=138.6109>\kappa \left({\mathbf{A}}^{\left(1\right)}\right)=127.8287$, indicating that the constraint block is ill-conditioned, most likely, because of the imposed data- (rather than theory-) based definition of stationary points, which rendered Step 2 infeasible (${\widehat{\mathbf{x}}}_M^{*}={\widehat{\mathbf{x}}}_M^{\left(1\right)}$) (limiting the fit quality).

Finally, in the case of a LPRLS/QPRLS numerical example, this text simulates potentially ill-posed LP (QP) problems, similar to the ones addressed in Whiteside et al. [93] and Blair [76], Blair [94]: a solution ${\widehat{\mathbf{x}}}^{*}⋛0$ in its classical LP (QP) notation ${\mathbf{x}}^{*}=arg{max}_{\mathbf{x}\in {\mathbb{R}}^p}{\mathbf{c}}^{\top }\mathbf{x}$ (${\mathbf{x}}^{*}=arg{min}_{\mathbf{x}\in {\mathbb{R}}^p}\left\{\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Q}\mathbf{x}+{\mathbf{c}}^{\top }\mathbf{x}\right\}$), $\mathrm{s}.\mathrm{t}.\mathbf{A}\mathbf{x}\le \mathbf{b}\equiv {\mathbf{A}}_{ub}\mathbf{x}\le {\mathbf{b}}_{ub}\wedge {\mathbf{A}}_{eq}\mathbf{x}={\mathbf{b}}_{eq}$ — where $\mathbf{Q}\in {\mathbb{R}}^{p\times p}$ is a symmetric positive definite matrix, $\mathbf{c}\in {\mathbb{R}}^p$, $\mathbf{A}=\left[\begin{array}{c}{\mathbf{A}}_{ub}\\ {\mathbf{A}}_{eq}\end{array}\right]\in {\mathbb{R}}^{(m_{ub}+m_{eq})\times p}$, and $\mathbf{b}=\left[\begin{array}{c}{\mathbf{b}}_{ub}\\ {\mathbf{b}}_{eq}\end{array}\right]\in {\mathbb{R}}^{m_{ub}+m_{eq}}$ — is estimated from two randomly generated (coefficient) matrices ${\mathbf{A}}_{ub}\equiv {\mathcal{X}}_1$, ${\left({\mathcal{X}}_1\right)}_{i,j}\stackrel{\mathrm{iid}}{\sim }\mathcal{N}(0,1)$, and ${\mathbf{A}}_{eq}\equiv {\mathcal{X}}_2$, ${\left({\mathcal{X}}_2\right)}_{i,j}\stackrel{\mathrm{iid}}{\sim }\mathcal{N}(0,1)$, and two (right-hand-side) vectors ${\mathbf{b}}_{ub}\equiv {\mathcal{X}}_3$, ${\left({\mathcal{X}}_3\right)}_i\stackrel{\mathrm{iid}}{\sim }\mathcal{N}(0,1)$, and ${\mathbf{b}}_{eq}\equiv {\mathcal{X}}_4$, ${\left({\mathcal{X}}_4\right)}_i\stackrel{\mathrm{iid}}{\sim }\mathcal{N}(0,1)$, via CLSP assuming $\mathbf{C}=\left[\begin{array}{c}{\mathbf{A}}_{ub}\\ {\mathbf{A}}_{eq}\end{array}\right]$, $\mathbf{S}={\mathbf{Q}}_S^L\left[\begin{array}{c}{\mathbf{S}}_s\subseteq {\mathbf{I}}_{m_{ub}}\\ \mathbf{0}\end{array}\right]{\mathbf{Q}}_S^R$, and $\mathbf{b}=\left[\begin{array}{c}{\mathbf{b}}_{ub}\\ {\mathbf{b}}_{eq}\end{array}\right]$, where ${\mathbf{Q}}_S^L$ and ${\mathbf{Q}}_S^R$ are permutation matrices, while omitting $\mathbf{Q}$ and $\mathbf{c}$. The Python 3 code, based on numpy and pylppinv modules, installable by executing pip install numpy pylppinv==1.3.0, for $m_{ub}=50$, $m_{eq}=25$, and $p=500$, with a non-iterated ($r=1$), unique ($\alpha >0$), and MNBLUE ($\alpha =1$) (two-step) CLSP solution (analogously consistent with the ones from cases (a) and (b) for APs (TMs) and the one from CMLS (RPs)), is implemented in Listing  (e.g., for a Jupyter Notebook).

Listing 4. Simulated numerical example for a underdetermined and potentially infeasible LP (QP) problem.

The nullity (i.e., underdeterminedness) of the CLSP design matrix, ${\mathbf{A}}^{\left(1\right)}$, is $dim(ker\left({\mathbf{A}}^{\left(1\right)}\right))=n-rank\left({\mathbf{A}}^{\left(1\right)}\right)=475$, corresponding to 86.36% of the total unknowns, accompanied by a greater (relatively high) error, $\mathrm{NRMSE}=12.247449$, moderate condition-weighted diagnostic intervals ${\widehat{\mathbf{x}}}_M^{*}\in [-1.2919,1.3918]$, and only the right-sided Monte Carlo-based t-test for the mean of the NRMSE (on a sample of 30 NRMSE values obtained by substituting $\mathbf{b}$ with ${\left({\mathbf{b}}^{\left(\tau \right)}\right)}_i\stackrel{\mathrm{iid}}{\sim }\mathcal{N}(0,1)$) suggesting consistency with expectations (i.e., with the $\mathrm{p}-\mathrm{value}\approx 1$). Similarly, in terms of numerical stability, $\kappa \left({\mathbf{C}}_{\mathrm{canon}}\right)=2.2168=\kappa \left({\mathbf{A}}^{\left(1\right)}\right)=2.2168$, indicating that the constraint block is well-conditioned despite a strongly rectangular ${\mathbf{A}}^{\left(1\right)}\in {\mathbb{R}}^{75\times 550}$. Hence, CLSP can be efficiently applied to AP (TM), CMLS (RP), and LPRLS/QPRLS special cases using Python 3, and the reader may wish to experiment with the code in Listings – by relaxing its uniqueness (setting $\alpha =0$) or the MNBLUE characteristic (setting $\alpha \in [0,1)$).

8. Discussion and Conclusions

The proposed modular two-step Convex Least Squares Programming (CLSP) estimator attempts to supersede historical efforts developed mainly in the 1960s–1990s to address ill-posed problems in both LS and LP (QP) — namely, the mentioned Wolfe [84], Stoer [80], Waterman [81], Dax [85], Osborne [78], George and Osborne [79], Barnes et al. [86], Donoho and Tanner [77], Grafarend and Awange [82], Bei et al. [87], and Qian et al. [83] — by correcting (1) the estimation, based on a Moore-Penrose pseudoinverse ${\mathbf{A}}^{†}$ (which has already been applied in an ad hoc manner in the field, particularly in econometrics in the 2010s; consult Pereira-López et al. [95] and Bolotov [16]) and a Bott-Duffin inverse ${\mathbf{A}}_{\mathbf{Z}}^{†}$ (based on the literature research, this text may be the first to consider a constrained generalized inverse to generalize an ${\mathbf{A}}^{†}$-based estimator) with the help of a (2) regularization-based correction [4] (pp. 477–479), achieving greater fit quality in terms of ensuring strict satisfaction of problem constraints (previously infeasible under an ${\mathbf{A}}^{†}$-based solution), while yielding a (potentially) unique and MNBLUE solution (by setting $\alpha =1$). At the same time, CLSP is well grounded in the algorithmic logic of the works in question, such as Wolfe [84] and Dax [85], as further formalized by, e.g., Osborne [78]. However, the primary strength of the framework in question lies in its “generality”: in contrast to the (individual) problem-specific algorithms, such as the ones developed by Übi [88], Übi [89], Übi [90], Übi [91], Übi [92] — routines combining (NN)LS, LP, and QP within a single framework to solve both well- and ill-posed convex optimization problems, — CLSP provides a unified estimator with a proper methodology for numerical stability and goodness-of-fit assessment, centered around RMSA and NRMSE. A secondary strength is its practical “versatility”: through at least three classes of special cases — APs (TMs), CMLS (RPs), and LPRLS/QPRLS, — it allows estimation of (rectangular) input-output tables (national, inter-country, or satellite), structural matrices (such as trade, country-product, investment, or migration data), and previously infeasible textbook-definition models of economic variables and/or LP (QP) problems. Furthermore, the inclusion of ${\mathbf{A}}_{\mathbf{Z}}^{†}$ allows probabilistic and Bayesian applications, since its subspace restrictions can be interpreted in terms of random variables, prior distributions, as well as Monte Carlo-based inference. To conclude, CLSP has already been fully implemented in Python 3 (pyclsp, pytmpinv, and pylppinv) and is to be ported into R and Stata in the future, making it a readily available alternative to purely theoretical or not freely accessible algorithms or proprietary packages, and applicable in a real-world setting (such as econometric research, replacing [15,16,95] with the two-step corrected estimation).

Appendix A.1

Let ${\mathbf{A}}^{\left(r\right)}\in {\mathbb{R}}^{m\times n}$ be a real matrix, and let $\mathbf{Z}\in {\mathbb{R}}^{n\times n}$ be a symmetric idempotent matrix (i.e., ${\mathbf{Z}}^2=\mathbf{Z}={\mathbf{Z}}^{\top }$) defining a subspace restriction. Then, the Bott-Duffin inverse ${\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}\in {\mathbb{R}}^{n\times m}$ is given by the formula ${\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}={\left(\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}\mathbf{Z}\right)}^{†}\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }$, as defined in Equation (5), where the projected Gram matrix $\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}\mathbf{Z}$ is computable via its SVD:

$$\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}\mathbf{Z}=\mathbf{U}\Sigma {\mathbf{V}}^{\top }=\mathbf{U}\Sigma {\mathbf{U}}^{\top }$$

(A1)

where $\mathbf{U},\mathbf{V}\in {\mathbb{R}}^{n\times n}$ are orthogonal (since the projected Gram matrix is symmetric positive semidefinite, $\mathbf{U}=\mathbf{V}$), and $\Sigma =diag({\sigma }_1,\dots ,{\sigma }_{\mathrm{rank}},0,\dots ,0)\in {\mathbb{R}}^{n\times n}$ contains the singular values. Hence, for ${\Sigma }^{†}=diag({\sigma }_1^{-1},\dots ,{\sigma }_{\mathrm{rank}}^{-1},0,\dots ,0)$, $\mathbf{U}{\Sigma }^{†}{\mathbf{U}}^{\top }$ substitutes ${\left(\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}\mathbf{Z}\right)}^{†}$:

$${\left({\mathbf{A}}_{\mathbf{Z}}^{\left(r\right)}\right)}^{†}={\left(\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }{\mathbf{A}}^{\left(r\right)}\mathbf{Z}\right)}^{†}\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }=\left(\mathbf{U}{\Sigma }^{†}{\mathbf{U}}^{\top }\right)\mathbf{Z}{\left({\mathbf{A}}^{\left(r\right)}\right)}^{\top }$$

(A2)

Appendix A.2

Let $rank\left(\left[\mathbf{C}\mathbf{S}\right]\right)\le rank\left(\mathbf{C}\right)+rank\left(\mathbf{S}\right)$ in the analysis of numerical stability of the first estimate ${\widehat{\mathbf{z}}}^{\left(r\right)}$, as defined in Equation (2). Then a decomposition of ${\mathbf{B}}^{\left(r\right)}$ reveals how the interaction of constraints $\mathbf{C}$ and slack structure $\mathbf{S}$ forms the curvature of the solution space:

$$\begin{array}{cc}\hfill {\mathbf{B}}^{\left(r\right)}& ={\mathbf{A}}^{\left(r\right)}{\left[\begin{array}{cc}\mathbf{C}& \mathbf{S}\end{array}\right]}^{†}\hfill \\ \hfill & =f\left(\left[\begin{array}{cc}\mathbf{C}& \mathbf{S}\end{array}\right],g\left(\mathbf{M},\mathbf{b}-{\mathbf{A}}^{(r-1)}{\widehat{\mathbf{z}}}^{(r-1)}\right)\right)·{\left[\begin{array}{cc}{\mathbf{C}}^{\top }\mathbf{C}& {\mathbf{C}}^{\top }\mathbf{S}\\ {\mathbf{S}}^{\top }\mathbf{C}& {\mathbf{S}}^{\top }\mathbf{S}\end{array}\right]}^{†}\left[\begin{array}{c}{\mathbf{C}}^{\top }\\ {\mathbf{S}}^{\top }\end{array}\right]\hfill \end{array}$$

(A3)

While the block ${\mathbf{C}}^{\top }\mathbf{C}$ captures the contribution of the constraints and ${\mathbf{S}}^{\top }\mathbf{S}$ reflects the curvature introduced by slack variables, the off-diagonal blocks ${\mathbf{C}}^{\top }\mathbf{S}$ and ${\mathbf{S}}^{\top }\mathbf{C}$ encode cross-interactions between constraints and slacks: if these off-diagonal blocks vanish, the effects are orthogonal and additive (in the column space sense since ${\mathbf{C}}^{\top }\mathbf{S}=0$ implies $\mathcal{R}\left(\mathbf{C}\right)\perp \mathcal{R}\left(\mathbf{S}\right)$); otherwise, $\mathbf{C}$ is orthogonalized with respect to $\mathcal{R}\left(\mathbf{S}\right)$, which is effected by $(\mathbf{I}-\mathbf{S}{\mathbf{S}}^{†})\mathbf{C}$. Hence, the presence of $\mathbf{S}$ in $\left[\mathbf{C}\mathbf{S}\right]$ can improve or worsen the ill-posedness of the canonized problem: namely, ${\mathbf{C}}^{\top }\mathbf{S}\ne 0$ tends to worsen the condition number $\kappa \left(\left[\mathbf{C}\mathbf{S}\right]\right)$.

Appendix B.1

Listing B.1. Stata code from the Monte Carlo experiment dedicated to mapping the distribution of NRMSE.