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

[Tur]

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

Maslov, спасибо за функцию и книгу.

На каждый совет, как известно, требуется еще десять советов чтобы он не остался втуне.

[Александрович]

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

При минимизации суммы min(e1+e2+e3+e4) как получить искомую прямую (m, n)? Какая здесь теория?

Теория, а точнее алгоритмы нахождения коэффициентов изложены в "Метод наименьших модулей" Myдров В.И., Кушко В.Л.

[alekcey]

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

[b][color=#3333FF]

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

Каждая прямая (шесть чисел) дает три таких уравнения. Если всего четыре прямых...

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

[Tur]

Спасибо всем за участие.
Задачу со сферой я решил как советовал мне Алексей К.
Даны $N$ прямых в пространстве, близких к касательным к сфере. Каждая прямая дана точкой и вектором. Точки $(x_i, y_i, z_i)$, соответствующие вектора $(m_i, n_i, p_i)$, i = 1 – N.
Центр сферы $M(X,Y,Z)$, радиус $R$. $r_i$ - расстояние от i-той прямой до центра сферы. (формула 1)

$$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}$$

Требуется минимизировать функцию (формула 2) $F(X,Y,Z,R)=\sum[r_i(X,Y,Z)-R]^2$
Вычислим частные производные( формулы 3)

$$\dfrac{\partial r_i}{\partial X}=\dfrac{X(n_i^2+p_i^2)-m_in_iY-m_ip_iZ-(n_i^2+p_i^2)x_i+m_in_iy_i+m_ip_iz_i}{r_i}$$

$$\dfrac{\partial r_i}{\partial Y}=\dfrac{-m_in_iX+(m_i^2+p_i^2)Y-n_ip_iZ+m_in_ix_i-(m_i^2+p_i^2)y_i+n_ip_iz_i}{r_i}$$

$$\dfrac{\partial r_i}{\partial Z}=\dfrac{-m_ip_iX-n_ip_iY+(m_i^2+n_i^2)Z+m_ip_ix_i+n_ip_iy_i-(m_i^2+n_i^2)z_i}{r_i}$$

Пусть нулевое приближение $M_0(X_0,Y_0,Z_0), R_0$. Заменим функцию ( 2 ) на новую функцию ( 4 ) линейную относительно переменных $\Delta X, \Delta Y, \Delta Z, \Delta R$

