Make CompressedRowSparseMatrix sensitive to StorageType

Even though we added support for storing the upper and
lower triangular parts of symmetric matrices in
CompressedRowSparseMatrix. RightMultiply, LeftMultiply
and SquaredColumnNorm were not modified to account for this.

This CL changes their implementation and adds thorough
tests.

Also methods that cannot work correctly with symmetric
storage now CHECK and fail, because that indicates
programmer error.

Change-Id: I76288472c8bac98db7376a79bdb6259e346ef2b7

internal/ceres/compressed_row_sparse_matrix.cc

@@ -133,6 +133,26 @@ void AddRandomBlock(const int num_rows,
   }
 }
 
+void AddSymmetricRandomBlock(const int num_rows,
+                             const int row_block_begin,
+                             std::vector<int>* rows,
+                             std::vector<int>* cols,
+                             std::vector<double>* values) {
+  for (int r = 0; r < num_rows; ++r) {
+    for (int c = r; c < num_rows; ++c) {
+      const double v = RandNormal();
+      rows->push_back(row_block_begin + r);
+      cols->push_back(row_block_begin + c);
+      values->push_back(v);
+      if (r != c) {
+        rows->push_back(row_block_begin + c);
+        cols->push_back(row_block_begin + r);
+        values->push_back(v);
+      }
+    }
+  }
+}
+
 }  // namespace
 
 // This constructor gives you a semi-initialized CompressedRowSparseMatrix.
