Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

BUG: least_squares with method='trf' with initial params at bound and solution outside of bounds results in division by zero warning #19102

Closed
JoepVanlier opened this issue Aug 21, 2023 · 3 comments
Labels
defect A clear bug or issue that prevents SciPy from being installed or used as expected scipy.optimize
Milestone

Comments

@JoepVanlier
Copy link

Describe your issue.

Not sure if this is considered a regression, but since 1.11.2, using method='trf' on scipy.optimize.least_squares results in a warning about a division by zero when the initial guess is initially exactly at the bound while the solution is outside of the bounds.

What I expected would be no warning with the resulting optimized value being the value at the bound.

Reproducing Code Example

import warnings
import scipy
import numpy as np

warnings.filterwarnings("error")

x = np.arange(0.0, 10.0, 1.0)
def f(a):
    return - a * x - x

scipy.optimize.least_squares(f, -0.5, bounds=(-0.5, 5.0), gtol=1e-6, method='trf')

Error message

---------------------------------------------------------------------------
RuntimeWarning                            Traceback (most recent call last)
Cell In[17], line 11
      8 def f(a):
      9     return - a * x - x
---> 11 scipy.optimize.least_squares(f, -0.5, bounds=(-0.5, 5.0), gtol=1e-6, method='trf')

File /opt/homebrew/Caskroom/miniconda/base/envs/pylake/lib/python3.11/site-packages/scipy/optimize/_lsq/least_squares.py:935, in least_squares(fun, x0, jac, bounds, method, ftol, xtol, gtol, x_scale, loss, f_scale, diff_step, tr_solver, tr_options, jac_sparsity, max_nfev, verbose, args, kwargs)
    931     result = call_minpack(fun_wrapped, x0, jac_wrapped, ftol, xtol, gtol,
    932                           max_nfev, x_scale, diff_step)
    934 elif method == 'trf':
--> 935     result = trf(fun_wrapped, jac_wrapped, x0, f0, J0, lb, ub, ftol, xtol,
    936                  gtol, max_nfev, x_scale, loss_function, tr_solver,
    937                  tr_options.copy(), verbose)
    939 elif method == 'dogbox':
    940     if tr_solver == 'lsmr' and 'regularize' in tr_options:

File /opt/homebrew/Caskroom/miniconda/base/envs/pylake/lib/python3.11/site-packages/scipy/optimize/_lsq/trf.py:123, in trf(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale, loss_function, tr_solver, tr_options, verbose)
    119     return trf_no_bounds(
    120         fun, jac, x0, f0, J0, ftol, xtol, gtol, max_nfev, x_scale,
    121         loss_function, tr_solver, tr_options, verbose)
    122 else:
--> 123     return trf_bounds(
    124         fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale,
    125         loss_function, tr_solver, tr_options, verbose)

File /opt/homebrew/Caskroom/miniconda/base/envs/pylake/lib/python3.11/site-packages/scipy/optimize/_lsq/trf.py:234, in trf_bounds(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale, loss_function, tr_solver, tr_options, verbose)
    232 v, dv = CL_scaling_vector(x, g, lb, ub)
    233 v[dv != 0] *= scale_inv[dv != 0]
--> 234 Delta = norm(x0 * scale_inv / v**0.5)
    235 if Delta == 0:
    236     Delta = 1.0

RuntimeWarning: divide by zero encountered in divide

SciPy/NumPy/Python version and system information

