Производительность загрузки и хранения разреженных данных с помощью numpy

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

Потребность возникла из-за обработки таких данных: N «изображений» (32x32) по четыре канала каждое. Изображения содержат в среднем ~ 5% ненулевых значений. Поскольку N становится очень большим, хранение всех изображений в ОЗУ неэффективно. Таким образом, сохраняется только количество ненулевых записей, их индексы и значения, а также исходная форма.

Вот пример того, как это можно реализовать:

import numpy as np


def disassemble_data(data):
    # get some dense data and make it sparse
    lengths = np.count_nonzero(data, axis=(1, 2, 3))
    idxs = np.flatnonzero(data)
    vals = data.ravel()[idxs]

    return lengths, idxs, vals, data.shape


def assemble_data(lengths, idxs, vals, shape):
    # get some sparse data and make it dense
    data = np.zeros(shape)

    lower_idx = 0
    for length in lengths:
        upper_idx = lower_idx + length
        data.ravel()[idxs[lower_idx:upper_idx]] = vals[lower_idx:upper_idx]
        lower_idx = upper_idx
    return data


# Create some dummy data 
my_data = np.random.uniform(0, 1, (10, 4, 32, 32))
my_data[my_data > 0.05] = 0

# make data sparse and then dense again
my_new_data = assemble_data(*disassemble_data(my_data))

# assert that this actually works
assert np.allclose(my_data, my_new_data)

Теперь мы можем сразу увидеть преимущество: данные уплотняются изображение за изображением. Это позволяет нам загружать весь набор данных в ОЗУ и генерировать плотные изображения по запросу с помощью генератора:

def image_generator(lengths, idxs, vals, shape):          
        idxs %= np.prod(shape[1:])

        lower_idx = 0
        for length in lengths:
            upper_idx = lower_idx + length
            data = np.zeros(shape[1:])
            data.ravel()[idxs[lower_idx:upper_idx]] = vals[lower_idx:upper_idx]
            lower_idx = upper_idx
            yield data

Кроме того, также можно создавать пакеты изображений:

def image_batch_generator(lengths, idxs, vals, shape, batch_size):
    idxs %= np.prod((batch_size, *shape[1:]))
    lengths = np.sum(lengths.reshape(-1, batch_size), axis=1)

    lower_idx = 0
    for length in lengths:
        upper_idx = lower_idx + length
        data = np.zeros((batch_size, *shape[1:]))
        data.ravel()[idxs[lower_idx:upper_idx]] = vals[lower_idx:upper_idx]
        lower_idx = upper_idx
        yield data

Это довольно удобный подход для моих нужд. И все же мне было интересно, можно ли это ускорить.

Например. Я видел, что набор элементов numpys работает быстрее, чем прямое назначение (согласно документы). Но это работает только для одного элемента, а не для индексных массивов.

Есть ли другие подходы? Я совсем не знаком с cython и им подобным, поэтому буду рад получить несколько советов!

Вы только что заново изобрели scipy.sparse? Единственная разница в том, что он фактически поддерживает вычисления, преобразованные в плотные по мере необходимости.

Marat 09.04.2018 04:11

Это немного зависит от вашего варианта использования. Если вы не хотите использовать разреженные матрицы только в качестве «алгоритма сжатия», но работать только с плотными матрицами, используя набор данных HDF5, сжатие также может быть хорошим и более эффективным способом сделать это. (Это должно обеспечить пропускную способность выше 600 МБ / с) stackoverflow.com/a/48997927/4045774

max9111 09.04.2018 10:03

Я знаю, что есть scipy.sparse, но представленные там форматы не соответствуют тому, который я использую. scipy.sparse предназначен для матриц, а этот метод обобщается на тензоры любой формы. Также большим отличием этого метода является возможность получения срезов разреженного тензора, поэтому нет необходимости сначала раздувать весь тензор. Есть ли лучший способ добиться этого с помощью scipy.sparse? Я использовал HDF5 со сжатием (lz4). При использовании только половины пространства декомпрессия эмпирически занимает больше времени и занимает все ядро ​​процессора. Даже распараллеливание более 8 ядер было медленнее.

Obay 09.04.2018 16:43
Почему в Python есть оператор "pass"?
Почему в Python есть оператор "pass"?
Оператор pass в Python - это простая концепция, которую могут быстро освоить даже новички без опыта программирования.
Некоторые методы, о которых вы не знали, что они существуют в Python
Некоторые методы, о которых вы не знали, что они существуют в Python
Python - самый известный и самый простой в изучении язык в наши дни. Имея широкий спектр применения в области машинного обучения, Data Science,...
Основы Python Часть I
Основы Python Часть I
Вы когда-нибудь задумывались, почему в программах на Python вы видите приведенный ниже код?
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
Алиса и Боб имеют неориентированный граф из n узлов и трех типов ребер:
Оптимизация кода с помощью тернарного оператора Python
Оптимизация кода с помощью тернарного оператора Python
И последнее, что мы хотели бы показать вам, прежде чем двигаться дальше, это
Советы по эффективной веб-разработке с помощью Python
Советы по эффективной веб-разработке с помощью Python
Как веб-разработчик, Python может стать мощным инструментом для создания эффективных и масштабируемых веб-приложений.
1
3
219
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

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

Я немного увеличил скорость ваших функций, дав около 1,4 ГБ / с для сжатия и 1,2 ГБ / с для распаковки, что совсем неплохо. С h5py (blosclz) я смог достичь только около 450 МБ / с, но записал данные также на диск.

Улучшенный разреженный алгоритм

import numpy as np
import numba as nb

