Skip to content

BUG: Fix linalg.norm to handle empty matrices correctly. #28343

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

Merged
merged 6 commits into from
Mar 12, 2025

Conversation

carlosgmartin
Copy link
Contributor

Fixes linalg.norm to correctly handle empty matrices, which currently raises a ValueError sometimes. For example:

$ py -c "import numpy as np; print(np.linalg.norm(np.zeros([2, 0]), ord=1))"
ValueError: zero-size array to reduction operation maximum which has no identity

$ py -c "import numpy as np; print(np.linalg.norm(np.zeros([0, 2]), ord=float('inf')))"
ValueError: zero-size array to reduction operation maximum which has no identity

This resolves the following issues:

@eendebakpt
Copy link
Contributor

I tested this locally and works fine. Also the types are stable:

np.linalg.norm(numpy.array([], dtype=np.float16), ord=np.inf)

has the same output type (np.float16) as the input.

@carlosgmartin Could you add some unit tests for the cases that have been fixed?

@carlosgmartin
Copy link
Contributor Author

@eendebakpt Done.

Copy link
Contributor

@eendebakpt eendebakpt left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks good to me.

@carlosgmartin
Copy link
Contributor Author

@eendebakpt Done.

@eendebakpt
Copy link
Contributor

eendebakpt commented Feb 18, 2025

@carlosgmartin Thanks for the updates. Failing tests seem unrelated, so I approved. The PR needs approval of a numpy dev to be merged though.

@seberg
Copy link
Member

seberg commented Feb 18, 2025

@charris do you have a very clear picture of what the right thing is?

@carlosgmartin you seem to have just extended this with from the issue for initial=inf without any explanation, please provide one if you have.
Without thinking about it further, the one thing that makes sense to me is to that an empty matrix always has a norm of zero, but I am not 100% sure zero is safe for all norms.

Also ping @melissawm, maybe you have this type of knowledge more available than me. If there are some norms where an empty matrix default value isn't pretty clear, I would be fine with adding a default= kwarg (but I it may be that the need for this is so small that I am not sure we will agree on it being worthwhile).
(NaN might also be a reasonable default for some in that case.)

EDIT: I think this may be an interesting example/reason for why maximum/minimum indeed shouldn't define -inf/inf as a default.

@melissawm
Copy link
Member

So - all of these norms should indeed be 0 on zero matrices. However, some of these norms are the result of computations e.g. "the sum of singular values". In that case it doesn't seem like assert_equal is the right choice for these tests. I haven't looked that closely at the implementation though!

On the other hand, the initial issue was about empty matrices, not zero-valued matrices. So I'm unsure if this solves that issue.

@seberg
Copy link
Member

seberg commented Feb 18, 2025

@melissawm all of the tests work with empty matrices (so zero or empty would be the same)?
As you said, the question is really about what the right thing is for empty matrices, and that is what the PR modifies.

EDIT: Ahhh, zero matrix => zero norm, I suppose also means that empty matrix == empty zero matrix => zero norm. Unless there is some argument to the contrary?

@melissawm
Copy link
Member

Oh i see, sorry I looked to quickly. Yes I'd say the default should be 0 for all by following what LAPACK does in *LANGE.

@charris
Copy link
Member

charris commented Feb 18, 2025

Mathematically, raising an error for the norm of an empty matrix makes the most sense. What is the use case for doing otherwise?

@carlosgmartin
Copy link
Contributor Author

carlosgmartin commented Feb 18, 2025

@seberg Suppose we are working in a complete lattice $L$.

The identity element of sup (join) $\vee$ is the least element (bottom) $\bot$ of $L$. Thus $\vee \varnothing = \bot$.
The identity element of inf (meet) $\wedge$ is the greatest element (top) $\top$ of $L$. Thus $\wedge \varnothing = \top$.

In the case of norms, $L = [0, \infty] \subset \mathbb{R}$ (the non-negative extended real numbers). Notice the presence of absolute values inside the expressions involving sup and inf. This is why it makes sense, in this context, to define the empty sup as $0$ and empty inf as $\infty$: These are the identity elements of sup and inf, respectively, in that $L$. If $L$ were $[-\infty, \infty]$ instead, then the empty sup would be $-\infty$.

Furthermore, the matrix $p$-norm is defined by

