WIP: typestable factorizations

[skariel]

There are several more factorization functions that can be made type-stable. If this is positively accepted I'l gladly change them too :)

[andreasnoack]

👍 However, we need to figure out if we want use types or instances for this, see #9452. I think I prefer the type version.

[simonster]

Is a stagedfunction really necessary here? Couldn't you do:

copy_oftype{T}(A::AbstractArray{T}, ty::Type{T}) = copy(A)
copy_oftype{T,N}(A::AbstractArray{T,N}, ty) = convert(AbstractArray{ty,N}, A)
qrfact{T}(A::StridedMatrix{T}) = qrfact!(copy_oftype(A, typeof(one(T)/norm(one(T))))

[andreasnoack]

The solution we find better should be used for all the promotion methods in the linear algebra code.

[simonster]

Since this is likely to be useful outside of Base, perhaps we could have:

copy{T}(::Type{T}, x::T) = copy(x)
copy{T}(::Type{T}, x) = convert(T, x)

Although the second method may be too broad, since convert to a new type is not guaranteed to copy everything that copy would.

[andreasnoack]

As I read the discussion, we should use the type version. It would be great if you could also prepare a similar change for the Cholesky as it also has a pivoted version.

[ViralBShah]

Before merging - it would be nice to squash the commits, and have an explanation on type stability in the commit message.

I faintly recollect a discussion on the mailing list too, but couldn't find it. If so, could someone who knows, link it here?

[skariel]

Sure, I'll do the following:

1. use copy_oftype instead of stagedfunction
2. type stability also for Cholesky factorization
3. squash the commits into a single one

There is some discussion on this type instability in the docs, is this what you mean, see here:
http://julia.readthedocs.org/en/latest/stdlib/linalg/?highlight=qrfact#Base.qrfact
(see note in blue below...)

[andreasnoack]

@skariel Great. I like @simonster's suggestion of extending copy to allow a type specification.

I think @ViralBShah is thinking of this discussion

[simonster]

This is in the Julia 0.3 release notes. I think we need a separate entry for these changes under Julia 0.4. Also it would be good to have deprecations in deprecated.jl.

[skariel]

lufact has a reversed default pivot argument ie pivot=true. So separating this into types means defining a type NoPivot and a cooresponging constant nopivot. Is this the way to go? the other options would be to use a different default, pivot=false like other functions do. But then this a a change in api that would require many changes all over the place. I actually started doing this but maybe I'll just revert.

What do you think?

[skariel]

there is no need actually for copy_oftype since there is no dependence of return type on value for this. So this is not a type instability. Only interesting are cases where return type depends on value of pivot. I'm changing this now

[simonster]

@skariel I'm pretty sure you need copy_oftype because Julia still considers the branch that's not taken in type inference (even though it could be inferred to be dead; see #5560), and if S != T the two branches give different types, e.g.:

julia> f{T}(A::StridedMatrix{T}) = (S = typeof(one(T)/norm(one(T)));S != T ? convert(AbstractMatrix{S},A) : copy(A));

julia> code_typed(f, (Matrix{Int},))
1-element Array{Any,1}:
 :($(Expr(:lambda, Any[:A], Any[Any[:S,:_var0,:_var1,:_var2,:_var3,:_var4,:_var6,:_var7,:_var9,:_var10,:_var11,:_var12],Any[Any[:A,Array{Int64,2},0],Any[:S,Type{Float64},18],Any[:_var0,Int64,18],Any[:_var1,Bool,18],Any[:_var2,Type{Int64},18],Any[:_var3,Type{Float64},18],Any[:_var4,Int64,18],Any[:_var6,Int64,18],Any[:_var7,Int64,18],Any[:_var9,Int64,18],Any[:_var10,Int64,18],Any[:_var11,Int64,18],Any[:_var12,Array{Int64,2},18]],Any[]], :(begin  # none, line 1:
        _var0 = norm(1)::Int64
        S = typeof((top(box))(Float64,(top(div_float))((top(box))(Float64,(top(sitofp))(Float64,1)),(top(box))(Float64,(top(sitofp))(Float64,_var0::Int64)))))::Type{Float64}
        _var3 = S::Type{Float64}
        _var2 = T::Any
        _var1 = _var3::Type{Float64} == _var2::Type{Int64}::Bool
        unless (top(box))(Bool,(top(not_int))(_var1::Bool)) goto 0
        _var4 = (top(arraysize))(A::Array{Int64,2},1)::Int64
        _var7 = _var4::Int64
        _var6 = (top(arraysize))(A::Array{Int64,2},2)::Int64
        return copy!((top(ccall))(:jl_alloc_array_2d,$(Expr(:call1, :(top(apply_type)), :Array, Float64, 2)),$(Expr(:call1, :(top(tuple)), :Any, :Int, :Int)),Array{Float64,2},0,_var7::Int64,0,_var6::Int64,0)::Array{Float64,2},A::Array{Int64,2})::Array{Float64,2}
        0: 
        _var9 = (top(arraysize))(A::Array{Int64,2},1)::Int64
        _var11 = _var9::Int64
        _var10 = (top(arraysize))(A::Array{Int64,2},2)::Int64
        _var12 = (top(ccall))(:jl_alloc_array_2d,$(Expr(:call1, :(top(apply_type)), :Array, Int64, 2)),$(Expr(:call1, :(top(tuple)), :Any, :Int, :Int)),Array{Int64,2},0,_var11::Int64,0,_var10::Int64,0)::Array{Int64,2}
        return copy!(_var12::Array{Int64,2},1,A::Array{Int64,2},1,(top(arraylen))(A::Array{Int64,2})::Int64)::Array{Int64,2}
    end::Union(Array{Int64,2},Array{Float64,2})))))

[simonster]

We could define !(::Pivot) = NoPivot() and use !pivot instead of nopivot.

[skariel]

yes, I see your point. I like the !pivot suggestion, lets use this

[skariel]

About copy_oftype -- it won't work for more elaborate cases I think, like this one for e.g.:

function cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0)
    S = promote_type(typeof(chol(one(T))),Float32)
    S <: BlasFloat && return cholfact!(convert(AbstractMatrix{S}, A), uplo, pivot = pivot, tol = tol)
    pivot && throw(ArgumentError("pivot only supported for Float32, Float64, Complex{Float32} and Complex{Float64}"))
    S != T && return cholfact!(convert(AbstractMatrix{S}, A), uplo)
    return cholfact!(copy(A), uplo)
end

Here (unless I'm mistaken of course :) only a staged function would help. So maybe better to use them for these cases?

EDIT -- I mean, maybe better to handle all cases equally ie using a stagedfunction

[andreasnoack]

@skariel I think the cholfact example you gave would go into at least two new definitions when using the Pivot type instead of a keyword. Therefore I think than @simonster's suggestion is still enough, but I might be wrong.

In #9452, the preference seemed to be for the type version, so I think we should use Pivot instead of the pivot instance.

Would Pivot{true} and Pivot{false} be a better option here?

[simonster]

@skariel Adopting @andreasnoack's suggestion of Pivot{false} and Pivot{true}, I would handle that as:

cholfact!{T}(A::StridedMatrix{T}, uplo::Symbol=:U, ::Type{Pivot{false}}) = # non-pivoting alg goes here
cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, ::Type{Pivot{false}}) = # pivoting alg goes here
cholfact!{T}(A::StridedMatrix{T}, uplo::Symbol=:U, ::Type{Pivot{false}}) =
    throw(ArgumentError("pivot only supported for Float32, Float64, Complex{Float32} and Complex{Float64}"))
cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U, pivot::Union(Type{Pivot{false}},Type{Pivot{true}})=Pivot{false}) =
    cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, pivot)

