Обложка: Задача на собеседовании: провести прямую через набор точек

Задача на собеседовании: провести прямую через набор точек

Дмитрий Хизбуллин
Дмитрий Хизбуллин

ведущий инженер в исследовательском центре Huawei

Будучи энтузиастом в области беспилотных автомобилей, я проходил собеседования в ряд компаний на позицию разработчика программного обеспечения. Чаще всего компании, занимающиеся беспилотными автомобилями, разрабатывают машину, которая может проехать по местности и хорошо ориентируется в 3D окружении. Следовательно, фирмы беспокоятся о том, насколько хорош кандидат в вычислительной геометрии, и в частности, как легко и быстро он/она может закодить решение. Один из вопросов на собеседовании меня особенно заинтересовал. Задача формулируется следующим образом:

Представьте, что ваш автономный автомобиль едет по дороге прямо, и его 2D позиции определяются при помощи ГНСС (Глобальной Навигационной Спутниковой Системы, например GPS) не очень точно, а с определенной погрешностью. Требуется найти наиболее вероятные положения автомобиля такие, что они лежат на одной прямой.

Эту задачу, как и практически любую другую, можно решать несколькими способами. В этой статье мы рассмотрим 3 решения на языке Python и сравним их: посмотрим на достоинства и недостатки. С математической точки зрения задание формулируется как задача линейной регрессии: найти аппроксимацию набора точек при помощи прямой. С программисткой стороны нам понадобится подходящая библиотека по манипулированию векторами и матрицами. Мы будем использовать NumPy для этой цели.

Подход 1: Метод наименьших квадратов

На первый взгляд, нам подойдет стандартное решение линейной регрессии: метод наименьших квадратов.

Нам дан массив 2D координат, представляющих собой пары точек (x, y). Мы можем считать x за аргумент, а y — за значение функции. Тогда, проконсультировавшись с Википедией, мы можем рассчитать коэффициенты прямой при помощи следующей формулы.

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

def least_squares(points: np.ndarray, axis: Optional[Any] = None) \
      -> np.ndarray:
   """
   Функция для аппроксимации массива точек прямой, основанная на
   методе наименьших квадратов.

   :param points: Входной массив точек формы [N, 2]
   :return: Numpy массив формы [N, 2] точек на прямой
   """

   x = points[:, 0]
   y = points[:, 1]
   # Для метода наименьших квадратов нам нужно, чтобы X был матрицей,
   # в которой первый столбец - единицы, а второй - x координаты точек
   X = np.vstack((np.ones(x.shape[0]), x)).T
   normal_matrix = np.dot(X.T, X)
   moment_matrix = np.dot(X.T, y)
   # beta_hat это вектор [перехват, наклон], рассчитываем его в
   # в соответствии с формулой.
   beta_hat = np.dot(np.linalg.inv(normal_matrix), moment_matrix)
   intercept = beta_hat[0]
   slope = beta_hat[1]
   # Теперь, когда мы знаем параметры прямой, мы можем
   # легко вычислить y координаты точек на прямой.
   y_hat = intercept + slope * x
   # Соберем x и y в единую матрицу, которую мы собираемся вернуть
   # в качестве результата.
   points_hat = np.vstack((x, y_hat)).T

   return points_hat

Посмотрим, что мы получили на выходе.

Выглядит в точности так, как мы и хотели, — линия прямая и максимально подогнана под набор точек. Однако у метода наименьших квадратов есть один существенный недостаток — он не работает в случае, если траектория вертикальная. В таком случае наклон прямой обращается в бесконечность, и программа или упадет, или выдаст NaN — «не число». Таким образом, данный подход не может претендовать на продуктовый уровень.

Подход 2: RANSAC

RANSAC — метод достижения консенсуса на основе случайных выборок, предложенный в 1981 году.

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

где p — вероятность успеха, ω — доля выбросов, n — размер суппорта. В нашем случае размер суппорта равен двум точкам для прямой в двумерном пространстве. Для прямой в трехмерном пространстве n было бы тоже равно двум, а для трехмерной плоскости в 3D — уже трем точкам. В примере кода количество итераций получилось 16, что является вполне небольшим значением. Я предлагаю читателю далее следовать за комментариями в коде.

