Как подтолкнуть цикл for к numpy

У меня есть следующий фрагмент кода, который делает именно то, что я хочу (это часть метода кригинга). Но проблема в том, что он идет слишком медленно, и я хочу знать, есть ли какая-нибудь возможность перевести цикл for в numpy? Если я вытолкну numpy.sum и использую там аргумент оси, он немного ускорится, но, видимо, это не узкое место. Любые идеи о том, как я могу нажать forloop на numpy, чтобы ускорить его, или другие способы ускорить его?)

# n = 2116
print GRZVV.shape  # (16309, 2116)
print GinvVV.shape  # (2117, 2117) 
VVg = numpy.empty((GRZVV.shape[0]))

for k in xrange(GRZVV.shape[0]):
    GRVV = numpy.empty((n+1, 1))
    GRVV[n, 0] = 1
    GRVV[:n, 0] = GRZVV[k, :]
    EVV = numpy.array(GinvVV * GRVV)  # GinvVV is numpy.matrix
    VVg[k] = numpy.sum(EVV[:n, 0] * VV)

Я разместил размеры матрицы ndarrays n, чтобы прояснить некоторые вещи.

редактировать: форма ВВ 2116


person usethedeathstar    schedule 23.09.2013    source источник
comment
Если VV.shape == (16309,), как вы можете умножить его на EVV[:n, 0], который имеет форму (n,)?   -  person askewchan    schedule 23.09.2013
comment
Возможно, в последней строке вашего цикла должно быть EVV[:n, 0] * VV[k], что, похоже, и предполагает ответ @Jaime.   -  person askewchan    schedule 23.09.2013
comment
@askewchan форма была 2116, я немного запутался, я отредактировал ее,   -  person usethedeathstar    schedule 23.09.2013


Ответы (2)


Вы можете сделать следующее вместо вашего цикла по k (время выполнения ~ 3 с):

tmp = np.concatenate((GRZVV, np.ones((16309,1),dtype=np.double)), axis=1)
EVV1 = np.dot(GinvVV, tmp.T)
#Changed line below based on *askewchan's* recommendation
VVg1 = np.sum(np.multiply(EVV1[:n,:],VV[:,np.newaxis]), axis=0)
person Joel Vroom    schedule 23.09.2013
comment
Это дает тот же результат, что и код @usethedeathstar, и работает на моей машине в 15 раз быстрее. - person askewchan; 23.09.2013
comment
Нет необходимости в вызове плитки, так как np.multiply вещает. Измените его на: VVg1 = np.sum(np.multiply(EVV1[:n,:],VV[:,np.newaxis]), axis=0) для небольшого ускорения. - person askewchan; 23.09.2013
comment
+1 Гораздо быстрее, чем np.einsum для больших массивов, не уверен, что понимаю почему... - person Jaime; 23.09.2013
comment
небольшая проблема: VVg1 по-прежнему является numpy.matrix вместо ndarray? это влияет на np.multiply или около того? Я предполагаю/надеюсь, что нет, или я должен после np.dot сделать приведение к numpy.array на EVV1? - person usethedeathstar; 25.09.2013
comment
Еще один побочный вопрос: в чем разница между [:,numpy.newaxis] и использованием [:,None]? И почему numpy.multiply вместо того, чтобы просто использовать * ? И ускорение в 12 раз (что ХОРОШО) - person usethedeathstar; 25.09.2013
comment
Разница между np.newaxis и None --› stackoverflow.com/questions/944863/ - person Joel Vroom; 25.09.2013
comment
np.multiply использовался для поэлементного умножения в np.matrix (stackoverflow.com/questions/4151128/) - person Joel Vroom; 25.09.2013
comment
@joelVroom, так что в основном я могу снова передать его в ndarray и покончить с этим после numpy.dot? - person usethedeathstar; 26.09.2013

Вы в основном берете каждую строку GRZVV, добавляете 1 в конце, умножаете ее на GinvVV, а затем складываете все элементы в векторе. Если бы вы не делали «добавить 1», вы могли бы сделать все это без циклов, как:

VVg = np.sum(np.dot(GinvVV[:, :-1], GRZVV.T), axis=-1) * VV

или даже:

VVg = np.einsum('ij,kj->k', GinvVV[:, :-1], GRZVV) * VV

Как нам справиться с этим дополнительным 1? Что ж, результирующий вектор, полученный в результате умножения матриц, будет увеличен на соответствующее значение в GinvVV[:, -1], а когда вы их все сложите, значение будет увеличено на np.sum(GinvVV[:, -1]). Таким образом, мы можем просто вычислить это один раз и добавить ко всем элементам возвращаемого вектора:

VVg = (np.einsum('ij,kj->k', GinvVV[:-1, :-1], GRZVV) + np.sum(GinvVV[:-1, -1])) * VV

Приведенный выше код работает, если VV является скаляром. Если это массив формы (n,), то будет работать следующее:

GinvVV = np.asarray(GinvVV)
VVgbis = (np.einsum('ij,kj->k', GinvVV[:-1, :-1]*VV[:, None], GRZVV) +
          np.dot(GinvVV[:-1, -1], VV))
person Jaime    schedule 23.09.2013
comment
С редактированием @usethedeathstar VV.shape теперь равно 2116, поэтому оно не транслируется в вашем решении (поскольку Wg.shape равно 16309) - person askewchan; 23.09.2013
comment
Подумал, что это скаляр. С полным массивом, не сворачивая массив до конца, как в ответе Джоэла, это намного быстрее, чем для больших массивов, описанное выше. - person Jaime; 23.09.2013