Аппроксимация по абсолютным разностям. Проще МНК.

[ewert]

Алексей К. в сообщении #334549 писал(а):

Я бы искал первое приближение,

Да вся беда в том, что поди его найди, это первое. Я бы, не заморачиваясь, безмятежно искал бы минимум суммы квадратов от разностей квадратов же расстояний методом градиентного спуска. Периодически пытаясь перейти к методу Ньютона (конечноразностному, конечно -- ну нафик, ещё и честные частные производные выписывать). А там уж -- как повезёт.

Tur в сообщении #334546 писал(а):

, а может ее и нет, а есть два горба у функции

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

[Tur]

ewert, Вы правы, вот график

Не могли бы Вы хотя бы кратко пояснить алгоритм в данном случае градиентного спуска. Я конечно о задаче со сферой.
Алексей К.,

Цитата:

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

К сожалению не понимаю. Вот исходная формула
$r_1 = \sqrt{[(x_0-x_1)n_1-(y_0-y_1)m_1]^2+[(y_0-y_1)p_1-(z_0-z_1)n_1]^2+[(z_0-z_1)m_1-(x_0-x_1)p_1]^2}$
Обозначим
$A_1(x_0,y_0) = (x_0-x_1)n_1-(y_0-y_1)m_1$
$B_1(y_0,z_0) = (y_0-y_1)p_1-(z_0-z_1)n_1$
$C_1(x_0,z_0) = (z_0-z_1)m_1-(x_0-x_1)p_1$
тогда $r_1 = \sqrt{A_1^2+B_1^2+C_1^2}$

Исходная сумма $S(x_0,y_0,z_0,R) = \sum\limits_i^n(r_i^2-R^2)^2$. Дифференцируем ее по $x_0$
$dS/dx_0=2(A_1^2+B_1^2+C_1^2-R^2)(2A_1n_1-2C_1p_1)+...+2(A_n^2+B_n^2+C_n^2-R^2)(2A_nn_n-2C_np_1)$
$dS/dx_0=4[(A_1^2+B_1^2+C_1^2-R^2)(A_1n_1-C_1p_1)+...+(A_n^2+B_n^2+C_n^2-R^2)(A_nn_n-C_np_n)]$
(почему происходит перенос на новую строку?)
Что Вы предлагаете дальше? Это уравнение третьей степени и таких будет четыре уравнения. Объясните мне и дальше я сам попробую все это развернуть.

[Алексей К.]

Я пока объясню только про простенькое, про перенос. Это просто считается, что формула шире ширины страницы.
Заключайте длинную формулу в двойные доллары: $$ ... $$.

А по делу --- вечером. Т.е. Вы решились на эти итерации-линеаризации?

[Tur]

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

[Алексей К.]

Во-первых, на всякий случай отмечу, что я в этих делах не профессионал. Занимался когда-то, лет 20 назад, подобными задачками. Тогда форума dxdy.ru ещё не было, и всё приходилось делать самому. Тогда я занимался допусковым контролем, и МНК был основным инструметном. Т.е. отфрезеровали --- отМНКачили. Обточили --- отМНКачили. Сейчас я как бы строю дороги, и ничего не приходится минимизировать. Обходимся преобразованием Мёбиуса.

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

На форуме научили интегрировать по частям. Я воспользуюсь этим методом, и буду отвечать по частям. Отправлю для начала эту преамбулу, и приступлю к объяснению линеаризации. Блин, там формулы ТеХом наколачивать, как я этого не люблю! Но ещё больше не люблю, когда не ТеХом.

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

[ewert]

Алексей К. в сообщении #335211 писал(а):

Ещё надо будет обсудить проблему первого приближения.

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

Метод же градиентного спуска -- при всей своей медленности гарантированно надёжен и гарантированно спускает нас к точкам, начиная с которых "линеаризация" будет уже разумна. Да и что значит "медленно"; четыре переменных -- тьфу.

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

[Алексей К.]

Упростим малость обозначения. Пусть искомый центр сферы называется $X,Y,Z$, радиус --- $R$. Индексы типа $X_0$ будем пользовать для очередного приближения. Очередной итерации. Их будет немного.
Ваша формула из сообщения #333979 примет вид: расстояние от $i$-той прямой, проходящей через точку $M_i (x_i, y_i, z_i)$ с направляющим нормализованным вектором $(m_i, n_i, p_i)$ до точки $M_0(X, Y, Z)$ определяется так

$$r_i = \sqrt{[(X-x_i)n_i-(Y-y_i)m_i]^2+[(Y-y_i)p_i-(Z-z_i)n_i]^2+[(Z-z_i)m_i-(X-x_i)p_i]^2}$$

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

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

(Оффтоп)

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

-- 25 июн 2010, 23:36 --

ewert в сообщении #335217 писал(а):