def ransac(points: np.ndarray,
          min_inliers: int = 4,
          max_distance: float = 0.15,
          outliers_fraction: float = 0.5,
          probability_of_success: float = 0.99) -> Optional[np.ndarray]:
   """
   RANdom SAmple Consensus метод нахождения наилучшей
   аппроксимирующей прямой.

   :param points: Входной массив точек формы [N, 2]
   :param min_inliers: Минимальное количество не-выбросов
   :param max_distance: максимальное расстояние до поддерживающей прямой,
                        чтобы точка считалась не-выбросом
   :param outliers_fraction: Ожидаемая доля выбросов
   :param probability_of_success: желаемая вероятность, что поддерживающая
                                  прямая не основана на точке-выбросе
   :param axis: Набор осей, на которых рисовать график
   :return: Numpy массив формы [N, 2] точек на прямой,
            None, если ответ не найден.
   """

   # Давайте вычислим необходимое количество итераций
   num_trials = int(math.log(1 - probability_of_success) /
                    math.log(1 - outliers_fraction**2))

   best_num_inliers = 0
   best_support = None
   for _ in range(num_trials):
       # В каждой итерации случайным образом выбираем две точки
       # из входного массива и называем их "суппорт"
       random_indices = np.random.choice(
           np.arange(0, len(points)), size=(2,), replace=False)
       assert random_indices[0] != random_indices[1]
       support = np.take(points, random_indices, axis=0)

       # Здесь мы считаем расстояния от всех точек до прямой
       # заданной суппортом. Для расчета расстояний от точки до
       # прямой подходит функция векторного произведения.
       # Особенность np.cross в том, что функция возвращает только
       # z координату векторного произведения, а она-то нам и нужна.
       cross_prod = np.cross(support[1, :] - support[0, :],
                             support[1, :] - points)
       support_length = np.linalg.norm(support[1, :] - support[0, :])
       # cross_prod содержит знаковое расстояние, поэтому нам нужно
       # взять модуль значений.
       distances = np.abs(cross_prod) / support_length

       # Не-выбросы - это все точки, которые ближе, чем max_distance
       # к нашей прямой-кандидату.
       num_inliers = np.sum(distances < max_distance)
       # Здесь мы обновляем лучший найденный суппорт
       if num_inliers >= min_inliers and num_inliers > best_num_inliers:
           best_num_inliers = num_inliers
           best_support = support

   # Если мы успешно нашли хотя бы один суппорт,
   # удовлетворяющий всем требованиям
   if best_support is not None:
       # Спроецируем точки из входного массива на найденную прямую
       support_start = best_support[0]
       support_vec = best_support[1] - best_support[0]
       # Для расчета проекций отлично подходит функция
       # скалярного произведения.
       offsets = np.dot(support_vec, (points - support_start).T)
       proj_vectors = np.outer(support_vec, offsets).T
       support_sq_len = np.inner(support_vec, support_vec)
       projected_vectors = proj_vectors / support_sq_len
       projected_points = support_start + projected_vectors

   return projected_points

Что мы получим, запустив код?

Полученная прямая хорошо аппроксимирует набор точек. Черными квадратами отмечены точки, которые были выбраны в качестве суппорта.

В чем недостатки метода RANSAC? В том, что он не гарантирует успешного результата. Вероятность того, что генератор случайных чисел будет выбирать исключительно выбросы или просто неудачные точки, например очень близкие друг к другу, достаточно мала. Тем не менее, в применениях, требующий высокой надежности, например, в беспилотных автомобилях, RANSAC может не подойти. Также RANSAC не подойдет там, где доступ к генератору случайных или псевдослучайных чисел затруднен. Помимо того у RANSAC есть много параметров, которые надо подбирать под задачу. В частности, max_distance довольно сложно оценить, не потратив длительное время на изучение интенсивности шума.

Подход 3: PCA

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

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

def pca(points: np.ndarray, axis: Optional[Any] = None) -> np.ndarray:
   """
   Метод главных компонент (PCA) оценки направления
   максимальной дисперсии облака точек.

   :param points: Входной массив точек формы [N, 2]
   :param axis: Набор осей, на которых рисовать график
   :return: Numpy массив формы [N, 2] точек на прямой
   """

   # Найдем главные компоненты.
   # В первую очередь нужно центрировать облако точек, вычтя среднее
   mean = np.mean(points, axis=0)
   centered = points - mean
   # Функция вычисления собственных значений и векторов np.linalg.eig
   # требует ковариационную матрицу в качестве аргумента.
   cov = np.cov(centered.T)
   # Теперь мы можем посчитать главные компоненты, заданные
   # собственными значениями и собственными векторами.
   eigenval, eigenvec = np.linalg.eig(cov)
   # Мы хотим параметризовать целевую прямую в координатной системе,
   # заданной собственным вектором, собственное значение которого
   # наиболее велико (направление наибольшей вариативности).
   argmax_eigen = np.argmax(eigenval)
   # Нам понадобятся проекции входных точек на наибольший собственный
   # вектор.
   loc_pca = np.dot(centered, eigenvec)
   loc_maxeigen = loc_pca[:, argmax_eigen]
   max_eigenval = eigenval[argmax_eigen]
   max_eigenvec = eigenvec[:, argmax_eigen]
   # Ре-параметризуем прямую, взяв за начало отрезка проекции
   # первой и последней точки на прямую.
   loc_start = mean + max_eigenvec * loc_maxeigen[0]
   loc_final = mean + max_eigenvec * loc_maxeigen[-1]
   linspace = np.linspace(0, 1, num=len(points))
   # Получаем позиции точек, которые идут с одинаковым интервалом,
   # таким образом удаляя шум измерений и вдоль траектории движения.
   positions = loc_start + np.outer(linspace, loc_final - loc_start)

   return positions

Запустив код, получаем следующий график.

PCA также справился с задачей! При этом этот метод лишен недостатка метода наименьших квадратов: PCA отлично аппроксимирует вертикальные линии. Также PCA всегда даст результат в отличие от вероятностного RANSAC. Однако, в отличие от RANSAC метод PCA чувствителен к сильным выбросам в локациях точек.

Заключение

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

Из трех решений мой личный фаворит — метод главных компонент PCA.

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

Спасибо за внимание и удачи на собеседованиях!