ENH: add singularity warnings

pyblp/economies/economy.py

@@ -10,9 +10,9 @@
 from ..configurations.formulation import Formulation, Absorb
 from ..configurations.iteration import Iteration
 from ..primitives import Container
-from ..utilities.algebra import precisely_identify_collinearity, precisely_identify_psd
+from ..utilities.algebra import precisely_identify_collinearity, precisely_identify_singularity, precisely_identify_psd
 from ..utilities.basics import (
-    Array, Bounds, Error, Groups, RecArray, StringRepresentation, format_table, get_indices, output, warn
+    Array, Bounds, Error, Groups, RecArray, StringRepresentation, format_number, format_table, get_indices, output, warn
 )
 
 
@@ -202,8 +202,22 @@ def _detect_collinearity(self, added_exogenous: bool) -> None:
                 )
 
     @staticmethod
-    def _detect_psd(matrix: Array, name: str) -> None:
-        """Detect whether a matrix is PSD."""
+    def _detect_singularity(matrix: Array, name: str) -> None:
+        """Detect any singularity issues in a matrix."""
+        singular, successful, condition = precisely_identify_singularity(matrix)
+        common_message = "To disable singularity checks, set options.singular_tol = numpy.inf."
+        if not successful:
+            warn(f"Failed to compute the condition number of {name} while checking for singularity. {common_message}")
+        if singular:
+            prefix = "nearly " if condition < np.inf else ""
+            warn(
+                f"Detected that {name} is {prefix}singular with condition number {format_number(condition).strip()}. "
+                f"{common_message}"
+            )
+
+    @staticmethod
+    def _require_psd(matrix: Array, name: str) -> None:
+        """Require that a matrix is PSD."""
         psd, successful = precisely_identify_psd(matrix)
         common_message = "To disable PSD checks, set options.psd_atol = options.psd_rtol = numpy.inf."
         if not successful:

pyblp/economies/problem.py