$$F(\Delta X,\Delta Y,\Delta Z,\Delta R)=\sum[r_i(X_0,Y_0,Z_0)+r_i'_X(X_0,Y_0,Z_0)\Delta X+r_i'_Y(X_0,Y_0,Z_0)\Delta Y+r_i'_Z(X_0,Y_0,Z_0)\Delta Z-(R_0+\Delta R)]^2$$

где (формулы 5)

$$r_i(0)=r_i(X_0,Y_0,Z_0) \ r_i'_X(0)=r_i'_X(X_0,Y_0,Z_0) \ r_i'_Y(0)=r_i'_Y(X_0,Y_0,Z_0) \ r_i'_Z(0)=r_i'_Z(X_0,Y_0,Z_0)$$

значение частных производных $r_i$ в точке нулевого приближения
Новый шаг для нового приближения (формулы 6) $X_1=X_0+\Delta X \ Y_1=Y_0+\Delta Y \ Z_1=Z_0+\Delta Z \ R_1=R_0+\Delta R$
Введем новые обозначения (формулы 7)

$$u=\Delta X \ v=\Delta Y \ w=\Delta Z \ t=\Delta R \ a_i=r_i'_X(X_0,Y_0,Z_0) \ b_i=r_i'_Y(X_0,Y_0,Z_0) \ c_i=r_i'_Z(X_0,Y_0,Z_0) \ d_i=r_i(X_0,Y_0,Z_0)$$

В этих обозначениях функция 4 примет вид (формула 8)

$$F(u,v,z,t)=\sum[r_i(X_0,Y_0,Z_0)+r_i'_X(X_0,Y_0,Z_0)u+r_i'_Y(X_0,Y_0,Z_0)v+r_i'_Z(X_0,Y_0,Z_0)w-(R_0+t)]^2$$

Минимизируем F, берем частные производные по u, v, w, t, приравниваем их к нулю, решаем систему четырех линейных уравнений с четырьмя неизвестными и т.о. находим новый шаг. Функция F сходиться за три такта. На сколько я понимаю это не градиентный метод наискорейшего спуска. Но вопрос мой сейчас в другом. Чтобы найти нулевое приближение мне понадобилось решить одну задачу. О ней в следующем посту, там же будет и вопрос.

-- Пн июл 12, 2010 18:26:35 --

На плоскости даны n точек $(x_i, y_i)$. Требуется найти прямую минимизирующую сумму квадратов расстояний от точек до этой прямой. Эта задача решается простым школьным методом так

Но мне на этом примере хотелось бы понять метод главных компонент. Вот теория http://ru.wikipedia.org/wiki/%D0%9C%D0% ... 0%BD%D1%82 Неясно как вот эти формулы записать для моего случая.

$$dist^2(x_i, L_k)=||x_i-a_0-\sum\limits_{i=1}^ka_j(a_j,x_i-a_0)||^2\\$$

$$dist^2(x_i,L_k)=\sum\limits_{i=1}^n(x_i_l-a_0_l-\sum\limits_{j=1}^ka_j_l\sum\limits_{q=1}^na_j_q(x_j_q-a_0_q))^2$$

Например даны три точки (1, 2), (3, 1), (4, 3). Как записать эти формулы?

$k = 1 \ L_1 = a_0+\beta_1a_1 \ n=3$

[Tur]

Цитата из той же википедии

2) находим первую главную компоненту как решение задачи;
$$a_1=argmin\Bigg(\sum\limits_{i=1}^m||x_i-a_1(a_1,x_i)||^2\Bigg)$

Что это значит? Как это выглядит при подстановке моих трех точек?
С $a_0$ тоже непонятно. Речь идет о векторах, но что такое нулевой вектор?

[Алексей К.]

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

На плоскости даны n точек $(x_i, y_i)$. Требуется найти прямую минимизирующую сумму квадратов расстояний от точек до этой прямой. Эта задача решается простым школьным методом так...
Но мне на этом примере хотелось бы понять метод главных компонент.

Что такое метод главных компонент я, признаться, не знаю. А может и знаю, просто не знаю, что это он.
Известное мне решение этой задачи (3D) я постил на форуме. Название "главные компоненты" ему очень идёт. Найду ссылочку, и сюда же вставлю.

-- 13 июл 2010, 16:00 --

Здесь же, где и про окружность.

-- 13 июл 2010, 16:15 --

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

[AKM]

А почему $y=kx$, а не $y=kx+b$? Как-то совсем неприкольно...

[Tur]

AKM, можно и так, но дело ведь не в этом. Нет ничего проще центрирования данных. Мне надо на простом примере понять метод. Вот есть три точки (вектора) X1(1, 2) X2(3, 1) X3(4, 2.5)

Надо найти прямую (ее проект изображен красным) минимизирующую сумму квадратов расстояний от нее до точек. Искать ее будем в виде $L = a_0 + t*a_1$. Находим нулевой вектор $a_0(2.67,1.83)$ показан белым. Вектор $a_1(u,v)$ показан голубым. Всего то надо найти u и v, причем $u^2+v^2 = 1$
Чтобы понять формулу $$a_1=argmin\Bigg(\sum\limits_{i=1}^m||x_i-a_1(a_1,x_i)||^2\Bigg)$ пришлось решить школьную задачку. Даны две точки А и B или два вектора $a(a1,a2) , x(x1,x2)$ Надо найти расстояние от В до ОА, т.е. длину ВС или модуль вектора b(b1,b2).

Нормализуем вектор $a(1,4)$, получим вектор $a_1(1/\sqrt{17}, 4/\sqrt{17}).$ Скалярное произведение $(a_1,x)=|a_1|*|x|*cos(\angle BOC) = 4/\sqrt{17}*1+1/\sqrt{17}*2 =6/\sqrt{17} $ Т.к. $|c| = |x|*cos(\angle BOC)$ умножив результат скалярного произведения на нормализованный $a_1$ получим вектор $\overrightarrow{OC}=c(c1,c2) = a_1(a_1,x)=(4,1)*6/\sqrt{17}=c(24/17,6/17)$.
Т.к. $\overrightarrow{BC}=\overrightarrow{OC}-\overrightarrow{OB}$ , то $\overrightarrow{BC}=\overrightarrow{b}=x-a_1(a_1,x)=(1,2)-(24/17,6/17)=(7/17,28/17)$ $|\overrightarrow{BC}|=\sqrt((7/17)^2+(28)^2)=1.6977$
Все действия в точности по формуле $||x_i-a_1(a_1,x_i)||^2$

Итак, центрируем данные, т.е. из Х1,Х2,Х3 вычитаем $a_0(2.67,1.83)$ получаем три новых вектора $X1(-1.67,0.17), X2(0.33,-0.83), X3(1.33,0.67)$ Требуется найти координаты $a_1(u,v)$

Обозначим для общности $X1(x1,y1),X2(x2,y2),X3(x3,y3)$ Учитывая что $||x_i-a_k(a_k,x_i)||^2=||x_i||^2-(a_k,x_i)^2$ раскроем сумму
$$a_1=argmin\Bigg(\sum\limits_{i=1}^m||x_i-a_1(a_1,x_i)||^2\Bigg)$

$$F_a_1=||X1||^2-(ux_1+vy_1)^2+||X2||^2-(ux_2+vy_2)^2+||X3||^2-(ux_3+vy_3)^2]$$

