`schurfact!` fails on `StridedMatrix`

[Jutho]

I noticed that the following fails:

B = randn(100,100);
A = view(B,1:50,1:50)
schurfact!(A)

The reason is that the Schur immutable assumes that T and Z are of the same type:

immutable Schur{Ty<:BlasFloat, S<:AbstractMatrix} <: Factorization{Ty}
    T::S
    Z::S
    values::Vector
    Schur(T::AbstractMatrix{Ty}, Z::AbstractMatrix{Ty}, values::Vector) = new(T, Z, values)
end

which is not the case, because T will be stored in the input matrix A (in this case a SubArray) whereas Z will be newly allocated in LAPACK.gees! using similar on A, which produces a regular Array.

I could simply fix this by widening the definition of the Schur type (similar to how Eigen is currently defined), but wanted to consult the linear algebra/lapack experts (my apologies for pinging @andreasnoack , @kshyatt , @tkelman , @jiahao ) to see if that makes sense and whether any other plans or changes are in the pipeline that I could contribute to.

For example, the type of values is abstract (Vector) so I could also change this (like it is in Eigen). However, the output of schurfact will still be type unstable because of this line in the output of LAPACK.gees!

all(wi .== 0) ? wr : complex(wr, wi)

Would it be that bad to just always return complex values, with the advantages of having type stability? I noticed the current behaviour also appears in the various function related to eigfact etc.

Another question/remark: I noticed that throughout LAPACK, new arrays are sometimes allocated with Array{$elty} and sometimes with similar. Given that the input array is always a StridedMatrix, I don't think similar can produce anything else than an Array (I checked for SubArray and for SharedArray).

With some input, I would be happy to prepare a PR.

[tkelman]

Would it be that bad to just always return complex values, with the advantages of having type stability?

I've wanted to make this change for a while. x-ref #13362 where even the number of outputs differs between real and complex, which I think is confusing. I haven't revisited that PR in a long time though.

[andreasnoack]

There are at least two somewhat separate problems here. I can see two possible solutions for the first problem of the types of T and Z. Each with some issues attached.

1. We could throw some extra type parameters at Schur but generally, it makes things simpler with fewer type parameters.
2. The other solution is to define something like

Base.convert{T,N,S<:AbstractArray,I,J}(::Type{SubArray{T,N,S,I,J}}, A::Array{T,N}) = view(A, indices(A)...

which would simply wrap Z in a SubArray when constructing the Schur. This is a little inefficient but shouldn't matter much for larger problems.

Regarding the element type, there might be a tradeoff depending on problem size. Type stability is important for smaller problems whereas avoiding complex arithmetic could be a saving for larger problems. However, from a random perspective, real-valued eigenvalues become pretty rare as the matrix size grows so maybe we should go ahead and make the change to always return complex results from Schur and general Eigen. However, then we definitely want split Eigen into two such that the Hermitian version can remain real.

I think that the similar -> Array change would be fine for the LAPACK wrappers although the change is also kind of unnecessary.

[Jutho]

Thanks for the input. For now, I will restrict to just making schurfact! work for subarrays using either your solution 1 or 2.

The rest should clearly be part of a larger plan that might require further discussion.

More type parameters in Schur would make it analogous to Eigen and would have less of a global impact than defining a new convert method, but I am fine with both solutions.

Anyway, I am currently writing my own pure julia schurfact! (for now just for matrices which are already in Hessenberg form); is there any interest in this?

[andreasnoack]

Eigen was actually one of the main reasons I mentioned the number of type parameters. It has annoyed me a couple of times that so many were needed. It might still be the right solution, though.

I've also been playing around with a generic implemention. I can't guarantee for the quality because it requires a lot of testing which brings me to your question. In theory, it would be great to have a generic schurfact but in contrast to the finite-step linear algebra routines (Cholesky, LU, QR) the spectral like routines will require much more testing. At the same time, we are having discussions about making Base smaller and migrating modules to a standard library so I think this functionality might be better in a package. At least until the standard library discussion has settled.

[Jutho]

Ok, if I get my version of schur working I'll take it to a discussion in your LinearAlgebra.jl package. As to the problem at hand, one reason I was asking about the return type of similar in the lapack routines is that, if the result is anyway always Matrix, one could for example define

immutable Schur{Ty<:BlasFloat, S<:AbstractMatrix} <: Factorization{Ty}
    T::S
    Z::Matrix{Ty}
    values::Vector
    Schur(T::AbstractMatrix{Ty}, Z::AbstractMatrix{Ty}, values::Vector) = new(T, Z, values)
end

such that at least the issue for Z is resolved without new type parameters.

[andreasnoack]

I'm not sure. Even though it would probably work for all relevant cases right now, it is not impossible to think of a distributed array schurfact that could use the same type.