How can I efficiently duplicate a large amount of objects?
Remove repeated operator calls
Most likely this is a veiled duplicate of
Python performance with Blender operators
Removing 17,000 operator calls will help... ie removing bpy.ops.object.duplicate_move using low level API calls
# copy the object
copy = ob.copy()
# copy and link data, prob not needed here
mesh_copy = copy.data.copy()
copy.data = mesh_copy
# set the local location
copy.location = (0, 0, 0)
# and / or set the global translation
copy.matrix_world.translation = (0, 0, 0) # scene origin

Link to scene in &lt;= 2.7x
scene.objects.link(ob)

Link to collection in &gt;= 2.80
context.collection.objects.link(ob)

In 2.8 for example move_obj can become, (Note there is no need to link them to collection here, could do it later)
def dupe_move_ob(coll, ob, vec):
    copy = ob.copy()
    copy.matrix_world = Matrix.Translation(vec) @ copy.matrix_world
    coll.objects.link(copy)
    return copy

Returns  a reference to the newly created and moved object.
Similarly would change other methods to pass objects
def rotate_object(ob, rot_mat):

Dupli
For 17, 000 objects will still be slowish. As commented consider using some form of dupli for this, eg dupliverts, make one mesh, add verts at all the locations calculated, set their normals (via the rotation) and one object will be duplicated for each vertex aligned to the vertex normal.
Make the new mesh object the parent of dupli object and under   Properties &gt; Object &gt; Instancing choose Verts.
To illustrate a script to duplivert a cone to a subsurfed suzanne.

The vertex normals of suzanne are determined by the faces. If we construct a vert only mesh, then set v.normal = (x, y, z) for each vertex
import bpy
from math import radians
from mathutils import Matrix

context = bpy.context

bpy.ops.mesh.primitive_monkey_add()
monkey = context.object
monkey.modifiers.new(name=&quot;Foo&quot;, type='SUBSURF')

bpy.ops.mesh.primitive_cone_add()
cone = context.object
cone.scale *= 0.1
cone.parent = monkey
cone.data.transform(Matrix.Translation((0, 0, 1)))
cone.data.transform(Matrix.Rotation(radians(-90), 4, 'X'))

monkey.instance_type = 'VERTS'
monkey.use_instance_vertices_rotation = True