BUG/ENH?: scipy.linalg.norm: overflow/underflow when `axis` is provided

[mdhaber]

Describe your issue.

I'm working on adding an axis argument to scipy.stats.pearsonr. It uses scipy.linalg.norm instead of (e.g.) np.linalg.norm because scipy.linalg.norm tends to avoid premature overflow.

from scipy import linalg
x = [-5e210, 5e210, 3e200, -3e200]
linalg.norm(x)
# 7.071067811865475e+210

Unfortunately, that advantage is lost when axis is not None.

linalg.norm(x, axis=0)  # warning and inf

Similarly, premature underflow can occur when the argument has small magnitude elements and axis is not None.

I can work around it if need be (e.g. manually scale), but can linalg.norm avoid premature under/overflow regardless of axis?

Reproducing Code Example

from scipy import linalg
x = [-5e210, 5e210, 3e200, -3e200]
linalg.norm(x, axis=0)

Error message

RuntimeWarning: overflow encountered in multiply s = (x.conj() * x).real

SciPy/NumPy/Python version and system information

import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info); scipy.show_config()
1.13.0.dev0+1376.8360bac 1.26.0 sys.version_info(major=3, minor=11, micro=6, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /Users/matthaberland/miniforge3/envs/scipy-dev/include
    lib directory: /Users/matthaberland/miniforge3/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: /Users/matthaberland/miniforge3/envs/scipy-dev/lib/pkgconfig
    version: 0.3.24
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /Users/matthaberland/miniforge3/envs/scipy-dev/include
    lib directory: /Users/matthaberland/miniforge3/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: /Users/matthaberland/miniforge3/envs/scipy-dev/lib/pkgconfig
    version: 0.3.24
  pybind11:
    detection method: pkgconfig
    include directory: /Users/matthaberland/miniforge3/envs/scipy-dev/include
    name: pybind11
    version: 2.11.1
Compilers:
  c:
    args: -ftree-vectorize, -fPIC, -fPIE, -fstack-protector-strong, -O2, -pipe, -isystem,
      /Users/matthaberland/miniforge3/envs/scipy-dev/include, -D_FORTIFY_SOURCE=2,
      -isystem, /Users/matthaberland/miniforge3/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,/Users/matthaberland/miniforge3/envs/scipy-dev/lib, -L/Users/matthaberland/miniforge3/envs/scipy-dev/lib,
      -ftree-vectorize, -fPIC, -fPIE, -fstack-protector-strong, -O2, -pipe, -isystem,
      /Users/matthaberland/miniforge3/envs/scipy-dev/include, -D_FORTIFY_SOURCE=2,
      -isystem, /Users/matthaberland/miniforge3/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, /Users/matthaberland/miniforge3/envs/scipy-dev/include,
      -D_FORTIFY_SOURCE=2, -isystem, /Users/matthaberland/miniforge3/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,/Users/matthaberland/miniforge3/envs/scipy-dev/lib, -L/Users/matthaberland/miniforge3/envs/scipy-dev/lib,
      -ftree-vectorize, -fPIC, -fPIE, -fstack-protector-strong, -O2, -pipe, -stdlib=libc++,
      -fvisibility-inlines-hidden, -fmessage-length=0, -isystem, /Users/matthaberland/miniforge3/envs/scipy-dev/include,
      -D_FORTIFY_SOURCE=2, -isystem, /Users/matthaberland/miniforge3/envs/scipy-dev/include
    name: clang
    version: 15.0.7
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.4
  fortran:
    args: -march=armv8.3-a, -ftree-vectorize, -fPIC, -fno-stack-protector, -O2, -pipe,
      -isystem, /Users/matthaberland/miniforge3/envs/scipy-dev/include
    commands: /Users/matthaberland/miniforge3/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,/Users/matthaberland/miniforge3/envs/scipy-dev/lib, -L/Users/matthaberland/miniforge3/envs/scipy-dev/lib,
      -march=armv8.3-a, -ftree-vectorize, -fPIC, -fno-stack-protector, -O2, -pipe,
      -isystem, /Users/matthaberland/miniforge3/envs/scipy-dev/include
    name: gcc
    version: 12.3.0
  pythran:
    include directory: ../../../miniforge3/envs/scipy-dev/lib/python3.11/site-packages/pythran
    version: 0.14.0
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: /Users/matthaberland/miniforge3/envs/scipy-dev/bin/python
  version: '3.11'

[lucascolley]

np.linalg.norm is used whenever axis is provided:

scipy/scipy/linalg/_misc.py

Lines 177 to 178 in fdf3b90

	# fall back to numpy in every other case
	return np.linalg.norm(a, ord=ord, axis=axis, keepdims=keepdims)

[ilayn]

On top of what @lucascolley mentioned, NumPy does not use the safe norm from LAPACK but does its own unscaled dot product hence bypasses the safety guards. If you don't have too many cols/rows in the axis a for loop is an easy way out without too much of python overhead.

[ilayn]

Should we fix 2-norm working with axis? I can cook up something in C or Cython and nrm2 is not complicated.

Based on what we discussed over kron deprecation, not sure what it entails in terms of array api though since this is a central function and unlike kron we should support it and make it as general as possible, in my opinion. However if I cook something up in C or Cython not sure if it is even possible to support it.

[dschmitz89]

From what I understand, numpy should be able to use nrm2 as well as it is BLAS, not LAPACK. In the long run, would it not be best if this was fixed in numpy itself on top of our own efforts in scipy? norm is definitely much less involved than many other linalg functions from numpy such as the matrix decompositions but a fundamentally important function in my opinion.

[lucascolley]

that sounds sensible - see data-apis/array-api#213 for the array API perspective. We'd probably rewrite scipy.linalg.norm in terms of xp.linalg.vector_norm and xp.linalg.matrix_norm for the non-np codepath.

If we can upstream all efficiency gains from the SciPy implementation to those (new) NumPy functions, we'll be able to simplify it to just one codepath (well, there'll be a codepath for if the namespace doesn't implement the linalg extension, but that's orthogonal to the discussion here).

[dschmitz89]

Oh my, I typed too fast probably: openblas does not support the safe scaling yet apparently: OpenMathLib/OpenBLAS#4313. Until it is available there, probably difficult to use for both numpy and scipy but I would like to be convinced of the opposite by the BLAS gurus.

Xref numpy issue numpy/numpy#19097

[ilayn]

It doesn't have to use the BLAS code as all based on Anderson's work. We can write our own native version of nrm2 in C or Cython. Not sure if NumPy implements it since it is not using nrm2 even now.

[ilayn]

Overall judging by the discussions oat the array api, it seems to me the linalg parts are not well thought out.