Proposal: cleanup of dense and structured matrix types

[jiahao]

The goal of this proposal to replace the existing zoo of non-sparse matrix types (i.e. not AbstractSparse) with the following 9 types:

• 4 types representing different storage formats used by LAPACK:
  • Matrix (Conventional storage)
  • Banded{n,m,order} (Band storage: a m x n banded matrix with bandwidth (kl, ku) is stored in a (kl+ku+1) x n matrix if order == :row and a n x (kl+lu+1) matrix if order == :col)
  • Packed (Packed storage)
  • RFP (rectangular full packed; see LAWNS 199)
• 4 types representing optional symmetry labels (no symmetry label implies no symmetry) with an uplo label referencing whether the upper or lower part of the matrix is stored and sometimes an isunit label referring to whether the diagonal is all ones or not.
  • Hermitian{T,AbstractMatrix{T},UpLo}
  • Hessenberg{T,AbstractMatrix{T},UpLo}
  • Symmetric{T,AbstractMatrix{T},UpLo}
  • Triangular{T,AbstractMatrix{T},UpLo,IsUnit}
  • Trapezoid{T,AbstractMatrix{T},UpLo,IsUnit} (essentially a generalization of Triangular to rectangular matrices)

which I argue will be sufficient to wrap the vast majority of LAPACK functionality going forward.

My thoughts about this proposal:

1. This proposal reduces the number of existing structured matrix types from 9 to appropriate compositions of 3 symmetry label types (Hermitian, Symmetric, Triangular) and 2 storage types (Banded, RFP).

   Current	Proposed
   Bidiagonal{T}	Banded{1,0,:col} or {0,1,:col}
   Diagonal{T}	Banded{0,0,:col/:row}
   Hermitian{T,S<:AbstractMatrix{T},UpLo}	no change
   SymTridiagonal{T}	Symmetric{Banded{1,0,:col}}
   SymmetricRFP{T}	Symmetric{RFP{T}}
   Symmetric{T,S<:AbstractMatrix{T}}	no change
   TriangularRFP{T}	Triangular{RFP{T}}
   Triangular{T,S<:AbstractMatrix{T},UpLo,IsUnit}	no change
   Tridiagonal{T}	Banded{1,1,:col}

2. Redefining Bidiagonal, Tridiagonal and SymTridiagonal in terms of Banded provides a consistent interface for storing and retrieving super-/sub-/diagonals.

3. The Banded matrix type is closed under the ring of matrix addition and multiplication, which is not true of Union(Diagonal,Bidiagonal,Tridiagonal,SymTridiagonal). A cute side-effect of the algebraic closure is that the type parameters m and n specifying the number of sup-/sub-diagonals obey a (max-plus) tropical semiring algebra which mirrors the underlying matrix algebra:

   *(A::Banded{T,ma,na,order}, B::Banded{T,mb,nb,order}) -> Banded{T,ma+mb,na+nb,order}
   +(A::Banded{T,ma,na,order}, B::Banded{T,mb,nb,order}) -> Banded{T,max(ma,mb),max(na,nb),order}

   You can exploit part of the algebra of m,n by defining + using diagonal dispatch and a new promotion rule:

     promote_type(::Type{Banded{Ta,ma,na,order}, ::Type{Banded{Tb,mb,nb,order}}}) =
        Banded{promote_type(Ta,Tb), max(ma,mb), max(na,nb), order}
     +{T,m,n,order}(A::Banded{T,m,n,order}, B::Banded{T,m,n,order}) -> Banded{T,m,n,order)

   and similarly for -.

4. The full proposal allows representations of all 29 LAPACK matrix types (as of v3.5.0) by appropriate compositions of the aforementioned 9 types, and tuples where necessary to represent generalized problems:

   LAPACK code	Current Julia type	Proposed typealias
   BB	-	don't wrap (part of the cosine-sine decomposition)
   BD	Bidiagonal{:L or :U}	Bidiagonal{:uplo} is a typealias for Banded{1,0,:col} or {0,1,:col}
   DI	Diagonal	Diagonal is a typealias for Banded{0,0,:col) (the order parameter is mostly irrelevant)
   GB	-	Banded{n,m}
   GE(, OR, UN)	Matrix	-
   GG	(Matrix, Matrix)	-
   GT	Tridiagonal	Tridiagonal is a typealias for Banded{1,1,:col}
   HB(, PB)	-	Hermitian{Banded{n,m}}
   HE(, PO, PS)	Hermitian{T, Matrix{T}}	-
   HG	-	(Hessenberg, Triangular)
   HP	-	Hermitian{Packed}
   HS	-	Hessenberg{Matrix}
   SB(, PB)	-	Symmetric{Banded{n,0}}
   SF(, PF)	SymmetricRFP	Symmetric{RFP}
   SP	-	Symmetric{Packed}
   ST, PT	SymTridiagonal	Symmetric{Banded{1,0,:col}}
   SY, PO	Symmetric{T, Matrix{T}}	-
   TB	-	Triangular{Banded{n,0} or {0,n}}
   TF	TriangularRFP	Triangular{RFP}
   TG	(Triangular, Triangular)	-
   TP(, PP, OP, UP)	-	Triangular{Packed}
   TR	Triangular{T, Matrix{T}}	-
   TZ	-	Trapezoid{Matrix}

