inv(F::SVD) specialized method

[carstenbauer]

I realized that inv for a SVD currently uses the generic inv(F::Factorization): https://github.com/JuliaLang/julia/blob/25c33e40c5abf5a43d8e45832a889755d6423f38/stdlib/LinearAlgebra/src/factorization.jl#L50

With a simple specialized implementation like

myinv(F::SVD) = (Diagonal(1 ./ F.S) * F.Vt)' * F.U'

I get the following speed up

julia> x = rand(100,100);

julia> F = svd(x);

julia> @btime inv($F);
  187.700 μs (13 allocations: 313.02 KiB)

julia> @btime myinv($F);
  98.100 μs (6 allocations: 157.30 KiB)

julia> isapprox(inv(F), myinv(F))
true

Am I missing something or should this be added to svd.jl?

[StefanKarpinski]

Shouldn't inv(::SVD) give you another SVD object which represents the inverse?

[carstenbauer]

I guess that's a matter of definition. It does currently give you the inverse matrix as for all the other factorizations. In my own packages I define fact_inv for this purpose.

[StefanKarpinski]

Generally, having the SVD of a matrix is much more powerful than having the matrix itself; I would much prefer to keep the inverse factorized as long as it generally works in places where multiplied out inverse matrix would work, which an SVD object generally should.

[carstenbauer]

I would much prefer to keep the inverse factorized as long as it generally works in places where multiplied out inverse matrix would work, which an SVD object generally should

It currently doesn't. For example, neither can I multiply a SVD with something nor add something to it. SVD is not a <: AbstractMatrix and doesn't implement many (anything besides left division?) matrix like operations.

Also, getting the factorization of the inverse of a factorization can sometimes be harder (slower) than getting the inverse matrix. This isn't true for an SVD but should be for a QR if I'm not mistaken (please correct me if I am).

Nonetheless, it makes sense to me to also have a function that calculates the inverse matrix of a factorization efficiently.

[dkarrasch]

In the case of easily invertible factorizations, I agree with @StefanKarpinski that we should keep factorization structure as long as possible. If the user wants a Matrix, then simply construct one. For the F defined in the OP, I get the following.

julia> @btime inv($F);
  86.660 μs (13 allocations: 313.02 KiB)

Now consider this:

myinv2(F::SVD) = SVD(F.V, inv.(F.S), F.U') # perhaps we should check here if !iszero(F.S[end]) or so
julia> @btime myinv2($F);
  98.675 ns (4 allocations: 960 bytes)

Now, constructing the matrix:

myinv3(F::SVD) = Matrix(myinv2(F))
@btime myinv3($F);
  56.079 μs (13 allocations: 235.66 KiB)

And for the original myinv

@btime myinv($F);
  47.879 μs (6 allocations: 157.30 KiB)

If inv returned the matrix right away, one would loose the intermediate decomposition. And, after all, myinv is pretty much nothing but the Matrix constructor for an SVD:

https://github.com/JuliaLang/julia/blob/efd794e199994237e1278a69b38ddad9f949b6b5/stdlib/LinearAlgebra/src/svd.jl#L463-L466

[carstenbauer]

Just to be clear, I do agree that keeping the factorized structure is generally superior.

In the case of easily invertible factorizations, I agree with @StefanKarpinski that we should keep factorization structure as long as possible.

It strikes me as odd to switch between returning a <: Factorization or a matrix depending on the type of factorization given. Especially because factorizations do not behave like matrices.

If the user wants a Matrix, then simply construct one.

As you showed, there is a performance difference, about double the number of allocations and about 17 % slowdown for a SVD (and it would be worse for other factorizations). We could try to argue that, in most performance critical cases, one shouldn't construct the full inverse matrix anyways. But I'm not sure whether this is enough to make the case.

What about being explicit about it and having two distinct functions, maybe inv -> AbstractMatrix and invfact -> Factorization?

[dkarrasch]

Hm, I think I agree, especially since the current behavior of inv(::Factorization) is to return the inverse matric. May I suggest an (arguably) improved version?

julia> myinv(F::SVD) = (F.S .\ F.Vt)' * F.U'
myinv (generic function with 1 method)

julia> @btime myinv($F);
  43.203 μs (4 allocations: 156.41 KiB)

[StefanKarpinski]

In any case, I think we should open a separate issue for the question of whether we should return factorization object in a future Julia version and just make this change as a straightforward performance improvement.

[carstenbauer]

Alright, I'll make a PR. (#32126)