Весьма нетривиальная оптимизация

[_genius_]

В общем изрядно поиздевавшись над "Леннардом Джонсом" пришел к следующему выражению:
$h(d) = \frac{1}{d}-(\frac{0.1}{d})^{10}$
Как мне кажется обрыв достаточно резкий чтобы отбить всяческое желание сближать точки ближе чем $0.1$. Завтра проверю.

[_genius_]

vpb в сообщении #1357281 писал(а):

$$ f(\Phi)=\sum_{i=1,\ldots,N;\ j=1,\ldots,M} h(d(x_i, \Phi y_j)) (f(x_i)-g(y_j)) ? $$

где $x_i$ --- это точки первой поверхности, $M$ штук, $y_j$ --- точки второй, $N$ штук, $d(x,y)$ --- обычное расстояние между точками в трехмерном пространстве, $\Phi$ означает то самое движение (а $\Phi y$ соответственно --- результат применения движения к точке $y$) ?

vpb в сообщении #1357558 писал(а):

Есть такая возможность ускорить вычисление $f(\Phi)$. Рассмотрим функции от $x$ и $y$:

$$ p(x)=\sum_{j=1}^M h(d(x,y_j))g(y_j)\,, \qquad q(y)=\sum_{i=1}^N h(d(x_i,y)) f(x_i)\,.$$

Тогда легко видеть, что

$$f(\Phi)=\sum_{j=1}^Mq(\Phi y_j)-\sum_{i=1}^Np(\Phi^{-1}x_i)$$

В общем я только сейчас понял ещё один неприятный момент - хоть я и писал ранее "разность энергий близлежащих (одна от первой поверхности, другая от второй) точек", на деле надо учитывать и когда наборот.
Т.е. на деле мне нужна не разность, а модуль разности (f и g нормированы так чтобы лежать в интервале $\pm 0.5$)

$$ f(\Phi)=\sum_{i=1,\ldots,N;\ j=1,\ldots,M} h(d(x_i, \Phi y_j)) \lvert f(x_i)-g(y_j) \rvert ? $$

При этом с $f(x_i)$ и $g(y_j)$ можно делать все что вздумается (сдвигать/масштабировать), лишь бы сохранялся принцип - маленькое f хочет к большому g, а маленькое g хочет к большому f.

vpb, возможно ли все ещё ускорить вычисления при таком раскладе?

[_genius_]

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

[vpb]

Подумать надо, однако !

[_genius_]

Так и не смог раскрыть честно, поэтому решил просто прикинуть как оно должно выглядеть чтобы было красиво:

$$ p(x)=\sum_{j=1}^M h(d(x,y_j))g(y_j)\,, \qquad q(y)=\sum_{i=1}^N h(d(x_i,y)) f(x_i)\ $$

$$ pp(x)=\sum_{j=1}^M h(d(x,y_j))g^2(y_j)\,, \qquad qq(y)=\sum_{i=1}^N h(d(x_i,y)) f^2(x_i)\ $$

$$ f(\Phi)=\sum_{j=1}^Mqq(\Phi y_j)+\sum_{i=1}^Npp(\Phi^{-1}x_i)-\sum_{j=1}^Mg(y_j)q(\Phi y_j)-\sum_{i=1}^Nf(x_i)p(\Phi^{-1}x_i) $$

Только что опробировал с $h(d) = \frac{1}{d}-(\frac{0.1}{d})^{10}$ и результаты хоть и адекватные (нет постоянных вылетов за границы, несходимости и т.д.), но противоречат ожиданиям - теперь второе тело отчаянно стремится залезть в первое.
Завтра попробую с какой-нибудь экспонентой в качестве $h$. Просто не могу быть уверен, проблема в моем $h$ или в интуитивном нахождении $f$.

[_genius_]

Ах да, $\sum_{j=1}^Mg(y_j)q(\Phi y_j) $ и $ \sum_{i=1}^Nf(x_i)p(\Phi^{-1}x_i) $ равны друг другу вплоть до 5го знака, очевидно что незначительная ошибка исходит из интерполяции. Это подверждает что все считается как я думаю. Вопрос лишь в том правильно ли это.

(Оффтоп)

upd: в поисках ошибок нахожу даже то, чего нет

-- Чт дек 20, 2018 11:59:42 --

В общем сегодня ещё попробую $h(d) = \frac{100}{d}-(\frac{1}{d})^{10}$, т.е. теперь пик в 1, а не в 0.1. Теперь внутри первого тела точно не должно хватить места для второго.

[_genius_]

