scikit-learn
diff --git a/‎sklearn/linear_model/_glm/_newton_solver.py
Lines changed: 7 additions & 2 deletions b/‎sklearn/linear_model/_glm/_newton_solver.py
Lines changed: 7 additions & 2 deletions
diff --git a/‎sklearn/linear_model/_glm/tests/test_glm.py
Lines changed: 5 additions & 5 deletions b/‎sklearn/linear_model/_glm/tests/test_glm.py
Lines changed: 5 additions & 5 deletions
@@ -974,6 +974,9 @@ def inner_solve(self, X, y, sample_weight):
             atol=eta * norm_G / (self.A_norm * self.r_norm),
             btol=self.tol,
             maxiter=max(n_samples, n_features) * n_classes,  # default is min(A.shape)
+            # default conlim = 1e8, for compatible systems 1e12 is still reasonable,
+            # see LSMR documentation
+            conlim=1e12,
             show=self.verbose >= 3,
         )
         # We store the estimated Frobenius norm of A and norm of residual r in
@@ -988,6 +991,8 @@ def inner_solve(self, X, y, sample_weight):
             conda,
             normx,
         ) = result
+        if self.verbose >= 2:
+            print(f"  Inner iterations in LSMR = {itn}")
         if self.coef.dtype == np.float32:
             self.coef_newton = self.coef_newton.astype(np.float32)
         if not self.linear_loss.base_loss.is_multiclass:
@@ -1010,7 +1015,7 @@ def inner_solve(self, X, y, sample_weight):
         if self.iteration == 1:
             return
         # Note: We could detect too large steps by comparing norm(coef_newton) = normx
-        # with norm(gradient) o with the already available condition number of A, e.g.
+        # with norm(gradient) or with the already available condition number of A, e.g.
         # conda.
         if istop == 7:
             self.use_fallback_lbfgs_solve = True
@@ -1033,7 +1038,7 @@ def inner_solve(self, X, y, sample_weight):
                 msg
                 + "It will now resort to lbfgs instead.\n"
                 "This may be caused by singular or very ill-conditioned Hessian "
-                " matrix. "
+                "matrix. "
                 "Further options are to use another solver or to avoid such situation "
                 "in the first place. Possible remedies are removing collinear features"
                 "of X or increasing the penalization strengths.",
 
@@ -1093,8 +1093,8 @@ def test_solver_on_ill_conditioned_X(
         np.exp(X_orig @ np.ones(X_orig.shape[1])), size=X_orig.shape[0]
     ).astype(np.float64)
 
-    tol = 1e-7
-    model = PoissonRegressor(solver=solver, alpha=0.0, tol=tol)
+    tol = 1e-8
+    model = PoissonRegressor(solver=solver, alpha=0.0, tol=tol, max_iter=200)
 
     # No warning raised on well-conditioned design, even without regularization.
     with warnings.catch_warnings():
@@ -1122,8 +1122,8 @@ def test_solver_on_ill_conditioned_X(
 
     # Construct another ill-conditioned problem by scaling of features.
     X_ill_conditioned = X_orig.copy()
-    X_ill_conditioned[:, 0] *= 1e-4
-    X_ill_conditioned[:, 1] *= 1e4
+    X_ill_conditioned[:, 0] *= 1e-6
+    X_ill_conditioned[:, 1] *= 1e2  # too large X may overflow in link function
     # Make sure that it is ill conditioned >=> large condition number.
     assert np.linalg.cond(X_ill_conditioned) > 1e5 * np.linalg.cond(X_orig)
 
@@ -1144,7 +1144,7 @@ def test_solver_on_ill_conditioned_X(
     if test_loss:
         # Without penalty, scaling of columns has no effect on predictions.
         ill_cond_deviance = mean_poisson_deviance(y, reg.predict(X_ill_conditioned))
-        if solver in ("lbfgs", "newton-cholesky", "newton-lsmr"):
+        if solver in ("lbfgs", "newton-cholesky"):
             pytest.xfail(
                 f"Solver {solver} does not converge but does so without warning."
             )