Example for Estimating the Fundamental Matrix

[RoyiAvital]

The Fundamental Matrix, $\boldsymbol{F} \in \mathbb{R}^{3 \times 3}$, is the matrix which given a set of correspondence between 2 images, $\left\{ \left( \boldsymbol{p}_{i}, \boldsymbol{q}_{i} \right) \right\}_{i = 1}^{n}$ , obeys:

$$ \boldsymbol{p}_{i}^{\top} \boldsymbol{F} \boldsymbol{q}_{i} = 0, \; \forall i $$

The Fundamental Matrix must have Rank 2 and by definition is defined up to a scale.

The above can be solved by:

$$ \arg \min_{\boldsymbol{F}} \sum_{i = 1}^{n} {\left( \boldsymbol{p}_{i}^{\top} \boldsymbol{F} \boldsymbol{q}_{i} \right)}^{2}, \; \text{ subject to } \; \boldsymbol{F} \in \mathcal{R}_{2}^{3}, \; {\left\| \boldsymbol{F} \right\|}_{F} = 1 $$

Where $\mathcal{R}_{2}^{3}$ is the set of a $3 \times 3$ matrices with rank 2: $\mathcal{R}_{2}^{3} = \left\{ \boldsymbol{A} \in \mathbb{R}^{3 \times 3} \mid \text{rank} \left( \boldsymbol{A} \right) = 2 \right\}$.

I wonder if such manifold exist?

This is opened after discussion with @kellertuer .

Remark
An extension with solving it with the objective $\arg \min_{\boldsymbol{F}} \sum_{i = 1}^{n} \left| \boldsymbol{p}_{i}^{\top} \boldsymbol{F} \boldsymbol{q}_{i} \right|$ can make it more robust and even a paper worth.

[kellertuer]

This looks interesting. The Frobenius constraint $\lVert F\rVert_{\mathrm{F}}=1$ would yield a “sphere” basically (it is actually a sphere of you vectorise and consider it in $\mathbb R^9$). We have https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere/#Manifolds.ArraySphere for that constrained alone and yes that is a manifold.

The rank constraint $\mathcal R_2^3$ yields this manifold https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/fixedrankmatrices/

My feeling says that should be a (topological) manifold, it might even inherit a metric from the embedding. But from that sentence – that is neither implemented yet (also not in Manopt/Matlab), nor am I aware of literature that would talk about geodesic, exp/log, or retractions. Do you have literature where this manifold is considered?

But the robust case looks even more interesting, since the non smooth case on manifolds is what I usually work on – though then you have to come up with alternatives and can not do “just gradient descent”.

[RoyiAvital]

Indeed the classic algorithm to solve it, called the Direct Linear Transformation (DLT) Algorithm transform the problem into solving $\arg \min_{\boldsymbol{f}} {\left\| \boldsymbol{A} \boldsymbol{f} \right\|}_{2}^{2}, \; \text{ subject to } ; {\left\| \boldsymbol{f} \right\|}_{2} = 1$.

I am not an expert on Manifold Optimization.
Currently I grasp it as a generalization of the Projected Gradient Descent.
Yet instead of only having a correction step of projection one must also take care of the actual gradient step.

Searching for Manifold based solutions yielded sparse results. Not many chose this direction.
Yet there is a similar matrix called the Essential Matrix with similar properties. So it might be a good start.

• An Efficient Solution to Non Minimal Case Essential Matrix Estimation.
• Roberto Tron, Kostas Daniilidis - The Space of Essential Matrices as a Riemannian Quotient Manifold.
• Liran Goshen, Ilan Shimshoni - Recovery of Epipolar Geometry as a Manifold Fitting Problem.
• Gijs Dubbelman, Leo Dorst, Henk Pijls - Manifold Statistics for Essential Matrices.
• Saket Anand, Sushil Mittal, Peter Meer - Robust Estimation for Computer Vision Using Grassmann Manifolds.
• Czech Technical University in Prague - Three Dimensional Computer Vision (TDV) - Sampson Error for Fundamental Matrix Manifold.
• Fang Bai, Adrien Bartoli - The Proxy Step Size Technique for Regularized Optimization on the Sphere Manifold. They solve on the Sphere Manifold and use the Nuclear Norm instead of the Rank constraint. They also have a MATLAB Code.

