switched singulrity detecction to condition number of hessian

Author: davidsebfischer

batchglm/train/numpy/base_glm/estimator.py

@@ -265,25 +265,31 @@ def iwls_step(
         if isinstance(delta_theta, dask.array.core.Array):
             delta_theta = delta_theta.compute()
 
-        with np.errstate(all="ignore"):
-            if isinstance(a, dask.array.core.Array):
-                # Have to use a workaround to solve problems in parallel in dask here. This workaround does
-                # not work if there is only a single problem, ie. if the first dimension of a and b has length 1.
-                if a.shape[0] != 1:
-                    delta_theta[:, idx_update] = dask.array.map_blocks(
-                        np.linalg.solve, a, b[:, :, None], chunks=b[:, :, None].shape
-                    ).squeeze().T.compute()
-                else:
+        if isinstance(a, dask.array.core.Array):
+            # Have to use a workaround to solve problems in parallel in dask here. This workaround does
+            # not work if there is only a single problem, ie. if the first dimension of a and b has length 1.
+            if a.shape[0] != 1:
+                invertible = np.where(dask.array.map_blocks(
+                    np.linalg.cond, a, chunks=a.shape
+                ).compute() < 1 / sys.float_info.epsilon)[0]
+                delta_theta[:, idx_update[invertible]] = dask.array.map_blocks(
+                    np.linalg.solve, a[invertible], b[invertible, :, None],
+                    chunks=b[invertible, :, None].shape
+                ).squeeze().T.compute()
+            else:
+                if np.linalg.cond(a.compute(), p=None) < 1 / sys.float_info.epsilon:
                     delta_theta[:, idx_update] = np.expand_dims(
                         np.linalg.solve(a[0], b[0]).compute(),
                         axis=-1
                     )
-            else:
-                delta_theta[:, idx_update] = np.linalg.solve(a, b).T
-        linalg_errors = np.isnan(delta_theta[0, :])
-        if np.any(linalg_errors):
-            print("caught %i linalg errors" % np.sum(linalg_errors))
-            delta_theta[:, linalg_errors] = 0.
+                    invertible = np.array([0])
+                else:
+                    invertible = np.array([])
+        else:
+            invertible = np.where(np.linalg.cond(a, p=None) < 1 / sys.float_info.epsilon)[0]
+            delta_theta[:, idx_update[invertible]] = np.linalg.solve(a[invertible], b[invertible]).T
+        if invertible.shape[0] < len(idx_update):
+            print("caught %i linalg singular matrix errors" % len(idx_update) - invertible.shape[0])
         # Via np.linalg.lsts:
         #delta_theta[:, idx_update] = np.concatenate([
         #    np.expand_dims(np.linalg.lstsq(a[i, :, :], b[i, :])[0], axis=-1)
@@ -516,8 +522,10 @@ def finalize(self):
         """
         # Read from numpy-IRLS estimator specific model:
         self._hessian = - self.model.fim.compute()
-        with np.errstate(all="ignore"):
-            self._fisher_inv = np.linalg.inv(- self._hessian)
+        fisher_inv = np.zeros_like(self._hessian)
+        invertible = np.where(np.linalg.cond(self._hessian, p=None) < 1 / sys.float_info.epsilon)[0]
+        fisher_inv[invertible] = np.linalg.inv(- self._hessian[invertible])
+        self._fisher_inv = fisher_inv
         self._jacobian = np.sum(np.abs(self.model.jac.compute() / self.model.x.shape[0]), axis=1)
         self._log_likelihood = self.model.ll_byfeature.compute()
         self._loss = np.sum(self._log_likelihood)