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

[mdhaber]

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

[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]

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.

[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.