Метод же градиентного спуска -- при всей своей медленности гарантированно надёжен и гарантированно спускает нас к точкам, начиная с которых "линеаризация" будет уже разумна.

ewert, Корна я оставил на работе (хи-квадраты там завелись, я их тоже эабыл), интернетов не люблю. Т.е. не смогу где-то выучить про спуск, и научить автора. Сам, наверное, читает-изучает. А я буду пилить о том, что знаю.

ewert в сообщении #335217 писал(а):

Понимаю, что лень -- это нхорошо, но что поделаешь)

Если б только лень. А жара какая! gris, заразка, из речки небось не вылезает.

[Алексей К.]

Итак, мы хотим минимизировать функцю $F(X,Y,Z,R)=\sum(r_i-R)^2$ (суммирование везде от i=1 до N). Предлагался также вариант $\sum(r_i^2-R^2)^2$, который как бы избавляет нас от радикала, и очень здорово срабатывает в случае поиска наилучшей МНК-окружности (и, кстати, сферы, но не с заданными прямыми, а с заданными отфрезерованными точками). Поскольку мы собираемся линеаризовывать, нам до лампочки этот радикал, и можно основываться на первом варианте, $F(X,Y,Z,R)=\sum[r_i(X,Y,Z)-R]^2$. Кабы это $r_i(X,Y,Z)$ было линейно по искомым переменным, задачка бы легко решилась. Но оно, сволочь, у нас нелинейно (в обоих вариантах). Поэтому мы делаем его как бы линейным. А именно, если нам как-то удалось найти начальное приближение $X_0,\ldots,R_0)$, даже весьма грубое, то мы можем завести функцию

$$g(\Delta X,\Delta Y,\Delta Z,\Delta R)= r(X_0,Y_0,Z_0) +r'_X(X_0,Y_0,Z_0) \Delta X + r'_Y(X_0,Y_0,Z_0) \Delta Y+ r'_Z(X_0,Y_0,Z_0) \Delta Z-(R_0+\Delta R).$$

Теперь наши новые неизвестные --- эти самые дэльты.
Именно по ним всё линейно. Именно по ним мы берём производные $\dfrac{\partial g}{\partial(\Delta X)}$ итд. Может даже удобнее для них завести специальные буковки, например, $u,v,w$ вместо $\Delta X,\Delta Y,\Delta Z$. А коэффициенты типа ${{r_i^\prime}}_{Y}(X_0,Y_0,Z_0)$ известны, т.к. зависят только от известных величин (параметров i-той прямой и выбранного приближения). Мы решаем линейную задачу МНК, находим эти самые $\Delta X,\Delta Y,\Delta Z,\Delta R$, поправляем наше начальное приближение $X_1=X_0+\Delta X,\; Y_1=\ldots$, и с этим новым начальным приближением перерёшиваем задачу (следующая итерация). Если начальное приближение было приличным, то новые $\Delta X,\Delta Y,\Delta Z,\Delta R$ будут уже совсем маленькими. Итерации очень быстро сходятся.

Я понятно объяснил?

[Tur]

Понятно в первом приближении. Почему так можно делать? Как называется и где изложена эта теория? Решить бы подобную задачу с одной или двумя переменными чтобы все стало окончательно ясно. Какой геометрический смысл у $r'_X(X_0,Y_0,Z_0)\Delta X ?$
Начальное приближение я думаю что смогу получить. Эта задача фрагмент другой задачи. Касательные проведены только к сферическому углу градусов в 60. Причем точки касания лежат близко к двум дугам на сфере.

[Алексей К.]

Боюсь, разочарую Вас

Tur в сообщении #335451 писал(а):

Почему так можно делать? Как называется и где изложена эта теория?

Точно знаю, что нигде об этом не читал. Как бы изучал в ВУЗе математику, и это решение было естественным следствием. Ну а зачем я тогда функции в ряд разлагал, дифференциалы функции, в т.ч. нескольких переменных изучал? Как бы для того чтобы если потребуется --- применить. Даже не уверен, что это где-то описано. А с другой стороны --- слова "линеаризация функционала" --- неужели сам придумал? А почему бы и нет? (Кажется, ewert где-то повыше их закавычил.)
Почему так можно делать? Наверное, в студенчестве на контрольной доказал бы. Сейчас всё труднее составить доказательство того, что представляется очевидным.

Цитата:

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

Ну, это просто. Можно сгенерить точки для функции типа $e^{bx}$ и отфитировать. Можно решить какю-нибудь традиционную линейную задачку, типа провести МНК-прямую, но линейность проигнорировать. Как бы мы её не заметили.

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

Цитата:

Какой геометрический смысл у $r'_X(X_0,Y_0,Z_0)\Delta X ?$

Ну как бы естественные комноненты полного дифференциала функции нескольких переменных... Да и для одной: $y(x_0+\Delta x)\approx y(x_0)+y'_x(x_0)\Delta x$.

