I need to minimize a cost function for a large number (1000s) of different inputs. Obviously, this can be implemented by looping over scipy.optimize.minimize or any other minimization routine. Here is an example:
import numpy as np
import scipy as sp
def cost(x, a, b):
return np.sum((np.sum(a * x.reshape(a.shape), axis=1) - b)**2)
a = np.random.randn(500, 40)
b = np.array(np.arange(500))
x = []
for i in range(a.shape[0]):
res = sp.optimize.minimize(cost, np.zeros(40), args=(a[None, i], b[None, i]))
x.append(res.x)
It finds x[i, :] that minimize cost for each a[i, :] and b[i], but this is very slow. I guess looping over minimize causes considerable overhead.
A partial solution is to solve for all x simultaneously:
res = sp.optimize.minimize(cost, np.zeros_like(a), args=(a, b))
This is even slower than the loop. minimize does not know that elements in x are group-wise independent. So it computes the full hessian although a block-diagonal matrix would be sufficient, considering the problem structure. This is slow and overflows my computer's memory.
Is there any way to inform minimize or another optimization function about the problem structure so that it can solve multiple indepentent optimizations in a single function call? (Similar to certain options supported by Matlab's fsolve.)