Repository: scipy/scipy
Pull request: BUG: make `L-BFGS-B` optimizer work with single precision gradient
#### Reference issue
Closes #18730 

#### What does this implement/fix?
Casts the gradient array to double precision as Fortran expects a double precision array. I am not sure if this is the best place for it.


Problem statement:
BUG: `scipy.optimize.minimize` fails when `dtype=float32` Jacobian is passed
### Describe your issue.

`scipy.optimize.minimize` only accept `dtype=float64` Jacobian (the called Fortran routine takes `float64/double`). If fed with `dtype=float32` Jacobian, the failing message doesn't help users debug too much.

`dtype=float32` gradient is quite common for users from deep learning community. Is that possible to do some extra chck in the `scipy.optimize.minimize` function, e.g. `if jac.dtype != np.float64: raise ValueError('Jacobian must be float64')` or simply `jac = jac.astype(np.float64)`? Many thx.

### Reproducing Code Example

```python
import numpy as np
import scipy.optimize
N0 = 4
# return value from hf0/hf1: function value, jacobian
hf0 = lambda x: (np.sum(x**2), (2*x).astype(np.float64))
hf1 = lambda x: (np.sum(x**2), (2*x).astype(np.float32))
x0 = np.random.default_rng().uniform(0, 1, size=N0).astype(np.float64) #float32/float64 doesn't matter here
method = 'L-BFGS-B' #BFGS #same
theta_optim = scipy.optimize.minimize(hf0, x0, jac=True, method=method) #pass
# theta_optim = scipy.optimize.minimize(hf1, x0, jac=True, method=method) #fail
```


### Error message

