BUG: scipy.optimize.linprog gives optimized solution which does not respect linear equation

[Markus-Weiand]

Describe your issue.

in the code example below the solution fond by linprog does not respect the linear equation even though it claims otherwise.
If I run the code snippet below, I get the following output

norm residual x on input equation is 0.021185053416350418
norm residual x given by linprog is 0.0

Reproducing Code Example

import numpy as np # version 1.26.4
import scipy       # version 1.12.0

if __name__ == '__main__':
    boundaries=[(10000.0, 9010000.0), (0.0, None), (10000.0, None), (0.0, 84.62623413258109), (10000.0, None), (10000.0, None), (10000.0, None), (10000.0, None), (10000.0, None), (10000.0, None), (10000.0, None), (10000.0, None), (10000.0, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None), (None, None)]
    eq_entries=[-1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 0.001, -0.001, 3.7337777768059636e-10, 3.7337777768059636e-10, 1.0, -1.0, 0.001, -0.001, 3.7337777768059636e-10, 3.7337777768059636e-10, 1.0, -1.0, 0.001, -0.001, 3.7337777768059636e-10, 3.7337777768059636e-10, 1.0, -1.0, 0.001, -0.001, 3.7337777768059636e-10, 3.7337777768059636e-10, 1.0, -1.0, 0.001, -0.001, 3.7337777768059636e-10, 3.7337777768059636e-10, 1.0, -1.0, 0.001, -0.001, 3.7337777768059636e-10, 3.7337777768059636e-10, 1.0, -1.0, 0.001, -0.001, 3.7337777768059636e-10, 3.7337777768059636e-10, 1.0, -1.0, 0.001, -0.001, 3.7337777768059636e-10, 3.7337777768059636e-10, 1.0, -1.0, 0.001, -0.001, 3.7337777768059636e-10, 3.7337777768059636e-10, 1.0, -1.0, 0.001, -0.001, 3.7337777768059636e-10, 3.7337777768059636e-10, 1.0, -1.0]
    eq_indizes=[0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 11, 11, 12, 12, 12, 12, 13, 13, 14, 14, 14, 14, 15, 15, 16, 16, 16, 16, 17, 17, 18, 18, 18, 18, 19, 19, 20, 20, 20, 20, 21, 21, 22, 22, 22, 22, 23, 23, 24, 24, 24, 24, 25, 25, 26, 26, 26, 26, 27, 27, 28, 28, 28, 28, 29, 29, 30, 30, 30, 30, 31, 31]
    eq_vars=[15, 14, 17, 16, 19, 18, 21, 20, 23, 22, 25, 24, 27, 26, 29, 28, 31, 30, 13, 1, 0, 32, 3, 14, 13, 4, 0, 4, 0, 32, 31, 2, 12, 2, 12, 16, 15, 5, 4, 5, 4, 18, 17, 6, 5, 6, 5, 20, 19, 7, 6, 7, 6, 22, 21, 8, 7, 8, 7, 24, 23, 9, 8, 9, 8, 26, 25, 10, 9, 10, 9, 28, 27, 11, 10, 11, 10, 30, 29, 12, 11, 12, 11]
    eq_values=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 9000000.0, 0.0, 0.006587392118285457, -5032.197406716549, 0.0041860502789104696, -7918.93439542944, 0.0063205763583549035, -5244.625751707402, 0.006053760598424349, -5475.7793929428, 0.005786944838493795, -5728.248403917573, 0.0055201290785632405, -6005.123623538355, 0.005253313318632687, -6310.123825488683, 0.004986497558702133, -6647.763714796453, 0.004719681798771578, -7023.578908071522, 0.004452866038841024, -7444.431798646482]
    coefficients=[0.0, 0.0, 0.0, -0.011816666666666668, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    np_eq_entries=np.asarray(eq_entries,dtype=np.float64)
    np_eq_indizes=np.asarray(eq_indizes,dtype=np.int32)
    np_eq_vars=np.asarray(eq_vars,dtype=np.int32)

    a_eq=scipy.sparse.csr_array((np_eq_entries,(np_eq_indizes, np_eq_vars)),shape=(32,33))
    b_eq=np.asarray(eq_values,dtype=np.float64)
    c=np.asarray(coefficients,dtype=np.float64)
    assert len(b_eq)==32
    assert len(c)==33
    assert len(boundaries)==33
    result = scipy.optimize.linprog(c, A_ub=None, b_ub=None, A_eq=a_eq, b_eq=b_eq, bounds=boundaries)
    assert result.status==0
    x=result.x
    b_x=a_eq @ x
    r_x=b_x-b_eq
    n_r_x=np.linalg.norm(r_x)
    r=result.eqlin.residual
    n_r=np.linalg.norm(r)
    print("norm residual x on input equation is %s"%n_r_x)
    print("norm residual x given by linprog is %s"%n_r)

Error message

None

SciPy/NumPy/Python version and system information

1.12.0 1.26.4 sys.version_info(major=3, minor=9, micro=13, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /c/opt/64/include
    lib directory: /c/opt/64/lib
    name: openblas
    openblas configuration: USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= SKYLAKEX MAX_THREADS=2
    pc file directory: c:/opt/64/lib/pkgconfig
    version: 0.3.21.dev
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /c/opt/64/include
    lib directory: /c/opt/64/lib
    name: openblas
    openblas configuration: USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= SKYLAKEX MAX_THREADS=2
    pc file directory: c:/opt/64/lib/pkgconfig
    version: 0.3.21.dev
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.11.1
Compilers:
  c:
    commands: cc
    linker: ld.bfd
    name: gcc
    version: 10.3.0
  c++:
    commands: c++
    linker: ld.bfd
    name: gcc
    version: 10.3.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.8
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 10.3.0
  pythran:
    include directory: C:\Users\runneradmin\AppData\Local\Temp\pip-build-env-rikezz0v\overlay\Lib\site-packages/pythran
    version: 0.15.0
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: windows
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: windows
Python Information:
  path: C:\Users\runneradmin\AppData\Local\Temp\cibw-run-ee1litsm\cp39-win_amd64\build\venv\Scripts\python.exe
  version: '3.9'

[mdhaber]

This seems to be a scaling issue. Scaling the left and right hand side like:

    i = np.abs(b_eq) > 0
    a_eq[i, :] = a_eq[i, :] / b_eq[i, np.newaxis]
    b_eq[i] = b_eq[i]/b_eq[i]

Results in

norm residual x on input equation is 3.238717519617626e-13

I'm not sure what is going on with the norm residual provided by linprog. I'd rather not dig until we see if the problem is fixed with an update.

@HaoZeke are you able to check this out with a more recent version of HiGHS?

[Markus-Weiand]

as the correct residuum is a bit large, the problem seems most likely to be that linprog didn't solve the problem or solved the wrong problem. I now use google ortools, which didn't have the problem.

[fancidev]

Does it help by tweaking the various tolerance parameters to the function? It seems the default tolerances (1e-7 to 1e-8) are rather “tolerant”.

In general, I find it a hard exercise as an average user to pick a tolerance for many numerical routines of this kind. it would be ideal if these routines specify a default tolerance that is “as accurate as possible”, or provide clear guidance on how to pick the tolerances. At least, elaborate the exact definition of the tolerance in the doc (just “feasibility tolerance” is not clear enough).

[Markus-Weiand]

@fancidev: You are trying to solve my optimization problem. The problem here is a bug with linprog. It may be ok if linprog doesn't find a solution, it may be ok if it finds a bad (almost) solution with a residuum>0, but claiming to have found a solution with residuum 0 when in fact the solution given has a quite large residuum is simple a bug.
I'm now using google ortools and can work with these, but I think people should also be able to work with scipy ....