$h(d) = \frac{100}{d}-(\frac{1}{d})^{10}$ оправдала ожидания, теперь наложения практически нет (если использовать Пауэлла и Нелдер-Мида, ЛБФГС да и любой другой метод использующий производные как мне кажется здесь не справляется с задачей), однако чисто на глаз мне кажется что найденное положение отнюдь не самое лучшее, хотя и неплохое.
Попробую ещё сдвинуть пик например в 3. Экспонента теперь точно отменяется: в ней нет смысла, т.к. автоматическое предотвращение наложения перевешивает все минусы.

[_genius_]

$h(d) = \frac{1}{d^3}-(\frac{1}{d})^{10}$
Поигравшись со стартовыми углами поворота (раньше они были нулевые) и стартовыми сдвигами убедился что увы наложение очень даже возможно. Не знаю почему так получается. Как мне казалось с таким резким спадом отсутствие наложений должно быть гарантировано. Хотя быть может он слишком резкий и поэтому любой дальнейший сдвиг приводит к ухудшению ситуации и как следствие метод застревает.
https://i.ibb.co/dt4kPWF/untitled.png

[vpb]

_genius_
Как я раньше Вам писал в ЛС, я более-менее решил вопрос, как ускорить вычисления с новой функцией (содержащей $|f(x_i)-g(y_j)|$). Но аккуратно это написать требует некоторого времени. Имейте терпение. Начало, однако, я уже написал,
см.ниже.

Вы пишете что-то странное:

_genius_ в сообщении #1362603 писал(а):

Так и не смог раскрыть честно, поэтому решил просто прикинуть как оно должно выглядеть чтобы было красиво:

$$ pp(x)=\sum_{j=1}^M h(d(x,y_j))g^2(y_j)\,, \qquad qq(y)=\sum_{i=1}^N h(d(x_i,y)) f^2(x_i)\ $$

$$ f(\Phi)=\sum_{j=1}^Mqq(\Phi y_j)+\sum_{i=1}^Npp(\Phi^{-1}x_i)-\sum_{j=1}^Mg(y_j)q(\Phi y_j)-\sum_{i=1}^Nf(x_i)p(\Phi^{-1}x_i) $$

Но, собственно, к задаче то это не имеет отношения ! Это очень напоминает то, как школьник (причем бестолковый), "решая" задачу, комбинирует данные в условии числа каким-то случайным образом, типа складывает метры с килограммами. В общем, подозреваю, что Вы в логике задачи путаетесь.

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

Сначала несколько общих слов.
1) Произошла накладка в обозначениях, которую я и сам не заметил: употреблял букву $f$ в двух разных смыслах. Пусть то, что было $f(\Phi)$, обозначается $E(\Phi)$.
2) Разумно будет у функции $h$ поменять знак, так как это на самом деле потенциальная энергия: $h(d)=(0,1/d)^{10}-(1/d)$. Кроме того, в задачах оптимизации обычно ведется речь о нахождении минимума какой-либо функции, а не максимума. Имея же дело с потенциальной энергией, сам бог велел искать минимум.
3) Функции $f$ и $g$ сдвинем, из соображений удобства. Считаем, что они принимают значения в отрезке $I=[0;1]$.

Итак, пусть $x_i$, $i=1,\ldots,N$, $y_j$, $j=1,\ldots,M$, $f(x),g(y)\in[0;1]$ --- как раньше, $h(d)=(0,1/d)^{10}-(1/d)$,

$$ E(\Phi)=\sum_{i,j} h(d(x_i,\Phi y_j))\cdot |f(x_i)-g(y_j)|, $$

$$ p(x)=\sum_{j=1}^M h(d(x,y_j)) g(y_j)\, \qquad q(y)=\sum_{i=1}^N h(d(x_i,y)) f(x_i)\,.$$

Собственно, функции $p(x)$, $q(y)$ как таковые нам нужны не будут. Нужны же будут некоторые их модификации.

Разобьем отрезок $I$ на несколько частей, скажем, на три равных: $I_1=[0; 1/3]$, $I_2=[1/3; 2/3]$, $I_3=[2/3; 1]$. Также для дальнейших вычислений, заметим, удобно будет перенумеровать точки $x_i$ в порядке возрастания функции $f(x_i)$, и аналогично поступить с $y_j$.

Основная мысль --- разбить множество всех пар $(i,j)$ на девять частей, согласно тому, в какой из отрезков $I_{1,2,3}$ попадают $f(x_i)$ и $g(y_j)$. А именно, пусть

$$ E_{ab}(\Phi)=\sum'_{i,j} h(d(x_i,\Phi y_j))\cdot |f(x_i)-g(y_j)|, $$

