Привет, Хабр! Я недавно начал свой путь в 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)


  1. haqreu
    19.08.2023 16:50
    +4

    Ваша реализация не позволит работать с большими (например, миллион на миллион) матрицами, которых в датасайнс очень много. Да и с маленькими будет работать медленно и численно неустойчиво. Пользуйтесь сайпаем, хотя полезно иметь представление, что у него под капотом, да. Но первые шаги в это проще на бумаге делать...


  1. sci_nov
    19.08.2023 16:50
    +2

    Транспонировать можно и так:

    # Возвращает транспонированную матрицу, т.е. Y = X^T
    def transpose(X):
        return list(map(list, zip(*X)))


    1. economist75
      19.08.2023 16:50
      +3

      И так:

      X.T


      1. sci_nov
        19.08.2023 16:50

        В numpy?


        1. economist75
          19.08.2023 16:50

          Да, конечно в нем. Очень сложно представить на практике работу с матрицами вне numpy/pyarrow. Академический, учебный кейс - тот же: вращать матрицы нужно в том, в чем они вращаются на практике.


          1. sci_nov
            19.08.2023 16:50

            Я делал проект учебный, для студентов, там норм использовать вложенные списки. Плюс там целочисленная арифметика. Я поэтому сразу запостил свой вариант транспонирования :)


  1. vadimr
    19.08.2023 16:50
    +5

    Серьёзная работа с матрицами для data science - это сложная вещь, тут на таком уровне делать нечего.

    В Питоне есть специальные пакеты для этого, использующие эффективное представление матриц и качественные численные алгоритмы.

    Что касается операции транспонирования матрицы, то она вообще не должна иметь физической реализации в реальной работе. Это просто переобозначение индексов.


  1. sci_nov
    19.08.2023 16:50
    +5

    А разве определитель обратной матрицы не равен обратному определителю прямой?


    1. akhalat
      19.08.2023 16:50
      +3

      Там и у исходной неправильно вычислен


      1. sci_nov
        19.08.2023 16:50
        +2

        )))


  1. Daddy_Cool
    19.08.2023 16:50
    +5

    Я python не знаю, но осуждаю. Писать в чистом виде математику на python это кажется самое вредное, что только можно. Если всё то же самое сделать с использованием библиотек - было бы уже как-то полезно.
    P.S.
    Можно я напишу статью как выводить "Hello, word!" на Си?


    1. mikko_kukkanen
      19.08.2023 16:50
      -2

      Мне было бы интересно какая будет разница если выводить "Hello, world!" на C# и С++. Когда я писал на С, этих "диалектов" еще не было. А разобраться, что под капотом в библиотеках - прекрасное начало долгого пути.


    1. Tzimie
      19.08.2023 16:50

      Julia !!!


  1. 9982th
    19.08.2023 16:50

    del


  1. rdmr
    19.08.2023 16:50

    Действительно зачем зная про нампай, математикий аппарат тут довольно простой.


  1. namedictor
    19.08.2023 16:50

    Эх, алгебра 1 курс, вспоминаю как через метод гаусса считали этот кошмар вещественночисленный


  1. net_racoon
    19.08.2023 16:50

    Простите, совсем не дата саинтист и даже не разраб. Но разве такие вещи не должны на GPU считать как-то? Или там совсем неподходящая структура данных?