Find the closest converging point of a group of vectors

I am trying to find the point that is closest to a group of vectors.

For context, the vectors are inverted rays emitted from center of aperture stop after exiting a lens, this convergence is meant to locate the entrance pupil.

The backward projection of the exiting rays, while not converging at one single point due to spherical aberration, is quite close to converging toward a point, as illustrated in the figure below.

(For easier simulation the positive z is pointing down)

I believe that to find the closest point, I would be finding a point that has the shortest distance to all these lines. And wrote the method as follow:

def FindConvergingPoint(position, direction):

    A = np.eye(3) * len(direction) - np.dot(direction.T, direction)
    b = np.sum(position - np.dot(direction, np.dot(direction.T, position)), axis=0)

    return np.linalg.pinv(A).dot(b)

For the figure above and judging visually, I would have expected the point to be something around [0, 0, 20]

However, this is not the case. The method yielded a result of [ 0., 188.60107764, 241.13690715], which is far from the converging point I was expecting.

Is my algorithm faulty or have I missed something about python/numpy?

Attached are the data for the vectors:

position = np.array([
    [0, 0, 0],
    [0, -1.62, 0.0314],
    [0, -3.24, 0.1262],
    [0, -4.88, 0.2859],
    [0, -6.53, 0.5136],
    [0, -8.21, 0.8135],
    [0, -9.91, 1.1913],
    [0, -11.64, 1.6551],
    [0, -13.43, 2.2166],
    [0, -15.28, 2.8944],
    [0, -17.26, 3.7289]
])

direction = np.array([
    [0, 0, 1],
    [0, 0.0754, 0.9972],
    [0, 0.1507, 0.9886],
    [0, 0.2258, 0.9742],
    [0, 0.3006, 0.9537],
    [0, 0.3752, 0.9269],
    [0, 0.4494, 0.8933],
    [0, 0.5233, 0.8521],
    [0, 0.5969, 0.8023],
    [0, 0.6707, 0.7417],
    [0, 0.7459, 0.6661]
])

I'm assuming the answer to my question in the comment is that you have symmetry such that the "closest point" must lie on the z-axis. Furthermore, I'm assuming you are somewhat flexible about the notion of "closest".

First, let's remove the zeroth point from position and direction, since that will pass through all points on the z-axis (and cause numerical problems).

Next, let's find the distance from each position along direction that brings us back to the z-axis.

# transpose and remove the x-axis for simplicity
position = position.T[1:]
direction = direction.T[1:]

# more carefully normalize direction vectors
direction /= np.linalg.norm(direction, axis=0)

# find the distance from `position` along `direction` that
# brings us back to the z-axis.
distance = position[0] / direction[0]

# All these points should be on the z-axis
point_on_z = position - distance * direction
# array([[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
#          0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
#        [21.45665199, 21.380772  , 21.34035527, 21.23103513, 21.09561354,
#         20.89001607, 20.608748  , 20.26801397, 19.79193392, 19.14234148]])

Indeed, the y-coordinate of all these points is 0, and the z-coordinate is nearly 20 for all of them, as you expected. If you are flexible about your notion of "closest", the mean of the z-coordinates will minimize the sum of square of distances to the point where the rays intersect the z-axis (but not necessarily the minimum distances between your point and the rays).

This might not solve exactly the problem you were hoping to solve, but hopefully you will find it "useful". If it's not spot on, let me know what assumptions I've made that I shouldn't, and we can try to find a solution to the refined question.