<jupyter_start><jupyter_text>Метод эллипсоидовЭтот метод был придуман в 1979 году и является первым полиномиальным для задачи линейного программирования. Метод эллипсоидов является вариацией более старого метода разделяющих плоскостей (cutting plane method). Для задачи минимизации функции $f$ без ограничений этот метод имеет следующий вид:* Изначально у нас есть многогранник $P_0$, который гарантированно содержит минимум $f$ и начальное приближение $x_0\in P_0$.* На итерации $k$   * найти такой вектор $g_k$, что выполняется$$x^*\in P_k\cap\{x~|~g_k^T(x-x_k)\leq 0\}$$  * Взять$$P_{k+1}=P_k\cap\{x~|~g_k^T(x-x_k)\leq 0\}, ~x_{k+1}\in P_{k+1}$$Для выпуклых функций в качестве $g_k$ обычно выбирается как субградиент $f(x_k)$, из определения субградиента$$f(y)\geq f(x_k)+g_k^T(y-x_k)$$это означает, что если $g_k^T(y-x_k)>0$, то $f(y)> f(x_k)$, из чего следует, что $x^*\in\{x~|~g_k^T(x-x_k)\leq 0\}$. Точка $x_{k+1}$ может быть выбрана произвольным образом, по сути в большей степени определяет где стоит сделать отсечение на следующей итерации. В целом выбор точки следует производить таким образом, чтобы в худшем/среднем случае мы отсекали как можно больше, чтобы сходимость была быстрее, таким образом в целом стоит выбирать новое приближение как некоторый центр многогранника. Здесь есть несколько возможных подходов:* Центр масс$$x_k=\frac{\int_{P_k}xdx}{\int_{P_k}dx}$$к сожалению вычисление $x_k$ обычно является более сложной задачей, чем исходная. В теории однако этот метод дает очень хорошую оценку сходимости$$\mathbf{vol}P_{k+1}\leq (1-1/e)\mathbf{vol}P_k$$вне зависимости от размерности пространства.* Центр эллипсоида максимального размера: вычисляется решением вспомогательной выпуклой задачи, дает оценку$$\mathbf{vol}P_{k+1}\leq (1-1/n)\mathbf{vol}P_k$$* Аналитический центр: при $P_k=\{x~|~a_ix\leq b_i\}$$$\DeclareMathOperator*\argmin{argmin}x_k=\argmin_{x}-\sum_{i=1}^m\log(b_i-a_ix)$$ Проблема этого метода во многом обусловлена тем, что в большинстве случаев стоимость итерации растет со временем, при этом дополнительные трудности возникают из-за необходимости найти хороший центр для следующего приближения. Обе этих проблемы решаются в методе эллипсоидов путем ухудшения скорости сходимости.Метод эллипсоидов работает по тому же принципу, что и метод разделающих плоскостей, но теперь множество $P_k$ является эллипсоидом, который в общем случае задается следующим образом$$P_k=\{x~|~(x-x_k)^TE(x-x_k)\leq 1 \}, E\succeq 0$$Это автоматически решает проблему выбора центра эллипсоида, которым является точка $x_k$. В силу того, что $P_k\cap\{x~|~g_k^T(x-x_k)\leq 0\}$ не является элипсоидом, в качестве $P_{k+1}$ берется эллипсоид минимального объема, содержащий полуэлипсоид $P_k\cap\{x~|~g_k^T(x-x_k)\leq 0\}$, который задается по следующим формулам:$$\begin{array}{rl}x_{k+1}&=x_k-\frac{1}{n+1}E_kg_k\\\hat{g}_k&=\frac{1}{\sqrt{g_kE_kg_k}}g_k \\E_{k+1}=&\frac{n^2}{n^2-1}\left(E_k-\frac{2}{n+1}E_k\hat{g}_k\hat{g}_k^TE_k\right)\end{array}$$Для этого метода выполняется оценка$$\mathbf{vol}P_{k+1}<\exp\left(-\frac{1}{2n}\right)\mathbf{vol}P_k$$что делает его линейным по точности приближения. На практике однако этот алгоритм не прижился из-за того, что в большинстве случаев заметно проигрывал симплекс-методу и, появившемуся позже, методу внутренней точки. Метод эллипсоидов в задачах с ограничениямиДля начала заметим, что если нам известно, что текущее приближение $x_k$ хуже некоторого $x_l$, которое мы видели до этого момента, то вместо полуплоскости $$g_k^T(y-x_k)\leq 0$$можно взять гиперплоскость$$g_k^T(y-x_k)\leq f(x_l)-f(x_k)$$как и раньше, если для $y$ лежит в оставшейся части, то$$f(y)\geq f(x_k)+g_k^T(y-x_k)> f(x_l)\geq f(x^*)$$что все еще дает нам гарантию, что мы отсекаем ту часть, которая не содержит точки минимума. Такое отсечение называется глубоким разрезом и в применении к методу эллипсоидов имеет следующий вид$$\begin{array}{rl}\alpha_k&=\frac{f(x_k)-f(x_{\min})}{\sqrt{g_kE_kg_k}} \\\hat{g}_k&=\frac{1}{\sqrt{g_kE_kg_k}}g_k \\x_{k+1}&=x_k-\frac{1+\alpha_k n}{n+1}E_kg_k\\E_{k+1}=&\frac{n^2}{n^2-1}(1-\alpha_k)^2\left(E_k-\frac{2(1+\alpha_k n)}{(n+1)(1+\alpha_k)}E_k\hat{g}_k\hat{g}_k^TE_k\right)\end{array}$$здесь мы используем глубокий разрез для того, чтобы улучшить скорость сходимости, но схожая техника может быть использована, чтобы решать задачу с ограничениями. Идея проста: если на в какой-то момент центр эллипсоида нарушает некоторое ограничение $h_i(x)\leq 0$, то отсечение эллипсоида происходит по субградиенту функции $h_i$ со значением $h_i(x_k)$.<jupyter_code>def ellipsoid_step_LP(x_k, P_k, direction, A, b, prev_best):
    num = -1
    for i in range(b.shape[0]):
        if A[i].dot(x_k) > b[i]:
            num = i
            break
            
    if num != -1:
        #print('Infeasible')
        g_k = A[num]
        alpha = (A[num].dot(x_k) - b[num]) / np.sqrt(g_k.dot(P_k.dot(g_k)))
        g_k = g_k / np.sqrt(g_k.dot(P_k.dot(g_k)))
        best = prev_best
    else:
        #print('Feasible')
        g_k = direction
        alpha = max((direction.dot(x_k) - prev_best), 0) / np.sqrt(g_k.dot(P_k.dot(g_k)))
        g_k = g_k / np.sqrt(g_k.dot(P_k.dot(g_k)))
        best = min(prev_best, direction.dot(x_k))
        
    x_new = x_k - (1.0 + 2.0 * alpha) / 3.0 * P_k.dot(g_k)
    #print(x_new)
    #print(np.outer(P_k.dot(g_k), P_k.dot(g_k)))
    P_new = 4.0 / 3.0 * (1 - alpha ** 2) * (P_k - 2.0 * (1.0 + 2.0 * alpha) / 3.0 / (1.0 + alpha) * np.outer(P_k.dot(g_k), P_k.dot(g_k)))  
    #print('Eig', np.linalg.eig(P_k))
    #print('Eig', np.linalg.eig(P_new))
    return x_new, P_new, best
    
