Использование scipy.fsolve на матрице - PullRequest
Новая
       95

Использование scipy.fsolve на матрице

0 голосов
/ 20 января 2019

После помощи, которую мне дали здесь Я пытался реализовать ее в своем скрипте, но мне не удалось запустить ее умно.

Мне нужно использовать этот алгоритм для каждого пикселя изображения 4072x3080, это занимает около 1:30 для всего процесса, поэтому я попытался как-то форсировать его, но я получаю эту ошибку: 

ValueError                                Traceback (most recent call last)
<ipython-input-12-99c1f41dbba7> in <module>()
----> 1 res = scipy.optimize.fsolve(func, x0=np.ones((K.shape[0], K.shape[1])), args=(f[:,None], g[:,None], K))

/*/python2.7/site-packages/scipy/optimize/minpack.pyc in fsolve(func, x0, args, fprime, full_output, col_deriv, xtol, maxfev, band, epsfcn, factor, diag)
    146                'diag': diag}
    147 
--> 148     res = _root_hybr(func, x0, args, jac=fprime, **options)
    149     if full_output:
    150         x = res['x']

/*/python2.7/site-packages/scipy/optimize/minpack.pyc in _root_hybr(func, x0, args, jac, col_deriv, xtol, maxfev, band, eps, factor, diag, **unknown_options)
    212     if not isinstance(args, tuple):
    213         args = (args,)
--> 214     shape, dtype = _check_func('fsolve', 'func', func, x0, args, n, (n,))
    215     if epsfcn is None:
    216         epsfcn = finfo(dtype).eps

/*/python2.7/site-packages/scipy/optimize/minpack.pyc in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
     25 def _check_func(checker, argname, thefunc, x0, args, numinputs,
     26                 output_shape=None):
---> 27     res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
     28     if (output_shape is not None) and (shape(res) != output_shape):
     29         if (output_shape[0] != 1):

<ipython-input-7-911c817cb57d> in func(x, f, g, K)
      1 def func(x, f, g, K):
----> 2     return np.sum(f * np.exp(-g*x), axis=0) - K
      3 
      4 
      5 def derivative(x, f, g, K):

ValueError: operands could not be broadcast together with shapes (13551616,) (4072,3328)

Это код, который я пробовал: 

def func(x, f, g, K):
    return np.sum(f * np.exp(-g*x), axis=0) - K


def derivative(x, f, g, K):
    return np.sum(-g*f * np.exp(-g*x), axis=0)

+ 

res = scipy.optimize.fsolve(func, x0=np.ones((K.shape[0], K.shape[1])), args=(f[:,None], g[:,None], K))

f и g являются массивами (47,), где K является (4072, 3328) изображением

В остальном медленный процесс идет таким образом: (но этот процесс все равно не работает с производной. 

res = np.ones((mbn.shape[0],mbn.shape[1]))
for i_x in range(0,mbn.shape[0]):
    if i_x%10 == 0:
        print i_x/4070 
    for i_y in range(0,mbn.shape[1]):
        res[i_x,i_y] = scipy.optimize.fsolve(func, x0=1, args=(f[:], g[:], K[i_x,i_y]) )

Это будет ошибкой, если я попытаюсь использовать медленный метод с производным 

    ---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
ValueError: object of too small depth for desired array

---------------------------------------------------------------------------
error                                     Traceback (most recent call last)
<ipython-input-8-3587dcccfd93> in <module>()
      4         print i_x/4070
      5     for i_y in range(0,mbn.shape[1]):
----> 6         res[i_x,i_y] = scipy.optimize.fsolve(func, fprime=derivative, x0=1, args=(f[:], g[:], K[i_x,i_y]) )

/*/python2.7/site-packages/scipy/optimize/minpack.pyc in fsolve(func, x0, args, fprime, full_output, col_deriv, xtol, maxfev, band, epsfcn, factor, diag)
    146                'diag': diag}
    147 
--> 148     res = _root_hybr(func, x0, args, jac=fprime, **options)
    149     if full_output:
    150         x = res['x']

