Розв’язання задачі лінійної регресії за допомогою швидкого перетворення Хафа
Друзі, розглянемо нині завдання лінійної регресії у присутності викидового (некорельованого з сигналом) шуму. Це завдання часто виникає при обробці зображень (напр., при колірній сегментації [1]), у тому числі акустичних [2]. У випадках, коли координати випадкових величин можна грубо дискретизувати, а розмірність завдання низька (2-3), крім стандартних методів робастної регресії, можна скористатися швидким перетворенням Хафа (БПХ) [3]. Спробуємо порівняти цей останній метод за точністю та стійкістю з «класичними».
Використання БПХ для лінійної регресії
Завдання лінійної регресії на площині полягає у відновленні лінійної залежності між двома змінними, заданими у вигляді множини пар (x, y). Задавшись деяким рівнем дискретизації координат, можна відобразити цю множину на однобітному або цілісному зображенні (у першому випадку ми відзначаємо лише факт наявності у вихідних даних точки з приблизно такими координатами, у другому — ще й їхнє число). Фактично йдеться про двовимірну гістограму вихідних даних. Таким чином, неформально завдання може бути зведена до пошуку на зображенні прямої, яка найкраще описує зображений розподіл точок. У обробці зображень у подібних випадках використовується перетворення Хафа.
Перетворення Хафа є дискретним аналогом перетворення Радона і ставить у відповідність кожної прямої на зображенні суму яскравостей пікселів уздовж неї (тобто одночасно обчислює всілякі суми уздовж дискретних прямих). Можна ввести розумну дискретизацію прямих по зсувах і нахилам так, щоб паралельні прямі дискретні щільно упаковували площину, а зображення, що виходять з однієї точки на одному краї, прямі розходилися по нахилу на протилежному краї націла кількість пікселів. Тоді таких дискретних прямих квадраті n 2 буде приблизно 4 * n 2 . Для цієї дискретизації існує алгоритм швидкого обчислення перетворення Хафа з асимптотикою O(n 2 *log n). Цей алгоритм є близьким аналогом алгоритму швидкого перетворення Фур'є, що добре паралелізується і не вимагає жодних операцій, крім складання. У роботі [3] можна прочитати трохи більше, крім того, там пояснюється, чому перетворення Хафа від згладженого гауссовским фільтром зображення взагалі можна застосовувати в задачі лінійної регресії. Тут ми продемонструємо стійкість цього методу.
Тестові дані
Поставимо тепер питання, який розподіл ми очікуємо у вхідних даних. В задачі акустичної локації можна вважати, що у нас є рівномірно розподілені точки вздовж прямолінійного відрізка траєкторії руху судна, спотворені нормальним адитивним координатним шумом. У разі аналізу кольорових кластерів розподіл уздовж відрізка невідомий і часто дуже химерний. Для простоти розглядатимемо виключно двомірний гауссовий розподіл із високим ексцентриситетом. При досить велику дисперсію вздовж прямої і малої дисперсії в поперечному напрямку властивості побудованого розподілу не сильно відрізняються від «зашумленого відрізка».
Алгоритм генерації тестових розподілів може бути описаний такими трьома кроками:
- випадково вибираються «центр» і нахил (моделюється відрізок прямої);
- генерується задану кількість точок із зазначеним розподілом (і вибраним центром);
- додається однорідний викидний шум.
Стандартне відхилення нормального розподілу вздовж прямої становило одну четверту довжини відрізка, що моделюється (при такому виборіприблизно 95% крапок потрапляють всередину), а довжина відрізка бралася рівною половині лінійного розміру квадратного зображення (що становить 512 пікселів). Відношення дисперсій дорівнювало 10.0, а кількість точок розподілу і частка викидного шуму варіювалися. Приклад такого роду розподілу зображено малюнку 1.
Щоб порівнювати різні алгоритми, треба вміти оцінювати, наскільки добре було знайдено відрізок у кожному разі. Тонкість у тому, що аналізовані алгоритми шукають не відрізок, а пряму. Тому братимемо «ідеальний» відрізок (A, B) і знаходитимемо максимум з двох відстаней від його кінців до експериментально знайденої прямої. Якщо позначити через Q(L, AB) помилку в параметрах експериментальної прямої L, це визначення можна записати так: Q(L, AB) = max(r(L, A), r(L, B)).

Мал. 1. Зразок розподілу з адитивним та викидним шумами. Позначками показані кінці і центр відрізка, що моделюється.
Залежність помилки методу БПХ від розміру вікна згладжування
Як було зазначено вище, регресія з допомогою БПХ передбачає згладжування. Очевидно, що розмір вікна згладжування має істотно впливати на якість роботи алгоритму. Досліджуємо залежність помилки знаходження прямої від розміру вікна. Для кожного розміру деякого набору була згенерована серія зі ста зображень, потім для кожного обчислювалася пряма і фіксувалася помилка. Отриманий результат представлений на малюнку 2, де вісь абсцис задає відношення параметра sH, що згладжує, до стандартного відхилення нормальної («поперечної»), що становить координатного шуму sL. Крапками на малюнку позначаються середні помилки, а вгору та вниз від них відкладаються середньоквадратичні відхилення. Можна помітити, що результат є найкращим уу разі, коли sH близько до sL, тобто. коли розмір вікна не дуже відрізняється від «товщини» відрізка.