упростив получим функцию от u и v для минимизации

$$F_a_1=||X1||^2+||X2||^2+||X3||^2-[u^2(x_1^2+x_2^2+x_3^2)+2uv(x_1y_1+x_2y_2+x_3y_3)+v^2(y_1^2+y_2^2+y_3^2)]$$

или $F_a_1=Au^2+Bv^2+2Cuv+D$
Вот цитата из той же статьи в википедии

Цитата:

Векторы главных компонент для задач о наилучшей аппроксимации и о поиске ортогональных проекций с наибольшим рассеянием — это ортонормированный набор собственных векторов эмпирической ковариационной матрицы C

Будьте добры вопрос: какова связь между собственными векторами ковариационной матрицы и минимизацией $F_a_1$ ?

Алексей К. по ссылке про окружность похоже не то.

[ewert]

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

Будьте добры вопрос: какова связь между собственными векторами ковариационной матрицы и минимизацией $F_a_1$ ?

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

[Алексей К.]

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

AKM, можно и так, но дело ведь не в этом. Нет ничего проще центрирования данных.

Tur,
Ваш ответ AKMу я проинтерпретировал так: Вы решали задачу для уже центрированных данных, т.е. перешли предварительно в систему центра тяжести точек $(x_i,y_i)$. И Вам понятно, что при этом добавленное AKMом $b$

AKM в сообщении #338978 писал(а):

А почему $y=kx$, а не $y=kx+b$?

будет равно нулю. До сих пор не особо было ясно, что Вы этим пользуетесь и это осознаёте.

Также замечу, что это Ваше "нет ничего проще центрирования" на самом деле штука или просто важная, или архиважная. Важна она хотя бы с точки зрения точности вычислений: так в суммах типа $\sum x_i^2$ нужны и работают в итоге именно некоторые малые отконения от "идеальных" иксов, которые просто затеряются, если все иксы кучкуются где-то там возле миллиона, с разбросом плюс-минус десять, с отклонениями плюс-минус один. Но эти плюс-минус однушки не затеряются, если этот $\bar x=\dfrac1N\sum x_i=$миллион посадить в ноль: плюс-минус единички (ошибки) будут уже фигурировать в этих суммах на фоне десяток.
Архи- возникает, видимо, от того, что в этих всех ковариционных марицах наверняка тоже фигурируют не $\sum x_i y_i$, а $\sum (x_i-\bar x)(y_i- \bar y)$. Я плохо всё это помню и не буду говорить в терминах ковариаций.

