Я знаю что есть модификации алгоритмов для SVD и Eig, но я пока хочу для начала написать стандартные реализации этих алгоритмов чтобы сравнить их работу с numpy.linalg.svd и numpy.linalg.eig, используя подходы с numpy.memmap и cudamat.
Проблема как раз в том, что эти ф-ии numpy всего лишь обёртка вокруг lapack, вот я и интересуюсь возможно ли используя только ф-ии обёртки вокруг blas из numpy реализовать lapack (пусть это будет даже чуть медленней) на numpy.По логике это возможно,но понять это сходу у меня не получается, т.к. сорцы lapack и blas довольно таки сложные для понимания.
Смысл использовать MPI есть только если у меня будет кластер или несколько компьютеров, а учитывая, что я с этим не работал займет очень много времени.
Что мешает докупить память?
ну сейчас у меня 8Гб, максимум можно поставить 16Гб, если менять мат. плату, то за приемлемые деньги можно поставить 32-64Гб я так думаю, но модификация железа это экстенсивный путь, вопрос скорее в правильном алгоритме который умеет работать с иерархией памяти hdd<->sdd<->ram<->cache т.е. по сути мы имеем многоуровневое кэширование в более медленной памяти
тут кстати есть еще интересное сравнение, что можно передавать данные быстрее сжав их или же надо алгоритм который работает с данными кусками, т.е. занят только определенный объем памяти, это так же полезно и для вычислений на gpu, т.к. там явно память ограниченна.
Например есть блочное перемножение матриц, не надо загружать все матрицы в память и по идее если мы правильно выберем размер блока, то может быть даже быстрее, чем если бы мы делали перемножение матриц используя memmaped array.
p.s. говорят что np.dot работает медленнее в некоторых случаях из-за копирования хотя по сути это обёртка над blas dgemm как я понимаю и это зависит от того как хранятся матрицы row major или column major т.е. С order или Fortran order.
вот тест np.dot
http://stackoverflow.com/questions/24016207/how-order-of-numpy-array-influence-on-multiplication-speed а вот например реализация PCA с помощью memmap, но работать будет только если кол-во сэмплов большое, а кол-во переменных достаточно мало чтобы np.linalg.eig не упало (как раз это меня и побудило переписать np.linalg.eig). Пришлось еще переписать np.cov .
Смысл всего этого только в том, чтобы алгоритм не упал при нехватке памяти и для теста, а так на практике нахождение матрицы ковариации и потом проекция при больших размерах массива требует много времени и нужен другой алгоритм.
def cov_mat(fmat):
#if fmat centered then 2*(A.T*A)/(n-1) covariance matrix
M= fmat.shape[0]
N= fmat.shape[1]
fcov= np.memmap('cov.npy', dtype='float32', mode='w+', shape=(N,N))
#assinging don't copy?
# fmat_tr= np.memmap('A_tr.npy', dtype='float32', mode='w+', shape=(N,M))
fmat_tr= fmat.T
np.dot(fmat_tr,fmat,out=fcov)
fcov*= 2/(N-1)
return fcov
#works ok x64
#np.linalg.eig ok only if number of dimensions is small
#assume M >> N
#fdata is memmaped array M x N
def pca_mm(fdata,k):
M= fdata.shape[0]
N= fdata.shape[1]
print fdata.shape
#get mean
mean= np.mean(fdata,axis=0)
# print mean.shape
# print mean
#M x N
fdata_c= np.memmap('data_c.npy', dtype='float32', mode='w+', shape=(M,N))
# fdata_c= fdata-mean
for i in xrange(M):
fdata_c[i,:]= fdata[i,:]-mean
# print fdata_c.shape
# print data_c
#N x N
#calculate covariance matrix
#t0= time.time()
fcov= cov_mat(fdata_c)
# print (time.time()-t0)
# covData=np.cov(fdata_c,rowvar=0)
# print fcov
# print covData
# print fcov.shape
#we store eval,evec in memory
#t0= time.time()
eigenvalues, eigenvectors = np.linalg.eig(fcov)
# print (time.time()-t0)
# print eigenvalues.shape # N long
# print eigenvectors.shape # N x N
#print eigenvalues
#print eigenvectors
#sort and get k largest eigenvalues
idx = eigenvalues.argsort()[-k:][::-1]
# print idx
eigenvalues = eigenvalues[idx] # k long
eigenvectors = eigenvectors[:,idx] # N x k
# print eigenvalues.shape
# print eigenvectors.shape
# print eigenvalues
# print eigenvectors
#projection
fpr= np.memmap('pr.npy', dtype='float32', mode='w+', shape=(M,k))
np.dot(fdata_c,eigenvectors,out=fpr) # (M N) * (N k) = (M k)
# print fpr.shape
#reconstruction
frec= np.memmap('rec.npy', dtype='float32', mode='w+', shape=(M,N))
np.dot(fpr, eigenvectors.T,out=frec) #(M k) * (N k).T = (M N)
# print frec.shape
#reconstruction error
print np.allclose(fdata_c,frec) # true if matrices are equal
Еще вопрос что лучше использовать для PCA np.linalg.eig или numpy.linalg.svd? говорят, что numpy.linalg.svd более численно стабилен, а что по памяти?