Мал. 2. Залежність помилки методу Хафа від розміру вікна згладжування.
Порівняння з МНК та його роботними модифікаціями
Найбільш очевидним, хоч і неправильним (у нашому випадку) методом лінійної регресії є, звичайно ж, метод найменших квадратів (МНК), який для ортогональної регресії збігається з методом головних компонентів (англ. PCA). Свого часу для випадку викидного шуму його допилили до «МНК з ітеративним перерахуванням терезів», причому до всіх підручників увійшли варіанти з ваговою функцією Х'юбера та біквадратною ваговою функцією [4].
Порівняємо метод БПХ з простим МНК, а також з МНК з перерахуванням ваги з обома ваговими функціями. Параметри функцій візьмемо, як рекомендують теоретики: kH = 1345 для Хьюбера і kB = 4685 для біквадратної функції.
На малюнку 3 видно, як середня помилка кожного з чотирьох методів залежить від густини викидного шуму µ - відношення кількості викидних шумових точок до загального числа пікселів на зображенні (щільність викидного шуму відкладена по осі абсцис, а середня помилка методу представлена на осі ординат). Кожен замір проводився на ста різних згенерованих зображень. З малюнка очевидно, що метод, заснований на застосуванні БПХ, виявився найстійкішим до додавання шуму викиду, хоча при малій кількості викидних точок він поступається іншим методам в точності. До того ж, цей метод виявляється стабільнішим, тобто рідше помиляється катастрофічно (вгору і вниз від кожної точки на малюнку відкладено чверть середньоквадратичного відхилення помилки).

Мал. 3. Залежність середньої помилки від щільності викидного шуму у БПХ та МНК-методів.
Порівняння з медіаннимиметодами
Дивно, але про по-справжньому хороші методи робастної регресії рідко розповідають у ВНЗ. Але ми з вами, дорогий читачу, знаємо про чарівний метод Тейла-Сена [5], і тому обмежитися порівнянням з роботним МНК було б безчесно. У класичному алгоритмі Тейла-Сена нахил прямої обчислюється як медіана серед нахилів прямих, проведених через усі пари точок вибірки. Крім цього використовується модифікований варіант, в якому спочатку для кожної точки розподілу вибиралася медіана нахилів всіх прямих, що проходять через цю та якусь іншу точку розподілу, після чого з отриманих значень вибиралася медіана («медіана медіан»). У обох випадках визначення параметра зсуву прямий береться медіана серед значень i — m xi>, де m – знайдений описаними способами кутовий коефіцієнт прямий.
На малюнку 4 показано порівняння тепер із цими методами. У рамках кожної серії, як і раніше, проводилося сто експериментів, позначення теж не змінилися. Видно, що в розглянутому інтервалі густини викидного шуму середня помилка методу БПХ є найменшою і не зростає, тоді як методи Тейла-Сіна мають пропорційне зростання помилки. Також видно, що метод БПХ, як і раніше, має найменший розкид помилки.

Мал. 4. Залежність середньої помилки від щільності викидного шуму у БПХ та медіанних методів.
Порівняння методів за наявності структурованого шуму
Окремий інтерес представляє випадок, коли викид шум має власну просторову структуру. Для прикладу додамо до вже розглянутих видів шумів некорельовану з сигналом шумову компоненту з двовимірним круговим нормальним розподілом. Приклад комбінації всіх трьох видів шумів наведено малюнку 5. Прямий показанийідеальна відповідь.

Мал. 5. Приклад відрізка, зашумленого адитивним, а також рівномірним та структурованим викидами.
Подивимося тепер залежність помилки від щільності «круглого» структурованого шуму при фіксованих параметрах викидного і адитивного шумів. Як і раніше, для кожного значення щільності згенеровано сто зображень. Результати наведено малюнку 6. Помітно, що методи розділилися втричі групи, причому дивним чином класичні робастные методи мало від МНК. Видно, що «зрив» у всіх МНК-методів відбувається практично відразу, медіанні методів стабільно працюють при набагато більших забрудненнях даних, але теж рано чи пізно перестають працювати, тоді як у методу БПХ різкого зриву немає зовсім.

Мал. 6. Залежність середньої помилки від густини структурованого шуму.
Обговорення
Вирішення задачі лінійної регресії за допомогою комбінації гаусового згладжування та перетворення Хафа є стійким до одночасно присутніх адитивного координатного та викидного шумів. Порівняння з іншими поширеними робастними методами показало, що запропонований алгоритм є більш стійким, коли щільність шуму викидів досить висока. Слід зазначити при цьому, що метод застосовується лише в маломірному випадку (алгоритм БПХ може бути узагальнений на будь-яку розмірність, але обсяг пам'яті, необхідний для кодування завдання у вигляді зображення розмірності вище трьох - лякає). Крім того, навіть у двовимірному випадку використання БПХ буде невигідним, якщо точок мало, а потрібна відносна координатна точність висока. Це з тим, що з інших методів обчислювальна складність залежить від кількості точок, а й у БПХ — від числа дискретпростору. Втім, за малої кількості точок можливе обчислення перетворення Хафа «грубою силою». Здавалося б, так не вийде, оскільки заздалегідь ми хочемо згладити зображення, проте ще Ізраїлю Мойсейовичу Гельфанду було відомо, що гаусове згладжування у певному сенсі перестановкове з перетворенням Радона. Використовуючи цей трюк, асимптотику методу в наближенні малого числа точок можна поліпшити до O(n 2 ), проте n, як і раніше, означає лінійний розмір зображення, так що майже напевно це все одно занадто погано. Зате за великої кількості точок метод БПХ криє решту методів як бик — вівцю. Такі справи.
Література
Хардкорна конфа за С++. Ми запрошуємо лише профі.