add new vecnorm(x, p), generalizing and replacing normfro

[stevengj]

As discussed on the mailing list, this adds a new vecnorm(x,p) function (generalizing and replacing normfro) which computes a p-norm of any iterable container as if it were a vector.

Along the way, I also made a few other improvements, e.g. the generic-p case no longer allocates temporary arrays, accumulation is done in double precision for Float16 or Float32 arrays (but the result is still returned as Float16 or Float32), and cheapnorm in quadgk is eliminated in favor of vecnorm.

Also fixes #6056 (makes norm type-stable).

cc: @toivoh, @StefanKarpinski

[stevengj]

Also, note that I changed the various norm functions to accept only p::Real, rather than p::Number. Non-Real values for p made no sense to me whatsoever.

[JeffBezanson]

normfro should also be removed from exports.jl.

[stevengj]

Doesn't it have to be in exports because it still exists, albeit in deprecated form?

[StefanKarpinski]

No, the @deprecate macro automatically exports it for you (there's no point in having a private deprecation, after all).

[stevengj]

Ah, okay. Will fix.

[JeffBezanson]

Really nice.

[pao]

Can we get a NEWS.md snippet for the behavior change to norm() with a row vector as well?

[andreasnoack]

It is not necessary to do the conversion here. It should be done in vecnorm. This function is type stable without the conversion.

[stevengj]

It's just as easy to do the conversion here, and things are a bit more uniform if all the vecnormX functions return the same type.

[andreasnoack]

Actually, I considered the type instability of norm a feature, because I like that the 1- and inf-norms return integers for integer input, but I recognise the type inference problem. I've just realised that the old norm2(AbstractVector) suffered from that. However, I don't think it is necessary do the conversion before vecnorm.

@pao I did that yesterday.

[toivoh]

It can definitely be nice to get the 1-norm or inf-form of an integer array as an integer, but I believe that the type instability might be a real performance issue. Maybe we should add norm1, norminf, vecnorm1, and vecnorminf, for these cases? Though it feels a bit ugly that we should need to add four functions for it...

[stevengj]

I really don't see the practical utility of integer-valued norm functions. Vectors of integers don't form a Banach space... where would an integer-valued norm be useful?

[StefanKarpinski]

I really don't see the practical utility of integer-valued norm functions. Vectors of integers don't form a Banach space... where would an integer-valued norm be useful?

Full agree with this. Returning integer norms doesn't make any sense to me.

[StefanKarpinski]

Of all people, I think @stevengj is the last one we have to worry about not updating the docs.

[stevengj]

(I will add a NEWS item, but I usually add a NEWS item after the patch is merged.)

[andreasnoack]

Do vectors of floats form a Banach space? What has completeness to do with this? Are integers less real than floats? And what about the minus infinity norm?

[stevengj]

Vectors of floats are a reasonable approximation for vectors of reals, and are well-accepted as such. Vectors of integers are not. (Yes, integers are still subsets of the reals, but if you view them as a mere subset rather than as a sub-space then there is no need for the return type to be an Integer.) Of course, this doesn't mean their "norm" is not useful, but it eliminates all the usual algebraic justifications for a norm function — you need to find some other need for it in order to justify its existence.

Yes, we compute several things that are not norms. I was particularly skeptical of the utility of the p < 0 norms. But at least they are continuous functions for a certain subset of the reals, and don't require any extra code to compute (except for p=-Inf) so I suppose we might as well support them.

[pao]

Of all people, I think @stevengj is the last one we have to worry about not updating the docs.

The change to document isn't a part of this patch (it was already changed), so it seemed worth noting that as a separate item, particularly because it's a non-deprecable (deprecatable?) change.

[lindahua]

It is true that from a math standpoint that norms only make sense in Banach space.

However, allowing vectors of integers as input would be quite handy especially in interactive sessions. For example, I think it would be reasonable to expect the following to work:

vecnorm([1,2,3])

This may be interpreted as treating [1,2,3] as a real vector and computing the vector norm. I don't see this would violate the math correctness much.

[StefanKarpinski]

That's not the argument – the question is whether the 1-norm of an integer vector should be integral or floating point. And @stevengj and I are arguing that it doesn't make much sense to get an integral norm value even though in the case of the 1-norm the result happens to always be an integer.

[lindahua]

Ok. Then I agree with using floating point value for the norm.

[toivoh]

I actually do think that 1- and inf-norms make sense for integer matrices as well. Integer matrices are closed under matrix addition and multiplication, these vector norms give integer values on integer vectors, as do the induced norms on integer matrices, which in turn give exact bounds on the maximal amplification that an integer matrix will cause in the respective norms. Coincidence or not, there is clearly some mathematical structure here.

That said, integer 1- and inf-norms are probably a niche enough application, and it's easy to roll your own.

[stevengj]

@toivoh, much of the concept of "linear algebra" goes out the window for integer matrices and vectors, because integer vectors are not a vector space (over the integers) because integers do not form a field (with the usual + and *). Yes, you can still define certain operations, but the practical utility of those operations is greatly limited compared to the situation where you have a vector space (and especially a complete normed vector space, i.e. a Banach space).

The question is, are there any real-world applications of integer-valued L1 norms in which returning a Float64 norm value is not sufficient? Moreover, do they outweigh the performance advantages of having norm(x,p) and vecnorm(x,p) be type-stable? So far, the answer seems to be a resounding no.

[toivoh]

Oh, I definitely agree that norm and vecnorm should be type-stable. The
only question is whether we should expose specialized 1- and inf-norms as
well. I personally think that that would be nice, but I'm fine with not
doing it.

[stevengj]

The question of having different norm functions with various names was discussed back in #1875.

[StefanKarpinski]

In the context of wanting an integer answer for integer inputs, while what you're computing may happen to be the same as the 1-norm or the inf-norm, in such cases, I'd say you are doing abs-sum or abs-max reduction, rather than taking a norm. Exposing those computations as norms is weird.

[toivoh]

Yes, but that was before the type stability issue was brought up. Anyway,
please carry on.

[andreasnoack]

@stevengj If you insist that the input type represents the field of the vector space then think of the Julia type hate most: Rational{Int}.

However, my math arguments won't make the type inference problem go away.

[stevengj]

@andreasnoackjensen, I don't want to restrict the input type, I want to make the output type floating-point. We should certainly accept rational (and integer) vectors as input, but it is insane for the L1 norm (for instance) to try to return a rational output (= near-certain overflow for fixed-precision types, and exponential complexity for Rational{BigInt}).

Even when we can practically compute a rational or integer result from a norm function (e.g. for the Linf norm of a rational vector or the L1 norm of an integer vector), I don't think that we should, because the additional code required is simply not matched by any practical benefit. So far, not a single practical use of an integer- or rational-valued norm (vs. float-valued) has been raised.

(PS. While vectors over the rationals are indeed a vector space, the are not a Banach space because they are not complete.)

[StefanKarpinski]

Be it resolved that norm returns an appropriate floating-point type.

[stevengj]

So is the patch okay to merge?

[andreasnoack]

I think it is unfortunate to promote when there is no reason to do it. For norm(x,p) there is a type inference problem but not for vecnormInf or vecnorm1. No mathematical sound argument has been raised for promoting to floating point. What in the definition of a norm excludes integer valued output? You pretend that integers are not real numbers and I simply don't agree in that. I don't think that because 0.5+0.5 is a floating point 1+1 should be a floating point.

A cool thing about Julia is that you can write very generic code such that you can define a new Number type and compute with it. For the generic LU, the number type only has to be stable under division and for the generic QR it has to stable under normalisation with the Euclidian norm. Then someone can define a new type and just use the LU factorisation on that or another new type and do QR and eigenvalue calculations. I can't see that it makes the code more complicated, (maybe not assuming commutativity makes it slightly more complicated). The only complicated thing is to make the right promotion.

It would have costed you ten seconds of work to acknowledge that point and move the promotion to floating point from vecnormInf and vecnorm1 to vecnorm. Instead we have wasted out time on mentioning such things as Banach spaces, which complete irrelevant to this issue.

[stevengj]

You should never call them elsewhere because they don't even handle the isempty iterable case (since this way the isempty check occurs only once rather than in 5 places). Putting the float call in vecnormInf saves calling it (or similar complications) not just in vecnorm but in vecnorm2 and vecnormp. Given that, you might as well put float in vecnorm1 too for uniformity's sake, not to mention that the sum is done in floating point to prevent overflows.

I don't see the point in duplicating code, even if it is a few simple calls, and risking overflows to return a result no one wants from functions no one will call.

[alanedelman]

At the risk of jumping in without reading the whole thread --- I still hate the name
vecnorm for normfro -- because it forgets that normfro is indeed a matrix norm

i'm okay with both existing, and no depracating

my two cents