I am trying to write a function that determines if a set of vertices, edges and faces is a pure simplicial complex.

A pure simplicial complex is a set where all facets have the same degree, a facet is a simplex that is not contained in a simplex of higher degree. So for example, a single vertex is a facet. A set of vertices is a pure simplicial complex because all the facets are of degree one and so on.

In this case a simplicial complex is up to order 2 (triangles) so it is expressed as 2 matrices.

An $|edges|\times|vertices|$ adjacency matrix and a $|face|\times|edges|$ adjacency matrix. If you are curious about the details here's a cool resource on simplicial complexes in the context of CS: https://www.cs.cmu.edu/~kmcrane/Projects/DDG/

I represent a subset of the complex as 3 vectors, a vertex vector, an edge vector and a face vector, where every entry is 0 or 1. 1 means the element is in the subset, 0 means it isn't.

Attempted solution

So to determine if a subset is a pure simplicial complex I thought of the following approach. let's say we want to find all the vertex facets first.

0 out all the edge rows in the edge adjacency matrix that are not in the subset. This is achieved by the product $\text{sub_edge_adjacency} = D(edges) \times \text{edge_adjacency}$ where $D$ is an operator that generates a diagonal matrix from a vector. Then multiply the above by a vector filled with 1's of the same dimension as the number of edges and multiply from the right $ones * \text{subedge_adjacency}$. Which is equivalent to just adding all the rows into a single vector...

From then on...

I am not going to keep explaining the algorithm as it should be obvious this is convoluted. I suspect there is a much nicer way to check if a subset is a pure simplicial complex. Essentially I am trying to code the 8th programming problem of chapter 2 in the resource I linked.

I tested my monstrosity and it works, but I am wondering if someone can come up with a cleaner algorithm.


To clarify, the algorithm should be constraint to using mostly, or even exclusively, linear algebra operations. This means that most information should be obtained through the use of matrix products, additions, and queries of the properties of the matrices, like the number of zeroes, the trace of the diagonal etc...

Refined implementation

I managed to simplify the logic quite a bit, but I want to see if it is possible to get a more straighforward implementation:

using SBM = Eigen::SparseVector<bool>;
std::array<SBM, 3> Closure(
    const Eigen::SparseMatrix<bool>& edge_matrix,
    const Eigen::SparseMatrix<bool>& face_matrix,
    const std::array<Eigen::SparseVector<bool>, 3>& simplices)
    const auto face_vertex_adjacency = face_matrix * edge_matrix;
    const auto& vertices = simplices[0];
    const auto& edges = simplices[1];
    const auto& faces = simplice

    SBM closure1_edges = edges + face_matrix.transpose() * faces;
    SBM closure1_vertices = vertices + edge_matrix.transpose() * edges;

    return {closure1_vertices, closure1_edges, faces};
bool IsComplex(
    const Eigen::SparseMatrix<bool>& edge_matrix,
    const Eigen::SparseMatrix<bool>& face_matrix,
    const std::array<Eigen::SparseVector<bool>, 3>& simplices)
    const auto& vertices = simplices[0];
    const auto& edges = simplices[1];
    const auto& faces = simplices[2];

    auto[closure_vertices, closure_edges, closure_faces] =
        Closure(edge_matrix, face_matrix, {vertices, edges, faces});

        closure_vertices.isApprox(vertices) &&
        closure_edges.isApprox(edges) &&

int IsPureComplex(
    const Eigen::SparseMatrix<bool>& edge_matrix,
    const Eigen::SparseMatrix<bool>& face_matrix,
    const std::array<SBM, 3>& simplices)
    const SBM& verts = simplices[0];
    const SBM& edges = simplices[1];
    const SBM& faces = simplices[2];

    if(!IsComplex(edge_matrix, face_matrix, simplices)) return -1;

    SBM point_faces = edges.transpose() * edge_matrix;
    const bool has_point_facets = SBM((verts - point_faces).pruned()).nonZeros();

    SBM line_faces = faces.transpose() * face_matrix;
    const bool has_line_facets = SBM((edges - line_faces).pruned()).nonZeros();
    const bool has_triangle_facets = faces.nonZeros();

    if(int(has_point_facets) + int(has_line_facets) + int(has_triangle_facets) != 1)
        return -1;

    return 3 * has_triangle_facets + 2 * has_line_facets + has_point_facets;

I added all the code so that this is copy pasteable, but the only important function is the last one IsPureComplex. The above works and is better than what I originally had, but I wonder if it is possible to make it even simpler. That is to say, shorter and with less cyclomatic complexity.


First, suppose there is at least one triangle in the set. Now the complex is pure if and only if all vertices are contained in some triangle, because any edge that is a facet contains at least one vertex. This is simple to check: iterate over all triangles and mark all vertices contained in them. If there are no unmarked vertices, the complex is pure. Otherwise, each unmarked vertex lies inside a facet that is not a triangle.

If your set does not contain a triangle, you can use the same approach to determine whether there is facet that is not an edge.

  • $\begingroup$ Iterations should not be necessary this should be achievable exclusively through matrix multiplication. The reason for that is, the major advantage of expressing meshes this way is that you get a major speedup in expressability as complex operations are reduced to matrix products. My current solution operates that way but I am very convinced it can be simplified. $\endgroup$
    – Makogan
    Oct 3 '20 at 15:54
  • $\begingroup$ I see. Perhaps you could edit your question to clarify what the exact constraints on your algorithm should be, so that it is clearer what sort of answer you are looking for. $\endgroup$
    – Discrete lizard
    Oct 3 '20 at 17:28
  • $\begingroup$ You are correct, my bad. I have added the constraints $\endgroup$
    – Makogan
    Oct 3 '20 at 17:32
  • $\begingroup$ @Makogan I'm not exactly sure what operations you have available on matrices and vectors, but if you can count the number of 0's in a vector, then it seems you can implement this idea by multiplying both matrices (or maybe their transpose) with a all-one vector and checking whether the resulting vector has at least one 0 entry. $\endgroup$
    – Discrete lizard
    Oct 3 '20 at 17:53
  • $\begingroup$ I added some sample code to make the question easier to understand and answer. $\endgroup$
    – Makogan
    Oct 3 '20 at 18:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.