Как оптимизировать построение очень большого массива с помощью numpy или scipy - PullRequest
Новая
       29

Как оптимизировать построение очень большого массива с помощью numpy или scipy

0 голосов
/ 14 июня 2019

Я обрабатываю очень большие наборы данных на Python 64 бит и мне нужна помощь для оптимизации моего кода интерполяции.

Я привык использовать numpy, чтобы избежать циклов, но здесь есть 2 цикла, которые я не могунайти способ избежать.

Основная проблема также в том, что размер массивов, которые мне нужно вычислить, выдает ошибку памяти, когда я использую numpy, поэтому я переключился на scipy разреженные массивы, которые работают, но занимают слишком много временичтобы вычислить 2 левых цикла ...

Я пытался итеративно построить свою матрицу, используя функцию numpy.fromfunction, но она не запустилась, потому что размер массива слишком велик.

У меня естьЯ уже прочитал много постов о построении больших массивов, но запрашиваемые массивы были слишком просты по сравнению с тем, что я должен построить, поэтому решения здесь не работают.

Я не могу уменьшить размер данныхустановить, поскольку это облако точек, которое я уже выложил плиткой размером 10x10.

Вот мой код интерполяции: 

z_int = ss.dok_matrix((x_int.shape))

n,p = x_obs.shape
m = y_obs.shape[0]
a1 = ss.coo_matrix( (n, 3), dtype=np.int64 )
a2 = ss.coo_matrix( (3, 3), dtype=np.int64 )
a3 = ss.dok_matrix( (n, m))
a4 = ss.coo_matrix( (3, n), dtype=np.int64)

b = ss.vstack((z_obs, ss.coo_matrix( (3, 1), dtype=np.int64 ))).tocoo()

a1 = ss.hstack((ss.coo_matrix(np.ones((n,p))), ss.coo_matrix(x_obs), ss.coo_matrix(y_obs)))

shape_a3 = a3.shape[0]

for l in np.arange(0, shape_a3):
    for c in np.arange(0, shape_a3) :
        if l == c:
            a3[l, c] = rho
        else:
            a3[l, c] = phi(x_obs[l] - x_obs[c], y_obs[l] - y_obs[c])


a4 = a1.transpose()

a12 = ss.vstack((a1, a2))
a34 = ss.vstack((a3, a4))
a = ss.hstack((a12, a34)).tocoo()


x = spsolve(a, b)


for i in np.arange(0, z_int.shape[0]):
    for j in np.arange(0, z_int.shape[0]):
        z_int[i, j] = x[0] + x[1] * x_int[i, j] + x[2] * y_int[i, j] + np.sum(x[3:] * phi(x_int[i, j] - x_obs, y_int[i, j] - y_obs).T)


return z_int.todense()

, где dist () - функцияэто вычисляет расстояние, и phi выглядит следующим образом: 

return dist(dx, dy) ** 2 * np.log(dist(dx, dy))

Мне нужен код, чтобы работать быстрее, и я знаю, что он может быть написан очень плохо, но я хотел бы узнать, как написать что-то более оптимизированное дляулучшить свои навыки кодирования.

1 Ответ

0 голосов
/ 14 июня 2019

Этот код сложен для понимания и понятен медленно. Итерация на разреженных матрицах даже медленнее, чем итерация на плотных массивах. Я почти хотел бы, чтобы вы начали с небольшого рабочего примера с использованием плотных массивов, прежде чем беспокоиться о том, чтобы он работал для большого случая. Я не собираюсь пробовать всестороннее исправление или ускорение, просто грызу тут и там.

Это первое a1 создание ничего не делает для вас (кроме траты времени). Python не является скомпилированным языком, где вы определяете тип переменных в начале. a1 после второго назначения - разреженная матрица, потому что это то, что hstack создано, а не из-за предыдущего coo назначения. 

a1 = ss.coo_matrix( (n, 3), dtype=np.int64 )
...
a1 = ss.hstack((ss.coo_matrix(np.ones((n,p))), ss.coo_matrix(x_obs), ss.coo_matrix(y_obs)))

Инициализация матриц dok, zint и a3 правильная, потому что вы выполняете итерацию для заполнения значений. Но мне нравится видеть такую ​​инициализацию ближе к циклу, а не назад наверху. Я бы использовал lil вместо dok, но я не уверен, что это быстрее. 

for l in np.arange(0, shape_a3):
    for c in np.arange(0, shape_a3) :
        if l == c:
            a3[l, c] = rho
        else:
            a3[l, c] = phi(x_obs[l] - x_obs[c], y_obs[l] - y_obs[c])

Тест l==c определяет основную диагональ. Существуют способы изготовления диагональных матриц. Но похоже, что вы устанавливаете все элементы a3. Если так, зачем использовать более медленный разреженный подход?

Что такое phi Требуются ли скалярные входы? x_obs[:,None]-x_obs должен давать массив (n, n) напрямую.

Что производит spsolve? x, разреженный или плотный. Из вашего использования в цикле z_int это выглядит как 1d плотный массив. Похоже, вы устанавливаете все значения z_int.

Если phi принимает массив (n, n), я думаю 

x[0] + x[1] * x_int[i, j] + x[2] * y_int[i, j] + np.sum(x[3:] * phi(x_int[i, j] - x_obs, y_int[i, j] - y_obs).T)

x[0] + x[1] * x_int + x[2] * y_int +  np.sum(x[3:]) * phi(x_int-x_obs, y_int-y_obs).T)
...