diff --git a/internal/ceres/block_jacobi_preconditioner_benchmark.cc b/internal/ceres/block_jacobi_preconditioner_benchmark.cc
index edd431d..b5c1613 100644
--- a/internal/ceres/block_jacobi_preconditioner_benchmark.cc
+++ b/internal/ceres/block_jacobi_preconditioner_benchmark.cc
@@ -85,19 +85,17 @@
   auto jacobian = CreateFakeBundleAdjustmentJacobian(
       kNumCameras, kNumPoints, kCameraSize, kPointSize, kVisibility, prng);
 
-  CompressedRowSparseMatrix jacobian_crs(
-      jacobian->num_rows(), jacobian->num_cols(), jacobian->num_nonzeros());
-  jacobian->ToCompressedRowSparseMatrix(&jacobian_crs);
+  auto jacobian_crs = jacobian->ToCompressedRowSparseMatrix();
   Preconditioner::Options preconditioner_options;
   ContextImpl context;
   preconditioner_options.context = &context;
   preconditioner_options.num_threads = static_cast<int>(state.range(0));
   context.EnsureMinimumThreads(preconditioner_options.num_threads);
-  BlockCRSJacobiPreconditioner p(preconditioner_options, jacobian_crs);
+  BlockCRSJacobiPreconditioner p(preconditioner_options, *jacobian_crs);
 
-  Vector d = Vector::Ones(jacobian_crs.num_cols());
+  Vector d = Vector::Ones(jacobian_crs->num_cols());
   for (auto _ : state) {
-    p.Update(jacobian_crs, d.data());
+    p.Update(*jacobian_crs, d.data());
   }
 }
 
@@ -154,19 +152,17 @@
   std::mt19937 prng;
 
   auto jacobian = BlockSparseMatrix::CreateRandomMatrix(options, prng);
-  CompressedRowSparseMatrix jacobian_crs(
-      jacobian->num_rows(), jacobian->num_cols(), jacobian->num_nonzeros());
-  jacobian->ToCompressedRowSparseMatrix(&jacobian_crs);
+  auto jacobian_crs = jacobian->ToCompressedRowSparseMatrix();
   Preconditioner::Options preconditioner_options;
   ContextImpl context;
   preconditioner_options.context = &context;
   preconditioner_options.num_threads = static_cast<int>(state.range(0));
   context.EnsureMinimumThreads(preconditioner_options.num_threads);
-  BlockCRSJacobiPreconditioner p(preconditioner_options, jacobian_crs);
+  BlockCRSJacobiPreconditioner p(preconditioner_options, *jacobian_crs);
 
-  Vector d = Vector::Ones(jacobian_crs.num_cols());
+  Vector d = Vector::Ones(jacobian_crs->num_cols());
   for (auto _ : state) {
-    p.Update(jacobian_crs, d.data());
+    p.Update(*jacobian_crs, d.data());
   }
 }
 BENCHMARK(BM_BlockCRSJacobiPreconditionerUnstructured)
diff --git a/internal/ceres/block_random_access_sparse_matrix.cc b/internal/ceres/block_random_access_sparse_matrix.cc
index bd9f1ba..992bd01 100644
--- a/internal/ceres/block_random_access_sparse_matrix.cc
+++ b/internal/ceres/block_random_access_sparse_matrix.cc
@@ -53,53 +53,46 @@
   CHECK_LE(blocks.size(), std::numeric_limits<std::int32_t>::max());
 
   const int num_cols = NumScalarEntries(blocks);
+  const int num_blocks = blocks.size();
 
