The solution is to place the vertex coordinate vectors in the matrix columns i.e.
M = [v_0, v_1, ..., v_n] = [[x_0, x_1, ...., x_n],
[y_0, y_1, ...., y_n],
[z_0, y_1, ...., y_n],
[ 1 , 1 , ...., 1 ]]
Note that we are using homogenous coordinates, hence vectors need to be 4 dimensional where the first three elements are the coordinates (x,y,z), and the fourth is either 1 (if we are multiplying position vectors) or 0 if they are directional vectors (e.g. normals, however, certain considerations need to be taken in if they are directional unit vectors).
So the datablock for M could be constructed as:
M = np.ones([4, len(obj.data.vertices)], dtype=np.float32)
And the columns could be filled in some manor e.g. a possible simple solution is:
for i in len(obj.data.vertices):
M[0:3, i] = obj.data.vertices[i]
Note, there are many ways to fill the 2D array and you might be able to make a better solution (fyi mine may not work as i didn’t test it). For example, your current ‘coords’ array could be transposed and appended with a 4:th row filled with ones, something like:
M[0:3,:] = coords.T
I’m pretty sure above would yield correct result, and i think is more efficient compared to the same solution below (my knowledge of numpy data views is limited):
M[0:3,:] = np.transpose(coords)
Finally, to multiply the data you just use:
world_coords = np.matmul(obj.matrix_world, M)
which will give you the coordinates in world space as defined by the 4x4 matrix: ‘obj.matrix_world’
in the same format as M i.e.
world_coords = [v_0, v_1, ... , v_n]
To get the i:th 3D coordinate vector just do:
v_i = world_coords[0:3, i]
From here i think i hinted at a solution to your current problem, so depending on your next step you could remove the fourth row, creating a 3xn matrix/array. If you need to apply more transformations after modifying the data, keeping the fourth dimension could be preferable as 4x4 transformations could be applied without modyfying the array once again, and could be more efficient. I’m unsure of the extent of the performance gains if numpy is only used to transform the data, as the gain would probably be related to how efficient you can get a matrix on the form of M.