где штрих означает, что сумма берется по парам $(i,j)$, для которых $f(x_i)\in I_a$, $g(y_j)\in I_b$, $a,b=1,2,3$. Ясно, что $E(\Phi)=\sum_{a,b=1}^3 E_{ab}(\Phi)$.

Суммы $E_{aa}(\Phi)$ вычисляются прямо. А при $a\ne b$ можно сэкономить. Например, если $a=1$, $b=2$, то $|f(x_i)-g(y_j)|=g(y_j)-f(x_i)$, значит сумма есть

$$ E_{12}(\Phi)=\sum'_{i,j} h(d(x_i,\Phi y_j))\cdot (g(y_j)-f(x_i)) = $$

$$ = -\sum_{\{j\mid g(y_j)\in I_2\}} q_1(\Phi y_j) + \sum_{\{i\mid f(x_i)\in I_1\}} p_2(\Phi^{-1} x_i), $$

где

$$ q_1(y)=\sum_{\{i\mid f(x_i)\in I_1\}} h(d(x_i,y)) f(x_i)\,, \qquad p_2(x)=\sum_{\{j\mid g(y_j)\in I_2\}}h(d(x,y_j)) g(y_j)\,.$$

Аналогично можно поступить для вычисления сумм $E_{ab}(\Phi)$ в остальных случаях с $a\ne b$.

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

(продолжение следует)

[vpb]

Теперь мы изложенное выше обобщим и улучшим.
Пусть $s\geq2$. Разобьем отрезок $I$ на $s$ частей $I_1,\ldots,I_s$, где $I_t=[(t-1)/s; t/s]$. Определим

$$ p_a(x)=\sum_{\{j\mid g(y_j)\in I_a\}} h(d(x,y_j))g(y_j)\,,\qquad q_a(y)=\sum_{\{i\mid f(x_i)\in I_a\}} h(d(x_i,y))f(x_i)\,, \qquad a=1,\ldots,s.$$

Также определим функции

$$ \widetilde p_1(x)=p_2(x)+p_3(x)+\ldots+p_s(x)\,, \qquad \widetilde p_2(x)=-p_1(x)+p_3(x)+\ldots+p_s(x)\,, $$

$$ \widetilde p_3(x)=-p_1(x)-p_2(x)+p_4(x)+\ldots+p_s(x)\,. $$

короче,

$$ \widetilde p_a(x)=\sum_{l=1}^s{\rm sgn}(l-a) p_l(x)\,, $$

где ${\rm sgn}(m)=-1$ при $m<0$, $1$ при $m>0$, $0$ при $m=0$.
Аналогично положим

$$ \widetilde q_a(y)=\sum_{l=1}^s{\rm sgn}(l-a) q_l(y)\,.$$

Сумму для $E(\Phi)$ разобьем на два слагаемых: $E(\Phi)=E_1(\Phi)+E_2(\Phi)$, где $E_1(\Phi)$ --- сумма всех слагаемых, отвечающих тем парам $(i,j)$, для которых $f(x_i)$ и $g(y_j)$ лежат в разных отрезках $I_a$, а $E_2(\Phi)$ --- тем парам, дляя которых они лежат в одном $I_a$. Понятно, что если значения $f(x_i)$ и $g(y_j)$ распределены по отрезку
$I$ примерно равномерно, то во вторую сумму попадает примерно $1/s$-я доля всех слагаемых, а все остальные (т.е. почти все, при $s>>1$) --- в первую.

Сумму $E_2(\Phi)$ вычисляем непосредственно.

Затем рассмотрим вычисление суммы $E_1(\Phi)$. Для $1\leq i\leq M$ пусть $a_i$ --- номер того из отрезков $I_a$, в который попадает $f(x_i)$. Аналогично определим $b_j$ тем, что $g(y_j)\in I_{b_j}$.

Выпишем формулу для $E_1(\Phi)$ (обдумать и проверить ее правильность оставляю Вам самому):

$$ E_1(\Phi)=\sum_{i=1}^N\widetilde p_{a_i}(\Phi^{-1}x_i)\ +\ \sum_{j=1}^M\widetilde q_{b_j}(\Phi y_j). $$

Применение этой формулы даже для небольших $s$ дает выигрыш в скорости вычисления функции $E(\Phi)$.

Теперь рассмотрим вопрос о том, каковы вычислительные ресурсы (время и память), нужные для вычисления $E_1(\Phi)$. Мы покажем, что время, по существу, примерно то же, что и при вычислении по (неверной) формуле

$$ E(\Phi)=\sum_i p(\Phi^{-1}x_i) -\sum_j q(\Phi y_j)\, $$

