Введение
Инверсия матрицы — это математическая операция, которая, как следует из названия, включает в себя получение «обратной матрицы» для данной матрицы. Обратная матрица — это матрица, которая при умножении на исходную матрицу дает единичную матрицу. Важно отметить, что не все матрицы имеют обратную; матрицы, имеющие обратную, называются «обратимыми матрицами» или «несингулярными матрицами».
Процесс обращения матриц требует больших вычислительных ресурсов, в первую очередь из-за характера используемых математических операций. Стандартный алгоритм обращения матриц имеет временную сложность O(n³), что означает, что время вычислений быстро увеличивается с размером матрицы (увеличивается кубически с количеством компонентов). Кроме того, инверсия матриц включает в себя выполнение многочисленных вычислений с плавающей запятой, что может привести к ошибкам округления и числовой нестабильности, если ими не управлять должным образом.
Несмотря на эти проблемы, инверсия матриц играет фундаментальную роль во многих областях науки и техники. Например, в животноводстве матричная инверсия используется при анализе оценок параметров в смешанных линейных моделях, которые составляют основу генетической селекции. Эти модели используются для прогнозирования генетических характеристик животных с целью определения лучших особей для размножения.
Помимо животноводства, инверсия матриц также имеет решающее значение во многих других областях, включая, помимо прочего, решение систем линейных уравнений, сетевой анализ, эконометрику, системную оптимизацию, машинное обучение и обработку сигналов. Таким образом, несмотря на свои вычислительные затраты, обращение матриц является важным инструментом в арсенале современной науки и техники.
Инверсия матриц, хотя и требует больших вычислительных затрат, может быть оптимизирована с помощью ряда альтернативных методов, которые значительно снижают вычислительные затраты, особенно для больших матриц. Такие методы, как LU-разложение, разложение Холецкого, QR-разложение, метод Гаусса-Жордана и метод определителей, обеспечивают более эффективные стратегии обращения матриц. Эти методы включают в себя разложение исходной матрицы на более простые формы или применение ряда элементарных операций для приведения матрицы к более удобной для манипулирования форме, что делает вычисление обратной матрицы более эффективным и менее подверженным ошибкам округления. Во многих ситуациях вместо прямого обращения матрицы более эффективно решать соответствующие системы линейных уравнений, что можно сделать с помощью таких методов, как LU-разложение или метод Гаусса-Жордана. Такие альтернативные подходы играют решающую роль в практической реализации матричных вычислений, способствуя прогрессу в широком спектре научных и технологических областей.
Перейдём к Питону...
Во-первых, давайте создадим функцию для генерации обратимых матриц и определим нашу вычислительную стоимость, принимая во внимание время и пиковое потребление памяти, необходимое для инвертирования матрицы.
Мы можем создать эти матрицы довольно просто. Давайте сгенерируем случайную матрицу (A) и умножим ее на ее транспонирование (A.T). Это создаст положительно определенную симметричную матрицу.
import numpy as np # Generate a random matrix A = np.random.rand(1000,1000) # Create a positive definite symmetric matrix Matrix = np.dot(A, A.T) print(Matrix)
Итак, у нас есть:
[[346.5979911 265.38275116 256.26402957 ... 245.80078252 267.13380138 256.84901862] [265.38275116 360.55081534 263.97122819 ... 261.22509249 277.64630251 271.39898847] [256.26402957 263.97122819 336.49803262 ... 244.84151498 262.82901393 252.33494426] ... [245.80078252 261.22509249 244.84151498 ... 324.99095213 251.39806184 252.99577792] [267.13380138 277.64630251 262.82901393 ... 251.39806184 353.25654529 260.15549937] [256.84901862 271.39898847 252.33494426 ... 252.99577792 260.15549937 335.02079887]]
И как мы будем измерять время и пиковое потребление памяти? Мы будем использовать модули time и memory_profiler и вычислим стоимость инвертирования матрицы размером 1000 строк и столбцов напрямую с помощью функции np.linalg.inv(). из numpy, чтобы протестировать нашу методологию. Вот пример на Python:
Сначала мы импортируем соответствующие библиотеки:
import numpy as np import time from memory_profiler import memory_usage
И мы создаем функцию для инвертирования матрицы на основе numpy вместе с пакетами, которые мы будем использовать для измерения вычислительной стоимости:
def invert_matrix(Matrix):
# Measure or use initial memory
mem_before = memory_usage(-1, interval=0.1, timeout=1)[0]
# Measure the start time
start_time = time.time()
try:
A_inv = np.linalg.inv(A)
except np.linalg.LinAlgError:
print("Matriz é singular e não pode ser invertida.")
return None
# Measure the end time
end_time = time.time()
# Measure ultimate memory usage
mem_after = memory_usage(-1, interval=0.1, timeout=1)[0]
# Calculate and print runtime
print("Execution time: ", end_time - start_time, "seconds")
# Calculate and print memory usage
print("Memory usage: ", mem_after - mem_before, "MiB")
return A_inv
# Example of use
A_inv = invert_matrix(A)
if A_inv is not None:
print(A_inv)
Здесь мы используем функцию np.linalg.inv для нашего теста (которая не является прямой инверсией матрицы).
Результат:
Execution time: 0.22894835472106934 segundos Memory usage: 8.4453125 MiB [[-0.22440157 -0.11734082 0.2651931 ... 0.05743279 -0.02913667 0.03216334] [ 0.32815818 0.19427521 -0.22804691 ... -0.12614379 -0.05541012 -0.02195372] [ 0.02748701 0.008787 0.09213851 ... -0.04512599 -0.00655794 -0.04441288] ... [ 0.32767353 0.11321216 -0.23788117 ... 0.01607856 -0.07020798 0.0057955 ] [ 0.12156325 0.06088216 -0.07164596 ... 0.01916564 -0.01757929 0.00743931] [ 0.76110064 0.33654407 -0.64459896 ... -0.33153027 -0.08015695 -0.00390345]]
Здесь у нас есть время 0,22 секунды и вычислительные затраты 8,44 МБ для инвертирования нашей матрицы 1000 x 1000, которая работает правильно. Теперь перейдем к тестируемым методам.
1. Метод разложения ЛУ
Разложение LU разбивает матрицу на две матрицы — нижнюю треугольную матрицу и верхнюю треугольную матрицу. Этот метод полезен при работе с большими матрицами, так как делает вычисления более эффективными.
Мы можем использовать этот метод через функцию библиотеки scipy scipy.linalg.lu_factor().
import numpy as np
import scipy.linalg
import time
from memory_profiler import memory_usage
def invert_matrix_lu(Matrix):
n = Matrix.shape[0]
# Measure initial memory usage
mem_before = memory_usage(-1, interval=0.1, timeout=1)[0]
# Measure the start time
start_time = time.time()
try:
lu, piv = scipy.linalg.lu_factor(A)
A_inv = np.zeros((n, n))
for i in range(n):
e = np.zeros((n,))
e[i] = 1.
A_inv[:, i] = scipy.linalg.lu_solve((lu, piv), e)
except np.linalg.LinAlgError:
print("Matriz é singular e não pode ser invertida.")
return None
# Measure the end time
end_time = time.time()
# Measure final memory usage
mem_after = memory_usage(-1, interval=0.1, timeout=1)[0]
# Calculate and print runtime
print("Execution time: ", end_time - start_time, "seconds")
# Calculate and print memory usage
print("Memory usage: ", mem_after - mem_before, "MiB")
return A_inv
# Example of use
A_inv = invert_matrix_lu(A)
if A_inv is not None:
print(A_inv)
В этом коде функция invert_matrix_lu() выполняет LU-разложение матрицы A с помощью функции scipy.linalg.lu_factor(). Затем он решает несколько систем уравнений вида Ax = e, где e — единичный вектор, для вычисления столбцов обратной матрицы.
2. Метод разложения Холецкого
Этот метод в основном применяется к симметричным и положительно определенным матрицам. Матрица разлагается на нижнюю треугольную матрицу и транспонирование этой нижней треугольной матрицы.
Здесь мы снова будем использовать библиотеку scipy, но функцией будет scipy.linalg.cholesky().
import numpy as np
import scipy.linalg
import time
from memory_profiler import memory_usage
def invert_matrix_cholesky(Matrix):
n = Matrix.shape[0]
# Measure initial memory usage
mem_before = memory_usage(-1, interval=0.1, timeout=1)[0]
# Measure the start time
start_time = time.time()
try:
# Doing the Cholesky Decomposition
c = scipy.linalg.cholesky(Matrix)
# Inverting the c matrix
c_inv = scipy.linalg.inv(c)
# Getting the inverse of A
A_inv = np.dot(c_inv.T, c_inv)
except np.linalg.LinAlgError:
print("Matrix is singular and cannot be inverted, or is not suitable for Cholesky decomposition.")
return None
# Measure the end time
end_time = time.time()
# Measure final memory usage
mem_after = memory_usage(-1, interval=0.1, timeout=1)[0]
# Calculate and print runtime
print("Execution time: ", end_time - start_time, "seconds")
# Calculate and print memory usage
print("Memory usage: ", mem_after - mem_before, "MiB")
return A_inv
# Example of use
A_inv = invert_matrix_cholesky(Matrix)
if A_inv is not None:
print(A_inv)
Процесс обращения матриц с помощью этого кода происходит следующим образом:
c = scipy.linalg.cholesky(Matrix): Эта строка выполняет разложение Холецкого на входной матрице. Разложение Холецкого разлагает входную матрицу на нижнюю треугольную матрицу C и ее транспонирование.c_inv = scipy.linalg.inv(c): Эта строка вычисляет обратную нижнюю треугольную матрицу 'c', полученную в результате разложения Холецкого. Обращение треугольной матрицы вычислительно более эффективно, чем обращение произвольной матрицы.A_inv = np.dot(c_inv.T, c_inv): Эта строка вычисляет обратную исходной матрице Matrix. Свойство разложения Холецкого состоит в том, что обратная исходная матрица равна произведению обратной верхней треугольной матрицы (которая является транспонированной c_inv) и обратной нижней треугольной матрицы. c_inv. Функция np.dot() используется для выполнения матричного умножения.
3. Метод QR-разложения
Разложение QR разлагает матрицу на ортогональную матрицу (q) и верхнюю треугольную матрицу (r). Этот метод полезен для решения систем линейных уравнений, вычисления собственных значений и собственных векторов, среди других приложений.
Здесь мы будем использовать функцию scipy.linalg.qr() из той же библиотеки.
import numpy as np
import scipy.linalg
import time
from memory_profiler import memory_usage
def invert_matrix_qr(Matrix):
n = Matrix.shape[0]
# Measure initial memory usage
mem_before = memory_usage(-1, interval=0.1, timeout=1)[0]
# Measure the start time
start_time = time.time()
try:
# Doing the QR decomposition
q, r = scipy.linalg.qr(Matrix)
# Inverting the r matrix
r_inv = scipy.linalg.inv(r)
# Getting the inverse of Matrix
A_inv = np.dot(r_inv, q.T)
except np.linalg.LinAlgError:
print("Matrix is singular and cannot be inverted, or is not suitable for QR decomposition.")
return None
# Measure the end time
end_time = time.time()
# Measure ultimate memory usage
mem_after = memory_usage(-1, interval=0.1, timeout=1)[0]
# Calculate and print runtime
print("Execution time: ", end_time - start_time, "seconds")
# Calculate and print memory usage
print("Memory usage: ", mem_after - mem_before, "MiB")
return A_inv
# Example of use
A_inv = invert_matrix_qr(Matrix)
if A_inv is not None:
print(A_inv)
Процесс обращения матриц с помощью этого кода происходит следующим образом:
q, r = scipy.linalg.qr(A): Эта строка выполняет QR-разложение входной матрицы A. Функция scipy.linalg.qr() возвращает две матрицы: 'q', которая является ортогональной матрицей, и 'r', которая является верхней треугольной матрицей.r_inv = scipy.linalg.inv(r): Эта строка вычисляет обратную матрицу R, полученную в результате QR-разложения. Поскольку «r» — верхняя треугольная матрица, эта операция более эффективна, чем обращение произвольной матрицы.A_inv = np.dot(r_inv, q.T): Эта строка вычисляет обратную исходную матрицу A. Согласно свойствам разложения QR, обратная матрица A является произведением обратной матрицы 'r' и транспонированной матрицы 'q'. Функция np.dot() используется для выполнения матричного умножения.
4. Метод Гаусса-Жордана
Метод Гаусса-Жордана используется для нахождения обратной матрицы путем приведения данной матрицы к тождественной форме. Хотя он может быть менее эффективным, чем некоторые другие упомянутые методы, это прямой метод нахождения обратной матрицы.
Для этого метода вычисления могут стать более сложными, и мы будем использовать только библиотеку numpy:
import numpy as np
import time
from memory_profiler import memory_usage
def invert_matrix_gauss_jordan(Matrix):
n = Matrix.shape[0]
# Measure initial memory usage
mem_before = memory_usage(-1, interval=0.1, timeout=1)[0]
# Measure the start time
start_time = time.time()
try:
# Creates an augmented matrix [Matrix|I]
augmented_matrix = np.hstack((Matrix, np.identity(n)))
for i in range(n):
# Normalizing line i
augmented_matrix[i] = augmented_matrix[i] / augmented_matrix[i, i]
for j in range(n):
if i != j:
# Subtracting from row j the row i times the element j, i of the matrix
augmented_matrix[j] -= augmented_matrix[i] * augmented_matrix[j, i]
# Separating the inverse matrix from the augmented matrix
A_inv = augmented_matrix[:, n:]
except np.linalg.LinAlgError:
print("Matrix is singular and cannot be inverted.")
return None
# Measure the end time
end_time = time.time()
# Measure final memory usage
mem_after = memory_usage(-1, interval=0.1, timeout=1)[0]
# Calculate and print runtime
print("Execution time: ", end_time - start_time, "seconds")
# Calculate and print memory usage
print("Memory usage: ", mem_after - mem_before, "MiB")
return A_inv
# Example of use
A_inv = invert_matrix_gauss_jordan(Matrix)
if A_inv is not None:
print(A_inv)
Вот краткое объяснение кода:
augmented_matrix = np.hstack((A, np.identity(n))): Эта строка создает расширенную матрицу путем объединения исходной матрицы A и единичной матрицы I размера n x n рядом (по горизонтали). Расширенная матрица представляет собой систему уравнений AX = I, где X — обратная матрица, которую мы хотим найти.- Первый цикл
for i in range(n):перебирает каждую строку матрицы.
augmented_matrix[i] = augmented_matrix[i] / augmented_matrix[i, i]: Эта строка нормализует i-ю строку матрицы путем деления всей строки на диагональный элемент. Это гарантирует, что диагональный элемент станет равным 1.- Второй цикл
for j in range(n):повторяет каждую строку снова.
if i != j: Если строка j отличается от строки i (мы не хотим изменять строку, которую мы только что нормализовали)… augmented_matrix[j] -= augmented_matrix[i] * augmented_matrix[j, i]: мы вычитаем из строки j строку i, умноженную на элемент в позиции (j, i). Это гарантирует, что все элементы за пределами диагонали станут равными 0, что приведет к получению матрицы в сокращенной эшелонированной форме строки Гаусса-Жордана.A_inv = augmented_matrix[:, n:]Эта линия отделяет обратную матрицу от расширенной матрицы. В конце процесса Гаусса-Жордана расширенная матрица будет иметь вид [I|X], где I — единичная матрица, а X — обратная матрица. Следовательно, мы можем получить обратную матрицу, взяв все столбцы после n-го столбца.
Это алгоритм нахождения обратной матрицы, который применяет ряд элементарных операций для сведения заданной матрицы к единичной матрице.
5. Метод Монтанте
Метод Монтанте (также известный как алгоритм Барейса) представляет собой разновидность метода исключения Гаусса, используемого для решения систем линейных уравнений. Как и метод Гаусса-Жордана, метод Монтанте включает в себя операции со строками и столбцами для преобразования матрицы в форму, в которой решение можно легко идентифицировать (например, диагональную форму или форму эшелона строк). Следовательно, вычислительные затраты быстро возрастают с размером матрицы.
Однако преимущество метода Монтанте перед методом Гаусса-Жордана заключается в том, что он позволяет избежать операций деления, что может помочь уменьшить ошибки округления и сделать метод более стабильным в численном отношении. Но он до сих пор обычно не используется в крупномасштабных численных вычислениях из-за его кубической временной сложности.
Здесь мы снова будем использовать только Numpy:
import numpy as np
import time
from memory_profiler import memory_usage
def invert_matrix_montant(Matrix):
n = Matrix.shape[0]
# Measure the start time
mem_before = memory_usage(-1, interval=0.1, timeout=1)[0]
# Measure the start time
start_time = time.time()
# Creates an augmented matrix [Matrix|I]
augmented_matrix = np.hstack((Matrix, np.identity(n)))
try:
# Loop over each line
for k in range(n):
pivot = augmented_matrix[k, k]
augmented_matrix[k] = augmented_matrix[k] / pivot
for i in range(n):
if i != k:
factor = augmented_matrix[i, k]
augmented_matrix[i] -= factor * augmented_matrix[k]
# Separating the inverse matrix from the augmented matrix
A_inv = augmented_matrix[:, n:]
except np.linalg.LinAlgError:
print("Matrix is singular and cannot be inverted.")
return None
# Measure the end time
end_time = time.time()
# Measure final memory usage
mem_after = memory_usage(-1, interval=0.1, timeout=1)[0]
# Calculate and print runtime
print("Execution time: ", end_time - start_time, "seconds")
# Calculate and print memory usage
print("Memory usage: ", mem_after - mem_before, "MiB")
return A_inv
# Example of use
A_inv = invert_matrix_montant(Matrix)
if A_inv is not None:
print(A_inv)
for k in range(n):— это внешний цикл, который перебирает каждую строку матрицы.pivot = augmented_matrix[k, k]— Эта строка выбирает опорный элемент, который является диагональным элементом строки k.augmented_matrix[k] = augmented_matrix[k] / pivot— Эта линия делит всю строку k на опорный элемент, тем самым нормализуя опорный элемент до 1.for i in range(n):— это внутренний цикл, который перебирает все строки матрицы.if i != k:— проверяет, отличается ли текущая строка i от строки k. Мы не хотим изменять строку k, которую мы только что нормализовали.factor = augmented_matrix[i, k]— Эта строка получает коэффициент, который будет использоваться для исключения элемента в строке i и столбце k.augmented_matrix[i] -= factor * augmented_matrix[k]— Эта строка вычитает строку k, умноженную на коэффициент, из строки i. Это приводит к нулевому элементу в позиции (i, k), потому что мы вычитаем элемент (i, k) из строки k (где элемент (k, k) равен единице).A_inv = augmented_matrix[:, n:]— Наконец, после обработки всех строк расширенная матрица представляет собой единичную матрицу, объединенную с обратной исходной матрицей. Эта строка отделяет обратную матрицу от расширенной матрицы.
Давайте сравним методы…
Коды установлены, теперь, чтобы иметь представление о прогрессии, давайте инвертируем матрицы размером 2500, 5000, 7500 и 10000 строк и столбцов.
Начнем со времени вычислений для 4 протестированных методов:
Рисунок 1 —Время в секундах для инверсии матрицы.
Таблица 1 — Время в секундах для инверсии матрицы.
Мы хорошо видим, что метод Холецкого является самым быстрым в обращении матриц, и здесь мы видим не маленькую разницу, а скорее огромный разрыв в скорости обращения. Учитывая матрицу 10 000 x 10 000, метод Холецкого занял менее 5% времени по сравнению с методом Гаусса-Джордана и только половину времени по сравнению со вторым по скорости методом (QR-разложение). Однако является ли этот метод тем, который потребляет наименьшее количество оперативной памяти?
Рисунок 2 – Потребление оперативной памяти для инверсии матрицы.
Таблица 2 — Потребление оперативной памяти для инверсии матриц.
Здесь мы видим, что медленные методы Гаусса-Жордана и Монтанте потребляют почти половину памяти по сравнению с методом QR-разложения, а быстрый метод Холецкого не впечатляет по потреблению памяти. При рассмотрении обоих показателей вычислительных затрат, которые мы использовали для анализа затрат и результатов, мы можем выделить метод разложения LU и даже метод Холецкого.
Заключение
Результаты здесь могут различаться в зависимости от того, что вы делаете, но я надеюсь, что это небольшое сравнение прольет свет на ваши исследования и поможет вам выбрать лучший метод для вашего исследования. Например, если вам нужна скорость и имеется отличное оборудование, выберите метод Холецкого. Однако, если ваше оборудование ограничено и вы не возражаете против длительного ожидания, методы Гаусса-Жордана и Монтанте могут быть хорошим вариантом.
Важно отметить, что хотя эти методы точны при обращении матриц, они могут быть не самыми эффективными при работе с очень большими матрицами. В случае больших объемов данных вы можете рассмотреть методы, которые жертвуют некоторой точностью в обмен на меньшее время и потребление памяти, такие как итерационные или приближенные методы для обращения матриц.
