Python-сообщество

Допустим имеем обёртку вокруг BLAS на python (т.е. имеем numpy), можем ли руками написать ф-ии из LAPACK?(по минимуму мне надо eig, svd.)
я знаю про numpy.linalg.svd и numpy.linalg.eig так же там можно посмотреть сорцы, которые показывают в случае svd

The decomposition is performed using LAPACK routine _gesdd

implemented using the _geev LAPACK routines

References
———-
G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, FL,
Academic Press, Inc., 1980, Various pp.

mrgloom
переписать стандартные функции из numpy (svd и eig) с использованием numpy.memmap

mrgloom
и требуется написать только высокоуровневые алгоритмы.

Отредактировано doza_and (Июнь 23, 2014 19:47:27)

Я знаю что есть модификации алгоритмов для SVD и Eig, но я пока хочу для начала написать стандартные реализации этих алгоритмов чтобы сравнить их работу с numpy.linalg.svd и numpy.linalg.eig, используя подходы с numpy.memmap и cudamat.

Проблема как раз в том, что эти ф-ии numpy всего лишь обёртка вокруг lapack, вот я и интересуюсь возможно ли используя только ф-ии обёртки вокруг blas из numpy реализовать lapack (пусть это будет даже чуть медленней) на numpy.По логике это возможно,но понять это сходу у меня не получается, т.к. сорцы lapack и blas довольно таки сложные для понимания.

Смысл использовать MPI есть только если у меня будет кластер или несколько компьютеров, а учитывая, что я с этим не работал займет очень много времени.

Что мешает докупить память?

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

Отредактировано mrgloom (Июнь 24, 2014 10:20:28)

mrgloom
Еще вопрос что лучше использовать для PCA np.linalg.eig или numpy.linalg.svd?

mrgloom
вопрос скорее в правильном алгоритме

Отредактировано doza_and (Июнь 24, 2014 23:29:19)

я пытаюсь написать свою версию метода PCA

для теста можно использовать MNIST
но надо чтобы это работало на mnist8m
на текущий момент я тестирую на матрице ~1kk x 1k .

я так понимаю разговор о данных имеет смысл только в том случае, если известно что они разрежены, т.е. использовать методы которые принимают sparse matrix.

Для больших симметричных матриц обычно используют метод Арнольди.

Отредактировано mrgloom (Июнь 25, 2014 09:52:06)