python - nD "cube" from ranges

Question

I have a mixed integer problem. I need to minimize a function, which is a weighted least square regression, the weights being dependent on the regression (iteratively reweighted least square). 7 parameters define my piecewise regression. I need to find a local minima around a first guess.

I tried to write the problem in gekko, but I somehow find it very difficult to implement. After many tries, I stopped at "negative DOF".

Anyway, I decided to brute force the problem. It works, but it's slow. I build a cube (itertools) around my working point in 7D and calculate the weighted square errors at each of the 3^7 points. I have boundaries for each dimension, and sometimes my working point is on one of the faces of my 7D domain. Technically, I have 2^p * 3^(7-p) points. I now have a list of all the values, find the minimum, move my working point there and restart building a cube, excluding all the points that I have already calculated in the previous loop steps.

Now I want to accelerate it by calculating the gradient at my working point and move faster (skip a step or two in my loop). np.gradient will require a 7d array in order to perform correctly.

Given a point, and 7 ranges around that point, how to make a 7D array in an efficient way? How to make an image of this array with the values of my function?

Please don't say 7 for loops.

question from:https://stackoverflow.com/questions/65847978/nd-cube-from-ranges

Answer 1

Regardless of whether your function is vectorized, you can use an approach with np.indices like this:

base_grid = np.indices(7 * (3,), sparse=False) - 1

This produces an array of all the combinations of -1, 0, 1 that you need. np.meshgrid does something similar, but the arrays will be separated into a tuple, which is inconvenient.

At each iteration, you modify the grid with your step (scale) and offset:

current_grid = base_grid * scale + offset

If your function is vectorized, you call it directly, the grid is 7 3x3x3x3x3x3x3 arrays. If it accepts seven inputs, just use star expansion.

If your function is not vectorized, you can still step along the corresponding elements in a single loop, not seven loops, using np.nditer:

with np.nditer([current_grid, None],
        op_axes=[list(range(1, current_grid.ndim)), None]) as it:
    for x, y in it:
        y[:] = f(*x)
    j = it.operands[1]