[Tur]

Вероятно это градиентный метод. Не нашел пока в сети простого и ясного примера. Вот простая программа на матлабе. Задаю пять точек близких к экспоненте, задаю нулевое приближение, шаг, вычисляю целевую функцию celfunc, изменяю b и снова считаю целевую функцию. Она уменьшается. Все верно пока? Но как узнать прибавлять db или вычитать?

function grad_exp_bx
e = 2.718281828459;
M1.X = 1; M1.Y = 1;
M2.X = 2; M2.Y = 3;
M3.X = 3; M3.Y = 4;
M4.X = 4; M4.Y = 8;
M5.X = 5; M5.Y = 12;
x = [M1.X M2.X M3.X M4.X M5.X]; 
y = [M1.Y M2.Y M3.Y M4.Y M5.Y];
axis([0 5.5 0 12]); grid on; hold on; plot(x, y, '.');

bo = 0.47; % нулевое приближение
yo = e.^(bo*x); 
plot(x, yo, '-g'); % зеленая линия нулевого приближения
celfunc= sum((yo - y).^2); % целевая функция
db = 0.01; % шаг
b = bo;
y1 = e.^(b*x) + b*e.^(b*x)*db;
celfunc1 = sum((y1 - y).^2);
b = bo + db;
y2 = e.^(b*x) + b*e.^(b*x)*db;
celfunc2 = sum((y2 - y).^2);
???????
end

[Алексей К.]

Я Вашего языка, пардон, не знаю (как, наверное, и Вы моего).
Вот, наколотил для проверки своих выкладок, для $e^{bx},\quad b=0.2,\quad b_0=0$.
Ноль --- ну очень плохое начальное приближение (сравните логарифмы), сходится на 6-й итерации.
Ну и можно там поправить цифирьки во первых строках.
Вы же умеете вьючить *.ps файлы?

1. %! 
2. /B  0.2 def 
3. /Bi 0.0 def   % Начальное приближение; можно изменить, например, на 0.1, -0.1 итп. 
4. /Niter 7 def   
5. /top 10 def 
6.  
7. /Exp {2.718281828 exch exp} bind def 
8. /Data [0 1 10 {dup exch B mul Exp [ 3 1 roll]} for] def  
9.  
10. 50 50 scale 1 1 translate /Helvetica findfont .3 scalefont setfont .02 setlinewidth 
11. -1 0 moveto 11 0 lineto 0 -1 moveto 0 11 lineto stroke 
12.  
13. Data {gsave 1 setlinecap .1 setlinewidth aload pop moveto 0 0 rlineto stroke grestore} forall 
14.   0 1 moveto .1 .1 10 {dup Bi mul Exp lineto} for stroke 
15.  
16. 0 1 Niter {/I exch def 
17.   0 1 moveto .1 .1 10.5 {dup Bi mul Exp dup 3 1 roll lineto 10 gt {exit} if} for  
18.   I 2 string cvs -.1 .1 rmoveto show stroke 
19.   .5 top moveto (Iter )show I 2 string cvs show %.2 0 rmoveto  
20.   (,  B=) show Bi 12 string cvs show 
21.  
22.   /top top .5 sub def 
23.   0 0 
24.   Data {aload pop /y exch def dup /x exch def Bi mul Exp /ebx exch def 
25.         x ebx mul dup mul add exch 
26.         x ebx mul y ebx sub mul add exch} forall 
27.   div Bi add /Bi exch def 
28.  
29.   I 0 eq {1 0 0}{currentrgbcolor 3 1 roll} ifelse setrgbcolor   
30. } for 

Сами же формулки выпишу (или впишу) чуть попозже.

-- 28 июн 2010, 15:23 --

Tur в сообщении #335732 писал(а):

Но как узнать прибавлять db или вычитать?

Прибавлять, конечно. Мы же ищем db из условия, чтобы новое $b_{i+1}=b_i\text{\color{blue}\bf\Large+}db$ было не хуже старого. Если ему понадобится отняться, оно само сделается отрицательным.

[Алексей К.]

Hack attempt!

[Tur]

Вот это хорошо, сразу полегчало. В первом уравнении производная по $b_0$ умножается на шаг $\Delta b$
Попробую все это запрограммировать чуть позднее сегодня надеюсь.

Вот нашел книжку http://matlab.exponenta.ru/optimiz/book_2/index.php
Там в методах можно еще переменный шаг сделать или вторую производную подключить.
Может быть еще одну задачу разберу с двумя переменными.

[Maslov]

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

Можно, например, вот здесь почитать: Гилл, Мюррэй, Райт. Практическая оптимизация (Мир, 1985)

В матлабе есть функция LSQNONLIN (Optimization Toolbox) для поиска решения нелинейной задачи МНК методами Гаусса-Ньютона и Левенберга-Марквардта.