Возникла необходимость аппроксимировать полученые данные эллипсом с минизацией квадратов геометрических (не алгебраических!) рсстояний от полученых точек до эллипса. Нашёл много быстрых реализаций для алгебраических расстояний, но нужно именно геометрические минимизировать (поскольку вектора X и Y имеют равную погрешность). В качестве затравки к алгоритму можно использовать результат минимизации алгебраических расстояний так что использование метода Ньютона-Гауса (или других, чувствительных к начальному приближению) вполне пригодно. Может кто-либо встречался с какой-нибудь статьёй на эту тему? или есть готовые библиотеки для таких манипуляций?

Ах как интересно! А в каком направлении вы собираетесь мерить (минимизируемые) расстояния в этой задаче? Или спрошу об этом иначе - как вы узнаете точку эллипса, которая соответствует конкретной точке данных?

В класическом методе наименьших квадрартов, например, когда прямую линию по экспериментальным точкам проводят, то отклонения измеряют по вертикали (ординате). Их и минимизируют в смысле наименьших квадратов. А когда эллипс, то как? Можно, по-прежнему, по вертикали, то тогда на краях плохо будет. Можно измерять разницу в полярных координатах, по вектору из начала координат. Но где тут это начало? В точке пересечения осей эллипса или в одном из его фокусов? Сплошные непонятки... А может быть из точки надо опускать перпендикуляр на поверхность эллипса? Но тогда очень уж сложно получится.

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

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

Что такое геометрическое расстояние (ну и алгебраическое заодно)?
В МНК минимизируют функцию правдоподобия.
Обычно правдоподобие сводят к минимизации квадрата Евклидовой нормы.
Квадрат выбирают:
1. Сильно выпуклая функция
2. Можно дифференцировать

Евклидова (L2) - можно дифференцировать.
Кроме того обычно именно в L2 решение существует (например в l1, l2, C и пр. его может просто не существовать для этой же задачи).

Норма - это длина вектора.
Метрика (расстояние между двумя точками) определяют через норму.
Все это чистая геометрия.
Дискретное представление эллипса можно представить как N_мерный вектор.
Искомое решение на этой же сетке - это другой N-мерный вектор.
Вы (с помощью МНК) минимизируете квадрат расстояния ||Zo-Zp||^2.
Что вы имеете в виду под геометрическим и алгебраическим расстоянием?

Про эллипс подумалось:
Если вам надо [a,b,xo,yo], то это в рамках МНК тоже вектор (размерность 4).

Я имел ввиду геометрическое расстояние в терминах аналитической геометрии, где геометрическим расстоянием от точки до кривой есть длина перпендикуляра от точки до кривой. То, что задача получается не очень простая мне и так ясно, поэтому спрашиваю о готовом её решении, а не решаю сам.


Спасибо за введение в аналитическую геометрию, однако это не является вопросом темы. И если уж придираться - то в МНК минимизируют квадрат некой нормы, которая может являться оптимальной либо не оптимальной в зависимости от условий задачи и выбор которой во многом определяет эффективность метода.
Что касается геометрического расстояния и алгебраического - термины взяты из зарубежной литературы на эту тему, например:
http://www.ulb.ac.be/assoc/bms/Bulletin/sup962/gander.pdf

Если подробнее про расстояния, то алгебраическое это когда норма, которую следует минимизировать является функцией (F(xi)-yi)^2 т.е. предполагается, что по одной из осей погрешность отсчета много меньше чем по другой. Однако, это не всегда верно. и чтобы уравнять оси в "правах" необходимо использовать целевую функцию в смысле геометрического расстояния: (xi-x0)^2+(yi-y0)^2. Для нахождения x0,y0 можно построить нормаль к эллипсу через точку xi,yi и найти её точку пересечения с кривой. Для определённости: в моих рассуждениях использую прямоугольные декартовы координаты. При первых оценках метода у меня всё свелось к уравнению 4й степени, однако, возможно, эту проблему можно обойти. Поэтому я и интересуюсь, может кто-то уже нашёл метод решения в литературе.
Заранее спасибо!

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

Алгебраическая постановка - это линеаризованная модель для коник.
У каждой линеаризации есть область где она работает (область определения).
Авторы распространили эту область на все пространство и удивляются хреновости....

