Поиск:

Ответ в темуСоздание новой темы Создание опроса
> библиотеки для матричных вычислений 
:(
    Опции темы
mrgloom
Дата 14.11.2012, 13:43 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Опытный
**


Профиль
Группа: Участник
Сообщений: 829
Регистрация: 8.6.2011

Репутация: нет
Всего: нет



Суть такова, что мне надо делать PCA для данных которые не влезают в память.
собсвенно получается, что либо нужна функция которая находит собсвенные векторы собсвенные значения не загружая всё в память.
либо аналогичная реализация SVD.(т.к. вроде есть реализация PCA через SVD).


попробовал библиотеку http://arma.sourceforge.net по дефолту используется blas lapack, но может подключать еще и другие например Intel MKL и еще какие то.

ну на тесте 

Код

mat A = randu<mat>(10000,200*200);
cx_vec eigval;
сx_mat eigvec;
eig_gen(eigval,eigvec,A,'r'); 


всё благополучно упало, ссылаясь на нехватку памяти.


вообщем существуют ли в природе такие библиотеки?
PM MAIL   Вверх
W4FhLF
Дата 14.11.2012, 19:44 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


found myself
****


Профиль
Группа: Участник Клуба
Сообщений: 2831
Регистрация: 2.12.2006

Репутация: 5
Всего: 121



1. Какая размерность данных?
2. Нужны все компоненты или только доминантные?



--------------------
"Бог умер" © Ницше
"Ницше умер" © Бог
PM ICQ   Вверх
Pavia
Дата 14.11.2012, 20:16 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Опытный
**


Профиль
Группа: Участник
Сообщений: 418
Регистрация: 6.12.2008

Репутация: 11
Всего: 12



W4FhLF, 
Нужна не просто библиотека. А что-бы там была защита от дурака.  (mrgloom, без обид)
А тоже бывает задашь параметр 10 000, а прога рассчитана на 10. Конечно можно использовать жёсткий диск, но время расчёта может оказаться очень долгим. Нужно чтобы библиотека говорила столько памяти будет использовано и сколько времени будет считаться.

Добавлено через 5 минут и 25 секунд
mrgloom, 
eig_gen более 100 не потянет по точности.
PM MAIL   Вверх
W4FhLF
Дата 14.11.2012, 22:02 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


found myself
****


Профиль
Группа: Участник Клуба
Сообщений: 2831
Регистрация: 2.12.2006

Репутация: 5
Всего: 121



Цитата(Pavia @  14.11.2012,  20:16 Найти цитируемый пост)
W4FhLF, 
Нужна не просто библиотека. А что-бы там была защита от дурака.  (mrgloom, без обид)
А тоже бывает задашь параметр 10 000, а прога рассчитана на 10. Конечно можно использовать жёсткий диск, но время расчёта может оказаться очень долгим. Нужно чтобы библиотека говорила столько памяти будет использовано и сколько времени будет считаться.


У меня была однажды задача, где я считал SVD матрицы (dense matrix причём) размером 60000х60000. Считал в обычном matlab. Бывают такие задачи. Другой вопрос, требует ли задача mrgloom'a в действительности таких ресурсов. И прежде, чем бросаться в поиски каких-то экзотических библиотек нужно осмыслить алгоритм. 

Для PCA считается ков. матрица, которая симметричная и positive semidefinite, соответственно SVD будет давать точно такой же результат, как и EVD. При том, что SVD сам по себе дороже, как в плане памяти, так и по-времени. Рассчёт ков. матрицы для EVD не требует загрузки в память данных целиком, можно использовать накопительный алгоритм. Поэтому я и спросил про размерность данных, т.к. если просто много реализаций, а размерность относительно небольшая, проблемы никакой нет и задача решается просто.

Если же размерность большая, тогда да, есть проблема. В таком случае, если требуются только доминантные компоненты (скажем первые 100-200), тогда надо смотреть в сторону iterative PCA (гугл в помощь), которые позволяют строить ортогональный базис в подпространстве Крылова без использования дорогих EVD. В целом iterative PCA всегда оправдан в случае необходимости использовать только доминантные компоненты.  



Это сообщение отредактировал(а) W4FhLF - 14.11.2012, 22:12


--------------------
"Бог умер" © Ницше
"Ницше умер" © Бог
PM ICQ   Вверх
mrgloom
Дата 15.11.2012, 10:30 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Опытный
**


Профиль
Группа: Участник
Сообщений: 829
Регистрация: 8.6.2011

Репутация: нет
Всего: нет



размерность допустим 100,000 х 1,000 , но вопрос даже не в этом,
а в том чтобы например программа работала на компьютере с 1Гб оперативной памяти.

почему падает Armadillo? возможно потому что ХР х32 или потому что не умеет использовать файл подкачки?
eig_gen вроде как выводит еще для подсчитанных значений погрешности при выводе в консоль?

Если не удастся подобрать алгоритм, то хорошо бы повесить всё это на систему распределения памяти в винде,
матлаб написал вроде out of memory, но видимо ему надо поставить своп на диске побольше, вообщем еще попробую.

теперь про алгоритмическую часть:
Из ответа я так понял, что если сэмплов много, а размерность относительно малая, то это не проблема?
cov в матлабе не справится, надо самому переписывать?

про iterative PCA уже слышал, мне как раз надо первые k максимальных собственных значений,
 только тут непонятен вопрос с точностью, т.е. будет ли давать такой же ответ как обычный алгоритм.
PM MAIL   Вверх
Pavia
Дата 15.11.2012, 11:00 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Опытный
**


Профиль
Группа: Участник
Сообщений: 418
Регистрация: 6.12.2008

Репутация: 11
Всего: 12



mrgloom, 
Цитата(mrgloom @  15.11.2012,  10:30 Найти цитируемый пост)
размерность допустим 100,000 х 1,000 , но вопрос даже не в этом,

400 МБайт или 800 МБайт в зависимости от типа
Цитата(mrgloom @  14.11.2012,  13:43 Найти цитируемый пост)
(10000,200*200);

1600 МБайт 3200МБайт.
XP может выделить непрерывный кусок 500МБайт максимум 1500МБайт. 
Вот оно и падает. 

Цитата

Из ответа я так понял, что если сэмплов много, а размерность относительно малая, то это не проблема?

Да матрицу ковариации считаешь она будет MxM  где M размерность. N число сэмплов не участвует.

Цитата

, мне как раз надо первые k максимальных собственных значений

Для поиска k максимальных собственных значений есть свои алгортмы. Про точность непомню. 
PM MAIL   Вверх
W4FhLF
Дата 15.11.2012, 11:07 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


found myself
****


Профиль
Группа: Участник Клуба
Сообщений: 2831
Регистрация: 2.12.2006

Репутация: 5
Всего: 121



Цитата(mrgloom @  15.11.2012,  10:30 Найти цитируемый пост)
размерность допустим 100,000 х 1,000 , но вопрос даже не в этом,


То есть 100000 реализаций в пространстве размерностью 1000?

Цитата(mrgloom @  15.11.2012,  10:30 Найти цитируемый пост)
Из ответа я так понял, что если сэмплов много, а размерность относительно малая, то это не проблема?


Если семплов 100000, а размерность 1000, то для рассчёта ков. матрицы и последующего ортогонального базиса потребуется ~3*1000^2 * 8 = 24 мегабайт памяти. 

Цитата(mrgloom @  15.11.2012,  10:30 Найти цитируемый пост)
 только тут непонятен вопрос с точностью, т.е. будет ли давать такой же ответ как обычный алгоритм. 


От спектра матрицы зависит. В подавляющем большинстве задач точности этой достаточно.

Добавлено через 9 минут и 22 секунды
При этом т.н. накопительный алгоритм рассчёта ков. матрицы позволяет добавлять новые реализации и обновлять ков. матрицу без пересчёта с нуля. Для этого достаточно хранить накопленные суммы. 

Если дин. диапазон значений большой, то данные можно центрировать и нормировать. 



--------------------
"Бог умер" © Ницше
"Ницше умер" © Бог
PM ICQ   Вверх
mrgloom
Дата 8.11.2013, 12:47 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Опытный
**


Профиль
Группа: Участник
Сообщений: 829
Регистрация: 8.6.2011

Репутация: нет
Всего: нет



Написал блочное перемножение для больших матриц на питоне, которое не требует много RAM,используя hdf5, единственное непонятно как его оптимизировать под конкретный компьютер и будет ли это работать лучше чем тот же своп системы(но во всяком случае более контролируемо и прозрачно), в питоне например используя numpy.memmap были проблемы(и вообще получается это более ограниченнывй подход завязанный на систему), а как работает матлаб со свопом надо еще протестировать.


Код

import numpy as np
import tables
import time
n_row=1000
n_col=1000
n_batch=100
def test_hdf5_disk():
    rows = n_row
    cols = n_col
    batches = n_batch
    #settings for all hdf5 files
    atom = tables.Float32Atom() #if store uint8 less memory?
    filters = tables.Filters(complevel=9, complib='blosc') # tune parameters
    Nchunk = 4*1024  # ?
    chunkshape = (Nchunk, Nchunk)
    chunk_multiple = 1
    block_size = chunk_multiple * Nchunk
   
    fileName_A = 'carray_A.h5'
    shape_A = (n_row*n_batch, n_col)  # predefined size
    h5f_A = tables.open_file(fileName_A, 'w')
    A = h5f_A.create_carray(h5f_A.root, 'CArray', atom, shape_A, chunkshape=chunkshape, filters=filters)
    for i in range(batches):
        data = np.random.rand(n_row, n_col)
        A[i*n_row:(i+1)*n_row]= data[:]
    rows = n_col
    cols = n_row
    batches = n_batch
    fileName_B = 'carray_B.h5'
    shape_B = (rows, cols*batches)  # predefined size
    h5f_B = tables.open_file(fileName_B, 'w')
    B = h5f_B.create_carray(h5f_B.root, 'CArray', atom, shape_B, chunkshape=chunkshape, filters=filters)
    sz= rows/batches
    for i in range(batches):
        data = np.random.rand(sz, cols*batches)
        B[i*sz:(i+1)*sz]= data[:]
    fileName_C = 'CArray_C.h5'
    shape = (A.shape[0], B.shape[1])
    h5f_C = tables.open_file(fileName_C, 'w')
    C = h5f_C.create_carray(h5f_C.root, 'CArray', atom, shape, chunkshape=chunkshape, filters=filters)
    sz= block_size
    t0= time.time()
    for i in range(0, A.shape[0], sz):
        for j in range(0, B.shape[1], sz):
            for k in range(0, A.shape[1], sz):
                C[i:i+sz,j:j+sz] += np.dot(A[i:i+sz,k:k+sz],B[k:k+sz,j:j+sz])
    print (time.time()-t0)
    h5f_A.close()
    h5f_B.close()
    h5f_C.close()

PM MAIL   Вверх
  
Ответ в темуСоздание новой темы Создание опроса
Правила форума "Алгоритмы"

maxim1000

Форум "Алгоритмы" предназначен для обсуждения вопросов, связанных только с алгоритмами и структурами данных, без привязки к конкретному языку программирования и/или программному продукту.


Если Вам понравилась атмосфера форума, заходите к нам чаще! С уважением, maxim1000.

 
1 Пользователей читают эту тему (1 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема | Алгоритмы | Следующая тема »


 




[ Время генерации скрипта: 0.0450 ]   [ Использовано запросов: 20 ]   [ GZIP включён ]


Реклама на сайте     Информационное спонсорство

 
По вопросам размещения рекламы пишите на vladimir(sobaka)vingrad.ru
Отказ от ответственности     Powered by Invision Power Board(R) 1.3 © 2003  IPS, Inc.