@@ -251,10 +271,60 @@ void CompressedRowSparseMatrix::RightMultiply(const double* x,
   CHECK_NOTNULL(x);
   CHECK_NOTNULL(y);
 
-  for (int r = 0; r < num_rows_; ++r) {
-    for (int idx = rows_[r]; idx < rows_[r + 1]; ++idx) {
-      y[r] += values_[idx] * x[cols_[idx]];
+  if (storage_type_ == UNSYMMETRIC) {
+    for (int r = 0; r < num_rows_; ++r) {
+      for (int idx = rows_[r]; idx < rows_[r + 1]; ++idx) {
+        const int c = cols_[idx];
+        const double v = values_[idx];
+        y[r] += v * x[c];
+      }
     }
+  } else if (storage_type_ == UPPER_TRIANGULAR) {
+    // Because of their block structure, we will have entries that lie
+    // above (below) the diagonal for lower (upper) triangular matrices,
+    // so the loops below need to account for this.
+    for (int r = 0; r < num_rows_; ++r) {
+      int idx = rows_[r];
+      const int idx_end = rows_[r + 1];
+
+      // For upper triangular matrices r <= c, so skip entries with r
+      // > c.
+      while (r > cols_[idx] && idx < idx_end) {
+        ++idx;
+      }
+
+      for (; idx < idx_end; ++idx) {
+        const int c = cols_[idx];
+        const double v = values_[idx];
+        y[r] += v * x[c];
+        // Since we are only iterating over the upper triangular part
+        // of the matrix, add contributions for the strictly lower
+        // triangular part.
+        if (r != c) {
+          y[c] += v * x[r];
+        }
+      }
+    }
+  } else if (storage_type_ == LOWER_TRIANGULAR) {
+    for (int r = 0; r < num_rows_; ++r) {
+      int idx = rows_[r];
+      const int idx_end = rows_[r + 1];
+      // For lower triangular matrices, we only iterate till we are r >=
+      // c.
+      for (; idx < idx_end && r >= cols_[idx]; ++idx) {
+        const int c = cols_[idx];
+        const double v = values_[idx];
+        y[r] += v * x[c];
+        // Since we are only iterating over the lower triangular part
+        // of the matrix, add contributions for the strictly upper
+        // triangular part.
+        if (r != c) {
+          y[c] += v * x[r];
+        }
+      }
+    }
+  } else {
+    LOG(FATAL) << "Unknown storage type: " << storage_type_;
   }
 }
 
@@ -262,10 +332,15 @@ void CompressedRowSparseMatrix::LeftMultiply(const double* x, double* y) const {
   CHECK_NOTNULL(x);
   CHECK_NOTNULL(y);
 
-  for (int r = 0; r < num_rows_; ++r) {
-    for (int idx = rows_[r]; idx < rows_[r + 1]; ++idx) {
-      y[cols_[idx]] += values_[idx] * x[r];
+  if (storage_type_ == UNSYMMETRIC) {
+    for (int r = 0; r < num_rows_; ++r) {
+      for (int idx = rows_[r]; idx < rows_[r + 1]; ++idx) {
+        y[cols_[idx]] += values_[idx] * x[r];
+      }
     }
+  } else {
+    // Since the matrix is symmetric, LeftMultiply = RightMultiply.
+    RightMultiply(x, y);
   }
 }
 
@@ -273,11 +348,58 @@ void CompressedRowSparseMatrix::SquaredColumnNorm(double* x) const {
   CHECK_NOTNULL(x);
 
   std::fill(x, x + num_cols_, 0.0);
-  for (int idx = 0; idx < rows_[num_rows_]; ++idx) {
-    x[cols_[idx]] += values_[idx] * values_[idx];
+  if (storage_type_ == UNSYMMETRIC) {
+    for (int idx = 0; idx < rows_[num_rows_]; ++idx) {
+      x[cols_[idx]] += values_[idx] * values_[idx];
+    }
+  } else if (storage_type_ == UPPER_TRIANGULAR) {
+    // Because of their block structure, we will have entries that lie
+    // above (below) the diagonal for lower (upper) triangular
+    // matrices, so the loops below need to account for this.
+    for (int r = 0; r < num_rows_; ++r) {
+      int idx = rows_[r];
+      const int idx_end = rows_[r + 1];
+
+      // For upper triangular matrices r <= c, so skip entries with r
+      // > c.
+      while (r > cols_[idx] && idx < idx_end) {
+        ++idx;
+      }
+
+      for (; idx < idx_end; ++idx) {
+        const int c = cols_[idx];
+        const double v2 = values_[idx] * values_[idx];
+        x[c] += v2;
+        // Since we are only iterating over the upper triangular part
+        // of the matrix, add contributions for the strictly lower
+        // triangular part.
+        if (r != c) {
+          x[r] += v2;
+        }
+      }
+    }
+  } else if (storage_type_ == LOWER_TRIANGULAR) {
+    for (int r = 0; r < num_rows_; ++r) {
+      int idx = rows_[r];
+      const int idx_end = rows_[r + 1];
+      // For lower triangular matrices, we only iterate till we are r >=
+      // c.
+      for (; idx < idx_end && r >= cols_[idx]; ++idx) {
+        const int c = cols_[idx];
+        const double v2 = values_[idx] * values_[idx];
+        x[c] += v2;
+        // Since we are only iterating over the lower triangular part
+        // of the matrix, add contributions for the strictly upper
+        // triangular part.
+        if (r != c) {
+          x[r] += v2;
+        }
+      }
+    }
+  } else {
+    LOG(FATAL) << "Unknown storage type: " << storage_type_;
   }
 }
-
 void CompressedRowSparseMatrix::ScaleColumns(const double* scale) {
   CHECK_NOTNULL(scale);
 
@@ -301,6 +423,7 @@ void CompressedRowSparseMatrix::ToDenseMatrix(Matrix* dense_matrix) const {
 void CompressedRowSparseMatrix::DeleteRows(int delta_rows) {
   CHECK_GE(delta_rows, 0);
   CHECK_LE(delta_rows, num_rows_);
+  CHECK_EQ(storage_type_, UNSYMMETRIC);
 
   num_rows_ -= delta_rows;
   rows_.resize(num_rows_ + 1);
@@ -324,6 +447,7 @@ void CompressedRowSparseMatrix::DeleteRows(int delta_rows) {
 }
 
 void CompressedRowSparseMatrix::AppendRows(const CompressedRowSparseMatrix& m) {
+  CHECK_EQ(storage_type_, UNSYMMETRIC);
   CHECK_EQ(m.num_cols(), num_cols_);
 
   CHECK((row_blocks_.empty() && m.row_blocks().empty()) ||
@@ -481,15 +605,24 @@ CompressedRowSparseMatrix* CompressedRowSparseMatrix::Transpose() const {
 }
 
 CompressedRowSparseMatrix* CompressedRowSparseMatrix::CreateRandomMatrix(
-    const CompressedRowSparseMatrix::RandomMatrixOptions& options) {
+    CompressedRowSparseMatrix::RandomMatrixOptions options) {
   CHECK_GT(options.num_row_blocks, 0);
   CHECK_GT(options.min_row_block_size, 0);
   CHECK_GT(options.max_row_block_size, 0);
   CHECK_LE(options.min_row_block_size, options.max_row_block_size);
-  CHECK_GT(options.num_col_blocks, 0);
-  CHECK_GT(options.min_col_block_size, 0);
-  CHECK_GT(options.max_col_block_size, 0);
-  CHECK_LE(options.min_col_block_size, options.max_col_block_size);
+
+  if (options.storage_type == UNSYMMETRIC) {
+    CHECK_GT(options.num_col_blocks, 0);
+    CHECK_GT(options.min_col_block_size, 0);
+    CHECK_GT(options.max_col_block_size, 0);
+    CHECK_LE(options.min_col_block_size, options.max_col_block_size);
+  } else {
+    // Symmetric matrices (LOWER_TRIANGULAR or UPPER_TRIANGULAR);
+    options.num_col_blocks = options.num_row_blocks;
+    options.min_col_block_size = options.min_row_block_size;
+    options.max_col_block_size = options.max_row_block_size;
+  }
+
   CHECK_GT(options.block_density, 0.0);
   CHECK_LE(options.block_density, 1.0);
 
@@ -504,12 +637,17 @@ CompressedRowSparseMatrix* CompressedRowSparseMatrix::CreateRandomMatrix(
     row_blocks.push_back(options.min_row_block_size + delta_block_size);
   }
 
-  // Generate the col block structure.
-  for (int i = 0; i < options.num_col_blocks; ++i) {
-    // Generate a random integer in [min_col_block_size, max_col_block_size]
-    const int delta_block_size =
-        Uniform(options.max_col_block_size - options.min_col_block_size);
-    col_blocks.push_back(options.min_col_block_size + delta_block_size);
+  if (options.storage_type == UNSYMMETRIC) {
+    // Generate the col block structure.
+    for (int i = 0; i < options.num_col_blocks; ++i) {
+      // Generate a random integer in [min_col_block_size, max_col_block_size]
+      const int delta_block_size =
+          Uniform(options.max_col_block_size - options.min_col_block_size);
+      col_blocks.push_back(options.min_col_block_size + delta_block_size);
+    }
+  } else {
+    // Symmetric matrices (LOWER_TRIANGULAR or UPPER_TRIANGULAR);
+    col_blocks = row_blocks;
   }
 
   vector<int> tsm_rows;
@@ -533,15 +671,31 @@ CompressedRowSparseMatrix* CompressedRowSparseMatrix::CreateRandomMatrix(
     for (int r = 0; r < options.num_row_blocks; ++r) {
       int col_block_begin = 0;
       for (int c = 0; c < options.num_col_blocks; ++c) {
+        if (((options.storage_type == UPPER_TRIANGULAR) && (r > c)) ||
+            ((options.storage_type == LOWER_TRIANGULAR) && (r < c))) {
+          col_block_begin += col_blocks[c];
+          continue;
+        }
+
         // Randomly determine if this block is present or not.
         if (RandDouble() <= options.block_density) {
-          AddRandomBlock(row_blocks[r],
-                         col_blocks[c],
-                         row_block_begin,
-                         col_block_begin,
-                         &tsm_rows,
-                         &tsm_cols,
-                         &tsm_values);
+          // If the matrix is symmetric, then we take care to generate
+          // symmetric diagonal blocks.
+          if (options.storage_type == UNSYMMETRIC || r != c) {
+            AddRandomBlock(row_blocks[r],
+                           col_blocks[c],
+                           row_block_begin,
+                           col_block_begin,
+                           &tsm_rows,
+                           &tsm_cols,
+                           &tsm_values);
+          } else {
+            AddSymmetricRandomBlock(row_blocks[r],
+                                    row_block_begin,
+                                    &tsm_rows,
+                                    &tsm_cols,
+                                    &tsm_values);
+          }
         }
         col_block_begin += col_blocks[c];
       }
@@ -559,7 +713,7 @@ CompressedRowSparseMatrix* CompressedRowSparseMatrix::CreateRandomMatrix(
           kDoNotTranspose);
   (*matrix->mutable_row_blocks()) = row_blocks;
   (*matrix->mutable_col_blocks()) = col_blocks;
-  matrix->set_storage_type(CompressedRowSparseMatrix::UNSYMMETRIC);
+  matrix->set_storage_type(options.storage_type);
   return matrix;
 }
 

internal/ceres/compressed_row_sparse_matrix.h

@@ -91,9 +91,7 @@ class CompressedRowSparseMatrix : public SparseMatrix {
   // double the peak memory usage.
   //
   // The storage type is set to UNSYMMETRIC.
-  CompressedRowSparseMatrix(int num_rows,
-                            int num_cols,
-                            int max_num_nonzeros);
+  CompressedRowSparseMatrix(int num_rows, int num_cols, int max_num_nonzeros);
 
   // Build a square sparse diagonal matrix with num_rows rows and
   // columns. The diagonal m(i,i) = diagonal(i);
@@ -108,7 +106,6 @@ class CompressedRowSparseMatrix : public SparseMatrix {
   virtual void LeftMultiply(const double* x, double* y) const;
   virtual void SquaredColumnNorm(double* x) const;
   virtual void ScaleColumns(const double* scale);
-
   virtual void ToDenseMatrix(Matrix* dense_matrix) const;
   virtual void ToTextFile(FILE* file) const;
   virtual int num_rows() const { return num_rows_; }
@@ -160,8 +157,7 @@ class CompressedRowSparseMatrix : public SparseMatrix {
   //
   // Caller owns the result.
   static CompressedRowSparseMatrix* CreateBlockDiagonalMatrix(
-      const double* diagonal,
-      const std::vector<int>& blocks);
+      const double* diagonal, const std::vector<int>& blocks);
 
   // Options struct to control the generation of random block sparse
   // matrices in compressed row sparse format.
@@ -178,14 +174,22 @@ class CompressedRowSparseMatrix : public SparseMatrix {
   // block which are distributed normally.
   struct RandomMatrixOptions {
     RandomMatrixOptions()
-        : num_row_blocks(0),
+        : storage_type(UNSYMMETRIC),
           min_row_block_size(0),
           max_row_block_size(0),
           num_col_blocks(0),
           min_col_block_size(0),
           max_col_block_size(0),
-          block_density(0.0) {
-    }
+          block_density(0.0) {}
+
+    // Type of matrix to create.
+    //
+    // If storage_type is UPPER_TRIANGULAR (LOWER_TRIANGULAR), then
+    // create a square symmetric matrix with just the upper triangular
+    // (lower triangular) part. In this case, num_col_blocks,
+    // min_col_block_size and max_col_block_size will be ignored and
+    // assumed to be equal to the corresponding row settings.
+    StorageType storage_type;
 
     int num_row_blocks;
     int min_row_block_size;
@@ -206,7 +210,7 @@ class CompressedRowSparseMatrix : public SparseMatrix {
   //
   // Caller owns the result.
   static CompressedRowSparseMatrix* CreateRandomMatrix(
-      const RandomMatrixOptions& options);
+      RandomMatrixOptions options);
 
  private:
   static CompressedRowSparseMatrix* FromTripletSparseMatrix(

internal/ceres/compressed_row_sparse_matrix_test.cc

@@ -99,35 +99,6 @@ class CompressedRowSparseMatrixTest : public ::testing::Test {
   scoped_ptr<CompressedRowSparseMatrix> crsm;
 };
 
-TEST_F(CompressedRowSparseMatrixTest, RightMultiply) {
-  CompareMatrices(tsm.get(), crsm.get());
-}
-
-TEST_F(CompressedRowSparseMatrixTest, LeftMultiply) {
-  for (int i = 0; i < num_rows; ++i) {
-    Vector a = Vector::Zero(num_rows);
-    a(i) = 1.0;
-
-    Vector b1 = Vector::Zero(num_cols);
-    Vector b2 = Vector::Zero(num_cols);
-
-    tsm->LeftMultiply(a.data(), b1.data());
-    crsm->LeftMultiply(a.data(), b2.data());
-
-    EXPECT_EQ((b1 - b2).norm(), 0);
-  }
-}
-
-TEST_F(CompressedRowSparseMatrixTest, ColumnNorm) {
-  Vector b1 = Vector::Zero(num_cols);
-  Vector b2 = Vector::Zero(num_cols);
-
-  tsm->SquaredColumnNorm(b1.data());
-  crsm->SquaredColumnNorm(b2.data());
-
-  EXPECT_EQ((b1 - b2).norm(), 0);
-}
-
 TEST_F(CompressedRowSparseMatrixTest, Scale) {
   Vector scale(num_cols);
   for (int i = 0; i < num_cols; ++i) {
@@ -188,7 +159,6 @@ TEST_F(CompressedRowSparseMatrixTest, AppendAndDeleteBlockDiagonalMatrix) {
   scoped_ptr<CompressedRowSparseMatrix> appendage(
       CompressedRowSparseMatrix::CreateBlockDiagonalMatrix(
           diagonal.get(), row_and_column_blocks));
-  LOG(INFO) << appendage->row_blocks().size();
 
   crsm->AppendRows(*appendage);
 
@@ -408,6 +378,223 @@ TEST(CompressedRowSparseMatrix, FromTripletSparseMatrixTransposed) {
   }
 }
 
+typedef ::testing::tuple<CompressedRowSparseMatrix::StorageType> Param;
+
+std::string ParamInfoToString(testing::TestParamInfo<Param> info) {
+  if (::testing::get<0>(info.param) ==
+      CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
+    return "UPPER";
+  }
+
+  if (::testing::get<0>(info.param) ==
+      CompressedRowSparseMatrix::LOWER_TRIANGULAR) {
+    return "LOWER";
+  }
+
+  return "UNSYMMETRIC";
+}
+
+class RightMultiplyTest : public ::testing::TestWithParam<Param> {};
+
+TEST_P(RightMultiplyTest, _) {
+  const int kMinNumBlocks = 1;
+  const int kMaxNumBlocks = 10;
+  const int kMinBlockSize = 1;
+  const int kMaxBlockSize = 5;
+  const int kNumTrials = 10;
+
+  for (int num_blocks = kMinNumBlocks; num_blocks < kMaxNumBlocks;
+       ++num_blocks) {
+    for (int trial = 0; trial < kNumTrials; ++trial) {
+      Param param = GetParam();
+      CompressedRowSparseMatrix::RandomMatrixOptions options;
+      options.num_col_blocks = num_blocks;
+      options.min_col_block_size = kMinBlockSize;
+      options.max_col_block_size = kMaxBlockSize;
+      options.num_row_blocks = 2 * num_blocks;
+      options.min_row_block_size = kMinBlockSize;
+      options.max_row_block_size = kMaxBlockSize;
+      options.block_density = std::max(0.5, RandDouble());
+      options.storage_type = ::testing::get<0>(param);
+      scoped_ptr<CompressedRowSparseMatrix> matrix(
+          CompressedRowSparseMatrix::CreateRandomMatrix(options));
+      const int num_rows = matrix->num_rows();
+      const int num_cols = matrix->num_cols();
+
+      Vector x(num_cols);
+      x.setRandom();
+
+      Vector actual_y(num_rows);
+      actual_y.setZero();
+      matrix->RightMultiply(x.data(), actual_y.data());
+
+      Matrix dense;
+      matrix->ToDenseMatrix(&dense);
+      Vector expected_y;
+      if (::testing::get<0>(param) ==
+          CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
+        expected_y = dense.selfadjointView<Eigen::Upper>() * x;
+      } else if (::testing::get<0>(param) ==
+                 CompressedRowSparseMatrix::LOWER_TRIANGULAR) {
+        expected_y = dense.selfadjointView<Eigen::Lower>() * x;
+      } else {
+        expected_y = dense * x;
+      }
+
+      ASSERT_NEAR((expected_y - actual_y).norm() / actual_y.norm(),
+                  0.0,
+                  std::numeric_limits<double>::epsilon() * 10)
+          << "\n"
+          << dense
+          << "x:\n"
+          << x.transpose() << "\n"
+          << "expected: \n" << expected_y.transpose() << "\n"
+          << "actual: \n" << actual_y.transpose();
+    }
+  }
+}
+
+INSTANTIATE_TEST_CASE_P(
+    CompressedRowSparseMatrix,
+    RightMultiplyTest,
+    ::testing::Values(CompressedRowSparseMatrix::LOWER_TRIANGULAR,
+                      CompressedRowSparseMatrix::UPPER_TRIANGULAR,
+                      CompressedRowSparseMatrix::UNSYMMETRIC),
+    ParamInfoToString);
+
+class LeftMultiplyTest : public ::testing::TestWithParam<Param> {};
+
+TEST_P(LeftMultiplyTest, _) {
+  const int kMinNumBlocks = 1;
+  const int kMaxNumBlocks = 10;
+  const int kMinBlockSize = 1;
+  const int kMaxBlockSize = 5;
+  const int kNumTrials = 10;
+
+  for (int num_blocks = kMinNumBlocks; num_blocks < kMaxNumBlocks;
+       ++num_blocks) {
+    for (int trial = 0; trial < kNumTrials; ++trial) {
+      Param param = GetParam();
+      CompressedRowSparseMatrix::RandomMatrixOptions options;
+      options.num_col_blocks = num_blocks;
+      options.min_col_block_size = kMinBlockSize;
+      options.max_col_block_size = kMaxBlockSize;
+      options.num_row_blocks = 2 * num_blocks;
+      options.min_row_block_size = kMinBlockSize;
+      options.max_row_block_size = kMaxBlockSize;
+      options.block_density = std::max(0.5, RandDouble());
+      options.storage_type = ::testing::get<0>(param);
+      scoped_ptr<CompressedRowSparseMatrix> matrix(
+          CompressedRowSparseMatrix::CreateRandomMatrix(options));
+      const int num_rows = matrix->num_rows();
+      const int num_cols = matrix->num_cols();
+
+      Vector x(num_rows);
+      x.setRandom();
+
+      Vector actual_y(num_cols);
+      actual_y.setZero();
+      matrix->LeftMultiply(x.data(), actual_y.data());
+
+      Matrix dense;
+      matrix->ToDenseMatrix(&dense);
+      Vector expected_y;
+      if (::testing::get<0>(param) ==
+          CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
+        expected_y = dense.selfadjointView<Eigen::Upper>() * x;
+      } else if (::testing::get<0>(param) ==
+                 CompressedRowSparseMatrix::LOWER_TRIANGULAR) {
+        expected_y = dense.selfadjointView<Eigen::Lower>() * x;
+      } else {
+        expected_y = dense.transpose() * x;
+      }
+
+      ASSERT_NEAR((expected_y - actual_y).norm() / actual_y.norm(),
+                  0.0,
+                  std::numeric_limits<double>::epsilon() * 10)
+          << "\n"
+          << dense
+          << "x\n"
+          << x.transpose() << "\n"
+          << "expected: \n" << expected_y.transpose() << "\n"
+          << "actual: \n" << actual_y.transpose();
+    }
+  }
+}
+
+INSTANTIATE_TEST_CASE_P(
+    CompressedRowSparseMatrix,
+    LeftMultiplyTest,
+    ::testing::Values(CompressedRowSparseMatrix::LOWER_TRIANGULAR,
+                      CompressedRowSparseMatrix::UPPER_TRIANGULAR,
+                      CompressedRowSparseMatrix::UNSYMMETRIC),
+    ParamInfoToString);
+
+class SquaredColumnNormTest : public ::testing::TestWithParam<Param> {};
+
+TEST_P(SquaredColumnNormTest, _) {
+  const int kMinNumBlocks = 1;
+  const int kMaxNumBlocks = 10;
+  const int kMinBlockSize = 1;
+  const int kMaxBlockSize = 5;
+  const int kNumTrials = 10;
+
+  for (int num_blocks = kMinNumBlocks; num_blocks < kMaxNumBlocks;
+       ++num_blocks) {
+    for (int trial = 0; trial < kNumTrials; ++trial) {
+      Param param = GetParam();
+      CompressedRowSparseMatrix::RandomMatrixOptions options;
+      options.num_col_blocks = num_blocks;
+      options.min_col_block_size = kMinBlockSize;
+      options.max_col_block_size = kMaxBlockSize;
+      options.num_row_blocks = 2 * num_blocks;
+      options.min_row_block_size = kMinBlockSize;
+      options.max_row_block_size = kMaxBlockSize;
+      options.block_density = std::max(0.5, RandDouble());
+      options.storage_type = ::testing::get<0>(param);
+      scoped_ptr<CompressedRowSparseMatrix> matrix(
+          CompressedRowSparseMatrix::CreateRandomMatrix(options));
+      const int num_cols = matrix->num_cols();
+
+      Vector actual(num_cols);
+      actual.setZero();
+      matrix->SquaredColumnNorm(actual.data());
+
+      Matrix dense;
+      matrix->ToDenseMatrix(&dense);
+      Vector expected;
+      if (::testing::get<0>(param) ==
+          CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
+        const Matrix full = dense.selfadjointView<Eigen::Upper>();
+        expected = full.colwise().squaredNorm();
+      } else if (::testing::get<0>(param) ==
+                 CompressedRowSparseMatrix::LOWER_TRIANGULAR) {
+        const Matrix full = dense.selfadjointView<Eigen::Lower>();
+        expected = full.colwise().squaredNorm();
+      } else {
+        expected = dense.colwise().squaredNorm();
+      }
+
+      ASSERT_NEAR((expected - actual).norm() / actual.norm(),
+                  0.0,
+                  std::numeric_limits<double>::epsilon() * 10)
+          << "\n"
+          << dense
+          << "expected: \n" << expected.transpose() << "\n"
+          << "actual: \n" << actual.transpose();
+    }
+  }
+}
+
+INSTANTIATE_TEST_CASE_P(
+    CompressedRowSparseMatrix,
+    SquaredColumnNormTest,
+    ::testing::Values(CompressedRowSparseMatrix::LOWER_TRIANGULAR,
+                      CompressedRowSparseMatrix::UPPER_TRIANGULAR,
+                      CompressedRowSparseMatrix::UNSYMMETRIC),
+    ParamInfoToString);
+
+
 // TODO(sameeragarwal) Add tests for the random matrix creation methods.
 
 }  // namespace internal