Привет, Хабр! Я недавно начал свой путь в data science и хочу поделиться свои опытом в написание алгоритмов для работы с матрицами, я планирую активно пополнять свой репозиторий новыми функциями (понимаю что можно все сделать с нампаем).
Транспонирование матрицы
Во-первых, что такое транспонирование, это операция в последствие которой строки и столбцы меняются местами.
На фото наглядный пример поворота, это довольна простая операция, которую я реализовал в два цикла.
def transpose_matrix(matrix: list[list]) -> list[list]:
transposed_matrix = [[0 for i in range(len(matrix))] for i in range(len(matrix[0]))]
for i in range(len(matrix)):
for j in range(len(matrix[0])):
transposed_matrix[j][i] = matrix[i][j]
return transposed_matrix
Матрицы представляют вложенными списками, поэтому функция принимает список в списке. Переменная transposed_matrix хранит в себе нулевую матрицу с размерами идентичной нашей только с поворотом в 90 градусов, затем происходит двойная итерация в которой элементы меняются местами, в результате выходит поворот.
Ранг матрицы
Ранг матрицы — наивысший из порядков всевозможных ненулевых миноров этой матрицы.
def rank_of_matrix(matrix: list[list]) -> int:
rank = min(len(matrix), len(matrix[0]))
row_index = 0
for i in range(len(matrix[0])):
found_nonzero = False
for j in range(row_index, len(matrix)):
if matrix[j][i] != 0:
found_nonzero = True
matrix[row_index], matrix[j] = matrix[j], matrix[row_index]
break
if found_nonzero:
for j in range(row_index + 1, len(matrix)):
factor = matrix[j][i] / matrix[row_index][i]
for k in range(i, len(matrix[0])):
matrix[j][k] -= matrix[row_index][k] * factor
row_index += 1
return rank
В самом начале в переменная rank инициализируется минимальным значением между количеством строк и количеством столбцов в матрице. Это определяется тем, что ранг матрицы не может быть больше, чем количество строк или столбцов, row_index равна нулю, так как это будет индексом строки, с котором мы будем работать на каждом шаге. Так же проходимся двумя циклами и ищем первый ненулевой элемент в текущем столбце. Если такой элемент найден, он меняется местами. Это делается для того, чтобы разместить ненулевой элемент на позиции (
row_index, i)
. Если найден ненулевой элемент в текущем столбце на позиции (
row_index, i)
, производится процесс приведения матрицы к ступенчатому виду. Для этого вычисляется коэффициент factor, равный matrix[j][i] / matrix[row_index][i], и затем вычитается matrix[row_index][k] * factor из всех элементов строки j, начиная с позиции i. Это приводит к обнулению всех элементов ниже matrix[row_index][i] в столбце i. Увеличиваем row_index на 1, чтобы перейти к следующей строке и продолжить процесс приведения к ступенчатому виду. По окончании алгоритма возвращается rank, которая представляет собой максимальное количество линейно независимых строк или столбцов в матрице.
Инверсия матрицы
Обратная матрица - это матрица, которая умножается на исходную матрицу таким образом, что их произведение дает единичную матрицу. Другими словами, если у нас есть матрица A и ее обратная матрица обозначается как A^-1.
def inverse_matrix(matrix: list[list]) -> list[list]:
augmented_matrix = [
[
matrix[i][j] if j < len(matrix) else int(i == j - len(matrix))
for j in range(2 * len(matrix))
]
for i in range(len(matrix))
]
for i in range(len(matrix)):
pivot = augmented_matrix[i][i]
if pivot == 0:
raise ValueError("Matrix is not invertible")
for j in range(2 * len(matrix)):
augmented_matrix[i][j] /= pivot
for j in range(len(matrix)):
if i != j:
scalar = augmented_matrix[j][i]
for k in range(2 * len(matrix)):
augmented_matrix[j][k] -= scalar * augmented_matrix[i][k]
inverse = [
[augmented_matrix[i][j] for j in range(len(matrix), 2 * len(matrix))]
for i in range(len(matrix))
]
return inverse
Переменная augmented_matrix представляет собой расширенную матрицу, добавив единичную матрицу справа. Происходит итерация, получаем значение pivot, которое является текущим диагональным элементом. Если pivot равна 0, вызывается исключение, так как матрица не обратима. Далее делаем нормирование строки, делим все элементы строки на pivot. Это делается для того, чтобы текущий диагональный элемент стал равным 1. Внутренний цикл итерируется по всем столбцам расширенной матрицы. Вычитаем из текущей строки другие строки, умноженные на значение элемента scalar. После завершения внутреннего цикла, матрица будет приведена к верхнетреугольному виду. Создается матрица inverse, которая содержит элементы справа от вертикальной черты в augmented_matrix. Возвращается матрица inverse, которая представляет собой обратную матрицу исходной матрицы.
Заключение
На этом все, это моя первая статья в целом, не знаю что с этого выйдет, но очень хотелось написать и поделиться реализацией. Я привел три алгоритма по обработке матриц, в моем репозитории вы можете посмотреть другие. Пока что я только пробую себя в науке о данных, в дальнейшем хочу производить различный анализы и делиться ими на сайте.
Спасибо за выделенное время! Надеюсь было интересно, пишите свои мысли, критикуйте, буду рад исправиться и писать лучшие статьи.
Комментарии (17)
sci_nov
19.08.2023 16:50+2Транспонировать можно и так:
# Возвращает транспонированную матрицу, т.е. Y = X^T def transpose(X): return list(map(list, zip(*X)))
economist75
19.08.2023 16:50+3И так:
X.T
sci_nov
19.08.2023 16:50В numpy?
economist75
19.08.2023 16:50Да, конечно в нем. Очень сложно представить на практике работу с матрицами вне numpy/pyarrow. Академический, учебный кейс - тот же: вращать матрицы нужно в том, в чем они вращаются на практике.
sci_nov
19.08.2023 16:50Я делал проект учебный, для студентов, там норм использовать вложенные списки. Плюс там целочисленная арифметика. Я поэтому сразу запостил свой вариант транспонирования :)
vadimr
19.08.2023 16:50+5Серьёзная работа с матрицами для data science - это сложная вещь, тут на таком уровне делать нечего.
В Питоне есть специальные пакеты для этого, использующие эффективное представление матриц и качественные численные алгоритмы.
Что касается операции транспонирования матрицы, то она вообще не должна иметь физической реализации в реальной работе. Это просто переобозначение индексов.
Daddy_Cool
19.08.2023 16:50+5Я python не знаю, но осуждаю. Писать в чистом виде математику на python это кажется самое вредное, что только можно. Если всё то же самое сделать с использованием библиотек - было бы уже как-то полезно.
P.S.
Можно я напишу статью как выводить "Hello, word!" на Си?mikko_kukkanen
19.08.2023 16:50-2Мне было бы интересно какая будет разница если выводить "Hello, world!" на C# и С++. Когда я писал на С, этих "диалектов" еще не было. А разобраться, что под капотом в библиотеках - прекрасное начало долгого пути.
namedictor
19.08.2023 16:50Эх, алгебра 1 курс, вспоминаю как через метод гаусса считали этот кошмар вещественночисленный
net_racoon
19.08.2023 16:50Простите, совсем не дата саинтист и даже не разраб. Но разве такие вещи не должны на GPU считать как-то? Или там совсем неподходящая структура данных?
haqreu
Ваша реализация не позволит работать с большими (например, миллион на миллион) матрицами, которых в датасайнс очень много. Да и с маленькими будет работать медленно и численно неустойчиво. Пользуйтесь сайпаем, хотя полезно иметь представление, что у него под капотом, да. Но первые шаги в это проще на бумаге делать...