#We can use uint16 on (4,32,32), since max. idx<2**16
@nb.jit()
def to_sparse_data_uint16(data):
  data_flat=data.reshape(-1)

  idx=np.empty(data.size,dtype=np.uint16)
  data_out=np.empty(data.size,dtype=data.dtype)

  ii=0
  for i in range(data_flat.shape[0]):
    if (data_flat[i]!=0):
      idx[ii]=i
      data_out[ii]=data_flat[i]
      ii+=1

  return idx[0:ii], data_out[0:ii], data.shape


def to_dense_data(idx,data,shape):
  length=np.prod(shape)
  data_out=np.zeros(length,dtype=data.dtype)

  data_out[idx]=data

  return data_out.reshape(shape)


########################
#do you really need float64 here?
images = np.array(np.random.uniform(0, 1, (100000, 4, 32, 32)),dtype=np.float32)
images[images > 0.05] = 0.

res=[]
t1=time.time()
for i in range(100000):
  res.append(to_sparse_data_uint16(images[i,:,:,:]))

print(time.time()-t1)

t1=time.time()
for i in range(100000):
  data=to_dense_data(res[i][0],res[i][1],res[i][2])

print(time.time()-t1)

Пример HDF5

import numpy as np
import tables #register blosc
import h5py as h5
import h5py_cache as h5c
import time

# I assume here that you don't need float64 for images..
# 1650MB Testdata

images = np.array(np.random.uniform(0, 1, (100000, 4, 32, 32)),dtype=np.float32)
images[images > 0.05] = 0.

#Write data (32,7 GB uncompressed)
hdf5_path='Test.h5'
f = h5c.File(hdf5_path, 'w',chunk_cache_mem_size=1024**2*100) #200 MB cache size
dset_images = f.create_dataset("images", shape=(20*100000, 4, 32, 32),dtype=np.float32,chunks=(1000, 4, 32, 32),compression=32001,compression_opts=(0, 0, 0, 0, 9, 1, 1), shuffle=False)

t1=time.time()
#Don't call h5py to often, this will lead to bad performance
for i in range(20):
    dset_images[i*100000:(i+1)*100000,:,:,:]=images

f.close()
print(time.time()-t1)

print("MB/s: " + str(32700/(time.time()-t1)))

Интересно увидеть numba. Можно ли это применить и к декомпрессии? Предположение о том, что индексы малы, неверно, поскольку 2 ** 16 = 65536, к сожалению, слишком мало для набора данных из нескольких миллионов изображений. Хороший улов на float32!

Obay 09.04.2018 21:19

В примере с декомпрессией я обращаюсь к массиву только один раз (с ограничением памяти). Также обратите внимание, что накладные расходы на вызов функции и добавление списка сильно ухудшают производительность. Сжатие занимает всего (0,65 с, список добавляет еще 0,4 с). Ad uint16: количество изображений здесь не проблема, поскольку idx начинается с нуля на каждом изображении. Если изображения не увеличиваются, этого должно быть достаточно.

max9111 09.04.2018 23:30

Кажется, я что-то упускаю, разве цикл не увеличивает i до data.size? Итак, i = data.size-1 на последней итерации. Затем i присваивается idx [ii]. Так что это приведет к переполнению (или усечению?), Насколько я могу судить!

Obay 09.04.2018 23:52

Тип данных рекламного изображения: обычно используются типы данных для изображений uint8, а для некоторых изображений - uint16, для которых датчик не может превышать 10–14 бит на канал. Обычно 16 бит, потому что это легче обрабатывать, но это также можно оптимизировать ... Проверьте, что здесь действительно достаточно. Пожалуйста, также очистите шаблон доступа и, возможно, не очень случайные ненулевые места и значения между ними (это также может привести к распространенным алгоритмам сжатия без потерь на переднем плане). Ваши изображения также могут быть (или нет) коррелированы с другими. Не существует самого быстрого \ лучшего алгоритма сжатия без знания точных данных.

max9111 09.04.2018 23:52

Изображения принимаются несколькими матрицами фотоумножителей и должны быть откалиброваны, поэтому остаются дробные интенсивности пикселей (то есть количество фотоэлектронов) с высоким динамическим диапазоном. Вы совершенно правы насчет корреляции. Однако я уже пробовал использовать сжатое хранилище фрагментов HDF5 с LZ4 (из пакета 'hdf5plugin'), это было примерно в 20 раз медленнее, чем разреженный формат. Я попробую установить blosc и размер блока.

Obay 10.04.2018 00:04

В вашем примере каждое изображение (4x32x32), поэтому idx не может выходить за рамки 4 * 32 * 32-1 = 4095 (и ii не за 4096). Как сказано выше, uint16 подходит только в том случае, если np.prod (dims_of_single_image) <2 ** 16. Пожалуйста, извините, у вас, очевидно, нет 3 каналов ... Использование python-blosc для прямого сжатия фрагментов изображений также возможно и может достигать скорости более 2-3 ГБ / с на четырехъядерном процессоре и намного выше, чем на мощной рабочей станции. Я могу завтра добавить пример. github.com/Blosc/python-blosc

max9111 10.04.2018 00:07

Ааа, я вижу. Я полностью неправильно прочитал ваш код, извините! Ваш код предназначен для сжатия каждого изображения по отдельности и сохранения их в списке. Мой код берет все изображения сразу и сохраняет их все в одном большом массиве. Затем я просто сохраняю length, idxs, vals, data.shape с numpy.savez в одном файле и получаю невероятно высокую скорость чтения. Я добавлю тест позже.

Obay 10.04.2018 00:14

Другие вопросы по теме