$$\|\mathbf{A}\|_p = \sup_{\mathbf{x} \in \mathbb{R}_d : \|\mathbf{x}\|_p \leq 1} \|\mathbf{A} \mathbf{x}\|_p = \sup_{\mathbf{x} \in \mathbb{R}_d : \|\mathbf{x}\|_p \leq 1} \left\|\sum_{i \in [n]} \mathbf{e}_i \sum_{j \in [d]} A_{ij} x_j\right\|_p = \sup_{\mathbf{x} \in \mathbb{R}_d : \|\mathbf{x}\|_p \leq 1} \|\mathbf{0}_n\|_p$$

where the last equality holds if $n = 0$ or $d = 0$.

The norm of a zero vector is, by the definition of a norm, zero. (Note that ord < 1 is not strictly a norm, as noted in the docs.)

The fact that the norm of a zero matrix (of which every empty matrix is an example) is zero is also explicitly stated here.

@melissawm

some of these norms are the result of computations e.g. "the sum of singular values". In that case it doesn't seem like assert_equal is the right choice for these tests.

The additional tests are for when $A$ is empty, for which no actual arithmetic is involved, and for which the result is (and should be) always exactly 0 or inf, so exact comparison makes sense in this context.

@charris

Mathematically, raising an error for the norm of an empty matrix makes the most sense.

Hard disagree. Mathematically, the answer is well-defined. That'd be like saying that an empty sum or product should raise an error rather than return 0 or 1 respectively.

What is the use case for doing otherwise?

This is the example that inspired me to create this PR. Other examples might be in the issues I linked to.

@charris
Copy link
Member

charris commented Feb 18, 2025

That'd be like saying that an empty sum or product should raise an error rather than return 0 or 1 respectively.

Yes, that is a useful convention for sums and product because of how they are used. But a norm, while it may use sums and products, is not just a sum, it has a meaning. Does the height and width of non-existent object make sense? You can return numbers for both, but what purpose would it serve? If your use case is to answer the question "will it fit in my closet", then returning zero for both makes sense. That is why I am asking about the use case.

@carlosgmartin
Copy link
Contributor Author

carlosgmartin commented Feb 18, 2025

@charris Geometrically, $\|A\|_p$ measures the longest $p$-radius of the image of the $p$-norm unit ball under $A$. This is used, for example, to compute Lipschitz constants of objective functions for optimization purposes.

@jakevdp
Copy link
Contributor

jakevdp commented Feb 18, 2025

@charris I think it's helpful to point out that in many cases NumPy already returns 0 for an empty norm; e.g.

In [1]: import numpy as np

In [2]: np.__version__
Out[2]: '2.2.2'

In [3]: np.linalg.norm(np.zeros((0, 0)))
Out[3]: np.float64(0.0)

In [4]: np.linalg.matrix_norm(np.zeros((10, 0, 0)))
Out[4]: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])

So the contribution here is not so much about adding new behavior, but rather about making already existing behavior more consistent.

@seberg
Copy link
Member

seberg commented Feb 18, 2025

I don't think you can argue with the supremum/infinum being 0 or all inf on the set of all possible values for the norm.

