vectorizing a function of an array of arrays

Question

I am new to vectorizing code, and I am really psyched about how much faster everything is, but I can't get the high speed out of this particular piece of code...

Here is the housing class...

class GaussianMixtureModel:
    def __init__(self, image_matrix, num_components, means=None):

        self.image_matrix = image_matrix
        self.num_components = num_components
        if(means is None):
            self.means = np.zeros(num_components)
        else:
            self.means = np.array(means)
        self.variances = np.zeros(num_components)
        self.mixing_coefficients = np.zeros(num_components)

And here is what I've got so far that works:

def likelihood(self):
    def g2(x):

#N =~ 5
#self.mixing_coefficients = 1D, N items
#self.variances = 1D, N items
#self.means = 1D, N items

        mc = self.mixing_coefficients[:,None,None]
        std = self.variances[:,None,None] ** 0.5
        var = self.variances[:,None,None]
        mean = self.means[:,None,None]
        return np.log((mc*(1.0/(std*np.sqrt(2.0*np.pi)))*(np.exp(-((x-mean)**2.0)/(2.0*var)))).sum())

    f = np.vectorize(g2)

#self.image_matrix =~ 400*700 2D matrix

    log_likelihood = (f(self.image_matrix)).sum()
    return log_likelihood

And here is what I've got that gives a strange result (note that self.image_matrix is an nxn matrix of a grayscale image):

def likelihood(self):
    def g2():

#N =~ 5
#self.mixing_coefficients = 1D, N items
#self.variances = 1D, N items
#self.means = 1D, N items
#self.image_matrix = 1D, 400x700 2D matrix

        mc = self.mixing_coefficients[:,None,None]
        std = self.variances[:,None,None] ** 0.5
        var = self.variances[:,None,None]
        mean = self.means[:,None,None]
        return np.log((mc*(1.0/(std*np.sqrt(2.0*np.pi)))*(np.exp(-((self.image_matrix-mean)**2.0)/(2.0*var)))).sum())

    log_likelihood = (g2()).sum()
    return log_likelihood

However, the second version is really fast compared to the first (which takes almost 10 seconds...and speed is really important here, because this is part of a convergence algorithm)

Is there a way to replicate the results of the first version and the speed of the second? (And I'm really not familiar enough with vectorizing to know why the second version isn't working)

What are the dimensions of the other arrays involved? self.mixing_coefficients, self.variances, and self.means? (And where is self coming from?) — user2357112
– user2357112, Commented Nov 14, 2015 at 23:15
I'll edit the answer to include the dims, but the overall class is huge...i'll try to get the important stuff — Chris
– Chris, Commented Nov 14, 2015 at 23:26
In your second version you only seem to be using the top left element in your image matrix (self.image_matrix[0,0]) — ali_m
– ali_m, Commented Nov 14, 2015 at 23:53

user2357112 · Accepted Answer · 2015-11-14 23:53:52Z

2

The second version is so fast because it only uses the first cell of self.image_matrix:

return np.log((mc*(1.0/(std*np.sqrt(2.0*np.pi)))*(np.exp(-((self.image_matrix[0,0]-mean)**2.0)/(2.0*var)))).sum())
#                                                                            ^^^^^

This is also why it's completely wrong. It's not actually a vectorized computation over self.image_matrix at all. Don't try to use its runtime as a point of comparison; you can always make wrong code faster than right code.

By eliminating the use of np.vectorize, you can make the first version much faster, but not as fast as the wrong code. The sum inside the log simply needs the appropriate axis specified:

def likelihood(self):
    def f(x):
        mc = self.mixing_coefficients[:,None,None]
        std = self.variances[:,None,None] ** 0.5
        var = self.variances[:,None,None]
        mean = self.means[:,None,None]
        return np.log((mc*(1.0/(std*np.sqrt(2.0*np.pi)))*(np.exp(-((x-mean)**2.0)/(2.0*var)))).sum(axis=0))

    log_likelihood = (f(self.image_matrix)).sum()

This can be further simplified and optimized in a few ways. For example, the nested function can be eliminated, and multiplying by 1.0/whatever is slower than dividing by whatever, but eliminating np.vectorize is the big thing.

answered Nov 14, 2015 at 23:53

user2357112

286k32 gold badges490 silver badges569 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

Chris Over a year ago

no, that was an error in copying over...I corrected it...but I think the axis will work-thank you!

Chris Over a year ago

and...that is good to know about vectorize...the speed increase from changing the axis is exactly what I was expecting. Thanks again

Collectives™ on Stack Overflow

vectorizing a function of an array of arrays

1 Answer 1

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related