Adding matrix norms to the `stdlib_linalg` module.

[loiseaujc]

Adding matrix norms to the `stdlib_linalg` module.

[loiseaujc]

We probably could. Julia is more exhaustive than SciPy for this. Having norm and opnorm allows them to have both the element wise p-norm (for abitrary p) and the induced one (for p=1, 2 and infinity). SciPy is restricted to the standard operator induced norms + the Frobenius one. I think it mostly depends on how exhaustive we want to be (and how often the non-Frobenius element wise norms are actually being used which I think is not much).

The more I think about it, the more it feels like the Julia way may actually cover more use-cases than are actually used in practice. I still like the difference between norm and opnorm though just to recall the user that these are fundamentally different mathematical objects. Having the two interfaces might also make it more transparent when a 2D array represent an actual matrix or a collection of vectors instead of having to play with the axis keyword.

[loiseaujc]

I've used norm and opnorm in the discussion because that is what is being used in Julia. opnorm indeed stands for operator norm. While using norm for vectors is pretty obvious, I don't have a strong opinion regarding opnorm (although I don't have a better alternative to propose).

[loiseaujc]

I took some time to think about all the different norms (both for vectors and matrices) I've ever needed in my scientific computing adventure over the past 15 years or so. Here is a pretty standard list (+ whatever scipy/Julia are offering).

Vector norms

• norm(x, 1): 1-norm of a vector (i.e. sum of absolute values),
• norm(x, 2): 2-norm of a vector (i.e. standard Euclidean norm),
• norm(x, inf): $\infty$-norm of a vector (i.e. maximum absolute value),
• norm(x, p): p-norm of a vector.

Among these four, the Euclidean 2-norm is by far the most popular and, more often than not, is what people mean when discussing the norm of a vector $\mathbf{x}$.

I've used the other three quite regularly as well but more from a convex optimization perspective. I hardly ever had to actually compute them but they're simple enough to implement and I believe are something everybody expects from any standard linear algebra module.

Matrix norms

• opnorm(A, 1): 1-norm of a matrix (i.e. maximum absolute column sum)
• opnorm(A, 2): 2-norm of a matrix (i.e. largest singular value)
• opnorm(A, inf): $\infty$-norm of a matrix (i.e. maximum absolute row sum)
• opnorm(A, "nuclear"): nuclear norm (i.e. sum of singular values)
• opnorm(A, "fro"): Frobenius norm.

Like everyone, I use the Frobenius norm all the time. Similarly, the 2-norm and the nuclear norm are quite extensively used in model order reduction. As before, the 1-norm and the $\infty$-norm, although I don't quite know where they are being used, are some things everybody would expect from a standard linalg module.

There are a bunch of other matrix-related norms that I have used or seen being used over the years but they are more of a niche thing. These include:

• schatten_norm(A, p): Schatten p-norm (i.e. the vector p-norm applied to the singular values of A).
• lognorm(A, 2): Log-norm of a matrix, also known as its numerical abcissa and defined as $\lambda_{\max} ( \frac{\mathbf{A} + \mathbf{A}^H}{2})$.

Among these two, the lognorm is probably the one I've used most often. It is related to the non-normality of the matrix and proves important when studying linear time invariant systems of the form $\dot{x} = Ax + Bu$.

---

Here are my thoughts (in no particular order) on the matter:

• Even if it goes against the scipy notation, I kind of prefer the norm/opnorm dichotomy as it makes norm(x, 2) and opnorm(x, 2) very clear from just reading the code and not having to check the declarations to see if x is rank-1 or rank-2 array.
• norm and opnorm should return the Euclidean and Frobenius norms by default, respectively.
• The schatten_norm and lognorm could be added although they wouldn't be very high in my priority list.
• As a starting point, the implementations should focus on general matrices. Specialized computations for say Hessenberg matrices, Circulant matrices etc could be included later on. Since they're all standard rank-2 arrays, we can't naturally dispatch based on the matrix type and, the simplest, would be to write a master subroutine with is_symmetric(A), is_hessenberg(A) etc dispatching to the specialized solvers. If you know of any better/more fortranic way, I'm all ears.
• On the development and mathematical sides, I would restrict norm to take only a rank-1 array as argument. This would make less code bloat to handle all the corner cases scipy is handling (i.e. whether the norm is applied on the whole array, only along the rows, or only the columns). On the other hand, this would force the user to write a for loop to compute the norm of a collection of vectors represented as a rank-2 array and would be a major departure from the numpy/scipy standard.

I tend to think that being strict about the rank-1 array argument for norm is a better scientific/programming practice than what scipy currently offers. It doesn't leave any room for interpretation. Forcing the user to write a for loop for a collection of vectors might also make the code more readable and less error-prone (although it may possibly cause some performance loss (?), I don't know).

