Consistency of return type of special matrices

[andreasnoack]

When I wrote functionality for the Triangular and Hermitian types, the general opinion seemed to be that the types should only be used for dispatch and that all methods should return a Matrix.

I have just noticed that methods for Diagonal return Diagonal e.g.

julia> inv(Diagonal(randn(4)))
4x4 Float64 Diagonal:
diag:  -1.54012  2.12275  2.08832  -0.976705

whereas

julia> inv(Triangular(triu(randn(4,4))))
4x4 Float64 Array:
 24.2769  -24.7589   31.9379    -2634.95  
  0.0       1.07824  -0.981735    107.485 
  0.0       0.0      -0.709239     54.6777
  0.0       0.0       0.0         -29.9744

Which solution should we go for?

[stevengj]

My feeling is that if the user is employing special matrix types for performance/memory reasons, then the user is entitled to receive those matrices back from operations that preserve special-matrix properties.

[andreasnoack]

I think I agree, but I forgot to mention that the original reason not to return special matrix types was that it would require many new methods to be written, i.e. ensure that all kinds of matrix multiplication is defined for the new types.

[vtjnash]

I haven't followed that discussion too closely, but is it possible to insert a superclass CompressedMatrix or MatrixLike (perhaps one of the existing ones already serves this purpose) and define generic convert methods on them for operations that call LAPACK? For many other operations, I imagine that indexing is easy to compute and offers roughly the same efficiency. I probably am not the best person to comment however.

[andreasnoack]

@vtjnash I think that is the point in JuliaLang/julia#2345. It has been discussed several times. I don't know if there are technical reasons why StridedArray is not yet an abstract type. I have tried to change the super class in jltypes.c but I cannot make it work.

[ViralBShah]

We should certainly try and resolve JuliaLang/julia#2345. Perhaps start a branch?

[stevengj]

It would be great to resolve JuliaLang/julia#2345, but special matrices like Triangular would still be different than a StridedMatrix. A StridedMatrix would denote any matrix with the same memory layout as a SubArray (and hence one which could be used directly with most BLAS/LAPACK routines).

A Triangular matrix would not be of this type. It's not clear to me how the proposed MatrixLike type would differ from AbstractMatrix, in fact. (More operations that work on AbstractMatrix would be great, of course.)

[ViralBShah]

Regarding the point in this issue, I agree with @stevengj 's suggestion to return the same special matrix type when possible. We can implement more operations on special matrices as we go along, and until then, users will have to manually convert to Matrix for operations that are not implemented.

Also, currently, we do try to implement as many operations on AbstractArray as possible with the current type hierarchy. I haven't reviewed this recently, but if you have specific cases in mind that can be made more generic, do file an issue. Fixing JuliaLang/julia#2345 should make it possible to make more of the array codebase generic.

[andreasnoack]

@stevengj I misunderstood the purpose of JuliaLang/julia#2345. Still, my feeling is that it would be useful to have a DenseAbstractArray. At least when you convert Integer (special) arrays to Float64 (special) arrays and hopefully also in other situations. For example the definition in matmul.jl line 65 could be changed to DenseAbstractArray and then cover the special matrices as well.

[stevengj]

How does a DenseAbstractArray differ from an AbstractArray? What methods could it support that a regular AbstractArray cannot?

It seems like only point of having a new abstract array class is if the new abstraction exposes more information about the storage or capabilities of the matrix, e.g. that it is laid out in memory with regular strides, or contiguously with transposed storage, or... I'm not sure what is being exposed in this case.

[andreasnoack]

AbstractArray covers distributed and sparse arrays so I think that a DenseAbstractArray do expose storage capabilities. In the example, a new array is constructed and it would not be clear which type it should be if the method at line 65 was defined for AbstractArray.

[stevengj]

What storage capabilities does a DenseAbstractArray have that AbstractArray does not? Line 65 of what file?

[JeffBezanson]

To me a DenseAbstractArray would be one that can use algorithms that touch every element, while for sparse arrays you'd try to touch only nonzeros.

[ViralBShah]

What if we just have AbstractArray, and in the case of AbstractSparseMatrix overload all those methods that in AbstractArray that touch all elements? That would avoid having one more type. That is essentially what we are doing already, and seems to work ok.

[JeffBezanson]

I agree we don't necessarily need the type, just saying what I thought it meant.

[andreasnoack]

@stevengj Line 65 of matmul.jl (I gave the reference first time). My idea was that DenseAbstractArray should be used for LAPACK calculations. It is possible that I don't understand what you mean by storage capabilities.

[stevengj]

