<jupyter_start><jupyter_text>Введение в численные методы оптимизации (Ю. Е. Нестеров Введение в выпуклую оптимизацию, гл. 1 $\S$ 1.1) 1. Обзор материала весеннего семестра 2. Постановка задачи 3. Общая схема решения 4. Сравнение методов оптимизации 5. Методы одномерной минимизации  Обзор материала весеннего семестраТакже на [странице курса](https://github.com/amkatrutsa/MIPT-Optspring-term).1. Методы решения задач **безусловной** оптимизации    - Одномерная минимизация (**уже сегодня!**)    - Градиентный спуск    - Метод Ньютона    - Квазиньютоновские методы    - Метод сопряжённых градиентов      - Опционально:        - Решение задачи наименьших квадратов        - Оптимальные методы и нижние оценки2. Методы решения задач **условной** оптимизации    - Методы проекции градиента и условного градиента    - Методы штрафных и барьерных функций    - Метод модифицированой функции Лагранжа  Организация работы в семестре1. Семинар и лекция раз в неделю2. Два задания в течение семестра3. Midterm в середине семестра4. Итоговая контрольная в конце семестра5. Экзамен в конце семестра (схема выставления оценки аналогична осеннему семестру)6. Миниконтрольные в начале каждого семинара7. Домашнее задание почти каждую неделю: $\TeX$ или Jupyter Notebook  Постановка задачи\begin{equation}\begin{split}& \min_{x \in S} f_0(x)\\\text{s.t. } & f_j(x) = 0, \; j = 1,\ldots,m\\& g_k(x) \leq 0, \; k = 1,\ldots,p\end{split}\end{equation}где $S \subseteq \mathbb{R}^n$, $f_j: S \rightarrow \mathbb{R}, \; j = 0,\ldots,m$, $g_k: S \rightarrow \mathbb{R}, \; k=1,\ldots,p$Все функции как минимум непрерывны. Важный факт: задачи **нелинейной** оптимизации в их самой общей форме являются **численно неразрешимыми**!  Аналитические результаты- Необходимое условие первого порядка: если $x^*$ точка локального минимума дифференцируемой функции $f(x)$, тогда $$f'(x^*) = 0$$- Необходимое условие второго порядка если $x^*$ точка локального минимума дважды дифференцируемой функции $f(x)$, тогда $$f'(x^*) = 0 \quad \text{и} \quad f''(x^*) \succeq 0$$- Достаточное условие: пусть $f(x)$ дважды дифференцируемая функция, и пусть точка $x^*$ удовлетворяет условиям$$f'(x^*) = 0 \quad f''(x^*) \succ 0,$$тогда $x^*$ является точкой строго локального минимума функции $f(x)$.**Замечание**: убедитесь, что Вы понимаете, как доказывать этирезультаты!  Особенности численного решения1. Точно решить задачу принципиально невозможно из-за погрешности машинной арифметики2. Необходимо задать критерий обнаружения решения3. Необходимо определить, какую информацию о задаче использовать  Общая итеративная схемаДано: начальное приближение $x$, требуемая точность $\varepsilon$.```pythondef GeneralScheme(x, epsilon):        while StopCriterion(x) > epsilon:                OracleResponse = RequestOracle(x)                UpdateInformation(I, x, OracleResponse)                x = NextPoint(I, x)            return x```  Вопросы1. Какие критерии остановки могут быть?2. Что такое оракул и зачем он нужен?3. Что такое информационная модель?4. Как вычисляется новая точка?  Критерии остановки1. Сходимость по аргументу: $$\| x_k - x^*  \|_2 < \varepsilon$$ 2. Сходимость по функции: $$\| f_k - f^* \|_2 < \varepsilon$$ 3. Выполнение необходимого условия $$\| f'(x_k) \|_2 < \varepsilon$$Но ведь $x^*$ неизвестна!Тогда\begin{align*}& \|x_{k+1} - x_k \| = \|x_{k+1} - x_k + x^* - x^* \| \leq \\& \|x_{k+1} - x^* \| + \| x_k - x^* \| \leq 2\varepsilon\end{align*}Аналогично для сходимости по функции, однако иногда можно оценить $f^*$! **Замечание**: лучше использовать относительные изменения этих величин! Например $\dfrac{\|x_{k+1} - x_k \|_2}{\| x_k \|_2}$  Что такое оракул?**Определение**: оракулом называют некоторое абстрактное устройство, которое отвечает на последовательные вопросы методаАналогия из ООП: - оракул - это виртуальный метод базового класса- каждая задача - производный класс- оракул определяется для каждой задачи отдельно согласно общему определению в базовом классе**Концепция чёрного ящика**1. Единственной информацией, получаемой в ходе работы итеративного метода, являются ответы оракула2. Ответы оракула являются *локальными*  Информация о задаче1. Каждый ответ оракула даёт **локальную** информацию о поведении функции в точке2. Агрегируя все полученные ответы оракула, обновляем информацию о **глобальном** виде целевой функции:    - кривизна    - направление убывания    - etc  Вычисление следующей точки$$x_{k+1} = x_{k} + \alpha_k h_k$$- **Линейный поиск**: фиксируется направление $h_k$ и производится поиск по этому направлению "оптимального" значения  $\alpha_k$- **Метод доверительных областей**: фиксируется допустимый размер *области* по некоторой норме $\| \cdot  \| \leq \alpha$ и *модель* целевой функции, которая хорошо её аппроксимирует в выбранной области.         Далее производится поиск направления $h_k$, минимизирующего модель целевой функции и не выводящего точку $x_k + h_k$ за пределы доверительной области  Вопросы1. Как выбрать $\alpha_k$?2. Как выбрать $h_k$?3. Как выбрать модель?4. Как выбрать область?5. Как выбрать размер области?     В курсе рассматривается только линейный поиск!     Однако несколько раз копцепция метода доверительных областей будет использована.  Как сравнивать методы оптимизации?Для заданного класса задач сравнивают следующие величины:1. Сложность    - аналитическая: число обращений к оракулу для решения задачи с точностью $\varepsilon$    - арифметическая: общее число всех вычислений, необходимых для решения задачи с точностью $\varepsilon$2. Скорость сходимости3. Эксперименты  Скорости сходимости 1. Сублинейная$$\| x_{k+1} - x^* \|_2 \leq C k^{\alpha},$$где $\alpha < 0$ и $ 0 < C < \infty$2. Линейная (геометрическая прогрессия)$$\| x_{k+1} - x^* \|_2 \leq Cq^k, $$где $q \in (0, 1)$ и $ 0 < C < \infty$3. Сверхлинейная $$\| x_{k+1} - x^* \|_2 \leq Cq^{k^p}, $$где $q \in (0, 1)$, $ 0  1$4. Квадратичная$$\| x_{k+1} - x^* \|_2 \leq C\| x_k - x^* \|^2_2, \qquad \text{или} \qquad \| x_{k+1} - x^* \|_2 \leq C q^{2^k}$$где $q \in (0, 1)$ и $ 0 < C < \infty$  Оптимальные методы: can we do better?- Доказывают нижние оценки скоростей сходимости для класса задач и методов фиксированного порядка- Предлагают методы, на которых эти нижние оценки достигаются $\Rightarrow$ доказана оптимальность - Ниже про значение теорем сходимостиОптимальным методам и нижним оценкам будет, возможно, посвящён отдельный семинар или часть домашнего задания.<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt

USE_COLAB = False
if not USE_COLAB:
    plt.rc("text", usetex=True)

import numpy as np
C = 10
alpha = -0.5
q = 0.9
num_iter = 7
sublinear = np.array([C * k**alpha for k in range(1, num_iter + 1)])
linear = np.array([C * q**k for k in range(1, num_iter + 1)])
superlinear = np.array([C * q**(k**2) for k in range(1, num_iter + 1)])
quadratic = np.array([C * q**(2**k) for k in range(1, num_iter + 1)])
plt.figure(figsize=(12,8))
plt.semilogy(np.arange(1, num_iter+1), sublinear, 
             label=r"Sublinear, $\alpha = -0.5$")
plt.semilogy(np.arange(1, num_iter+1), superlinear, 
             label=r"Superlinear, $q = 0.5, p=2$")
plt.semilogy(np.arange(1, num_iter+1), linear, 
             label=r"Linear, $q = 0.5$")
plt.semilogy(np.arange(1, num_iter+1), quadratic, 
             label=r"Quadratic, $q = 0.5$")
plt.xlabel("Number of iterations, $k$", fontsize=28)
plt.ylabel("Error rate upper bound", fontsize=28)
plt.legend(loc="best", fontsize=26)
plt.xticks(fontsize = 28)
_ = plt.yticks(fontsize = 28)<jupyter_output><empty_output><jupyter_text>Значение теорем сходимости (<NAME> в оптимизацию, гл. 1, $\S$ 6)1. Что дают теоремы сходимости     - класс задач, для которых можно рассчитывать на применимость метода (важно не завышать условия!)         - выпуклость         - гладкость     - качественное поведение метода         - существенно ли начальное приближение         - по какому функционалу есть сходимость     - оценку скорости сходимости         - теоретическая оценка поведения метода без проведения экспериментов         - определение факторов, которые влияют на сходимость (обусловленность, размерность, etc)         - иногда заранее можно выбрать число итераций для достижения заданной точности 2. Что **НЕ** дают теоремы сходимости    - сходимость метода **ничего не говорит** о целесообразности его применения    - оценки сходимости зависят от неизвестных констант - неконструктивный характер    - учёт ошибок округления и точности решения вспомогательных задач    **Мораль**: нужно проявлять разумную осторожность и здравый смысл!  Классификация задач1. Безусловная оптимизация    - целевая функция липшицева    - градиент целевой функции липшицев2. Условная оптимизация    - многогранник    - множество простой структуры    - общего вида  Классификация методов1. Методы нулевого порядка: оракул возвращает только значение функции $f(x)$2. Методы первого порядка: оракул возвращает значение функции $f(x)$ и её градиент $f'(x)$3. Методы второго порядка: оракул возвращает значение функции $f(x)$, её градиент $f'(x)$ и гессиан $f''(x)$.**Вопрос**: существуют ли методы более высокого порядка?1. Одношаговые методы $$x_{k+1} = \Phi(x_k)$$2. Многошаговые методы$$x_{k+1} = \Phi(x_k, x_{k-1}, ...)$$  Одномерная минимизация**Определение**. Функция $f(x)$ называется унимодальной на $[a, b]$, если существует такая точка $x^* \in [a, b]$, что - $f(x_1) > f(x_2)$ для любых $a \leq x_1 < x_2 < x^*$, и - $f(x_1) < f(x_2)$ для любых $x^* < x_1 < x_2 \leq b$.**Вопрос**: какая геометрия унимодальных функций?  Метод дихотомииИдея из информатики первого семестра: делим отрезок $[a,b]$ на две равные части пока не найдём минимум унимодальной функции.- $N$ - число вычислений функции $f$- $K = \frac{N - 1}{2}$ - число итерацийТогда$$|x_{K+1} - x^*| \leq \frac{b_{K+1} - a_{K+1}}{2} = \left( \frac{1}{2} \right)^{\frac{N-1}{2}} (b - a) \approx 0.5^{K} (b - a) $$<jupyter_code>def binary_search(f, a, b, epsilon, callback=None):
    c = (a + b) / 2.0
    while abs(b - a) > epsilon:
#         Check left subsegment
        y = (a + c) / 2.0
        if f(y) <= f(c):
            b = c
            c = y
        else:
#         Check right subsegment
            z = (b + c) / 2.0
            if f(c) <= f(z):
                a = y
                b = z
            else:
                a = c
                c = z
        if callback is not None:
            callback(a, b)
    return c
def my_callback(a, b, left_bound, right_bound, approximation):
    left_bound.append(a)
    right_bound.append(b)
    approximation.append((a + b) / 2.0)
import numpy as np

left_boud_bs = []
right_bound_bs = []
approximation_bs = []

callback_bs = lambda a, b: my_callback(a, b, 
            left_boud_bs, right_bound_bs, approximation_bs)

# Target unimodal function on given segment
f = lambda x: (x - 2) * x * (x + 2)**2 # np.power(x+2, 2)
# f = lambda x: -np.sin(x)
x_true = -2
# x_true = np.pi / 2.0
a = -3
b = -1.5
epsilon = 1e-8
x_opt = binary_search(f, a, b, epsilon, callback_bs)
print(np.abs(x_opt - x_true))
plt.figure(figsize=(10,6))
plt.plot(np.linspace(a,b), f(np.linspace(a,b)))
plt.title("Objective function", fontsize=28)
plt.xticks(fontsize = 28)
_ = plt.yticks(fontsize = 28)<jupyter_output>9.313225746154785e-10<jupyter_text>Метод золотого сеченияИдея: делить отрезок $[a,b]$ не на две равные насти, а в пропорции "золотого сечения".Оценим скорость сходимости аналогично методу дихотомии:$$|x_{K+1} - x^*| \leq b_{K+1} - a_{K+1} = \left( \frac{1}{\tau} \right)^{N-1} (b - a) \approx 0.618^K(b-a),$$где $\tau = \frac{\sqrt{5} + 1}{2}$.- Константа геометрической прогрессии **больше**, чем у метода дихотомии- Количество вызовов функции **меньше**, чем у метода дихотомии<jupyter_code>def golden_search(f, a, b, tol=1e-5, callback=None):
    tau = (np.sqrt(5) + 1) / 2.0
    y = a + (b - a) / tau**2
    z = a + (b - a) / tau
    while b - a > tol:
        if f(y) <= f(z):
            b = z
            z = y
            y = a + (b - a) / tau**2
        else:
            a = y
            y = z
            z = a + (b - a) / tau
        if callback is not None:
            callback(a, b)
    return (a + b) / 2.0
left_boud_gs = []
right_bound_gs = []
approximation_gs = []

cb_gs = lambda a, b: my_callback(a, b, left_boud_gs, right_bound_gs, approximation_gs)
x_gs = golden_search(f, a, b, epsilon, cb_gs)

print(f(x_opt))
print(f(x_gs))
print(np.abs(x_opt - x_true))<jupyter_output>6.93889390875399e-18
9.549014390504221e-18
9.313225746154785e-10<jupyter_text>Сравнение методов одномерной минимизации<jupyter_code>plt.figure(figsize=(10,6))
plt.semilogy(np.arange(1, len(approximation_bs) + 1), np.abs(x_true - np.array(approximation_bs, dtype=np.float64)), label="Binary search")
plt.semilogy(np.arange(1, len(approximation_gs) + 1), np.abs(x_true - np.array(approximation_gs, dtype=np.float64)), label="Golden search")
plt.xlabel(r"Number of iterations, $k$", fontsize=26)
plt.ylabel("Error rate upper bound", fontsize=26)
plt.legend(loc="best", fontsize=26)
plt.xticks(fontsize = 26)
_ = plt.yticks(fontsize = 26)
%timeit binary_search(f, a, b, epsilon)
%timeit golden_search(f, a, b, epsilon)<jupyter_output>22.6 µs ± 1.13 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
84.1 µs ± 6.64 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)<jupyter_text>Пример иного поведения методов$$f(x) = \sin(\sin(\sin(\sqrt{x}))), \; x \in [2, 60]$$<jupyter_code>f = lambda x: np.sin(np.sin(np.sin(np.sqrt(x))))
x_true = (3 * np.pi / 2)**2
a = 2
b = 60
epsilon = 1e-8
plt.plot(np.linspace(a,b), f(np.linspace(a,b)))
plt.xticks(fontsize = 28)
_ = plt.yticks(fontsize = 28)<jupyter_output><empty_output><jupyter_text>Сравнение скорости сходимости и времени работы методов  Метод дихотомии<jupyter_code>left_boud_bs = []
right_bound_bs = []
approximation_bs = []

