BUG: Do not materialize weighted data when all weights are 1

linearmodels/iv/model.py

@@ -282,13 +282,23 @@ def __init__(
         self._sqrt_weights = w = sqrt(self.weights.ndarray)
         # dependent variable
         self._y = self.dependent.ndarray
-        self._wy = self._y * w
         # model regressors
         self._x = c_[self.exog.ndarray, self.endog.ndarray]
-        self._wx = self._x * w
         # first-stage regressors
         self._z = c_[self.exog.ndarray, self.instruments.ndarray]
-        self._wz = self._z * w
+        if npall(w == 1.0):
+            # All of the weights are 1, as when weights is None or constant,
+            # and so scaling by the square root of the weights does not
+            # change the data. The weighted data is the data, and avoiding
+            # the copies substantially reduces the memory used when the
+            # data is large (GH 607).
+            self._wy = self._y
+            self._wx = self._x
+            self._wz = self._z
+        else:
+            self._wy = self._y * w
+            self._wx = self._x * w
+            self._wz = self._z * w
 
         self._has_constant = False
         self._regressor_is_exog = array(

linearmodels/tests/iv/test_model.py

@@ -17,6 +17,7 @@
 from linearmodels.iv import IV2SLS, IVGMM, IVGMMCUE, IVLIML
 from linearmodels.iv.model import _OLS
 from linearmodels.iv.results import compare
+from linearmodels.shared.exceptions import MissingValueWarning
 from linearmodels.shared.hypotheses import WaldTestStatistic
 from linearmodels.shared.utility import AttrDict
 
@@ -508,3 +509,202 @@ def test_first_state_f(data):
         indiv = _OLS(endog[:, i], rhs).fit(cov_type="unadjusted")
         stat = indiv.wald_test(null_rest, null_value)
         assert_allclose(diag.iloc[i].loc["f.stat"], stat.stat / nz)
+
+
+def test_weights_all_one_not_materialized(data):
+    # GH 607: weighted copies of the data are not stored when all of the
+    # weights are 1, which is the case when weights is None or constant
+    mod = IV2SLS(data.dep, data.exog, data.endog, data.instr)
+    assert mod._wy is mod._y
+    assert mod._wx is mod._x
+    assert mod._wz is mod._z
+
+    weights = np.full(data.nobs, 2.5)
+    mod = IV2SLS(data.dep, data.exog, data.endog, data.instr, weights=weights)
+    assert mod._wy is mod._y
+    assert mod._wx is mod._x
+    assert mod._wz is mod._z
+
+
+def test_weights_nonconstant_materialized(data):
+    rs = np.random.RandomState(12345)
+    weights = 0.5 + rs.random_sample(data.nobs)
+    mod = IV2SLS(data.dep, data.exog, data.endog, data.instr, weights=weights)
+    assert mod._wz is not mod._z
+    root_w = np.sqrt(mod.weights.ndarray)
+    assert_allclose(mod._wy, mod._y * root_w)
+    assert_allclose(mod._wx, mod._x * root_w)
+    assert_allclose(mod._wz, mod._z * root_w)
+
+
+def test_weights_all_one_results_unchanged(data):
+    # Estimation with no weights or constant weights is identical to
+    # explicitly weighting by the square root of the normalized weights
+    res = IV2SLS(data.dep, data.exog, data.endog, data.instr).fit()
+    weights = np.full(data.nobs, 3.0)
+    res_const = IV2SLS(
+        data.dep, data.exog, data.endog, data.instr, weights=weights
+    ).fit()
+    assert_allclose(res_const.params, res.params, rtol=0, atol=0)
+    assert_allclose(res_const.cov, res.cov, rtol=0, atol=0)
+
+
+def test_weights_float32_not_upcast_by_weighting():
+    # GH 607: float32 data with float32 weights produces the same estimates
+    # as the float64 data, up to the precision of the float32 inputs, and
+    # constant float32 weights do not materialize weighted copies
+    rs = np.random.RandomState(12345)
+    n = 500
+    dep = rs.standard_normal((n, 1))
+    exog = rs.standard_normal((n, 2))
+    endog = rs.standard_normal((n, 1))
+    instr = rs.standard_normal((n, 2))
+    weights = 0.5 + rs.random_sample(n)
+
+    res64 = IV2SLS(dep, exog, endog, instr, weights=weights).fit()
+    res32 = IV2SLS(
+        dep.astype(np.float32),
+        exog.astype(np.float32),
+        endog.astype(np.float32),
+        instr.astype(np.float32),
+        weights=weights.astype(np.float32),
+    ).fit()
+    assert_allclose(res32.params, res64.params, rtol=1e-5, atol=1e-6)
+
+    mod32 = IV2SLS(
+        dep.astype(np.float32),
+        exog.astype(np.float32),
+        endog.astype(np.float32),
+        instr.astype(np.float32),
+        weights=np.ones(n, dtype=np.float32),
+    )
+    assert mod32._wz is mod32._z
+
+
+def test_weights_float32_all_one_fit():
+    # GH 607: float32 data with all-one float32 weights are fitted through
+    # the unmaterialized path, and the estimates agree with the estimates
+    # of the same data stored as float64
+    rs = np.random.RandomState(12345)
+    n = 500
+    dep = rs.standard_normal((n, 1)).astype(np.float32)
+    exog = rs.standard_normal((n, 2)).astype(np.float32)
+    endog = rs.standard_normal((n, 1)).astype(np.float32)
+    instr = rs.standard_normal((n, 2)).astype(np.float32)
+    weights = np.ones(n, dtype=np.float32)
+
+    mod = IV2SLS(dep, exog, endog, instr, weights=weights)
+    assert mod._wy is mod._y
+    assert mod._wx is mod._x
+    assert mod._wz is mod._z
+    res = mod.fit()
+    res64 = IV2SLS(
+        dep.astype(np.float64),
+        exog.astype(np.float64),
+        endog.astype(np.float64),
+        instr.astype(np.float64),
+    ).fit()
+    assert_allclose(res.params, res64.params, rtol=1e-6, atol=1e-8)
+    assert_allclose(res.cov, res64.cov, rtol=1e-5, atol=1e-9)
+
+
+def test_weights_integer_inputs_fit():
+    # GH 607: integer data with constant integer weights are fitted through
+    # the unmaterialized path, and the estimates are identical to the
+    # estimates of the same data stored as float64
+    rs = np.random.RandomState(12345)
+    n = 500
+    dep = np.round(10 * rs.standard_normal((n, 1))).astype(np.int64)
+    exog = np.round(10 * rs.standard_normal((n, 2))).astype(np.int64)
+    endog = np.round(10 * rs.standard_normal((n, 1))).astype(np.int64)
+    instr = np.round(10 * rs.standard_normal((n, 2))).astype(np.int64)
+    weights = np.full(n, 2)
+
+    mod = IV2SLS(dep, exog, endog, instr, weights=weights)
+    assert mod._wy is mod._y
+    assert mod._wx is mod._x
+    assert mod._wz is mod._z
+    res = mod.fit()
+    res64 = IV2SLS(
+        dep.astype(np.float64),
+        exog.astype(np.float64),
+        endog.astype(np.float64),
+        instr.astype(np.float64),
+        weights=np.full(n, 2.0),
+    ).fit()
+    assert_allclose(res.params, res64.params, rtol=0, atol=0)
+    assert_allclose(res.cov, res64.cov, rtol=0, atol=0)
+
+
+def test_weights_all_one_matches_materialized_reference(data):
+    # GH 607: the estimates from the unmaterialized path agree with a
+    # reference computed directly from explicitly materialized weighted
+    # data, which is what the model stored before the copies were skipped
+    y = np.asarray(data.dep, dtype=np.float64)
+    x = np.c_[
+        np.asarray(data.exog, dtype=np.float64),
+        np.asarray(data.endog, dtype=np.float64),
+    ]
+    z = np.c_[
+        np.asarray(data.exog, dtype=np.float64),
+        np.asarray(data.instr, dtype=np.float64),
+    ]
+    n = y.shape[0]
+    root_w = np.ones((n, 1))
+    wy, wx, wz = y * root_w, x * root_w, z * root_w
+    pz = wz @ pinv(wz.T @ wz) @ wz.T
+    x_pz_x_inv = pinv(wx.T @ pz @ wx)
+    params_ref = x_pz_x_inv @ (wx.T @ pz @ wy)
+    resids_ref = wy - wx @ params_ref
+    s2 = float(np.squeeze(resids_ref.T @ resids_ref)) / n
+    cov_ref = s2 * x_pz_x_inv
+
+    res = IV2SLS(data.dep, data.exog, data.endog, data.instr).fit(
+        cov_type="unadjusted", debiased=False
+    )
+    assert_allclose(res.params, params_ref.ravel(), rtol=1e-8, atol=1e-10)
+    assert_allclose(
+        np.asarray(res.resids).ravel(), resids_ref.ravel(), rtol=1e-8, atol=1e-10
+    )
+    assert_allclose(res.cov, cov_ref, rtol=1e-8, atol=1e-12)
+
+
+def test_weights_single_non_one_materialized(data):
+    # GH 607: a single weight that is not 1 takes the materializing path,
+    # and the stored data is scaled by the square root of the weights
+    weights = np.ones(data.nobs)
+    weights[0] = 2.0
+    mod = IV2SLS(data.dep, data.exog, data.endog, data.instr, weights=weights)
+    assert mod._wy is not mod._y
+    assert mod._wx is not mod._x
+    assert mod._wz is not mod._z
+    root_w = np.sqrt(mod.weights.ndarray)
+    assert_allclose(mod._wy, mod._y * root_w)
+    assert_allclose(mod._wx, mod._x * root_w)
+    assert_allclose(mod._wz, mod._z * root_w)
+
+
+def test_weights_nan_weight_drops_observation():
+    # GH 607: a NaN weight is a missing value, so the observation is
+    # dropped, as it was before the weighted copies were skipped, and the
+    # remaining all-one weights use the unmaterialized path
+    rs = np.random.RandomState(12345)
+    n = 200
+    dep = rs.standard_normal((n, 1))
+    exog = np.column_stack([np.ones(n), rs.standard_normal(n)])
+    endog = rs.standard_normal((n, 1))
+    instr = rs.standard_normal((n, 2))
+    weights = np.ones(n)
+    weights[3] = np.nan
+
+    with pytest.warns(MissingValueWarning):
+        mod = IV2SLS(dep, exog, endog, instr, weights=weights)
+    res = mod.fit()
+    keep = np.ones(n, dtype=bool)
+    keep[3] = False
+    res_ref = IV2SLS(dep[keep], exog[keep], endog[keep], instr[keep]).fit()
+    assert res.nobs == n - 1
+    assert_allclose(res.params, res_ref.params, rtol=0, atol=0)
+    assert mod._wy is mod._y
+    assert mod._wx is mod._x
+    assert mod._wz is mod._z