Приблизительная двумерная поверхность с использованием scipy.minimize L_BFGS_B

Question

Мне нужно восстановить высшее поле f(x,y) -> z из векторного поля g(x,y) -> (a,b,c), в соответствии с набором заданных уравнений, но у меня есть проблема для вычисления градиента (я использую L_BFGS_B, потому что я думаю, что это единственный оптимизатор, который позволит мне не взрывать память).

У меня следующий конфликт:

• Изображения представляют собой двумерные массивы, градиент вычисляется вдоль каждой оси и требует двух значений на пиксель
• Поле высоты имеет только 1 переменную на пиксель (высоту)

Поэтому, после вычисления градиента, Сципи жалуется, что функция не может возвращать в качестве градиента массив с (img_x * img_y, 2) формой. Только 1D фигура с точно таким же количеством, что и массив полей высоты.

Я что-то упускаю здесь очевидное?

Подробная задача

Я определяю три переменные nx, ny, nz:

nx = g(x,y).x
ny = g(x,y).y
nz = g(x,y).z

Теперь для каждой точки (x, y) сетки мне нужно вычислить высоту точки (f (x, y)) в соответствии со следующими уравнениями:

nz * [f(x+1, y) - f(x,y)] == nx
nz * [f(x, y+1) - f(x,y)] == ny

Я попытался выразить это в функции потерь:

class Eval():
    def __init__(self, g):
        '''
        g: (x, y, 3)
        '''
        self.g = g

    def loss(self, x):
        depth = x.reshape(self.g.shape[:2])

        x_roll = np.roll(depth, -1, axis=0)
        y_roll = np.roll(depth, -1, axis=1)

        dx = depth - x_roll
        dy = depth - y_roll

        nx = self.g[:,:,0]
        ny = self.g[:,:,1]
        nz = self.g[:,:,2]

        loss_x = nz * dx - nx
        loss_y = nz * dy - ny

        self.error_loss = np.stack([loss_x, loss_y], axis=-1)

        total_loss = norm(self.error_loss, axis=-1)

        return np.sum(total_loss)

def grads(self, x):
    x_roll = np.roll(self.error_loss[:,:,0], -1, axis=0)
    y_roll = np.roll(self.error_loss[:,:,1], -1, axis=1)

    dx = self.error_loss[:,:,0] - x_roll
    dy = self.error_loss[:,:,1] - y_roll

    g_xy = np.stack([dx, dy], axis=-1)

    # g_xy has shape (x, y, 2)
    # BUT THIS MUST RETURN (x * y), not (x * y, 2) nor (x * y * 2)
    # WHAT SHOULD I RETURN HERE ? ||g_xy|| ?

И назовите это так

vector_field = ... # shape: 1024, 1024, 3
x0 = np.random.uniform(size=vector_field.shape[:2])
ev = Eval(vector_field)
result = optim.fmin_l_bfgs_b(ev.loss, x0, fprime=ev.grads)

Вопросы

• Есть ли способ выразить эти уравнения для их решения по-другому? Теоретически, они представляют собой чрезмерно ограниченную линейную систему, которая может быть решена с помощью метода наименьших квадратов, но я не могу найти, как переформулировать уравнения в виде Ax = b
• В противном случае, что я должен вернуть для моего градиента? или я вычисляю это правильно?