Matrix norms

[jiahao]

The matrix norm() function could use a redesign to reflect the many matrix norms that are used in numerical analysis. The usefully computable norms can be divided into two classes:

• Hölder p-norm of matrix for arbitrary p. This is the "usual" norm that is some sense the natural extension of vector p-norms to matrices. Currently, norm(A::Matrix, p::Number) throws an error if p is not 1, 2, or Inf. The computation of such p-norms is not usually presented explicitly, since it is NP-hard to compute; however, Higham's estimation algorithm appears to work well in practice. (doi:10.1007/BF01396242, Algorithm Pnorm, p. 548)
• Schatten p-norms, based on the singular values, which are simply norm(svdvals(M),p). This reduces to the Frobenius norm normfro() when p=2.

We could fruitfully introduce an extra symbolic argument to specify the desired type matrix norm, e.g. :fro, :holder, :schatten.

[stevengj]

I would tend to prefer just normschatten(M, p) rather than adding more symbolic arguments to norm (we went through this already in JuliaLang/julia#1875); we might as well add this since its implementation is trivial...although that is also an argument for not adding it. I'm skeptical that the Hölder norm is worth the effort of implementing (approximately) in Base.

[marcusps]

I am working on a package that is mostly aimed at a norm induced by a Schatten norm (i.e., completely bounded norm), which is non-trivial to compute, but I am throwing in the usual (easy to compute) Schatten norms in as well.

I have them defined as snorm(A::Matrix, p::Number), with aliases for trace norm (trnorm), Frobenius norm (fnorm), and nuclear norm (nucnorm).

I plan to make the package public this week.

I suppose Higham's algorithm for Holder p norms should make a good package.

[stevengj]

@marcusps, the usual convention in Julia would be to name the functions norm*, so that tab completion norm<TAB> would find them. (Although this would argue for normvec rather than our current vecnorm.)

[marcusps]

So vecnorm would be the exception? That is what I was emulating.

The norm* convention would get a bit confusing for Schatten norms (snorm -> norms), although I suppose normtr/normnuc (for p=1), normfro (for p=2), and normspec (for p=Inf) are clear enough.

[hayd]

Perhaps, similar to sort (where you pass an alg e.g. QuickSort/MergeSort) we could pass some norm-type to norm. Rather than have a handful of similarly named functions.

Perhaps vecnorm be deprecated in favour of something else.

[stevengj]

normschat or normschatten rather than norms or snorm? It probably doesn't make much sense to abbreviate too much.

[stevengj]

@hayd, see the discussion in JuliaLang/julia#1875 on passing an extra argument to norm.

I don't have any problem with renaming vecnorm and vecdot to normvec and dotvec or whatever; I don't have strong feelings about the name, but in retrospect beginning it with norm would be nice for tab completion. If we want to rename it, someone should submit a PR now before the 0.4 freeze.

[StefanKarpinski]

I have to say I'm not fond of the norm* naming convention. The tab-completion thing doesn't strike as sufficiently important to make the names totally unnatural. ?norm should still list all of them.

[hayd]

In JuliaLang/julia#1875 was about removing symbols as arguments, using types is really an alternative solution.

[StefanKarpinski]

Writing norm(L{2}, x) isn't so bad, and you could have const L2 = L{2} to make it norm(L2, x).

[marcusps]

norm(Lp,x) for vector and Holder p-norms, and norm(SchattenLp,A) for Schatten p-norms seem manageable. Should other aliases use specialized functions? I like trnorm, specnorm, cbnorm, etc for readability, but I can see how it borders on getting out of hand.

[stevengj]

I would hate to have L be reserved as a type name by Base. One-letter symbols are valuable.

I guess we could continue to use norm(x, p::Integer) for Lp (vector and Holder/induced) norms, but switch to norm(x, Frobenius) for normfro etc.

[StefanKarpinski]

As usual, we don't actually reserve these things – you can still use e and im as variables.

[KristofferC]

But not as a type?

[StefanKarpinski]

In any particular scope you can use the global version or use the name as a local, not both.

[tkelman]

I still think shadowing exported base names with local variable names is bad practice and should generally be avoided, and leads to some confusing errors once in a while (e.g. JuliaPackaging/WinRPM.jl#24 (comment)).

Could also dispatch on Val{:L2} now, which is a little ugly but not too terrible.

[marcusps]

I like the original suggestion in JuliaLang/julia#1875 to have other norms in a separate package, and then use the full names there. e.g.,

Norm.schatten(A,p)
Norm.frobenius(A)
Norm.nuclear(A)
Norm.cb(A)
Norm.diamond(A)

The ones that seem a bit odd are

Norm.operator(A)  # schatten(A,Inf)
Norm.trace(A) # schatten(A,1)
Norm.spectral(A) # == schatten(A,Inf) == norm(A,2)

but maybe there is an argument to drop all aliases, refer to them in the docs, and focus on the more general name, e.g., schatten vs trace, operators, nuclear, frobenius, etc. norms, and also one particular (distinctive) choice for each swath of aliases, e.g., nuclear(A) for all schatten(A,1).

It is a more verbose version of what @stevengj suggested, but I find it a bit more readable.

[oscardssmith]

Closing as we now have p-norms.

[stevengj]

Re-opening as we don't have the specific Hölder and Schatten matrix norms at issue here.

[StefanKarpinski]

Maybe good to open two separate issues, one for each norm and spell out the API that should be implemented?

[aravindh-krishnamoorthy]

Could we recap what exactly is missing here? I can quickly code the missing items then.

[stevengj]

It's not at all clear that we want any additional norms in Base or LinearAlgebra, as opposed to adding packages for more specialized norms.