Python 3.10.12 | packaged by conda-forge | (main, Jun 23 2023, 22:41:52) [Clang 15.0.7 ] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info); scipy.show_config()
1.12.0.dev0+1521.da82ac8 1.25.2 sys.version_info(major=3, minor=10, micro=12, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/include
    lib directory: /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/lib
    name: openblas
    openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=0 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK=0 NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP=1 VORTEX MAX_THREADS=128
    pc file directory: /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/lib/pkgconfig
    version: 0.3.23
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/include
    lib directory: /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/lib
    name: openblas
    openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=0 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK=0 NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP=1 VORTEX MAX_THREADS=128
    pc file directory: /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/lib/pkgconfig
    version: 0.3.23
  pybind11:
    detection method: pkgconfig
    include directory: /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/include
    name: pybind11
    version: 2.11.1
Compilers:
  c:
    args: -ftree-vectorize, -fPIC, -fPIE, -fstack-protector-strong, -O2, -pipe, -isystem,
      /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/include, -D_FORTIFY_SOURCE=2,
      -isystem, /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/include
    commands: arm64-apple-darwin20.0.0-clang
    linker: ld64
    linker args: -Wl,-pie, -Wl,-headerpad_max_install_names, -Wl,-dead_strip_dylibs,
      -Wl,-rpath,/opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/lib, -L/opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/lib,
      -ftree-vectorize, -fPIC, -fPIE, -fstack-protector-strong, -O2, -pipe, -isystem,
      /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/include, -D_FORTIFY_SOURCE=2,
      -isystem, /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/include
    name: clang
    version: 15.0.7
  c++:
    args: -ftree-vectorize, -fPIC, -fPIE, -fstack-protector-strong, -O2, -pipe, -stdlib=libc++,
      -fvisibility-inlines-hidden, -fmessage-length=0, -isystem, /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/include,
      -D_FORTIFY_SOURCE=2, -isystem, /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/include
    commands: arm64-apple-darwin20.0.0-clang++
    linker: ld64
    linker args: -Wl,-pie, -Wl,-headerpad_max_install_names, -Wl,-dead_strip_dylibs,
      -Wl,-rpath,/opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/lib, -L/opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/lib,
      -ftree-vectorize, -fPIC, -fPIE, -fstack-protector-strong, -O2, -pipe, -stdlib=libc++,
      -fvisibility-inlines-hidden, -fmessage-length=0, -isystem, /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/include,
      -D_FORTIFY_SOURCE=2, -isystem, /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/include
    name: clang
    version: 15.0.7
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.0
  fortran:
    args: -march=armv8.3-a, -ftree-vectorize, -fPIC, -fno-stack-protector, -O2, -pipe,
      -isystem, /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/include
    commands: /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/bin/arm64-apple-darwin20.0.0-gfortran
    linker: ld64
    linker args: -Wl,-pie, -Wl,-headerpad_max_install_names, -Wl,-dead_strip_dylibs,
      -Wl,-rpath,/opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/lib, -L/opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/lib,
      -march=armv8.3-a, -ftree-vectorize, -fPIC, -fno-stack-protector, -O2, -pipe,
      -isystem, /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/include
    name: gcc
    version: 12.3.0
  pythran:
    include directory: ../../../../../../opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/lib/python3.10/site-packages/pythran
    version: 0.13.1
Machine Information:
  build:
    cpu: aarch64
    endian: little
    family: aarch64
    system: darwin
  cross-compiled: false
  host:
    cpu: aarch64
    endian: little
    family: aarch64
    system: darwin
Python Information:
  path: /opt/homebrew/Caskroom/miniconda/base/envs/scipy-dev/bin/python3.10
  version: '3.10'
@JoepVanlier JoepVanlier added the defect A clear bug or issue that prevents SciPy from being installed or used as expected label Aug 21, 2023
@JoepVanlier JoepVanlier changed the title BUG: Trust region reflective option results in division by zero warning BUG: Trust region reflective option with initial guess at bound results in division by zero warning Aug 21, 2023
@JoepVanlier JoepVanlier changed the title BUG: Trust region reflective option with initial guess at bound results in division by zero warning BUG: least_squares with method='trf' with initial params at bound and solution outside of bounds results in division by zero warning Aug 21, 2023
@hwalinga
Copy link
Contributor

This is a duplicate of my #19103. In my issue, I pinpoint what change has caused this problem, what change to the source code would satisfy the initial problem as well as the current problem, and I show a test case in which there is even a RuntimeError instead of a RuntimeWarning.

@sweettyler
Copy link

Hit by this too. Hope it gets fixed quickly.

hwalinga added a commit to hwalinga/scipy that referenced this issue Aug 22, 2023
PR scipy#18896 removed make_strictly_feasible because it makes solutions
worse if they are close to the boundary.
However, the downstream algorithm that eventually solves the
optimization problem has some trouble with starting points that
are on the boundaries of a bounded optimization problem.
This either results in division by zero warnings, see scipy#19102,
or in failure of finding the correct solution at all,
which is added as a regression test in this commit.
This PR brings back make_strictly_feasible but with the
smallest possible step from the boundary a floating point allows
by setting rstep=0 on the function make_strictly_feasible,
which uses np.nextafter to find such a point, hereby respecting
the strictly feasible point condition and minimizing boundary effects.
See strictly feasible points in interior-point methods for more context
https://nmayorov.wordpress.com/2015/06/19/trust-region-reflective-algorithm/
hwalinga added a commit to hwalinga/scipy that referenced this issue Aug 23, 2023
PR scipy#18896 removed make_strictly_feasible because it makes solutions
worse if they are close to the boundary.
However, the downstream algorithm that eventually solves the
optimization problem has some trouble with starting points that
are on the boundaries of a bounded optimization problem.
This either results in division by zero warnings, see scipy#19102,
or in failure of finding the correct solution at all,
which is added as a regression test in this commit.
This PR brings back make_strictly_feasible but with the
smallest possible step from the boundary a floating point allows
by setting rstep=0 on the function make_strictly_feasible,
which uses np.nextafter to find such a point, hereby respecting
the strictly feasible point condition and minimizing boundary effects.
See strictly feasible points in interior-point methods for more context
https://nmayorov.wordpress.com/2015/06/19/trust-region-reflective-algorithm/
@JoepVanlier
Copy link
Author

JoepVanlier commented Aug 28, 2023

Confirmed as fixed by #19111

@tylerjereddy tylerjereddy added this to the 1.12.0 milestone Aug 28, 2023
tylerjereddy pushed a commit to tylerjereddy/scipy that referenced this issue Sep 21, 2023
PR scipy#18896 removed make_strictly_feasible because it makes solutions
worse if they are close to the boundary.
However, the downstream algorithm that eventually solves the
optimization problem has some trouble with starting points that
are on the boundaries of a bounded optimization problem.
This either results in division by zero warnings, see scipy#19102,
or in failure of finding the correct solution at all,
which is added as a regression test in this commit.
This PR brings back make_strictly_feasible but with the
smallest possible step from the boundary a floating point allows
by setting rstep=0 on the function make_strictly_feasible,
which uses np.nextafter to find such a point, hereby respecting
the strictly feasible point condition and minimizing boundary effects.
See strictly feasible points in interior-point methods for more context
https://nmayorov.wordpress.com/2015/06/19/trust-region-reflective-algorithm/
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
defect A clear bug or issue that prevents SciPy from being installed or used as expected scipy.optimize
Projects
None yet
Development

No branches or pull requests

5 participants