def ellipse_generating(P, center, x):
    d = x - center
    P_1 = np.linalg.inv(P)
    return d.dot(P_1.dot(d))
import matplotlib.pyplot as plt
import numpy as np
from interactive_visualization.animation_utils import animate_list
def fix_scaling(ax=None):
    if not ax:
        xlim = plt.xlim()
        ylim = plt.ylim()
        d1 = xlim[1] - xlim[0]
        d2 = ylim[1] - ylim[0]
        if d1 > d2:
            plt.ylim((ylim[0] - (d1-d2) / 2, ylim[1] + (d1-d2) / 2))
        else:
            plt.xlim((xlim[0] + (d1-d2) / 2, xlim[1] - (d1-d2) / 2))
    else:
        xlim = ax.get_xlim()
        ylim = ax.get_ylim()
        d1 = xlim[1] - xlim[0]
        d2 = ylim[1] - ylim[0]
        if d1 > d2:
            ax.set_ylim((ylim[0] - (d1-d2) / 2, ylim[1] + (d1-d2) / 2))
        else:
            ax.set_xlim((xlim[0] + (d1-d2) / 2, xlim[1] - (d1-d2) / 2))
def get_line(x1, x2):
    a = x1[1] - x2[1]
    b = x2[0] - x1[0]
    c = a * x1[0] + b * x1[1]
    return a, b, c

