Максиминное расположение точек на сфере

[B@R5uk]

мат-ламер в сообщении #1681600 писал(а):

В исходной постановке ведь речь идёт о системе точек на сфере. И ни о каком многограннике не говорится. В картинках, что я вижу у вас и у Haars много граней - четырёхугольники и пятиугольники. И не факт, что они всегда плоские.

А как ещё визуализировать топологию связей? Табличка длинных десятичных дробей N на 3 с координатами точек не очень информативна. Граф связности вершин стал новым для меня понятием (после открытия этой темы), над ним надо будет ещё поработать. А многогранники — это привычное дело (для меня во всяком случае, и процедура для них была давно написана).

мат-ламер в сообщении #1681600 писал(а):

Что касается проверки решения. Можно сделать несколько просчётов с разной начальной конфигурации.

Мультистарт — это вариант поиска новых решений "для отчаявшихся" (как и метод Монте-Карло). В зависимости от алгоритма оптимизации, вопрос об оптимальности находимого в каждом случае решения остаётся открытым.

[Евгений Машеров]

B@R5uk в сообщении #1681413 писал(а):

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

А не получится, что оптимальные решения для длин и углов разнятся?

[B@R5uk]

Евгений Машеров в сообщении #1681631 писал(а):

оптимальные решения для длин и углов разнятся?

Так ведь эти величины связаны между собой монотонной функцией:

B@R5uk в сообщении #1681413 писал(а):

Длина ребра и угол связаны формулой:

$$l=2\sin\frac{\alpha}{2}$$

Поскольку целевая функция является минимумом набора величин, то в экстремуме все минимальные величины будут равны и совершенно не важно, если их пересчитать по нелинейной формуле: равенство от этого не нарушится, экстремальность — тоже, лишь бы функция была монотонной. Вот если бы целевой функцией была какая-нибудь взвешенная сумма модулей или квадратов, тогда да, нелинейность бы всё испортила.

В этом плане Харс в своей статье пошёл ещё дальше: он минимизирует максимум скалярных произведений векторов. И это же гениально! Поскольку вектора единичные (указывают точку на единичной сфере) скалярное произведение в точности равно косинусу оптимизируемых углов:

$$\Bigl\langle\mathbf{R}\bigl(P,\;m\bigr),\;\mathbf{R}\bigl(P,\;n\bigr)\Bigr\rangle=\cos\alpha\bigl(P,\;m,\;n\bigr)$$

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

$$l_k=\Bigl\lVert\mathbf{R}\bigl(P,\;m_k\bigr)-\mathbf{R}\bigl(P,\;n_k\bigr)\Bigr\rVert^2=\bigl(x_{m_k}-x_{n_k}\bigr)^2+\bigl(y_{m_k}-y_{n_k}\bigr)^2+\bigl(z_{m_k}-z_{n_k}\bigr)^2$$

$$dl_k=2\bigl(x_{m_k}-x_{n_k}\bigr)\bigl(dx_{m_k}-dx_{n_k}\bigr)+2\bigl(y_{m_k}-y_{n_k}\bigr)\bigl(dy_{m_k}-dy_{n_k}\bigr)+2\bigl(z_{m_k}-z_{n_k}\bigr)\bigl(dz_{m_k}-dz_{n_k}\bigr)$$

будет

$$l'_k=-\Bigl\langle\mathbf{R}\bigl(P,\;m_k\bigr),\;\mathbf{R}\bigl(P,\;n_k\bigr)\Bigr\rangle=-x_{m_k}x_{n_k}-y_{m_k}y_{n_k}-z_{m_k}z_{n_k}$$

$$dl'_k=-x_{n_k}dx_{m_k}-x_{m_k}dx_{n_k}-y_{n_k}dy_{m_k}-y_{m_k}dy_{n_k}-z_{n_k}dz_{m_k}-z_{m_k}dz_{n_k}$$

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

