Гессиан и смешанная частная производная

[AlexMAS]

Здравствуйте.

Хотелось бы попросить Вашей консультации по следующему вопросу. Допустим, дана функция двух переменных D=f(x,y). Для этой функции в некоторой точке строится матрица Гессе H=[ Dxx Dxy ; Dxy Dyy ]. Компоненты Dxx и Dyy, если я правильно понимаю, определяют нормальные кривизны в направлении осей Ox и Oy соответственно. Какую кривизну определяют смешанные частные производные Dxy?

На форуме я встречал ветку, в которой пытались объяснить геометрический смысл смешанной производной функции двух переменных. Вот тут поправьте меня, если я не прав. Можно ли считать, что компоненты Dxy в матрице Гессе определяют нормальную кривизну в заданной точке в направлении dx·dy (то есть где-то по "диагонали" между Ox и Oy)?

Второй вопрос касается случаев, когда Dxy=0. Когда это возможно и что это дает? Иначе говоря, какой вывод я могу сделать, если я буду знать (или условлюсь), что компоненты Dxy матрицы H равны нулю?

Наконец, последнее, с чего все и началось. Я разбираюсь с алгоритмом поиска особых точек на изображении (SIFT) и в описании встретил соотношение вида:

Tr(H)^2 / Det(H)

где Tr(H) = Dxx + Dyy – след матрицы H, а Det(H) = Dxx * Dyy – (Dxy)^2 – определитель матрицы H. Мне не ясен смысл этого соотношения. Возможно, есть какая-то формула (опять же, связанная с кривизной поверхности), которую я упустил из виду?

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

Спасибо!

P.s. Ссылка на статью, которая вызвала вопросы: http://habrahabr.ru/blogs/algorithm/106302

[svv]

AlexMAS писал(а):

Компоненты Dxx и Dyy, если я правильно понимаю, определяют нормальные кривизны в направлении осей Ox и Oy соответственно.

Для кривизны важны еще и первые производные. Рассмотрим полуцилиндр $z=\sqrt{1-x^2}$, ось которого параллельна $Oy$. Нормальная кривизна сечения плоскостью $Oxz$ всюду одна и та же, а производная $\frac {\partial^2 z }{\partial x^2}$ разная.
Обратный пример: параболический цилиндр $z=x^2$. Производная $\frac {\partial^2 z }{\partial x^2}$ всюду равна $2$, но кривизна зависит от $x$ (легко заметить, что парабола $z(x)=x^2$ вблизи нуля более искривлена, чем вдали от нуля).
Поэтому лучше говорить не о нормальной кривизне, а о второй производной. Это не только $\frac {\partial^2 z }{\partial x^2}$ и $\frac {\partial^2 z }{\partial y^2}$, но и производная по произвольному направлению -- мы ее сейчас вычислим.

AlexMAS писал(а):

Можно ли считать, что компоненты Dxy в матрице Гессе определяют нормальную кривизну в заданной точке в направлении dx·dy (то есть где-то по "диагонали" между Ox и Oy)?

Пусть вектор $\mathbf{a}$, задающий направление дифференцирования в точке, получается поворотом вектора $\mathbf{e}_x$ против часовой стрелки на угол $\varphi$. Тогда будем характеризовать направление этим углом $\varphi$.
Найдем вторую производную по направлению $\varphi$. Введем систему координат $(\xi, \eta)$, повернутую против часовой стрелки на угол $\varphi$ относительно координат $(x, y)$:
$x = \xi \cos \varphi - \eta \sin \varphi$
$y = \xi \sin \varphi + \eta \cos \varphi$
Координатные линии $\xi$ имеют как раз направление $\varphi$, поэтому находим
$\frac {\partial} {\partial\xi} = \frac {\partial x}{\partial\xi} \frac {\partial} {\partial x} + \frac {\partial y} {\partial\xi} \frac {\partial} {\partial y} = \cos \varphi \frac {\partial} {\partial x} + \sin \varphi \frac {\partial} {\partial y}$,
$\frac{\partial^2 z}{\partial\xi^2} = (\cos \varphi \frac {\partial}{\partial x} + \sin \varphi \frac {\partial}{\partial y})^2 z = \cos^2 \varphi \frac {\partial^2 z }{\partial x^2} + 2 \cos \varphi \sin \varphi \frac {\partial^2 z }{\partial x\partial y} + \sin^2 \varphi \frac {\partial^2 z }{\partial y^2}$.
Это и есть вторая производная функции по направлению $\varphi$. Видно, что "задействованы" все вторые производные, и ни при каком $\varphi$ не остается только $\frac {\partial^2 z }{\partial x\partial y}$. Поэтому $\frac {\partial^2 z }{\partial x\partial y}$ не имеет смысла второй производной по какому-то "косому" направлению.
Пример: параболический цилиндр $z=(y-x)^2$. Значение $z$ не меняется в "диагональном" направлении $\varphi=\pi/4$, но $\frac {\partial^2 z }{\partial x\partial y}$ прекрасненько отлична от нуля.

