Оптимизировать решение A * x = B для матрицы трехдиагональных коэффициентов

У меня есть система уравнений в виде A*x = B где [A] – матрица трехдиагональных коэффициентов. Используя Решатель numpy.linalg.solve я могу решить систему уравнений для x.

См. Пример ниже о том, как я развиваю трехдиагональный [A] мартикс. вектор {B} и решить для x :

 # Solve system of equations with a tridiagonal coefficient matrix # uses numpy.linalg.solve # use Python 3 print function from __future__ import print_function from __future__ import division # modules import numpy as np import time ti = time.clock() #---- Build [A] array and {B} column vector m = 1000 # size of array, make this 8000 to see time benefits A = np.zeros((m, m)) # pre-allocate [A] array B = np.zeros((m, 1)) # pre-allocate {B} column vector A[0, 0] = 1 A[0, 1] = 2 B[0, 0] = 1 for i in range(1, m-1): A[i, i-1] = 7 # node-1 A[i, i] = 8 # node A[i, i+1] = 9 # node+1 B[i, 0] = 2 A[m-1, m-2] = 3 A[m-1, m-1] = 4 B[m-1, 0] = 3 print('A \n', A) print('B \n', B) #---- Solve using numpy.linalg.solve x = np.linalg.solve(A, B) # solve A*x = B for x print('x \n', x) #---- Elapsed time for each approach print('NUMPY time', time.clock()-ti, 'seconds') 

Поэтому мой вопрос касается двух разделов приведенного выше примера:

  1. Поскольку я имею дело с трехдиагональной матрицей для [A] , также называемой полосой матрицы, существует ли более эффективный способ решения системы уравнений вместо использования numpy.linalg.solve ?
  2. Кроме того, есть ли лучший способ создать трехдиагональную матрицу вместо использования for-loop ?

Вышеприведенный пример работает в Linux примерно через 0.08 seconds соответствии с time.clock() .

Функция numpy.linalg.solve работает нормально, но я пытаюсь найти подход, который использует тридиагональную форму [A] в надежде ускорить решение еще дальше, а затем применить этот подход к более сложному примеру.

Есть два моментальных улучшения производительности (1), которые не используют цикл, (2) используют scipy.linalg.solve_banded() .

Я бы написал код, более похожий

 import scipy.linalg as la # Create arrays and set values ab = np.zeros((3,m)) b = 2*ones(m) ab[0] = 9 ab[1] = 8 ab[2] = 7 # Fix end points ab[0,1] = 2 ab[1,0] = 1 ab[1,-1] = 4 ab[2,-2] = 3 b[0] = 1 b[-1] = 3 return la.solve_banded ((1,1),ab,b) 

Могут быть более элегантные способы построения матрицы, но это работает.

Используя %timeit в ipython исходный код занял 112 мс для m = 1000. Этот код занимает 2,94 мс для m = 10 000, на порядок больше, но все же почти на два порядка быстрее! У меня не было терпения ждать по исходному коду для m = 10 000. В большинстве случаев в оригинале может быть построение массива, я не тестировал это. Несмотря на это, для больших массивов гораздо эффективнее хранить ненулевые значения матрицы.

Существует матричный тип scipy.sparse.dia_matrix называемый scipy.sparse.dia_matrix который фиксирует структуру вашей матрицы (он будет хранить 3 массива, в «положениях» 0 (диагональ), 1 (выше) и -1 (ниже)). Используя этот тип матрицы, вы можете попробовать scipy.sparse.linalg.lsqr для решения. Если ваша проблема имеет точное решение, она будет найдена, иначе она найдет решение в смысле наименьших квадратов.

 from scipy import sparse A_sparse = sparse.dia_matrix(A) ret_values = sparse.linalg.lsqr(A_sparse, C) x = ret_values[0] 

Однако это может быть не совсем оптимальным с точки зрения использования триадиагональной структуры, может быть теоретический способ сделать это быстрее. То, что это преобразование делает для вас, сокращает расходы на умножение матрицы на существенное: используются только 3 полосы. Это, в сочетании с итеративным решателем lsqr должно уже дать ускорение.

Примечание . Я не предлагаю scipy.sparse.linalg.spsolve , потому что он преобразует вашу матрицу в формат csr . Однако замена lsqr на spsolve стоит попробовать, особенно потому, что spsolve может связывать UMFPACK , см. Соответствующий документ о spsolve . Кроме того, может быть интересно взглянуть на этот вопрос и ответ на stackoverflow, относящийся к UMFPACK

Вы можете использовать scipy.linalg.solveh_banded .

EDIT : вы НЕ МОЖЕТЕ использовать вышеупомянутое, так как ваша матрица не симметрична, и я подумал, что это так. Однако, как было упомянуто выше в комментарии, алгоритм Томаса отлично подходит для этого

 a = [7] * ( m - 2 ) + [3] b = [1] + [8] * ( m - 2 ) + [4] c = [2] + [9] * ( m - 2 ) d = [1] + [2] * ( m - 2 ) + [3] # This is taken directly from the Wikipedia page also cited above # this overwrites b and d def TDMASolve(a, b, c, d): n = len(d) # n is the numbers of rows, a and c has length n-1 for i in xrange(n-1): d[i+1] -= 1. * d[i] * a[i] / b[i] b[i+1] -= 1. * c[i] * a[i] / b[i] for i in reversed(xrange(n-1)): d[i] -= d[i+1] * c[i] / b[i+1] return [d[i] / b[i] for i in xrange(n)] 

Этот код не оптимизирован и не использует np , но если у меня (или у любого из других прекрасных людей) есть время, я отредактирую его так, чтобы он делал это. В настоящее время он составляет ~ 10 мс для m = 10000.