inv(::StaticMatrix) should always return StaticMatrix

[nbaum]

Currently, inv(SMatrix{4, 4}(rand(Float32, 4, 4))) returns an Array{Float32, 2}.

I think it should return a SMatrix{4, 4}.

[andyferris]

Currently only up to 3x3 returns SMatrix with specialised (fast) algorithms. So far things that aren't specialised return Matrix or MMatrix and this is one way users can explore what has been implemented faster than in Base.

Not that this PR is a bad idea at all, just I was hoping to either have a specialised algorithm (maybe just 4x4) or a way of avoiding the extra allocations within Julia to make it a bit speedier.

[andyferris]

By the way, this method was copied from Base.

[andyferris]

Other opinions? @c42f @SimonDanisch

[c42f]

FixedSizeArrays has inverses up to 4x4; looks like they're based on cofactors at a quick glance. I don't know off the top of my head (and couldn't immediately find) a definitive numerical method for small matrices.

As for this PR, I think it's good to make the output type SMatrix because it makes the mutablility of the result more consistent. 👍 from me.

[andyferris]

Let's see if we can get JuliaLang/LinearAlgebra.jl#370 fixed quickly, or else write our own workarounds. Then lufact would work properly and inv would fall back gracefully to MMatrix like eig does (rather than Array like inv does now).

@c42f Do you think preserving immutability is very important? Preserving the static size, yes, but from a performance perspective the problem is largely the allocation and GC which is the same either way (an allocation is done in the working code in Base). With an MMatrix you may avoid unnecessary copies to and around the stack...

[c42f]

Hmm, well type inference should be fine either way. Performance-wise, forcing a copy into an SMatrix will presumably be slightly slower as you say.

To me it's just more consistent at the API level if the output type is always an SMatrix. Currently:

m = blah()
invm = inv(m)
invm[1,1] = 1 # succeeds or fails, depending on implementation details of inv().

I dunno, the size dependent behaviour just seems a bit unfortunate.

[andyferris]

I dunno, the size dependent behaviour just seems a bit unfortunate.

Very true.

[andyferris]

At the moment this returns a SizedArray... is this sufficient @nbaum?

[c42f]

I think the underlying issue here has been solved for some time.