$\frac {\partial^2 z }{\partial x\partial y}$ имеет такой смысл: это скорость изменения $\frac {\partial z }{\partial x}$ в направлении $Oy$ = скорость изменения $\frac {\partial z }{\partial y}$ в направлении $Ox$. В "трехмерном" представлении это выглядит так. Возьмем на поверхности $ z=f(x, y)$ точку $(x_0, y_0, z_0)$, проведем через нее плоскость, параллельную $Oxz$, получится сечение -- кривая $z(x)=f(x, y_0)$. Угловой коэффициент касательной к ней в точке $x_0$ -- это $\frac {\partial z}{\partial x}$. Теперь будем двигать точку $y_0$, повторяя построение. Скорость изменения углового коэффициента -- и есть смешанная производная $\frac {\partial^2 z }{\partial x\partial y}$. Если угловой коэффициент сечения $z(x)=f(x, y_0)$ не меняется при изменении $y_0$ -- смешанная производная равна нулю.

Моё мнение: "трехмерные" дифференциально-геометрические аналогии были бы полезны, если бы координата $z=f(x, y)$ была в каком-то смысле равноценна координатам $x$, $y$. Тогда, например, имел бы смысл поворот координат или поверхности вокруг произвольной оси, и дифгеометрические инварианты работали бы в полную силу. А так -- у нас все-таки область в двумерном пространстве $(x, y)$, в которой задана функция $z=f(x, y)$. Тривиальное преобразование $z \mapsto kz$ изменит поверхность с постоянной кривизной (сферу) в поверхность с переменной кривизной (эллипсоид). Но оно не превратит функцию с постоянной второй производной в функцию с переменной второй производной.

[AlexMAS]

Большое спасибо за подробный ответ!

Однако остался вопрос, связанный с соотношением $Tr(H)^2/Det(H)$. Далее я приведу свои суждения по этому поводу. Пожалуйста, поправьте меня, если я где-то неправ...

Рассмотрим все ту же матрицу Гессе для заданной функции $z=f(x,y)$:

$$H=\left(\begin{array}{cccc} H_{xx} & H_{xy}\\ H_{xy} & H_{yy}\end{array}\right)$$

Задача состоит в том, чтобы определить, находится ли данная точка $(x_{0},y_{0})$ вблизи локального экстремума. Для решения задачи найдем главные нормальные кривизны функции $z=f(x,y)$ в точке $(x_{0},y_{0})$, произведя аппроксимацию заданной функции соприкасающимся параболоидом вида $z=ax^2+by^2$, где $a=H_{xx}$, $b=H_{yy}$. Вычислив $a$ и $b$ (определив компоненты $H_{xx}$ и $H_{yy}$ в заданной точке), получим максимальную и минимальную кривизны, по которым можем судить о том, какой вид имеет соприкасающаяся поверхность: если $a$ и $b$ - одного знака, то эллиптический параболоид (экстремум есть), если разного - гиперболический параболоид (экстремума нет).

Положим, что $a=rb$, тогда:
$\frac {Tr(H)^2}{Det(H)}=\frac {(a+b)^2}{(ab)}=\frac {(r+1)^2}{r}$

Так как при аппроксимации было принято, что $H_{xy}=0$, то:
$\frac {Tr(H)^2}{Det(H)} \le \frac {(r+1)^2}{r}$

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

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

[svv]

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