а памяти надо в $s$ раз больше.

Для простоты оценок будем считать, что $M=N$.

Пусть $L$ --- число узлов сетки, в которых табулируются $\widetilde p_a$ и $\widetilde q_b$. При вычислении по старой формуле табуляция функций $p$ и $q$ требует, очевидно, $O(ML)$ арифметических операций (вычисление каждого слагаемого $h(d(x,y_j))q(y_j)$ требует, очевидно, не более $C$ арифметических операций, где $C$ --- некоторая не очень большая постоянная (порядка нескольких десятков).) Объем же памяти для хранения значений $p$ и $q$ есть $O(L)$, очевидно (надо хранить значения $p$ и $q$
во всех узлах).

При вычислении "по новому" памяти надо в $s$ раз больше, так как нам надо хранить значения в тех же узлах $2s$ функций $\widetilde p_a$, $\widetilde q_b$.

Собственно вычисление распадается на две части:
(1) табулирование значений функций $\widetilde p_a$, $\widetilde q_b$ и
(2) применение найденной таблицы к вычислению функции $E_1(\Phi)$.
(Есть еще (2а) --- вычисление части $E_2(\Phi)$, которая рассматривается несколько позже.)

Часть (2) имеет ту же трудоемкость, что и в "неправильной" ситуации (это ясно, более-менее).

Для вычисления значений функций $p_a(x)$ в узлах мы вычисляем те же самые слагаемые $h(d(x,y_j))g(y_j)$, что и раньше. Однако, теперь эти слагаемые дают вклады в разные функции $p_a$. Но всего слагаемых столько же, сколько и ранее, а именно порядка $ML$. А стало быть, и сложность
та же. После того, как вычислены все числа $p_a(x)$, мы вычисляем все $\widetilde p_a(x)$. При этом, заметим, не обязательно каждый раз (для каждого $a$) вычислять сумму из $s-1$ слагаемого. Мы можем вычислить $\widetilde p_1(x)$, это требует $O(Ls)$ операций, а затем вычислить $\widetilde p_a(x)$ последовательно, используя соотношение $\widetilde p_a(x)=\widetilde p_{a-1}(x) -p_a(x)-p_{a-1}(x)$. Это опять потребует $O(Ls)$ операций. Таким образом, общая сложность табулирования функций $\widetilde p_a$ и $\widetilde q_b$ есть $O(ML+Ls)$. Что касается вычисления слагаемых $E_2(\Phi)$, то оно для каждого $\Phi$ имеет сложность порядка $M^2/s$ операций.

(Изложенное рекомендуется продумать и проверить. У меня, знаете ли, тоже бывают ошибки.)

[vpb]

vpb в сообщении #1362834 писал(а):

Но, собственно, к задаче то это не имеет отношения ! Это очень напоминает то, как школьник (причем бестолковый), "решая" задачу, комбинирует данные в условии числа каким-то случайным образом, типа складывает метры с килограммами. В общем, подозреваю, что Вы в логике задачи путаетесь.

В этом месте я был неправ (не въехал поначалу).

[_genius_]

Никакой научной тайны здесь нет, буду даже рад если кому-нибудь когда-нибудь наша дискуссия поможет. Поэтому смело переношу сообщение из ЛС в тему:

0) Если

$$ \sum_{j=1}^Mqq(\Phi y_j)+\sum_{i=1}^Npp(\Phi^{-1}x_i)-\sum_{j=1}^Mg(y_j)q(\Phi y_j)-\sum_{i=1}^Nf(x_i)p(\Phi^{-1}x_i) = \sum \sum h(d(x_i, \Phi y_j)) (f(x_i)-g(y_j))^2 $$

то в целом я не так уж и не прав. Так ли это, как Вы считаете? Если так, то быть может мне лучше использовать именно квадрат вместо модуля, т.к. я уже убедился что считается он довольно быстро.
1) Это нестрашно, меня это ни разу не смутило.
2) Все правильно, я собственно так и делал при оптимизации (т.е. менял знак), но не из каких-то принципиальных соображений, а просто потому что в библиотеке scipy есть только minimize.
3) Как Вам будет угодно.
4) С $h$ я ещё поиграюсь, вне всякого сомнения. Точек очень и очень много, поведение оптимизитора непредсказуемо, поэтому надо просто пробовать разные варианты.
Я в любом случае попробую Вашу новую идею для сравнения, но мне кажется что она может быть значительно медленнее, т.к. вычислений будет больше плюс дополнительные проверки кто больше, особенно если $f$ и $g$ лежат в одном интервале. Но насколько медленнее можно узнать только протестировав.