Потом решили сделать как я уже написал. Сделали дискретную модель эллипса (через cos(t), sin(t)).
Каждое изменение параметров искомой модели создает сове уникальное множетсво.
И после этого к этому множеству точек применяют минимизацию.
Замкнутая задача является полной и нелинейной.
Нового тут уже ничего не придумать. Все проблемы нелинейности известны. Книг море.

Вам надо внести априорную информацию, всю какая у вас есть, чтобы сжать область поиска.
А желательно свести эту область к компакту.

Только в том, что его нужно решить аналитически, затем сделать подстановку в 5 уравнений МНК(по числу неизвестных параметров) , взять соответствующие производные, а затем выразить все параметры через суммы. Если говорить проще - казалось, что задача достаточно стандартна чтобы быть решённой и не тратить на неё много времени. А ещё не ошибиться в многостраничных выкладках...

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



Мне казалось более правильным аналитически свести задачу к системе нелинейных уравнений и далее численными методами за несколько итераций (имея хорошую затравку) спуститься к требуемой точности.

Как способ задания уравнений определяет скорость решения и точность (я не говорю про численное представление)?
Для нелинейных уравнений нельзя использовать метод Гаусса (приведение к LU виду).

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

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

У коник очень много свойств на вписанные и описанные окружности, треугольники и касательные.
Можно их по изучать, возможно вы знаете о задаче больше, чем сейчас считаете.
Дома посмотрю название книги где много свойств коник.

Перпендикуляр к касательной - биссектриса угла от точки касания к полюсам.
Легче?

Акопян, Заславский "Геометрические свойства кривых второго порядка"


Мне все равно, поэтому и так не тяжело.

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

