5

My Setup: Python 2.7.4.1, Numpy MKL 1.7.1, Windows 7 x64, WinPython

Context:

I tried to implement the Sequential Minimal Optimization algorithm for solving SVM. I use maximal violating pair approach.

The problem:

In working set selection procedure i want to find maximum value of gradient and its index for elements which met some condition, y[i]*alpha[i]<0 or y[i]*alpha[i] 

#y - array of -1 and 1
y=np.array([-1,1,1,1,-1,1])
#alpha- array of floats in range [0,C]
alpha=np.array([0.4,0.1,1.33,0,0.9,0])
#grad - array of floats
grad=np.array([-1,-1,-0.2,-0.4,0.4,0.2])


GMaxI=float('-inf')       
GMax_idx=-1        
n=alpha.shape[0] #usually n=100000
C=4
B=[0,0,C]
for i in xrange(0,n):
    yi=y[i]  #-1 or 1
    alpha_i=alpha[i]

    if (yi * alpha_i< B[yi+1]): # B[-1+1]=0 B[1+1]=C
        if( -yi*grad[i]>=GMaxI):
            GMaxI= -yi*grad[i]
            GMax_idx = i

This procedure is called many times (~50000) and profiler shows that this is the bottleneck. It is possible to vectorize this code?

Edit 1: Add some small exemplary data

Edit 2: I have checked solution proposed by hwlau , larsmans and Mr E. Only solutions proposed Mr E is correct. Below sample code with all three answers:

import numpy as np
y=np.array([   -1,  -1,   -1,  -1,   -1,   -1,   -1,   -1])
alpha=np.array([0,  0.9,  0.4, 0.1, 1.33,    0,  0.9,    0])
grad=np.array([-3, -0.5,   -1,  -1, -0.2,   -4, -0.4, -0.3])
C=4
B=np.array([0,0,C])

#hwlau - wrong index and value
filter = (y*alpha < C*0.5*(y+1)).astype('float')
GMax_idx = (filter*(-y*grad)).argmax()
GMax = -y[GMax_idx]*grad[GMax_idx]

print GMax_idx,GMax

#larsmans - wrong index
neg_y_grad = (-y * grad)[y * alpha < B[y + 1]]
GMaxI = np.max(neg_y_grad)
GMax_ind = np.argmax(neg_y_grad)

print GMax_ind,GMaxI

#Mr E - correct result
BY = np.take(B, y+1)
valid_mask = (y * alpha < BY)
values = -y * grad
values[~valid_mask] = np.min(values) - 1.0

GMaxI = values.max()
GMax_idx = values.argmax()
print GMax_idx,GMaxI

Output (GMax_idx, GMaxI)
0 -3.0 
3 -0.2
4 -0.2

Conclusions

After checking all solutions, the fastest one (2x-6x) is solution proposed by @ali_m. However it requires to install some python packages: numba and all its prerequisites.

I have some trouble to use numba with class methods, so I create global functions which are autojited with numba, my solution look something like this:

from numba import autojit  

@autojit
def FindMaxMinGrad(A,B,alpha,grad,y):
    '''
    Finds i,j indices with maximal violatin pair scheme
    A,B - 3 dim arrays, contains bounds A=[-C,0,0], B=[0,0,C]
    alpha - array like, contains alpha coeficients
    grad - array like, gradient
    y - array like, labels
    '''
    GMaxI=-100000
    GMaxJ=-100000

    GMax_idx=-1
    GMin_idx=-1

    for i in range(0,alpha.shape[0]):

        if (y[i] * alpha[i]< B[y[i]+1]):
            if( -y[i]*grad[i]>GMaxI):
                GMaxI= -y[i]*grad[i]
                GMax_idx = i

        if (y[i] * alpha[i]> A[y[i]+1]):
            if( y[i]*grad[i]>GMaxJ):
                GMaxJ= y[i]*grad[i]
                GMin_idx = i

    return (GMaxI,GMaxJ,GMax_idx,GMin_idx)  

class SVM(object):

    def working_set(self,....):
        FindMaxMinGrad(.....)
Improve this question
4
  • You could use list comprehension to start with - but I doubt it'd make it execute much faster (if any faster at all). Maybe try looking at in-built Numpy functions to do some of the work. Commented Dec 19, 2013 at 12:49
  • Can you please add some input data of the right shape and dtype as Python code? e.g. you can either type it out by hand for a small value of n or generate random data of the right shape with NumPy's random module. Make it so that we can directly run this code without guessing the input and output values. Then you'll get a useful answer in no time. Commented Dec 19, 2013 at 12:56
  • 1
    I have edit my question and add some data. Commented Dec 19, 2013 at 13:15
  • Indeed, I had a bug in the argmax. Fixed that; still only three lines of code. Commented Dec 19, 2013 at 19:26

4 Answers 4

Reset to default
3

You can probably do quite a lot better than plain vectorization if you use numba to JIT-compile your original code that used nested loops.

import numpy as np
from numba import autojit