Вернёмся к задаче о поиске прямой $ux+vy+l=0$ ($u^2+v^2=1$), ближайшей к заданному набору точек $(x_i,y_i)$. Предположим для начала (что не суть важно), что точки вполне равномерно-симметрично расположены на некотором эллипсе, малость или сильно вытянутом. Вам, видимо, уже понятно, что после центрирования ответ будет $l=0$, и можно возиться с прямой $ux+vy=0$. Выше (на рисунке из Вашего сообщения post338784.html#p338784 ) Вы подменили её прямой $y=kx$. Найдены два экстремальных значения $k_{1,2}$. Так вот, эти два значения есть максимум и минимум. Наилучшая прямая и наихудшая. И $k_1k_2=-1$: прямые взаимно перпендикулярны. Это две полуоси нашего эллипса.

(Теперь мельком подумайте о 3D-варианте, об эллипсоиде, о трёх экстремумах..., и возвращайтесь назад, в 2D).

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

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

Алексей К. по ссылке про окружность похоже не то.

То. Тема про окружность, но ссылка обсуждает подзадачу с плоскостью-прямой. Просто вместо той матрицы я буду пользовать матрицу попроще, 2х2, без z.
А Вы, кстати, умеете искать экстремум функции $F(u,v)$ при дополнительном условии типа $u^2+v^2=1$ (метод неопр. множителей Лагрангжа)? А знаете, что такое собственный вектор матрицы? Даже совсем маленькой матрицы, 2х2, вполне сойдёт...

[Алексей К.]

Продолжаю...

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

Находим нулевой вектор $a_0(2.67,1.83)$ показан белым.

Ну да, этот вектор у Вас ровно в центр тяжести точек и лезет. С этим, похоже, всё ясно. Точки далее считаем центрироваными: $\sum x_i = \sub y_i = 0$.

Итак, ищем минимум (экстремум) функции $F(u,v)=\sum (u x_i + v y_i)^2$ при условии $u^2+v^2-1=0$. Мы, конечно, можем избавиться от этого условия, положив $u=-\sin\xi$, $v=\cos\xi$, сделав функцию одной переменной $\Phi(\xi)$, и решая итерациями-линеаризациями, как уже освоено для исходной задачи.
Но мы не будем так делать.

По методу Лагранжа, для отыскания условного экстремума мы должны экстремизировать функцию

$$L(u,v,\lambda)=F(u,v)+\lambda(u^2+v^2-1).$$

Уравнения $L'_{u,v,\lambda}=0$ принимают вид

$$ \begin{align*}u\sum x_i^2 + v\sum x_i y_i = \lambda u,\\ u\sum x_i y_i + v\sum y_i^2 = \lambda v,\\ u^2+v^2=1.\end{align*}$$

Находим из первых двух $u(\lambda),\:v(\lambda)$, и всё решается. Будет два решения.

Теперь нам надо увидеть в этом задачу на собственные вектора.

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

$$M= \left| \begin{array}{lll} \sum x_i^2 & \sum x_i y_i & \sum x_i z_i\\ \sum x_i y_i & \sum y_i^2 & \sum y_i z_i\\ \sum x_i z_i & \sum y_i z_i & \sum z_i^2 \end{array} \right| $$

(в системе центра тяжести). Ищутся они легко, особенно если программкой, а не формулой...
Собственные значения неотрицательны; соб. вектор, соответствующий максимальному с.з., даёт направление наилучшей МНК-прямой. Чтобы найти его, достаточно взять какое-то первое приближение $v_0$ и слегка поитерировать: $v_{i+1}=M v_i$...

У нас (суммы малость переобозначил)

$$ \left|\begin{array}{ll} S_{xx} & S_{xy}\\ S_{xy}& S_{yy}\end{array} \right| \cdot \left|\begin{array}{l}u\\v \end{array} \right| = \lambda \left|\begin{array}{l}u\\v \end{array} \right|, \text{~~~~что приводит к~уравнению~~~~} \mathrm{det}\left|\begin{array}{cc} S_{xx}-\lambda & S_{xy}\\ S_{xy}& S_{yy}-\lambda\end{array} \right| = 0. $$

Откликайтесь, ежели надо дальше колотить... а то я в этом не уверен...

[Tur]

Алексей К., спасибо за разъяснения, сейчас гораздо яснее. Не совсем ясен смысл лямбда в $L(u,v,\lambda)=F(u,v)+\lambda(u^2+v^2-1)$

Цитата:

Находим из первых двух $u(\lambda),\:v(\lambda)$, и всё решается.

Это тоже неясно, т.к. определитель системы ноль, т.е. первые два ур-я линейно зависимы. Я сначала нашел лямбда как максимальное собственное число, затем первое ур-е поделил на $u$ и выразил из него $v/u$, найдя т.о. направление вектора.
Следующая задача та же, только в пространстве. Точки центрированы, надо найти только вектор $a_1(u,v,w)$. В следующей формуле он ищется изначально нормализованным. $a_1=argmin\Bigg(\sum\limits_{i=1}^m||x_i-a_1(a_1,x_i)||^2\Bigg)$
Учитывая $||x_i-a_k(a_k,x_i)||^2=||x_i||^2-(a_k,x_i)^2$ запишем $F(u,v,w)=\sum{(ux_i+vy_i+wz_i)^2}$
Формально записываю как в предыдущем случае функцию для минимизации

$$L(u,v,w,\lambda)=F(u,v,w)+\lambda(u^2+v^2+w^2-1)$$

Уравнения $L'_{u,v,w\lambda}=0$ принимают вид

$$ \begin{align*}u(S_x_x-\lambda) + vS_x_y + wS_x_z = 0\\ uS_x_y + v(S_y_y-\lambda) + wS_y_z = 0\\ uS_x_z + vS_y_z + w(S_z_z-\lambda) = 0\\ u^2+v^2+w^2=1.\end{align*}$$

отсюда

$$\mathrm{det}\left|\begin{array}{ccc} S_{xx}-\lambda & S_{xy} & S_{xz}\\ S_{xy} & S_{yy}-\lambda & S_{yz}\\ S_x_z & S_y_z & S_z_z-\lambda\end{array} \right| = 0. $$

Решаем кубическое ур-е и находим лямбда, подставляем максимальную в систему и полагаем, например, $w = 1$. Получаем систему (например первые два) двух линейно независимых уравнений. Находим $u, v$ и нормализуем их с $w$. Все верно?

ewert,

Цитата:

Сумма квадратов расстояний -- это сумма квадратов проекций радиус-векторов всех точек на направление нормали. Которая довольно легко сводится к квадратичной форме той самой "ковариационной матрицы" на векторе нормали, нормированном на единицу.

Я и спрашивал как при минимизации $F(u,v)=\sum{(ux_i+vy_i)^2}$ придти к ковариационной матрице.

Цитата:

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

У меня почему то минимум получился при максимальном лямбда.

[ewert]

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

Я и спрашивал как при минимизации $F(u,v)=\sum{(ux_i+vy_i)^2}$ придти к ковариационной матрице.

Хорошо, детально:

$F(u,v)=\sum\limits_{i=1}^n(u^2x_i^2+2uvx_iy_i+v^2y_i^2)=$

$=u^2\sum\limits_{i=1}^nx_i^2+uv\sum\limits_{i=1}^nx_iy_i+vu\sum\limits_{i=1}^ny_ix_i+v^2\sum\limits_{i=1}^ny_i^2\equiv(C\vec n,\vec n)\;,$

где $\vec n=(u,v)$ -- это единичный вектор нормали к прямой и $C=\begin{pmatrix}\sum x_i^2 & \sum x_iy_i \\ \sum y_ix_i & \sum y_i^2 \end{pmatrix}$ -- ковариационная матрица (учитывая, конечно, что всё уже центрировано).

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

[Алексей К.]

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

Не совсем ясен смысл лямбда в $L(u,v,\lambda)=F(u,v)+\lambda(u^2+v^2-1)$

Надо посмотреть тему условный экстремум, метод множителей Лагранжа.

Tur писал(а):

Цитата:

Находим из первых двух $u(\lambda),\:v(\lambda)$, и всё решается.

Это тоже неясно, т.к. определитель системы ноль, т.е. первые два ур-я линейно зависимы.

А здесь я не доглядел и лопухнулся.

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