Погрузитесь во внутреннее устройство linalg.polyfit () от NumPy.

Сегодня программа линейной регрессии в машинном обучении может быть написана всего в несколько строк настолько легко, что мы не сможем оценить глубину прогресса в математике и компьютерных системах, которая делает эту простоту для нас. За мощными numpy и scipy библиотеками от нас скрывается огромное количество математических и числовых сложностей. В этой статье мы глубоко погрузимся в внутреннюю реализацию линейной регрессии, чтобы лучше понять, что именно происходит под покровом. В процессе мы откроем для себя множество элегантных методов линейной алгебры, сложности вызова функций Fortran из Python и т. Д.

В задаче линейной регрессии нам нужно решить Ax = b. Здесь (A, b) - это набор наблюдаемых данных, и нам нужно найти параметры x, которые соответствуют этим данным. Существует бесчисленное множество математических способов решения системы линейных уравнений. Нас интересуют внутренние математические методы, лежащие в основе пакетов машинного обучения, таких как numpy / scipy.

Начнем с примера использования numpy.linalg.polyfit (). (Вы можете скачать код Python на Github линейная регрессия с использованием numpy polyfit).

Вот визуализация результата. «*» - это наблюдаемые данные (или сгенерированные в нашем случае), а линия проводится по параметрам, которые были вычислены с помощью polyfit. Ясно, что polyfit хорошо справляется с линейной регрессией.

Большинство статей в Интернете останавливаются на этом и не углубляются. Однако глубокое погружение и попытки понять, что происходит, поможет нам оценить, как было вычислено решение. Итак, давайте углубимся в np.linalg.polyfit ().

Первое, что мы сделаем, это то, что np.linalg.polyfit () не полностью реализован в python. Суть его вычислений происходит в библиотеке, написанной с помощью оптимизированного кода FORTRAN! Посмотрим, как это так. Вызов polyfit () проходит через следующие блоки, как показано ниже. Я вырезал и вставил соответствующие вызовы из множества файлов, чтобы вы могли увидеть, как все это сшивается:

В блоке 2 вызов polyfit () создаст матрицу Вандермонда через вызов numpy.linalg.polyvander (), специальной матрицы, в которой столбцы расположены в геометрической прогрессии. Предполагая, что пользователь хочет подогнать под полином 5-й степени, строятся столбцы до x ⁵. Обратите внимание, что пользователю нужно только указать одномерный столбец со значениями x, а polyfit () генерирует матрицу Вандермонда. Тем самым он готовится к полиномиальной линейной регрессии для этого уравнения:

Исходная матрица A преобразовалась следующим образом:

Обратите внимание, что это прямоугольная матрица с m ›n (градус). В зависимости от количества наблюдений данных m может быть очень большим, тогда как n обычно является небольшим числом. Как только эта матрица настроена, начинается подготовка к вызову реализации матрицы в LAPACK (библиотеке линейной алгебры).

Библиотека LAPACK была написана на Fortran 90 еще в 1980-х и 1990-х годах, но регулярно обновляется. Последняя версия LAPACK с открытым исходным кодом - 3.9.0 - была выпущена в 2019 году. Она также оптимизирована такими поставщиками процессоров, как Intel, Apple и т. Д. Библиотека Fortran скомпилирована в динамическую библиотеку на основе архитектуры платформы и поставляется с операционная система. Кроме того, когда вы устанавливаете numpy, 2 общих объекта (.so) также идут вместе с ним, чтобы упростить вызовы C API. Чтобы просмотреть эти файлы, перейдите по ‹пути к вашей установке python› /lib/python3.7/site-packages/numpy/linalg/, и вы увидите следующее:

Здесь два важных файла:

  1. _umath_linalg.cpython-37m-darwin.so. : Этот файл преобразует структуры данных Python в структуры данных C, чтобы можно было вызвать интерфейс C для LAPACK.
  2. lapack_lite.cpython-37m-darwin.so: этот файл является интерфейсом C для подпрограмм LAPACK. Он вызывает процедуры LAPACK и преобразует возвращаемые значения обратно в структуры данных Python. Фактические процедуры LAPACK предоставляются ОС. Например, в Apple macOS версию, оптимизированную для платформы, можно найти по адресу /System/Library/Frameworks/vecLib.framework/libLAPACK.dylib

