diff --git a/src/lib/MeshFEMSparse/BlockCSCHessian.hh b/src/lib/MeshFEMSparse/BlockCSCHessian.hh
index 150cc15..0e0b68b 100644
--- a/src/lib/MeshFEMSparse/BlockCSCHessian.hh
+++ b/src/lib/MeshFEMSparse/BlockCSCHessian.hh
@@ -833,11 +833,18 @@ struct MESHFEM_EXPORT BlockCSCHessian final : public BlockToScalarPolicyDefault<
 
     // Call f(j, loc), passing the location in `Ax` of the diagonal entry
     // in scalar column j, for each j in 0..numScalarVars() - 1
+    // Note: empty block columns are skipped. Their diagonal entries are not
+    // stored (so there is no location to hand to `f`), and
+    // `diagBlockScalarLoc()` would compute an offset into the preceding
+    // column's storage for them. `j` still advances so callers see the correct
+    // scalar column indices. Callers that must touch every diagonal entry
+    // should call `assertAllDiagonalBlocksPresent` first.
     template<class F>
     void visitDiagonalScalarEntries(F &&f) const {
         _Index j = 0;
         for (_Index bj = 0; bj < n; ++bj) {
             auto cs = columnScanner(bj);
+            if (col_nnz(bj) == 0) { j += cs.colBlockSize(); continue; }
             _Index loc = cs.diagBlockScalarLoc();
             for (_Index c_j = 0; c_j < cs.colBlockSize(); ++c_j) {
                 f(j++, loc);
@@ -847,6 +854,14 @@ struct MESHFEM_EXPORT BlockCSCHessian final : public BlockToScalarPolicyDefault<
         }
     }
 
+    // Throw unless every block column has a diagonal block, which the
+    // diagonal-mutating operations below require: for a column without one
+    // there is no stored entry to write.
+    void assertAllDiagonalBlocksPresent(const char *op) const {
+        if (numDiagonalBlocks() < size_t(n))
+            throw std::runtime_error(std::string("BlockCSCHessian::") + op + ": matrix is missing diagonal blocks; insert them first (see BorderedSparseHessian::insertSparsityPatternDiagonalBlocksIfNeeded)");
+    }
+
     virtual Real trace() const override {
         Real result = 0;
         visitDiagonalScalarEntries([&result, this](size_t /* j */, _Index loc) { result += Ax[loc]; });
@@ -1319,14 +1334,17 @@ private:
     VarStructure m_vars;
 
     virtual void m_addDiag(const _Real *d) override {
+        assertAllDiagonalBlocksPresent("addDiag");
         visitDiagonalScalarEntries([d, this](size_t j, _Index loc) { Ax[loc] += d[j]; });
     }
 
     virtual void m_addDiag(_Real d) override {
+        assertAllDiagonalBlocksPresent("addDiag");
         visitDiagonalScalarEntries([d, this](size_t /* j */, _Index loc) { Ax[loc] += d; });
     }
 
     virtual void m_setDiag(_Real d) override {
+        assertAllDiagonalBlocksPresent("setDiag");
         visitDiagonalScalarEntries([d, this](size_t /* j */, _Index loc) { Ax[loc] = d; });
     }
 
diff --git a/src/lib/MeshFEMSparse/SparseMatrices.hh b/src/lib/MeshFEMSparse/SparseMatrices.hh
index e72b20d..18b98ae 100644
--- a/src/lib/MeshFEMSparse/SparseMatrices.hh
+++ b/src/lib/MeshFEMSparse/SparseMatrices.hh
@@ -1171,7 +1171,11 @@ struct CSCMatrix {
                 for (_Index block_j = 0; block_j < blockHsp.n; ++block_j) {
                     size_t gvar_j = block_j * N;
                     size_t numBlocks = blockHsp.Ap[block_j + 1] - blockHsp.Ap[block_j];
-                    bool hasDiagonal = AssumeDiagonalExists || (blockHsp.Ai[blockHsp.Ap[block_j + 1] - 1] == block_j);
+                    // Guard against empty block columns (possible, e.g., for contact
+                    // Hessians, where most vertices are collision-free): reading
+                    // Ai[Ap[block_j + 1] - 1] would be out of bounds, and a spurious
+                    // `hasDiagonal` would underflow `colSize` below.
+                    bool hasDiagonal = (numBlocks > 0) && (AssumeDiagonalExists || (blockHsp.Ai[blockHsp.Ap[block_j + 1] - 1] == block_j));
                     if (hasDiagonal) {
                         size_t colSize = (numBlocks - 1) * N + 1;
                         for (size_t c_j = 0; c_j < N; ++c_j)