@andreasnoackjensen, in order to pass an array directly to most LAPACK or BLAS routines, it needs to be stored in memory in certain ways; basically, it needs to be a StridedArray (and often at least one dimension needs to be contiguous, ideally the first dimension although a contiguous second dimension can be handled with some implicit transpositions).

It cannot be a triangular or tridiagonal array in compressed format (except for a few specialized LAPACK routines which only accept that compressed format). It certainly can't be some abstracted data structure that hides its memory layout. For such data structures, your only choice is to make a copy before passing to LAPACK or BLAS, but in that case you might as well support any arbitrary AbstractArray.

[andreasnoack]

@stevengj The main reason for the special matrix types in Julia is to exploit the special LAPACK routines. By the way, Triangular and Hermitian do not use compressed formats. I don't see why copying should be necessary as long as the element type is BlasFloat. For example there is no copying in

julia> Asym=A+A'
4x4 Float64 Array:
 -0.6        1.07975    1.54888   0.395639 
  1.07975    0.278261  -0.60896   1.69276  
  1.54888   -0.60896    1.19782  -1.84686  
  0.395639   1.69276   -1.84686  -0.0994641

julia> eigfact!(Symmetric(Asym))
Eigen{Float64,Float64}([-2.3701568194781486,-1.8149528758493196,1.6447505826965325,3.3169781703376313],4x4 Float64 Array:
  0.63831     0.425256     -0.633801  -0.100068
 -0.0518874  -0.691626     -0.5831     0.423031
 -0.558721   -0.000552329  -0.453427  -0.694431
 -0.52697     0.583791     -0.229552   0.573408)

julia> Asym
4x4 Float64 Array:
 -0.514571   1.86865    0.905495  -0.0902633
  1.07975    0.181334   0.707043  -0.386195 
  1.54888   -0.60896    1.20932    2.53631  
  0.395639   1.69276   -1.84686   -0.0994641

To ask a specific question, would it make sense to have a method

*(A::AbstractMatrix,B::AbstractMatrix) = *(convert(Matrix,A),convert(Matrix,B))

?
That would make it easier to cover arithmetic for the special matrix types. The right hand side could also be chosen to be DArray or SparseMatrixCSC but as I understand @ViralBShah the suggestion is that the default AbstractArray should be Array.

Finally, all the conversions to BlasFloat before the calls to LAPACK, i.e. eigfact!(A::StridedMatrix) = eigfact!(float(A)) could defined for AbstractMatrix right? Then I wouldn't need to write that line for each special matrix type.

[stevengj]

@andreasnoackjensen, I didn't realize Symmetric and Hermitian don't use compressed storage. In that case, you're right that it makes sense to be able to pass them to ordinary LAPACK and BLAS routines, and JuliaLang/julia#2345 applies. (It doesn't apply to Triadiagonal and other banded-matrix types, which do use compressed storage, however.)

A DenseArray (#987 ... though I prefer some more descriptive name like ColumnMajorArray) (would not be "abstract" since we only use "abstract" when the storage format is unknown) would be then any matrix type with contiguous column-major storage. But this is more specific than what is needed for LAPACK: for LAPACK and BLAS we only need to have the contiguous storage for the 1st (column) dimension (e.g. ColumnMajorSubArray since this almost exclusively occurs for subarrays of column-major arrays).

As discussed in JuliaLang/julia#987, we should really have a hierarchy of abstract types with known memory layout: ColumnMajorArray <: StridedColumnMajorArray <: StridedArray <: AbstractArray. And probably also a RowMajorArray <: RowMajorSubArray <: StridedArray for contiguous row-major arrays, since these may be returned by external libraries (e.g. NumPy is row-major by default) and can be passed to LAPACK and BLAS without copying via some trickery (some routines take an explicit transpose argument, while for others you can implicitly solve a transposed version of the problem).

(Regarding your other question, I would prefer that * on AbstractMatrix not make a copy of the matrices if it doesn't have to. We shouldn't make copies of potentially large data structures unless the user explicitly requests it.)

[ViralBShah]

I like the idea to have ColumnMajorArray. Note that Array is not implemented in julia and is in the C code, and I believe that the codegen also has this assumption about Array. Some reworking of these parts may be required to clean up our array hierarchy and move towards indexing returning views.

[jiahao]

The recent discussion in JuliaLang/julia#987 suggests that we are leaning in favor of preserving the special properties of the output matrix where possible.

[andreasnoack]

Effectively, we have now decided to preserve structure whenever possible. There might still be some missing cases, but they should be handled in individual issues as they show up.