callback_bs = lambda a, b: my_callback(a, b, 
            left_boud_bs, right_bound_bs, approximation_bs)

x_opt = binary_search(f, a, b, epsilon, callback_bs)
print(np.abs(x_opt - x_true))<jupyter_output>2.1968899233115735e-07<jupyter_text>Метод золотого сечения<jupyter_code>left_boud_gs = []
right_bound_gs = []
approximation_gs = []

cb_gs = lambda a, b: my_callback(a, b, left_boud_gs, right_bound_gs, approximation_gs)
x_gs = golden_search(f, a, b, epsilon, cb_gs)

print(np.abs(x_opt - x_true))<jupyter_output>2.1968899233115735e-07<jupyter_text>Сходимость<jupyter_code>plt.figure(figsize=(8,6))
plt.semilogy(np.abs(x_true - np.array(approximation_bs, dtype=np.float64)), label="Binary")
plt.semilogy(np.abs(x_true - np.array(approximation_gs, dtype=np.float64)), label="Golden")
plt.legend(fontsize=28)
plt.xticks(fontsize=28)
_ = plt.yticks(fontsize=28)
plt.xlabel(r"Number of iterations, $k$", fontsize=26)
plt.ylabel("Error rate upper bound", fontsize=26)<jupyter_output><empty_output><jupyter_text>Время работы<jupyter_code>%timeit binary_search(f, a, b, epsilon)
%timeit golden_search(f, a, b, epsilon)<jupyter_output>455 µs ± 35.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
416 µs ± 26.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
{'fraction_non_alphanumeric': 0.13357712765957447, 'fraction_numerical': 0.028656914893617023, 'mean_word_length': 4.027072192513369, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 36, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '458346', 'license': 'MIT', 'n_tokens_mistral': 6529, 'n_tokens_neox': 7022, 'n_words': 1567, 'path': 'Spring2017-2019/12-NumMethods/Seminar12.ipynb', 'repo_name': 'Jhomanik/MIPT-Opt', 'size': '15053'}