projection onto tangent space to unitary matrices

[Nikdwal]

I don't think that the current implementation of project(::OrthogonalMatrices, p, X) works as intended.
For any X in the tangent space, the projection should act as the identity operator, but it doesn't do that in this example:

using Manifolds
M = OrthogonalMatrices(2)
p = [1 1; -1 1] / √2
@assert is_point(M, p) # test passes
X = [0 1; -1 0]
@assert is_vector(M, p, X) # test passes
X_proj = project(M, p, X)
println(norm(M, p, X - X_proj)) # prints 0.414, not 0 or a small number

If the tangent space is T_p M = { pA | A^T = - A } and it is identified in code with the skew-symmetric matrices through the isomorphism A -> pA, then I think the correct projection should be A -> (A - A^T)/2.

[kellertuer]

The currently implemented – and documented version is

   project(M::OrthogonalMatrices{n}, p, X)
   project(M::Rotations{n}, p, X)
   project(M::UnitaryMatrices{n}, p, X)

  Orthogonally project the tangent vector X ∈ 𝔽^{n × n}, \mathbb F ∈ \{\mathbb R, \mathbb C\} to the tangent space of M at p,
  and change the representer to use the corresponding Lie algebra, i.e. we compute

      \operatorname{proj}_p(X) = \frac{pX-(pX)^{\mathrm{T}}}{2},

So currently, for OrthogonalMatrices we do not do the isomorphism, but sure your variant would (as the easiest way= be

project(SkewHermitianMatrices(2,ℝ), X)

which keeps X as is, yes.

I think we should switch to the isomorphism you mentioned – also because thats nicer for the SO(3) – though we then have to adapt exp. We seem to have changed the default when we unified Unitary/Orthogonal/SpecialOrthogonal/SpecialUnitary about a month ago, sorry for that.

I would like to wait for a short feedback from @mateuszbaran ; but the fix is quite easy; for all 4 it is just the one line at

Manifolds.jl/src/manifolds/GeneralUnitaryMatrices.jl

Line 678 in 8274392

project!(SkewHermitianMatrices(n, 𝔽), Y, p \ X)

and we just have to remove the p \ ;) and currently line 669 states that the change is performed. But you are right, not changing is mit be the more usual approach.

[kellertuer]

I checked our old discussion and the reason for this choice is the following: Since T_p\mathcal M is embedded we added the (not completely correctly documented) p \ X befor skew-symmetrizing.

The PR I opened now proposes to use pX as a representation for the embedded situation for the manifolds and the Lie algebra (skew hermitian) idea (without the ( p \ X) for the Lie groups. But the tests have to still be adapted and I would also document this more precisely.

And while this is not a bug I can see that this asymmetric usage might be irritating. Your example works correct if you first produce the “embedded tangent vector” pX, that is

using Manifolds
M = OrthogonalMatrices(2)
p = [1 1; -1 1] / √2
@assert is_point(M, p) # test passes
X = [0 1; -1 0]
@assert is_vector(M, p, X) # test passes
pX = p*X
X_proj = project(M, p, pX)
println(norm(M, p, X - X_proj)) # is now zero

[mateuszbaran]

I see that projections bite again 🙂 . The current design is that embedding and then projecting is supposed to be identity on manifolds. So your example should be:

julia> X_proj = project(M, p, embed(M, p, X))
2×2 Matrix{Float64}:
  0.0  1.0
 -1.0  0.0

julia> println(norm(M, p, X - X_proj))
0.0

Currently the entire differentiation code is written based on the assumption that project(M, p, embed(M, p, X)) is identity, not just project, so reworking it consistently would be a lot of work. Perhaps the easiest solution would be to add a new function that does the idempotent projection.

[Nikdwal]

Your explanation makes a lot of sense. Thanks.

[kellertuer]

So @mateuszbaran we should keep it as is? I spent a momentanen a PR but it is a little work to change that.

[kellertuer]

Your explanation makes a lot of sense. Thanks.

The one thing that might be irritating here is that both embed and project might change how tangents are represented. So that is why Mateusz explanations better than just claiming project does not change anything on tangent vectors.

[mateuszbaran]

So @mateuszbaran we should keep it as is? I spent a momentanen a PR but it is a little work to change that.

Yes, I think we should keep it as is for now, as changing it would be a lot of work and somewhat breaking. I'd suggest adding a new function that does what is currently handled by embed + project instead of #523 .

[kellertuer]

I accidentally closed this one yesterday then I wanted to close the PR; @Nikdwal thanks for noticing this, do you think we should add some documentation to make this more clear? Maybe even to the generic project/embed functions in ManifoldsBase?

I think the explanation from Mateusz would probably help to be documented somewhere central.