Need help vectorizing some Python code

Question

I have some code that I need help vectorizing. I want to convert the following to vector form, how can I? I want to get rid of the inner loop - apparently, it's possible to do so. X is an NxD matrix. y is a 1xD vector.

def foo(X, y, mylambda, N, D, epsilon): 
... 
    for j in xrange(D): 
        aj = 0 
        cj = 0 
        for i in xrange(N): 
            aj += 2 * (X[i,j] ** 2) 
            cj += 2 * (X[i,j] * (y[i] - w.transpose()*X[i].transpose() + w[j]*X[i,j])) 

...

If I call numpy.vectorize() on the function, it throws an error at runtime.

Complete code:

import scipy
import scipy.io
import numpy
from numpy import linalg
from scipy import *

def data(N, d, k, sigma, seed=12231):
    random.seed(seed)
    X = randn(N, d)
    wg = zeros(1 + d)
    wg[1:k + 1] = 10 * sign(randn(k))
    eps = randn(N) * sigma
    y = X.dot(wg[1:]) + wg[0] + eps
    return (y, X)


def foo(X, y, mylambda, n, D, epsilon):
    identityMatrix = numpy.matrix(numpy.identity(D))

    w = (X.transpose() * X + mylambda * identityMatrix).getI() * X.transpose() * y
    newweight = (X.transpose() * X + mylambda * identityMatrix).getI() * X.transpose() * y

    iterate = 1
    iteration = 0

    while iterate > 0 and iteration < 10000:
        iteration += 1
        iterate = 0
        maxerror = 0
        for j in xrange(D):
            aj = 0
            cj = 0
            for i in xrange(n):
                aj += 2 * (X[i,j] ** 2)
                cj += 2 * (X[i,j] * (y[i] - w.transpose()*X[i].transpose() + w[j]*X[i,j]))

            if cj < -mylambda:
                newweight[j,0] = (cj + mylambda)/ aj
            elif cj > mylambda:
                newweight[j,0] = (cj - mylambda)/ aj
            else:
                newweight[j,0] = 0

            if abs(newweight[j,0] - w[j,0]) > epsilon:
                iterate += 1
            if abs(newweight[j,0] - w[j,0]) > maxerror:
                maxerror = abs(newweight[j,0] - w[j,0])            
            w[j,0] = newweight[j,0]

N, D, k = 50, 75, 5
(y, X) = data(N, D, k, 1, 123)
X = numpy.matrix(X)
y = numpy.matrix(y).transpose()
foo(X, y, 1, N, D, 0.1)

Error: File "C:\Python27\lib\site-packages\numpy\lib\function_base.py", line 1572, in call return self._vectorize_call(func=func, args=vargs) File "C:\Python27\lib\site-packages\numpy\lib\function_base.py", line 1632, in _vectorize_call ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args) File "C:\Python27\lib\site-packages\numpy\lib\function_base.py", line 1596, in _get_ufunc_and_otypes outputs = func(*inputs) File "<stdin>", line 21, in simplelasso IndexError: 0-d arrays can only use a single () or a list of newaxes (and a single ...) as an index — user3210796
– user3210796, Commented Jan 18, 2014 at 21:09
Do you mean a[j] = 0 and a[j] += 2 * (X[i,j] **2) (and a similar change of bj to b[j]? Otherwise, in the outer loop you overwrite aj without ever using it. — Warren Weckesser
– Warren Weckesser, Commented Jan 18, 2014 at 21:30
As you can tell from the comments, there are plenty of people who desperately want to help you, but you're making it hard for them. Try updating your question with a short, self contained, correct example, I suspect you'll get much more useful answers instead of confused comments. — Bi Rico
– Bi Rico, Commented Jan 18, 2014 at 23:15

HYRY · Accepted Answer · 2014-01-19 00:46:43Z

3

You can replace:

aj = 0
cj = 0
for i in xrange(n):
    aj += 2 * (X[i,j] ** 2)
    cj += 2 * (X[i,j] * (y[i] - w.transpose()*X[i].transpose() + w[j]*X[i,j]))

with:

aj = 2*np.sum(X[:,j].T*X[:,j])
cj = 2*np.sum(np.multiply(X[:, j].T, (y.T - w.T*X.T + w[j] * X[:, j].T)))

answered Jan 19, 2014 at 0:46

HYRY

97.8k28 gold badges197 silver badges192 bronze badges

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

Collectives™ on Stack Overflow

Need help vectorizing some Python code

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related