С другой стороны, Харс в своей работе для оптимизации использует некий алгоритм LD_SLSQP (Sequential Least-Squares Quadratic Programming) из библиотеки NLopt library, предоставленной третьими лицами, и даже не потрудился вникнуть и расписать основные принципы его работы и трудоёмкость. При этом он растекается водой по дереву, обсуждая то, как ускорить и вообще улучшить процесс оптимизации, совершенно не привязывая свои действия к особенностям внутренней работы этого алгоритма оптимизации. (Не говоря даже о том, чтобы разработать свой метод под задачу). Читать это всё было подобно тому, как прийти на конференцию по проектированию автомобилей и попасть на лекцию про то, как и какой сабуфер лучше прикрутить в Жигули.

Некоторые решения автора в плане улучшения быстродействия очень сомнительны. Например, он убивает вращение всех точек как целого, фиксируя одну точку в северном полюсе сферы, а другую — в плоскости OXZ. Это приводит к тому, что когда алгоритм пытается оптимизировать рёбра, выходящие из этих двух вершин, то ему приходится вращать как целое систему всех остальных точек, вместо того, чтобы просто сдвинуть эти две фиксированные точки. У меня нет строго доказательства, что это плохо, но моя выработанная на практике интуиция говорит, что это резко ухудшает сходимость. Это равносильно искусственному внесению в целевую функцию особенностей типа функции Розенброка (может не с коэффициентом 100, а по-меньше на порядок, но тем не менее), что никак не может способствовать улучшению сходимости.

Однако, работа, проделанная автором по вычислению, оформлению и визуализации результатов, достойна похвалы. мат-ламер, ещё раз спасибо за ссылку.

[B@R5uk]

zgemm в сообщении #1681359 писал(а):

Исходим из того, что каждая грань икосаэдра заполняется равносторонними треугольниками, а потом полученную фигуру можно раздуть до того, что она ляжет на сферу. Назовем это дополненно-раздудый икосаэдр. Для заданного N ищем ближайший сверху такой дополненно-раздутый, и ищем на сколько он отличается от того, что надо нам. Далее находим две параллели близкие по значению одну в северном, другую в южном полушариях, на которых в сумме будет столько точек, как если отнять числа точек в этом раздуто-дополненного икосаэдре значение N. Удаляем точки на этих двух параллелях. Разрешаем каждой точке двигаться только вдоль ее меридиана, и пододвигаем так, чтобы дырки от удаленных точек схлопнулись. Запускаем минимизацию, уже без отжига обычным BFGSом.

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

$$\theta=\arccos\left(1-\frac{2n-1}{N}\right),\qquad\varphi=\pi\left(\sqrt{5}\pm1\right)n,\qquad 1\le n\le N$$

$$\begin{cases}x=\sin\theta\cos\varphi\\y=\sin\theta\sin\varphi\\z=\cos\theta\end{cases}$$

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

$$r=\sqrt{n/N},\qquad\qquad\varphi=\pi\left(1+\sqrt{5}\right)n,\qquad 1\le n\le N$$

$$\begin{cases}x=r\cos\varphi\\y=r\sin\varphi\end{cases}$$

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

[мат-ламер]

B@R5uk в сообщении #1681694 писал(а):

забавный вариант начального равномерного набрасывания точек на сферу

Ещё один забавный вариант равномерного набрасывания. В куб вписываем шар. Равномерно заполняем куб точками. Выбрасываем точки, которые не попали в шар. Проектируем точки на сферу - границу шара.

[B@R5uk]

Интересно, что в оптимальном расположении валентность точки (количество ближайших соседей) не может превышать пяти. Происходит это потому, что сферический равносторонний треугольник (в отличие от плоского) имеет угол при вершинах больший 60°. Это легко показать из сферической теоремы косинусов:

$$\cos a=\cos b\cos c+\sin b\sin c\cos A$$

$$\cos A=\frac{\cos a(1-\cos a)}{\sin^2a}=\frac{\cos a}{2\cos^2(a/2)}=\frac{\cos a}{1+\cos a}<\frac{1}{2}$$

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

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

