Factorization is not a subtype of AbstractArray

[simonbyrne]

This means that I can't use a CholeskyDense or LUDense as the invertible matrix inside the Woodbury class.

However, changing the line in base/linalg_dense.jl:

abstract Factorization{T}

to

abstract Factorization{T} <: AbstractMatrix{T}

causes:

julia> X = randn(5,5); A= X'*X;

julia> C = chold(A)
5x5 Float64 CholeskyDense:
Illegal instruction: 4

[simonbyrne]

Ah, it's just due to the show method. Is there any easy way to make it a subtype, but keep the old behaviour of show without defining methods for each factorization?

[ViralBShah]

We should probably define a show method for each factorization anyways. It also seems logical to make Factorization a subtype of AbstractArray.

Cc: @dmbates

[ViralBShah]

Could you submit a pull request?

[JeffBezanson]

I'm actually not sure whether a factorization is an array. Does it have elements you can access by indexing? AbstractArray functions assume that.

[ViralBShah]

Yes, you can do so with dense factorizations, and possibly even sparse factorizations if you try very hard. Factorizations would only meaningfully support a part of the Matrix behaviour.

[ViralBShah]

The other alternative would be to modify the routines in Woodbury to expect Union(Matrix,Factorization).

[JeffBezanson]

I'm sure it's possible, but is that how they're meant to be used? Often a better solution is just to remove the type restriction rather than make something an AbstractArray just so you can call the function you want.

[ViralBShah]

I think creating a Union(Matrix, Factorization) is the way to go rather than Any, so that the interface is still well-defined. This kind of thing will come up in iterative solvers as well. For example, there are cases when an iterative method expects either a matrix or an operator (which is just a function) - and it wouldn't then make sense to make every Function to be an AbstractArray as well.

[timholy]

I'm actually not sure whether a factorization is an array. Does it have elements you can access by indexing?

No. A factorization is a a concept, not a representation, a point which is clearest when you realize that the factorization of different matrix types will have different representations. Moreover, take an example like LU decomposition. While a dense matrix's LU decomposition might be stored as a single chunk of memory of the same size as the matrix, it's actually representing both the L and U factors. I wouldn't know how ref should work in such cases.

I guess they're of type AbstractAbstractArray :-).

[toivoh]

Would it be fair to say that a Factorization is a linear operator that supports applying its inverse using \ and /?
I e conceptually Matrix <: Factorization <: LinearOperator (and perhaps concretely with multiple inheritance)

[timholy]

Yes, I think that's the true meaning.

@ViralBShah , I agree that Union is the way to go here.

[simonbyrne]

I guess I didn't quite understand the definition of AbstractArray. The manual is somewhat unclear on the topic, saying:

Operations on AbstractArray objects are defined using higher level operators and functions, in a way that is independent of the underlying storage class. These operations are guaranteed to work correctly as a fallback for any specific array implementation.

So does this mean that all operations (ref, +, *, etc.) should be defined for any subtype of AbstractArray? It is true that ref would be somewhat awkward, but a Factorization could always be converted to an Array if no native method was defined.

On a related note, perhaps Woodbury itself should be a Factorization?

[JeffBezanson]

AbstractArray is independent of specific representations, but it is still understood to be some kind of container for keeping elements in a grid. Examples of different representations are sparse and distributed.

If you want to do Array-style operations on a Factorization, then explicitly converting it to an Array first might be a good interface, since this has a cost and isn't a normal use of Factorizations.

There is also an issue open about the array type hierarchy in general, which needs work. We might need abstract types both more general and more specific than AbstractArray.

[dmbates]

The factors generic creates a tuple of matrices constituting the factorization from a Factorization object. Even with this generic, however, there are still cases where the result contains, e.g. permutation vectors, rather than the corresponding permutation matrix.

[simonbyrne]

AbstractArray is independent of specific representations, but it is still understood to be some kind of container for keeping elements in a grid.

Okay, that makes sense.

We might need abstract types both more general and more specific than AbstractArray.

I think that might be the best approach. So as to avoid AbstractAbstractArray, how about MetaArray?:

Array <: AbstractArray <: MetaArray
LU <: Factorization <: MetaArray
Woodbury <: MetaArray

[StefanKarpinski]

MetaArray doesn't really mean anything to me intuitively. I think we need to figure out what the essential properties are and that will tell us what it's called. A name like ArrayLike might be ok since factorizations are kind of like arrays but not quite arrays.

[simonbyrne]

True, meta is not quite right: quasi would be the more suitable latin prefix.

[nolta]

I like @toivoh's LinearOperator suggestion (although i would call it LinearTransformation).

[StefanKarpinski]

Yes, that's a good name too. It would be nice to see if we can figure out some other things that are in the class that are neither factorizations nor matrices. This is verging on the need for multiple inheritance since matrices are linear operators but tensors generally aren't.

[ViralBShah]

LinearOperator is a term I like, and generalizes very nicely. I was thinking of Operator, but that was too general. This can also include functions and such as we go forward. For example, a stencil could be a LinearOperator that is passed as a function to iterative methods, or any other place where such transforms are needed, without needing to explicitly doing vectorized operations manually every time to apply a stencil.

[StefanKarpinski]

Are stencils really linear operators in the mathematical sense?

[ViralBShah]

@simonbyrne Does the Woodbury issue still exist? Otherwise, we have really not had issues with Factorization objects, and we could close this.

[simonbyrne]

The direct issue seems resolved, so I'll close this issue, but there is still the larger issue of where linear operators fit in the type hierarchy.

[ViralBShah]

Yes the larger issue still remains.