@@ -490,7 +490,8 @@ def solve(
                 micro_moment_covariances = np.c_[np.asarray(micro_moment_covariances, options.dtype)]
                 if micro_moment_covariances.shape != (moments.MM, moments.MM):
                     raise ValueError(f"extra_micro_moments must be a square {moments.MM} by {moments.MM} matrix.")
-                self._detect_psd(micro_moment_covariances, "micro_moment_covariances")
+                self._require_psd(micro_moment_covariances, "micro_moment_covariances")
+                self._detect_singularity(micro_moment_covariances, "micro_moment_covariances")
 
         # choose whether to do an initial update
         if initial_update is None:
@@ -521,12 +522,15 @@ def solve(
             M = self.MD + self.MS + moments.MM
             if W.shape != (M, M):
                 raise ValueError(f"W must be a square {M} by {M} matrix.")
-            self._detect_psd(W, "W")
+            self._require_psd(W, "W")
+            self._detect_singularity(W, "W")
         else:
-            W, successful = precisely_invert(scipy.linalg.block_diag(
+            S = scipy.linalg.block_diag(
                 self.products.ZD.T @ self.products.ZD / self.N,
                 self.products.ZS.T @ self.products.ZS / self.N,
-            ))
+            )
+            self._detect_singularity(S, "the 2SLS weighting matrix")
+            W, successful = precisely_invert(S)
             if not successful:
                 raise ValueError("Failed to compute the 2SLS weighting matrix. There may be instrument collinearity.")
 

pyblp/economies/simulation.py

@@ -503,7 +503,7 @@ def __init__(
         if self.xi is None and self.omega is None and self.K3 > 0:
             covariance = correlation * np.sqrt(xi_variance * omega_variance)
             covariances = np.array([[xi_variance, covariance], [covariance, omega_variance]], options.dtype)
-            self._detect_psd(covariances, "the covariance matrix from xi_variance, omega_variance, and correlation")
+            self._require_psd(covariances, "the covariance matrix from xi_variance, omega_variance, and correlation")
             xi_and_omega = state.multivariate_normal([0, 0], covariances, self.N, check_valid='ignore')
             self.xi = xi_and_omega[:, [0]].astype(options.dtype)
             self.omega = xi_and_omega[:, [1]].astype(options.dtype)

pyblp/options.py

@@ -64,6 +64,14 @@
     To always attempt to compute classic inverses first, set ``pyblp.options.pseudo_inverses = False``. If a classic
     inverse cannot be computed, an error will be displayed, and a pseudo-inverse may be computed instead.
 
+weights_tol : `float`
+    Tolerance for detecting integration weights that do not sum to one, which is by default ``1e-10``. In most setups
+    weights should essentially sum to one, but for example with importance sampling they may be slightly different.
+    Warnings can be disabled by setting this to ``numpy.inf``.
+singular_tol : `float`
+    Tolerance for detecting singular matrices, which is by default ``1 / numpy.finfo(options.dtype).eps``. If a matrix
+    has a condition number larger than this tolerance, a warning will be displayed. To disable singularity checks, set
+    ``pyblp.options.singular_tol = numpy.inf``.
 collinear_atol : `float`
     Absolute tolerance for detecting collinear columns in each matrix of product characteristics and instruments:
     :math:`X_1`, :math:`X_2`, :math:`X_3`, :math:`Z_D`, and :math:`Z_S`.
@@ -75,10 +83,6 @@
     The default absolute tolerance is ``1e-14``. To disable collinearity checks, set
     ``pyblp.options.collinear_atol = pyblp.options.collinear_rtol = 0``.
 
-weights_tol : `float`
-    Tolerance for detecting integration weights that do not sum to one, which is by default ``1e-10``. In most setups
-    weights should essentially sum to one, but for example with importance sampling they may be slightly different.
-    Warnings can be disabled by setting this to ``numpy.inf``.
 collinear_rtol : `float`
     Relative tolerance for detecting collinear columns, which is by default also ``1e-14``.
 psd_atol : `float`
@@ -109,5 +113,6 @@
 finite_differences_epsilon = _np.sqrt(_np.finfo(dtype).eps)
 pseudo_inverses = True
 weights_tol = 1e-10
+singular_tol = 1 / _np.finfo(dtype).eps
 collinear_atol = collinear_rtol = 1e-14
 psd_atol = psd_rtol = 1e-8

pyblp/results/problem_results.py

@@ -418,6 +418,7 @@ def __init__(
 
             # update the weighting matrix
             S_for_weights = self._compute_S(progress.moments, W_type, center_moments)
+            self.problem._detect_singularity(S_for_weights, "the estimated covariance matrix of GMM moments")
             self.updated_W, W_errors = compute_gmm_weights(S_for_weights)
             self._errors.extend(W_errors)
 
@@ -430,6 +431,10 @@ def __init__(
                 S_for_covariances = S_for_weights
                 if se_type != W_type or center_moments:
                     S_for_covariances = self._compute_S(progress.moments, se_type)
+                    self.problem._detect_singularity(
+                        S_for_covariances,
+                        "the estimated covariance matrix of GMM moments used for computing standard errors",
+                    )
 
                 # if this is the first step, an unadjusted weighting matrix needs to be used when computing unadjusted
                 #   covariances so that they are scaled properly

pyblp/utilities/algebra.py

@@ -21,6 +21,19 @@ def compute_condition_number(x: Array) -> float:
         return np.nan
 
 
+def precisely_identify_singularity(x: Array) -> Tuple[Array, bool, Array]:
+    """Compute the condition number of a matrix to identify whether it is nearly singular."""
+    singular = False
+    successful = True
+    condition = np.nan
+    if options.singular_tol < np.inf:
+        condition = compute_condition_number(x)
+        successful = not np.isnan(condition)
+        singular = successful and condition > options.singular_tol
+
+    return singular, successful, condition
+
+
 def precisely_identify_collinearity(x: Array) -> Tuple[Array, bool]:
     """Compute the QR decomposition of a matrix and identify which diagonal elements of the upper diagonal matrix are
     within absolute and relative tolerances.