`sqrt(Hermitian(A))` errors (due to fallback) for complex `A`

[ericphanson]

I hope this isn't a duplicate; I couldn't find an issue though. Here's an example:

julia> A = @SMatrix(rand(2, 2)) + im*@SMatrix(rand(2,2))
2×2 SArray{Tuple{2,2},Complex{Float64},2,4}:
 0.911161+0.30874im   0.152706+0.0236032im
  0.55638+0.107542im  0.175589+0.100683im 

julia> B= A*A'
2×2 SArray{Tuple{2,2},Complex{Float64},2,4}:
 0.949411+0.0im        0.569345+0.0625584im
 0.569345-0.0625584im  0.362093+0.0im      

julia> sqrt(B)
2×2 SArray{Tuple{2,2},Complex{Float64},2,4}:
 0.859886+0.0im       0.455524+0.050052im
 0.455524-0.050052im  0.389982+0.0im     

julia> sqrt(Hermitian(B))
ERROR: setindex!(::SArray{Tuple{2,2},Complex{Float64},2,4}, value, ::Int) is not defined.
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] setindex!(::SArray{Tuple{2,2},Complex{Float64},2,4}, ::Float64, ::Int64) at /Users/eh540/.julia/packages/StaticArrays/VyRz3/src/indexing.jl:3
 [3] macro expansion at /Users/eh540/.julia/packages/StaticArrays/VyRz3/src/indexing.jl:51 [inlined]
 [4] _setindex!_scalar at /Users/eh540/.julia/packages/StaticArrays/VyRz3/src/indexing.jl:39 [inlined]
 [5] setindex! at /Users/eh540/.julia/packages/StaticArrays/VyRz3/src/indexing.jl:35 [inlined]
 [6] sqrt(::Hermitian{Complex{Float64},SArray{Tuple{2,2},Complex{Float64},2,4}}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/symmetric.jl:832
 [7] top-level scope at none:0

As you can see, it works without the Hermitian wrapper, and it also works in the real case. I will just not use the Hermitian wrapper, but I thought maybe the problem should be noted somewhere.

[andyferris]

Basically we haven't provided overloads for all the functions in LinearAlgebra for all the relevant types... it's typically a matter of chipping away at this.

[fipelle]

+1

It would also be great to add support for Symmetric. I am working on a new package for forecasting and having it would be great to speed up many calculations.

[c42f]

We'd be really happy to merge some code to make Symmetric (and Hermitian) work better in various cases. To give some quick hints on how this could be fixed (we accept PRs ;-) ):

• Firstly, we probably want to define sqrt(m::StaticMatrixLike{T}) to dispatch to an internal _sqrt implementation rather than dispatching on m::StaticMatrix.
• Then need to figure out how to implement _sqrt for the various types. There's a few general strategies which could be used:
  • Convert to a drop_sdims(MMatrix(M)) (or maybe just a normal Matrix) and pass to Base to do the hard work (see, eg, the svd implementation). On the way back out, use similar_type to wrap the results in an appropriate static array type, preserving the static array size and type where possible. Or try the option of using Base.invoke to force use of the AbstractMatrix fallback here.
  • Alternatively, copy the implementation from stdlib LinearAlgebra, modifying where necessary to make it work for statically sized arrays. If you want to go really fast for small sizes possibly create an entirely unrolled version (see, eg, the lu() implementation). This has the advantage of allowing optimizations for small sizes, but the disadvantage that we share less code with Base and being a lot more work, especially in tuning. Try this option if you want an adventure!
  • For extra internet points, create specialized optimized versions for tiny matrices (2x2, 3x3).