Steepest Descent Algorithm - Maneiras de optimizar?

[JoaoRodrigues]

Caros,

Já lá vai um tempo que posto aqui algo. Andei ocupado a escrever um pequeno script que dado um número de pontos, optimiza a distribuição dos mesmo numa esfera de forma a que a distância entre pontos adjacentes seja sempre a mesma.

A ideia é primeiro gerar os pontos na esfera (uso o Golden Spiral Algorithm para isso) e depois minimizar uma função (E=1/r, onde r é a distância entre 2 pontos) até à convergência. A minimização fica a cargo de um algoritmo chamado Steepest Descent. Basicamente, a cada ciclo da minimização, são calculadas forças a actuar em cada ponto (força é a derivada da energia numa determinada direcção, x y ou z) e depois move-se o ponto de acordo com estas forças.

Isto funciona bastante bem. Para 4 pontos temos um belo tetraedro por exemplo, o que é a distribuição óptima. Para os restantes também funciona porreiro. E é rápido, para poucos pontos (digamos, N < 20 ). Não preciso de usar isto para N maior que 20 sinceramente, ou os casos são (muito) raros. Mas estou curioso em ver se consigo optimizar o código, tirando partido de mais funções de Numpy ou Scipy (sem ser a scipy.optimizer obviamente..).

Ideias? Se for preciso dou explicações mais detalhadas sobre cada passo.

O código está aqui: https://gist.github.com/4698472

Abraços

João

[ffunenga]

Deixo aqui duas possiveis optimizações que eu tentaria fazer à função calc_forces:

def calc_forces(system):

forces = [ [0,0,0] for particle in system ]
# Em vez de [0,0,0] utilizava arrays de numpy de forma a poder melhorar a
# operação que ocorre em cada iteração dos ciclos aninhados abaixo:
#	 forces = [ np.zeros(3) for particle in system ]
nparticles = len(system)

# Cautela com estes dois ciclos aninhados. Para *nparticles* elevado
# (>1000) isto é lento.
for i in xrange(nparticles):
	#(...)
	for j in xrange(i+1, nparticles):
		#(...)
		# As próximas linhas resultam da alteração feita na lista forces
		r = particle_i-particle_j
		f = numpy.linalg.norm(r)**-3
		__t = r*f
		forces[i] += __t
		forces[j] -= __t

def calc_forces(system):
n = len(system)
forces = np.array([np.zeros(3) for particle in range(n)])

idxs = np.arange(n)
idxi,idxj = np.meshgrid(idxs, idxs.T)
__triu = np.triu_indices(n, 1)
i = idxj[__triu]
j = (idxi-idxj)[__triu]

r = system[i]-system[j]
__t = np.array([r[i]*np.linalg.norm(ri)**-3 for ri in r])

#(...) aplicar na array forces as somas/subtrações nos indices i e j

return forces

[motherFFH]

Falta-me o scipy para correr o código e também não conheço o algoritmo. Tendo em conta isto, a única optimização que posso sugerir é remover algumas invariantes dos loops: o r*f como dito acima, e outra que parece importante, optimizar o steepest_descent() para remover cálculos duplicados.

Nessa função os cálculos de #3.1 e # 3.2 só variam se o valor de coordinates tiver mudado, o que só acontece se se entrar em #3.3 ou mais tarde em coordinates = old_coordinates. A estrutura do código pode ser revista para estes valores serem apenas recalculados se o valor de coordinates tiver mudado e também aproveitar os valores previamente calculados de old_coordinates. Também o "new_energy = calc_coulomb(coordinates)" logo em #3.4 devia ser feito à saída do if do #3.3 e antes ser inicializado a cur_energy.

Concretizando um pouco:

def steepest_descent(coordinates, max_steps, step_size, min_step_size, min_energy_delta, verbose=True):
#(...)
from collections import deque
## queue com duas posições apenas: old, cur. item da queue: tuple(energy, forces, force_normal)
results_cache = deque([], 2)
#(...)
# dentro do for, usar os resultados da cache

	if len(results_cache) > 0:
		t = results_cache[-1]
		cur_energy = t[0]
		forces = t[1]
		force_normal = t[2]
		results_cache.append(t) # simula old_collections = collections a nível de resultados guardados
	else:
		# 3.1 Calculate energy
		cur_energy = calc_coulomb(old_coordinates)

		# 3.2 Calculate forces
		forces = calc_forces(coordinates)
		force_normal = numpy.linalg.norm(forces)
		results_cache.append((cur_energy, forces, force_normal))
#(...)
	# 3.3 Move particles
	new_energy = cur_energy ## iniciar aqui new_energy
	if force_normal > 0:
		for p in xrange(nparticles):
			directions = (forces[p]/force_normal)*step_size
			coordinates[p] += directions
			coordinates[p] /= numpy.linalg.norm(coordinates[p])
		new_energy = calc_coulomb(coordinates) ## só calcula nova se coordinates tiverem mudado
		results_cache.clear() ## e faz reset da cache nesse caso
#(...)
	# 3.4 Check new energy and accept move if favourable + accelerate in that direction
	if new_energy >= cur_energy:
		step_size /= 2
		if step_size < min_step_size:
			print "-- Step Size Converged (%1.5f < %1.5f)" %(step_size, min_step_size)
			break
		coordinates = old_coordinates
		## se coordinates reverterem para o valor original, anula a cache referente a coordinates
		## em results_cache ficará a cache de old_coordinates, caso exista...
		if len(results_cache) > 0:
			results_cache.pop()