Источник: Nuances of Programming
Наука о данных и разложение матриц
Специалистам по данным стоит хорошо знать несколько разложений матриц, потому что они помогают находить методы для актуальных вычислений и оценки результатов для моделей и алгоритмов, с которыми мы работаем. В некоторых случаях конкретной формой разложения является алгоритм (например, в методе главных компонент и сингулярном разложении). И во всех случаях разложение матриц помогает развивать интуицию и способность к анализу.
QR—разложение — одно из самых полезных. Ему находится множество важных применений в науке о данных, статистике и анализе данных. Одно из подобных применений — вычисление решения задачи наименьших квадратов.
Содержание статьи
- Повторение проблемы наименьших квадратов.
- Описание с QR-разложения.
- Решение задачи наименьших квадратов с помощью QR-разложения.
- Реализация вычисления QR-разложения на R и C++ и сравнение реализаций.
Задача наименьших квадратов
QR-разложение позволяет вычислить решение задачи наименьших квадратов. Подчеркиваю, именно вычислить, поскольку обычный метод наименьших квадратов дает нам закрытое решение в форме нормальных уравнений. Это замечательно, но, если нам нужно найти актуальное числовое решение, этот метод не подойдет.
Вспомним задачу наименьших квадратов. Нужно решить уравнение ниже:
Проблема состоит в том, что не существует решения для β, потому что обычно, если у нас больше наблюдений, чем переменных, X не имеет обратного значения, следовательно, вычисление ниже невозможно:
Вместо этого попробуем найти некоторое β̂, решающее уравнение неидеально, но с минимально возможной ошибкой. Один из способов— минимизировать следующую целевую функцию, являющуюся функцией от β̂.
Минимизация этой суммы квадратов отклонений и дает имя задаче наименьших квадратов. Взятие производных по β̂ и приравнивание к нулю приведут к нормальным уравнениям и предоставят решение в замкнутой форме.
Это один из способов. Но можно использовать линейную алгебру. И вот здесь QR-разложение подойдет как нельзя лучше.
QR-разложение
Для начала давайте опишем это разложение. QR-разложение позволяет отобразить матрицу как произведение двух отдельных матриц Q и R.
Q — ортогональная матрица, а R — треугольная матрица.
Это означает, что:
Так как R квадратная, до тех пор, пока диагональные элементы не равны нулю, она также обратима. Если столбцы X линейно независимы, это всегда будет верным. Хотя, если в данных есть коллинеарность, все же будут возникать проблемы. Тем не менее — и в этом суть QR-разложения — прямоугольная и необратимая X может быть выражена как две обратимые матрицы! И вот это уже имеет смысл.
Решение задачи наименьших квадратов с помощью QR-разложения.
Теперь, зная, что из себя представляет QR-разложение, решим задачу наименьших квадратов следующим образом:
Следовательно:
Это означает, что все, что нужно сделать — это найти матрицу, обратную R, транспонировать Q и вычислить их произведение. Мы получим коэффициенты обычного метода наименьших квадратов. Нам даже не нужно вычислять ковариационную матрицу и обратную ей, что происходит в решении обычным методом наименьших квадратов.
Реализация QR-разложения
Чтобы найти QR-коэффициенты матрицы, сначала найдем Q, используя процесс Грама-Шмидта, затем просто умножим исходную матрицу на транспонированную Q, чтобы найти R. Далее применим QR-разложение, используя функции, внедренные в R и C++. Позднее мы заглянем внутрь этих функций, чтобы рассмотреть весь процесс в деталях.
Вычисление коэффициентов
Я использую две функции: myQRR и myQRCpp, применяющие процесс Грама-Шмидта для QR-разложения. Одна функция написана на R, а другая на C++, она загружается в среду R через Rcpp. В конце статьи я сравню их производительность.
Начнем с маленького примера, в котором смоделируем y и X, а затем решим уравнение, используя QR-разложение. Также мы сможем провести двойную проверку QR-разложения — работает ли оно и возвращает ли смоделированный X. Вот смоделированные переменные ответа:
Вот данные, которые мы используем для определения коэффициентов наименьших квадратов. В нашем распоряжении 3 переменные:
Теперь, чтобы найти Q и R, используем myQRCpp.
- 1. Видим, что R действительно верхнетреугольная матрица.
2. Убеждаемся в ортогональности Q.
3. И в том, что QR действительно возвращает исходную матрицу X.
Теперь вычисляем актуальные коэффициенты.
Чтобы проверить, верен ли ответ, сравним вычисленное значение β̂ с тем, что выдает функция lm.
Мы получили в точности то же решение, что и для рассчитанных коэффициентов.
Реализация QR-разложения
Процесс Грама-Шмидта
Процесс Грама-Шмидта — это метод вычисления ортогональной матрицы Q, которая состоит из ортогональных или независимых единичных векторов и занимает такое же пространство, что и матрица X.
- В качестве начального шага алгоритм включает в себя выбор вектора столбца X, скажем, x1 = u1.
- Затем находим вектор, ортогональный u1, проецируя на него следующий столбец X, скажем, x2, и вычитая из него проекцию u2 = x2 − proj u1x2. Теперь у нас есть набор из двух ортогональных векторов.
- Следующий шаг — проделать то же самое, но вычитая сумму проекций на каждый вектор из набора ортогональных векторов uk.
Выразим это следующим образом:
Получив полный набор ортогональных векторов, просто делим каждый вектор на его нормаль и помещаем их в матрицу:
Зная Q, легко вычисляем R:
Реализация в R и C++
Конечно, в R существуют встроенные функции, осуществляющие QR-разложение. Так как алгоритм Грама-Шмидта является итеративным по своей природе, я решил реализовать его на C++ , являющемся хорошим инструментом для подобных задач, и сравнить его с подобной функцией R. Вот как выглядит моя версия R:
Очень точно. Внутри цикла for есть цикл while, и вызываемая проекционная функция также является функцией, написанной на R.
А вот как выглядит моя версия C++. Логика по сути та же, за исключением того, что есть другой цикл for, нормализующий ортогональные столбцы.
В дополнение к двум функциям выше у меня есть третья, идентичная R, только она называется projC, а не proj. Я назвал ее myQRC (projC написана на C++, а proj — на R). В противном случае у нас есть чистая C++ функция myQRCpp и чистая R функция myQR.
Чтобы сравнить, как быстро эти три функции выполняют QR-разложение, я поместил их в функцию QR_comp, которая вызывает функции и измеряет время выполнения каждой с одним и тем же аргументом матрицы.
Их производительность можно сравнить по сетке из случайных матриц n на m. Эти матрицы генерируются при вызове функции QR_comp.
Наконец, визуализируем.
Очевидно, что чем больше задействован C++ , тем быстрее вычисление QR-разложения. Функция C++ решает его менее, чем за минуту, для матриц размером до 250 столбцов на 3000 строк или 600 столбцов на 500. R-функция в 2–3 раза медленнее.
Заключение
QR — это просто разложение матрицы, а метод наименьших квадратов просто одно из применений QR. Надеюсь, обсуждение выше демонстрирует, насколько важна и полезна линейная алгебра для науки о данных.
Читайте также:
Читайте нас в телеграмме и vk
Перевод статьи Ben Denis Shaffer: QR Matrix Factorization