Repository navigation
Conversation
Add the keyword solver to IV2SLS, IVLIML and their from_formula. solver="qr" tests the rank of [exog, endog] and [exog, instruments] after scaling the columns to unit length, with a relative tolerance of 1e-9, and computes the k-class estimates, their covariance and the kappa of LIML from QR decompositions, as R's ivreg does. The results do not depend on the units of the variables, and their rounding error grows with the condition number of the scaled data, not with its square. The default is solver="legacy", which computes as before. If solver is not given, a FutureWarning is issued for data that solver="qr" finds rank deficient, and the error for columns that are only rank deficient because of their units suggests solver="qr". linearmodels 8.0 will use "qr" by default. The covariance estimators, proj/annihilate and the Cragg-Donald and Kleibergen-Paap statistics take the same keyword, and the first stage results and the tests of IVResults use the solver of the model. Builds on the unit-length rank check of #850. Fixes case 2 of #823 (kappa depends on the units with 2 or more endogenous regressors) with solver="qr".
Rename solver="legacy" to solver="pinv", so that both values name the numerical method. The default, None, still computes as "pinv" and warns.
The covariance estimators accepted solver positionally after extra_df, and the private _first_stage_rank and _rank did too. Every function that takes solver now requires it as a keyword, and a test checks the public ones.
…heck The default rank check computed the singular values of the columns scaled to unit length for every model, to decide whether to warn, which made creating a model up to twice as slow. The rank is now computed from the singular values of the columns as they are, as matrix_rank computes it, and the scaled columns are only tested if the condition number is at least 1e9 / (2 sqrt(ncol)). Scaling the columns cannot increase the condition number by more than sqrt(ncol) (van der Sluis), so the decision does not change. With solver="qr", fit decomposes the instruments once for kappa, the estimate and the covariance, OLS (kappa = 0) does not decompose the instruments, and the covariance uses one QR decomposition of x when z is x.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #870 +/- ##
========================================
Coverage 99.68% 99.69%
========================================
Files 201 202 +1
Lines 38176 38950 +774
Branches 2442 2521 +79
========================================
+ Hits 38057 38831 +774
Misses 51 51
Partials 68 68
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
solver="pinv" (bit-for-bit the same results): - fit computes pinv(z) once for kappa, the estimate and the covariance, which each computed it before, and kappa reuses it for x1 if x1 is z - the covariance reuses pinv(z) @ x, x'z pinv(z) x and V - the terms that are multiplied by exactly 0 when kappa is 0 or 1, such as the x'x of 2SLS, are not computed (the inputs are finite) - the kernel covariance computes the projected regressors once solver="qr": - the covariance reuses the factors of the estimate, and computes U from R instead of decomposing x again for OLS - kappa < 1 applies the reflections of the QR of x to y without forming Q - the rank check and kappa use QR without forming Q, and proj uses QR and only falls back to the SVD if the scaled columns are rank deficient
Sharing pinv(z) between kappa, the estimate and the covariance kept it alive while kappa decomposed x1, which raised the peak memory of fit by up to 19%. Kappa now computes the residuals on x1 first, so that their decomposition is freed before the one of z is computed and kept. The residuals are independent, so the results are the same.
With kappa other than 1, the k-class estimator fits a residual of y - Xb, whose error grows with the square of the condition number. The Hadamard data only gave exact estimates with kappa = 1.25 because the rounding error of their products cancelled in the order of the sums of OpenBLAS on x86. It did not on Windows ARM64 (6.4e-7 and 1.0e-2), and does not on x86 with the rows permuted. Fixed kappa is now tested on data without the error, where the error grows with the condition number, with the rows in both orders. The tolerance of the Fuller test, which has a residual, follows its bound.
With solver="pinv", kappa only uses the cross product of the residuals on x1, which is now computed before z is decomposed, so that the residuals are freed. This removes the last 2-4% of extra peak memory of fit compared with main; the results are the same. Also narrow the optional basis of kclass_factors for mypy.
The keyword of the numerical method of IV2SLS and IVLIML is now method, as in statsmodels, with the values "pinv" and "qr". The internal attribute is _numerical_method, since _method holds the name of the estimator. If method is not given, a FutureWarning says that the default changes from "pinv" to "qr" in linearmodels 8.0, that method="qr" gives the new method now and that method="pinv" keeps the current one without the warning. For data that method="qr" finds rank deficient, the warning also says so. Only the models that the user creates warn: the models of the diagnostics and the covariance estimators get the method of the model. The test suite ignores the warning in pyproject.toml; the tests of the warning use pytest.warns. The kappa tests that turn every warning into an error ignore it too.
The IV covariance estimators, estimate_parameters, cragg_donald, cragg_donald_f, kleibergen_paap, kleibergen_paap_f and FirstStageResults default to method=None, which is "pinv". If they are called directly from outside linearmodels, a FutureWarning says that the default will be "qr" in linearmodels 8.0. They do not warn when linearmodels calls them: the models pass their own method, and AbsorbingLS and IVGMM use the default without a warning. The tests of linearmodels count as callers outside linearmodels. The docstrings of the models and of fit explain that with "pinv" the estimates lose digits well below the condition number of 1e9 at which "qr" rejects the data, about 16 - 2 log10(c) correct digits against 16 - log10(c).
The tests pass method="pinv" wherever they create IV2SLS or IVLIML or call a lower-level function that takes method, except the tests of the warning, which use pytest.warns. Where the estimator is a parameter that can also be IVGMM, IVGMMCUE, IVJIVE or a panel model, the new test helper method_kwargs adds method only for those that take it. pyproject.toml no longer ignores the warning, and the kappa tests turn every warning into an error again. The suite passes with the warning promoted to an error.
The IV examples pass method="qr", so that they do not show the warning, and the basic examples have a new cell on the two numerical methods: what each does, the trade-offs (accuracy and units, the stricter rank check, the exact results of 7.1 with "pinv") and why "qr" is recommended.
The 7.2 notes explain the digits that "pinv" loses below the tolerance of "qr", why IVGMM, IVGMMCUE and IVJIVE do not have method, and that the lower-level functions warn only when they are called directly.
Codecov found added lines that no test runs. Those that no input can reach, or that change nothing, are removed: - _rank_and_condition_number: the fallback to matrix_rank. NumPy 1 computes the tolerance as (smax * m) * eps and NumPy 2 as smax * (m * eps), which are the same number since eps is a power of 2, unless smax * m overflows. - _estimate_kappa: the checks of the keys of shared, which is always empty when kappa is estimated. - kclass_solve: the check that the factors have t, which every caller computes, and the guard for no regressors, since solve_triangular without trans handles empty arrays. - proj: the separate branch for fewer rows than columns. The QR decomposition then has a square Q, which is the projection if x has full row rank. - The OLS shortcut of the covariance estimators tests the shape in array_equal. The guards for no regressors in kclass_factors and in the covariance are kept, with comments: qr_multiply fails if x has no columns, and solve_triangular with trans="T" fails with empty arrays in SciPy 2.0.dev.
- The rank tests with a constant and a full set of dummies as exog, which proj projects on with the SVD. "qr" agrees with pinv, and Kleibergen-Paap does not change without the redundant dummy. Without the SVD, Cragg-Donald was 307.7 instead of 303.9. - A model with fewer observations than regressors is rank deficient with both methods. - estimate_parameters and the covariance estimators with no regressors return empty arrays with both methods. - The covariance estimators agree when the model is just identified, so that z has the shape of x, which tests that the OLS shortcut is only used if z is x. - proj agrees with pinv for full rank, rank deficient and wide x. method_kwargs no longer catches errors of inspect.signature, which cannot fail for the IV models and functions.
bc96124 removed it after checking SciPy 1.18 and 2.0.dev, where solve_triangular returns an empty array. SciPy 1.13 and earlier raise LinAlgError for empty arrays, which made test_no_regressors fail in the minimums (SciPy 1.8), NumPy 1.26 (1.11) and Cython coverage (1.13) jobs. Checked with SciPy 1.13.1 and 1.11.4, where all 255 tests of test_numerical_method.py now pass.
bashtage
marked this pull request as ready for review
November 7, 2026 10:31
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Adds an opt-in QR method to
IV2SLSandIVLIMLfor 7.2. It supersedes #850, which is closed. The default computation does not change in 7.2; aFutureWarningannounces that the default becomes"qr"in 8.0.What changes
IV2SLS,IVLIMLand theirfrom_formulatake a keyword-onlymethod, as in statsmodels, with values"pinv"and"qr"(not case sensitive).method="qr":[exog, endog]and[exog, instruments]after scaling the columns to unit length. They are rank deficient if the condition number is at least 1e9 (relative tolerance 1e-9).ivreg. For kappa < 1 it is a stacked least-squares problem. For kappa > 1 (LIML) it usesR'GRwithG = I - (kappa-1)F'F. The cross products are never formed.R⁻¹, with scores in orthonormal coordinates, so robust, clustered and kernel SEs keep the cond·eps accuracy.σ_min(M₁W R⁻¹)²). This removes case 2 of IV: LIML kappa with an endogenous regressor spanned by the instruments, kappa's unit dependence, scale-dependent instrument rank check #823: the old kappa changed with the units when there were 2 or more endogenous regressors.method="pinv"gives today's results bit for bit: rawmatrix_rank,pinv, inverted cross products and the eigenvalue kappa. It no longer repeats work within a fit:pinv(z)is computed once instead of three times, and the covariance reuses the estimator's products. Terms multiplied by exactly zero at kappa 0 or 1 (such as 2SLS'sx'x) are skipped. That is exact because the rank check rejects non-finite data. With these changes the default path takes about half the time of main.None, behaves as"pinv", plus:FutureWarningwhenevermethodis not given. It says that the default changes from"pinv"to"qr"in linearmodels 8.0, thatmethod="qr"gives the new method now, and thatmethod="pinv"keeps the current one and silences the warning. For data that"qr"would refuse (the cond ≥ 1e9 band, where"pinv"estimates have few or no correct digits), the warning also says so. It points at the caller's line, also throughfrom_formulaandlinearmodels.formula.iv_2sls/iv_liml. Models built internally by the diagnostics pass the model's method and don't warn.method="qr".estimate_parameters,cragg_donald(_f),kleibergen_paap(_f)andFirstStageResultstakemethodtoo, with defaultNone(="pinv"). Called directly from outside linearmodels, they issue the same kind ofFutureWarning. A frame check stays silent when linearmodels itself calls them: the models pass their own method, and AbsorbingLS and IVGMM use the default. The first-stage results and the tests ofIVResults(Sargan, Basmann, Wooldridge, Durbin, Wu-Hausman, Anderson-Rubin, C) use the model's method.fitdocstrings and the 7.2 notes explain the accuracy below the 1e9 tolerance. With"pinv"about 16 − 2·log10(c) digits are correct, e.g. 6 at c = 1e5 and none at 1e8; with"qr"about 16 − log10(c).IVGMM,IVGMMCUEandIVJIVEdon't havemethod.method="qr"(they run without warnings), andiv_basic-exampleshas a new cell on the two methods, their trade-offs, and why"qr"is recommended.IVGMM,IVGMMCUE,IVJIVEandAbsorbingLSkeep today's unscaled rank check and computation. GMM inverts Z'Z and needs its own follow-up.Evidence
Exact references (tests): the k-class estimates and all four 2SLS covariances are exact on Hadamard-built data, with the references computed with
fractions. Fixed kappas use data without an error term, since a residual adds a cond² term for any method. Each case also runs with permuted rows, so that the result can't depend on the BLAS summation order."qr"params"qr"SEs"pinv"paramsLinAlgErrorScale invariance (tests): rescaling
y, an exog, an endog or an instrument by 2^±60 leaves the"qr"params and SEs bit-identical after mapping back, for all four covariance types and for 2SLS, LIML, Fuller and kappa=0.5. So are all the diagnostics listed above, kappa included. The one exception is the automatic kernel bandwidth, which sums score columns and so depends on units by definition, with"pinv"too.Simulated data, 7,560 data sets from the BUG: Test the rank of the IV inputs after scaling the columns #850 study. Scored against a 120-digit reference at tolerance 1e-6. A refusal is "needed" if the exact estimate is undefined or QR can't get within 1e-6. The grid is deliberately heavy in collinear and extreme-scale data, so compare the rows rather than reading the shares as frequencies.
"qr"The
"qr"row is the same for the 2SLS params, all four SE types and LIML; the main row varies by at most 0.2 points. Median error with collinear regressors:"qr"params"qr"robust SEDefault path = main, bit for bit: 10,484 numbers from 477 fits compare equal with exact equality. The fits cover 2SLS, LIML, Fuller, kappa 0/0.5, OLS, IVGMM and AbsorbingLS; unadjusted, robust, debiased, one- and two-way clustered, and automatic and fixed kernel covariances; weighted and unweighted data; scaled, collinear and just-identified designs; and every test statistic and the first-stage results. The same holds on all 7,560 simulated data sets.
Large models: seconds (best of 1 to 3 runs), as main → branch default /
"qr", with 1 BLAS thread. With 12 threads the ratios are similar: default 0.43–0.69× main,"qr"0.33–0.59×. The"qr"estimates are within 1.2e-13 of main's on these well-conditioned data.pytest -m "not slow" linearmodelsgives 23,350 passed, 2,473 skipped and 47 xfailed. That is main's 23,095 plus the 255 new tests inlinearmodels/tests/iv/test_numerical_method.py, with the same skips and the same 992 warnings (all from pyhdfe's scipy.sparse use). The tests passmethodexplicitly everywhere except in the tests of the warning, and the suite also passes with the default-method warning promoted to an error. Where a test's estimator is a parameter that can also be IVGMM, JIVE or a panel model, the helperlinearmodels.tests._utility.method_kwargsaddsmethodonly for those that take it.proj), fewer observations than regressors, and no regressors inestimate_parametersand the covariance estimators. Branches that no input reaches were removed.Decided
method(as in statsmodels), with values"pinv"and"qr", and both methods stay for now.methodis not given. The lower-level functions warn only when they are called directly.methodexplicitly.Builds on the unit-length rank check proposed by inchang-ing in #850.
🤖 Generated with Claude Code