sqrtm is too slow for matrices of size 1000

[Vivek Kulkarni]

Hi,
I am implementing spectral clustering for my course work and am using the sqrtm function to find the square root of a matrix . But its far too slow, my matrix is of size (1258,1258). Any suggestions on how I can speed things up or some other function that scipy supports which can be used.

I am implementing the algorithm as described here:
http://books.nips.cc/papers/files/nips14/AA35.pdf

Finding the square root of D^(-1) is way too slow for D^-1 of size (1258,1258).


Thanks.
Vivek.

[Paul Blelloch]

I was interested so I tried this and in my case it just crashed Python.  I think that the crash might be related to an issue that posted earlier related to numpy 1.7.1.  In fact any of the linalg functions in scipy 0.12.0 seem to cause Python to crash on my computer.

[Paul Blelloch]

I don't have an answer, but I can confirm that the algorithm seems to be fairly slow.  I compared a sqrtm or a 1258x1258 random (non-symmetric) matrix in Scipy 0.12.0 and Matlab 2011b on the same machine (weak laptop).  Matlab was 38.5 seconds.  Scipy was 367 seconds.  I did look briefly at the code in both Scipy and Matlab and they both appear to be using an algorithm by by Nicholas J. Higham that's based on a Schur decomposition. I checked the elapsed time for the Schur decomposition itself, and it wasn't very different (~9 seconds on Scipy and 7 seconds in Matlab).  In fact, the Scipy version first returns the real Schur decomposition and then converts that to complex, which only seems to take about 3.5 seconds.  The difference would appear to be somewhere else in the algorithm.  Clearly it can be improved, but that would require some digging.

On Saturday, April 13, 2013 8:47:25 AM UTC-7, Vivek Kulkarni wrote:

[Paul Blelloch]

I looked at the rest of the algorithm, and Scipy and Matlab appear to be doing nearly the same thing.  In Scipy there are three nested for loops as follows:

for j in range(n):

     R[j,j] = sqrt(T[j,j])

     for i in range(j-1,-1,-1):

         s = 0

         for k in range(i+1,j):

             s = s + R[i,k]*R[k,j]

             R[i,j] = (T[i,j] - s)/(R[i,i] + R[j,j])

While in Matlab it's two nested loops as follows:

    for j=1:n
        R(j,j) = sqrt(T(j,j));
        for i=j-1:-1:1
            k = i+1:j-1;
            s = R(i,k)*R(k,j);
            R(i,j) = (T(i,j) - s)/(R(i,i) + R(j,j));
        end
    end

The difference is that Matlab is replacing the inner loop by a Vector multiply of a slice of the matrix.  This is where the speed difference may be.  So if you're feeling adventurous you could try to do the same thing.

On Saturday, April 13, 2013 8:47:25 AM UTC-7, Vivek Kulkarni wrote:

[Paul Blelloch]

I think that I found the problem.  It was in not recognizing that the inner for loop is actually a dot product.  If I replace the following lines of code:

s = 0

for k in range(i+1,j):

    s = s + R[i,k]*R[k,j]

with

s = np.dot(R[i,(i+1):j],R[(i+1):j,j])

Run time decreases from 367 to 15.6 seconds.  My guess is that you could get considerable further speedup, but I'm pleased with the 15.6 seconds.  If you copy the sqrtm function from scipy and make that change I think that you'll see considerable improvement.

On Saturday, April 13, 2013 8:47:25 AM UTC-7, Vivek Kulkarni wrote: