Gaussian blurring a function defined over a surface mesh I have a scalar function defined over the vertices of a surface mesh. I want to compute an approximate (and generalized, I suppose) "Gaussian blur/convolution" of this function over the surface.

I can imagine for each vertex, taking an average of the function at vertices in it's (multi-hop) neighborhood weighted according to the Gaussian of their Euclidean distance (in R3) from the current vertex multiplied by the local area they represent. This approximation would be good if the surface isn't too curved at the scale of the Gaussian kernel. However it would still be computationally prohibitive.

Is there a more efficient convolving algorithm, perhaps based on some kind of iterated "flow" or partial exchange through edges of the mesh?

The mesh can have either triangular or quadrilateral faces - whichever is more convenient in answering this question. But the edge lengths are not constant, so local neighborhood geometry differs between vertices.
A few different approaches

I'll consider a few variations on your specific request, since you mention efficiency and I suspect your specific request may be the least efficient. I'll also suggest ways of improving the efficiency without varying from your intended approach, so you can weigh up the alternatives.

Blurring the volume instead of the surface

If you want the distance metric to be 3D Euclidean distance instead of 2D Euclidean distance within the surface, then you could perform the blur on a regular 3D grid to which the scalar function you have in mind has been applied. Then you can use the final result of your Difference of Gaussians to calculate the scalar values at the vertices of your irregular mesh. This avoids having to take into account the mesh shape for the bulk of the calculation.

The 3D grid is likely to have a much larger number of vertices than the 2D mesh, but they will all be equally spaced and the large kernel blur can be achieved by repeated application of a small kernel blur taking into account only the 6 nearest neighbours, which will always be at a constant distance away. This approach involves potentially more calculation, but the ease with which the regular grid could be GPU accelerated may appeal.

This will give a different result from performing the blur on the mesh vertices using 3D Euclidean distance. For example, the 3D grid approach will be affected by distinctive regions of the 3D scalar function that are near but not on the mesh. This may be desirable or not depending on your specific purpose.

Using 2D distance instead of 3D distance

If you find that you need the distance metric to be 2D Euclidean within the surface, then you can get a good approximation to a larger kernel Gaussian blur by repeatedly applying a smaller kernel Gaussian blur. If there is not too much variation in the edge lengths within your mesh you may be able to choose a kernel size which allows for only including vertices one edge away at each iteration. This allows for only using single edge lengths to calculate the size of the contribution of a vertex, rather than calculating a 2D multi-edge distance.

3D distance using the surface without the volume

If you need the calculation to be precisely as described in your question, being calculated within the mesh rather than within the surrounding volume, but also using the 3D Euclidean distance, then using nearest neighbours and several iterations will not work. Unless the mesh is near to flat, the repeated application of a nearest neighbour blur will result in an approximation to the 2D Euclidean distance case, since the values will only be able to bleed from vertex to vertex, not directly along the shortest path as they would in a single pass. This will give less spread than would be achieved by a single pass that calculates the 3D distance to a vertex 10 edges away. (I have used 10 edges since you mention a 10-hop in your comment on the question.)

Implementing the blur in a single pass will mean calculating the 3D Euclidean distance between every vertex and every other vertex within a 10 edge radius. This will be expensive, but perfectly possible. Since you mention efficiency, consider that there are some redundancies you can eliminate provided you have sufficient available memory.

The two blurs that you produce prior to taking the Difference of Gaussians will use the same set of 3D distances up to the edge radius of the smaller kernel blur. If you can save these then you only need to calculate them once, rather than once per blur.

Also, each distance will be used twice per blur - once in each direction as the length from vertex A to vertex B is the same as the length from B to A. Caching/memoising these distances will avoid calculating them twice.

Effects from arbitrarily many edges away

If the surface curves such that some vertices which are many edges away are still near enough to affect each other in 3D distance, then rather than considering vertices within a certain number of edges away, you may need to consider all vertices within a certain 3D radius, regardless of how long the path via edges. In this case you can consider fewer vertices by using space partitioning, choosing a specific method which suits the mesh.

If you don't want parts of the surface which approach each other to influence each other, then you probably want a 2D distance metric rather than the 3D one.

If you have a wide range of different edge lengths then you may find the same problem of not being able to define a set number of edges to traverse, even if the mesh is fairly flat. Again you may need to define a 3D distance instead of a number of edges, and consider all vertices that lie within that radius.