Widen (inv)quad(!) signatures to allow StridedArrays

[ararslan]

Should help with JuliaStats/Distributions.jl#594.

[ararslan]

Based on my understanding of the logic here, widening the signatures to AbstractVector should be pretty harmless. Please let me know if that isn't the case or if there's another reason why we shouldn't go ahead with this.

[jmxpearson]

I tested it with a view, as in the automated test, and works for me.

One question, though: the quad! and invquad! methods still seem to require Matrix and StridedMatrix. I know that in some cases, we want this restriction to call BLAS, but how reasonable is it to push that requirement further down into the low-level methods?

Perhaps more to the point: what method does * call in quad when we don't have isa(x, StridedVecOrMat)?

[jmxpearson]

Is there a method *(a::PDMat, x::AbstractVector) or only *(a::PDMat, x::StridedVecOrMat)?

I get

julia> methods(*, (PDMat, AbstractVector))
# 1 method for generic function "*":
*(a::PDMats.PDMat, x::Union{Union{Base.ReshapedArray{T,1,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:DenseArray, DenseArray{T,1}, SubArray{T,1,A,I,L} where L} where I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex, Int64, Range{Int64}},N} where N} where A<:Union{Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:DenseArray where N where T, DenseArray}, Union{Base.ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:DenseArray, DenseArray{T,2}, SubArray{T,2,A,I,L} where L} where I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex, Int64, Range{Int64}},N} where N} where A<:Union{Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:DenseArray where N where T, DenseArray}} where T) in PDMats at /Users/jmxp/.julia/v0.6/PDMats/src/pdmat.jl:39

[ararslan]

Ah, you're right. It looks like the definition of * is

*(a::PDMat, x::StridedVecOrMat) = a.mat * x

where the mat field is <:AbstractMatrix. So I suppose one could end up with a multiplication that doesn't make sense even with x::StridedVecOrMat depending on what a.mat is. To that end, perhaps the signature for * should be widened as well, and we can just let the type mismatch case fall back to a MethodError in *?

[jmxpearson]

I'd like to have @andreasnoack weigh in, but that seems reasonable. It would let people know that * and \ are the minimal interface they'd need to define for everything else to work.

In this case, you'd also want to widen quad! and invquad! everywhere, I think.

[ararslan]

Looking through the internals of quad! and invquad! for the various types, all of the things that are called look safe to allow AbstractArrays, so that shouldn't be a problem. Thanks for your help here, @jmxpearson!

[jmxpearson]

Sounds good!

[ararslan]

I tested it with a view, as in the automated test, and works for me.

Do you mean with this change or without it?

Perhaps more to the point: what method does * call in quad when we don't have isa(x, StridedVecOrMat)?

Currently it's a MethodError (at least for PDMats), but with this change it falls back to a multiplication of the mat field with the other input, as in my inline comment above.

[andreasnoack]

Allowing StridedArray is sufficient to fix the in distributions, right? If so, I think we should stick with that. Indexing is now more complicated to get right for AbstractArrays since you'd have to consider the non-one-based indexing case. Using AbstractArrays is also much more likely to cause ambiguity problems.

[ararslan]

Allowing StridedArray is sufficient to fix the in distributions, right?

No, because in Distributions we're using a view of the matrix to avoid copying each column, and SubArray is not a subtype of StridedArray.

[andreasnoack]

No, because in Distributions we're using a view of the matrix to avoid copying each column, and SubArray is not a subtype of StridedArray.

julia> view(randn(2,2), :, 1) isa StridedVector
true

julia> view(randn(2,2), [1,2], 1) isa StridedVector
false

[ararslan]

...oh

[ararslan]

Thanks for setting me straight, @andreasnoack. I'll revert the Abstract changes and just widen the signatures for ScalMat and PDiagMat to Strided. I should be able to get to this later today.

[ararslan]

How do things look now, @andreasnoack?