Unexpected result setting a shear matrix as transform
Mesh.transform

Can transform mesh directly by a matrix with Mesh.transform(matrix)  The matrix must be a 4x4.  Use in object mode. Update the mesh to see result.

Python console code, C = bpy.context

&gt;&gt;&gt; mat = Matrix.Shear('XY', 4, (3,3))
&gt;&gt;&gt; mat
Matrix(((1.0, 0.0, 3.0, 0.0),
        (0.0, 1.0, 3.0, 0.0),
        (0.0, 0.0, 1.0, 0.0),
        (0.0, 0.0, 0.0, 1.0)))

&gt;&gt;&gt; C.object.data.transform(mat)
&gt;&gt;&gt; C.object.data.update()



  But this is transforms the vertices for good. In the modal operator, I
  need to apply it to the original positions of the sphere on each
  update. And the matrix is not necessarily invertable, so I can't just
  do the inverse beforehand, which makes it tricky


Shear matrix is not orthogonal, this is why it is not seen as an object matrix, and only in edit mode.

&gt;&gt;&gt; mat.is_orthogonal
False


AFAIK will need to apply to verts to shear.

Some ways to copy mesh coords, swap out a mesh copy

&gt;&gt;&gt; me = C.object.data
&gt;&gt;&gt; me2 = C.object.data.copy()
&gt;&gt;&gt; me2.transform(mat)
&gt;&gt;&gt; C.object.data = me2 # see transform
&gt;&gt;&gt; C.object.data = me # back to original
&gt;&gt;&gt; bpy.data.meshes.remove(me2) # remove it when done


Foreach set and get

load the 3d vert coordinates into a ravelled array.  Have used numpy as the matrix algebra could be performed on it reshaped, although will find Mesh.transform is quick too.

&gt;&gt;&gt; import numpy as np
&gt;&gt;&gt; data = np.zeros(len(me.vertices) * 3)
&gt;&gt;&gt; me.vertices.foreach_get("co", data)
&gt;&gt;&gt; data
array([ 0.00000000e+00,  1.95090324e-01,  9.80785251e-01, ...,
        1.51397259e-08,  1.95089921e-01, -9.80785310e-01])


The data is ravelled, here it is reshaped to display the coordinates 

&gt;&gt;&gt; np.reshape(data, (-1, 3))
array([[ 0.00000000e+00,  1.95090324e-01,  9.80785251e-01],
       [ 0.00000000e+00,  3.82683456e-01,  9.23879504e-01],
       [ 0.00000000e+00,  5.55570245e-01,  8.31469595e-01],
       ...,
       [ 1.52975474e-07,  7.07106352e-01, -7.07106769e-01],
       [ 7.84696610e-08,  3.82682920e-01, -9.23879623e-01],
       [ 1.51397259e-08,  1.95089921e-01, -9.80785310e-01]])


Transform the data, update to see result

&gt;&gt;&gt; me.transform(mat)
&gt;&gt;&gt; me.update()


write back, update to see result.

&gt;&gt;&gt; me.vertices.foreach_set("co", data)
&gt;&gt;&gt; me.update()


Another OTT method would be to make a shape key and transform that.