# probably also want this one...
cholfact!{T}(A::StridedMatrix{T}, uplo::Symbol=:U) = cholfact!(A, uplo, Pivot{false})

which gives the same error message for types for which pivoting is unsupported for both cholfact and cholfact!.

[edit: fixed method signature for cholfact to use a union type]

[skariel]

@simonster, I don't think your version is equivalent to the original one. For e.g. if using T==Type{BigInt} and pivot==false then the original version would dispatch cholfact!(convert(AbstractMatrix{S}, A), uplo) while if T==Type{Float32} it would dipatch cholfact!(convert(AbstractMatrix{S}, A), uplo, tol = tol). Watch the tolerance parameter, there is no distinction betwen the two cases in your version...

(in the T==Type{BigInt} case s==Type{BigFloat})

[simonster]

@skariel You can define:

cholfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, pivot::Union(Type{Pivot{false}},Type{Pivot{true}})=Pivot{false}; tol=0.0) =
    cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, pivot; tol=tol)
cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U, pivot::Union(Type{Pivot{false}},Type{Pivot{true}})=Pivot{false}) =
    cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, pivot)

Alternatively, everything could accept a tol parameter and cholfact! for non-BlasFloats could warn/error if it's non-zero.

[skariel]

@simonster the behavior is still not the same:

consider the original version given
pivot==false,
T==Type{Rational{Int64}} and
tol>0
it dispatches cholfact!(convert(AbstractMatrix{S}, A), uplo, pivot = pivot, tol = tol)
where S==Type{Float32}.

In your version, since Rational{Int64} is not a subtype of BlasFloat this would dispatch without specifying tol. This would cause a hard to catch bug.

This function has a complex dispatch on types which is hard (or impossible) to express with current multiple-dispatch semantics, this is likly why it was implemented in a type-unstable way before stagedfunctions were available

[simonster]

You're right. But I believe you can still get the current behavior by dispatch:

cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U, pivot::Union(Type{Pivot{false}},Type{Pivot{true}})=Pivot{false}; tol=0.0) =
    _cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, pivot, tol)
_cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo, pivot, tol) = cholfact!(A, uplo, pivot; tol=tol)
_cholfact!(A, uplo, pivot, tol) = cholfact!(A, uplo, pivot)

Obviously it would be good to have test coverage for this to make sure it actually behaves properly. And we may still want to do something besides silently drop non-zero tol parameters for non-BlasFloats.

I'm not sure using a staged function would actually make this much simpler. @andreasnoack would know more about the history of the current implementation.

[skariel]

Yes, that seems to work! I guess I'm still getting used to multiple dispatch and just didn't think of breaking a single function into three levels of dispatch. But it makes sense, it's equivalent to if equals else on types and that's all we have here. So I'll use this of course

[andreasnoack]

@skariel Right now when pivot=false, the tol argument is just ignored at a later point. I don't think it is a huge problem, but having a way of throwing an error when the user it trying to set tol for normal Cholesky could be a possibility.

The reason the dispatch is so convolved right now is partly that the tolerance only makes sense for pivoted Cholesky and partly that non pivoted Cholesky has a generic fall back and pivoted Cholesky is LAPACK types only.

There also used to be a separate cholpfact method for the pivoted version, but we decided to clean up the namespace a bit and only later realized that the type stability issue was too costly.

[andreasnoack]

@skariel Bump.

[skariel]

I'm here, I'll have time tomorrow, will try to push then

[skariel]

@andreasnoack I think this is it, all factorizations should be stable, all tests are passing. Is this solution good enough? If so ill start updating the metadata (documentation, news, etc.)

[skariel]

@andreasnoack sorry false alarm... I found a few things that still need doing. I'll try to finish it tomorrow

[skariel]

@andreasnoack it's good now, can you please check the interface is what we want, If so I'll just update the docs and that should be it

Thanks!

[skariel]

I think using Val{T} is good because it is a general solution. I just checked and having pivot a keyword argument seems no problem for type stability see attached image. So I suggest just changing Pivot{T} to Val{T} and merge. One open question: how should the union Val{ture} and Val{false} be named?

[skariel]

(closed by mistake)

[simonster]

@skariel That test case has a default argument (which is equivalent to an ordinary positional argument for type inference purposes), not a keyword argument, and the call to @code_typed is within the function rather than outside of it. This is not type-stable:

julia> f(; n::Union(Float64,Int)=1) = 1+n;

julia> g(x) = f(n=x);

julia> @code_typed g(1)
1-element Array{Any,1}:
 :($(Expr(:lambda, Any[:x], Any[Any[],Any[Any[:x,Int64,0]],Any[]], :(begin  # none, line 1:
        return (top(kwcall))(call,1,:n,x::Int64,f,(top(ccall))(:jl_alloc_array_1d,$(Expr(:call1, :(top(apply_type)), :Array, Any, 1)),$(Expr(:call1, :(top(tuple)), :Any, :Int)),Array{Any,1},0,2,0)::Array{Any,1})::Union(Float64,Int64)
    end::Union(Float64,Int64)))))

Also, I tried building this branch and I see:

julia> Base.return_types(qrfact, (Matrix{Float64},))
1-element Array{Any,1}:
 Union(Base.LinAlg.QRCompactWY{Float64,Array{Float64,2}},Base.LinAlg.QRPivoted{Float64,Array{Float64,2}})

julia> Base.return_types(cholfact, (Matrix{Float64},))
1-element Array{Any,1}:
 Union(Base.LinAlg.Cholesky{Float64,S<:AbstractArray{T,2},UpLo},Base.LinAlg.CholeskyPivoted{Float64,S<:AbstractArray{T,2}})

[simonster]

It seems this needs to be split into two identical methods for the different Pivot arguments to avoid a method ambiguity warning:

_cholfact{T<:BlasFloat}(A::StridedMatrix{T}, pivot::Type{Pivot{true}}, uplo::Symbol=:U; tol=0.0) =
    cholfact!(A, uplo, pivot = pivot, tol = tol)
_cholfact{T<:BlasFloat}(A::StridedMatrix{T}, pivot::Type{Pivot{false}}, uplo::Symbol=:U; tol=0.0) =
    cholfact!(A, uplo, pivot = pivot, tol = tol)

[skariel]

yep you're right. I'll change Pivot{T} to Val{T}, and the pivot keyword argument to a positional one. Also I'll drop the PivotUnion in favor of an in-place union this way we get one thing less to name

[skariel]

more type instability:

function (\)(A::StridedMatrix, B::StridedVecOrMat)
    m, n = size(A)
    if m == n
        if istril(A)
            return istriu(A) ? \(Diagonal(A),B) : \(LowerTriangular(A),B)
        end
        istriu(A) && return \(UpperTriangular(A),B)
        return \(lufact(A),B)
    end
    return qrfact(A,eltype(A)<:BlasFloat?Val{true}:Val{false})\B
end

[skariel]

I mean, all this can be fixed, why even allow this? is there a single instance in the standard lib where type instability is actually needed? I'd like to see that

[skariel]

Tests seem to fail because of a bug, I've filed a new issue #9938

[skariel]

@simonster this is ready. I'm just waiting for #9938 to squash

[StefanKarpinski]

@skariel: I mean, all this can be fixed, why even allow this? is there a single instance in the standard lib where type instability is actually needed? I'd like to see that

Why is Julia a dynamic language rather than a static one? Kind of a deep question to sneak into a comment on a pull request. Largely because dynamic languages are easier to use. A static language assigns types to all expressions based on the structure of your program. This means that there are rules that are part of the language – pretty complicated ones, inevitably – for how types get assigned to expressions. If you write a program that doesn't obey those rules, then it won't compile, let alone run. This means that programmers have to learn those rules before they can write programs that will run. In dynamic languages, on the other hand, there are no such rules, which is a massive reduction in the learning curve before you can use the language. You can write a Julia program that is wildly type unstable and it will work. If you wonder why it's not as fast as you would like, then there's more to learn, but at least you can get a program working in the first place.

[StefanKarpinski]

I also believe that \ is precisely one of those cases where type instability is intentional.

[andreasnoack]

@skariel I thought for a long time that one of branches in S==T ? typeA : typeB were eliminated when possible which is part of the reason why there are so many type instabilities in the linalg code. The other reason is that, in some cases, we have decided to be type unstable, e.g. when we removed the cholpfact function in favor of keyword in cholfact we were fully aware of the type stability issue. I think we underestimated the performance impact of this so I this it is good that some if this is reverted by your fixes here.

There are a couple of problematic cases:

• The factorize function has been mentioned a couple of times as an example of a function which cannot be type stable.
• In cholmod.jl there will be a function that can read files in the MatrixMarket format. It dynamically determines if the elements are real or complex, the integer size needed and if the matrix is symmetric/Hermitian.
• The eigen functions for non-symmetric problems. Here we could make it type stable, but it would require us to always return complex results which we might don't want to do.

However, I don't see a reason why Matrix{T}\VecOrMat{S} should be type unstable. Though, there is a type stability problem for SubArray rhs, which is something we'll have to consider quite generally in the linalg code.

I think it is important to be specific about for which types there is a type instability. Many Julia function have a very broad signature and therefore it will probably be type unstable for some weird input type that doesn't really matter.

[skariel]

Yeah I already had a long discussion of this in the forums. Sorry for all the noise :)

@StefanKarpinski I enjoy daily the interactive workflow that Julia enables
@andreasnoack Dynamically loading files is a stong point in favor of type-instability. It could be done type-stable at the price of a more convoluted solution say using some ifs etc. which might be OK for library developpers

Anyway this is done and ready to merge (no new warnings, test passing, etc.)

[skariel]

(closed by mistake)

[timholy]

Though, there is a type stability problem for SubArray rhs, which is something we'll have to consider quite generally in the linalg code.

This sounds important; more details, please?

[andreasnoack]

@timholy

