More matrix property checks

[ghbrown]

Description

I'd been tossing around this proposal for a while, but was spurred into action by @zoziha's recent issue which proposed is_square and is_symmetric checks.

This proposal suggests built-in matrix property checks for common properties like triangularity, hermiticity, etc. These kinds of checks are not difficult to write inline with a main program, but having them as a one line call with an obvious name would make any calling code shorter and more readable, and I believe such checks are common enough to warrant inclusion in the standard library.

In addition to first two proposed by zoziha, I am proposing a few extras. Below is my full list for now, but additions are welcome:

• is_square
• is_symmetric
• is_skew_symmetric
• is_triangular (upper, lower, diagonal (upper and lower), and false none results)
• is_hermitian
• is_hessenberg (upper, lower, tridiagonal (upper and lower), and false/none results)
  etc.

Prior Art

• Numpy does not appear to have an is_symmetric or any other of the proposed checks (after light digging)
• MatLab has an issymmetric that uses an option to also work as a skew symmetric test, it takes a matrix and returns a bool
• Julia has issym (here and see links for other is* tests), isdiag, istriu, istril, ishermitian, all of which take a matrix and return a bool

API Discussion
is_square, is_symmetric, is_skew_symmetric, and is_hermitian will probably be four separate tests that all take a matrix and return a bool*
*The one exception I might suggest is perhaps is_square having a form like is_square(A,n) that also sets n=dim(A). Personally I have found this extremely useful when checking inputs to my own liinear algebra routines for methods or decompositions in which only square matrices are permitted. I believe this form of is_square could be snuck into the interface.

The case of triangular and Hessenberg tests will perhaps require a bit more discussion. We could only provide the "exact" tests for upper and lower (and possibly diag like Julia) and let the user construct the logic "on their side", in which case our tests would take a matrix and return a bool. Alternatively, our functions could take a matrix and return a letter for one of the four possibilities. For triangularity, I imagine the outputs being something like 'u': upper, 'l': lower, 'd': diag, 'n': not triangular. Alternatively we could use 'b': both which would generalize for both triangular and Hessenberg but rely on the user to understand what it implies (for example 'b' for Hessenberg implies tridiagonal). I'm partial to the string returns, but I recognize that returning a character/string may not be the expected output of a function called is_*.

For all these checks, I think the key will be to find a solution that doesn't require too many unique stdlib functions for a user to remember, doesn't require the caller to do much of the work/logic/thinking on their side, while being easy and obvious in style/format/usage.

Once we decide on the API, I'm happy to implement these checks myself or split them among whoever is interested.

[zoziha]

These routines you mentioned will be useful and promising code implementations.
keurfonluu/Forlab(MIT license) has three implementations that hope to inspire you: is_member, is_square, is_symmertic. (see API)

[ghbrown]

Thanks for the link. I'll probably do some speed tests on different methods for some of the checks once we decide on API, but having a reference implementation is great.

As for is_member, I'm not sure if it fits in this proposal, but we can and should add it at some point.

[ejovo13]

This is a great proposal and would add more core features that are necessary to start filling out stdlib a little bit more. For efficiency's sake, I say go for it! Start implementing them. I'll be happy to review (although I don't have write access) when you've written some code and have working prototypes!

Godspeed!

[nshaffer]

How often does one need these tests? In my own work, I have certainly tested if a matrix is square, but that is also the simplest of these. I have also tested if matrices are symmetric/hermitian, but there the context was checking result of a calculation, where I was willing to accept small errors. The main use case I can think of for checking that a matrix is exactly symmetric would be if I wanted to pass the matrix in to a procedure that assumes it is symmetric, like the LAPACK dsy* family of routines. I think the main use case of is_triangular would be similar.

So I personally would not find much use for these routines, except maybe for debugging. However, I should not assume that my experience is typical. Maybe these kinds of checks are useful to other people? Or maybe there are use cases I have overlooked?

[ejovo13]

I also would not find much personal use for these routines but I think prior art and the rare case where you might have to call an LAPACK routine and you want to make sure that your input is in accord with the subroutines parameters. Who knows, maybe one day stdlib will incorporate clean LAPACK bindings and these functions could be used as checks!

Regardless, I think that these are nice, straightforward algorithms to implement. I also think it's great to encourage new contributions whenever possible, and this proposal seems legitimate and useful to me.

[nshaffer]

Rereading my comment, I sound too pointed. To be clear, I think these are fine functions to have in stdlib. My opinion on the API:

pure logical function is_{square,symmetric,skew_symmetric,hermitian,diagonal}(a)
  real, intent(in) :: a(:,:)
end function

pure logical function is_{triangular,hessenberg}(a, uplo)
  real, intent(in) :: a(:,:)
  character, intent(in) :: uplo ! should be 'u', 'U', 'l', or 'L'
end function

I am on the fence about intent(out) arguments for some of these. For instance, is_square(a, n) could optionally return the dimension as discussed upthread. Or is_symmetric(a, i, j) could optionally return the indices where it first detected that a(i,j) /= a(j,i). If we go that route, it should be implemented a separate functions under a single generic name, e.g.

pure logical function is_square_1(a)
  real, intent(in) :: a(:,:)
end function

impure logical function is_square_2(a, n)
  real, intent(in) :: a(:,:)
  integer, intent(out) :: n
end function

Allowing for returned arguments would make these functions more useful for debugging, but it also makes the implementation more complicated. So I slightly prefer not doing that for now to keep the PR size reasonable.

[ghbrown]

Thanks for the suggestions and pointers. This issue was to gauge use cases and interest/usefulness, so open discussion is more than welcome. I certainly don't want to propose and implement something useless and bloat stdlib.

@nshaffer If inexact checks are more useful in practice, perhaps they could go in another PR (if people want them). But yes, the plan was for these checks to be exact for now. I don't think your experience of use cases for these functions is atypical, but as you say, for prior art and convenience I think they are worth it. My experience may see this functions more than average though.

I prefer your API to the two options I presented at first, especially since it doesn't require the user to store character variables their side of the function. I'll give the implementation a shot with that format.

[zoziha]

I see that the is* routines in numpy are called logic functions (see numpy logic functions). I feel that although is_square and the like act on matrices, they also belong to logic functions.
I found that we have a need for numpy.isclose/allclose function. I am considering implementing it in stdlib, I think we can create a stdlib_logic module (stdlib_logic.f90) to include them?

[ghbrown]

Based on the entries listed in your link, I'm not so sure that the matrix is_* functions would fit in with other "logic functions", even though they obviously return truthy values. The problem of where to put these is made more difficult since Numpy doesn't have the matrix checks built in (that I could find).

Either place is fine with me, though, and moving it over should be pretty simple if we decide where to put after implementation. I should be able to start my PR this weekend with couple checks implemented. I'll put them in stdlib_linalg.f90 for now, since that makes the most sense to me. Glad you implemented the close functions, they'll be good to have.

[ghbrown]

My PR was merged quite a while ago, so I'm closing the issue to remove clutter. Feel free to reopen if needed.