Второй аргумент связан с тем, что рано или поздно минимальные длины рёбер приблизятся к оптимальным и надо будет уже задумываться о том, а не оказались ли мы вблизи локального минимума с необходимой точностью и пора остановить расчёт. Здесь надо будет проверять переполненную систему на возможность положительного решения, что является задачей линейного программирования. И только при отсутствии такого решения можно заключить, что локальный минимум достигнут. Тут ещё, кроме всего прочего, есть большой вопрос о том, как учитывать желаемую точность в процессе решения этой задачи ЛП. И решать её надо, потому что в качестве ответа кроме двух вариантов "достигли минимум" и "есть направление спуска" задача ЛП даёт ещё ответ "находимся в седловой точке, направление неизменности значения целевой функции такое-то". Именно этот третий вариант будет особенно чувствителен к погрешностям положений точек. Но правильный учёт этого случая обязателен, так как есть возможность выйти из седловой точки и продолжить улучшение значения целевой функции.

Могут быть и другие нюансы с выбором рёбер на этапе оптимизации, когда их длина заметно разнится.

[zgemm]

B@R5uk в сообщении #1681694 писал(а):

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

согласен, что всякие спирали проще программировать, только вот спираль потом плохо сходится и надо отжигом дожимать, а мой метод уже отжигом не дожимался, хотя я не спорю, что отжига тоже можно не дождаться. В моем методе почти все точки уже изначально лежат в вершинах почти равностороннего треугольника. Мой метод не работает, если надо точное число точек, а оно не удовлетворяет формуле $5k+2$, где $k$ - целое. Для маленького числа точек около двадцати и меньше я этот метод вообще не использовал, мне всегда тысячи и миллионы точек нужны были.

-- 10.04.2025, 23:49 --

Про минимизацию и выбор ребер. Вы Делоне стройте, и по ним выбирайте те ребра, что получились. Сложно программировать, но в общем тоже линейно, если предисторию помнить.

[B@R5uk]

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

$$2N-2$$

В то же время каждое ребро соединяет две точки, поэтому количество связей на точку равно половине её валентности. Это означает, что в среднем в ФФФ из вершины выходит 4 ребра, не смотря на то, что для фиксации положения точки требуется всего 3 ребра.

В связи с этим, я понял, что в очередной раз лоханулся, когда выдвинул гипотезу о том, что куб является седловой точкой:

B@R5uk в сообщении #1681597 писал(а):

Расположение куб для 8 точек является седловой точкой (я в лоб ещё не проверял, но скорее всего "градиент" целевой функции в ней равен нулю).

Куб имеет 12 рёбер, что на 2 ребра меньше, чем требуемые 14 для получения ФФФ из 8 точек. Дифферециал dL компонент целевой функции в формуле

$$\left(\begin{matrix}Q_R\bigl(P\bigr)\\Q_S\bigl(P\bigr)\\Q_E\bigl(P,\;H\bigr)\end{matrix}\right)dP=\left(\begin{matrix}0_{N+3}\\dL\end{matrix}\right),\qquad dL\ge 0$$

можно сделать любым (допустимым) и она будет разрешима с по крайней мере одной степенью свободы (возможно больше), так как матрица слева недоопределена (24 столбца, 23 строки).

Встаёт вопрос: существует ли (и как выглядит, если да) расположение точек, которое является седловым для задачи Таммеса? Под седловым я имею в виду такое расположение P, что уравнение выше разрешимо относительно дифференциала dP координат, но только при нулевом дифференциале dL. А то может моё беспокойство о седловых точках в этой задаче вообще не обосновано.

мат-ламер в сообщении #1681695 писал(а):

Равномерно заполняем куб точками. Выбрасываем точки, которые не попали в шар. Проектируем точки на сферу - границу шара.

Не понятно, что именно вы имеете в виду под выделенными словами. Если случайным образом — то это не годится, потому что я имел в виду детерминированный алгоритм. Если же в вершинах кристалла простой кубической сингонии, то будут существовать точки с пропорциональными координатами, которые спроектируются в одну точку сферы. Это противоречит требованию равномерного распределения по площади сферы. Если две точки могут совпасть, то ничто не запрещает им оказаться сколь угодно близко. Можно, конечно, попытаться починить это выбором некоторого достаточно тонкого сферического слоя, из которого будет делаться проекция (вместо всего шара), но это целый дополнительный подгонный параметр, в то время как результируемое количество точек и так весьма неоднозначно связано с величиной шага решётки. Эти проблемы остаются (и усложняются), если попытаться добавить смещение кристаллической решётки точек относительно центра сферы.

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

