Covariance algorithm for rank deficient Jacobian

[Kyle]

Hi Everyone,

I'm working on bundle adjustement and want to check covariance between estimated parameters.

Ceres supports 3 different algorithms: DENSE_SVD, EIGEN_SPARSE_QR and SUITE_SPARSE_QR.

However, they can't handle my problem with rank deficient and large sparse matrices Jacobian.

E0217 18:34:30.221156 10993 covariance_impl.cc:501] Jacobian matrix is rank deficient. Number of columns: 6253 rank: 6240

F0217 18:34:30.222707 10993 application.cpp:1568] Check failed: covariance.Compute(covariance_blocks_, &problem)

I'd really to know if there are ways to do for large sparse deficient cases.

Thanks for all your help.

Kyle

[Sameer Agarwal]

Kyle,

Two things.

1. Your problem seems to have a larger rank deficiency than the usual gauge freedom 7 - degrees of freedom that bundle adjustment problems have. 

2. Unfortunately, right now there is no sparse rank deficient covariance estimation in Ceres. It requires a sparse SVD algorithm the last time I checked. There is some hope of a complete orthogonal decomposition doing the trick, but I have not looked into it with any depth.

Sameer

--
You received this message because you are subscribed to the Google Groups "Ceres Solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ceres-solver...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/8ea9865d-f692-489f-a498-eb07d9de03aa%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Kyle]

Hi Sameer,

Thanks for your rely.

Could you tell me more details about the first one: gauge freedom?

I try to understand why the Jacobian (in my problem) is rank deficient.

According to Ceres' API reference, there are two types: structural and numerical rank deficiency.

I can ignore the first one because I'm sure all parameter blocks are used and variable.

Moreover, Euler's angles are used to parameterize SO(3).

So, we can ignore the reason arises from overparameterization.

Are there other reasons?

Again, thanks for this great project.

Kyle

[Sameer Agarwal]

Kyle,

A traditional bundle adjustment problem has six degrees of freedom. Because you can scale, rotate and translate the entire reconstruction (including the cameras) without changing the reprojection error. scale + rotation + translation = 1 + 3 + 3 = 7 degrees of freedom. Therefore the Jacobian for bundle adjustment has a null space of rank 7. Another name for this rank deficiency or indeterminacy is Gauge Freedom. This has nothing to do with the fact that all parameters in your problem are used by some residual block.

For more on this I recommend reading "Bundle Adjustment: A Modern Synthesis" or Hartley & Zisserman's book.

Sameer

--

You received this message because you are subscribed to the Google Groups "Ceres Solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ceres-solver...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/f32dd642-6dfe-4af2-9d8a-02d19d644c95%40googlegroups.com.

[Kyle]

 Hi Sameer,

The lastest stable release Eigen 3.2.8 http://eigen.tuxfamily.org/dox/classEigen_1_1SparseQR.html

and SuiteSparse v4.5.1 http://faculty.cse.tamu.edu/davis/suitesparse.html

can handle rank-deficient large sparse matrices.

Could you tell me how I can use these algorithms in Ceres for covariance estimation?

Thanks for your help.

Kyle

[Sameer Agarwal]

Kyle,

As far a I can tell they still cannot be used for actually inverting the matrix the way we want to, but I maybe wrong.

Please file an issue on github and I will take at a look at it as soon as I can.

Sameer

--

You received this message because you are subscribed to the Google Groups "Ceres Solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ceres-solver...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/f1e67e79-624f-46f5-8be1-0a64cbb78cc2%40googlegroups.com.

[Kyle]

Hi Sameer,

I fixed 7 independent parameters in order to obtain minimal parametrization and did 2 experiments with SUITE_SPARSE_QR.

First, used SuiteSparse package 4.2.1 in the official package respository (Debian/jessie) https://packages.debian.org/source/jessie/suitesparse

with the following code in covariance_impl.cc (without change)

647  if (rank < cholmod_jacobian.ncol) {

648    LOG(ERROR) << "Jacobian matrix is rank deficient. "

649               << "Number of columns: " << cholmod_jacobian.ncol

650               << " rank: " << rank;

651     free(permutation);

652    cholmod_l_free_sparse(&R, &cc);

653    cholmod_l_finish(&cc);

654    return false;

655  }

-> Return: Jacobian matrix is rank deficient

Second, used the SuiteSparse 4.5.1  http://faculty.cse.tamu.edu/davis/suitesparse.html

with a small modification in covariance_impl.cc

if (rank < cholmod_jacobian.ncol) {

    LOG(WARNING) << "Jacobian matrix is rank deficient. "

               << "Number of columns: " << cholmod_jacobian.ncol

               << " rank: " << rank;

  }

-> It worked very well, without any warning.

It implied that Jacobian matrix is full deficient.

A+

Kyle

[Kyle]

Second, used the SuiteSparse 4.5.1  http://faculty.cse.tamu.edu/davis/suitesparse.html

with a small modification in covariance_impl.cc

if (rank < cholmod_jacobian.ncol) {

    LOG(WARNING) << "Jacobian matrix is rank deficient. "

               << "Number of columns: " << cholmod_jacobian.ncol

               << " rank: " << rank;

  }

-> It worked very well, without any warning.