5. It is necessary to have two different representations of banded matrices since two memory layouts are required. The LAPACK storage format for band storage stores a m x n banded matrix with bandwidth (kl, ku) as (kl+ku+1) x n matrix; however, LAPACK's routines for tridiagonal/bidiagonal matrices expect to be passed the super-/sub-/diagonals as separate vectors. We cannot store these in the correct memory order using just Banded{m,n,:row}, but it is possible to do so with Banded{m,n,:col}.

   1. One could separate out Banded{m,n,:row} and Banded{m,n,:col) as two separate types (say Banded and BandedCol, but if we do that then there is no unique representation of Diagonal; both Banded{0,0} and BandedCol{0,0} would have equivalent memory mappings. One possible solution is to make the Diagonal constructor consistently create one of them. For functions computing on Diagonal matrices, it would be possible to make Diagonal a typealias of Union(Banded{0,0},BandedCol{0,0}). However, the Diagonal constructor still needs to preferentially select one or the other. For this reason it seems like Banded{m,n,:row/:col} is a better choice.
   2. The Symmetric, Hermitian, Hessenberg and Triangular labels modify the storage requirements; composing these types with a Banded(m,n) matrix makes either m or n semantically meaningless and used only to describe the storage of irrelevant elements.

6. I've specifically ignored the need to annotate matrices with spectral properties like positive semi-/definiteness. There is no need to introduce a special matrix type for the purposes of dispatching onto LAPACK routines of the P? type, since they can all be dispatched from Cholesky factorization objects.

   1. It is unclear, however, if there might be uses other than LAPACK dispatch for which such a label would be nesessary and there is no other reason to annotate a Julia matrix with a spectral label. (See JuliaStats/PDMats.jl ?)
   2. It's easy to add a PosDef or PosSemiDef type along the lines of Hermitian, for example, but then the composition of the two types of labels like Symmetric{PosDef{M}} or PosDef{Symmetric{M}} produces new types that are distinguishable by their order of composition, even though the order of labeling ought not to make a difference.

7. The OP, UP, OR and UN matrix types are special matrix types resulting from QR factorizations and are never used in any other context. Dispatching onto their appropriate LAPACK routines can be handled entirely from within factorization objects like QR* and Hessenberg and there should be no need to expose them as new matrix types. Similarly, BB is only ever used as part of the representation of the result of the cosine-sine decomposition, and can be dispatched on entirely through a hypothetical new CS factorization object.

8. This proposal turns matrix types like Hessenberg and Tridiagonal into immutable fixed size matrices, since constructing them as special cases of Banded{m,n} means that the type parameters m and n are completely sufficient to determine the size of the matrix, which means that the matrices must have known dimensions at the time of construction. This could be a major breaking change, although I don't have a good feel for how extensively these special matrix types are used in the wild. The type parameters give a minimum bound on the size of the matrix but are insufficient to determine the matrix dimensions exactly.

9. Ideally, UpLo and order should be enums with just two allowed values each ((:U, L) and (:row, :col)) but it seems much clearer to represent them with symbols instead of Bools since we don't have enums (yet).

10. I haven't worked out the details yet, but I suspect that this reorganization of the matrix types will lead to correspondingly cleaner matrix Factorization objects and type hierarchy.

The proposal here is the result of extensive discussion with @simonbyrne and @andreasnoack.

[tkelman]

Yes please. This would be a major improvement re: #8001. Doesn't address the brokenness of sparse at all, but one step at a time.

Banded will cover many of the dense structured types currently in lapack as special cases, and the algebra of banded matrix operations will be cool to see done right. For some reason I thought it had already been done in Julia, but apparently not.

For many of the reductions and binary operators I think adding a banded type will necessitate rethinking how the generic fallback implementations are done, in order for those operations not to be painfully slow for narrow-banded matrices. This will likely expose many of the same issues as currently exist for sparse matrices. We should perhaps consider an iterator-based protocol for more generic operations on matrices. That way the fallback for the iterator protocol could be a dense 1:n and we'd have the same behavior as today on full, dense matrices. But types like banded and sparse could provide specialized iterators to operate on each column or row.

[simonbyrne]

+ a lot

I found this nice guide to the LAPACK types.

The only note I really have is that Hessenberg HS is actually stored as a full matrix, not a banded one.

[jiahao]

I think adding a banded type will necessitate rethinking how the generic fallback implementations are done, in order for those operations not to be painfully slow for narrow-banded matrices.

True. Anecdotally from @andreasnoack's experiments it seems like the only real gains in going through LAPACK and BLAS are in routines that can take advantage of BLAS3 (matrix-matrix) kernels. It's difficult to see how BLAS3 can be used for banded algebra.

The only note I really have is that Hessenberg HS is actually stored as a full matrix, not a banded one.

Thanks, fixed.

[stevengj]

@jiahao, I don't understand this comment:

It is necessary to have two different representations of banded matrices since two memory layouts are required. The LAPACK storage format for band storage stores a m x n banded matrix with bandwidth (kl, ku) as (kl+ku+1) x n matrix; however, LAPACK's routines for tridiagonal/bidiagonal matrices expect to be passed the super-/sub-/diagonals as separate vectors

Since a (kl+ku+1) x n matrix has contiguous columns (assuming column-major storage), can't you pass separate "vectors" for the super-/sub-/diagonals simply by pointer arithmetic?

[jiahao]

Assuming I've read the diagram correctly, banded storage stores nonzeros column by column, so the sub-/super-/diagonals would be stored as vectors of stride kl+ku+1. Last I checked, the bidiagonal and tridiagonal routines do not support matrices represented as vectors of stride other than 1.

[stevengj]

Ah, I see; I misread it. This is annoying, although I suppose we could transpose the data on the fly when it is needed.

[andreasnoack]

My pull request #9779 relates to this issue.

The way I have been thinking about the uplo type parameter is that is useful when the output type of some methods on the matrix type depend on whether uplo=:L or uplo=:U, e.g. Triangular{:L}+Triangular{:L}=Triangular{:L}butTriangular{:L}+Triangular{:U}=Matrix`.

As far as I can see, this is not the case for Hermitian and Symmetric and therefore it is fine for those types to carry the uplo information as a field.

I think it is okay to make the assumption that a Hessenberg matrix is lower Hessenberg as I haven't seen an application of lower Hessenberg matrices. Hence, the type doesn't need an uplo parameter.

If it is necessary to introduce a trapezoid matrix type, I think we should follow the approach in #9779.

This will result in more types compared to the present solution, but it is easier to work with the versions without the parameters in a way that insures good type inference.

[mschauer]

When working with FixedArrays I noted that currently - and this is also what you propose - Diagonal is a storage and not a symmetry property (for example Diagonal{T} is parametrized with eltype, different from e.g. Triangular{T, AbstractMatrix{T}}.) A different possibility is to have a parametrization Diagonal{T,AbstractVector{T}} and have Diagonal behave closer to Triangular matrices.

[jiahao]

@mschauer this issue was focused on LAPACK, where Diagonal only has one representation. Are there any circumstances where Diagonal matrices are stored in more general AbstractVectors? I suppose Diagonal of SubArray could be made to work, but it has to be strided correctly in order for LAPACK to work with that data. Otherwise we'd have to dispatch onto more general fallback methods, which is possible but would complicate reasoning about performance.

[mschauer]

I was thinking especially about fixed-size diagonal matrices. But also the diagonals of subarrays you mention is a nice example. Further the 2d-Laplacian operator induces a tridiagonal matrix of tridiagonal matrices of which for example the Cholesky factor is a bidiagonal matrix of bidiagonal matrices. These these structured matrices of structured matrices are a recurring pattern. As example take the discrete Poisson equation, which could be directly solved by the abstract fallback version of the tridiagonal matrix solution algorithm.

[oscardssmith]

It's really a shame that we didn't get this for 1.0. Is there any way we could still do this without breaking everyone's code?