[мат-ламер]

B@R5uk в сообщении #1681778 писал(а):

Не понятно, что именно вы имеете в виду под выделенными словами. Если случайным образом — то это не годится, потому что я имел в виду детерминированный алгоритм.

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

[B@R5uk]

Опять я ерунду выше написал, и что-то никто меня не поправит. При вращении одной грани куба относительно противоположной движение точек происходит в точности перпендикулярно удлиняемым ребрам. Это значит, что величина их удлинения имеет второй порядок малости по величине смещения. Это значит, что первые производные длины по смещению нулевые. В рассуждении выше у меня изъян:

B@R5uk в сообщении #1681778 писал(а):

Дифферециал dL компонент целевой функции в формуле

$$\left(\begin{matrix}Q_R\bigl(P\bigr)\\Q_S\bigl(P\bigr)\\Q_E\bigl(P,\;H\bigr)\end{matrix}\right)dP=\left(\begin{matrix}0_{N+3}\\dL\end{matrix}\right),\qquad dL\ge 0$$

можно сделать любым

Я его выделил. Несмотря на то, что матрица системы недоопределена (число строк меньше числа столбцов), это ещё не значит, что у неё сами строки не могут быть линейно зависимыми (ранг матрицы меньше числа строк). В этом случае компоненты dL будут друг от друга зависеть, то есть его нельзя взять произвольным.

С кубом именно так и происходит. Я даже программку составил для того, чтобы в этом убедиться на практике (ниже). Матрица для расчёта dL после перехода к независимым дифференциальным координатам имеет одну линейно зависимую строку. Более точнее, сумма элементов в каждом столбце равна нулю. Из этого следует, что то же самое выполнено и для dL, и если в dL есть хотя бы одна положительная координата, то обязательно найдётся и отрицательная. То есть нельзя из этого положения улучшить целевую функцию, однако есть направления (подпространство дифференциалов размерности 2), в котором целевая функция в первом порядке остаётся неизменной. Это и есть седловая точка для данной задачи.

Вообще, если внимательно подумать, то можно эту ситуацию обобщить на произвольный случай. Пусть некоторое подмножество точек образует жёсткую фигуру и соединяется с основным каркасом отрезками, продолжения которых пересекаются в одной точке (примеры на рисунке ниже). Это может быть и равносторонний треугольник, и пара, и что-нибудь по-сложнее. Тогда можно начать вращать это подмножество вершин вокруг точки пересечения продолжений. Движение внешних вершин подмножества будет перпендикулярно соответствующим отрезкам, соединяющим подмножество с внешним каркасом, поэтому они будут удлиняться, но это удлинение будет величиной второго порядка малости по величине поворота. Если конструкция в общем, за исключением этого поворота является жёсткой (никакие другие вершины нельзя сдвинуть, чтобы улучшить целевую функцию), то какое расположение точек является седловым для задачи Таммеса.

То есть седловая точка не является таким уж необычным явлением. Встаёт, правда, вопрос, на сколько вероятно, что алгоритм оптимизации на неё наткнётся. Тут есть два момента. Первый заключается в том, что конфигурации выше весьма надуманы. Интуитивно кажется, что чтобы к ним прийти, алгоритм должен обладать определёнными предпочтениями в плане выбора новых положений точек на каждом шаге. Но опять же, это гипотетически, а мои гипотезы меня в последнее время подводят. Второй момент связан с тем, что седловая точка требует от координат определённой "когерентности", если можно так сказать. Я не знаю, о всех случаях, когда нетривиальное решение матричного неравенства выше может давать в правой части только нули, но несомненно, что между координатами должна быть определённая взаимосвязь. Это так же следует из геометрического толкования: чтобы продолжения отрезков пересекались в одной точке, координаты соответствующих вершин должны быть взаимосвязаны. В компьютере числа представляются конечной точностью, поэтому эта взаимосвязь обязана выполнятся не идеально, а лишь в заданном отрезке. То есть объём области с выполнением приближённого равенства к объёму всего пространства параметров не является нулевым. Другое дело, что это отношение очень мало, и вероятность можно прикинуть как

