Jacobian matrix - Analytic derivatives

[Kyle]

Hi everybody,

I'm newbie. I'm trying to use ceres with analytic derivatives in order to solve bundle ajusterment problem.

I know someone ask me why I don't use AutoDiffCostFunction but I've some reason to use analytic method.

I read an example (https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/ellipse_approximation.cc) to understand how to use analytic Jacobian.

Context: I have 2 cameras jointly moving, 30 poses and 15k point 3D.

Each observation, i.e. point 2d, corresponds to a triple of a camera, a pose and a point 3D which are identified by id_camera, id_pose, id_point.

15 parameters for each camera, 6 parameters for each pose, and 3 for each point 3D.

Vector<double*> parameters_blocks = <cam0> <cam1> <pose0> ... <pose29> <point3D_0> ... <point3D_14999>

I used linear solver = SPARSE_NORMAL_CHOLESKY.

Could anyone help me how to fill Jacobian matrix?

I tried to do it 

" if (jacobians ==NULL) return true;

// For camera

for (int i = 0; i < 2; ++i) {

  if (jacobians[i] != NULL) {

     ceres::MatrixRef(jacobian[i], 2, 15).setZero();

     if (i = id_camera) {

       // fill jacobian[i][j] = ... ; j = 0 -> 29

     }

  } 

}

// for pose

for (int k = 0; k < 30; ++k) {

  int i = k + 2; // 2 cameras

  if (jacobians[i] != NULL) {

     ceres::MatrixRef(jacobian[i], 2, 6).setZero();

     if (i = id_camera) {

       // fill jacobian[i][j] = ... ; j = 0 -> 11

     }

  } 

}

// for point3D

...

"

but I got an error:

"CHLMOD error: problem too large"

and "Linear solver failed due to unrecoverable non-numeric causes."

Something wrong but I don't know where.

Help me, please!

[Kyle]

Oups, sorry for the mistake

"/ for pose

for (int k = 0; k < 30; ++k) {

  int i = k + 2; // 2 cameras

  if (jacobians[i] != NULL) {

     ceres::MatrixRef(jacobian[i], 2, 6).setZero();

     if (i = id_pose) { // modified here

[William Rucklidge]

Are you making one huge cost function that you're passing >15000 parameter blocks to? It's not clear if that's what you're doing or not... if it is, don't. Break it down into one cost function per observation. This lets Ceres see the block structure of the Jacobian.

-wjr

--
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/7ed202ab-5dbf-47f0-9036-45cc9e8bc37d%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

[Thanh-Tin NGUYEN]

Thank you for your rely.

I'll told you more details

I defined a 

class ReprojectionError : public CostFunction {

public:

   ReprojectionError(const int num_cameras, const int num_poses, const int num_points,

                             const double Obs_x, const double Obs_y,

                             const int camera_index, const int pose_index, const int point_index) : 

                            // set value in private {

     for (int i = 0; i<num_cameras_; ++i)         mutable_parameter_block_sizes()->push_back(15);

     for (int i=0; i < num_poses_; ++i)         mutable_parameter_block_sizes()->push_back(6);

     for (int i=0; i < num_points_; ++i)         mutable_parameter_block_sizes()->push_back(3);

     set_num_residual(2);

  }

  virtual bool Evaluate (double const* const* parameters, // <cam0> <cam1> <pose0> ... <pose29> <point3D_0> ... <point3D_14999>

                                 double* resuduals, double** jacobians) const {

     // Calculate residuals

    residuals[0] = ...;

    residuals[1] = ...;

   

    // Calculate jacobians

    if (jacobians ==NULL) return true;

    // For camera

    for (int i = 0; i < 2; ++i) {

       if (jacobians[i] != NULL) {

          ceres::MatrixRef(jacobian[i], 2, 15).setZero();

          if (i = id_camera) {

              // fill jacobian[i][j] = ... ; j = 0 -> 29

          }

        } 

     }

    // for pose

    for (int k = 0; k < 30; ++k) {

        int i = k + 2; // 2 cameras

        if (jacobians[i] != NULL) {

            ceres::MatrixRef(jacobian[i], 2, 6).setZero();

            if (i = id_camera) {

                // fill jacobian[i][j] = ... ; j = 0 -> 11

            }

        } 

     }

    // for point3D

     ...

  }

  private:

     const int num_cameras_;     const int num_poses_;      const int num_points_;

     const double Obs_x_;      const double Obs_y_;

     const int camera_index_; const int pose_index_; const int point_index_;

};

This is cost function per observation.

@William, do you mean it's a bad idea?

In your opinion, what should I do?

Thanks.

--
You received this message because you are subscribed to a topic in the Google Groups "Ceres Solver" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/ceres-solver/zeJs4Ek4zwA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to ceres-solver...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/CABsh3u5hfHnHJ5PK5YYm-9OR%2Bjv-_G%2B5dA2zA3S7ze9vZ3QOYw%40mail.gmail.com.

[William Rucklidge]

Yes, the way you're doing it is a bad idea. As far as Ceres can tell, each set of 2 residuals is a function of all (looks like) 45210 parameters. Sure, almost all of the Jacobian entries are zero - but Ceres can't know that in advance, and so it has to make a full-size Jacobian (45210x30000) and can't exploit the sparsity structure. And then bad things happen.

What you should be doing instead is having each ReprojectionError push only three entries onto its parameter block list:

- its camera

- its pose

- its point

That way Ceres knows in advance what the block sparsity structure of the Jacobian will be; this creates a 2x15, a 2x6, and a 2x3 block, whereas your current approach creates a 2x45210 block.

-wjr

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/CAL06NGOd454DfnHjRhD-j4M7bwaR1o_4s%2BFQtAF1b8tNNFsr%3Dw%40mail.gmail.com.

[Kyle]

Yeah, my code works. :D

You saved my life.

Thanks for your help.