Allocation for TagentSpaceAtPoint fails

[kellertuer]

With

k = 5
N = 20
M = 30
M = FixedRankMatrices(M,N,k)
p = rand(M)
TpM = TangentSpaceAtPoint(M,p)

We can easily generate a random tangent vector with

rand(M; vector_at=p)

which is a UMVTVector as expected, but it seems for rand on the tangent space the allocation of the result is wrong, calling

rand(TpM)

fails because

 allocate_result(TpM,rand)

creates a 30x20 matrix (from the embedding?) instead of an UMVTVector. I could patch this for exactly this one case, but is there maybe a more general way to do this?

[mateuszbaran]

No, we don't currently have a more general way of handling this. I tried designing something better some time ago but it was convoluted and and didn't really help that much.

[kellertuer]

I did not yet even mention my solution ;) What would you propose to solve this?

[mateuszbaran]

function Random.rand(M::AbstractManifold; kwargs...)
    return rand(Random.default_rng(), M; kwargs...)
end

function Random.rand(rng::AbstractRNG, M::TangentSpaceAtPoint; vector_at=nothing)
   return rand(rng, M.fiber.manifold; vector_at=M.point)
end

I think this would be the easiest solution (first method needs to go to ManifoldsBase.jl and the second one to Manifolds.jl)

[kellertuer]

Why the first one is necessary I do not understand.

[mateuszbaran]

It is required because right now we have "first go to allocating, then add default RNG" instead of "first add default RNG, then go to allocating". The second sequence makes it easier to overload non-mutating rand (you don't have to write both rng and non-rng methods).

[kellertuer]

I actually would prefer the first allocation part and was hoping we could just define allocate_result(::TangentSpaceAtPoint, rand)?

[mateuszbaran]

OK, but then you also have to implement allocate_result(::FixedRankMatrices, ::typeof(rand)). Currently FixedRankMatrices implements both mutating and non-mutating rand. So it's doable but a bit more work.

[kellertuer]

I will have to check how to distinguish the tangent vector or point randoms, but I feel fixing allocation is the way that fits better into what we usually did until now.

[mateuszbaran]

Ah, right, this is why rand is like that for FixedRankMatrices. You can't distinguish point and tangent allocation through the existing allocation system. There would be even more work because allocate_result(M, rand) calls would have to actually tell allocate_result what to allocate (point of tangent vector).

[kellertuer]

hm, then I am not sure how to fix this here, because we would either break our usual way to allocate first or it is not possible. If we break it – it is really only the two functions above and not more?

[kellertuer]

Ah besides that if allocate_result would get the kwargs it could actually decide what to allocate. but that would mean also quite some redesign of said function.

[mateuszbaran]

If we break it – it is really only the two functions above and not more?

It fixes the original issue but I haven't run all tests to see if it doesn't break anything.

Ah besides that if allocate_result would get the kwargs it could actually decide what to allocate. but that would mean also quite some redesign of said function.

Maybe not that much, something like this should do:

function Random.rand(M::AbstractManifold; kwargs...)
    if haskey(kwargs, :vector_at)
        pX = allocate_result(M, rand, kwargs[:vector_at])
    else
        pX = allocate_result(M, rand)
    end
    rand!(M, pX; kwargs...)
    return p
end

And then allocate_result(M, rand, p) would mean "allocate tangent vector" and pX = allocate_result(M, rand) would be point allocation.

EDIT: fixed p.

[kellertuer]

That reads for me like a solution that fits our usual approaches. And in principle one could even implement allocate_result(M, rand, p) = zero_vector(M,p) or such. Nice. And this seems to be something we should do in ManifoldsBase then.

[kellertuer]

Ah I do not understand how in your code the magic p appears in line 3 and why that should be the one we return.

[mateuszbaran]

Sorry, I've fixed the p.

[kellertuer]

Ok :) One could also just specify the keyword and pass it on instead of checking for its existence; I feel that is a bit more readable. I will fix this in base then soon.

[kellertuer]

hm, the idea of introducing the zero vector does not work in practice due to an ambiguity with the decorator pattern that I can not resolve (line 91 in DecoratorPattern), but without that (and doing that for every manifold instead) would work.

...or not. I am back to ?rand(TpM) not working then. But I'll open a PR and maybe we can resolve the ambiguity there.