Давайте продолжим, как происходит поток вызовов.

  1. Основная точка входа в вычисление решения методом наименьших квадратов - np.linalg.lstsq (), как показано в блоке 3.
  2. В блоке 4 на приведенной выше диаграмме np.linalg.lstsq () вызывает интерфейс C в _umath_linalg.lstsq_m ().
  3. Блоки 5–8 - это шаблоны для создания интерфейса python-C. Например, в блоке 6 call_ @ lapack_func @ становится call_sgelsd (), определенным в блоке 7. В блоке 7 LAPACK (@lapack_func @) (…) становится вызовом sgelsd_ (…). Это внешняя функция, реализованная в libLAPACK.dylib, предоставленном Apple.

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

  1. QR или LQ факторизация
  2. Полная ортогональная факторизация
  3. Разложение по сингулярным значениям (SVD) факторизация
  4. Разделяй и властвуй СВД

numpy.linalg.lstsq () решил использовать методы SVD "разделяй и властвуй". Откуда нам знать? Посмотрите в блоке 6, и вы найдете эти функции, которые указаны для генерации шаблонов. Итак, выбранные функции - SGELSD, CGELSD, ZGELSD и DGELSD.

Как работает СВД разделяй и властвуй? Высокоуровневое описание алгоритма приведено в блоке 9 и является частью комментариев к коду, но вы можете посмотреть все детали в исходной статье. Он показывает, как именно вычисляется решение методом наименьших квадратов. Вот аннотация статьи:

Чтобы прочитать вышеупомянутую статью, вот основные концепции линейной алгебры, которые необходимо знать: 1) Разложение по сингулярным значениям (SVD) 2) QR-факторизация 3) Псевдообратная матрица Мура-Пенроуза. На YouTube есть несколько хороших видеороликов по линейной алгебре, но я бы порекомендовал Лекции Гиберта Стрэнга по линейной алгебре (MIT). Не менее сказочен его учебник. (Я фиксирую все основные математические концепции, необходимые для машинного обучения здесь). Из документа мы видим, что текущая реализация метода наименьших квадратов, то есть SGELSD, в LAPACK в 9 раз быстрее, чем его предшественник для двудиагональных матриц! По этой причине SGELSD на первом этапе преобразует входную матрицу в двухдиагональную матрицу, чтобы мы могли извлечь выгоду из ускорения. После вычисления SVD вычислить решение методом наименьших квадратов несложно. Мы вычисляем псевдообратную матрицу Мура-Пенроуза, а затем выводим из нее решение методом наименьших квадратов. Последний шаг также подробно описан в статье.

Таким образом, мы проследили полный путь numpy.linalg.polyfit () вплоть до основного кода Fortran. Мы обнаружили, что он использует важные методы численной линейной алгебры, такие как разложение по сингулярным числам. Мы также просмотрели исходную статью, в которой были предложены алгоритмы, которые использует numpy.linalg.polyfit (), и обнаружили, что эти алгоритмы значительно быстрее, чем предыдущие реализации. Для беспрепятственного взаимодействия кода Fortran с Python была разработана тщательная реализация интерфейсов C и Python API, чтобы процесс не вылетел из строя. Чтобы сделать это надежным, требуются высококачественные методы разработки программного обеспечения. Вдобавок поставщики ОС предоставили тщательно настроенные библиотеки, чтобы сделать реализацию LAPACK более эффективной для своей процессорной архитектуры. Таким образом, достижения в области математики, разработки программного обеспечения, компьютерных архитектур и открытого исходного кода сделали numpy.linalg.polyfit () возможным и быстрым.

Надеюсь, вам понравилось глубокое погружение в nump.linalg.polyfit ().

Использованная литература:

  1. LAPACK - Пакет линейной алгебры, Руководство пользователя LAPACK, Функции наименьших квадратов LAPACK
  2. Исходный код Numpy на Github
  3. Справочник разработчика для Intel® Math Kernel Library 2020 - C Поставщики процессоров, такие как Intel, предлагают оптимизированные библиотеки для LAPACK и других математических библиотек.
  4. Библиотека ядра Intel Math
  5. VecLib: реализация LAPACK для Apple macOS
  6. Вы можете скачать исходный код LAPACK здесь, чтобы увидеть файлы Fortran.
  7. Исходный код SGELSD Fortran
  8. Эффективное вычисление разложения сингулярных чисел с приложениями к задачам наименьших квадратов.