lufact(), qr() and svd() use similar() but expect Array

[andyferris]

I have tried to make StaticArrays.jl compatible with the functions in Base.LinAlg but I get hung up on some inconsistencies.

The functions lufact(), qr() and svd() all use similar() to construct a (mutable) array for working with and returning output. In StaticArrays that would be e.g. the MMatrix type for a matrix. However, internal types like Base.LinAlg.SVD have fields which must be of type Array.

I consider this a bug since, according to the docs, similar should "Create an uninitialized mutable array with the given element type and size" and is not bound to return Array. We should widen the array field types in SVD and LU using an extra type parameter like in Eigen, and do something to fix the problem with qr() (also, there may be other functions with similar problems).

Alternatively we could use Array constructors where they are required, but I think this is not preferable. Julia is meant to be a generic language.

Examples:

               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.5.1-pre+2 (2016-09-20 03:34 UTC)
 _/ |\__'_|_|_|\__'_|  |  Commit f0d40ec (13 days old release-0.5)
|__/                   |  x86_64-linux-gnu

julia> using StaticArrays

julia> m = @SMatrix [1.0 2.0; 3.0 4.0]
2×2 StaticArrays.SMatrix{2,2,Float64,4}:
 1.0  2.0
 3.0  4.0

julia> eig(m) # successfully returns MVector and MMatrix (not Array)
([-0.372281,5.37228],
[-0.824565 -0.415974; 0.565767 -0.909377])

julia> lufact(m)
ERROR: MethodError: no method matching Base.LinAlg.LU{Float64,StaticArrays.MMatrix{2,2,Float64,4}}(::StaticArrays.MMatrix{2,2,Float64,4}, ::StaticArrays.MVector{2,Int64}, ::Int64)
Closest candidates are:
  Base.LinAlg.LU{Float64,StaticArrays.MMatrix{2,2,Float64,4}}{T,S<:AbstractArray{T,2}}(::AbstractArray{T,2}, ::Array{Int64,1}, ::Int64) at linalg/lu.jl:10
  Base.LinAlg.LU{Float64,StaticArrays.MMatrix{2,2,Float64,4}}{T}(::Any) at sysimg.jl:53
 in lufact!(::StaticArrays.MMatrix{2,2,Float64,4}, ::Type{Val{true}}) at ./linalg/lu.jl:20
 in lufact(::StaticArrays.SMatrix{2,2,Float64,4}) at ./linalg/lu.jl:72