Я говорил именно о численных решениях т.к. аналитических готовых я не нашёл в литературе. Система нелинейных уравнений решается методом Гаусса-Ньютона без особых проблем.
Кстати, в реальной жизни даже аналитическое решение, представленное в том или ином виде может иметь различные показатели по скорости и точности(!), поскольку, например, DSP могут иметь некоторые заложенные в них аппаратные примитивы и, опираясь на них, решение можно привести к более скоростному вычислению (тождественные преобразования, рекурсивные методы и т.д.) а ещё, в форматах чисел с плавающей запятой, может не выполняться свойство коммутативности: (a++c != a+(b+c) - соответственно и точность может быть разной в зависимости от последовательности действий.
Спасибо за помощь. Буду дальше копаться, если не удастся найти - придётся самому решать.

Задача для случая когда оси симметрии эллипса параллельны осям или расположены произвольно?
Аналитически решаются оба случая.

Интересная задача.
Но прежде чем пытаться искать методы оптимизации, стоит осмыслить:
- случаи, когда задача не имеет решений;
- случаи, когда задача имеет более одного решения.
Нетрудно придумать простейшие примеры, которые иллюстрируют по крайней мере, последний вариант.
Что будет делать алгоритм, если решений бесконечное множество?
Очевидно, следует задачу как-то дополнительно ограничить.

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

Это ассоциативность. Не знаю, как смайлик "плавающая запятая", но смайлик "ухмылки" явно не ассоциативен...

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

Я не виду основы для спора.
Для данных без шума тоже приходится искать решение, это наверное 90% вычислительных задач.
Забудьте про СКО. Минимум квадрата нормы (евклидовой) невязки. В нее может быть замешана мат статистика, а можно не мешать.
Смотрите на задачу с точки зрения функана.

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

Определиться с областью сходимости - внести АПРИОРНУЮ информацию.

Для двух точек нормальный алгоритм автоматически вернет вернет решение z=0, решение с нулевой нормой, которое по определению должно входить в область определения оператора.

Гарантия глобального минимума в нелинейной задаче - это предмет исследования для каждой задачи. Это зависит от решателя и типа задач.
Цитировать Остапа Бендера или просто классические работы по этому поводу нет смысла.

Самое простое по статье сделать алгоритм (все описано, кроме глобальной сходимости) и сверять его со своим самописным.
Алгоритм общий, проблемы его общие и понятные.
Если автор по другому решит задачу и захочет обсудить, то тогда продолжим.

Если есть вопросы, то продолжать можно в личке, так как это уже не тема ЦОС.

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

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

В моей задаче имеются точки на плоскости. часть из них хорошо описываются эллипсом, а другая часть (их может быть достаточно много, но меньше 30% от полезных) - шумовая, которые могут в свою очередь образовывать контуры, а могут быть хаотично рассыпаны.

Если по теме, то нашёл 2 метода:
1. пересчитывать точки в канонический базис текущего эллипса и в нём всё решать
http://www.geometrictools.com/Documentatio...aresFitting.pdf
2. Проводить из точки кратчайшее расстояние до эллипса и минимизировать квадрат этих расстояний.
http://ieeexplore.ieee.org/xpl/login.jsp?t...number%3D813057

К сожалениею, не имею полного доступа к материалам ссылки 2. Если у кого-нибудь есть - поделитесь плз. полным pdf.


Оси произвольно расположены.
Можно ссылочку на статью или книгу с аналитическим решением такого случая?

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

Если отбросить возможную обусловленность, то может быть "This problem is more dicult than that of fitting circles" будет простой на каждом шаге.
На выходных попробую проверить корректность перехода для данных с шумом.

Sci-Hub представляет по ссылке вашу статью, забирайте, если успеете

sci-hub.org/pdfcache/79cab939b0b85e2f48f3e2ae8692a600.pdf

Что касается самой нелинейной оптимизации - конечно гарантии сходимости к глобальному минимуму быть не может.
Но с очень большой вероятностью сходимость будет иметь место, если использовать не градиентные методы спуска,
а групповые по типу PSO или ещё более продвинутые "генетические" алгоритмы этого типа.
PSO совершенно элементарный алгоритм для реализации, но при грамотном задании области поиска 5-мерную задачу на современных компьютерах он успешно переберет с вероятностью почти 1.0
http://en.wikipedia.org/wiki/Particle_swarm_optimization

Успел! Спасибо большое!

Я даже не буду оспаривать математическую "ценность" подобных измышлений...

Упростите для начала задачу. Представьте, что данные раньше измерялись с погрешностями sigma_1 и sigma_2 для координат x_1 и x_2, а после того, как в измерительную систему вселился гремлин, x_1 измеряется с погрешностью 0, а x_2 с погрешностью sigma_1 + sigma_2. Задайте аналитически произвольный эллипс. Внесите погрешности (исказите) в случайный набор точек на эллипсе. Найдите параметры "измеренного" эллипса по МНК, считая, что измерительная система населена гремлином.

Если полученное решение не слишком далеко от истинного (задайте метрику), то пользуйтесь упрощённой моделью.

Если решение "далеко", увеличьте количество измерений. Если и это не помогает, то пересмотрите модель - считайте честно овал ошибки.

Статью наврятли смогу привести, но суть в следующем:

Вашу задачу давным давно я решал следующим образом:

Эллипс на плоскости аппроксимируется алгебраической кривой заданной в неявной функцией.
Используем для этого уравнение второго порядка.
F(x,y)=a11*x^2 + a12*xy+ ...+a33 =0 (6)
где: aij = aji - коэффициенты соответствующей квадратичной формы с матрицей A:
Подставим в (6) измеренные координаты (xk,yk) и получим
F(xk,yk)=a11*xk^2 + a12*xkyk+ ... = dk
В случае малости расстояния от точки до эллипса, получаемая невязка δk также будет мала, что
позволяет определить ее как расстояние от точки (xk,yk) до кривой . Минимизируем
невязки для n точек методом наименьших квадратов:
Ф(aij) = SUM(dk^2) -> 0
для чего решим систему из шести уравнений:
dФ/daij = 0 для i,j,=1..3
Получаемую систему удобно записать в матричном виде:
Xu=0 (7) ,
где X – симметричная матрица, равная:
и вектор искомых коэффициентов
u=[a11,....a33]
Получившуюся систему можно решить численным методом. Найдем собственные значения для матрицы X любым из известных методов [3]. Минимальному собственному значению будет соответствовать собственный вектор , который и будет являться решением системы (7), т.е представлять набор коэффициентов эллипса.

Попробуйте например, с помощью Матлаба, это легко.