Blender cameras convention
I'm assuming that this is not a math.se question,  but a blender.se question, therefore [frame challenge] I'd suggest to use the builtin types and methods.

If you want to know why your math doesn't work, Blender.SE might not be the right place.

I will add some more tips to restructure your code, but take them with grain of salt, as I am no software dev.


I would also discourage the excessive scoping of types from imports at the start of functions. I can see at a glance that np.sin is a numpy function, but have a hard time otherwise.
The x, y, z tuple in the sphere_golden_spiral function was barely readable.   


  PEP: Limit all lines to a maximum of 79 characters

Separate the lookat function from a "position" functions. Right now the lookat function produces a matrix containing the position. This is not very intuitive.
Use enumerate if you want both index and object. (Not range(len(...)).)
If you want to multiply your numpy array times 10, don't multiply the individual elements, multiply them all at once.
Replace look_at with the built-in to_track_quat.
Use the built-in Vector, Quaternion and Matrix types.


A modified script.

import bpy
import numpy as np
import mathutils

def sphere_golden_spiral(num_pts=1000):
    indices = np.arange(0, num_pts, dtype=float) + 0.5

    phi = np.arccos(1 - 2*indices/num_pts)
    theta = np.pi * (1 + 5**0.5) * indices

    x = np.cos(theta) * np.sin(phi)
    y = np.sin(theta) * np.sin(phi)
    z = np.cos(phi)
    points = np.stack((x,y,z),axis=-1)
    return points

def cams_sphere(num_pts):
    points = sphere_golden_spiral(num_pts) * 10

    for i, point in enumerate(points):
        bpy.ops.object.camera_add()
        cam = bpy.context.object
        cam.name='Traj_sphere_%i'%i

        cam.location = mathutils.Vector(point)
        cam.rotation_mode = 'QUATERNION'
        cam.rotation_quaternion = cam.location.to_track_quat('Z', 'Y')

cams_sphere(100)


Obviously, you can use the longer (possibly more stable) version with a matrix multiplication.

        bpy.ops.object.camera_add()
        cam = bpy.context.object
        cam.name='Traj_sphere_%i'%i

        location = mathutils.Vector(point)
        rotation = location.to_track_quat('Z', 'Y')

        mat_loc = mathutils.Matrix.Translation(location)
        mat_rot = rotation.to_matrix().to_4x4()

        cam.matrix_world = mat_loc * mat_rot