How to fit axis and radius of 3D cylinder？

Once I get some 3-D point coordinates, what algorithm do I use to fit an optimal cylindrical and get the direction vector and radius of the central axis？

My previous idea was to divide a cylinder into layers, and as the number of layers increased, the figure formed by the dots got closer to the cylinder, but in this case I couldn't get an exact radius of the cylinder. (The central axis is obtained by fitting the center of the circle through each layer)

Here is a MCVE to regress axis defined by p0 and p1 (vectors) and radius R (scalar).

First we create some dummy dataset:

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
from scipy.spatial.transform import Rotation

def cylinder(n=60, m=20, r=2, h=5):
    t = np.linspace(0, m * 2 * np.pi, m * n)
    z = np.linspace(0, h, m * n)
    x = r * np.cos(t)
    y = r * np.sin(t)
    return np.stack([x, y, z]).T

X = cylinder()
rot = Rotation.from_rotvec(np.array([-1, 2, 0.5]))
x0 = np.array([1., 2., 0.])
    
X = rot.apply(X)
X = X + x0

Which create a generic use case including origin shift.

Now it is sufficient to write down the geometric equation (see equation 10) as residuals and minimize it by least squares.

def residuals(p, xyz):
    return np.power(np.linalg.norm(np.cross((xyz - p[0:3]), (xyz - p[3:6])), axis=1) / np.linalg.norm((p[3:6] - p[0:3])), 2) - p[6] ** 2

p, res = optimize.leastsq(residuals, x0=[0., 0., 0., 1., 1., 1., 0.], args=(X,), full_output=False)

Which in this case returns:

# array([ -1.8283916 ,  -1.65918186,   3.29901757,  # p0
#         20.31455462,  26.98786514, -22.52837088,  # p1
#          2.        ])                             # R

Graphically it leads to: