SciPy, оптимизация

  • Tutorial

SciPy (произносится как сай пай) — это пакет прикладных математических процедур, основанный на расширении Numpy Python. С SciPy интерактивный сеанс Python превращается в такую же полноценную среду обработки данных и прототипирования сложных систем, как MATLAB, IDL, Octave, R-Lab и SciLab. Сегодня я хочу коротко рассказать о том, как следует применять некоторые известные алгоритмы оптимизации в пакете scipy.optimize. Более подробную и актуальную справку по применению функций всегда можно получить с помощью команды help() или с помощью Shift+Tab.


Введение


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


Итак, модуль scipy.optimize включает в себя реализацию следующих процедур:


  1. Условной и безусловной минимизации скалярных функций нескольких переменных (minim) с помощью различных алгоритмов (симплекс Нелдера-Мида, BFGS, сопряженных градиентов Ньютона, COBYLA и SLSQP)
  2. Глобальной оптимизации (например: basinhopping, diff_evolution)
  3. Минимизация остатков МНК (least_squares) и алгоритмы подгонки кривых нелинейным МНК (curve_fit)
  4. Минимизации скалярной функций одной переменной (minim_scalar) и поиска корней (root_scalar)
  5. Многомерные решатели системы уравнений (root) с использованием различных алгоритмов (гибридный Пауэлла, Левенберг-Марквардт или крупномасштабные методы, такие как Ньютона-Крылова).

В этой статье мы рассмотрим только первый пункт из всего этого списка.


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


Функция minim из пакета scipy.optimize предоставляет общий интерфейс для решения задач условной и безусловной минимизации скалярных функций нескольких переменных. Чтобы продемонстрировать ее работу, нам понадобится подходящая функция нескольких переменных, которую мы будем по-разному минимизировать.


Для этих целей прекрасно подойдет функция Розенброка от N переменных, которая имеет вид:


$f\left(\mathbf{x}\right)=\sum_{i=1}^{N-1}[100\left(x_{i+1}-x_{i}^{2}\right)^{2}+\left(1-x_{i}\right)^{2}]$


Несмотря на то, что функция Розенброка и ее матрицы Якоби и Гессе (первой и второй производной соответственно) уже определены в пакете scipy.optimize, определим ее самостоятельно.


import numpy as np

def rosen(x):
    """The Rosenbrock function"""
    return np.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0, axis=0)

Для наглядности отрисуем в 3D значения функции Розенброка от двух переменных.


Код для отрисовки
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter

# Настраиваем 3D график
fig = plt.figure(figsize=[15, 10])
ax = fig.gca(projection='3d')

# Задаем угол обзора
ax.view_init(45, 30)

# Создаем данные для графика
X = np.arange(-2, 2, 0.1)
Y = np.arange(-1, 3, 0.1)
X, Y = np.meshgrid(X, Y)
Z = rosen(np.array([X,Y]))

# Рисуем поверхность
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()


Зная заранее, что минимум равен 0 при $x_i = 1$, рассмотрим примеры того, как определить минимальное значение функции Розенброка с помощью различных процедур scipy.optimize.


Симплекс-метод Нелдера-Мида (Nelder-Mead)


Пусть имеется начальная точка x0 в 5-мерном пространстве. Найдем ближайшую к ней точку минимума функции Розенброка с помощью алгоритма симплекса Nelder-Mead (алгоритм указан в качестве значения параметра method):


from scipy.optimize import minimize
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='nelder-mead',
    options={'xtol': 1e-8, 'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 339
         Function evaluations: 571
[1. 1. 1. 1. 1.]

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


Метод Пауэлла


Другим алгоритмом оптимизации, в котором вычисляются только значения функций, является метод Пауэлла. Чтобы использовать его, нужно установить method = 'powell' в функции minim.


x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='powell',
    options={'xtol': 1e-8, 'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 19
         Function evaluations: 1622
[1. 1. 1. 1. 1.]

Алгоритм Бройдена-Флетчера-Голдфарба-Шанно (BFGS)


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


Найдем производную от функции Розенброка в аналитическом виде:


$ \frac{\partial f }{\partial x_j} = \sum\limits_{i=1}^N 200(x_i - x_{i-1}^2)(\delta_{i,j} - 2x_{i-1, j}) - 2(1 - x_{i-1}) \delta_{i-1,j} = $


$ = 200(x_j - x_{j-1}^2) - 400x_j(x_{j+1} - x_j^2) - 2(1-x_j) $


Это выражение справедливо для производных всех переменных, кроме первой и последней, которые определяются как:


$ \frac{\partial f }{\partial x_0} = -400 x_0(x_1 - x_0^2) - 2(1 - x_0), $


$ \frac{\partial f }{\partial x_{N-1}} = 200(x_{N-1} - x_{N-2}^2). $


Посмотрим на функцию Python, которая вычисляет этот градиент:


def rosen_der (x):
    xm = x [1: -1]
    xm_m1 = x [: - 2]
    xm_p1 = x [2:]
    der = np.zeros_like (x)
    der [1: -1] = 200 * (xm-xm_m1 ** 2) - 400 * (xm_p1 - xm ** 2) * xm - 2 * (1-xm)
    der [0] = -400 * x [0] * (x [1] -x [0] ** 2) - 2 * (1-x [0])
    der [-1] = 200 * (x [-1] -x [-2] ** 2)
    return der

Функция вычисления градиента указывается в качестве значения параметра jac функции minim, как показано ниже.


res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 25
         Function evaluations: 30
         Gradient evaluations: 30
[1.00000004 1.0000001  1.00000021 1.00000044 1.00000092]

Алгоритм сопряженных градиентов (Ньютона)


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


$ f\left(\mathbf{x}\right)\approx f\left(\mathbf{x}_{0}\right)+\nabla f\left(\mathbf{x}_{0}\right)\cdot\left(\mathbf{x}-\mathbf{x}_{0}\right)+\frac{1}{2}\left(\mathbf{x}-\mathbf{x}_{0}\right)^{T}\mathbf{H}\left(\mathbf{x}_{0}\right)\left(\mathbf{x}-\mathbf{x}_{0}\right)$


где $\mathbf{H}\left(\mathbf{x}_{0}\right)$ является матрицей вторых производных (матрица Гессе, гессиан).
Если гессиан положительно определен, то локальный минимум этой функции можно найти, приравняв нулевой градиент квадратичной формы к нулю. В результате получится выражение:


$ \mathbf{x}_{\textrm{opt}}=\mathbf{x}_{0}-\mathbf{H}^{-1}\nabla f $


Обратный гессиан вычисляется с помощью метода сопряженных градиентов. Пример использования этого метода для минимизации функции Розенброка приведен ниже. Чтобы использовать метод Newton-CG, необходимо задать функцию, которая вычисляет гессиан.
Гессиан функции Розенброка в аналитическом виде равен:


$ H_{i,j} = \frac{\partial^2 f }{\partial x_i x_j} = 200(\delta_{i,j} - 2x_{i-1} \delta{i-1,j} - 400x_i(\delta_{i+1,j} - 2x_i \delta{i,j}) - 400 \delta_{i,j} (x_{i+1} - x_i^2) + 2\delta_{i,j} = $


$ = (202 + 1200x_i^2 - 400x_{i+1})\delta_{i,j} - 400x_i \delta_{i+1,j} - 400x_{i-1}\delta_{i-1,j} $


где $i, j \in \left[ 1, N-2 \right]$ и $i, j \in \left [0, N-1 \right]$, определяют матрицу $ N \times N $.


Остальные ненулевые элементы матрицы равны:


$ \frac{\partial^2 f }{\partial x_0^2} = 1200x_0^2 - 400x_1 +2 $


$ \frac{\partial^2 f }{\partial x_0 x_1} = \frac{\partial^2 f }{\partial x_1 x_0} = -400x_0 $


$ \frac{\partial^2 f }{\partial x_{N-1} x_{N-2}} = \frac{\partial^2 f }{\partial x_{N-2} x_{N-1}} = -400x_{N-2} $


$ \frac{\partial^2 f }{\partial x_{N-1}^2} = 200x $


Например, в пятимерном пространстве N = 5, матрица Гессе для функции Розенброка имеет ленточный вид:


$ \tiny \mathbf{H}=\begin{bmatrix} 1200x_{0}^{2}-400x_{1}+2 & -400x_{0} & 0 & 0 & 0\\ -400x_{0} & 202+1200x_{1}^{2}-400x_{2} & -400x_{1} & 0 & 0\\ 0 & -400x_{1} & 202+1200x_{2}^{2}-400x_{3} & -400x_{2} & 0\\ 0 & & -400x_{2} & 202+1200x_{3}^{2}-400x_{4} & -400x_{3}\\ 0 & 0 & 0 & -400x_{3} & 200\end{bmatrix} $


Код, который вычисляет этот гессиан вместе с кодом для минимизации функции Розенброка с помощью метода сопряженных градиентов (Ньютона):


def rosen_hess(x):
    x = np.asarray(x)
    H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)
    diagonal = np.zeros_like(x)
    diagonal[0] = 1200*x[0]**2-400*x[1]+2
    diagonal[-1] = 200
    diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
    H = H + np.diag(diagonal)
    return H

res = minimize(rosen, x0, method='Newton-CG', 
               jac=rosen_der, hess=rosen_hess,
               options={'xtol': 1e-8, 'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 24
         Function evaluations: 33
         Gradient evaluations: 56
         Hessian evaluations: 24
[1.         1.         1.         0.99999999 0.99999999]

Пример с определением функции произведения гессиана и произвольного вектора


В реальных задачах вычисление и хранение всей матрицы Гессе может потребовать значительных ресурсов времени и памяти. При этом фактически нет необходимости задавать саму матрицу Гессе, т.к. для процедуры минимизации нужен только вектор, равный произведению гессиана с другим произвольным вектором. Таким образом, с вычислительной точки зрения намного предпочтительней сразу определить функцию, которая возвращает результат произведения гессиана с произвольным вектором.


Рассмотрим функцию hess, которая принимает вектор минимизации в качестве первого аргумента, а произвольный вектор — как второй аргумент (наряду с другими аргументами минимизируемой функции). В нашем случае вычислить произведение гессиана функции Розенброка с произвольным вектором не очень сложно. Если p — произвольный вектор, то произведение $H(x) \cdot p$ имеет вид:


$ \mathbf{H}\left(\mathbf{x}\right)\mathbf{p}=\begin{bmatrix} \left(1200x_{0}^{2}-400x_{1}+2\right)p_{0}-400x_{0}p_{1}\\ \vdots\\ -400x_{i-1}p_{i-1}+\left(202+1200x_{i}^{2}-400x_{i+1}\right)p_{i}-400x_{i}p_{i+1}\\ \vdots\\ -400x_{N-2}p_{N-2}+200p_{N-1}\end{bmatrix}. $


Функция, вычисляющая произведение гессиана и произвольного вектора, передается как значение аргумента hessp функции minimize:


def rosen_hess_p(x, p):
    x = np.asarray(x)
    Hp = np.zeros_like(x)
    Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1]
    Hp[1:-1] = -400*x[:-2]*p[:-2]+(202+1200*x[1:-1]**2-400*x[2:])*p[1:-1] \
    -400*x[1:-1]*p[2:]
    Hp[-1] = -400*x[-2]*p[-2] + 200*p[-1]
    return Hp

res = minimize(rosen, x0, method='Newton-CG',
               jac=rosen_der, hessp=rosen_hess_p,
               options={'xtol': 1e-8, 'disp': True})

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 24
         Function evaluations: 33
         Gradient evaluations: 56
         Hessian evaluations: 66

Алгоритм доверительной области (trust region) сопряженных градиентов (Ньютона)


Плохая обусловленность матрицы Гессе и неверные направления поиска могут привести к тому, что алгоритм сопряженных градиентов Ньютона может быть неэффективным. В таких случаях предпочтение отдается методу доверительной области (trust-region) сопряженных градиентов Ньютона.


Пример с определением матрицы Гессе:


res = minimize(rosen, x0, method='trust-ncg',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 20
         Function evaluations: 21
         Gradient evaluations: 20
         Hessian evaluations: 19
[1. 1. 1. 1. 1.]

Пример с функцией произведения гессиана и произвольного вектора:


res = minimize(rosen, x0, method='trust-ncg', 
                jac=rosen_der, hessp=rosen_hess_p, 
                options={'gtol': 1e-8, 'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 20
         Function evaluations: 21
         Gradient evaluations: 20
         Hessian evaluations: 0
[1. 1. 1. 1. 1.]

Методы Крыловского типа


Подобно методу trust-ncg, методы Крыловского типа хорошо подходят для решения крупномасштабных задач, поскольку в них используется только матрично-векторные произведения. Их суть в решении задачи в доверительной области, ограниченной усеченным подпространством Крылова. Для неопределенных задач лучше использовать этот метод, так как он использует меньшее количество нелинейных итераций за счет меньшего количества матрично-векторных произведений на одну подзадачу, по сравнению с методом trust-ncg. Кроме того, решение квадратичной подзадачи находится более точно, чем методом trust-ncg.
Пример с определением матрицы Гессе:


res = minimize(rosen, x0, method='trust-krylov',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 19
         Function evaluations: 20
         Gradient evaluations: 20
         Hessian evaluations: 18

print(res.x)

    [1. 1. 1. 1. 1.]

Пример с функцией произведения гессиана и произвольного вектора:


res = minimize(rosen, x0, method='trust-krylov',
               jac=rosen_der, hessp=rosen_hess_p,
               options={'gtol': 1e-8, 'disp': True})

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 19
         Function evaluations: 20
         Gradient evaluations: 20
         Hessian evaluations: 0

print(res.x)

    [1. 1. 1. 1. 1.]

Алгоритм приближенного решения в доверительной области


Все методы (Newton-CG, trust-ncg и trust-krylov) хорошо подходят для решения крупномасштабных задач (с тысячами переменных). Это связано с тем, что лежащий в их основе алгоритм сопряженных градиентов подразумевает приближенное нахождение обратной матрицы Гессе. Решение находится итеративно, без явного разложения гессиана. Поскольку требуется определить только функцию для произведение гессиана и произвольного вектора, этот алгоритм особенно хорош для работы с разреженными (ленточными диагональными) матрицами. Это обеспечивает низкие затраты памяти и значительную экономию времени.


В задачах среднего размера затраты на хранение и факторизацию гессиана не имеют решающего значения. Это значит, что можно получить решение за меньшее количество итераций, разрешив подзадачи области доверия почти точно. Для этого некоторые нелинейные уравнения решаются итеративно для каждой квадратичной подзадачи. Такое решение требует обычно 3 или 4 разложения Холецкого матрицы Гессе. В результате метод сходится за меньшее количество итераций и требует меньше вычислений целевой функции, чем другие реализованные методы доверительной области. Этот алгоритм подразумевает только определение полной матрицы Гессе и не поддерживает возможность использовать функцию произведения гессиана и произвольного вектора.


Пример с минимизацией функции Розенброка:


res = minimize(rosen, x0, method='trust-exact',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})
res.x

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 13
         Function evaluations: 14
         Gradient evaluations: 13
         Hessian evaluations: 14

array([1., 1., 1., 1., 1.])

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


Источник: https://docs.scipy.org/doc/scipy/reference/

Поделиться публикацией
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее
Реклама

Комментарии 19

    +1
    Спасибо за разбор! А есть ли какой-то более подробный критерий/гайд по выбору алгоритма оптимизации для конкретной задачи? Скажем, в зависимости от крутизны градиентов, локальных минимумов, периодичности и тп?
    PS
    у вас Гессиан в методе Ньютона не отображается правильно
      +1

      Посмотрите здесь.

        +1
        Есть опыт, что более важным фактором является качество и безглючность реализации.

        Лично я многие годы пользуюсь оптимизациями из библиотеки NAG для Fortran. Есть вариант для C, в частности, этим вариантом Maple пользуется. Говорят, они и для Python выпустили пакет, но лично я его не использовал.
          0
          О, круто, спасибо! Правда жалко, что проприетарно и платно…
            0
            Таки да, NAG не очень подходит для юношей с горящим взором, которые, прямо таки сразу, желают внести свою лепту в теорию и/или практику вычислительной математики. Тем да, проприетарность и невозможность совместного допиливания ломом через колено, как бы не очень.

            Но для остепенившихся товарищей, которые уже разобрались, что у них в математике может получиться, чего, увы, уже никогда. Или для тех, кому нужно реальные прикладные задачи решать, надёжно и в прогнозируемые сроки. NAG он для таких.
              0
              Ну да, я так и понял:) Тут скорее вопрос, что при покупке для института, например, должно быть хорошее обоснование, почему их, а не тот же scipy.optimize. Но тут вы правы, с опытом приходит и понимание, наверное, так что подожду, пусть полежит в закладках.
                0
                Как то уж слишком однобоко у Вас получается. Существует множество организаций, для которых покупка коммерческих математических программ экономически просто нецелесообразна. А практические задачки решать тоже надо, хотя бы и учебные.
                Добавим сюда основные преимущества python — интерактивность и простой синтаксис, и получим отличный бесплатный швейцарский нож. Он не заменит пилораму, но для подручной работы намного удобней.
                  +1
                  Капитализм, конечно, загнивает, и коммунизм, несомненно, с неизбежностью наступит, и OpenSource часть его (на полном серьёзе). Но надо понимать, что SciPy версии 1.0 появился всего пару лет назад. Для сравнения Linux 1.0 появился 25 лет назад и только последнее время Linux смог претендовать на нишу более менее надёжного швейцарского ножа.

                  «экономически просто нецелесообразна», т.е. машины на поездить или для посчитать целесообразны, а программное обеспечение нет. Странно немного.

                  Нет, ну если учить вычислительным методам, то можно и студентов знакомить с потрохами и процессом разработки той же SciPy/NumPy. Я, конечно, ретроград, и предпочел бы их учить на примере исходных кодов NAG для Fortran, или, на крайний случай, NAG для C. Ибо откуда они ещё узнают, как это должно нормально и надёжно работать?!

                  «интерактивность и простой синтаксис» — а как бы у MATLAB/Maple/IRAF интерактивность повыше будет. Не?

                  Учитывая основной недостаток Питона — нелюбовь к многоядерным/многопроцессорным архитектурам, всё таки, увы, Питон не Ява. И SciPy/NumPy, и NAG for Python, осваивают его медленно используя нативный код. По крайней мере, астрономические приложения на Питоне, типа PHOEBE (да, да, с использованием SciPy/NumPy) требуют существенного энтузиазма и всё равно:
                  PHOEBE 1.0 (legacy) should be used for reliable trustable science results and for cases that do not require the precision or additional physics introduced by PHOEBE 2. PHOEBE 1.0 (legacy) is still significantly faster than PHOEBE 2.

                  Примечание: PHOEBE 1.0 (legacy) — Питон обёртка FORTRAN кода образца (98..2007), стабильная версия Питон обёртки 2008 года. PHOEBE 2 — продукт почти 15-летних танцев с Питон, потом плюс NymPy, сейчас ещё и SciPy. А доверия всё нет, и нет, ещё и тормозит гад не по детски.

                  Ниша работы с Питоном, мне так кажется, перспективные сетевые вычисления в облаках. Т.е. нормальная жизнь — через 10..15 лет, а пока основной бонус — участие в процессе освоения. Впрочем, могу ошибаться.
                    0
                    «экономически просто нецелесообразна», т.е. машины на поездить или для посчитать целесообразны, а программное обеспечение нет. Странно немного.

                    Любой институт. Если научная группа не занимается постоянной и сложной оптимизацией, средств scipy обычно хватает с головой, особенно для студентов.
                    Если нет задачи научиться численным методам, а нужно обработать данные или построить простую модель в ограниченный срок (скажем, полгода-год написания магистерской работы у студентов) — питон в целом идеальный вариант.

                    «интерактивность и простой синтаксис» — а как бы у MATLAB/Maple/IRAF интерактивность повыше будет. Не?

                    Jupyter/ipython + любая IDE делают интерактивность питона вполне на уровне Матлаба.

                    Но в целом насчет скорости разработки вы правы.
                      0
                      … в ограниченный срок (скажем, полгода-год написания магистерской работы у студентов) ...
                      Вот в этом месте может скрываться засада. Все ж под Богом глюков в программном обеспечении ходим.
                      … сложной оптимизацией, средств scipy обычно хватает с головой...
                      Что есть сложная? Чем дольше живу, тем более вычислительно сложные модели (целевые функции) в т.ч. и у студентов/аспирантов. А расчёт вычислительной модели на Pythone/SciPy/NumPy, конечно как повезёт, но бывает и на порядок-два тормознее, чем на C/FORTRAN/MATLAB/Maple.

                      Плюс, на мой вкус, недоделанный и, главное, не отлаженный код самоих методов оптимизаций и прочих обвязок вплоть до matplotlib. Правда, я ж сам с SciPy имею дело только когда у товарищей и/или товарок проблемы, возможно, умные люди умеют как-то приспосабливаться? Кто знает?
                        0
                        А расчёт вычислительной модели на Pythone/SciPy/NumPy, конечно как повезёт, но бывает и на порядок-два тормознее, чем на C/FORTRAN/MATLAB/Maple.

                        Сложно сказать, хорошо написанный код на numpy в anaconda (т.е., скажем, с numba) субъективно работает не медленнее Матлаба. C и Fortran сложно, конечно, побить, но зависит от задач, иногда питон и побыстрее будет. Вот можно бенчмарки посмотреть, интереса ради.

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

                        Насчет кода методов не знаю, тут я совсем не специалист, да и сравнить не с чем. А про обвязки, кажется, все как раз очень неплохо обстоит сейчас: pandas для работы с данными, seaborn для визуализации, интерактивные графики с jupyther в matplotlib. В итоге воркфлоу для обработки данных получается быстрым, удобным и читабельным на выходе: получающийся в итоге ноутбук ipython на мой вкус вообще лучший вариант для научных задач, особенно если нужно делиться с кем-то.
                          +1
                          Ну и кстати, питон удобен разнообразием: у нас лабах цифровые контроллеры экспериментов работают на питоне, туда же встроен базовый анализ + документация, с выводом этого через ноутбуки jupyter в лабораторные книги/логи эксперимента. На питоне же написаны стандартные скрипты обработки данных и построения графиков/визуализации. В итоге новому студенту не нужно учить несколько разных языков, особенно специфичных, как фортран, а достаточно знать основы питона, и можно проводить полный цикл эксперимента, от собственно снятия данных до анализа/построения модели и визуализации. А вывод ноутбуков по сути можно прямо вставлять в финальные отчеты.
                            0
                            Где изолента не держит — сварке нечего делать! Или! Если Вам не удается отремонтировать что-нибудь при помощи изоленты, значит у вас мало изоленты!

                            Но при обучении студентов желательно научить их пользоваться не только изолентой. Так что институт/лаборатория должна иметь сварочный аппарат, той или иной системы, и учить им пользоваться ;)

                            Если говорить за оптимизацию, то «средство оптимизации» вызывает «ваше ядро». Написание ядер (мат.модель) на Питоне:
                            • Зашоривает мозги, ибо без векторизации в Питоне всё хреново, поэтому целевые функции, которые не выражаются или трудно выражаются в векторзованном виде, изначально идут лесом;
                            • Структура задачи нелинейной оптимизации, даже нелинейного МНК/правдоподобия и т.п., несколько отличаются от линейных, квадратичных и прочих фурье-преобразований. Так что для целевых функций среднего и небольшого размера при сложном рельефе все ваши тесты производительности игнорируются и всё тормозит без всякого стыда и совести;
                            • Конечно 99% мат. моделей уходят в утиль в течении пары лет, но 1% выживает. И с мат.моделями из 80-х проблем не имел, а вот, бывает, встретится/притащут что-то на Бейсике/Excel/AXUM. Перепишешь на C/Fortran — ошибок по насажаешь, к гадалке не ходи. Обертку сделаешь — научную достоверность сохранишь, но… И с Питон через 20...30 лет тот же геморрой ожидается;
                              +1
                              Давайте еще раз повторим. C/Fortran — это очень хорошо для промышленных приложений.
                              У python есть свои сформированные ниши, например www.kaggle.com/c/two-sigma-financial-news/kernels — практически все исследования на python.
                              Мы уже достаточно подробно обсудили сильные и слабые стороны инструмента scipy. Что Вы еще хотите сказать?
                              PS
                              Где изолента не держит — сварке нечего делать!
                              Где русский инженер использует изоленту и водку, британский возьмет скотч и скотч ;-)
                        +2
                        Учитывая основной недостаток Питона — нелюбовь к многоядерным/многопроцессорным архитектурам

                        А у Матлаба что ли всё хорошо с многопоточностью/многопроцессорностью? Там вообще нет многопоточности на уровне m-языка, а мультипроцессинг убогий через parfor и MATLAB Parallel Server, с огромными накладными расходами на передачу данных и запуск процессов матлаба. numpy/scipy внутри себя освобождают gil и работают многопоточно.


                        SciPy не реализует ядро оптимизации, он использует древние проверенные фортрановские библиотеки, например NetLib, MINPACK и т. д. Их же использует Matlab. Думаю, вам как любителю фортрана знакомы эти названия.


                        На счет версий — это полная ерунда. Есть библиотеки, которые по 20 лет с версией 0.X и что же теперь ими пользоваться нельзя? Просто есть такая система версионирования. Называется ZeroVer 0-based Versioning.


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

                          0
                          В матлабе начина с 2017а с загрузкой процессора стало сильно лучше, но она на уровне матричных операций и индексирования. Ну и введение GPU arrays временами пригождается.
                      0
                      Спасибо за подробный комментарий, Вы все правильно говорите, единственно позволю себе некоторые уточнения.
                      «экономически просто нецелесообразна», т.е. машины на поездить или для посчитать целесообразны, а программное обеспечение нет. Странно немного.

                      В маленькой конторе затраты на покупку полноценной коммерческой версии ПО могут просто не окупиться из-за небольшой выручки. Или в госконторах бывают проблемы с финансированием или банальная бюрократия. Всяко бывает.
                      Я, конечно, ретроград, и предпочел бы их учить на примере исходных кодов NAG для Fortran

                      Учить студентов Фортрану конечно можно, но начинать полезнее все-таки по принципу чем попроще, тем лучше.
                      «интерактивность и простой синтаксис» — а как бы у MATLAB/Maple/IRAF интерактивность повыше будет. Не?

                      Лично мне работа в Anaconda нравится намного больше MATLAB, но это конечно субъективно.
                      Для сравнения Linux 1.0 появился 25 лет назад и только последнее время Linux смог претендовать на нишу более менее надёжного швейцарского ножа.

                      Технологии совместной разработки сейчас получше, чем 25 лет назад, да и ядро Linux посложнее scipy будет. В любом случае, продукт не растет сам по себе, нужны усилия сообщества, нужно рассказывать пользователям, что есть альтернатива ворованному ПО.
                        0
                        В любом случае, продукт не растет сам по себе, нужны усилия сообщества
                        Собственное я ж это с самого начала и сказал: «NAG не очень подходит для юношей с горящим взором, которые, прямо таки сразу, желают внести свою лепту в теорию и/или практику вычислительной математики». А Вы писали — однобоко ;)
                        Учить студентов Фортрану конечно можно, но начинать полезнее все-таки по принципу чем попроще, тем лучше.
                        Проще что? Проще писать хорошие вычислительные модели на C/FORTRAN? Или проще быстрее начать и позже решить?

                        Как бы Питон авторы предназначали совсем для иного — сетевые сервисы, тот же DJANGO, администрирование, работа со сложно организованными текстовыми данным сравнительно небольшого объёма. Ещё визуалировать, рассортировать наблюдательные данные он ещё туда-сюда. Хотя с картинками от того же НАСА/LROC на сотни мегабайт он уже киснет и плачет.

                        Я вот смотрю на некоторые проекты и наблюдаю непонятное. Типа на одном процессоре всё медленно, считать не очень умеем, с многими процессорами не дружим. А давайте использовать MPI, использовать много вычислительных узлов. Мало того, заодно для преодолений ограничений Питона будем относится к соседним ядрам своего процессора как к удалённым узлам.

                        Не знаю, не знаю, пока Питон находится в парадигме своего основателя Гвидо ван Россума — «90% программ однопроцессорные и, по большому счёту, однопоточные». Всё будет странно устроено, к гадалке не ходи. Мне так кажется.
                +1
                Очень хороший вопрос. К сожалению, готовых методик выбора алгоритма оптимизации я не встречал, а у самого пока мало опыта в решении подобных задач. Так что пока остается только перебор и интуиция.

                Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                Самое читаемое