Support for symmetric & upper/lower triangular matrices

[dgrnbrg]

In my matrix programming, I encounter both upper/lower triangular matrices (from Cholesky decomps) and symmetric matrices (representing covariances and other positive definite forms). Since I am interesting in adding support for indexing APIs to matrices, I'd like to add another enum & variant to packed matrices, in order to mark whether their "other half" is a mirror image or all zeros.

Would you accept this enum? This would additionally affect the implementation of the count nonzero functions.

[IvanUkhov]

Thanks for the issue!

Then Packed will start to deviate from LAPACK’s definition of a packed matrix. If I’m not mistaken, they don’t even have the notion of a symmetric matrix; they always work explicitly with triangular matrices when symmetric matrices are expected. However, I do agree that it’d be nice to model symmetric matrices explicitly. Perhaps, instead of augmenting Packed, one could add a new type that would internally contain an instance of Packed?

[dgrnbrg]

I think that Packed would still cover LAPACK's definition. For instance, the family of LAPACK functions SSPTRD | CHPTRD | DSPTRD | ZHPTRD (for more packed representations on the symmetric eigenproblems, see here for more packed symmetric problems) all treat the packed triangular format as a symmetric matrix, while the xPOTRF functions treat the packed triangular format as having zeros on the other side. (Yay consistency!) By adding another variant field to the Packed representation, the type system could help us catch mis-applications of the format.

Do you think that's a good enough reason to include the variant field?

[IvanUkhov]

Yes, this is exactly what I meant by “they always work explicitly with triangular matrices [even] when [technically] symmetric matrices are expected.” Anyway, we can try what you propose; I don’t really have a strong opinion here. About the implementation, do you want to

1. change Variant’s variants to something like SymmetricLower, SymmetricUpper, TriangularLower, and TriangularUpper

2. introduce a new enum field to Packed with two variants Symmetric and Triangular, or

3. introduce a new boolean field to Packed that tells if it’s symmetric or triangular?

The last one is ugly, and the second needs a good name, and it’ll probably require renaming the existing enum field variant and its type Variant.

[dgrnbrg]

I agree that the 3rd approach is ugly.

I initially considered the first approach, but there are many cases when we only care about whether the matrix is triangular or symmetric but not upper or lower (for instance, to validate the input to a function). If it's alright with you, I'd like to take approach 2, with the following enums:

enum Representation {
    Upper, Lower
}

enum Kind {
  Symmetric, Triangular
}

I just chose the names arbitrarily--luckily, the enum variants are what will show up all over the code, not the enum names :)

Do you like this approach?

[IvanUkhov]

There is another complication: Hermitian matrices. The additional dimension that the issue proposes to introduce doesn’t really fit well into this story of Packed. They’ve really thought it through and defined the packed storage the way they did. It has just enough information to cover the three cases without introducing any irrelevance and confusion; it’s the least common denominator. I’m again leaning toward what I originally wrote: perhaps special cases of packed matrices should be represented by their own types leveraging Packed, in which case even the type system will be helping to avoid potential misuses of the APIs and to eliminate additional checks.

[dgrnbrg]

I see what you mean re: Hermitian matrices--they would be a 3rd variant of the Kind enum.

To me, then, there's a greater abstraction issue: in order to build an indexable API, we need to be able to return values that are not in the range covered by the packed representation. Thus, the packed matrix couldn't have an index API, but also the Matrix trait is incorrect for many situations, since the number of nonzeros depends on the matrix representation.

What do you think about stripping that functionality out of the Matrix trait, and then making wrapper types for the other variants?

[IvanUkhov]

Thus, the packed matrix couldn't have an index API.

Since it’s more relevant to #3, I’ve replied to this comment there. Please have a look.

The Matrix trait is incorrect for many situations, since the number of nonzeros depends on the matrix representation.

I’m not sure what you mean. Can you please elaborate? When one implements Matrix, one knows how to count the number of the nonzero elements that are present in the matrix at hand.