Аппроксимация $ax^2+by^2$ (если понимать формулу буквально) не будет хороша. Формула жестко фиксирует положение этой фигуры (сдвиг параллельно осям не допускается) и её ориентацию (поворот в плоскости $Oxy$ не допускается). Если Вы хотите иметь аппроксимацию второго порядка (аппроксимирующая функция имеет в данной точке "правильные" значение, первые и вторые производные), то двух параметров $a$, $b$ будет недостаточно.

Пусть нас интересует точка $(x_0, y_0)$.

1. Переместим туда начало координат, но, чтобы не плодить новых обозначений, назовем новые координаты опять $x$ и $y$ (ну, или будем считать, что в рассматриваемой точке сразу было начало координат). Итак, рассматриваем окрестность точки $(0, 0)$.

2. Разложим функцию $z=f(x,y)$ в окрестности $(0, 0)$ в ряд по степеням $x$ и $y$:
$z=\sum\limits_{m=0}^{\infty}\sum\limits_{n=0}^{\infty} c_{mn} x^m y^n$.

3. Замечаем, что вклад в смешанную производную в нуле дает только слагаемое $c_{11}xy$. Действительно,
$\frac {\partial^2 z}{\partial x \partial y} = \sum\limits_{m=1}^{\infty}\sum\limits_{n=1}^{\infty} mn \,c_{mn}\,x^{m-1} y^{n-1}$.
И если хоть один из индексов $m$, $n$ больше $1$, в нуль обращается $x^{m-1} y^{n-1}$. А суммирование здесь с $1$, так как смешанная производная от $c_{00}+c_{10}x+c_{01}y$ равна нулю. Остается только $m=n=1$.

4. Но слагаемое $c_{11}xy$ всегда можно устранить поворотом системы координат вокруг начала. Если требуется, доказательство (оно совершенно "лобовое") за мной.
Производим нужный поворот, обозначая новые координаты опять через $x$ и $y$, чтобы не плодить обозначений. Путаницы не будет.

5. Замечаем, что в нашей точке вклад в производную $\frac {\partial^2 z}{\partial x^2}$ дает только слагаемое $c_{20}x^2$, а вклад в производную $\frac {\partial^2 z}{\partial y^2}$ дает только слагаемое $c_{02}y^2$.

6. Теперь, так как нас интересуют вторые производные, мы уже совершенно законно оставляем в сумме только слагаемые $c_{20}x^2=ax^2$ и $c_{02}y^2=by^2$.

Итак, $c_{11}xy$ и тем самым смешанная производная исчезли не потому, что мы просто выбрали такую аппроксимацию (это была бы очень плохая аппроксимация), а потому, что мы совершили нужный поворот и наши главные оси стали параллельны осям координат. Если допустить нужный поворот параболоида, а также нужные сдвиги его по $x$, $y$ и $z$, параболоид сможет прилегать к поверхности в окрестности точки гораздо плотнее.

AlexMAS писал(а):

что в итоге проверяет указанное неравенство я так и не понял

И я не знаю. Можно заметить, что $\operatorname{Tr} H$ -- это лапласиан (часто используемый при обработке изображений). Вы хорошо описали, на что он указывает (фраза про эллиптический и гиперболический параболоиды). Возведение его в квадрат устраняет зависимость от направления вверх/вниз.
И можно предположить, что деление на $\det H$, делая результат безразмерным, устраняет влияние преобразования $z \mapsto kz$ на результат -- что-то вроде нормировки.

[svv]

$\frac {(\operatorname{tr}H)^2}{\det H}=\operatorname{tr}H \operatorname{tr}H^{-1}$
Красивая инвариантная форма. Но что это даёт?

[AlexMAS]

Здравствуйте, спасибо за замечания!

По поводу соотношения $Tr(H)^2/Det(H)$, мне кажется, что вы правы. Я долго копался в первоисточниках алгоритма SIFT в поисках истины, пока не увидел ссылку на статью «A combined corner and edge detector» (Harris and Stephens, 1988), из которой, со слов автора, взята идея указанного соотношения.

В этой статье для детектирования углов и линий Harris предлагает определить меру, которая является функцией, зависящей только от $a$ и $b$, где $a$ и $b$ – главные нормальные кривизны заданной функции в исследуемой точке. Далее дается определение меры:

$R(a,b)=Det(H)-k*Tr(H)^2$, где $Det(H)=a*b$, $Tr(H)=a+b$.

Если $R>0$, то точка является угловой (эллиптическая точка), если $R<0$, то точка принадлежит линии (гиперболическая точка), если же значение $R$ близко к нулю, то в окрестности точки заданная функция практически не меняется (точка уплощения).

Есть подозрения, что автор алгоритма SIFT ввел подобную эмпирическую меру $Tr(H)^2/Det(H)$ и, введя соотношение $Tr(H)^2/Det(H)<(r+1)^2/r$, производит отсев ненужных точек. На практике, действительно, все работает. Другой вопрос – хотелось бы дать этому соотношению (неравенству) хоть какое-то математическое обоснование.

P.s. Вообще, матрица Гессе очень интересная и хитрая вещь :) К виду $\frac {(\operatorname{tr}H)^2}{\det H}=\operatorname{tr}H \operatorname{tr}H^{-1}$ я бы, наверняка, не пришел самостоятельно. Но что это может значить? Может быть какой-нибудь инвариант? Например, относительно вращения (?). Об этом автор также упоминает...

[svv]

Да, $\operatorname{tr} H \operatorname{tr} H^{-1}$ инвариантно относительно линейных преобразований координат с постоянными коэффициентами, сюда относятся и вращения.

Посмотрел статью Харриса и Стивенса. Я так понял, они имеют в виду следующее. Вообще-то критерий должен использовать собственные числа матрицы $\alpha$ и $\beta$. Но поскольку в двумерном случае полную информацию о собственных числах дают $\operatorname{tr}$ и $\operatorname{det}$, а вычисляются они проще, то полезно выразить критерий в терминах именно $\operatorname{tr}$ и $\operatorname{det}$. Что они и делают.

[zluka]

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

[AlexMAS]

Разложение функции в ряд Тейлора необходимо для аппроксимации функции $D(x,y,\sigma)$, чтобы, во-первых, определить какие точки, полученные на этапе первичной выборки, действительно являются экстремальными, а какие нет; во-вторых, эта же функция используется для интерполяции (уточнения) положения найденных экстремумов (ведь на этапе оценки мы получили только приблизительное положение точек).

Чтобы найти экстремум функции производную функции $D'(x,y,\sigma)$ приравнивают к нулю. При получении производной дифференциалы высших порядков игнорируют, считая что они вносят незначительный вклад. Из этого уравнения получают вектор смещения - смещение экстремума относительно оценки, полученной на первом этапе. И так до тех пор, пока смещение не будет удовлетворять заданной точности.

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

-- Пт сен 23, 2011 09:38:53 --

@zluka Для простоты можно рассмотреть случай производной функции одной переменной:

$f(x)=f(x_0)+f'(x_0)(x-x_0)+0.5f''(x_0)(x-x_0)^2+\varepsilon(x-x_0)$

где $\varepsilon(x-x_0)$ очень близка к нулю, $x_0$ - предполагаемое положение экстремума. Тогда:

$f'(x)=f'(x_0)+0.5f''(x_0)(2x-2x_0)=f'(x_0)+f''(x_0)(x-x_0)$

Приравняв $f'(x)$ к нулю, получим смещение до экстремума:

(1): $dx_e=(x-x_0)=-(f''(x_0))^{-1}f'(x_0)$

Теперь нужно найти значение функции в точке экстремума $f(x_e)=f(x_0+dx_e)$:

$f(x_e)=f(x_0)+f'(x_0)(x_0+dx_e-x_0)+0.5f''(x_0)(x_0+dx_e-x_0)^2=$
$=f(x_0)+f'(x_0)dx_e+0.5f''(x_0){dx_e}^2$

Теперь, если сюда подставить выражение (1) получим:

$f(x_e)=f(x_0)+f'(x_0)dx_0-0.5f'(x_0)dx_e$

или:

(2): $f(x_e)=f(x_0)+0.5f'(x_0)dx_e$

[zluka]

Спасибо. Я просто боялся что нужно дифференцировать все 3 члена и решать матричное квадратное уравнение.