$$p=p_0\left(\frac{\delta}{l}\right)^{K-2}$$

δ — это машинная точность, l — длина ребра и K — число рёбер, продолжения которых пересекаются. Matlab, например, для определения ранга матрицы по умолчанию использует величину

$$\delta=2^{-52}N\sigma$$

где N — размер наибольшей из сторон матрицы, а σ — её наибольшее сингулярное значение. На первый взгляд вероятность напороться пренебрежима, а с другой — я не уверен, что вместо длины ребра не надо подставить изменение этой длины на шаге алгоритма. Тогда вероятность перестанет казаться такой уж маленькой и встанет вопрос: а на каком минимальном удалении от локального экстремума можно ещё встретить седловую точку? Ну и хорошим тоном в программировании является учёт всех особых случаев.

[B@R5uk]

При анализе связей в задаче часто появляется необходимость решать матричное уравнение

$$AX=Y$$

Как правило оно однородное (все элементы Y равны нулю), но тем не менее. Тут возможны различные ситуации, связанные с тем, как соотносятся между собой ранк матрицы A, ранк расширенной матрицы, количество строк и количество столбцов матрицы A. Когда решение существует, его можно выразить в виде

$$X=BY+CZ$$

с некоторыми "волшебными" матрицами B и C и столбцом независимых (свободных) переменных Z. Последние два объекта могут отсутствовать, если исходная матрица квадратная и невырожденная или если она переопределённая (число строк больше числа столбцов), но система в целом совместна. Как учат в курсе линейной алгебры, от исходной системы к решению (с искомыми матрицами B и C) можно естественным образом перейти, если воспользоваться методом Гаусса (для матрицы A, дополненной справа единичной матрицей) и построить приведённую ступенчатую форму по строкам (reduced row echelon form, функция rref в Matlab). Матрица B при этом получается из единичной. Дальше свободные переменные переносятся в правую часть, где коэффициенты перед ними образуют матрицу C после "подстыковки" снизу единичной матрицы с правильным знаком. Её столбцы образуют базис решения, но они не являются в определённом смысле наилучшими. В идеале хотелось бы, чтобы он был ортонормированным. Предпочтительность такого базиса связана с тем, что его матрица очень хорошо обусловлена (лучше в принципе уже нельзя), что важно для точности численных вычислений, и вообще является полезным свойством. Можно, конечно, применить к результату ортогонализацию Грамма-Шмидта или ещё какой-нибудь алгоритм (функция orth в Matlab). Но существует более прямолинейный метод решения системы через сингулярное разложение матрицы (функция svd).

Меня сейчас будет интересовать случай, когда число строк M матрицы A меньше числа столбцов N (как следствие, она имеет нетривиальное ядро, с нахождением базиса в котором связано вычисление матрицы C). При этом либо уравнение совместно, либо однородно. Сингулярное разложение матрицы действительных чисел — это представление в виде

$$A=U\Sigma V^T$$

где матрицы U и V являются матрицами поворота (и, возможно, включают в поворот отражение, то есть их определитель может быть равен −1). Это означает, что они квадратные и ортогональные (с единичным по модулю определителем):

$$UU^T=U^TU=I,\qquad\qquad VV^T=V^TV=I$$

Их размер разный: размер U совпадает с числом строк M матрицы A, а размер V — с числом столбцов N. Размер матрицы Σ совпадает с размером A, она заполнена нулями, кроме главной диагонали, на которой стоят сингулярные значения матрицы A (неотрицательные числа) в убывающем порядке. Её рангом r является число этих значений, превышающих некоторый порог, который обычно принимается равным произведению первого сингулярного значения на малую величину, то есть величина δ из поста выше. В зависимости от задачи порог можно увеличить, но уменьшать его нет смысла, потому что он связан с конечностью машинной точности. В любом случае, с учётом разложения решаемое уравнение примет вид:

