scipy: BUG: #18481 made BFGS minimization less accurate

Describe your issue.

#18481 appears to have made BFGS minimization less accurate. We noticed this in the JAX upstream nightly tests, which began failing after #18481 was merged: https://github.com/google/jax/actions/workflows/upstream-nightly.yaml

I suspect the dtype changes in approx_derivative are to blame.

Reproducing Code Example

import numpy as np
import scipy

def func(p):
    x, y = p
    return (x ** 2 + y - 11.) ** 2 + (x + y ** 2 - 7.) ** 2

x0 = np.ones(2, dtype='float32')
result = scipy.optimize.minimize(func, x0, method='BFGS')

print(f"{scipy.__version__=}")
print(result)
np.testing.assert_allclose(result.x, [3, 2])

Error message

Output before #18481:

scipy.__version__='1.11.1'
  message: Optimization terminated successfully.
  success: True
   status: 0
      fun: 2.288511533827582e-15
        x: [ 3.000e+00  2.000e+00]
      nit: 8
      jac: [-1.907e-08  2.738e-08]
 hess_inv: [[ 1.612e-02 -9.549e-03]
            [-9.549e-03  3.517e-02]]
     nfev: 39
     njev: 13

Output after #18481:

scipy.__version__='1.12.0.dev0+1440.64be173'
  message: Desired error not necessarily achieved due to precision loss.
  success: False
   status: 2
      fun: 9.628441344453949e-06
        x: [ 3.000e+00  2.000e+00]
      nit: 6
      jac: [ 6.255e-04  2.977e-04]
 hess_inv: [[ 1.615e-02 -9.191e-03]
            [-9.191e-03  3.561e-02]]
     nfev: 291
     njev: 70

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-1-d421a22eaf1d> in <cell line: 13>()
     11 result = scipy.optimize.minimize(func, x0, method='BFGS')
     12 print(result)
---> 13 np.testing.assert_allclose(result.x, [3, 2])

/usr/lib/python3.10/contextlib.py in inner(*args, **kwds)
     77         def inner(*args, **kwds):
     78             with self._recreate_cm():
---> 79                 return func(*args, **kwds)
     80         return inner
     81 

/usr/local/lib/python3.10/dist-packages/numpy/testing/_private/utils.py in assert_array_compare(comparison, x, y, err_msg, verbose, header, precision, equal_nan, equal_inf, strict)
    795                                 verbose=verbose, header=header,
    796                                 names=('x', 'y'), precision=precision)
--> 797             raise AssertionError(msg)
    798     except ValueError:
    799         import traceback

AssertionError: 
Not equal to tolerance rtol=1e-07, atol=0

Mismatched elements: 2 / 2 (100%)
Max absolute difference: 0.00049761
Max relative difference: 0.00016587
 x: array([2.999502, 1.999956])
 y: array([3, 2])

SciPy/NumPy/Python version and system information

The error is seen on the nightly scipy version:

$ pip install -i https://pypi.anaconda.org/scientific-python-nightly-wheels/simple --pre -U scipy