Unless you can argue that 0 is the correct answer for all (sensible) subsets of matrices, however chosen, having a default return is just wrong. That is exactly the same reason why maximum/minimum raise an error on empty arrays (The user might be working only with values within a specific subset, i.e. range, certainly inf is a special value so it isn't actually included in any sensible set to begin with!).
For example:

  • The set of matrices of any shape that contains only 0s always has 0 norm. So, IMO, clearly 0 is the only possible choice that might be right.
  • But say, the set of a matrix containing a single constant value? If something is defined via a sum([]) then I argue say that is OK, because the sum is const * N and for N=0 that is zero. In general, sum([]) is well defined as 0, so the extension makes sense to me; although I am happy if someone disagrees w.r.t. matrix norm.
    (we ignore matrices containing float inf here, because it is not a mathematical real value.)
  • For some norms, other sets of matrices might be interesting, although I suspect just the "constant value" may be enough.

@carlosgmartin
Copy link
Contributor Author

@seberg The mathematical definition of a norm forces the norm of an empty matrix to be zero. That particular fact has nothing to do with subsets one is working in.

@seberg
Copy link
Member

seberg commented Feb 18, 2025

How can it be that you have assert_equal(np.linalg.matrix_norm(x, ord=-2), np.inf)?

@carlosgmartin
Copy link
Contributor Author

@seberg ord < 1 aren't norms. They need to be considered separately.

According to the docs, ord=-2 is the smallest singular value. The infimum of an empty set is the greatest possible value, which is infinity.

@jakevdp
Copy link
Contributor

jakevdp commented Feb 18, 2025

According to the docs, ord=-2 is the smallest singular value. The infimum of an empty set is the greatest possible value, which is infinity.

This point is similar to the discussion in #5032, where the decision was to not return the supremum/infimum in the case of empty min/max. Given the consensus there, I think we should not add such behavor to norm when ord implies a min or max reduction.

@carlosgmartin
Copy link
Contributor Author

@jakevdp You mean for the cases that aren't norms, right?

@carlosgmartin
Copy link
Contributor Author

If it'll expedite merging, I can restrict this PR to only the actual norms, which are uncontroversial and what I'm currently most interested in.

@seberg
Copy link
Member

seberg commented Feb 19, 2025

@carlosgmartin here is the point:

  1. IMO (and I think everyone else agrees), returning anything but 0 is definitely incorrect. Clearly we can't merge something we know is wrong.
  2. 0 is clearly the only possible choice and it is definitely very often a reasonable and useful choice. But the bullet-proof argument would be that is the correct choice in (almost) all situations (rather than known to be useful in yours).
    E.g. for min/max -inf/inf are very reasonable (numerical) and in many situations correct choices. But this PR is a good example of why it would be a bad default since here 0 rather than -inf is correct.

FWIW, I suspect it is in fine, and probably completely correct, to define it to be 0 for all proper norms. But right now the most helpful argument was Melissa pointing to prior art.

It may be true that this is mathematical convention or drops out of the definition of the norm. But while 0 being the only possible (and definitely often correct) choice, none of the references I found until now explicitly state anything about empty matrices.

@carlosgmartin
Copy link
Contributor Author

@seberg

none of the references I found until now explicitly state anything about empty matrices.

Carl de Boor. An empty exercise. ACM Signum Newsletter 25 (4), 2-6, 1990.

any norm of an empty matrix is zero

@carlosgmartin
Copy link
Contributor Author

I've edited the PR to restrict changes to the proper norms.

@seberg
Copy link
Member

seberg commented Feb 28, 2025

We discussed this briefly and we decided to give this a try. But it would be good to have release note maybe "change" just in case it surprises someone. (You can find info in doc/release/upcoming_changes/README.md.

@seberg seberg added the 56 - Needs Release Note. Needs an entry in doc/release/upcoming_changes label Feb 28, 2025
@carlosgmartin
Copy link
Contributor Author

@seberg Done.

Copy link
Member

@ngoldbaum ngoldbaum left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry for popping in after others have already looked this over. I think the release note could use a tweak but otherwise don't have a problem with merging this.

Co-authored-by: Nathan Goldbaum <nathan.goldbaum@gmail.com>
@@ -0,0 +1 @@
* In all cases, ``np.linalg.norm``, ``np.linalg.vector_norm``, and ``np.linalg.matrix_norm`` now return zero for arrays with a shape of zero along at least one axis. Previously, NumPy would raises errors or return zero depending on the shape of the array.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
* In all cases, ``np.linalg.norm``, ``np.linalg.vector_norm``, and ``np.linalg.matrix_norm`` now return zero for arrays with a shape of zero along at least one axis. Previously, NumPy would raises errors or return zero depending on the shape of the array.
* ``np.linalg.norm``, ``np.linalg.vector_norm``, and ``np.linalg.matrix_norm`` now returns zero for empty arrays for most norms chosen by the ``ord`` parameter (empty arrays have at least one axis of size zero). Previously, NumPy raised an error or returned zero depending on the shape of the array.

@ngoldbaum maybe you can have a quick think (suggest-apply-merge)? Or just merge as is if you prefer...

Actually, maybe the clearest is to just list the changed ones:

* The vector norm ``ord=inf`` and the matrix norms ``ord={1, 2, inf, 'nuc'}``
  now always returns zero for empty arrays
  [...]
  This affects `np.linalg.norm`, `np.linalg.vector_norm`, and `np.linalg.matrix_norm`.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, I tweaked your suggestion a teeny bit.

@mattip
Copy link
Member

mattip commented Mar 12, 2025

Merging. We can continue to improve the release note as needed.

Edit: rerunning a failing 32-bit CI test, it failed to start. If it passes I will merge.

@mattip mattip merged commit 0dbee7e into numpy:main Mar 12, 2025
69 checks passed
@github-project-automation github-project-automation bot moved this from Awaiting a code review to Completed in NumPy first-time contributor PRs Mar 12, 2025
@mattip
Copy link
Member

mattip commented Mar 12, 2025

Thanks @carlosgmartin

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
Development

Successfully merging this pull request may close these issues.

8 participants