It implied that Jacobian matrix is full deficient.

Sorry, I mean that jacobian matrix is full rank. 

[Sameer Agarwal]

Kyle,

Are you implying that SuiteSparseQR is falsely indicating that certain matrices are rank deficient?

--

You received this message because you are subscribed to the Google Groups "Ceres Solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ceres-solver...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/d9c3f325-e5a8-4688-a692-8b6ace74d06c%40googlegroups.com.

[Kyle]

Sameer,

yes, I'm.

I don't know which version of SuiteSparseQR is wrong.

The first returns rank deficient, but the latter indicates full rank.

[Sameer Agarwal]

Maybe worth sharing your jacobian matrix with Prof. Tim Davis.

--

You received this message because you are subscribed to the Google Groups "Ceres Solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ceres-solver...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/139ab6fb-cc06-46d7-bbbf-7e7f8a7a36a4%40googlegroups.com.

[Kyle]

On Tuesday, March 22, 2016 at 1:18:36 PM UTC+1, Kyle wrote:

with a small modification in covariance_impl.cc

if (rank < cholmod_jacobian.ncol) {

    LOG(WARNING) << "Jacobian matrix is rank deficient. "

               << "Number of columns: " << cholmod_jacobian.ncol

               << " rank: " << rank;

  }

 

Apologize for my mistake.

In this case, LOG(WARNING) didn't work well (it didn't return the rank deficient message).

The last question, in ceres, does the function SolveRTRWithSparseRHS(...) in compressed_col_sparse_matrix_utils.h work with a upper trapezoidal R matrix?

Thanks.

[Sameer Agarwal]

I do not think so.

--

You received this message because you are subscribed to the Google Groups "Ceres Solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ceres-solver...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/21f32a52-f1d0-4520-a73e-2660fba9a511%40googlegroups.com.

[Kyle]

Hi Sameer and Cereser,

In bundle adjustment problem, we fix independent parameters to obtain a minimal parametrization.

In theory, rank(J) = n_cols(J).

However, in practical, numerical rank deficiency also occurs unexpectedly (my case).

For a large sparse problem, the good strategy is use SPARSE_SCHUR solver.

Moreover, we only need to have the covariance for camera parameters.

If we use SPARSE_SCHUR, we have a reduced camera matrix S.

Normally, S is full rank and small than H (Hessian matrix).

Importantly, for camera parameters, covariance matrix C = sigma^2 * S^{-1}.

My question: Is there any way to access matrix S in SPARSE_SCHUR solvers?

(I tried to read Ceres's code, but it's too complicated to me.)

Thanks for all your helps.

[Sameer Agarwal]

Kyle,

The short answer is that the schur complement matrix is not accessible to the user.

Also even if it were, it does not get around the rank deficiency problem. Since if the part of the covariance matrix dealing with points is rank deficient, we need to compute the pseudo inverse there, which we do not.

Sameer

--
You received this message because you are subscribed to the Google Groups "Ceres Solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ceres-solver...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/4079edf4-967c-491a-bbd6-bce62152441c%40googlegroups.com.

[Kyle]

Hi Sameer,

In the minimal parametrization case, matrix Hessain is full rank.

If H is of the block matrix form. https://en.wikipedia.org/wiki/Block_matrix_pseudoinverse , see section Block matrix inversion

H = [A B      => H^{-1} = [   S_D^{-1}               -A^{-1} B S_A^{-1}

       C D]                        -D^{-1} C S_D^{-1}        S_A^{-1}          ]

where A is a block sparse matrix corresponding to the cameras, Schur complements of D, S_D = A - B D^{-1} C

We can calculate the camera's covariance matrix by S_D^{-1}.

I checked my program with small sized problem. It's ok, rank(J) = num_col(J).

But with a large sized problem, num_col(J) - rank(J) is in [10,100]. I don't know why.

(I checked the quality of bundle input, the point 3D isn't too far the camera's pose).

That is why I try to find other way  to calculate covariance matrix (only for camera's parameter).

Finally, thank you so much.

[Sameer Agarwal]

Kyle,

S_D maybe rank deficient or D maybe rank deficient. Most likely what you are seeing is illconditioning in your problem manifesting itself as a rank deficiency.

You are better off looking into ways to better condition your problem.

Sameer

--

You received this message because you are subscribed to the Google Groups "Ceres Solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ceres-solver...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/9dc3a184-3c2a-426b-903f-0cc1ebf6ab0b%40googlegroups.com.

[Kyle]

Sameer,

S_D is calculated, like what you did in SPARSE_SCHUR.

When we have S_D, the pseudo inverse of S_D can be computed using Eigen's JacobiSVD.

A large problem becomes a smaller problem, more accurate and reliability.

That is my opinion.

Kyle.

[Sameer Agarwal]

Yes but the schur eliminator does not use SVD right now, so for us to be able to give you a schur matrix correctly that schur eliminator code will need to be modified, besides some significant API changes.
Sameer

--

You received this message because you are subscribed to the Google Groups "Ceres Solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ceres-solver...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/928336f7-2f45-4eb3-a0d7-a68bfd9d37f3%40googlegroups.com.

[Kyle]

I see.

Thanks for your helps.

Kyle