diff --git a/include/core/ceres_wrapper/ceres_manager.h b/include/core/ceres_wrapper/ceres_manager.h
index f124325950a06f74642cff38b6d0c053958cab8e..914002ef1d180aa7e5b85716a84334f5ce944e8a 100644
--- a/include/core/ceres_wrapper/ceres_manager.h
+++ b/include/core/ceres_wrapper/ceres_manager.h
@@ -84,6 +84,8 @@ class CeresManager : public SolverManager
 
         void check();
 
+        const Eigen::SparseMatrixd computeHessian() const;
+
     protected:
 
         std::string solveImpl(const ReportVerbosity report_level) override;
diff --git a/include/core/state_block/state_block.h b/include/core/state_block/state_block.h
index 33aa6fcb4b221ea1dc23f24b21f4a0fb12cbfc0d..ab3bb7147f3c81f57a3436216412274394acc04a 100644
--- a/include/core/state_block/state_block.h
+++ b/include/core/state_block/state_block.h
@@ -74,6 +74,10 @@ public:
          **/
         Eigen::VectorXd getState() const;
 
+        /** \brief Returns the state vector data array
+         **/
+        double* getStateData();
+
         /** \brief Sets the state vector
          **/
         void setState(const Eigen::VectorXd& _state, const bool _notify = true);
@@ -283,6 +287,11 @@ inline void StateBlock::resetLocalParamUpdated()
     local_param_updated_.store(false);
 }
 
+inline double* StateBlock::getStateData()
+{
+    return state_.data();
+}
+
 }// namespace wolf
 
 #endif
diff --git a/src/ceres_wrapper/ceres_manager.cpp b/src/ceres_wrapper/ceres_manager.cpp
index ecf57e7d7e29b4b8ccfecf09f8c39d8eacdafdbb..0d8cba6f97cb7e2c2cdf067e7ebcb65d8ac79c10 100644
--- a/src/ceres_wrapper/ceres_manager.cpp
+++ b/src/ceres_wrapper/ceres_manager.cpp
@@ -443,6 +443,77 @@ void CeresManager::check()
     }
 }
 
+void insertSparseBlock(const Eigen::MatrixXd& ins, Eigen::SparseMatrixd& original, const unsigned int& row, const unsigned int& col)
+{
+    for (auto ins_row = 0; ins_row < ins.rows(); ins_row++)
+        for (auto ins_col = 0; ins_col < ins.cols(); ins_col++)
+            original.insert(row+ins_row, col+ins_col) = ins(ins_row,ins_col);
+
+    original.makeCompressed();
+}
+
+const Eigen::SparseMatrixd CeresManager::computeHessian() const
+{
+    Eigen::SparseMatrixd H;
+    Eigen::SparseMatrixd A;
+    // fac_2_residual_idx_
+    // fac_2_costfunction_
+    // state_blocks_local_param_
+    int ix_row = 0;  // index of the factor/measurement in the 
+    for (auto fac_res_pair: fac_2_residual_idx_){
+        FactorBasePtr fac_ptr = fac_res_pair.first;
+        ix_row += fac_ptr->getSize();
+
+        // retrieve all stateblocks data related to this factor
+        std::vector<const double*> fac_states_ptr;
+        for (auto sb : fac_res_pair.first->getStateBlockPtrVector()){
+            fac_states_ptr.push_back(sb->getStateData());
+        }
+        
+        // retrieve residual value, not used afterwards in this case since we just care about jacobians
+        Eigen::VectorXd residual(fac_ptr->getSize());
+        // retrieve jacobian in group size, not local size
+        std::vector<Eigen::MatrixXd> jacobians;
+        fac_ptr->evaluate(fac_states_ptr, residual, jacobians);
+
+        // Retrieve the block row sparse matrix of this factor
+        Eigen::SparseMatrixd A_block_row(fac_ptr->getSize(), A.cols());
+        for (auto i = 0; i < jacobians.size(); i++)
+        {
+            StateBlockPtr sb = fac_ptr->getStateBlockPtrVector()[i];
+            if (!sb->isFixed())
+            {
+                assert((state_blocks_local_param_.find(sb) != state_blocks_local_param_.end()) && "factor involving a state block not added");
+                assert((A.cols() >= sb->getLocalSize() + jacobians[i].cols()) - 1 && "bad A number of cols");
+
+                // insert since A_block_row has just been created so it's empty for sure
+                if (sb->hasLocalParametrization()){
+                    // if the state block has a local parameterization, we need to right multiply by the manifold element / tangent element jacobian
+                    Eigen::MatrixXd J_manif_tang(sb->getSize(), sb->getLocalSize());
+                    Eigen::Map<Eigen::MatrixRowXd> J_manif_tang_map(J_manif_tang.data(), sb->getSize(), sb->getLocalSize());
+                    Eigen::Map<const Eigen::VectorXd> sb_state_map(sb->getStateData(), sb->getSize());
+
+                    sb->getLocalParametrization()->computeJacobian(sb_state_map, J_manif_tang_map);
+                    insertSparseBlock(jacobians[i] * J_manif_tang, A_block_row, 0, sb->getLocalSize());  // (to_insert, matrix_to_fill, row, col)
+                }
+                else {
+                    insertSparseBlock(jacobians[i], A_block_row, 0, sb->getLocalSize());
+                }
+            }
+        }
+        
+        // fill the weighted jacobian matrix block row
+        A.block(ix_row, 0, fac_ptr->getSize(), A.cols());
+
+    }
+
+    // compute the hessian matrix from the weighted jacobian 
+    H = A.transpose() * A;
+
+    return H;
+}
+
+
 } // namespace wolf
 #include "core/solver/solver_factory.h"
 namespace wolf {