I would like to compute a list of random camera positions and corresponding euler rotations on the sphere where at each locations, the camera should be pointing to the object in the center.

Welcome to the board. Do note that these sorts of general math utility questions are better served at Blender Stack Exchange. Higher traffic, more eyeballs, oriented more toward general questions - and many developers here answer there as well. This is a site about specifics of the Blender API. I see that you did post at BSE (how to get camera locations and euler rotations for 3d points hemisphere surface?). I suggest that you consider it your main go-to site for broad-topic Blender python questions.

# Create a camera and Empty by whatever means that suits.
# Choose a radius to suit
RADIUS = 15
# Translate matrix along +Z
xlate = Matrix()
xlate[2][3] = RADIUS
# Camera to position on the hemisphere.
cam = D.objects['Camera']
# Not needed for solution - Just visualizes target
empty = D.objects['Empty']
# randomly choose some orientation; scale to radius
v0 = RADIUS * noise.random_unit_vector()
# Fetch tracking quaternion to align camera to orientation
quat0 = v0.to_track_quat('Z', 'Y')
# Translate camera length of radius, orient camera with tracking quad
cam.matrix_world = quat0.to_matrix().to_4x4() @ xlate

v0 is an arbitrarily oriented sphere radius. By ‘sphere radius’, I mean an axis where I arbitrarily choose one endpoint on the surface of the sphere (RADIUS * noise.random_unit_vector()) and the other endpoint is assumed to be at the center of the sphere ((0, 0, 0)). Generate as many as these to suit; the example just illustrates one sample. For each radius, invoke its to_track_quat() to get a quaternion to align the camera to the radius. The 'Z', 'Y' parameters are so chosen to point the camera in the +Z axis direction, toward the center, and orient the top of the camera in the +Y direction. Compose with a translation matrix to move the camera from the center of the sphere to the surface of the sphere. See screenshots.