MarkWieczorek commented 1 month ago

Describe your issue.

In scipy 1.12, a python version of the NNLS (non-negative least squares) routine replaced the tried and tested fortran version. Unfortunately, the behavior of the new routine (which claims to be faster) is sometimes inconsistent with the old routine, and the new routine gives no indication of when the results are unreliable.

Description of problem

I am trying to solve the problem

min $$||Ax - b||$$

with the constraint that all elements of x are non-negative. With A having a shape (n, m), a maximum of n elements of x will be non-zero.

For my problem, the matrix A has shape (1074, 1257), the minimum and maximum values of the elements of A are -2.4193018654733958e-21, and 3.943646586725562e-21, and the minimum and maximum values of b are -1.4210278995247266e-08 and 1.217756814911735e-08.

calling (using scipy 1.14.1)

m, rnorm = nnls(A, b)

yields a solution vector m that is identically zero.

Since A and b were expressed in SI units of Teslas, I decided to convert to nT instead

m, rnorm = nnls(A*1.e9, b*1.e9)

This returns a solution vector that has 102 non-zero elements.

Noting that the elements of A were small, I then tried to normalize A and b by the maximum value of A

m, rnorm = nnls(A/A.max(), b/A.max())

This returns a solution vector that has 298 non-zero elements.

Next, I compared with the old Fortran based routine in scipy 1.11. For this, all three tests above returned the exact same answer with 298 non-zero elements in the solution vector.

Conclusion

As written, the 1.12+ version of nnls is not consistent with the 1.11 Fortran version.
The python version can be made to work by predividing A and b by A.max(). There might be a better way, as I am not that familiar with how things work under the hood.
If the python version can not be fixed, I would suggest giving the user the option to use either the fortran or python version. This would at least give confidence that the python-version results are correct, which they are not at the present time for my problem.

Related to https://github.com/scipy/scipy/issues/21484 Ping @ilayn

Reproducing Code Example

import numpy as np
from scipy.optimize import nnls

data = np.load('nnls.npz')
A = data['A']
b = data['b']

print(A.min(), A.max())
print(b.min(), b.max())

m, rnorm = nnls(A, b)
print(sum(m>0), rnorm)

m, rnorm = nnls(A*1.e9, b*1.e9)
print(sum(m>0), rnorm/1.e9)

m, rnorm = nnls(A/A.max(), b/A.max())
print(sum(m>0), rnorm*A.max())

Error message

There is no error message

SciPy/NumPy/Python version and system information

In [2]: import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.versi
   ...: on_info); scipy.show_config()
