*! version 1.1.0 20aug2025 I I Bolotov program define lppinvl2, eclass byable(recall) version 16.0 /* The Linear Programming via Regularized Least Squares (LPPinv) is a two- stage estimation method that reformulates linear programs as structured least-squares problems. Based on the Convex Least Squares Programming (CLSP) framework, LPPinv solves linear inequality, equality, and bound constraints by (1) constructing a canonical constraint system and computing a pseudoinverse projection, followed by (2) a second pseudoinverse estimation to improve fit, coherence, and achieve Ridge- style regularization (based on the pure Stata/Mata CLSP implementation). LPPinv is intended for underdetermined and ill-posed linear problems, for which standard solvers fail. Author: Ilya Bolotov, MBA, Ph.D. Date: 20 August 2025 */ tempname vars b c s // replay last result if replay() { if _by() { error 190 } cap conf mat e(z) if _rc { di as err "results of lppinvl2 not found" exit 301 } clspl2 // print summary exit 0 } // syntax syntax /// [anything] [if] [in], [ /// AUB(string) AEQ(string) BUB(string) BEQ(string) /// LOWERbound(numlist miss) UPPERbound(numlist miss) /// REPLACEvalue(numlist miss max=1) r(passthru) Z(passthru) /// RCOND(passthru) TOLerance(real -1) ITERationlimit(int 1) noFINAL /// CONDTOLerance(passthru) * /// ] // adjust and preprocess options for the LPPinvL2 Mata code if ustrregexm( `"`anything'"', "^estat", 1) { _estat `=ustrregexrf(`"`anything'"', "^estat", "", 1)', `options' exit 0 } else if ! ustrregexm( `"`anything'"', "^(non|free)",1) { di as error "invalid mode; specify {bf:nonneg} or {bf:free}" exit 198 } if ("`aub'`aeq'" == "") { di as err "at least one of aub() or aeq() is required" exit 198 } if ("`aub'" != "") { cap unab varlist : `aub' if _rc cap conf mat `aub' else loc aub `varlist' if _rc { di as err "aub: `aub'" di as err "invalid varlist or matrix" exit 198 } } if ("`aeq'" != "") { cap unab varlist : `aeq' if _rc cap conf mat `aeq' else loc aeq `varlist' if _rc { di as err "aeq: `aeq'" di as err "invalid varlist or matrix" exit 198 } } if ("`bub'`beq'" == "") { di as err "at least one of bub() or beq() is required" exit 198 } if ("`bub'" != "") { cap unab varlist : `bub', max(1) loc rc = _rc if _rc cap conf mat `bub' else loc bub `varlist' if _rc & `rc' != 103 { di as err "bub: `bub'" di as err "invalid varname or matrix" exit cond(`rc' == 103, 103, 198) } else if `rc' == 103 error 103 } if ("`beq'" != "") { cap unab varlist : `beq', max(1) loc rc = _rc if _rc cap conf mat `beq' else loc beq `varlist' if _rc & `rc' != 103 { di as err "beq: `beq'" di as err "invalid varname or matrix" exit cond(`rc' == 103, 103, 198) } else if `rc' == 103 error 103 } if ("`lowerbound'" != "") { numlist " `lowerbound'", miss mata: lowerbound = colshape(strtoreal(tokens("`r(numlist)'")),1) } else mata: lowerbound = J(1,1,.) if ("`upperbound'" != "") { numlist " `upperbound'", miss mata: upperbound = colshape(strtoreal(tokens("`r(numlist)'")),1) } else mata: upperbound = J(1,1,.) if (`iterationlimit' < 1) { di as err "r() and iterationlimit() must be at least 1" exit 198 } if ! (`tolerance' == -1 | (`tolerance' >= 0 & `tolerance' < 1)) { di as err "tolerance() must be -1 or lie in [0,1)" exit 198 } sca `vars' = "" // markout [if] [in] foreach name in `slackvars' `model' `brow' `bcol' `bval' { cap conf var `name' if ! _rc sca `vars' = `vars' + " `name'" } tempvar touse mark `touse' `if' `in' markout `touse' `= `vars'' mata: A1 =J(0,0,.) if "`aub'" != "" { mata: if (! hasmissing( _st_varindex(tokens("`aub'" )))) /// A1 = st_data(., _st_varindex(tokens("`aub'" )), /// "`touse'" );; mata: if (J(0,0,.) != st_matrix("`aub'" )) /// A1 = st_matrix("`aub'" );; } mata: A2 =J(0,0,.) if "`aeq'" != "" { mata: if (! hasmissing( _st_varindex(tokens("`aeq'" )))) /// A2 = st_data(., _st_varindex(tokens("`aeq'" )), /// "`touse'" );; mata: if (J(0,0,.) != st_matrix("`aeq'" )) /// A2 = st_matrix("`aeq'" );; } mata: b1 =J(0,1,.) if "`bub'" != "" { mata: if (! hasmissing( _st_varindex(tokens("`bub'" )))) /// b1 = st_data(., _st_varindex(tokens("`bub'" )), /// "`touse'" );; mata: if (J(0,0,.) != st_matrix("`bub'" )) /// b1 = st_matrix("`bub'" );; } mata: b2 =J(0,1,.) if "`beq'" != "" { mata: if (! hasmissing( _st_varindex(tokens("`beq'" )))) /// b2 = st_data(., _st_varindex(tokens("`beq'" )), /// "`touse'" );; mata: if (J(0,0,.) != st_matrix("`beq'" )) /// b2 = st_matrix("`beq'" );; } mata: replace_value = strtoreal( "`replacevalue'" ); mata: st_local( "tolerance", `tolerance' >= 0 ? /// strofreal(`tolerance') : /// strofreal(sqrt(epsilon(1))) ); loc iterationlimit = cond(`iterationlimit' < ., `iterationlimit', 50) // estimate with the help of the LPPinvL2 Mata code mata: if ( rows(A1) != rows(b1)) { /// errprintf("A_ub and b_ub must have the same number of " + /// "rows: " + strofreal(rows(A1)) + " vs " + /// strofreal(rows(b1)) ); /// exit(198); /// };; /// if ( rows(A2) != rows(b2)) { /// errprintf("A_eq and b_eq must have the same number of " + /// "rows: " + strofreal(rows(A2)) + " vs " + /// strofreal(rows(b2)) ); /// exit(198); /// };; /// if (cols(A1) & cols(A2) & cols(A1) != cols(A2)) { /// errprintf("A_ub and A_eq must have the same number of " + /// "columns: " + strofreal(cols(A1)) + " vs " + /// strofreal(cols(A2)) ); /// exit(198); /// };; /// if (ustrregexm(`"`anything'"', "^non",1) & /// "`lowerbound'`upperbound'"!= "" & /// (any(lowerbound :<0) | any(upperbound :<0))) { /// errprintf("Negative lower or upper bounds are not " + /// "allowed in linear programs\n" ); /// exit(198); /// };; /// if (rows(upperbound)> 1 & rows(upperbound)!=(n=max((cols(A1), /// cols(A2) )))| /// rows(lowerbound)> 1 & rows(lowerbound)!= n ) { /// errprintf("Bounds length ("+strofreal(rows(lowerbound)) + /// "," +strofreal(rows(upperbound)) + /// ") does not match number of variables " + /// strofreal(n) + "\n" ); /// exit(198); /// };; h = (rows(upperbound) == 1 ? J(n,1, upperbound) /// : upperbound); /// l = ustrregexm(`"`anything'"', "^non",1) /// ? J(n,1, 0) /// : (rows(lowerbound) == 1 ? J(n,1, lowerbound) /// : lowerbound); /// b = select((b=b1\b2\h\l), (nempty=b :<.)); /// if (b == J(0,1,.) ) { /// errprintf("RHS error: bub() and beq() are mis" + /// "specified\n" ); /// exit(198); /// };; C =((A1 != J(0,0,.) ? A1 : J(0,cols(A2),.)) \ /// (A2 != J(0,0,.) ? A2 : J(0,cols(A1),.)) \ /// J(2,1, I(n)))[selectindex( /// nempty), .]; /// if (C == J(0,0,.) ) { /// errprintf("RHS error: aub() and aeq() are mis" + /// "specified\n" ); /// exit(198); /// };; S =((S=I(rows(A1))\J(rows(A2),rows(A1),0)), /// J(rows(S ),2*n, 0) \ /// J(n,cols(S), 0), I(n),J(n,n,0) \ /// J(n,cols(S)+n,0),-I(n) )[selectindex( /// nempty), .]; /// if (rows(C ) != rows(S )| rows(C ) != rows(b )) { /// errprintf("Row mismatch: C="+strofreal(rows(C )) + /// ", S="+strofreal(rows(S )) + /// ", b="+strofreal(rows(b ))+"\n" ); /// exit(198); /// };; st_matrix("`b'", b ); /// st_matrix("`c'", C ); /// st_matrix("`s'", select(S, colsum(S ) :!= 0)); mata: st_local( "c", st_matrix("`c'") != J(0,0,.) ? /// " const( `c' )" : "" ); /// st_local( "s", st_matrix("`s'") != J(0,0,.) ? /// " slack( `s' )" : "" ); eret loc cmdline lppinvl2 `0' eret loc cmd lppinvl2 eret loc estat_cmd lppinvl2 estat eret loc title LPPinv estimation qui clspl2 general, b(`b') `c' `s' `r' `z' `rcond' tol(`tolerance') /* */ iter(`iterationlimit') `final' `condtolerance' /* */ nested /* replace out-of-bound values with replace_value */ if ( "`lowerbound'`upperbound'"!= "") /// mata: x = colshape(st_matrix("e(x)"), 1); /// idx = selectindex((x :< l :- `tolerance') :| /// (x :> h :+ `tolerance')); /// x[rows(idx) ? idx : J(0,1,.)] = /// J(rows(idx),1,replace_value); /// st_matrix("e(x)", x ); /* display the result summary */ clspl2 // print summary foreach s in b1 b2 A1 A2 x b C S idx n h l lowerbound upperbound /// replace_value U hit r { cap mata: mata drop `s' // clean up } end program define _estat version 15.1 // assert last result cap conf mat e(z) if _rc { di as err "results of lppinvl2 not found" exit 301 } // syntax syntax [anything], [*] clspl2 estat `=ustrregexrf(`"`anything'"',"^estat","",1)', `options' end * Mata code *************** version 15.1 loc RS real scalar loc RV real colvector loc RM real matrix mata: mata set matastrict on `RM' _uniqrows_tagindex(`RM' X, `TM' idx) { /* -- M = `RM' _uniqrows_tagindex(M, tmp) - - unique rows with first-occurrence indices returns: (uniqrows(X)) (k x p) */ `RS' r `RM' U `RV' hit U = uniqrows(X) idx = J(0, 1, .) for (r = 1; r <= rows(U); r++) { hit = selectindex(rowsum(X :!= U[r,.]) :== 0) idx = idx \ hit[1] } return(U) } end