Как реализовать линейную регрессию

В предыдущем уроке мы использовали реализацию scikit learn для вычисления линейной регрессии методом наименьших квадратов на наборе данных о диабете. Сегодня мы рассмотрим, как написать собственную реализацию.

Начало работы

from sklearn import datasets, linear_model, metrics
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import PolynomialFeatures
import math, scipy, numpy as np
from scipy import linalg

np.set_printoptions(precision=6)

data = datasets.load_diabetes()

feature_names=['age', 'sex', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6']

trn,test,y_trn,y_test = train_test_split(data.data, data.target, test_size=0.2)

trn.shape, test.shape

# ((353, 10), (89, 10))

def regr_metrics(act, pred):
    return (math.sqrt(metrics.mean_squared_error(act, pred)), 
    metrics.mean_absolute_error(act, pred))

Как это реализовано в sklearn

Как это делает sklearn? Проверяя исходный код, можно увидеть, что в плотном случае он вызывает scipy.linalg.lstqr, который вызывает метод LAPACK:

Опции — это 'gelsd', 'gelsy', 'gelss'. По умолчанию 'gelsd' является хорошим выбором, но 'gelsy' быстрее на многих задачах. 'gelss' используется по историческим причинам. Обычно он медленнее, но использует меньше памяти.

'gelsd' использует SVD и метод деления пополам;
'gelsy' использует QR-разложение;
'gelss' использует SVD.

Scipy разреженная линейная регрессия методом наименьших квадратов

Мы не будем подробно описывать разреженную версию метода наименьших квадратов. Если вам интересно, обратитесь к следующей информации:

Scipy разреженный метод наименьших квадратов использует итерационный метод, называемый двойной диагонализацией Паге и Сайндса.

Исходный код Scipy для разреженной линейной регрессии методом наименьших квадратов: предварительная обработка — это другой способ уменьшить количество итераций. Если возможно эффективно решить связанную систему M*x = b, где M приближает A некоторым полезным способом (например, M-A имеет низкий ранг или его элементы относительно элементов A малы), то LSQR может сходиться быстрее в системе A*M(inverse)*z = b. Затем можно восстановить x путём решения M * x = z.

Если A симметричен, то не следует использовать LSQR! Альтернативные варианты — симметричный сопряжённый градиентный метод (cg) и/или SYMMLQ. SYMMLQ — это реализация сопряжённого градиента для симметричных матриц, которая подходит для любого симметричного A и сходится быстрее, чем LSQR. Если A положительно определён, существуют другие реализации сопряжённого градиента, которые требуют меньше работы на каждой итерации (но требуют того же количества итераций).

`linalg.lstsq`

sklearn добавляет для нас константу (потому что для изучаемой нами прямой y перехват может быть не равен 0). Теперь нам нужно сделать это вручную:

trn_int = np.c_[trn, np.ones(trn.shape[0])]
test_int = np.c_[test, np.ones(test.shape[0])]

Поскольку linalg.lstq позволяет нам указать, какой метод LAPACK мы хотим использовать, давайте попробуем их и проведём некоторые временные сравнения:

%timeit coef, _, _, _ = linalg.lstsq(trn_int, y_trn, lapack_driver="gelsd")

# 290 µs ± 9.24 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit coef, _, _, _ = linalg.lstsq(trn_int, y_trn, lapack_driver="gelsy")

# 140 µs ± 91.7 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit coef, _, _, _ = linalg.lstsq(trn_int, y_trn, lapack_driver="gelss")

# 199 µs ± 228 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Простой метод

Вспомните, что мы ищем $\hat x$ , чтобы минимизировать:

$\big\vert\big\vert Ax - b \big\vert\big\vert_2$

Другой способ подумать об этом заключается в том, что нас интересует место, наиболее близкое к b в пространстве, охватываемом A. Это проекция b на A. Поскольку $b - A \hat x$ должно быть перпендикулярно пространству, охватываемому A, мы видим:

$A^T (b - A\hat{x}) = 0$

Используя $A^T$ , потому что мы хотим умножить A и $b - A\hat{x}$ на каждую колонку.

Это даёт нам нормальное уравнение:

$x = (A^TA)^{-1}A^T b$

def ls_naive(A, b):
     return np.linalg.inv(A.T @ A) @ A.T @ b

%timeit coeffs_naive = ls_naive(trn_int, y_trn)

# 45.8 µs ± 4.65 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

coeffs_naive = ls_naive(trn_int, y_trn)
regr_metrics(y_test, test_int @ coeffs_naive)

# (57.94102134545707, 48.053565198516438)

Нормальное уравнение (Cholesky)

Нормальное уравнение:

$A^TA x = A^T b$

Если у A полный ранг, псевдообратная матрица $(A^TA)^{-1}A^T$ является квадратной эрмитовой положительно определённой матрицей. Стандартным методом решения такой системы является разложение Холецкого, которое находит верхнюю треугольную матрицу R, удовлетворяющую $A^TA = R^TR$ .

Следующие шаги основаны на алгоритме Trefethen 11.1:

A = trn_int

b = y_trn

AtA = A.T @ A
Atb = A.T @ b

Предупреждение: Для Cholesky Numpy и Scipy по умолчанию используют разные верхние/нижние треугольники.

R = scipy.linalg.cholesky(AtA)

np.set_printoptions(suppress=True, precision=4)
R

'''
array([[  0.9124,   0.1438,   0.1511,   0.3002,   0.2228,   0.188
``` **Перевод текста на русский язык:**

Проверка нашего разложения:

```py
np.linalg.norm(AtA - R.T @ R)

# 4.5140158187158533e-16

Рисунок A^T A x = A^T b, R^T R x = A^T b, R^T w = A^T b, Rx = w (img/tex-82513c6fe20ef419fc103891436353ab.gif).

w = scipy.linalg.solve_triangular(R, Atb, lower=False, trans='T')

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

np.linalg.norm(R.T @ w - Atb)

# 1.1368683772161603e-13

coeffs_chol = scipy.linalg.solve_triangular(R, w, lower=False)

np.linalg.norm(R @ coeffs_chol - w)

# 6.861429794408013e-14

def ls_chol(A, b):
    R = scipy.linalg.cholesky(A.T @ A)
    w = scipy.linalg.solve_triangular(R, A.T @ b, trans='T')
    return scipy.linalg.solve_triangular(R, w)

%timeit coeffs_chol = ls_chol(trn_int, y_trn)

# 111 µs ± 272 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

coeffs_chol = ls_chol(trn_int, y_trn)
regr_metrics(y_test, test_int @ coeffs_chol)

# (57.9410213454571, 48.053565198516438)
## QR разложение

Рисунок Ax = b, A = Q R, Q R x = b, Rx = Q^T b (img/tex-d45f0a6200f855ffbed7dbbf8dc6094b.gif).

```py
def ls_qr(A,b):
    Q, R = scipy.linalg.qr(A, mode='economic')
    return scipy.linalg.solve_triangular(R, Q.T @ b)

%timeit coeffs_qr = ls_qr(trn_int, y_trn)

# 205 µs ± 264 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

coeffs_qr = ls_qr(trn_int, y_trn)
regr_metrics(y_test, test_int @ coeffs_qr)

# (57.94102134545711, 48.053565198516452)
### SVD

Рисунок Ax = b, A = U \Sigma V, \Sigma V x = U^T b, \Sigma w = U^T b, x = V^T w (img/tex-4b62075d74bf4d4e462ede239a6dea12.gif).

SVD даёт псевдообратную матрицу.

```py
def ls_svd(A,b):
    m, n = A.shape
    U, sigma, Vh = scipy.linalg.svd(A, full_matrices=False)
    w = (U.T @ b)/ sigma
    return Vh.T @ w

%timeit coeffs_svd = ls_svd(trn_int, y_trn)

# 1.11 ms ± 320 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit coeffs_svd = ls_svd(trn_int, y_trn)

# 266 µs ± 8.49 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

coeffs_svd = ls_svd(trn_int, y_trn)
regr_metrics(y_test, test_int @ coeffs_svd)

# (57.941021345457244, 48.053565198516687)
### Минимальное квадратичное отклонение регрессии с использованием методов случайного Sketching

[Линейный Sketching](http://researcher.watson.ibm.com/researcher/files/us-dpwoodru/journal.pdf) (Woodruff)

+ Выборка r×n случайной матрицы S, r << n
+ Вычисление SA и Sb
+ Нахождение точного решения Ax = b для SA x = Sb
### Сравнение времени

```py
import timeit
import pandas as pd

def scipylstq(A, b):
    return scipy.linalg.lstsq(A,b)[0]

row_names = ['Normal Eqns- Naive',
             'Normal Eqns- Cholesky', 
             'QR Factorization', 
             'SVD', 
             'Scipy lstsq']

name2func = {'Normal Eqns- Naive': 'ls_naive', 
             'Normal Eqns- Cholesky': 'ls_chol', 
             'QR Factorization': 'ls_qr',
             'SVD': 'ls_svd',
             'Scipy lstsq': 'scipylstq'}

*В тексте запроса присутствуют фрагменты кода на языке Python, которые были сохранены без изменений.* **Почему Cholesky быстрее:**

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

**QR-разложение: пример оптимального случая:**

Далее в тексте приводится пример использования QR-разложения для решения системы линейных уравнений. Автор сравнивает результаты решения с использованием QR-разложения и других методов, таких как метод нормальных уравнений, SVD (сингулярное разложение) и scipy lstsq. Результаты показывают, что QR-разложение является эффективным методом для решения линейных систем уравнений.

**Низкий ранг:**

Затем автор рассматривает случай, когда матрица имеет низкий ранг. Он сравнивает различные методы решения линейных систем с матрицами низкого ранга и приходит к выводу, что метод нормальных уравнений и Cholesky могут быть нестабильными для таких матриц.

**Сравнение:**

Автор сравнивает время выполнения и точность различных методов для разных размеров матриц и делает вывод, что для большинства случаев метод QR-разложения является наиболее предпочтительным.

**Выводы:**

На основе проведённого анализа автор приходит к заключению, что метод Cholesky может быть наиболее быстрым для определённых типов задач, но он ограничен симметричными положительно определёнными матрицами. Для большинства практических задач рекомендуется использовать метод QR-разложения.

OSCHINA-MIRROR/wizardforcel-fastai-num-linalg-v2-zh

Как реализовать линейную регрессию

Начало работы

Как это реализовано в sklearn

Scipy разреженная линейная регрессия методом наименьших квадратов

`linalg.lstsq`

Простой метод

Нормальное уравнение (Cholesky)

Опубликовать ( 0 )

Введение

Обновления

Участники

Недавние действия

OSCHINA-MIRROR/wizardforcel-fastai-num-linalg-v2-zh .gitee-modal { min-width: 500px !important; } .gitee-modal .close { right: 0.6rem !important; }

Как реализовать линейную регрессию

Начало работы

Как это реализовано в sklearn

Scipy разреженная линейная регрессия методом наименьших квадратов

linalg.lstsq

Простой метод

Нормальное уравнение (Cholesky)

Опубликовать ( 0 )

Введение

Обновления

Участники

Недавние действия

OSCHINA-MIRROR/wizardforcel-fastai-num-linalg-v2-zh

`linalg.lstsq`