1.14.1 2.1.0 sys.version_info(major=3, minor=12, micro=5, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /home/mrak/miniforge3/envs/py312/include
    lib directory: /home/mrak/miniforge3/envs/py312/lib
    name: blas
    openblas configuration: unknown
    pc file directory: /home/mrak/miniforge3/envs/py312/lib/pkgconfig
    version: 3.9.0
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /home/mrak/miniforge3/envs/py312/include
    lib directory: /home/mrak/miniforge3/envs/py312/lib
    name: lapack
    openblas configuration: unknown
    pc file directory: /home/mrak/miniforge3/envs/py312/lib/pkgconfig
    version: 3.9.0
  pybind11:
    detection method: pkgconfig
    include directory: /home/mrak/miniforge3/envs/py312/include
    name: pybind11
    version: 2.13.4
Compilers:
  c:
    args: -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, /home/mrak/miniforge3/envs/py312/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1724327061663/work=/usr/local/src/conda/scipy-split-1.14.1,
      -fdebug-prefix-map=/home/mrak/miniforge3/envs/py312=/usr/local/src/conda-prefix,
      -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, /home/mrak/miniforge3/envs/py312/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1724327061663/work=/usr/local/src/conda/scipy-split-1.14.1,
      -fdebug-prefix-map=/home/mrak/miniforge3/envs/py312=/usr/local/src/conda-prefix,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, /home/mrak/miniforge3/envs/py312/include,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, /home/mrak/miniforge3/envs/py312/include
    commands: /home/conda/feedstock_root/build_artifacts/scipy-split_1724327061663/_build_env/bin/x86_64-conda-linux-gnu-cc
    linker: ld.bfd
    linker args: -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now,
      -Wl,--disable-new-dtags, -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,/home/mrak/miniforge3/envs/py312/lib,
      -Wl,-rpath-link,/home/mrak/miniforge3/envs/py312/lib, -L/home/mrak/miniforge3/envs/py312/lib,
      -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now, -Wl,--disable-new-dtags,
      -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,/home/mrak/miniforge3/envs/py312/lib,
      -Wl,-rpath-link,/home/mrak/miniforge3/envs/py312/lib, -L/home/mrak/miniforge3/envs/py312/lib,
      -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, /home/mrak/miniforge3/envs/py312/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1724327061663/work=/usr/local/src/conda/scipy-split-1.14.1,
      -fdebug-prefix-map=/home/mrak/miniforge3/envs/py312=/usr/local/src/conda-prefix,
      -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, /home/mrak/miniforge3/envs/py312/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1724327061663/work=/usr/local/src/conda/scipy-split-1.14.1,
      -fdebug-prefix-map=/home/mrak/miniforge3/envs/py312=/usr/local/src/conda-prefix,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, /home/mrak/miniforge3/envs/py312/include,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, /home/mrak/miniforge3/envs/py312/include
    name: gcc
    version: 13.3.0
  c++:
    args: -fvisibility-inlines-hidden, -fmessage-length=0, -march=nocona, -mtune=haswell,
      -ftree-vectorize, -fPIC, -fstack-protector-strong, -fno-plt, -O2, -ffunction-sections,
      -pipe, -isystem, /home/mrak/miniforge3/envs/py312/include, -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1724327061663/work=/usr/local/src/conda/scipy-split-1.14.1,
      -fdebug-prefix-map=/home/mrak/miniforge3/envs/py312=/usr/local/src/conda-prefix,
      -fvisibility-inlines-hidden, -fmessage-length=0, -march=nocona, -mtune=haswell,
      -ftree-vectorize, -fPIC, -fstack-protector-strong, -fno-plt, -O2, -ffunction-sections,
      -pipe, -isystem, /home/mrak/miniforge3/envs/py312/include, -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1724327061663/work=/usr/local/src/conda/scipy-split-1.14.1,
      -fdebug-prefix-map=/home/mrak/miniforge3/envs/py312=/usr/local/src/conda-prefix,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, /home/mrak/miniforge3/envs/py312/include,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, /home/mrak/miniforge3/envs/py312/include
    commands: /home/conda/feedstock_root/build_artifacts/scipy-split_1724327061663/_build_env/bin/x86_64-conda-linux-gnu-c++
    linker: ld.bfd
    linker args: -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now,
      -Wl,--disable-new-dtags, -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,/home/mrak/miniforge3/envs/py312/lib,
      -Wl,-rpath-link,/home/mrak/miniforge3/envs/py312/lib, -L/home/mrak/miniforge3/envs/py312/lib,
      -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now, -Wl,--disable-new-dtags,
      -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,/home/mrak/miniforge3/envs/py312/lib,
      -Wl,-rpath-link,/home/mrak/miniforge3/envs/py312/lib, -L/home/mrak/miniforge3/envs/py312/lib,
      -fvisibility-inlines-hidden, -fmessage-length=0, -march=nocona, -mtune=haswell,
      -ftree-vectorize, -fPIC, -fstack-protector-strong, -fno-plt, -O2, -ffunction-sections,
      -pipe, -isystem, /home/mrak/miniforge3/envs/py312/include, -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1724327061663/work=/usr/local/src/conda/scipy-split-1.14.1,
      -fdebug-prefix-map=/home/mrak/miniforge3/envs/py312=/usr/local/src/conda-prefix,
      -fvisibility-inlines-hidden, -fmessage-length=0, -march=nocona, -mtune=haswell,
      -ftree-vectorize, -fPIC, -fstack-protector-strong, -fno-plt, -O2, -ffunction-sections,
      -pipe, -isystem, /home/mrak/miniforge3/envs/py312/include, -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1724327061663/work=/usr/local/src/conda/scipy-split-1.14.1,
      -fdebug-prefix-map=/home/mrak/miniforge3/envs/py312=/usr/local/src/conda-prefix,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, /home/mrak/miniforge3/envs/py312/include,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, /home/mrak/miniforge3/envs/py312/include
    name: gcc
    version: 13.3.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.11
  fortran:
    args: -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, /home/mrak/miniforge3/envs/py312/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1724327061663/work=/usr/local/src/conda/scipy-split-1.14.1,
      -fdebug-prefix-map=/home/mrak/miniforge3/envs/py312=/usr/local/src/conda-prefix
    commands: /home/conda/feedstock_root/build_artifacts/scipy-split_1724327061663/_build_env/bin/x86_64-conda-linux-gnu-gfortran
    linker: ld.bfd
    linker args: -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now,
      -Wl,--disable-new-dtags, -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,/home/mrak/miniforge3/envs/py312/lib,
      -Wl,-rpath-link,/home/mrak/miniforge3/envs/py312/lib, -L/home/mrak/miniforge3/envs/py312/lib,
      -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now, -Wl,--disable-new-dtags,
      -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,/home/mrak/miniforge3/envs/py312/lib,
      -Wl,-rpath-link,/home/mrak/miniforge3/envs/py312/lib, -L/home/mrak/miniforge3/envs/py312/lib,
      -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, /home/mrak/miniforge3/envs/py312/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1724327061663/work=/usr/local/src/conda/scipy-split-1.14.1,
      -fdebug-prefix-map=/home/mrak/miniforge3/envs/py312=/usr/local/src/conda-prefix
    name: gcc
    version: 13.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/lib/python3.12/site-packages/pythran
    version: 0.16.1
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/mrak/miniforge3/envs/py312/bin/python
  version: '3.12'

ilayn commented 1 month ago

Do you have any means to build the current main by any chance? Otherwise please paste the array somewhere online like pastebin or a gist and I'll see if the result of the latest version (#21021) is in accordance with your expectations.

MarkWieczorek commented 1 month ago

Build scipy? No. Sorry. Not that much free time...

I'll upload the file somewhere today. However, given that this appears to be a numerical resolution issue, I wonder if you just took any matrix A and vector b and divided them by 1.e20 if that would be sufficient to reproduce this.

ilayn commented 1 month ago

Build scipy? No. Sorry. Not that much free time...

Same applies to us. Hence my question.

I wonder if you just took any matrix A and vector b and divided them by 1.e20 if that would be sufficient to reproduce this.

We are trying to find if the new code solves, not the old code fails. Otherwise we have a number of examples for the old code to fail and we fixed them. We should see if it also fixed your problem.

MarkWieczorek commented 1 month ago

nnls.npz.zip

Here are the data, and I updated the original post with simple script to reproduce this.

ilayn commented 1 month ago

On the current main I get

-2.4193018654733958e-21 3.943646586725562e-21
-1.4210278995247266e-08 1.217756814911735e-08
298 4.3438564761143895e-08
298 4.343856476114391e-08
298 4.3438564761143895e-08

So it looks like it is working properly, but as a recommendation, it is always better to scale your data because it improves the conditioning.

MarkWieczorek commented 1 month ago

To confirm:

This is a bug with version 1.12, 1.13 and 1.14. The bug will be fixed in version 1.15.

ilayn commented 1 month ago

That's correct.

scipy / scipy

BUG: optimize.nnls: precision problems #21630

Describe your issue.

Description of problem

Conclusion

Reproducing Code Example

Error message

SciPy/NumPy/Python version and system information