LogisticRegression fails in sklearn 1.1.1 with newton-cg solver when X only contains one predictor

### Describe the bug

I think this is related to https://github.com/scikit-learn/scikit-learn/pull/21808

As explained in the title, this bug appeared in version 1.1.1. `LogisticRegression(solver="newton-cg")` fails when `X` only contains a single predictor, unless using `fit_intercept=False`.

### Steps/Code to Reproduce

```python
import numpy as np
from sklearn.linear_model import LogisticRegression

y = np.array([1, 1, 0, 0, 1, 1, 0, 1])
X = np.array([[0.5, 0.65, 1.1, 1.25, 0.8, 0.54, 0.95, 0.7]]).T

# Newton-cg solver
LogisticRegression(solver="newton-cg").fit(X, y)  # FAILS in 1.1.1, works in previous versions
LogisticRegression(solver="newton-cg", fit_intercept=False).fit(X, y)  # Works in 1.1.1 and previous versions

# Default solver
LogisticRegression().fit(X, y)  # Works in 1.1.1 and previous versions
```

### Expected Results

No error is thrown

### Actual Results

```bash
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-1-bc267bc834b6> in <module>
      4 y = np.array([1, 1, 0, 0, 1, 1, 0, 1])
      5 X = np.array([[0.5, 0.65, 1.1, 1.25, 0.8, 0.54, 0.95, 0.7]]).T
----> 6 LogisticRegression(solver="newton-cg").fit(X, y)  # Fails
      7 LogisticRegression().fit(X, y)

~/.pyenv/versions/3.8.3/lib/python3.8/site-packages/sklearn/linear_model/_logistic.py in fit(self, X, y, sample_weight)
   1231             n_threads = 1
   1232 
-> 1233         fold_coefs_ = Parallel(n_jobs=self.n_jobs, verbose=self.verbose, prefer=prefer)(
   1234             path_func(
   1235                 X,

~/.pyenv/versions/3.8.3/lib/python3.8/site-packages/joblib/parallel.py in __call__(self, iterable)
   1039             # remaining jobs.
   1040             self._iterating = False
-> 1041             if self.dispatch_one_batch(iterator):
   1042                 self._iterating = self._original_iterator is not None
   1043 

~/.pyenv/versions/3.8.3/lib/python3.8/site-packages/joblib/parallel.py in dispatch_one_batch(self, iterator)
    857                 return False
    858             else:
--> 859                 self._dispatch(tasks)
    860                 return True
    861 

~/.pyenv/versions/3.8.3/lib/python3.8/site-packages/joblib/parallel.py in _dispatch(self, batch)
    775         with self._lock:
    776             job_idx = len(self._jobs)
--> 777             job = self._backend.apply_async(batch, callback=cb)
    778             # A job can complete so quickly than its callback is
    779             # called before we get here, causing self._jobs to

~/.pyenv/versions/3.8.3/lib/python3.8/site-packages/joblib/_parallel_backends.py in apply_async(self, func, callback)
    206     def apply_async(self, func, callback=None):
    207         """Schedule a func to be run"""
--> 208         result = ImmediateResult(func)
    209         if callback:
    210             callback(result)

~/.pyenv/versions/3.8.3/lib/python3.8/site-packages/joblib/_parallel_backends.py in __init__(self, batch)
    570         # Don't delay the application, to avoid keeping the input
    571         # arguments in memory
--> 572         self.results = batch()
    573 
    574     def get(self):

~/.pyenv/versions/3.8.3/lib/python3.8/site-packages/joblib/parallel.py in __call__(self)
    260         # change the default number of processes to -1
    261         with parallel_backend(self._backend, n_jobs=self._n_jobs):
--> 262             return [func(*args, **kwargs)
    263                     for func, args, kwargs in self.items]
    264 

~/.pyenv/versions/3.8.3/lib/python3.8/site-packages/joblib/parallel.py in <listcomp>(.0)
    260         # change the default number of processes to -1
    261         with parallel_backend(self._backend, n_jobs=self._n_jobs):
--> 262             return [func(*args, **kwargs)
    263                     for func, args, kwargs in self.items]
    264 

~/.pyenv/versions/3.8.3/lib/python3.8/site-packages/sklearn/utils/fixes.py in __call__(self, *args, **kwargs)
    115     def __call__(self, *args, **kwargs):
    116         with config_context(**self.config):
--> 117             return self.function(*args, **kwargs)
    118 
    119 

~/.pyenv/versions/3.8.3/lib/python3.8/site-packages/sklearn/linear_model/_logistic.py in _logistic_regression_path(X, y, pos_class, Cs, fit_intercept, max_iter, tol, verbose, solver, coef, class_weight, dual, penalty, intercept_scaling, multi_class, random_state, check_input, max_squared_sum, sample_weight, l1_ratio, n_threads)
    452             l2_reg_strength = 1.0 / C
    453             args = (X, target, sample_weight, l2_reg_strength, n_threads)
--> 454             w0, n_iter_i = _newton_cg(
    455                 hess, func, grad, w0, args=args, maxiter=max_iter, tol=tol
    456             )

~/.pyenv/versions/3.8.3/lib/python3.8/site-packages/sklearn/utils/optimize.py in _newton_cg(grad_hess, func, grad, x0, args, tol, maxiter, maxinner, line_search, warn)
    191         # Inner loop: solve the Newton update by conjugate gradient, to
    192         # avoid inverting the Hessian
--> 193         xsupi = _cg(fhess_p, fgrad, maxiter=maxinner, tol=termcond)
    194 
    195         alphak = 1.0

~/.pyenv/versions/3.8.3/lib/python3.8/site-packages/sklearn/utils/optimize.py in _cg(fhess_p, fgrad, maxiter, tol)
     86             break
     87 
---> 88         Ap = fhess_p(psupi)
     89         # check curvature
     90         curv = np.dot(psupi, Ap)

~/.pyenv/versions/3.8.3/lib/python3.8/site-packages/sklearn/linear_model/_linear_loss.py in hessp(s)
    344                 if self.fit_intercept:
    345                     ret[:n_features] += s[-1] * hX_sum
--> 346                     ret[-1] = hX_sum @ s[:n_features] + hessian_sum * s[-1]
    347                 return ret
    348 

ValueError: matmul: Input operand 0 does not have enough dimensions (has 0, gufunc core with signature (n?,k),(k,m?)->(n?,m?) requires 1)
```

### Versions

```shell
System:
    python: 3.8.3 (default, Jun 26 2020, 00:34:38)  [Clang 11.0.3 (clang-1103.0.32.62)]
executable: /Users/raphael/.pyenv/versions/3.8.3/bin/python3.8
   machine: macOS-10.16-x86_64-i386-64bit

Python dependencies:
      sklearn: 1.1.1
          pip: 22.1.2
   setuptools: 41.2.0
        numpy: 1.21.6
        scipy: 1.8.0
       Cython: 0.29.23
       pandas: 1.4.1
   matplotlib: 3.5.1
       joblib: 1.0.1
threadpoolctl: 2.1.0

Built with OpenMP: True
```

COMMENT:
Thanks for the detailed report, it made the bug easy to fix.

Fix in #23608

Labels: Bug