$$U\Sigma V^TX=Y$$

Перейдём к новым переменным:

$$\Sigma X'=Y',\qquad\qquad X=VX',\qquad\qquad Y'=U^TY$$

Матрица Σ имеет вид, где в верхнем левом углу стоит квадратная матрица с ненулевой главной диагональю, а все остальные элементы вокруг — нулевые:

$$\Sigma=\left(\begin{matrix}\Sigma_\square&0\\0&0\end{matrix}\right)$$

Уравнение совместно, если оно однородно или если на последних строках его матричной записи стоят только нули, в том числе в столбце Y'. То есть, если ранг r матрицы A (и матрицы Σ) меньше числа строк, то решение о совместности системы принимается на основе того, что последние элементы столбца Y' равны нулю с заданной/желаемой точностью, так как вычисления с плавающей точкой по своей природе неточные. В результате уравнение преобразуется

$$\Sigma X'=\left(\begin{matrix}\Sigma_\square&0\\0&0\end{matrix}\right)\left(\begin{matrix}X'_{top}\\Z\end{matrix}\right)=\left(\begin{matrix}Y'_{top}\\0\end{matrix}\right)=Y'$$

Нижние элементы столбца X' выбираются произвольно, это свободные переменные в решении (я сразу обозначаю их как столбец Z), а верхние — вычисляются по формуле:

$$X'_{top}=\Sigma_\square^{-1}Y'_{top}=\Sigma_\square^{-1}U_r^TY$$

нижний индекс r у матрицы U означает, что перед транспонированием берутся первые r столбцов этой матрицы, а остальные — отбрасываются. Здесь вычисление обратной матрицы Σ в силу её природы заключается лишь в обращении элементов главной диагонали. Столбец искомых переменный выражается так:

$$X=V\left(\begin{matrix}X'_{top}\\Z\end{matrix}\right)=V_rX'_{top}+V_{-r}Z$$

Здесь отрицательный индекс −r означает, что у матрицы V берутся только последние столбцы, которые дополняют её первые r до квадратной. Окончательно имеем:

$$X=V_r\Sigma_\square^{-1}U_r^TY+V_{-r}Z$$

$$X=BY+CZ,\qquad\qquad B=V_r\Sigma_\square^{-1}U_r^T,\qquad\qquad C=V_{-r}$$

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

$$-\delta<U_{-r}^TY<\delta$$

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

Надо заметить, что в программе Matlab функции null, orth и rank, имеющие непосредственно отношение к задаче, (а так же, возможно, некоторые другие, например norm) — все работают через функцию svd даже тогда, когда время выполнения можно было бы в 2 раза улучшить с другими алгоритмами. Это делается потому, что сингулярное разложение даёт более точные и стойкие результаты.

[B@R5uk]

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

(Оффтоп)

имеется 6 жёстких рёбер: 3 из них связывают точки между собой, а остальные 3 — каждую из них с внешними точками фигуры. Считая внешнюю фигуру неподвижной имеем 6 равенств и 7 неизвестных: 3 точки по 2 координаты и длину ребра. Система недоопределена, поэтому возможно движение этих трёх точек, увеличивающих длину ребра (оставляя их все равными между собой, но не с рёбрами остальной фигуры), до тех пор, пока одна или несколько других точек не упрутся новыми жёсткими рёбрами в какие-либо другие точки, увеличив таким образом число связей и сделав систему определённой (или переопределённой, но совместной).

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

С простым подсчётом рёбер и степеней свободы, правда, есть проблемы. Например, на картинке справа имеется 4 точки и 9 рёбер. То есть 9 связей и 9 неизвестных. Кажется, что система определена, но это не так: одна связь является следствием всех остальных. Мы это понимаем, потому что сморим на фигуру и видим, что она допускает вращение, улучшающее целевую функцию, но что бы компьютер это понял, скорее всего не обойтись без матричного анализа дифференциалов, но только не для всей фигуры, а для рассматриваемых точек. Конкретно для рисунка решение парадокса в том, что одна из связей, соединяющая внутреннюю точку с внешней, после поворота перестанет быть жёсткой. Которая именно, скорее всего, зависит от направления вращения.

Другое следствие заключается в том, что движение из седловой точки возможно в две стороны: по и против часовой стрелки. И как уже замечено, они могут привести к разным конечным конфигурациям. Это значит, что в фигуре, доставляющей задаче Таммеса локальный экстремум, можно искусственно выискивать подмножества точек, которые достаточно близки к седловой точке (если подмножество рассматривать рассматривать независимо, считая остальную фигуру фиксированной), а затем "перещёлкивать" их в противоположное от седловой точки положение. После релаксации всей конфигурации скорее всего окажемся в новом локальном экстремуме задачи. Он может оказаться как лучше, так и хуже предыдущего, но в любом случае, в нём тоже можно поискать подобного рода "дорогу" в ещё один соседний локальный экстремум. Изучая таким образом связи различных локально оптимальных решений между собой и переходя от одного к другому, шанс найти глобальный экстремум кажется более вероятным, чем рандомным мультистартом. И точно менее трудоёмко, хотя бы потому, что мы не будем рассматривать одни и те же сомнительные конфигурации бессчётное число раз. Близость локальных экстремумов так же уменьшает время вычисления положения новых с требуемой точностью.

Возможно, это и делается в статье китайцев, что я приводил ранее. Надо уже в конце-концов прочитать её и вникнуть в суть дела там.

[Alex Krylov]

Цитата:

Ещё один забавный вариант равномерного набрасывания.

Есть еще вариант получения рандомизированных решений через полуопределенную релаксацию (см. Semidefinte Programming итп) исходной невыпуклой оптимизационной задачи. Такие рандомизированне решения дают в среднем большее минимальное расстояние, нежели (абсолютно случайные) рандомизированные решения, полученные посредством генерирования двух случайных углов (см. сферическую систему координат). Однако, если затем этими рандомизированными решениями инициализировать локальный решатель (fmincon в Matlab), то никакого выигрыша нет. Получается примерно одно и то же.

[zgemm]

(Оффтоп)

B@R5uk в сообщении #1681842 писал(а):

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

странно, у меня это было в 1994 году, и, кстати, даже на экзамене попалось.

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

Если SVD численно слишком "тяжел", гуглите на ранг-ревеалинг, им тоже можно такие задачи решать не боясь, что где-то численная ошибка с неустойчивостью перехлестывается

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

[B@R5uk]

zgemm в сообщении #1682092 писал(а):

исходная матрица слегка не точно задана

О, она всегда неточно задана! Ну, кроме случаев, когда все элементы матрицы имеют вид

$$a=m\cdot2^n,\qquad 0\le m<2^{53},\qquad n\in\mathbb{Z}$$

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

zgemm в сообщении #1682092 писал(а):

что кривой алгоритм может дать численную ошибку существенно большую, чем длина мантиссы.

Я знаю про один такой алгоритм. Вычисление сингулярных чисел через квадратные корни из собственных значения ковариационной матрицы столбцов (или строк — в зависимости от того чего меньше) убьёт точность для малых чисел, потому что квадрату числа гораздо легче оказаться меньше величины машинной точности. Подобный расчёт для чисел, по величине сравнимых и меньших квадратного корня из машинной точности, не даст ни одной точной значащей цифры (они даже комплексными могут оказаться). Тем не менее, этот алгоритм находит свою нишу. При анализе данных зачастую не требуется считать полное разложение, достаточным является лишь небольшое количество первых (наибольших по величине) сингулярных значений. Если данные на столько велики, что исходная матрица не помещается в память целиком и невозможно использовать что-то более точное, то этот алгоритм является рабочей альтернативой (для расчёта ковариационной матрицы достаточно лишь подгрузить два вектора — это медленно, но не смертельно), и вполне себе используется даже сейчас (когда известны быстрые и более точные альтернативы).