<p dir="auto">Hi,<br>
I would like to find the best-fit axis of points that are on a cylindrical surface.  Seems that <code class="notranslate">scipy.linalg.svd</code> is the function to look for.</p>
<p dir="auto">So to test out, I decide to generate some points, function <code class="notranslate">makeCylinder</code>, from this thread <a href="https://stackoverflow.com/questions/22285994/how-to-generate-regular-points-on-cylindrical-surface" rel="nofollow">https://stackoverflow.com/questions/22285994/how-to-generate-regular-points-on-cylindrical-surface</a>, and estimate the axis.</p>
<p dir="auto">This is the code:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="    def makeCylinder(radius, length, nlength, alpha, nalpha, center, orientation): 
    
        # Load 
        from numpy import array, allclose, linspace, tile, vstack 
        from numpy import pi, cos, sin, arccos, cross 
        from numpy.linalg import norm 
        from angoli import rotMatrixAxisAngle 
        
        # Create the length array 
        I = linspace(0, length, nlength) 
        
        # Create alpha array avoid duplication of endpoints 
        if int(alpha) == 360: 
            A = linspace(0, alpha, num=nalpha, endpoint=False)/180.0*pi 
        else: 
            A = linspace(0, alpha, num=nalpha)/180.0*pi 
        
        # Calculate X and Y 
        X = radius * cos(A) 
        Y = radius * sin(A) 
        
        # Tile/repeat indices so all unique pairs are present 
        pz = tile(I, nalpha) 
        px = X.repeat(nlength) 
        py = Y.repeat(nlength) 
        
        # Points 
        points = vstack(( pz, px, py )).T 
        ## Shift to center 
        points += array(center) - points.mean(axis=0) 
        
        # Orient tube to new vector 
        ovec = orientation / norm(orientation) 
        cylvec = array([1,0,0]) 
        
        if allclose(cylvec, ovec): 
            return points 
        
        # Get orthogonal axis and rotation 
        oaxis = cross(ovec, cylvec) 
        rot = arccos(ovec.dot(cylvec)) 
        
        R = rotMatrixAxisAngle(oaxis, rot) 
        
        return points.dot(R) 

    from numpy.linalg import norm 
    from numpy.random import rand 
    from scipy.linalg import svd 
    
    for i in xrange(100): 
        orientation = rand(3)
        orientation[0] = 0 
        orientation /= norm(orientation) 
    
        # Generate sample points 
        points = makeCylinder(radius = 3.0, 
                              length = 20.0, nlength = 20, 
                              alpha = 360, nalpha = 30, 
                              center = [0,0,0], 
                              orientation = orientation) 
        # Least Square 
        uu, dd, vv = svd(points - points.mean(axis=0)) 
        asse = vv[0] 
    
        assert abs( abs(orientation.dot(asse)) - 1) &lt;= 1e-4, orientation.dot(asse) "><pre class="notranslate"><code class="notranslate">    def makeCylinder(radius, length, nlength, alpha, nalpha, center, orientation): 
    
        # Load 
        from numpy import array, allclose, linspace, tile, vstack 
        from numpy import pi, cos, sin, arccos, cross 
        from numpy.linalg import norm 
        from angoli import rotMatrixAxisAngle 
        
        # Create the length array 
        I = linspace(0, length, nlength) 
        
        # Create alpha array avoid duplication of endpoints 
        if int(alpha) == 360: 
            A = linspace(0, alpha, num=nalpha, endpoint=False)/180.0*pi 
        else: 
            A = linspace(0, alpha, num=nalpha)/180.0*pi 
        
        # Calculate X and Y 
        X = radius * cos(A) 
        Y = radius * sin(A) 
        
        # Tile/repeat indices so all unique pairs are present 
        pz = tile(I, nalpha) 
        px = X.repeat(nlength) 
        py = Y.repeat(nlength) 
        
        # Points 
        points = vstack(( pz, px, py )).T 
        ## Shift to center 
        points += array(center) - points.mean(axis=0) 
        
        # Orient tube to new vector 
        ovec = orientation / norm(orientation) 
        cylvec = array([1,0,0]) 
        
        if allclose(cylvec, ovec): 
            return points 
        
        # Get orthogonal axis and rotation 
        oaxis = cross(ovec, cylvec) 
        rot = arccos(ovec.dot(cylvec)) 
        
        R = rotMatrixAxisAngle(oaxis, rot) 
        
        return points.dot(R) 

    from numpy.linalg import norm 
    from numpy.random import rand 
    from scipy.linalg import svd 
    
    for i in xrange(100): 
        orientation = rand(3)
        orientation[0] = 0 
        orientation /= norm(orientation) 
    
        # Generate sample points 
        points = makeCylinder(radius = 3.0, 
                              length = 20.0, nlength = 20, 
                              alpha = 360, nalpha = 30, 
                              center = [0,0,0], 
                              orientation = orientation) 
        # Least Square 
        uu, dd, vv = svd(points - points.mean(axis=0)) 
        asse = vv[0] 
    
        assert abs( abs(orientation.dot(asse)) - 1) &lt;= 1e-4, orientation.dot(asse) 