/*/python2.7/site-packages/scipy/optimize/minpack.pyc in _root_hybr(func, x0, args, jac, col_deriv, xtol, maxfev, band, eps, factor, diag, **unknown_options)
    232         with _MINPACK_LOCK:
    233             retval = _minpack._hybrj(func, Dfun, x0, args, 1,
--> 234                                      col_deriv, xtol, maxfev, factor, diag)
    235 
    236     x, status = retval[0], retval[-1]

error: Result from function call is not a proper array of floats.

Ответы [ 3 ]

0 голосов
/ 21 января 2019

Обычно плохая идея преобразовать решение n скалярных задач поиска корней в векторизованную n -мерную задачу поиска корней.Особые решения, которые принимает скалярный решатель, обычно не являются однородными по всем входным вариантам, и любой общий решатель систем будет пытаться вычислить или приблизить якобиан.Несмотря на то, что мы знаем, что это диагональ этого параметра, может оказаться сложным сообщить об этом решателю.Кроме того, как уже было сказано, решения о размерах шагов, поисках строк и т. Д., Когда они приняты равномерно для векторизованной задачи, почти гарантированно будут неоптимальными для каждой отдельной скалярной задачи.

Наилучшая стратегия заключается в уменьшении количестваскалярных задач для решения или, по крайней мере, использования одномерной природы, чтобы получить хорошее приближение обратной функции, которая требует только очень мало итераций для полировки решения.

случай небольшого количества целочисленных значенийв K

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

Первая идея для более быстрого вычисления состоит в том, чтобы сначала определить диапазон значений K, вычислить значение x для каждого из фактически существующих значений, а затем распределить эти значения обратно по массиву изображений.,Например, если K принимает только целочисленные значения от 0 до 255, то создание массива x_vals, где func(x_vals[k],f,g,0)=k позволяет получить результирующий массив как x_vals[K] a la 

x_vals = np.array([0.1,0.2,0.3,0.4])
K = np.array([[1,2],[0,2],[3,3]]);
x=x_vals[K]; print x

array([[ 0.2,  0.3],
       [ 0.1,  0.3],
       [ 0.4,  0.4]])

случай относительно небольшого диапазона нецелых значений в K

Если K содержит большее количество (нецелых) значений в некоторых относительно (дозначения в f) небольшом диапазоне, это все еще может обеспечить большое улучшение для вычисления решений для некоторой выборки K_samples = np.linspace(np.amin(K), np.amax(K), 256), так что func (x_samples [k], f, g, 0) = K_samples [k]) `

Затем решение или, по крайней мере, очень хорошее приближение для дальнейшего итеративного уточнения получается с использованием интерполяции 

x = np.interp(K, K_samples, x_samples)

с использованием прямой оценки функции

Если векторы f и g имеют положительные значения, то сумма экспонент представляет собой монотонно падающую функцию, обратное значение которой можно получить с помощью простой справочной таблицы и ее линейной интерполяции.Для значений функции суммы экспонент можно получить 

sum(f)*exp(-max(g)*x) <= K <= sum(f)*exp(-min(g)*x)

, чтобы можно было вычислить диапазон для x как 

- log(max(K)/sum(f)) / max(g)  <=  x  <=  - log(min(K)/sum(f)) / min(g)

Используя эти верхнюю и нижнюю границы длясоздать массив x_samples, получая соответствующий K_samples путем прямой оценки func.Линейная интерполяция для получения аппроксимации обратной функции выполняется, как описано выше.

0 голосов
/ 21 января 2019

Вы можете векторизовать проблему, создав функцию, которая принимает N входов и N выходов, где N - количество пикселей.Это включает в себя сглаживание входного изображения и обработку его как одномерного массива.В этой настройке входы независимы от выходов и, следовательно, якобиан диагональн.Поскольку fsolve вычисляет полную аппроксимацию якобиана, вам в конечном итоге не хватит памяти (MemoryError).Вместо этого вы можете использовать scipy.optimize.root с method='diagbroyden', который использует аппроксимацию, отслеживая только диагональный якобиан: 

import numpy as np
import scipy.optimize as optimize

def func(x, f, g, K):
    return np.sum(f * np.exp(-g*x[:, None]), axis=1) - K

np.random.seed(123)

f = np.random.uniform(size=(47,))
g = np.random.uniform(size=f.shape)
img = np.random.uniform(size=(4072, 3328)).ravel()
K = func(img, f, g, 0)

res = optimize.root(func, method='diagbroyden', x0=0.5*np.ones(img.size), args=(f, g, K))
print('Success:', res.success)
print('Message:', res.message)
assert np.allclose(img, res.x)

С этим методом, однако вы не можете использовать преимущества аналитической производной, котораяможет быть рассчитан для вашей конкретной функции.

0 голосов
/ 20 января 2019

func в optimize.fsolve может принимать одномерный вектор, но не двумерные массивы. Таким образом, хотя K и x являются двумерными, для этого расчета нам необходимо преобразовать их в одномерные массивы. 

Kshape = K.shape
K = K.ravel()

Затем, после вызова optimize.fsolve, вы можете изменить результат, чтобы он снова был 2D: 

res = optimize.fsolve(func, x0=np.ones(K.shape).ravel(), args=(f, g, K))
res = res.reshape(Kshape)

Затем вы можете избежать двойных циклов for, записав вычисление следующим образом: 

import numpy as np
import scipy.optimize as optimize

np.random.seed(123)

def func(x, f, g, K):
    return np.sum(f * np.exp(-g*x[:, None]), axis=-1) - K


def derivative(x, f, g, K):
    return np.sum(-g*f * np.exp(-g*x[:, None]), axis=-1)


f = np.random.uniform(size=(47,))
g = np.random.uniform(size=f.shape)
K = np.random.uniform(size=(4072,3080))
Kshape = K.shape
K = K.ravel()

res = optimize.fsolve(func, x0=np.ones(K.shape).ravel(), args=(f, g, K))
res = res.reshape(Kshape)
print(res)

Обратите внимание, что g*x[:, None] использует широковещание для генерации двумерного массива формы (4072*3080, 47). 2D массив f * np.exp(-g*x[:, None]), которая также имеет форму (4072*3080, 47), затем суммируется по последней оси (то есть axis=-1). Это оставляет одномерный массив формы (4072*3080,). fsolve решает для x и возвращает одномерный массив формы (4072*3080,).

res = res.reshape(Kshape) меняет решение, чтобы иметь форму (4072, 3080).

...