Proposed change to mat64: drop support for arbitrary matrices

[Brendan Tracey]

The mat64 package provides:

1) A set of interfaces for a variety of matrix types (Matrix, Symmetric, Triangular)
2) Concrete types that implement those interfaces (Dense, SymDense, TriDense)
3) A set of interfaces for representing as the "underlying" matrix type (RawMatrixer, RawSymmetricer)

Methods and functions have interface arguments, allowing a number of different concrete types as inputs. This is beneficial, as it allows, say,

func (m *Dense) Mul(a, b Matrix)

a could be a Dense, b could be a SymDense, and then we could type-switch to call Dsymm to make the multiplication most efficient. Right now, we also nominally allow "arbitrary" matrix types, that is, matrix types which are not expressible as a BLAS matrix. There are use cases, for example representing a large matrix implicitly in some way. However, this additional support comes with a cost. First, we would like to have matrix operations as efficient as possible, and calling the "At()" methods directly is quite expensive. Thus we will have specific implementations for RawMatrix, RawSymmetric, but we will also need implementations where the Matrix cannot be represented. See for example the Add method on Dense, which has both the Dense version and the normal Matrix version. This sort of idea should be replicated in every method and function. Right now our support is haphazard; the symmetric methods currently require the matrix be a RawSymmetricer.

This extra support is not that big of a deal for simple functions like Add. However, we want to be able to provide good and fast implementations of advanced functions like Cholesky decomposition, QR factorization, etc. Most matrix implementations (numpy, matlab, etc.) call Lapack to compute these routines. Lapack itself builds off of BLAS, and we can recreate functions as necessary to support mat64. For example, I am close to having a working Cholesky decomposition. The problem is that the full Matrix cannot support coding in the Lapack style. Matrices are often overwritten in-place with their types changed (so a Symmetric input will be overwritten to be a Triangular output), which is not supported via the interface. Additionally, Lapack routines typically need to take matrix "views", which is not supported with the Matrix interface. These issues prevent us from just implementing a "Matrix" version of lapack (that uses At and Set instead of a direct access). Instead we'll need to code these by hand (as is currently done in QR). I think this will severely hamper our ability to provide advanced functions in the package, as we will need to implement both the efficient Lapack call and also figure out how to implement it by hand without using View.

 PROPOSAL:  we drop support for "arbitrary" matrix types, and instead require that the type be able to express itself in one of the "Raw" forms. That is, the interface for Matrix will stay the same, the function signatures will stay the same, but the type itself must implement RawMatrix, RawTriangular, or RawSymmetric (in the future expanding to the other types).

Benefits:
1) Significantly reduced coding burden
    - We are already inconsistent in how much support we provide
    - Much easier to implement routines for the new types
    - Dramatically easier to support advanced lapack-like routines
2) Significantly reduced testing burden
    - Our current test support mostly tests that the function works on the existing types, not the general types
3) It would be backward compatible to "add" support
    - If in the future we decide we want/need to support them, it's all added code

Drawbacks:
1) Less flexibility in the package  (don't allow implicitly defined, etc.)
  - I'm not sure how much demand there is out there, especially since specific cases will likely need special implementations anyway for speed


In short, I think this reduction in scope will help us provide a robust foundation to the package, as it reduces the implementation and testing burden. Once the package is full-featured, we can think about ways to expand support. There is still so much to do on mat64, and I'd rather us not have to worry about supporting uncommon cases yet.

Thoughts?

[Dan Kortschak]

I'm not particularly aligned with either decision direction. Though
there are some things to think about.

Comments in-line.

One of the things that needs to be considered is what happens when an
arithmetic method is given blas-based types that do not agree. At the
moment, this case is handled as the slow path (which is also the support
for the general case Matrix). Without this fallthrough, we need to come
up with a way to handle interactions between different blas-based types.
This might involve adding a wrapper type that is not-exported that acts
as a shim between these. We would still however have some fall-through
path.

> This extra support is not that big of a deal for simple functions like
> Add. However, we want to be able to provide good and fast
> implementations of advanced functions like Cholesky decomposition, QR
> factorization, etc. Most matrix implementations (numpy, matlab, etc.)
> call Lapack to compute these routines. Lapack itself builds off of
> BLAS, and we can recreate functions as necessary to support mat64. For
> example, I am close to having a working Cholesky decomposition. The
> problem is that the full Matrix cannot support coding in the Lapack
> style.

I don't think that we need to do this. At the moment, we have a Cholesky
that accepts a SymDense and an SVD that accepts a Dense - we require
that people call DenseCopyOf(m) if they are not providing an appropriate
type. For the more complex functions or functions where there is an
obvious reason to require a blas-based type, I think this approach is
fine.

> Matrices are often overwritten in-place with their types changed (so a
> Symmetric input will be overwritten to be a Triangular output), which
> is not supported via the interface. Additionally, Lapack routines
> typically need to take matrix "views", which is not supported with the
> Matrix interface. These issues prevent us from just implementing a
> "Matrix" version of lapack (that uses At and Set instead of a direct
> access). Instead we'll need to code these by hand (as is currently
> done in QR). I think this will severely hamper our ability to provide
> advanced functions in the package, as we will need to implement both
> the efficient Lapack call and also figure out how to implement it by
> hand without using View.

The way that lapack does this is intrinsically unsafe. If we want to
mimic it, unsafe will do that for us if we accept the comment above as
reasonable.

[Dan Kortschak]

One other thing to consider it how this impacts on a future sparse
matrix implementation.