```shell
File ~/micromamba/envs/cuda118/lib/python3.10/site-packages/scipy/optimize/_lbfgsb_py.py:350, in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, finite_diff_rel_step, **unknown_options)
    346 n_iterations = 0
    348 while 1:
    349     # x, f, g, wa, iwa, task, csave, lsave, isave, dsave = \
--> 350     _lbfgsb.setulb(m, x, low_bnd, upper_bnd, nbd, f, g, factr,
    351                    pgtol, wa, iwa, task, iprint, csave, lsave,
    352                    isave, dsave, maxls)
    353     task_str = task.tobytes()
    354     if task_str.startswith(b'FG'):
    355         # The minimization routine wants f and g at the current x.
    356         # Note that interruptions due to maxfun are postponed
    357         # until the completion of the current minimization iteration.
    358         # Overwrite f and g:

ValueError: failed in converting 7th argument `g' of _lbfgsb.setulb to C/Fortran array
```


### SciPy/NumPy/Python version and system information

<details>

```shell
1.10.1 1.24.3 sys.version_info(major=3, minor=10, micro=11, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /home/zhangc/micromamba/envs/cuda118/include
    lib directory: /home/zhangc/micromamba/envs/cuda118/lib
    name: blas
    openblas configuration: unknown
    pc file directory: /home/zhangc/micromamba/envs/cuda118/lib/pkgconfig
    version: 3.9.0
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /home/zhangc/micromamba/envs/cuda118/include
    lib directory: /home/zhangc/micromamba/envs/cuda118/lib
    name: lapack
    openblas configuration: unknown
    pc file directory: /home/zhangc/micromamba/envs/cuda118/lib/pkgconfig
    version: 3.9.0
Compilers:
  c:
    commands: /home/conda/feedstock_root/build_artifacts/scipy_1683719288579/_build_env/bin/x86_64-conda-linux-gnu-cc
    linker: ld.bfd
    name: gcc
    version: 11.3.0
  c++:
    commands: /home/conda/feedstock_root/build_artifacts/scipy_1683719288579/_build_env/bin/x86_64-conda-linux-gnu-c++
    linker: ld.bfd
    name: gcc
    version: 11.3.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 0.29.34
  fortran:
    commands: /home/conda/feedstock_root/build_artifacts/scipy_1683719288579/_build_env/bin/x86_64-conda-linux-gnu-gfortran
    linker: ld.bfd
    name: gcc
    version: 11.3.0
  pythran:
    include directory: ../../_h_env_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_place/lib/python3.10/site-packages/pythran
    version: 0.12.2
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
Python Information:
  path: /home/zhangc/micromamba/envs/cuda118/bin/python
  version: '3.10'
```

</details>

Hints from issue discussion:
Thank you for reporting @husisy. I can reproduce locally. Until we fix that (or not) on the fortran side, would you be interested in adding some input validations?
I think we should do any casting lower down, rather than creating a high level wrapper in `minimize`. e.g. BFGS is pure python, so we probably don't want to cast to 64 because we have a chance of using the 32 bit grad directly. In comparison, with LBFGSB we'll need to cast a 32 grad to 64 just before it's sent back into the fortran code.
To be clear, I was not proposing to do any casting but just to raise properly in Python with a nice message saying that this function does not (yet) supports 32 bits.

Commit message: TST: strengthen test by making it 2 dimensional

Patch to review:
diff --git a/scipy/optimize/_lbfgsb_py.py b/scipy/optimize/_lbfgsb_py.py
index a4d06954178f..7def7605985a 100644
--- a/scipy/optimize/_lbfgsb_py.py
+++ b/scipy/optimize/_lbfgsb_py.py
@@ -196,8 +196,8 @@ def fmin_l_bfgs_b(func, x0, fprime=None, args=(),
             'callback': callback,
             'maxls': maxls}
 
-    res = _minimize_lbfgsb(fun, x0, args=args, jac=jac, bounds=bounds,
-                           **opts)
+    res = _minimize_lbfgsb(fun, x0, args=args, jac=jac,
+                           bounds=bounds, **opts)
     d = {'grad': res['jac'],
          'task': res['message'],
          'funcalls': res['nfev'],
@@ -348,6 +348,8 @@ def _minimize_lbfgsb(fun, x0, args=(), jac=None, bounds=None,
     n_iterations = 0
 
     while 1:
+        # convert gradient to double precision for Fortran
+        g = g.astype(np.float64)
         # x, f, g, wa, iwa, task, csave, lsave, isave, dsave = \
         _lbfgsb.setulb(m, x, low_bnd, upper_bnd, nbd, f, g, factr,
                        pgtol, wa, iwa, task, iprint, csave, lsave,
diff --git a/scipy/optimize/tests/test_lbfgsb_setulb.py b/scipy/optimize/tests/test_lbfgsb_setulb.py
index 3d9fc254aade..d3ebb19e7e99 100644
--- a/scipy/optimize/tests/test_lbfgsb_setulb.py
+++ b/scipy/optimize/tests/test_lbfgsb_setulb.py
@@ -1,5 +1,5 @@
 import numpy as np
-from scipy.optimize import _lbfgsb
+from scipy.optimize import _lbfgsb, minimize
 
 
 def objfun(x):
@@ -114,3 +114,11 @@ def test_setulb_floatround():
 
         assert (x <= upper_bnd).all() and (x >= low_bnd).all(), (
             "_lbfgsb.setulb() stepped to a point outside of the bounds")
+
+
+def test_gh_issue18730():
+    def fun(x):
+        return np.sum(x**2), (2*x).astype(np.float32)
+
+    res = minimize(fun, x0=np.array([1., 1.]), jac=True, method="l-bfgs-b")
+    np.testing.assert_allclose(res.fun, 0., atol=1e-15)


Code review comments:
diff_hunk: @@ -196,8 +196,8 @@ def fmin_l_bfgs_b(func, x0, fprime=None, args=(),
             'callback': callback,
             'maxls': maxls}
 
-    res = _minimize_lbfgsb(fun, x0, args=args, jac=jac, bounds=bounds,
-                           **opts)
+    res = _minimize_lbfgsb(fun, x0, args=args, jac=jac,
+                           bounds=bounds, **opts)
path: scipy/optimize/_lbfgsb_py.py
text: @user1:
minor: best to avoid unnecessary style tweaks like this. If this PR would be ready to merge I wouldn't comment, but if you're pushing a follow-up commit anyway, would you mind reverting this change?

diff_hunk: @@ -114,3 +114,11 @@ def test_setulb_floatround():
 
         assert (x <= upper_bnd).all() and (x >= low_bnd).all(), (
             "_lbfgsb.setulb() stepped to a point outside of the bounds")
+
+
+def test_gh_issue18730():
+    def fun(x):
+        return np.sum(x**2), (2*x).astype(np.float32)
path: scipy/optimize/tests/test_lbfgsb_setulb.py
text: @user1:
I was considering whether this function is actually valid to pass, because it's a little odd at least to have different dtypes for the two return values. However, if the first return value is also `float32`, the failure mode remains the same - and then it's definitely a valid/normal function. So consider tweaking this test like so:
```python
def test_gh_issue18730():
    def fun(x):
        # A function that does things only in float32
        x = x.astype(np.float32)
        return np.sum(x**2), (2*x)
```