@perazz : On a side note, do you have a publicly available roadmap for the development of the stdlib_linalg module? I have seen some discussions on the fortran-lang discourse but nothing centralized. Translating/adapting everything scipy.linalg is offering is too much for a single person, even a one-man team. I suppose however that we all have bits and pieces scattered around our different code bases (e.g. I have some sketch for schur, expm, hankel, toeplitz, etc).

[loiseaujc]

• the Fortran standard has norm2 that works for both arrays and matrices with the optional dim specification. We should find a way to be sure our notation is not easily confused with this.

I shall be able to start working on this by the end of next week. I'll start with the baseline implementation and will see from there how we can improve and make sure it is not confounded with norm2.

• We should find a unique way to express the norm type. For example, could it be that we require a character input for all norms, instead of a numeric one? I'm saying so because I'd personally like more having to deal with something like norm(x, '2'), norm(x, 'L2'), norm(x, 'inf') rather than having to do norm(x, ieee_value(0.0, ieee_positive_inf))). So I think the code would be clear for both the user and the developer:

character(*), optional, intent(in) :: norm_type
select case (norm_type)
   case ('2','Euclidean')
      ...
   case ('inf','Inf','Infinity')
     ...
end select   

I like this.

• For the (later) specialized matrix norms, I would think that an optional additional flag like is_hessenberg=.true. or matrix_form='Hessenberg' is probably the most performant design, because we can't really check if the matrix complies to any special properties at all times, it would be too expensive.

Agreed. I prefer the matrix_form = "Hessenberg" rather than is_hessenberg = .true.. I may be naïve but the second would incur a lot of optional args (is_hessenberg, is_symmetric, is_triangular, etc) of which only one can be set to .true.. The matrix_form = string seems easier to handle and if not specified, opnorm defaults to the computation using the general matrix approach.

[perazz]

Great idea @loiseaujc. Rn I'm trying to tackle the decompositions (pseudo-inverse, Cholesky, QR) and I was planning to address norms and condition number next. So, your contribution would be very welcome!

As a way to separate between the two approaches to norm, would there be a way to use the same norm interface for all, but with different arguments maybe? I was loosely thinking that something like norm(A,2) vs. norm(A,'Frobenius') could work.

[jalvesz]

@loiseaujc I basically agree with everything you said! I just have one question regarding the naming opnorm, is it for "operator norm" as to say the matrix is an operator on a vector space? could there be a more "explicit" name to make the distinction? in any case, totally agree that it is good to have 2 well distinguished interfaces. If not, opnorm is just fine.

Regarding your query about the type of matrix, I think stdlib has in the linalg module https://stdlib.fortran-lang.org/page/specs/stdlib_linalg.html there are several checks such https://stdlib.fortran-lang.org/page/specs/stdlib_linalg.html#is_hermitian-checks-if-a-matrix-is-hermitian, https://stdlib.fortran-lang.org/page/specs/stdlib_linalg.html#is_hessenberg-checks-if-a-matrix-is-hessenberg among others.

I tend to think that being strict about the rank-1 array argument for norm is a better scientific/programming practice than what scipy currently offers. It doesn't leave any room for interpretation. Forcing the user to write a for loop for a collection of vectors might also make the code more readable and less error-prone (although it may possibly cause some performance loss (?), I don't know).

Agree, and I don't think that in Fortran you would have a performance loss (that would be true in python). If the implementations and interfaces are well designed the compilers might even be able to properly vectorize nicely. One thing to consider though, would be to have at least dimensions 2, 3 and 4 with explicit unrolled implementations and then a generic one for d>4.

[perazz]

Thanks @loiseaujc for the detailed comments. I very mostly agree with your ideas and I would summarize the consensus so far:

• let's have two separate functions for vector (norm) and matrix (opnorm) norms
• let's start from general matrix types, and then we can extend them later.

Here is some further thoughts I'm having about the design.

• the Fortran standard has norm2 that works for both arrays and matrices with the optional dim specification. We should find a way to be sure our notation is not easily confused with this.
• We should find a unique way to express the norm type. For example, could it be that we require a character input for all norms, instead of a numeric one? I'm saying so because I'd personally like more having to deal with something like norm(x, '2'), norm(x, 'L2'), norm(x, 'inf') rather than having to do norm(x, ieee_value(0.0, ieee_positive_inf))). So I think the code would be clear for both the user and the developer:

character(*), optional, intent(in) :: norm_type
select case (norm_type)
   case ('2','Euclidean')
      ...
   case ('inf','Inf','Infinity')
     ...
end select   

• For the (later) specialized matrix norms, I would think that an optional additional flag like is_hessenberg=.true. or matrix_form='Hessenberg' is probably the most performant design, because we can't really check if the matrix complies to any special properties at all times, it would be too expensive.

[perazz]

I prefer the matrix_form = "Hessenberg"

Absolutely agree, and it would be more in line with the norm_type argument