[kellertuer]

One could of course consider the fixed rank manifold and do a projected gradient thereon – or consider a constrained problem on that manifold. But whether that is faster than existing approaches, I am not sure.

We have a manifold called the essential manifold, but whether that is related, I am not sure.

Thanks for all the links, maybe I find time to take a look some time. My main problem is, that currently issues arrive faster than I find time working on any of them; then I prioritise first bugs and then features – usually both are a bit easier to approach than a complete new example.
One example that I recently started is the relation to Density Functional Theory (e.g. DFTK.jl) – but that for example took me 3 years to get into.

So this might be a good example to add, but we might need someone who has time to work on that; a main challenge for Manopt.jl being a single-person project and Manifolds.jl a 2-person project is indeed very simply time.
I also do not do this full time, but mainly in my free time (unless code is directly related to my research).

[RoyiAvital]

OK, I guess my assumption was that working on the intersection of 2 manifolds is a solved problem just like Projection onto Intersection of Convex Set is solved (Dykstra's Projection Algorithm).
Hence I thought I just need a pointer where is that in the project and will be able to do that my self.

[kellertuer]

Ah, sorry, no. Even the sphere is an intersection of Euclidean and a norm-1 ;) and you can not just take all of Rn and project.
So, Still all formulae change; and before one would have to carefully check even that that actually is a manifold.

[mateuszbaran]

Unit norm, rank 2, 3x3 matrices actually can be expressed using: ProduceManifold(Stiefel(3, 2), Stiefel(3, 2), Sphere(1), through the SVD decomposition:

julia> M1 = Stiefel(3, 2)
Stiefel(3, 2, ℝ)

julia> M2 = Sphere(1)
Sphere(1, ℝ)

julia> V = rand(M1)
3×2 Matrix{Float64}:
 -0.930444   0.168209
 -0.246902  -0.944251
  0.270765  -0.283011

julia> M1 = Stiefel(3, 2)
Stiefel(3, 2, ℝ)

julia> M2 = Sphere(1)
Sphere(1, ℝ)

julia> U, V = rand(M1, 2)
2-element Vector{Matrix{Float64}}:
 [-0.48916066925963686 -0.7241944100053974; -0.35986505354443676 -0.3400692509320744; 0.7944929092741633 -0.5999127945688708]
 [-0.055036094840222205 0.7541384505899712; -0.07497497845898773 -0.6566695815944115; -0.9956654964745 0.007762599460112649]

julia> Σ = rand(M2)
2-element Vector{Float64}:
  0.8576563201823266
 -0.5142233332427754

julia> U * Diagonal(Σ) * V'
3×3 Matrix{Float64}:
 0.303929  -0.213088    0.420604
 0.148864  -0.0916925   0.30866
 0.195142  -0.253664   -0.676054

Where U * Diagonal(Σ) * V' is how you reconstruct the fundamental matrix. Note that multiplication by U and V does not change the norm, so all conditions are satisfied. Now you can freely use any method in Manopt.jl to optimize your objective :). I'd suggest Levenberg-Marquardt for the least squares problem. For the nonsmooth optimization problem I'm not sure what's the best available solver -- maybe the proximal gradient approach, though it'd require some work to figure out the proximal operator.

There are also other robustification approached that might be interesting (like http://ceres-solver.org/nnls_tutorial.html ), though we don't have that in Manopt yet.

If you have a script that computes the objective for some data, I can help you set it up in Manopt.

[kellertuer]

Nice! I did not think of that, a very good observation and parametrisation.