julia> @code_warntype eye(3)\sub(ones(3), 1:3)
Variables:
  A::Array{Float64,2}
  B::SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1}
  #s1994::(Int64,Int64)
  #s1992::(Int64,Int64)
  m::Int64
  #s1991::(Int64,Int64)
  n::Int64
  #s1993::Int64
  _var2::Int64
  _var3::Array{Float64,1}
  _var0::Int64
  _var1::Int64
  _var4::Int64
  _var5::Int64
  _var6::Int64
  _var7::Int64

Body:
  begin  # linalg/dense.jl, line 377:
      _var0 = (top(arraysize))(A::Array{Float64,2},1)::Int64
      _var1 = (top(arraysize))(A::Array{Float64,2},2)::Int64
      #s1993 = 1
      _var4 = _var0::Int64
      _var5 = (top(box))(Int64,(top(add_int))(1,1))
      m = _var4::Int64
      #s1993 = _var5::Int64
      _var6 = _var1::Int64
      _var7 = (top(box))(Int64,(top(add_int))(2,1))
      n = _var6::Int64
      #s1993 = _var7::Int64 # line 378:
      unless m::Int64 === n::Int64::Bool goto 3 # line 379:
      unless istril(A::Array{Float64,2})::Bool goto 1 # line 380:
      unless istriu(A::Array{Float64,2})::Bool goto 0
      _var2 = (top(arraysize))(A::Array{Float64,2},1)::Int64
      _var3 = getindex(A::Array{Float64,2},diagind(_var2::Int64,(top(arraysize))(A::Array{Float64,2},2)::Int64,0)::StepRange{Int64,Int64})::Array{Float64,1}
      return $(Expr(:new, Base.LinAlg.Diagonal{Float64}, :(_var3::Array{Float64,1}))) \ B::SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1}::Union(Array{Float64,1},SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1})
      0: 
      ((top(getfield))((top(getfield))(Base,:LinAlg),:chksquare)::Any)(A::Array{Float64,2})::Int64
      return $(Expr(:new, Base.LinAlg.LowerTriangular{Float64,Array{Float64,2}}, :(A::Array{Float64,2}))) \ B::SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1}::Union(Array{Float64,1},SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1})
      1:  # line 382:
      unless istriu(A::Array{Float64,2})::Bool goto 2
      ((top(getfield))((top(getfield))(Base,:LinAlg),:chksquare)::Any)(A::Array{Float64,2})::Int64
      return $(Expr(:new, Base.LinAlg.UpperTriangular{Float64,Array{Float64,2}}, :(A::Array{Float64,2}))) \ B::SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1}::Union(Array{Float64,1},SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1})
      2:  # line 383:
      return __lufact#181__(true,A::Array{Float64,2})::Base.LinAlg.LU{Float64,Array{Float64,2}} \ B::SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1}::Union(Array{Float64,1},SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1})
      3:  # line 385:
      return (top(kwcall))(call,1,:pivot,Float64 <: BlasFloat::Bool,qrfact,(top(ccall))(:jl_alloc_array_1d,$(Expr(:call1, :(top(apply_type)), :Array, Any, 1)),$(Expr(:call1, :(top(tuple)), :Any, :Int)),Array{Any,1},0,2,0)::Array{Any,1},A::Array{Float64,2})::Union(Base.LinAlg.QRPivoted{Float64,Array{Float64,2}},Base.LinAlg.QRCompactWY{Float64,Array{Float64,2}}) \ B::SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1}::Union(Array{Float64,1},SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1})
  end::Union(Array{Float64,1},SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1})

[timholy]

Thanks, I understand better now. Not really an issue with SubArrays themselves (which is what I was worried about), but a genuine challenge in structuring the linear algebra code. For example, in this method the problem is it might return the SubArray A (if m or n is zero) or the Array C.
You could presumably always return C?

[andreasnoack]

Why are the underscore versions necessary? Couldn't the underscore versions just be used as the ordinary cholfact! methods instead of having cholfact! calling _cholfact!?

[skariel]

This would require changing the API as dispatching on Val{false} and Val{true} would require pivot to not be a default argument and hence change places with uplo. It would still require 3 functions, and since _cholfact are not exported anyway I preferred to not changed the API

What do you think?

[andreasnoack]

Ah. I hadn't thought about that. Thanks for the explanation. Let's stick with your solution. I'll merge.

[andreasnoack]

@timholy No. It has nothing to do with SubArrays. It is about the linear algebra methods. I dont' think they have been widely used with SubArrays yet so the type problems haven't been revealed. We should begin to include @inferred tests in the linalg test suite.