При социализме средства производства должны находиться в общественной собственности. Это позволит распределять ресурсы справедливо и эффективно в интересах всего общества.
Однако совместное потребление в капиталистических странах - это другое. Здесь речь идёт не о производстве, а о потреблении товаров и услуг. Люди делятся друг с другом вещами, которыми владеют, или совместно арендуют что-то для временного использования. Такие инициативы могут быть полезны для отдельных людей. Но они не меняют капиталистических отношений собственности и производства. Поэтому я бы не стал говорить, что совместное потребление имеет прямое отношение к социалистическим идеям. Это скорее попытки смягчить некоторые негативные стороны капитализма.
(C) ЛЕНИН, модель ИИ от группы «Свободное время».
Стандартная реализация языка Python написана на языке C. Это значит, что каждый объект Python — просто искусно замаскированная структура данных языка C, содержащая не только значение, но и другую информацию. Например переменная x=1000 в Python не просто целое число, а указатель на место в памяти, где хранится составная структура, содержащая несколько значений (счётчик ссылок, код типа, размер следующих элементов, фактическое значение). Это влечёт накладные расходы, что становится особенно заметно при манипуляции с большим количеством таких объектов. Поэтому, хотя в Python и есть возможность описывать матрицы с помощью стандартных типов данных, например списков, но на практике чаще используют сторонние библиотеки, среди которых наиболее часто используемой является NumPy.
NumPy — это фундаментальная библиотека Python для научных вычислений и анализа данных. Он обеспечивает эффективную реализацию многомерных массивов и матриц, а также обширную коллекцию математических функций высокого уровня для работы с этими массивами. Для решения задач линейной алгебры NumPy предоставляет структуру данных n-мерного массива, называемую ndarray, а также векторизованные реализации наиболее распространённых операций, которые позволяют нам эффективно работать с большими наборами матричных данных
В отличие от списков Python, библиотека NumPy ограничивается массивами, содержащими элементы одного типа. Зато, массивы NumPy могут явным образом описываться как многомерные. Также NumPy предоставляет для многих типов операций удобный интерфейс для компилируемой процедуры со статической типизацией (векторизованной операции). Это значительно ускоряет вычисления.
Перед тем как начать работу, NumPy необходимо установить (если это не было сделано раньше). Для этого в консоли нужно выполнить команду:
pip install numpy
Затем, в проекте нужно импортировать эту библиотеку numpy:
import numpy as np
np- это псевдоним, обращаясь к которому можно вызывать функции из библиотеку numpy.
Создание массивов и матриц
Основной структурой данных в NumPy является класс ndarray. Мы можем создавать новые ndarrays из списков или кортежей в Python:
vector = np.array([1, 2, 3])
matrix = np.array([(1, 2, 3), (4, 5, 6)])
Конструктор ndarray преобразует вложенные последовательности, такие как списки списков, в многомерные массивы.
Мы также можем создавать конкретные формы матриц с помощью функций нулей и единиц:
# 3x3 матрица, заполненная нулями
zero_matrix = np.zeros((3,3))
# 4x2 матрица, заполненная единицами
one_matrix = np.ones((4,2))
NumPy предлагает методы arange и linspace для генерации числовых диапазонов:
# Массив значений от 0 до 9
arange_array = np.arange(10)
# 10 чисел равномерно распределённых между 0 и 5
linspace_array = np.linspace(0, 5, 10)
Основные свойства и атрибуты массивов NumPy
Атрибут shape возвращает размеры массива, а dtype возвращает тип данных:
matrix = np.array([(1, 2, 3), (4, 5, 6)])
print(matrix.shape) # (2, 3)
print(matrix.dtype) # int64
Массивы NumPy имеют полезные атрибуты, предоставляющие информацию о данных:
ndim - Количество измерений массива;
size - Общее количество элементов в массиве;
itemsize - Размер хранения каждого элемента в байтах.
Например:
arr = np.arange(10)
print("arrange matrix: ", arr)
print("dimension matrix: ", arr.ndim) # 1
print("size matrix: ", arr.size) # 10
print("size items in matrix:", arr.itemsize) # 8 (64-bit system)
Мы также можем легко получить основные статистические свойства, такие как среднее значение, стандартное отклонение, дисперсия и т. д., для массивов NumPy:
Мы также можем легко получить основные статистические свойства, такие как среднее значение, стандартное отклонение, дисперсия и т. д., для массивов NumPy:
stats = np.random.normal(size=100)
print(stats.mean())
print(stats.std()) # standard deviation
print(stats.var()) # variance
Индексирование и нарезка массива
Массивы NumPy поддерживают знакомый синтаксис срезов Python для извлечения подмассивов:
arr = np.arange(10)
# Первые 5 элементов
print(arr[:5])
# Последние 3 элемента
print(arr[-3:])
Доступ к конкретным элементам можно получить через целочисленные индексы:
matrix = np.array([(1, 2, 3), (4, 5, 6)])
# Получить первую строку
print(matrix[0])
# Получить элемент в строке 1 и столбце 2
print(matrix[1, 2]) # 6
Булевы индексные маски также могут извлекать элементы на основе условий:
rand_ints = np.random.randint(0, 10, 10)
# Элементы больше 5
mask = rand_ints > 5
print(rand_ints[mask])
Основные матричные операции
Давайте посмотрим на некоторые основные математические операции между массивами и скалярами.
Сложение и вычитание
arr = np.arange(5)
print(arr + 5) # Прибавить 5 к каждому элементу
mat1 = np.ones((2, 3))
mat2 = np.ones((2, 3))
print(mat1 + mat2) # Поэлементное сложение
arr = np.arange(5)
print(arr - 2) # Вычитает 2 из каждого элемента
Поэлементное вычитание матрицы происходит аналогично.
Умножение
- Умножение матрицы на число:
arr = np.full((5,), 5)
print(arr * 10) # Умножение каждого элемента на 10
- Умножение матриц с использованием функции dot():
mat1 = [[1, 2], [3, 4]]
mat2 = [[5, 6], [7, 8]]
mat1 = np.array(mat1)
mat2 = np.array(mat2)
print(np.dot(mat1, mat2))
- Поэлементное умножение матриц с использованием *:
mat1 = np.ones((2, 3))
mat2 = np.full((2, 3), 2)
print(mat1 * mat2)
Деление
- Скалярное деление:
arr = np.full((5,), 10)
print(arr / 2) # Деление каждого элемента на 2
- Поэлементное деление матриц:
mat1 = np.full((3, 3), 6)
mat2 = np.full((3, 3), 3)
print(mat1 / mat2)
Возведение в степень.
- Скалярное возведение массива в степень:
arr = np.array([1, 2, 3])
print(2**arr) # [2, 4, 8]
- Поэлементное возведение в степень с помощью **:
mat = np.full((2, 2), 2)
print(mat**3)
Скалярное произведение
Скалярное произведение — это фундаментальная операция линейной алгебры, которая проецирует два вектора друг на друга. В NumPy мы можем вычислять скалярные произведения с помощью функции dot().
numpy.dot((vector_a, vector_b, out = None) : возвращает скалярное произведение векторов a и b.
v1 = np.array([1, 2])
v2 = np.array([3, 4])
print(np.dot(v1, v2))
Для матриц функция dot() выполняет стандартное умножение матриц:
mat1 = [[1, 2], [3, 4]]
mat2 = [[5, 6], [7, 8]]
mat1 = np.array(mat1)
mat2 = np.array(mat2)
print(np.dot(mat1, mat2))
# [[19 22]
# [43 50]]
Мы также можем вычислить скалярные произведения по определенным осям, используя np.tensordot():
mat1 = np.arange(8).reshape(2, 2, 2)
mat2 = np.arange(8).reshape(2, 2, 2)
print(np.tensordot(mat1, mat2, axes=((1, 0), (0, 1))))
Здесь мы взяли скалярное произведение по заданным комбинациям осей.
Также скалярное произведение матриц можно выполнить с помощью оператора @:
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
B = np.array([[9, 8, 7], [6, 5, 4], [3, 2, 1]])
G = A @ B print("@ Operator:\n", G)
Если первый аргумент является комплексным, для вычисления скалярного произведения используется комплексное сопряжение первого аргумента (здесь vdot() отличается от метода dot()). Он может обрабатывать многомерные массивы, но работает с ними как с плоским массивом.
Результат:
Как происходит вычисление:
Другие функции для произведения матриц:
matmul(): Матричное произведение двух массивов.
inner(): Внутреннее произведение двух массивов — это операция, которая вычисляет скалярное произведение между двумя массивами. Для векторов (одномерных массивов) внутреннее произведение вычисляется как сумма произведений соответствующих элементов массивов. Для многомерных массивов внутреннее произведение вычисляется как сумма произведений по последним осям массива.
outer(): Внешнее произведение двух векторов — это операция, которая создаёт матрицу путём умножения каждого элемента оного вектора на каждый элемент другого вектора. Результатом является матрица, размерность которой равна произведению размерностей исходных веткоров. Внешнее произведение может использоваться для вычисления ковариационной матрицы, создания новых признаков, вычисления сходства между объектами и разложения матрицы на более простые компоненты.
linalg.multi_dot(): Вычислите скалярное произведение двух или более массивов за один вызов функции, автоматически выбирая самый быстрый порядок вычислений.
tensordot(): Вычислите скалярное произведение тензора вдоль заданных осей для массивов >= 1-D.
einsum(): Оценивает соглашение Эйнштейна о суммировании операндов. Соглашение Эйнштейна — это специальный синтаксис, используемый в библиотеке NumPy для выполнения операций над массивами. Оно позволяет указывать, какие оси входных массивов должны быть перемножены и какие должны быть сохранены в результирующем массива. В основе находится идея суммирования по повторящимся индексам.
einsum_path(): Оценивает порядок сокращения с наименьшей стоимостью для выражения einsum, учитывая создание промежуточных массивов.
linalg.matrix_power(): Возводит квадратную матрицу в (целую) степень n.
kron(): Кронекеровское произведение двух массивов — это операция, которая создаёт новый массив путём комбинирования всех возможных пар элементов из двух исходных массивов. Результатом является новый массив, размер которого равен произведению размеров исходных массивов.
Пример:
Результат:
Транспонирование и изменение формы массивов
Транспонирование матрицы Транспонирование матрицы получается путем ее переворачивания по диагонали. Это меняет местами индексы строк и столбцов каждого элемента.
Метод transpose() возвращает представление массива с поменянными местами осями:
mat = np.random.rand(3, 5)
print(mat.transpose()) # Транспонированный вид
Также для транспонирования можно использовать атрибут :
I = A.T
print("Transpose:\n", I)
Чтобы изменить форму массива без изменения данных, используйте reshape():
vector = np.arange(10)
print(vector.reshape(5, 2)) # Reshape to 5x2 matrix
Мы также можем выполнить выравнивание до 1D с помощью метода Flatten() или перевести в непрерывную память с помощью ravel().
Чтобы добавить новую ось, используйте np.expand_dims():
vector = np.arange(5)
col_vector = np.expand_dims(vector, axis=1)
# [[0], # [1], # [2], # [3], # [4]]
Вычисление определителя матрицы
Определитель — это скалярное значение, которое предоставляет важную информацию о матрице, например ее обратимость. В NumPy мы можем найти определитель, используя np.linalg.det():
mat = np.array([[1, 2], [3, 4]])
print(np.linalg.det(mat)) # -2.0
Для неквадратных матриц det() вызывает исключение. Значение определителя полезно для обращения матрицы и вычисления собственных значений.
Нахождение обратной матрицы
Обратная квадратная матрица A обозначается как A^-1 и удовлетворяет следующим условиям: A * A^-1 = E. Где E — единичная матрица.
В NumPy мы можем вычислить обратное следующим образом:
from numpy.linalg import inv
mat = np.array([[1, 2], [3, 4]])
mat_inv = inv(mat)
print(mat_inv)
# [[-2. , 1. ], # [ 1.5, -0.5]]
Это находит уникальную обратную матрицу. Обратите внимание, что не все матрицы имеют обратную. Матрица, не имеющая обратной, называется сингулярной или вырожденной матрицей. Если вы попытаетесь найти обратную такую матрицу в NumPy, это вызовет LinAlgError.
Специальные функции
numpy.trace() : Возвращает сумму по диагоналям массива. Если a двумерный, возвращается сумма по его диагонали с заданным смещением, т.е. сумма элементов a[i,i+offset] для всех i.Если a имеет больше чем два измерения, то оси, заданные осью1 и осью2, используются для определения двумерных подмассивов, трассировки которых возвращаются. Форма результирующего массива такая же, как и у массива a, но с удаленными осью1 и осью2.
Результат:
numpy.linalg.norm(): Матричная или векторная норма. Используется для вычисления нормы вектора или нормы матрицы. Норма вектора — это мера его длины или величины, а норма матрицы — это мера ее размера или протяженности. Функция поддерживает различные типы норм, такие как евклидова норма, манхэттенская норма и максимальная норма. Конкретная норма, подлежащая вычислению, указывается с помощью параметра ord. По умолчанию функция вычисляет евклидову норму (ord=2). Функция принимает на вход объект, подобный массиву, и возвращает скалярное значение, представляющее норму входного массива.
numpy.linalg.cond(): Вычисление числа обусловленности матрицы. Используется для вычисления числа обусловленности матрицы. Число обусловленности — это мера того, насколько чувствительна матрица к изменениям входных данных. Это свидетельствует об устойчивости решения линейной системы уравнений, включающей матрицу. Функция поддерживает различные типы норм, такие как 1-норма, 2-норма и норма бесконечности. Конкретная норма, которая будет использоваться в расчете, указывается с помощью параметра p. По умолчанию функция вычисляет число условий с 2 нормами (p=2). Функция принимает на вход матрицу и возвращает скалярное значение, представляющее число обусловленности матрицы.
numpy.linalg.matrix_rank(): Принимает матрицу в качестве входных данных и возвращает скалярное значение, представляющее ранг матрицы. Ранг матрицы — это максимальное количество линейно независимых строк или столбцов в матрице. Он предоставляет информацию о размерности пространства столбцов или пространства строк матрицы.
numpy.linalg.cholesky(): Используется для вычисления разложения Холецкого эрмитовой положительно определенной матрицы. Разложение Холецкого факторизует матрицу в произведение нижней треугольной матрицы и ее сопряженного транспонирования. Функция принимает на вход матрицу и возвращает нижнюю треугольную матрицу. Разложение Холецкого полезно при решении линейных систем уравнений и создании случайных выборок из многомерного нормального распределения.
numpy.linalg.qr(): Используется для вычисления QR-факторизации матрицы. QR-факторизация разлагает матрицу на произведение ортогональной матрицы и верхней треугольной матрицы. Функция принимает на вход матрицу и возвращает ортогональную матрицу и верхнюю треугольную матрицу. QR-факторизация полезна при решении линейных задач наименьших квадратов, поиске собственных значений и собственных векторов, а также решении линейных систем уравнений.
numpy.linalg.svd(): Используется для вычисления разложения по сингулярным значениям (SVD) матрицы. SVD факторизует матрицу в произведение трех матриц: левой сингулярной матрицы, диагональной сингулярной матрицы и правой сингулярной матрицы. Функция принимает на вход матрицу и возвращает три матрицы. SVD полезен при решении линейных систем уравнений, поиске псевдообратной матрицы и выполнении анализа главных компонент (PCA).
Решение линейных систем
Мы можем использовать обратные матрицы для решения систем линейных уравнений.
Например, рассмотрим:
3x + 5y = 14
2x + 4y = 8
Это можно представить в матричной форме как:
A * x = b
Где:
A = np.array([[3, 5], [2, 4]])
x = np.array([x, y])
b = np.array([14, 8])
Теперь нужно найти х:
A = np.array([[3, 5], [2, 4]])
b = np.array([14, 8])
x = np.linalg.inv(A).dot(b) print(x)
Здесь мы вычислили обратную матрицу коэффициентов A и умножили ее на b, чтобы найти x.
Решение линейных систем с регулярной матрицей
Предположим, что у нас есть система линейных алгебраических уравнений, заданная формулой: Ax=b, где A∈Cn×n и b∈Cn.
Чтобы найти решение для x, мы можем использовать метод numpy.linalg.solve. Как мы наверняка знаем из занятий алгебры, точное решение существует тогда и только тогда, когда A представляет собой квадратную матрицу полного ранга (также называемую регулярной матрицей), которая также требуется для упомянутого метода решения.
Пример. Решите систему линейных уравнений, заданную формулой:
A = matrix('0 2 1; 3 -1 2; 1 -1 1')
b = matrix('2; -3; -3')
print(A)
# Убедитесь, что A является полноранговым
print('rank(A) =', linalg.matrix_rank(A))
matrix (rank(A) = n)
x = linalg.solve(A,b) # Получить решение для линейной системы A*x=b
print(x)
print(b-A*x) # Убедитесь, что x является допустимым решением для A*x=b
Собственные значения и собственные векторы
Предположим, у нас есть линейное преобразование, заданное квадратной матрицей A.
Тогда матрица A имеет скалярное собственное значение λ, связанное с ненулевым собственным вектором v, если:
Очевидная геометрическая интерпретация такова, что собственные векторы - это векторы, на которые не влияет данное преобразование с точки зрения вращения, а только растягивание (масштабирование) на соответствующий коэффициент λ. Еще один фундаментальный смысл связан с системами линейных дифференциальных уравнений, но это уже другая история.
Если матрица А имеет n линейно независимых собственных векторов, то матрицу A можно факторизовать как:
где D — диагональная матрица, содержащая собственные значения на диагонали. Столбцы U являются собственными векторами, что делает U унитарной матрицей. Обратите внимание, что таким способом можно факторизовать только диагонализируемые матрицы, что также называется спектральным разложением или собственным разложением.
Пример. Найдите собственные значения и собственные векторы матрицы.
В NumPy есть функция numpy.linalg.eig(a): эта функция используется для вычисления собственных значений и правых собственных векторов квадратного массива. Используем её в данном примере.
A = matrix('2 0 0; 0 3 4; 0 4 9')
λ, U = linalg.eig(A)
print('Eigenvalues:')
print(λ)
print('\nEigenvectors:')
print(U)
Убедитесь, что Av−λv=0 для каждой пары собственное значение-собственный вектор:
eps = 0
for i in range(A.shape[0]):
eps += sum(A*U[:,i] - λ[i]*U[:,i])
print(eps)
#1.6653345369377348e-16
Подтвердите это:
print( U * diag(λ) * U.T )
Собственные значения можно найти также с помощью np.linalg.eigvals():
mat = np.array([[1, 2], [3, 4]])
print(np.linalg.eigvals(mat))
Снова найдём собственные значения и собственные векторы матрицы с помощью np.linalg.eig():
evals, evecs = np.linalg.eig(mat)
print(evals) # Собственные значения
print(evecs) # Матрица с собственными векторами в виде столбцов
Столбцы матрицы собственных векторов соответствуют собственным векторам отсортированных собственных значений.
Собственные значения и собственные векторы полезны для разложения матрицы, анализа главных компонентов и многого другого.
В NumPy есть также функция numpy.linalg.eigh(a, UPLO='L').
Эта функция используется для возврата собственных значений и собственных векторов комплексной эрмитовой (сопряженной симметричной) или действительной симметричной матрицы. Возвращает два объекта, одномерный массив, содержащий собственные значения a и двумерный квадратный массив или матрица (в зависимости от типа входных данных) соответствующих собственных векторов (в столбцах).
Результат:
Эта операция имеет множество практических приложений в различных областях. В алгоритмах машинного обучения, таких как метод главных компонент (PCA) и факторный анализ, собственные значения и собственные векторы матрицы ковариации используются для нахождения главных компонент или скрытых факторов данных.
Разложение матрицы
NumPy предоставляет различные методы декомпозиции матрицы, которые полезны для машинного обучения и анализа данных.
LU-разложение
LU-разложение разбивает матрицу на нижнюю и верхнюю треугольные матрицы. Мы можем выполнить LU-декомпозицию в NumPy, с помощью linalg.lu() из библиотеки SciPy. Данная библиотека является расширением NumPy и использует её массивы и матрицы для работы с данными.
from scipy.linalg import lu
A = np.array([[1, 2], [3, 4]])
P, L, U = lu(A)
print(P)
print(L)
print(U)
Это также возвращает матрицу перестановок P для применения замены строк.
Разложение Холецкого
Разложение Холецкого используется для положительно определенных матриц. Метод np.linalg.cholesky() возвращает нижний треугольный коэффициент Холецкого:
A = np.array([[4, -2, 1], [-2, 6, -3], [1, -3, 7]])
L = np.linalg.cholesky(A)
print(L)
Где L*L.T = А.
Разложение по сингулярным значениям (SVD)
SVD факторизует матрицу на сингулярные векторы и сингулярные значения.
В NumPy:
A = np.array([[1, 0, 0, 0, 2], [0, 0, 3, 0, 0], [0, 0, 0, 0, 0], [0, 2, 0, 0, 0]])
U, s, Vh = np.linalg.svd(A) print(s)
Решение линейных систем методом наименьших квадратов
Рассмотрим систему линейных уравнений Ax=b, где A ε Cm×n, b ε Cm и m>n. Такая система называется переопределенной, потому что в ней больше уравнений, чем неизвестных, поэтому, как правило, она не имеет единственного решения. Однако мы можем искать решение методом наименьших квадратов, которое минимизирует евклидову норму остатков, то есть:
Решение методом наименьших квадратов для систем с матрицей полного ранга
В случае, когда А является матрицей полного ранга, то есть RankA=min(m,n), решение методом наименьших квадратов можно легко найти как решение линейной системы:
В таком случае
Решение методом наименьших квадратов для систем с матрицей любого ранга
Обсудим наиболее общий случай, когда матрица A
может быть прямоугольным, а его ранг может быть произвольным, то есть ранг A>0. В таком случае обратная эрмитова матрица не существует, поэтому нам нужно разработать другой метод решения.
Использование разложения по сингулярным значениям:
мы можем найти это:
Если мы обозначим:
мы видим, что
минимизируется тогда и только тогда, когда
Благодаря этому мы можем найти решение для x по методу наименьших квадратов как:
можно найти как:
Пример.
Найдите решение методом наименьших квадратов для переопределенной системы линейных уравнений, заданной формулой:
A = matrix('1 0 -1 2; 1 1 1 -1; 0 -1 -2 3; 5 2 -1 4; -1 2 5 -8')
b = matrix('-1; 2; -3; 1; 7')
n = A.shape[1]
r = linalg.matrix_rank(A)
print(A)
print('n =', n)
print('rank(A) =', r)
Мы видим, что ранг матрицы A равен 2, что меньше количества неизвестных, и это делает матрицу недостаточной по рангу.
# SVD-разложение матрицы A
U, σ, VH = linalg.svd(A, full_matrices=False)
V = VH.H
print('Singular values:\n', σ)
print('\nLeft-singular vectors:')
print(U)
print('\nRight-singular vectors:')
print(V)
# Псевдообратная обратная связь Мура – Пенроуза
sigma_inv = diag(hstack([1/σ[:r], zeros(n-r)]))
A_plus = V * sigma_inv * U.H
print('\nMoore–Penrose pseudoinverse of A:')
print(A_plus)
# Решение методом наименьших квадратов
x = A_plus * b
print('\nLeast-squares solution x:')
print(x)
# Ошибка решения ||b-A*x||
eps = linalg.norm(b-A*x)
print('\nError of the least-squares solution: ||b-A*x|| =', eps)
Итак, используя разложение SVD и псевдообратную матрицу Мура-Пенроуза, мы нашли решение методом наименьших квадратов:
Мы также можем использовать некоторые встроенные методы из пакета NumPy. Для вычисления псевдообратной матрицы Мура-Пенроуза существует метод numpy.linalg.pinv:
print('Moore–Penrose pseudoinverse by NumPy:')
print(linalg.pinv(A))
В NumPy есть даже метод numpy.linalg.lstsq для вычисления решения по методу наименьших квадратов для любой линейной системы, которая может быть недостаточно, хорошо или переопределена:
print('Least-squares solution by NumPy:')
print(linalg.lstsq(A,b,rcond=-1)[0])
Практические примеры использования
Линейная регрессия
Линейная регрессия - это метод в машинном обучении, который используется для:
- Прогнозирования значений зависимой переменной на основе значений независимых переменных. Это может быть полезно, например, для прогнозирования цен на недвижимость, спроса на товары и т.д.
- Оценки влияния каждой независимой переменной на зависимую переменную. Это помогает понять, какие переменные оказывают наибольшее влияние на результат.
- Оценки модели и её пригодности для прогнозирования. Одним из инструментов для оценки модели является коэффициент детерминации (R2), который показывает, насколько хорошо модель соответствует данным.
Рассмотрим пример программы, которая выполняет линейную регрессию для заданных данных и строит график линии регрессии.
Результат:
Линия регрессии — это прямая линия, которая наилучшим образом соответствует зависимости между независимой переменной (в данном случае x) и зависимой переменной (y). Она используется для предсказания значений зависимой переменной на основе независимой переменной. В данном случае, линия регрессии строится с использованием метода наименьших квадратов (МНК), который минимизирует сумму квадратов отклонений между фактическими и предсказанными значениями.
Пошаговое объяснение программы:
0. Импортируются необходимые библиотеки: numpy и matplotlib.pyplot.
- Создается массив x с помощью функции np.arange(0,9), который представляет собой последовательность чисел от 0 до 8.
- Создается массив A, который является двумерным массивом, содержащим x и массив из единиц. Это нужно для выполнения линейной регрессии.
- Задаются значения массива y, которые представляют зависимую переменную.
- С помощью функции np.linalg.lstsq(A.T, y)[0] выполняется линейная регрессия. Результат сохраняется в массив w.
- Создается массив line, который представляет собой линию регрессии, вычисленную на основе массива x и коэффициентов w.
- С помощью функции plt.plot(x, line, 'r-') строится график линии регрессии.
- С помощью функции plt.plot(x, y, 'o') строится график точек данных.
- С помощью функции plt.show() отображается график.
Интерполяция полиномов с помощью линейных систем
Вы можете использовать линейные системы для расчета коэффициентов полинома, чтобы эти полиномы включали некоторые определенные точки.
Например, рассмотрим полином второй степени y = P(x) = a₀ + a₁x + a₂x². Напомним, что при построении полинома второй степени вы получаете параболу, которая будет разной в зависимости от коэффициентов a₀, a₁ и a₂.
Теперь предположим, что вы хотите найти конкретный полином второй степени, включающий точки (x, y) (1, 5), (2, 13) и (3, 25). Как можно вычислить a₀, a₁ и a₂ так, чтобы P(x) включала эти точки в свою параболу? Другими словами, вы хотите найти коэффициенты многочлена на этом рисунке:
Для каждой точки, которую вы хотите включить в параболу, вы можете использовать общее выражение полинома, чтобы получить линейное уравнение. Например, взяв вторую точку (x=2, y=13) и учитывая, что y = a₀ + a₁x + a₂x², вы можете написать следующее уравнение:
import numpy as np
from scipy import linalg
A = np.array([[1, 1, 1], [1, 2, 4], [1, 3, 9]])
print(linalg.det(A))
#1.9999999999999996
Стоит отметить, что существование решения зависит только от A. Поскольку значение определителя не равно нулю, вы можете быть уверены, что для системы существует единственное решение. Вы можете решить эту задачу, используя метод обратной матрицы с помощью следующего кода:
b = np.array([5, 13, 25]).reshape((3, 1))
a = linalg.inv(A) @ b
print(a)
#array([[1.], [2.], [2.]])
Этот результат говорит вам, что a₀ = 1, a₁ = 2 и a₂ = 2 являются решением для системы. Другими словами, полином, включающий точки (1, 5), (2, 13) и (3, 25), имеет вид y = P(x) = 1 + 2x + 2x². Вы можете проверить решение для каждой точки, введя x и проверив, что P(x) равно y.
В качестве примера системы без какого-либо решения предположим, что вы пытаетесь интерполировать параболу с точками (x, y), заданными (1, 5), (2, 13) и (2, 25). Если вы внимательно посмотрите на эти числа, то заметите, что вторая и третья точки учитывают x = 2 и разные значения y, что делает невозможным найти функцию, включающую обе точки. Выполнив те же действия, что и раньше, вы получите следующие уравнения этой системы:
A = np.array([[1, 1, 1], [1, 2, 4], [1, 2, 4]])
print(linalg.det(A))
# 0.0
Вы можете заметить, что значение определителя равно нулю, а это означает, что система не имеет единственного решения. Это также означает, что обратной матрицы коэффициентов не существует. Другими словами, матрица коэффициентов сингулярна.
В зависимости от архитектуры вашего компьютера вместо нуля вы можете получить очень маленькое число. Это происходит из-за числовых алгоритмов, которые det() использует для вычисления определителя. В этих алгоритмах ошибки числовой точности делают этот результат не совсем равным нулю.
В общем, всякий раз, когда вы встречаете небольшое число, вы можете сделать вывод, что система не имеет единственного решения.
Вы можете попытаться решить линейную систему, используя метод обратной матрицы, используя следующий код:
b = np.array([5, 13, 25]).reshape((3, 1))
x = linalg.inv(A) @ b
print(x)
#LinAlgError: Traceback (most recent call last) <ipython-input-10-e6ee9b06a6fe> #in <module> ----> 1 x = linalg.inv(A) @ b LinAlgError: singular matrix
Поскольку система не имеет решения, вы получаете исключение, сообщающее, что матрица коэффициентов сингулярна.
Если в системе имеется более одного решения, вы получите аналогичный результат. Значение определителя матрицы коэффициентов будет равно нулю или очень мало, что указывает на то, что матрица коэффициентов снова является сингулярной.
В качестве примера системы с более чем одним решением вы можете попытаться интерполировать параболу, учитывая точки (x, y), заданные (1, 5), (2, 13) и (2, 13). Как вы могли заметить, здесь вы рассматриваете две точки в одной и той же позиции, что допускает бесконечное количество решений для a₀, a₁ и a₂.
Минимизация ошибки с помощью наименьших квадратов
Вы видели, что иногда не удается найти полином, который точно соответствовал бы набору точек. Однако обычно, когда вы пытаетесь интерполировать полином, вас не интересует точная подгонка. Вы просто ищете решение, которое аппроксимирует точки, обеспечивая минимальную возможную ошибку.
Обычно это тот случай, когда вы работаете с реальными данными. Обычно он включает в себя некоторый шум, вызванный ошибками, возникающими в процессе сбора, например неточностью или неисправностью датчиков, а также опечатками при вводе данных вручную.
Используя метод наименьших квадратов, можно найти решение для интерполяции полинома, даже если матрица коэффициентов сингулярна. Используя этот метод, вы будете искать коэффициенты полинома, который обеспечивает минимальную квадратичную ошибку при сравнении полиномиальной кривой с вашими точками данных.
На самом деле метод наименьших квадратов обычно используется для сопоставления полиномов с большими наборами точек данных. Идея состоит в том, чтобы попытаться разработать модель, отражающую наблюдаемое поведение.
Примечание. Если линейная система имеет единственное решение, то решение методом наименьших квадратов будет равно этому уникальному решению.
Например, вы можете разработать модель, позволяющую прогнозировать цены на автомобили. Для этого вы можете собрать некоторые реальные данные, включая цену автомобиля и некоторые другие характеристики, такие как пробег, год и тип автомобиля. Используя эти данные, вы можете разработать полином, моделирующий цену как функцию других характеристик, и использовать метод наименьших квадратов для поиска оптимальных коэффициентов этой модели.
Для решения задач наименьших квадратов scipy.linalg предоставляет функцию lstsq(). Чтобы увидеть, как это работает, рассмотрим предыдущий пример, в котором вы пытались подогнать параболу к точкам (x, y), заданным (1, 5), (2, 13) и (2, 25). Помните, что эта система не имеет решения, так как существуют две точки с одинаковым значением x.
Как и раньше, используя модель y = a₀ + a₁x + a₂x², вы получаете следующую линейную систему:
import numpy as np
from scipy import linalg
A = np.array([[1, 1, 1], [1, 2, 4], [1, 2, 4]])
b = np.array([5, 13, 25]).reshape((3, 1))
p, *_ = linalg.lstsq(A, b)
print(p)
# array([[-0.42857143], [ 1.14285714], [ 4.28571429]])
В этой программе вы настроили следующее:
Строки с 1 по 2: вы импортируете numpy как np и linalg из scipy, чтобы использовать linalg.lstsq().
Строки с 4 по 5: вы создаете матрицу коэффициентов A, используя массив NumPy с именем A, и вектор с независимыми членами b, используя массив NumPy с именем b.
Строка 7: Вы вычисляете решение задачи методом наименьших квадратов, используя linalg.lstsq(), который принимает в качестве входных данных матрицу коэффициентов и вектор с независимыми членами.
lstsq() предоставляет несколько фрагментов информации о системе, включая остатки, ранг и сингулярные значения матрицы коэффициентов. В этом случае вас интересуют только коэффициенты полинома для решения задачи по критериям наименьших квадратов, которые хранятся в p.
Как видите, даже если рассматривать линейную систему, не имеющую точного решения, lstsq() предоставляет коэффициенты, минимизирующие квадраты ошибок. С помощью следующего кода вы можете визуализировать полученное решение, построив параболу и точки данных:
import matplotlib.pyplot as plt
x = np.linspace(0, 3, 1000)
y = p[0] + p[1] * x + p[2] * x ** 2
plt.plot(x, y)
plt.plot(1, 5, "ro")
plt.plot(2, 13, "ro")
plt.plot(2, 25, "ro")
Эта программа использует matplotlib для построения результатов:
Строка 1: вы импортируете matplotlib.pyplot как plt, что типично.
Строки с 3 по 4: вы создаете массив NumPy с именем x со значениями от 0 до 3, содержащий 1000 точек. Вы также создаете массив NumPy с именем y с соответствующими значениями модели.
Строка 6: вы строите кривую параболы, полученной с помощью модели, заданной точками в массивах x и y.
Строки с 7 по 9: красным цветом («ro») вы рисуете три точки, использованные для построения модели.
Обратите внимание, как кривая, предоставленная моделью, пытается максимально точно аппроксимировать точки.
Помимо lstsq(), существуют и другие способы расчета решений методом наименьших квадратов с использованием SciPy. Одна из альтернатив — использование псевдоинверсии, которую вы изучите далее.
Получение решений методом наименьших квадратов с использованием псевдообратной задачи
Другой способ вычисления решения по методу наименьших квадратов — использование псевдообратной матрицы Мура-Пенроуза.
Вы можете думать о псевдообратной матрице как об обобщении обратной матрицы, поскольку она равна обычной обратной матрице, когда матрица не является сингулярной.
Однако, когда матрица сингулярна, что имеет место в линейных системах, у которых нет единственного решения, тогда псевдообратная вычисляет матрицу, которая обеспечивает наилучшее соответствие, что приводит к решению методом наименьших квадратов.
Используя псевдообратный метод, можно найти коэффициенты параболы, использованной в предыдущем примере:
import numpy as np
from scipy import linalg
A = np.array([[1, 1, 1], [1, 2, 4], [1, 2, 4]])
b = np.array([5, 13, 25]).reshape((3, 1))
A_pinv = linalg.pinv(A)
p2 = A_pinv @ b
print(p2)
# array([[-0.42857143], [ 1.14285714], [ 4.28571429]])
Этот код очень похож на код из предыдущего примера, за исключением выделенных строк:
Строка 7: Вы вычисляете псевдообратную матрицу коэффициентов и сохраняете ее в A_pinv.
Строка 9: Следуя тому же подходу, который используется для решения линейных систем с обратной матрицей, вы вычисляете коэффициенты уравнения параболы, используя псевдообратное, и сохраняете их в векторе p2.
Как и следовало ожидать, решение методом наименьших квадратов такое же, как и решение lstsq(). В этом случае, поскольку A является квадратной матрицей, pinv() предоставит квадратную матрицу тех же размеров, что и A, оптимизируя ее для наилучшего соответствия в смысле метода наименьших квадратов:
print(A_pinv)
# array([[ 1. , -0.14285714, -0.14285714], [ 0.5 , -0.03571429, -0.03571429], # [-0.5 , 0.17857143, 0.17857143]])
Однако стоит отметить, что вы также можете вычислять pinv() для неквадратных матриц, что обычно и происходит на практике.
Прогнозирование цен на автомобили с помощью метода наименьших квадратов
В этом примере вы собираетесь построить модель с использованием метода наименьших квадратов для прогнозирования цен на подержанные автомобили, используя данные из набора данных подержанных автомобилей. Этот набор данных представляет собой огромную коллекцию с 957 МБ списков транспортных средств с сайта craigslist.org, включая самые разные типы транспортных средств.
К сожалению у нас этот набор данных недоступен для скачивания. Поэтому поручим сгенерировать их Алисе из яндекса.
Данные копируем в блокнот и сохраняем в файл cars.csv. Формат csv удобен тем, что импорт и экспорт из него поддерживают многие программы и библиотеки для работы с данными. Мы в, частности, будем использовать библиотеку pandas, которая загрузит данные из этого файла и сформирует из них датафрейм. Из этого датафрейма затем можно легко получить массивы, пригодные для работы с NumPy.
Сделаем пример с машинами в Google Colab. Для начала работы в этой среде нужно иметь аккаунт Google. Также для выполнения примера понадобится Google Drive.
Если это есть, вводим ссылку: https://colab.research.google.com и попадаем на стартовую страницу:
Жмём «Создать блокнот» и получаем чистый рабочий лист проекта.
Слева расположена панель быстрого доступа к вкладкам. Первая снизу вкладка - «Файлы», жмём её, чтобы подключить свой Google диск к проекту:
Жмём на кнопку «Подключить диск» (третья слева).
В отдельной вкладке браузера заходим в Google Drive и загружаем файл cars.csv (нажав кнопку «Создать», а затем «Загрузить файлы»).
Возвращаемся на вкладку с Google Colab, жмём кнопку «Обновить» (вторая слева) на вкладке Файлы. И двойным щелчком открываем файл cars.csv для предварительного просмотра.
Если всё как надо, переходим к написанию программы. В Colab программа состоит из блоков, которые можно добавлять наводя курсор на первую свободную строку. На чистом листе нужно навести курсор немного ниже панели с командами «+Код» и «+Текст» (или нажать на «+Код»).
Для начала импортируем библиотеку pandas, которая будет использоваться для предварительной обработки набора данных.
import pandas as pd
cars_data = pd.read_csv("/content/drive/MyDrive/cars.csv")
print(cars_data.columns)
Запустим этот блок кода, нажав на кнопку с треугольником. Под блоком кода появится результат его выполнения.
Здесь мы загрузили файл cars.csv и создали из него датафрейм cars_data. Этот объект включает в себя атрибут с именем columns, который позволяет просматривать имена столбцов, включённых в данные. Также мы можем просматривать записи, используя атрибут iloc:
print(cars_data.iloc[0])
Чтобы использовать эти данные для построения модели наименьших квадратов, вам необходимо представить категориальные данные в числовом виде. В большинстве случаев категориальные данные преобразуются в набор фиктивных переменных, то есть переменных, которые могут принимать значения 0 или 1.
В качестве примера такого преобразования рассмотрим столбец топлива, который может принимать значение газ или дизельное топливо. Вы можете преобразовать этот категориальный столбец в фиктивный столбец с именем Fuel_gas, который принимает значение 1, если топливом является газ, и 0, если топливом является дизельное топливо.
Обратите внимание: вам понадобится только один фиктивный столбец для представления категориального столбца, который может принимать два разных значения. Аналогично, для категориального столбца, который может принимать N значений, вам понадобится N-1 фиктивных столбцов, поскольку одно из значений будет считаться значением по умолчанию.
В pandas вы можете преобразовать эти категориальные столбцы в фиктивные столбцы с помощью get_dummies():
Здесь вы создаете новый DataFrame с именем cars_data_dummies, который включает фиктивные переменные для столбцов, указанных в аргументе columns. Теперь вы можете проверить новые столбцы, включенные в этот DataFrame:
Теперь, когда вы преобразовали категориальные переменные в наборы фиктивных переменных, вы можете использовать эту информацию для построения своей модели. По сути, модель будет включать коэффициент для каждого из этих столбцов, кроме цены, которая будет использоваться в качестве выходных данных модели. Цена будет определяться взвешенной комбинацией других переменных, где веса задаются коэффициентами модели.
Однако принято учитывать дополнительный коэффициент, который представляет собой постоянное значение, добавляемое к взвешенной комбинации других переменных. Этот коэффициент известен как перехват, и вы можете включить его в свою модель, добавив к данным дополнительный столбец со всеми строками, равными 1:
cars_data_dummies["intercept"] = 1
Теперь, когда у вас есть все данные, вы можете сгенерировать массивы NumPy для построения модели с помощью scipy.linalg. Это то, что вы сделаете дальше.
Построение модели
Чтобы сгенерировать массивы NumPy для ввода в lstsq() или pinv(), вы можете использовать .to_numpy():
A = cars_data_dummies.drop(columns=["price"]).to_numpy()
b = cars_data_dummies.loc[:, "price"].to_numpy()
Матрица коэффициентов A представлена всеми столбцами, кроме цены. Вектор b с независимыми членами задается значениями, которые вы хотите спрогнозировать, в данном случае это столбец цен. Если установлены A и b, вы можете использовать lstsq(), чтобы найти решение по методу наименьших квадратов для коэффициентов:
fromscipy importlinalg
p, *_ = linalg.lstsq(A, b)
print(p)
Это коэффициенты, которые вам следует использовать для моделирования цены с точки зрения взвешенной комбинации других переменных, чтобы минимизировать квадратичную ошибку. Как вы видели, эти коэффициенты также можно получить, используя pinv() со следующим кодом:
p2 = linalg.pinv(A) @ b
print(p2)
Теперь, когда вы получили модель, вы будете использовать ее для прогнозирования цены автомобиля.
Прогнозирование цен
Используя модель, полученную методом наименьших квадратов, вы можете спрогнозировать цену автомобиля, представленного вектором со значениями каждой из переменных, используемых в модели:
print(cars_data_dummies.drop(columns=["price"]).columns)
Итак, 4-цилиндровый хэтчбек 2010 года выпуска, с АКПП, газовое топливо, пробег 50 000 миль, в хорошем состоянии, можно представить следующим вектором:
import numpy as np
car = np.array(
[2010, 500000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0]
Вы можете получить прогноз цены, рассчитав скалярное произведение вектора автомобиля и вектора p коэффициентов. Поскольку оба вектора представляют собой одномерные массивы NumPy, вы можете использовать @ для получения скалярного произведения:
predicted_price = p @ car
print (predicted_price)
В этом примере прогнозируемая цена хэтчбека составляет примерно 1698 долларов. Стоит отметить, что коэффициенты модели включают некоторую неопределенность, поскольку данные, используемые для получения модели, могут быть смещены, например, в сторону определенного типа автомобиля.
Кроме того, выбор модели играет большую роль в качестве оценок. Метод наименьших квадратов — один из наиболее часто используемых методов построения моделей, поскольку он прост и позволяет получить объяснимые модели. В этом примере вы увидели, как использовать scipy.linalg для построения таких моделей.
На этом всё. В следующий раз я намерен перейти к рассматрению математических библиотек Python в приложении к лингвистическим задачам.
Использованные источники:
1. Дж. Вандер Плас, "Python для сложных задач. Наука о данных и машинное обучение". С: Питер, 2018.
2. An In-Depth Guide to Linear Algebra in NumPy
https://llego.dev/posts/numpy-linear-algebra-guide/
3. Linear Algebra in Python: Matrix Inverses and Least Squares
https://realpython.com/python-linear-algebra/
4. scipy.linalg.lu
https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu.html
5. Linear algebra (numpy.linalg)
https://numpy.org/devdocs/reference/routines.linalg.html
6. Numpy | Linear Algebra
https://www.geeksforgeeks.org/numpy-linear-algebra/
7. Датасет для прогнозирования цен автомобилей:
https://drive.google.com/file/d/1paHFSwtzEPqb9B2GSAS0SG_xeyb4uwFN/view?usp=sharing