vertices = [(2.0, 2.0), (1.9, 3.0), (2.5, 4.0), (4.0, 4.2), (4.7, 3.5), (4.5, 1.5), (3.5, 1.0), (2.0, 2.0)]
A = []
b = []

for i in range(len(vertices) - 1):
    a_, b_, c_ = get_line(vertices[i], vertices[i + 1])
    A.append([a_, b_])
    b.append(c_)
A = np.array(A)
b = np.array(b)
direction = np.array([-2, -1]) # c
P = 3.3 * np.identity(2)
x_cent = np.array([3.5, 2.8])
best = x_cent.dot(direction)

my_colors = ['red', 'green', 'magenta']
_num_steps = 11
some_data_to_animate = [[x_cent.copy(), P.copy()]]
for k in range(_num_steps + 1):
    x_cent, P, best = ellipsoid_step_LP(x_cent, P, direction, A, b, best)
    some_data_to_animate.append([x_cent.copy(), P.copy()])
#Level contours
def ellipse_state(t, data):
    fig, ax = plt.subplots(figsize=(10, 10))
    delta = 0.025
    x = np.arange(1, 6, delta)
    y = np.arange(0.5, 6, delta)
    X, Y = np.meshgrid(x, y)
    Z = np.zeros_like(X)
    #print(X.shape, Y.shape)
    P = 3.3 * np.identity(2)
    x_cent = np.array([3.5, 2.8])
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            #print(ellipse_generating(P, x_cent, np.array([X[i][j], Y[i][j]])))
            Z[i][j] = ellipse_generating(P, x_cent, np.array([X[i][j], Y[i][j]]))
    CS = ax.contour(X, Y, Z, [1], colors=['green'])
    ax.plot([x for x, y in vertices], [y for x, y in vertices])

    x_cent, P = data[t // 2]
    odd = t % 2 == 0
    num = -1
    for j in range(len(b)):
        if np.array(A[j]).dot(np.array(x_cent)) > b[j]:
            num = j
            break
            
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            Z[i][j] = ellipse_generating(P, x_cent, np.array([X[i][j], Y[i][j]]))
    CS = ax.contour(X, Y, Z, [1], colors=['green'])

    
    if num != -1:
        for i in range(X.shape[0]):
            for j in range(X.shape[1]):
                Z[i][j] = A[num][0] * X[i][j] + A[num][1] * Y[i][j] - b[num]
        CS = ax.contour(X, Y, Z, [0], colors=['black'])
    else:
        for i in range(X.shape[0]):
            for j in range(X.shape[1]):
                Z[i][j] = direction[0] * (X[i][j] - x_cent[0]) + direction[1] * (Y[i][j] - x_cent[1])
        CS = ax.contour(X, Y, Z, [0], colors=['black'])
    
    ax.plot([x_cent[0]], [x_cent[1]], 'o', color='black')
    if not odd:
        x_cent, P = some_data_to_animate[t // 2 + 1]            
        for i in range(X.shape[0]):
            for j in range(X.shape[1]):
                Z[i][j] = ellipse_generating(P, x_cent, np.array([X[i][j], Y[i][j]]))
        CS = ax.contour(X, Y, Z, [1], colors=['red'])
    ax.plot([x for x, y in vertices], [y for x, y in vertices])
    fix_scaling(ax)
    ax.axis('off')
    plt.close(fig)
    return fig
animate_list([ellipse_state(t, some_data_to_animate) for t in range(len(some_data_to_animate))])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1724907063197026, 'fraction_numerical': 0.025703664365374402, 'mean_word_length': 3.3795348837209302, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 12, 'https://': 0, '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': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '321572', 'license': 'MIT', 'n_tokens_mistral': 4367, 'n_tokens_neox': 4619, 'n_words': 921, 'path': 'optimization_course/12_ellipsoid_method.ipynb', 'repo_name': 'ivanoleinik/interactive-visualization', 'size': '9415'}