</code></pre></div>
<p dir="auto">As you can see, I generate multiple cylinder whose axis is random (rand(3)).</p>
<p dir="auto">The funny thing is that <code class="notranslate">svd</code> returns an axis that is absolutely perfect if the first component of <code class="notranslate">orientation</code> is zero (<code class="notranslate">orientation[0] = 0</code>).</p>
<p dir="auto">If I comment this line the estimated axis is way off.</p>
<p dir="auto">Even using leastsq on a cylinder equation returns the same behavior:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="    def bestLSQ1(points): 
        from numpy import array, sqrt 
        from scipy.optimize import leastsq 
    
        # Expand 
        points = array(points) 
        x = points[:,0] 
        y = points[:,1] 
        z = points[:,2] 
    
        # Calculate the distance of each points from the center (xc, yc, zc) 
        # http://geometry.puzzles.narkive.com/2HaVJ3XF/geometry-equation-of-an-arbitrary-orientated-cylinder
        def calc_R(xc, yc, zc, u1, u2, u3): 
            return sqrt( (x-xc)**2 + (y-yc)**2 + (z-zc)**2 - ( (x-xc)*u1 + (y-yc)*u2 + (z-zc)*u3 )**2 ) 
    
        # Calculate the algebraic distance between the data points and the mean circle centered at c=(xc, yc, zc) 
        def dist(c): 
            Ri = calc_R(*c) 
            return Ri - Ri.mean() 
    
        # Axes - Minimize residu 
        xM, yM, zM = points.mean(axis=0) 
        # Calculate the center 
        center, ier = leastsq(dist, (xM, yM, zM, 0, 0, 1)) 
        xc, yc, zc, u1, u2, u3 = center 
        asse = u1, u2, u3 
    
        return asse 
"><pre class="notranslate"><code class="notranslate">    def bestLSQ1(points): 
        from numpy import array, sqrt 
        from scipy.optimize import leastsq 
    
        # Expand 
        points = array(points) 
        x = points[:,0] 
        y = points[:,1] 
        z = points[:,2] 
    
        # Calculate the distance of each points from the center (xc, yc, zc) 
        # http://geometry.puzzles.narkive.com/2HaVJ3XF/geometry-equation-of-an-arbitrary-orientated-cylinder
        def calc_R(xc, yc, zc, u1, u2, u3): 
            return sqrt( (x-xc)**2 + (y-yc)**2 + (z-zc)**2 - ( (x-xc)*u1 + (y-yc)*u2 + (z-zc)*u3 )**2 ) 
    
        # Calculate the algebraic distance between the data points and the mean circle centered at c=(xc, yc, zc) 
        def dist(c): 
            Ri = calc_R(*c) 
            return Ri - Ri.mean() 
    
        # Axes - Minimize residu 
        xM, yM, zM = points.mean(axis=0) 
        # Calculate the center 
        center, ier = leastsq(dist, (xM, yM, zM, 0, 0, 1)) 
        xc, yc, zc, u1, u2, u3 = center 
        asse = u1, u2, u3 
    
        return asse 

</code></pre></div>
<p dir="auto">This is the version that I am using:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
('0.19.0', '1.12.1', sys.version_info(major=2, minor=7, micro=13, releaselevel='final', serial=0))"><pre class="notranslate"><code class="notranslate">import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
('0.19.0', '1.12.1', sys.version_info(major=2, minor=7, micro=13, releaselevel='final', serial=0))
</code></pre></div>
<p dir="auto">Thanks</p>
<p dir="auto">My issue is about a possible memory leak in scipy.spatial.distance.pdist and scipy.spatial.distance.cdist</p>
<h4 dir="auto">Reproducing code example:</h4>

<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="import numpy as np
from scipy.spatial.distance import pdist
X = np.random.uniform(0, 100, size=10000).reshape(-1, 10)
for _ in range(10000):
    test = pdist(X/0.01, metric='sqeuclidean')
    "><pre class="notranslate"><code class="notranslate">import numpy as np
from scipy.spatial.distance import pdist
X = np.random.uniform(0, 100, size=10000).reshape(-1, 10)
for _ in range(10000):
    test = pdist(X/0.01, metric='sqeuclidean')
    
</code></pre></div>
<h4 dir="auto">Error Description:</h4>
<p dir="auto">On my Mac OSX 10.14.6 the memory consumption of the snipped above quickly rises into the Gigabytes (after a few seconds). If I remove the division from the pdist argument (see below), the problem vanishes. (I don't think this is intended behavior, as it also leads to excessive memory consumption in the Gaussian Process Regression implementation of Sklearn, where they have a similar line in the implementation of the RBF kernel)</p>
<p dir="auto">The code, which does not produce this problem (division removed):</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="import numpy as np
from scipy.spatial.distance import pdist
X = np.random.uniform(0, 100, size=10000).reshape(-1, 10)
for _ in range(10000):
    test = pdist(X, metric='sqeuclidean')"><pre class="notranslate"><code class="notranslate">import numpy as np
from scipy.spatial.distance import pdist
X = np.random.uniform(0, 100, size=10000).reshape(-1, 10)
for _ in range(10000):
    test = pdist(X, metric='sqeuclidean')
</code></pre></div>
<h4 dir="auto">Scipy/Numpy/Python version information:</h4>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="1.7.0 1.21.0 sys.version_info(major=3, minor=8, micro=5, releaselevel='final', serial=0)"><pre class="notranslate"><code class="notranslate">1.7.0 1.21.0 sys.version_info(major=3, minor=8, micro=5, releaselevel='final', serial=0)
</code></pre></div>