Ben's suggestions.

 * Store rows/cols explicitly for Matrix serialization.
 * Clarifications.

Author: akleeman

albatross/cereal/eigen.h

@@ -29,7 +29,12 @@ template <class Archive, typename _Scalar, int _Rows, int _Cols>
 inline void save(Archive &archive,
                  const Eigen::Matrix<_Scalar, _Rows, _Cols> &v) {
   size_type rows = static_cast<size_type>(v.rows());
+  size_type cols = static_cast<size_type>(v.cols());
   archive(cereal::make_size_tag(rows));
+  // Storing the rows and columns is redundant information but
+  // makes loading a lot easier.
+  archive(rows);
+  archive(cols);
   for (size_type i = 0; i < rows; i++) {
     Eigen::Matrix<_Scalar, _Cols, 1> row = v.row(i);
     archive(row);
@@ -38,28 +43,21 @@ inline void save(Archive &archive,
 
 template <class Archive, typename _Scalar, int _Rows, int _Cols>
 inline void load(Archive &archive, Eigen::Matrix<_Scalar, _Rows, _Cols> &v) {
-  size_type rows;
-  archive(cereal::make_size_tag(rows));
+  size_type rows_plus_two, rows, cols;
+  archive(cereal::make_size_tag(rows_plus_two));
+  archive(rows);
+  archive(cols);
+  assert(rows == rows_plus_two - 2);
   /*
    * In order to determine the size of a matrix, we have to first determine
    * how many rows, then inspect the size of the first row to get the
    * number of columns.
    */
-  if (rows > 0) {
-    Eigen::Matrix<_Scalar, _Rows, 1> first;
-    archive(first);
-    size_type cols = first.rows();
-    v.resize(rows, cols);
-    v.row(0) = first;
-
-    for (size_type i = 1; i < rows; i++) {
-      Eigen::Matrix<_Scalar, _Cols, 1> row;
-      archive(row);
-      v.row(i) = row;
-    }
-  } else {
-    // Serialized matrix is empty.
-    v.resize(0, 0);
+  v.resize(rows, cols);
+  for (size_type i = 0; i < rows; i++) {
+    Eigen::Matrix<_Scalar, _Cols, 1> row;
+    archive(row);
+    v.row(i) = row;
   }
 };
 

albatross/eigen/serializable_ldlt.h

@@ -41,9 +41,15 @@ inline void load_lower_triangle(Archive &archive,
   // TODO: understand why the storage size upon load is always augmented by two.
   storage_size -= 2;
   // We assume the matrix is square and compute the number of rows from the
-  // storage size.
-  double drows = (std::sqrt(static_cast<double>(storage_size * 8 + 1)) - 1) / 2;
-  cereal::size_type rows = static_cast<cereal::size_type>(drows);
+  // storage size using the quadratic formula.
+  //     rows^2 + rows - 2 * storage_size = 0
+  double a = 1;
+  double b = 1;
+  double c = -2. * static_cast<double>(storage_size);
+  double rows_as_double = (std::sqrt(b * b - 4 * a * c) - b) / (2 * a);
+  assert(rows_as_double - static_cast<Eigen::Index>(rows_as_double) == 0. &&
+         "inferred a non integer number of rows");
+  cereal::size_type rows = static_cast<cereal::size_type>(rows_as_double);
   v.resize(rows, rows);
   for (cereal::size_type i = 0; i < rows; i++) {
     for (cereal::size_type j = 0; j <= i; j++) {
@@ -80,16 +86,12 @@ class SerializableLDLT : public LDLT<MatrixXd, Lower> {
         this->transpositionsP() * Eigen::MatrixXd::Identity(n, n);
     this->matrixL().solveInPlace(inverse_cholesky);
 
-    Eigen::VectorXd sqrt_diag = this->vectorD();
-    for (Eigen::Index i = 0; i < n; i++) {
-      sqrt_diag[i] = 1. / std::sqrt(sqrt_diag[i]);
-    }
+    Eigen::VectorXd sqrt_diag = this->vectorD().array().sqrt().matrix();
 
     Eigen::VectorXd inv_diag(n);
     for (Eigen::Index i = 0; i < n; i++) {
       const Eigen::VectorXd col_i =
-          inverse_cholesky.col(i).cwiseProduct(sqrt_diag);
-      ;
+          inverse_cholesky.col(i).cwiseQuotient(sqrt_diag);
       inv_diag[i] = col_i.dot(col_i);
     }