@autojit
def jit_max_grad(y, alpha, grad, B):
    maxgrad = -inf
    maxind = -1
    for ii in xrange(alpha.shape[0]):
        if (y[ii] * alpha[ii] < B[y[ii] + 1]):
            g = -y[ii] * grad[ii]
            if g >= maxgrad:
                maxgrad = g
                maxind = ii
    return maxind, maxgrad

For comparison, here's Mr E's vectorized version:

def mr_e_max_grad(y, alpha, grad, B):

    BY = np.take(B, y+1)
    valid_mask = (y * alpha < BY)
    values = -y * grad
    values[~valid_mask] = np.min(values) - 1.0
    GMaxI = values.max()
    GMax_idx = values.argmax()
    return GMax_idx, GMaxI

Timing:

y = np.array([   -1,  -1,   -1,  -1,   -1,   -1,   -1,   -1])
alpha = np.array([0,  0.9,  0.4, 0.1, 1.33,    0,  0.9,    0])
grad = np.array([-3, -0.5,   -1,  -1, -0.2,   -4, -0.4, -0.3])
C = 4
B = np.array([0,0,C])

%timeit mr_e_max_grad(y, alpha, grad, B)
# 100000 loops, best of 3: 19.1 µs per loop

%timeit jit_max_grad(y, alpha, grad, B)
# 1000000 loops, best of 3: 1.07 µs per loop

Update: if you want to see what the timings look like on bigger arrays, it's easy to define a function that generates semi-realistic fake data based on your description in the question:

def make_fake(n, C=4):
    y = np.random.choice((-1, 1), n)
    alpha = np.random.rand(n) * C
    grad = np.random.randn(n)
    B = np.array([0,0,C])
    return y, alpha, grad, B

%%timeit y, alpha, grad, B = make_fake(100000, 4)
mr_e_max_grad(y, alpha, grad, B)
# 1000 loops, best of 3: 1.83 ms per loop

%%timeit y, alpha, grad, B = make_fake(100000, 4)
jit_max_grad(y, alpha, grad, B)
# 1000 loops, best of 3: 471 µs per loop
Improve this answer
Sign up to request clarification or add additional context in comments.

3 Comments

Thank you, very useful suggestion. I wonder how the timings look like when the arrays are bigger.
Actually this function is a method in a class 'SVM'. I have tried to use your approach, but decorator @autojit don't work. I'm not very familiar with numba and don't know how to add specific decorators for such function. The timings you shown are very encouraging. In addition in the same loop I have to find min value, so I think this solution is the best one.
Take a look at this: numba.pydata.org/numba-doc/0.6/doc/extension_types.html
3

I think this is a fully vectorized version

import numpy as np

#y - array of -1 and 1
y=np.array([-1,1,1,1,-1,1])
#alpha- array of floats in range [0,C]
alpha=np.array([0.4,0.1,1.33,0,0.9,0])
#grad - array of floats
grad=np.array([-1,-1,-0.2,-0.4,0.4,0.2])


BY = np.take(B, y+1)
valid_mask = (y * alpha < BY)
values = -yi * grad
values[~valid_mask] = np.min(values) - 1.0

GMaxI = values.max()
GMax_idx = values.argmax()
Improve this answer

1 Comment

Thank you for your time and effort, I will check and let you know about timing results in my context.
2

Here you go:

y=np.array([-1,1,1,1,-1,1])
alpha=np.array([0.4,0.1,1.33,0,0.9,0])
grad=np.array([-1,-1,-0.2,-0.4,0.4,0.2])
C=4

filter = (y*alpha < C*0.5*(y+1)).astype('float')
GMax_idx = (filter*(-y*grad)).argmax()
GMax = -y[GMax_idx]*grad[GMax_idx]

No benchmark tried, but it is pure numerical and vectorized so it should be fast.

Improve this answer

1 Comment

Nice and compact solution, I will check all three answers and let you know about timing results in my context.
2

If you change B from a list to a NumPy array, you can at least vectorize the yi * alpha_i< B[yi+1] and push the loop inwards:

GMaxI = float('-inf')      
GMax_idx = -1      
for i in np.where(y * alpha < B[y + 1])[0]:
    if -y[i] * grad[i] >= GMaxI:
        GMaxI= -y[i] * grad[i]
        GMax_idx = i

That should save a bit of time. Next up, you can vectorize -y[i] * grad[i]:

GMaxI = float('-inf')
GMax_idx = -1
neg_y_grad = -y * grad
for i in np.where(y * alpha < B[y + 1])[0]:
    if neg_y_grad[i] >= GMaxI:
        GMaxI= -y[i] * grad[i]
        GMax_idx = i

Finally, we can vectorize away the entire loop by using max and argmax on -y * grad, filtered by y * alpha < B[y + 1]:

neg_y_grad = (-y * grad)
GMaxI = np.max(neg_y_grad[y * alpha < B[y + 1]])
GMax_idx = np.where(neg_y_grad == GMaxI)[0][0]
Improve this answer

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.