-  // Count the number of scalar non-zero entries and build the layout
-  // object for looking into the values array of the
-  // TripletSparseMatrix.
-  int num_nonzeros = 0;
-  for (const auto& block_pair : block_pairs) {
-    const int row_block_size = blocks_[block_pair.first].size;
-    const int col_block_size = blocks_[block_pair.second].size;
-    num_nonzeros += row_block_size * col_block_size;
+  std::vector<int> num_cells_at_row(num_blocks);
+  for (auto& p : block_pairs) {
+    ++num_cells_at_row[p.first];
   }
-
+  auto block_structure_ = new CompressedRowBlockStructure;
+  block_structure_->cols = blocks;
+  block_structure_->rows.resize(num_blocks);
+  auto p = block_pairs.begin();
+  int num_nonzeros = 0;
+  // Pairs of block indices are sorted lexicographically, thus pairs
+  // corresponding to a single row-block are stored in segments of index pairs
+  // with constant row-block index and increasing column-block index.
+  // CompressedRowBlockStructure is created by traversing block_pairs set.
+  for (int row_block_id = 0; row_block_id < num_blocks; ++row_block_id) {
+    auto& row = block_structure_->rows[row_block_id];
+    row.block = blocks[row_block_id];
+    row.cells.reserve(num_cells_at_row[row_block_id]);
+    const int row_block_size = blocks[row_block_id].size;
+    // Process all index pairs corresponding to the current row block. Because
+    // index pairs are sorted lexicographically, cells are being appended to the
+    // current row-block till the first change in row-block index
+    for (; p != block_pairs.end() && row_block_id == p->first; ++p) {
+      const int col_block_id = p->second;
+      row.cells.emplace_back(col_block_id, num_nonzeros);
+      num_nonzeros += row_block_size * blocks[col_block_id].size;
+    }
+  }
+  bsm_ = std::make_unique<BlockSparseMatrix>(block_structure_);
   VLOG(1) << "Matrix Size [" << num_cols << "," << num_cols << "] "
           << num_nonzeros;
-
-  tsm_ =
-      std::make_unique<TripletSparseMatrix>(num_cols, num_cols, num_nonzeros);
-  tsm_->set_num_nonzeros(num_nonzeros);
-  int* rows = tsm_->mutable_rows();
-  int* cols = tsm_->mutable_cols();
-  double* values = tsm_->mutable_values();
-
-  int pos = 0;
-  for (const auto& block_pair : block_pairs) {
-    const int row_block_size = blocks_[block_pair.first].size;
-    const int col_block_size = blocks_[block_pair.second].size;
-    cell_values_.emplace_back(block_pair, values + pos);
-    layout_[IntPairToInt64(block_pair.first, block_pair.second)] =
-        std::make_unique<CellInfo>(values + pos);
-    pos += row_block_size * col_block_size;
-  }
-
-  // Fill the sparsity pattern of the underlying matrix.
-  for (const auto& block_pair : block_pairs) {
-    const int row_block_id = block_pair.first;
-    const int col_block_id = block_pair.second;
-    const int row_block_size = blocks_[row_block_id].size;
-    const int col_block_size = blocks_[col_block_id].size;
-    int pos =
-        layout_[IntPairToInt64(row_block_id, col_block_id)]->values - values;
-    for (int r = 0; r < row_block_size; ++r) {
-      for (int c = 0; c < col_block_size; ++c, ++pos) {
-        rows[pos] = blocks_[row_block_id].position + r;
-        cols[pos] = blocks_[col_block_id].position + c;
-        values[pos] = 1.0;
-        DCHECK_LT(rows[pos], tsm_->num_rows());
-        DCHECK_LT(cols[pos], tsm_->num_rows());
-      }
+  double* values = bsm_->mutable_values();
+  for (int row_block_id = 0; row_block_id < num_blocks; ++row_block_id) {
+    const auto& cells = block_structure_->rows[row_block_id].cells;
+    for (auto& c : cells) {
+      const int col_block_id = c.block_id;
+      double* const data = values + c.position;
+      layout_[IntPairToInt64(row_block_id, col_block_id)] =
+          std::make_unique<CellInfo>(data);
     }
   }
 }
@@ -126,35 +119,41 @@
 // Assume that the user does not hold any locks on any cell blocks
 // when they are calling SetZero.
 void BlockRandomAccessSparseMatrix::SetZero() {
-  ParallelSetZero(
-      context_, num_threads_, tsm_->mutable_values(), tsm_->num_nonzeros());
+  bsm_->SetZero(context_, num_threads_);
 }
 
 void BlockRandomAccessSparseMatrix::SymmetricRightMultiplyAndAccumulate(
     const double* x, double* y) const {
-  for (const auto& cell_position_and_data : cell_values_) {
-    const int row = cell_position_and_data.first.first;
-    const int row_block_size = blocks_[row].size;
-    const int row_block_pos = blocks_[row].position;
+  const auto bs = bsm_->block_structure();
+  const auto values = bsm_->values();
+  const int num_blocks = blocks_.size();
 
-    const int col = cell_position_and_data.first.second;
-    const int col_block_size = blocks_[col].size;
-    const int col_block_pos = blocks_[col].position;
+  for (int row_block_id = 0; row_block_id < num_blocks; ++row_block_id) {
+    const auto& row_block = bs->rows[row_block_id];
+    const int row_block_size = row_block.block.size;
+    const int row_block_pos = row_block.block.position;
 
-    MatrixVectorMultiply<Eigen::Dynamic, Eigen::Dynamic, 1>(
-        cell_position_and_data.second,
-        row_block_size,
-        col_block_size,
-        x + col_block_pos,
-        y + row_block_pos);
+    for (auto& c : row_block.cells) {
+      const int col_block_id = c.block_id;
+      const int col_block_size = blocks_[col_block_id].size;
+      const int col_block_pos = blocks_[col_block_id].position;
 
-    // Since the matrix is symmetric, but only the upper triangular
-    // part is stored, if the block being accessed is not a diagonal
-    // block, then use the same block to do the corresponding lower
-    // triangular multiply also.
-    if (row != col) {
+      MatrixVectorMultiply<Eigen::Dynamic, Eigen::Dynamic, 1>(
+          values + c.position,
+          row_block_size,
+          col_block_size,
+          x + col_block_pos,
+          y + row_block_pos);
+      if (col_block_id == row_block_id) {
+        continue;
+      }
+
+      // Since the matrix is symmetric, but only the upper triangular
+      // part is stored, if the block being accessed is not a diagonal
+      // block, then use the same block to do the corresponding lower
+      // triangular multiply also
       MatrixTransposeVectorMultiply<Eigen::Dynamic, Eigen::Dynamic, 1>(
-          cell_position_and_data.second,
+          values + c.position,
           row_block_size,
           col_block_size,
           x + row_block_pos,
diff --git a/internal/ceres/block_random_access_sparse_matrix.h b/internal/ceres/block_random_access_sparse_matrix.h
index a64a139..cb6f511 100644
--- a/internal/ceres/block_random_access_sparse_matrix.h
+++ b/internal/ceres/block_random_access_sparse_matrix.h
@@ -39,18 +39,18 @@
 #include <vector>
 
 #include "ceres/block_random_access_matrix.h"
+#include "ceres/block_sparse_matrix.h"
 #include "ceres/block_structure.h"
 #include "ceres/context_impl.h"
 #include "ceres/internal/disable_warnings.h"
 #include "ceres/internal/export.h"
 #include "ceres/small_blas.h"
-#include "ceres/triplet_sparse_matrix.h"
 #include "ceres/types.h"
 
 namespace ceres::internal {
 
 // A thread safe square block sparse implementation of
-// BlockRandomAccessMatrix. Internally a TripletSparseMatrix is used
+// BlockRandomAccessMatrix. Internally a BlockSparseMatrix is used
 // for doing the actual storage. This class augments this matrix with
 // an unordered_map that allows random read/write access.
 class CERES_NO_EXPORT BlockRandomAccessSparseMatrix
@@ -81,19 +81,19 @@
   // locked.
   void SetZero() final;
 
-  // Assume that the matrix is symmetric and only one half of the
-  // matrix is stored.
+  // Assume that the matrix is symmetric and only one half of the matrix is
+  // stored.
   //
   // y += S * x
   void SymmetricRightMultiplyAndAccumulate(const double* x, double* y) const;
 
   // Since the matrix is square, num_rows() == num_cols().
-  int num_rows() const final { return tsm_->num_rows(); }
-  int num_cols() const final { return tsm_->num_cols(); }
+  int num_rows() const final { return bsm_->num_rows(); }
+  int num_cols() const final { return bsm_->num_cols(); }
 
   // Access to the underlying matrix object.
-  const TripletSparseMatrix* matrix() const { return tsm_.get(); }
-  TripletSparseMatrix* mutable_matrix() { return tsm_.get(); }
+  const BlockSparseMatrix* matrix() const { return bsm_.get(); }
+  BlockSparseMatrix* mutable_matrix() { return bsm_.get(); }
 
  private:
   int64_t IntPairToInt64(int row, int col) const {
@@ -117,12 +117,8 @@
   using LayoutType = std::unordered_map<int64_t, std::unique_ptr<CellInfo>>;
   LayoutType layout_;
 
-  // In order traversal of contents of the matrix. This allows us to
-  // implement a matrix-vector product which is 20% faster than using the
-  // iterator in the Layout object instead.
-  std::vector<std::pair<std::pair<int, int>, double*>> cell_values_;
   // The underlying matrix object which actually stores the cells.
-  std::unique_ptr<TripletSparseMatrix> tsm_;
+  std::unique_ptr<BlockSparseMatrix> bsm_;
 
   friend class BlockRandomAccessSparseMatrixTest;
 };
diff --git a/internal/ceres/block_random_access_sparse_matrix_test.cc b/internal/ceres/block_random_access_sparse_matrix_test.cc
index 2e7d439..f8afc1e 100644
--- a/internal/ceres/block_random_access_sparse_matrix_test.cc
+++ b/internal/ceres/block_random_access_sparse_matrix_test.cc
@@ -91,12 +91,11 @@
         Matrix::Ones(blocks[row_block_id].size, blocks[col_block_id].size);
   }
 
-  const TripletSparseMatrix* tsm = m.matrix();
-  EXPECT_EQ(tsm->num_nonzeros(), num_nonzeros);
-  EXPECT_EQ(tsm->max_num_nonzeros(), num_nonzeros);
+  const BlockSparseMatrix* bsm = m.matrix();
+  EXPECT_EQ(bsm->num_nonzeros(), num_nonzeros);
 
   Matrix dense;
-  tsm->ToDenseMatrix(&dense);
+  bsm->ToDenseMatrix(&dense);
 
   double kTolerance = 1e-14;
 
diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc
index a60ee89..df2d9e1 100644
--- a/internal/ceres/block_sparse_matrix.cc
+++ b/internal/ceres/block_sparse_matrix.cc
@@ -58,6 +58,120 @@
     rows[c].cumulative_nnz = curr_nnz + rows[c - 1].cumulative_nnz;
   }
 }
+
+template <bool transpose>
+std::unique_ptr<CompressedRowSparseMatrix>
+CreateStructureOfCompressedRowSparseMatrix(
+    const double* values,
+    int num_rows,
+    int num_cols,
+    int num_nonzeros,
+    const CompressedRowBlockStructure* block_structure) {
+  auto crs_matrix = std::make_unique<CompressedRowSparseMatrix>(
+      num_rows, num_cols, num_nonzeros);
+  auto crs_cols = crs_matrix->mutable_cols();
+  auto crs_rows = crs_matrix->mutable_rows();
+  int value_offset = 0;
+  const int num_row_blocks = block_structure->rows.size();
+  const auto& cols = block_structure->cols;
+  *crs_rows++ = 0;
+  for (int row_block_id = 0; row_block_id < num_row_blocks; ++row_block_id) {
+    const auto& row_block = block_structure->rows[row_block_id];
+    // Empty row block: only requires setting row offsets
+    if (row_block.cells.empty()) {
+      std::fill(crs_rows, crs_rows + row_block.block.size, value_offset);
+      crs_rows += row_block.block.size;
+      continue;
+    }
+
+    // Compute number of non-zeros per row from number of non-zeros in row-block
+    int row_block_nnz;
+    if constexpr (transpose) {
+      // Transposed block structure comes with nnz filled-in
+      row_block_nnz = row_block.nnz;
+    } else {
+      // Non-transposed block structure has sequential layout, but nnz field of
+      // block structure is not filled. Difference between values position of
+      // the first and the last cells in row-block gives number of non-zero
+      // elements in row-block excluding the last cell.
+      const auto& front = row_block.cells.front();
+      const auto& back = row_block.cells.back();
+      row_block_nnz =
+          (back.position - front.position) +
+          block_structure->cols[back.block_id].size * row_block.block.size;
+    }
+    const int row_nnz = row_block_nnz / row_block.block.size;
+
+    // Row-wise setup of matrix structure
+    for (int row = 0; row < row_block.block.size; ++row) {
+      value_offset += row_nnz;
+      *crs_rows++ = value_offset;
+      for (auto& c : row_block.cells) {
+        const int col_block_size = cols[c.block_id].size;
+        const int col_position = cols[c.block_id].position;
+        std::iota(crs_cols, crs_cols + col_block_size, col_position);
+        crs_cols += col_block_size;
+      }
+    }
+  }
+  return crs_matrix;
+}
+
+template <bool transpose>
+void UpdateCompressedRowSparseMatrixImpl(
+    CompressedRowSparseMatrix* crs_matrix,
+    const double* values,
+    const CompressedRowBlockStructure* block_structure) {
+  auto crs_values = crs_matrix->mutable_values();
+  auto crs_rows = crs_matrix->mutable_rows();
+  const int num_row_blocks = block_structure->rows.size();
+  const auto& cols = block_structure->cols;
+  for (int row_block_id = 0; row_block_id < num_row_blocks; ++row_block_id) {
+    const auto& row_block = block_structure->rows[row_block_id];
+    const int row_block_size = row_block.block.size;
+    const int row_nnz = crs_rows[1] - crs_rows[0];
+    crs_rows += row_block_size;
+
+    if (row_nnz == 0) {
+      continue;
+    }
+
+    MatrixRef crs_row_block(crs_values, row_block_size, row_nnz);
+    int col_offset = 0;
+    for (auto& c : row_block.cells) {
+      const int col_block_size = cols[c.block_id].size;
+      auto crs_cell =
+          crs_row_block.block(0, col_offset, row_block_size, col_block_size);
+      if constexpr (transpose) {
+        // Transposed matrix is filled using transposed block-strucutre
+        ConstMatrixRef cell(
+            values + c.position, col_block_size, row_block_size);
+        crs_cell = cell.transpose();
+      } else {
+        ConstMatrixRef cell(
+            values + c.position, row_block_size, col_block_size);
+        crs_cell = cell;
+      }
+      col_offset += col_block_size;
+    }
+    crs_values += row_nnz * row_block_size;
+  }
+}
+
+void SetBlockStructureOfCompressedRowSparseMatrix(
+    CompressedRowSparseMatrix* crs_matrix,
+    CompressedRowBlockStructure* block_structure) {
+  const int num_row_blocks = block_structure->rows.size();
+  auto& row_blocks = *crs_matrix->mutable_row_blocks();
+  row_blocks.resize(num_row_blocks);
+  for (int i = 0; i < num_row_blocks; ++i) {
+    row_blocks[i] = block_structure->rows[i].block;
+  }
+
+  auto& col_blocks = *crs_matrix->mutable_col_blocks();
+  col_blocks = block_structure->cols;
+}
+
 }  // namespace
 
 BlockSparseMatrix::BlockSparseMatrix(
@@ -330,29 +444,47 @@
       transpose_bs->rows.data(),
       [](const CompressedRow& row) { return row.cumulative_nnz; });
 }
+std::unique_ptr<CompressedRowSparseMatrix>
+BlockSparseMatrix::ToCompressedRowSparseMatrixTranspose() const {
+  auto bs = transpose_block_structure_.get();
+  auto crs_matrix = CreateStructureOfCompressedRowSparseMatrix<true>(
+      values(), num_cols_, num_rows_, num_nonzeros_, bs);
 
-void BlockSparseMatrix::ToCompressedRowSparseMatrix(
+  SetBlockStructureOfCompressedRowSparseMatrix(crs_matrix.get(), bs);
+
+  UpdateCompressedRowSparseMatrixTranspose(crs_matrix.get());
+  return crs_matrix;
+}
+
+std::unique_ptr<CompressedRowSparseMatrix>
+BlockSparseMatrix::ToCompressedRowSparseMatrix() const {
+  auto crs_matrix = CreateStructureOfCompressedRowSparseMatrix<false>(
+      values(), num_rows_, num_cols_, num_nonzeros_, block_structure_.get());
+
+  SetBlockStructureOfCompressedRowSparseMatrix(crs_matrix.get(),
+                                               block_structure_.get());
+
+  UpdateCompressedRowSparseMatrix(crs_matrix.get());
+  return crs_matrix;
+}
+
+void BlockSparseMatrix::UpdateCompressedRowSparseMatrixTranspose(
     CompressedRowSparseMatrix* crs_matrix) const {
-  {
-    TripletSparseMatrix ts_matrix;
-    this->ToTripletSparseMatrix(&ts_matrix);
-    *crs_matrix =
-        *CompressedRowSparseMatrix::FromTripletSparseMatrix(ts_matrix);
-  }
-
-  int num_row_blocks = block_structure_->rows.size();
-  auto& row_blocks = *crs_matrix->mutable_row_blocks();
-  row_blocks.resize(num_row_blocks);
-  for (int i = 0; i < num_row_blocks; ++i) {
-    row_blocks[i] = block_structure_->rows[i].block;
-  }
-
-  int num_col_blocks = block_structure_->cols.size();
-  auto& col_blocks = *crs_matrix->mutable_col_blocks();
-  col_blocks.resize(num_col_blocks);
-  for (int i = 0; i < num_col_blocks; ++i) {
-    col_blocks[i] = block_structure_->cols[i];
-  }
+  CHECK(crs_matrix != nullptr);
+  CHECK_EQ(crs_matrix->num_rows(), num_cols_);
+  CHECK_EQ(crs_matrix->num_cols(), num_rows_);
+  CHECK_EQ(crs_matrix->num_nonzeros(), num_nonzeros_);
+  UpdateCompressedRowSparseMatrixImpl<true>(
+      crs_matrix, values(), transpose_block_structure_.get());
+}
+void BlockSparseMatrix::UpdateCompressedRowSparseMatrix(
+    CompressedRowSparseMatrix* crs_matrix) const {
+  CHECK(crs_matrix != nullptr);
+  CHECK_EQ(crs_matrix->num_rows(), num_rows_);
+  CHECK_EQ(crs_matrix->num_cols(), num_cols_);
+  CHECK_EQ(crs_matrix->num_nonzeros(), num_nonzeros_);
+  UpdateCompressedRowSparseMatrixImpl<false>(
+      crs_matrix, values(), block_structure_.get());
 }
 
 void BlockSparseMatrix::ToDenseMatrix(Matrix* dense_matrix) const {
diff --git a/internal/ceres/block_sparse_matrix.h b/internal/ceres/block_sparse_matrix.h
index ed3146f..55f1cc4 100644
--- a/internal/ceres/block_sparse_matrix.h
+++ b/internal/ceres/block_sparse_matrix.h
@@ -91,7 +91,20 @@
   void ScaleColumns(const double* scale,
                     ContextImpl* context,
                     int num_threads) final;
-  void ToCompressedRowSparseMatrix(CompressedRowSparseMatrix* matrix) const;
+
+  // Convert to CompressedRowSparseMatrix
+  std::unique_ptr<CompressedRowSparseMatrix> ToCompressedRowSparseMatrix()
+      const;
+  // Create CompressedRowSparseMatrix corresponding to transposed matrix
+  std::unique_ptr<CompressedRowSparseMatrix>
+  ToCompressedRowSparseMatrixTranspose() const;
+  // Copy values to CompressedRowSparseMatrix that has compatible structure
+  void UpdateCompressedRowSparseMatrix(
+      CompressedRowSparseMatrix* crs_matrix) const;
+  // Copy values to CompressedRowSparseMatrix that has structure of transposed
+  // matrix
+  void UpdateCompressedRowSparseMatrixTranspose(
+      CompressedRowSparseMatrix* crs_matrix) const;
   void ToDenseMatrix(Matrix* dense_matrix) const final;
   void ToTextFile(FILE* file) const final;
 
diff --git a/internal/ceres/block_sparse_matrix_test.cc b/internal/ceres/block_sparse_matrix_test.cc
index 5036aa2..f749df3 100644
--- a/internal/ceres/block_sparse_matrix_test.cc
+++ b/internal/ceres/block_sparse_matrix_test.cc
@@ -482,38 +482,73 @@
 TEST(BlockSparseMatrix, ToCRSMatrix) {
   {
     std::unique_ptr<BlockSparseMatrix> m = CreateTestMatrixFromId(0);
-    CompressedRowSparseMatrix m_crs(
-        m->num_rows(), m->num_cols(), m->num_nonzeros());
-    m->ToCompressedRowSparseMatrix(&m_crs);
+    auto m_crs = m->ToCompressedRowSparseMatrix();
     std::vector<int> rows_expected = {0, 2, 4, 7, 10};
     std::vector<int> cols_expected = {0, 1, 0, 1, 2, 3, 4, 2, 3, 4};
     std::vector<double> values_expected = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
     for (int i = 0; i < rows_expected.size(); ++i) {
-      EXPECT_EQ(m_crs.rows()[i], rows_expected[i]);
+      EXPECT_EQ(m_crs->rows()[i], rows_expected[i]);
     }
     for (int i = 0; i < cols_expected.size(); ++i) {
-      EXPECT_EQ(m_crs.cols()[i], cols_expected[i]);
+      EXPECT_EQ(m_crs->cols()[i], cols_expected[i]);
     }
     for (int i = 0; i < values_expected.size(); ++i) {
-      EXPECT_EQ(m_crs.values()[i], values_expected[i]);
+      EXPECT_EQ(m_crs->values()[i], values_expected[i]);
     }
   }
   {
     std::unique_ptr<BlockSparseMatrix> m = CreateTestMatrixFromId(1);
-    CompressedRowSparseMatrix m_crs(
-        m->num_rows(), m->num_cols(), m->num_nonzeros());
-    m->ToCompressedRowSparseMatrix(&m_crs);
+    auto m_crs = m->ToCompressedRowSparseMatrix();
     std::vector<int> rows_expected = {0, 4, 8, 9};
     std::vector<int> cols_expected = {0, 1, 3, 4, 0, 1, 3, 4, 2};
     std::vector<double> values_expected = {1, 2, 5, 6, 3, 4, 7, 8, 9};
     for (int i = 0; i < rows_expected.size(); ++i) {
-      EXPECT_EQ(m_crs.rows()[i], rows_expected[i]);
+      EXPECT_EQ(m_crs->rows()[i], rows_expected[i]);
     }
     for (int i = 0; i < cols_expected.size(); ++i) {
-      EXPECT_EQ(m_crs.cols()[i], cols_expected[i]);
+      EXPECT_EQ(m_crs->cols()[i], cols_expected[i]);
     }
     for (int i = 0; i < values_expected.size(); ++i) {
-      EXPECT_EQ(m_crs.values()[i], values_expected[i]);
+      EXPECT_EQ(m_crs->values()[i], values_expected[i]);
+    }
+  }
+}
+
+TEST(BlockSparseMatrix, ToCRSMatrixTranspose) {
+  {
+    std::unique_ptr<BlockSparseMatrix> m = CreateTestMatrixFromId(0);
+    auto m_crs_transpose = m->ToCompressedRowSparseMatrixTranspose();
+    std::vector<int> rows_expected = {0, 2, 4, 6, 8, 10, 10};
+    std::vector<int> cols_expected = {0, 1, 0, 1, 2, 3, 2, 3, 2, 3};
+    std::vector<double> values_expected = {1, 3, 2, 4, 5, 8, 6, 9, 7, 10};
+    EXPECT_EQ(m_crs_transpose->num_nonzeros(), cols_expected.size());
+    EXPECT_EQ(m_crs_transpose->num_rows(), rows_expected.size() - 1);
+    for (int i = 0; i < rows_expected.size(); ++i) {
+      EXPECT_EQ(m_crs_transpose->rows()[i], rows_expected[i]);
+    }
+    for (int i = 0; i < cols_expected.size(); ++i) {
+      EXPECT_EQ(m_crs_transpose->cols()[i], cols_expected[i]);
+    }
+    for (int i = 0; i < values_expected.size(); ++i) {
+      EXPECT_EQ(m_crs_transpose->values()[i], values_expected[i]);
+    }
+  }
+  {
+    std::unique_ptr<BlockSparseMatrix> m = CreateTestMatrixFromId(1);
+    auto m_crs_transpose = m->ToCompressedRowSparseMatrixTranspose();
+    std::vector<int> rows_expected = {0, 2, 4, 5, 7, 9, 9};
+    std::vector<int> cols_expected = {0, 1, 0, 1, 2, 0, 1, 0, 1};
+    std::vector<double> values_expected = {1, 3, 2, 4, 9, 5, 7, 6, 8};
+    EXPECT_EQ(m_crs_transpose->num_nonzeros(), cols_expected.size());
+    EXPECT_EQ(m_crs_transpose->num_rows(), rows_expected.size() - 1);
+    for (int i = 0; i < rows_expected.size(); ++i) {
+      EXPECT_EQ(m_crs_transpose->rows()[i], rows_expected[i]);
+    }
+    for (int i = 0; i < cols_expected.size(); ++i) {
+      EXPECT_EQ(m_crs_transpose->cols()[i], cols_expected[i]);
+    }
+    for (int i = 0; i < values_expected.size(); ++i) {
+      EXPECT_EQ(m_crs_transpose->values()[i], values_expected[i]);
     }
   }
 }
diff --git a/internal/ceres/cuda_sparse_matrix_test.cc b/internal/ceres/cuda_sparse_matrix_test.cc
index a12ad40..74e74fd 100644
--- a/internal/ceres/cuda_sparse_matrix_test.cc
+++ b/internal/ceres/cuda_sparse_matrix_test.cc
@@ -80,10 +80,8 @@
 
 TEST_F(CudaSparseMatrixTest, RightMultiplyAndAccumulate) {
   std::string message;
-  CompressedRowSparseMatrix A_crs(
-      A_->num_rows(), A_->num_cols(), A_->num_nonzeros());
-  A_->ToCompressedRowSparseMatrix(&A_crs);
-  CudaSparseMatrix A_gpu(&context_, A_crs);
+  auto A_crs = A_->ToCompressedRowSparseMatrix();
+  CudaSparseMatrix A_gpu(&context_, *A_crs);
   CudaVector x_gpu(&context_, A_gpu.num_cols());
   CudaVector res_gpu(&context_, A_gpu.num_rows());
   x_gpu.CopyFromCpu(x_);
diff --git a/internal/ceres/evaluation_benchmark.cc b/internal/ceres/evaluation_benchmark.cc
index 5dc1ef4..ae1a65d 100644
--- a/internal/ceres/evaluation_benchmark.cc
+++ b/internal/ceres/evaluation_benchmark.cc
@@ -128,13 +128,7 @@
   std::unique_ptr<CompressedRowSparseMatrix> CreateCompressedRowSparseJacobian(
       ContextImpl* context) {
     auto block_sparse = BlockSparseJacobian(context);
-    auto crs_jacobian = std::make_unique<CompressedRowSparseMatrix>(
-        block_sparse->num_rows(),
-        block_sparse->num_cols(),
-        block_sparse->num_nonzeros());
-
-    block_sparse->ToCompressedRowSparseMatrix(crs_jacobian.get());
-    return crs_jacobian;
+    return block_sparse->ToCompressedRowSparseMatrix();
   }
 
   const BlockSparseMatrix* BlockSparseJacobian(ContextImpl* context) {
diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc
index aedc73a..e49d5bb 100644
--- a/internal/ceres/schur_complement_solver.cc
+++ b/internal/ceres/schur_complement_solver.cc
@@ -299,32 +299,36 @@
   summary.termination_type = LinearSolverTerminationType::SUCCESS;
   summary.message = "Success.";
 
-  const TripletSparseMatrix* tsm =
+  const BlockSparseMatrix* bsm =
       down_cast<const BlockRandomAccessSparseMatrix*>(lhs())->matrix();
-  if (tsm->num_rows() == 0) {
+  if (bsm->num_rows() == 0) {
     return summary;
   }
 
-  std::unique_ptr<CompressedRowSparseMatrix> lhs;
   const CompressedRowSparseMatrix::StorageType storage_type =
       sparse_cholesky_->StorageType();
   if (storage_type ==
       CompressedRowSparseMatrix::StorageType::UPPER_TRIANGULAR) {
-    lhs = CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm);
-    lhs->set_storage_type(
-        CompressedRowSparseMatrix::StorageType::UPPER_TRIANGULAR);
+    if (!crs_lhs_) {
+      crs_lhs_ = bsm->ToCompressedRowSparseMatrix();
+      crs_lhs_->set_storage_type(
+          CompressedRowSparseMatrix::StorageType::UPPER_TRIANGULAR);
+    } else {
+      bsm->UpdateCompressedRowSparseMatrix(crs_lhs_.get());
+    }
   } else {
-    lhs = CompressedRowSparseMatrix::FromTripletSparseMatrixTransposed(*tsm);
-    lhs->set_storage_type(
-        CompressedRowSparseMatrix::StorageType::LOWER_TRIANGULAR);
+    if (!crs_lhs_) {
+      crs_lhs_ = bsm->ToCompressedRowSparseMatrixTranspose();
+      crs_lhs_->set_storage_type(
+          CompressedRowSparseMatrix::StorageType::LOWER_TRIANGULAR);
+    } else {
+      bsm->UpdateCompressedRowSparseMatrixTranspose(crs_lhs_.get());
+    }
   }
 
-  *lhs->mutable_col_blocks() = blocks_;
-  *lhs->mutable_row_blocks() = blocks_;
-
   summary.num_iterations = 1;
   summary.termination_type = sparse_cholesky_->FactorAndSolve(
-      lhs.get(), rhs().data(), solution, &summary.message);
+      crs_lhs_.get(), rhs().data(), solution, &summary.message);
   return summary;
 }
 
diff --git a/internal/ceres/schur_complement_solver.h b/internal/ceres/schur_complement_solver.h
index bb2cd88..9f6146b 100644
--- a/internal/ceres/schur_complement_solver.h
+++ b/internal/ceres/schur_complement_solver.h
@@ -186,6 +186,7 @@
   std::vector<Block> blocks_;
   std::unique_ptr<SparseCholesky> sparse_cholesky_;
   std::unique_ptr<BlockRandomAccessDiagonalMatrix> preconditioner_;
+  std::unique_ptr<CompressedRowSparseMatrix> crs_lhs_;
   Vector cg_solution_;
   Vector* scratch_[4] = {nullptr, nullptr, nullptr, nullptr};
 };
diff --git a/internal/ceres/spmv_benchmark.cc b/internal/ceres/spmv_benchmark.cc
index 2a6e679..4211784 100644
--- a/internal/ceres/spmv_benchmark.cc
+++ b/internal/ceres/spmv_benchmark.cc
@@ -183,21 +183,18 @@
   auto bsm_jacobian = CreateFakeBundleAdjustmentJacobian(
       kNumCameras, kNumPoints, kCameraSize, kPointSize, kVisibility, prng);
 
-  CompressedRowSparseMatrix jacobian(bsm_jacobian->num_rows(),
-                                     bsm_jacobian->num_cols(),
-                                     bsm_jacobian->num_nonzeros());
-  bsm_jacobian->ToCompressedRowSparseMatrix(&jacobian);
+  auto jacobian = bsm_jacobian->ToCompressedRowSparseMatrix();
 
   ContextImpl context;
   context.EnsureMinimumThreads(num_threads);
 
-  Vector x(jacobian.num_cols());
-  Vector y(jacobian.num_rows());
+  Vector x(jacobian->num_cols());
+  Vector y(jacobian->num_rows());
   x.setRandom();
   y.setRandom();
   double sum = 0;
   for (auto _ : state) {
-    jacobian.RightMultiplyAndAccumulate(
+    jacobian->RightMultiplyAndAccumulate(
         x.data(), y.data(), &context, num_threads);
     sum += y.norm();
   }
@@ -225,21 +222,18 @@
   std::mt19937 prng;
 
   auto bsm_jacobian = BlockSparseMatrix::CreateRandomMatrix(options, prng);
-  CompressedRowSparseMatrix jacobian(bsm_jacobian->num_rows(),
-                                     bsm_jacobian->num_cols(),
-                                     bsm_jacobian->num_nonzeros());
-  bsm_jacobian->ToCompressedRowSparseMatrix(&jacobian);
+  auto jacobian = bsm_jacobian->ToCompressedRowSparseMatrix();
 
   ContextImpl context;
   context.EnsureMinimumThreads(num_threads);
 
-  Vector x(jacobian.num_cols());
-  Vector y(jacobian.num_rows());
+  Vector x(jacobian->num_cols());
+  Vector y(jacobian->num_rows());
   x.setRandom();
   y.setRandom();
   double sum = 0;
   for (auto _ : state) {
-    jacobian.RightMultiplyAndAccumulate(
+    jacobian->RightMultiplyAndAccumulate(
         x.data(), y.data(), &context, num_threads);
     sum += y.norm();
   }
@@ -258,19 +252,16 @@
   // Perform setup here
   auto bsm_jacobian = CreateFakeBundleAdjustmentJacobian(
       kNumCameras, kNumPoints, kCameraSize, kPointSize, kVisibility, prng);
-  CompressedRowSparseMatrix jacobian(bsm_jacobian->num_rows(),
-                                     bsm_jacobian->num_cols(),
-                                     bsm_jacobian->num_nonzeros());
-  bsm_jacobian->ToCompressedRowSparseMatrix(&jacobian);
+  auto jacobian = bsm_jacobian->ToCompressedRowSparseMatrix();
 
-  Vector x(jacobian.num_rows());
-  Vector y(jacobian.num_cols());
+  Vector x(jacobian->num_rows());
+  Vector y(jacobian->num_cols());
   x.setRandom();
   y.setRandom();
   double sum = 0;
   for (auto _ : state) {
     // This code gets timed
-    jacobian.LeftMultiplyAndAccumulate(x.data(), y.data());
+    jacobian->LeftMultiplyAndAccumulate(x.data(), y.data());
     sum += y.norm();
   }
   CHECK_NE(sum, 0.0);
@@ -291,19 +282,16 @@
   std::mt19937 prng;
 
   auto bsm_jacobian = BlockSparseMatrix::CreateRandomMatrix(options, prng);
-  CompressedRowSparseMatrix jacobian(bsm_jacobian->num_rows(),
-                                     bsm_jacobian->num_cols(),
-                                     bsm_jacobian->num_nonzeros());
-  bsm_jacobian->ToCompressedRowSparseMatrix(&jacobian);
+  auto jacobian = bsm_jacobian->ToCompressedRowSparseMatrix();
 
-  Vector x(jacobian.num_rows());
-  Vector y(jacobian.num_cols());
+  Vector x(jacobian->num_rows());
+  Vector y(jacobian->num_cols());
   x.setRandom();
   y.setRandom();
   double sum = 0;
   for (auto _ : state) {
     // This code gets timed
-    jacobian.LeftMultiplyAndAccumulate(x.data(), y.data());
+    jacobian->LeftMultiplyAndAccumulate(x.data(), y.data());
     sum += y.norm();
   }
   CHECK_NE(sum, 0.0);
@@ -319,10 +307,8 @@
   ContextImpl context;
   std::string message;
   context.InitCuda(&message);
-  CompressedRowSparseMatrix jacobian_crs(
-      jacobian->num_rows(), jacobian->num_cols(), jacobian->num_nonzeros());
-  jacobian->ToCompressedRowSparseMatrix(&jacobian_crs);
-  CudaSparseMatrix cuda_jacobian(&context, jacobian_crs);
+  auto jacobian_crs = jacobian->ToCompressedRowSparseMatrix();
+  CudaSparseMatrix cuda_jacobian(&context, *jacobian_crs);
   CudaVector cuda_x(&context, 0);
   CudaVector cuda_y(&context, 0);
 
@@ -360,10 +346,8 @@
   ContextImpl context;
   std::string message;
   context.InitCuda(&message);
-  CompressedRowSparseMatrix jacobian_crs(
-      jacobian->num_rows(), jacobian->num_cols(), jacobian->num_nonzeros());
-  jacobian->ToCompressedRowSparseMatrix(&jacobian_crs);
-  CudaSparseMatrix cuda_jacobian(&context, jacobian_crs);
+  auto jacobian_crs = jacobian->ToCompressedRowSparseMatrix();
+  CudaSparseMatrix cuda_jacobian(&context, *jacobian_crs);
   CudaVector cuda_x(&context, 0);
   CudaVector cuda_y(&context, 0);
 
@@ -392,10 +376,8 @@
   ContextImpl context;
   std::string message;
   context.InitCuda(&message);
-  CompressedRowSparseMatrix jacobian_crs(
-      jacobian->num_rows(), jacobian->num_cols(), jacobian->num_nonzeros());
-  jacobian->ToCompressedRowSparseMatrix(&jacobian_crs);
-  CudaSparseMatrix cuda_jacobian(&context, jacobian_crs);
+  auto jacobian_crs = jacobian->ToCompressedRowSparseMatrix();
+  CudaSparseMatrix cuda_jacobian(&context, *jacobian_crs);
   CudaVector cuda_x(&context, 0);
   CudaVector cuda_y(&context, 0);
 
@@ -433,10 +415,8 @@
   ContextImpl context;
   std::string message;
   context.InitCuda(&message);
-  CompressedRowSparseMatrix jacobian_crs(
-      jacobian->num_rows(), jacobian->num_cols(), jacobian->num_nonzeros());
-  jacobian->ToCompressedRowSparseMatrix(&jacobian_crs);
-  CudaSparseMatrix cuda_jacobian(&context, jacobian_crs);
+  auto jacobian_crs = jacobian->ToCompressedRowSparseMatrix();
+  CudaSparseMatrix cuda_jacobian(&context, *jacobian_crs);
   CudaVector cuda_x(&context, 0);
   CudaVector cuda_y(&context, 0);
 
diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc
index d2e4ee0..18bb864 100644
--- a/internal/ceres/visibility_based_preconditioner.cc
+++ b/internal/ceres/visibility_based_preconditioner.cc
@@ -392,27 +392,33 @@
 // Compute the sparse Cholesky factorization of the preconditioner
 // matrix.
 LinearSolverTerminationType VisibilityBasedPreconditioner::Factorize() {
-  // Extract the TripletSparseMatrix that is used for actually storing
+  // Extract the BlockSparseMatrix that is used for actually storing
   // S and convert it into a CompressedRowSparseMatrix.
-  const TripletSparseMatrix* tsm =
-      down_cast<BlockRandomAccessSparseMatrix*>(m_.get())->mutable_matrix();
-
-  std::unique_ptr<CompressedRowSparseMatrix> lhs;
+  const BlockSparseMatrix* bsm =
+      down_cast<BlockRandomAccessSparseMatrix*>(m_.get())->matrix();
   const CompressedRowSparseMatrix::StorageType storage_type =
       sparse_cholesky_->StorageType();
   if (storage_type ==
       CompressedRowSparseMatrix::StorageType::UPPER_TRIANGULAR) {
-    lhs = CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm);
-    lhs->set_storage_type(
-        CompressedRowSparseMatrix::StorageType::UPPER_TRIANGULAR);
+    if (!m_crs_) {
+      m_crs_ = bsm->ToCompressedRowSparseMatrix();
+      m_crs_->set_storage_type(
+          CompressedRowSparseMatrix::StorageType::UPPER_TRIANGULAR);
+    } else {
+      bsm->UpdateCompressedRowSparseMatrix(m_crs_.get());
+    }
   } else {
-    lhs = CompressedRowSparseMatrix::FromTripletSparseMatrixTransposed(*tsm);
-    lhs->set_storage_type(
-        CompressedRowSparseMatrix::StorageType::LOWER_TRIANGULAR);
+    if (!m_crs_) {
+      m_crs_ = bsm->ToCompressedRowSparseMatrixTranspose();
+      m_crs_->set_storage_type(
+          CompressedRowSparseMatrix::StorageType::LOWER_TRIANGULAR);
+    } else {
+      bsm->UpdateCompressedRowSparseMatrixTranspose(m_crs_.get());
+    }
   }
 
   std::string message;
-  return sparse_cholesky_->Factorize(lhs.get(), &message);
+  return sparse_cholesky_->Factorize(m_crs_.get(), &message);
 }
 
 void VisibilityBasedPreconditioner::RightMultiplyAndAccumulate(
diff --git a/internal/ceres/visibility_based_preconditioner.h b/internal/ceres/visibility_based_preconditioner.h
index 84773f8..96c9721 100644
--- a/internal/ceres/visibility_based_preconditioner.h
+++ b/internal/ceres/visibility_based_preconditioner.h
@@ -194,6 +194,7 @@
 
   // Preconditioner matrix.
   std::unique_ptr<BlockRandomAccessSparseMatrix> m_;
+  std::unique_ptr<CompressedRowSparseMatrix> m_crs_;
   std::unique_ptr<SparseCholesky> sparse_cholesky_;
 };