julia> svd(m)
ERROR: MethodError: no method matching Base.LinAlg.SVD{T,Tr,M<:AbstractArray}(::StaticArrays.MMatrix{2,2,Float64,4}, ::StaticArrays.MVector{2,Float64}, ::StaticArrays.MMatrix{2,2,Float64,4})
Closest candidates are:
  Base.LinAlg.SVD{T,Tr,M<:AbstractArray}{T,Tr}(::AbstractArray{T,N}, ::Array{Tr,1}, ::AbstractArray{T,N}) at linalg/svd.jl:10
  Base.LinAlg.SVD{T,Tr,M<:AbstractArray}{T}(::Any) at sysimg.jl:53
 in #svdfact!#32(::Bool, ::Function, ::StaticArrays.MMatrix{2,2,Float64,4}) at ./linalg/svd.jl:19
 in (::Base.LinAlg.#kw##svdfact!)(::Array{Any,1}, ::Base.LinAlg.#svdfact!, ::StaticArrays.MMatrix{2,2,Float64,4}) at ./<missing>:0
 in #svdfact#33(::Bool, ::Function, ::StaticArrays.SMatrix{2,2,Float64,4}) at ./linalg/svd.jl:37
 in (::Base.LinAlg.#kw##svdfact)(::Array{Any,1}, ::Base.LinAlg.#svdfact, ::StaticArrays.SMatrix{2,2,Float64,4}) at ./<missing>:0
 in #svd#36 at ./linalg/svd.jl:57 [inlined]
 in svd(::StaticArrays.SMatrix{2,2,Float64,4}) at ./linalg/svd.jl:57

julia> qr(m)
ERROR: TypeError: typeassert: expected Array{Float64,2}, got StaticArrays.MMatrix{2,2,Float64,4}
 in (::Base.LinAlg.#kw##_qr)(::Array{Any,1}, ::Base.LinAlg.#_qr, ::StaticArrays.SMatrix{2,2,Float64,4}, ::Type{Val{false}}) at ./<missing>:0
 in qr(::StaticArrays.SMatrix{2,2,Float64,4}) at ./linalg/qr.jl:173

[andreasnoack]

So far, my thinking has been that the vectors of singular values and permutation would be okay to store in a Vector. The other parts are already parametric but we should try to make the vector fields parametric as well. It can become a little tedious, though.

[JeffreySarnoff]

This is one of some overly limiting type constraints that were introduced early in the development of Julia and then layered upon without relaxation. As a current, remediable, consequence, at least five packagings of number types and of number collective types now balk, unable to promulgate their designed-in willingness to just work seamlessly. Please deossify; embolden future kin to play well.

[Evizero]

I am trying to wrap my head around this, but it seems that a lot has changed since this issue was authored, because basically all error messages listed in the first post are now different on master.

To start of with a simple confusion. Let's say the type restriction for the vector is gone. How could the base function get the type of the vector given your matrix? similar only seems to work if I query it without additional parameters.

julia> m = @SMatrix [1.0 2.0; 3.0 4.0]
2×2 StaticArrays.SMatrix{2,2,Float64,4}:
 1.0  2.0
 3.0  4.0

julia> similar(m)
2×2 StaticArrays.MMatrix{2,2,Float64,4}:
 1.4822e-323   6.92545e-310
 6.92545e-310  6.92544e-310

julia> similar(m, 2)
2-element Array{Float64,1}:
 6.92545e-310
 6.92545e-310

The second aspect is that the methods you list expect a StridedMatrix as parameter, which I don't think you can opt-in outside Base (given that it is a typealias for a union) (Edit: I was wrong there).
It does seem to me that the implementations that don't outsource to LAPLACK could be more type-tolerant in that aspect though. For example: https://github.com/JuliaLang/julia/blob/master/base/linalg/lu.jl#L32 ?

[andreasnoack]

The second aspect is that the methods you list expect a StridedMatrix as parameter, which I don't think you can opt-in outside Base (given that it is a typealias for a union).

julia> type MyVec{T} <: DenseVector{T} end

julia> MyVec{Float64} <: StridedVector{Float64}
true

[Evizero]

Ah ok. Thanks for clearing that up.

So then to even potentially allow it to work with StaticArrays the following has to be made true in the package itself if I am not mistaken

julia> m = @SMatrix [1.0 2.0; 3.0 4.0]
2×2 StaticArrays.SMatrix{2,2,Float64,4}:
 1.0  2.0
 3.0  4.0

julia> typeof(m) <: DenseArray
false

[andyferris]

So then to even potentially allow it to work with StaticArrays the following has to be made true in the package itself if I am not mistaken

Right, I removed the DenseArray inheritance due to the fact that basically no method from Base could be made to work.

How could the base function get the type of the vector given your matrix? similar only seems to work if I query it without additional parameters.

This is specific to StaticArrays, since the size is a type parameter the similar method is only defined when the output size is type inferrable. I have a new Size type which acts as both a trait and similar to Val, so you can type similar(m, Size(2)) -> MVector{2}.

There are several other steps that need to be taken to get StaticArrays to work nicely with methods in Base - some of them require compiler improvements. However what I raised in the OP affects any type of AbstractMatrix defined outside base - even those with standard similar and so-on.

It does seem to me that the implementations that don't outsource to LAPLACK could be more type-tolerant in that aspect though. For example: https://github.com/JuliaLang/julia/blob/master/base/linalg/lu.jl#L32 ?

+1. Not sure why this kind of thing is there...

[andreasnoack]

Not sure why this kind of thing is there

What kind of thing?

[Evizero]

What I was trying to ask, but maybe sounded a bit too assertive (sorry, not intentional), is if the linked method really requires the matrix to be a StridedMatrix and not just a AbstractMatrix, since it only seems to use getindex and setindex and doesn't reason about the memory layout.

I would like to take a shot at this issue if it is still considered relevant.

I don't think I could solve all problems for StaticArrays though, but I came across some things I think could be worth doing while trying to figure out if this issue is something I could manage. For example this method here uses copy, while this method here uses similar, which makes a big difference for static arrays in terms of mutability. As a consequence the behaviour is fundamentally different depending on if the eltype is a float or an integer.

julia> copy(m)
2×2 StaticArrays.SMatrix{2,2,Float64,4}:
 1.0  2.0
 3.0  4.0

julia> similar(m)
2×2 StaticArrays.MMatrix{2,2,Float64,4}:
 6.93696e-310  6.93696e-310
 6.93696e-310  4.94066e-324

[andyferris]

What kind of thing?

What I was trying to ask, but maybe sounded a bit too assertive (sorry, not intentional), is if the linked method really requires the matrix to be a StridedMatrix and not just a AbstractMatrix, since it only seems to use getindex and setindex and doesn't reason about the memory layout.

Exactly.

And, yes, copy() of an immutable array at first glance doesn't sound like it would become mutable, but I guess the copy function only has utility for arrays which can be mutated (either the input or the output). Do you think copy(::AbstractArray) should always return an array which can be mutated?

EDIT: there is the Base.copymutable() function for that, so I think having copy return an immutable array should be allowed.

[andreasnoack]

The base linalg methods were written without any StaticArrays considerations. Mutability is generally assumed. Maybe it is possible to write the base methods such that they can work for StaticArrays. Generally, the methods in linalg are structured such that foo! does the actual linear algebra computation and foo just promotes sufficiently and makes sure that memory is fresh. It might be that StaticArrays will have to implement foo methods in order to get this right.

[andyferris]

Maybe it is possible to write the base methods such that they can work for StaticArrays.

Well, if we get something along the lines of #11902, then similar can make a Ref{SMatrix}, say, then we perform the mutating algorithm on that, and then finally we return something like publish(output) so that the Ref is removed from a static array (but publish(a::Array) = a).

(There was also a talk at Juliacon in the context of parallelism (distributed arrays) of the power of using a locally-mutating function to evaluate the output of functions which are then "published" to be considered more-or-less immutable. This "functional" approach to parallelism supposedly avoided issues with concurrency and so-on that arise when working with multiple workers writing and reading to mutable objects at the same time.)

It might be, that StaticArrays will have to implement foo methods in order to get this right.

That's the current approach, but it's a lot of work to duplicate all the functionality in Base!

[ViralBShah]

All the factorizations reported in OP work on StaticArrays now.

julia> using StaticArrays

julia> m = @SMatrix [1.0 2.0; 3.0 4.0]
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 1.0  2.0
 3.0  4.0

julia> using LinearAlgebra

julia> eigen(m)
Eigen{Float64, Float64, SMatrix{2, 2, Float64, 4}, SVector{2, Float64}}
values:
2-element SVector{2, Float64} with indices SOneTo(2):
 -0.3722813232690143
  5.372281323269014
vectors:
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 -0.824565  -0.415974
  0.565767  -0.909377

julia> lu(m)
StaticArrays.LU
L factor:
2×2 LowerTriangular{Float64, SMatrix{2, 2, Float64, 4}}:
 1.0        ⋅ 
 0.333333  1.0
U factor:
2×2 UpperTriangular{Float64, SMatrix{2, 2, Float64, 4}}:
 3.0  4.0
  ⋅   0.666667

julia> svd(m)
StaticArrays.SVD{Float64, SMatrix{2, 2, Float64, 4}, SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}([-0.40455358483375703 -0.9145142956773042; -0.9145142956773045 0.4045535848337568], [5.464985704219043, 0.3659661906262574], [-0.5760484367663209 -0.8174155604703631; 0.8174155604703631 -0.5760484367663209])

julia> qr(m)
StaticArrays.QR{SMatrix{2, 2, Float64, 4}, SMatrix{2, 2, Float64, 4}, SVector{2, Int64}}([-0.316227766016838 -0.9486832980505138; -0.9486832980505138 0.316227766016838], [-3.1622776601683795 -4.427188724235731; 0.0 -0.6324555320336751], [1, 2])