Решение разреженного Ax=b в scipy
Мне нужно решить Ax=b
где A
является матрицей, которая представляет метод конечных разностей для PDE. Типичный размер A
для двумерной задачи около (256^2) х (256^2), и она состоит из нескольких диагоналей. Следующий пример кода, как я строю A:
N = Nx*Ny # Nx is no. of cols (size in x-direction), Ny is no. rows (size in y-direction)
# finite difference in x-direction
up1 = (0.5)*c
up1[Nx-1::Nx] = 0
down1 = (-0.5)*c
down1[::Nx] = 0
matX = diags([down1[1:], up1[:-1]], [-1,1], format='csc')
# finite difference in y-direction
up1 = (0.5)*c
down1 = (-0.5)*c
matY = diags([down1[Nx:], up1[:N-Nx]], [-Nx,Nx], format='csc')
Добавление matX
а также matY
вместе получается четыре диагонали. Выше для дискретизации второго порядка. Для дискретизации четвертого порядка у меня есть восемь диагоналей. Если у меня есть вторая производная, то главная диагональ также отлична от нуля.
Я использую следующий код для фактического решения:
# Initialize A_fixed, B_fixed
if const is True: # the potential term V(x) is time-independent
A = A_fixed + sp.sparse.diags(V_func(x))
B = B_fixed + sp.sparse.diags(V_func(x))
A_factored = sp.sparse.linalg.factorized(A)
for idx, t in enumerate(t_steps[1:],1):
# Solve Ax=b=Bu
if const in False: #
A = A_fixed + sp.sparse.diags(V_func(x,t))
B = B_fixed + sp.sparse.diags(V_func(x,t))
psi_n = B.dot(psi_old)
if const is True:
psi_new = A_factored(psi_n)
else:
psi_new = sp.sparse.linalg.spsolve(A,psi_n,use_umfpack=False)
psi_old = psi_new.copy()
У меня есть пара вопросов:
- Какой лучший способ решить
Ax=b
в scipy? Я используюspsolve
вsp.sparse.linalg
библиотека, которая использует LU-декомпозицию. Я пытался использовать стандартsp.linalg.solve
для плотной матрицы, но это значительно медленнее. Я также попытался использовать все другие итерационные решатели вsp.sparse.linalg
библиотека, но они также медленнее (для 1000x1000 все они занимают пару секунд, по сравнению с менее чем полсекунды дляspsolve
, и мойA
скорее всего будет намного больше). Есть ли альтернативные способы сделать вычисления?
Изменить: проблема, которую я пытаюсь решить, на самом деле зависит от времени уравнения Шредингера. Если потенциальный член не зависит от времени, то, как и предполагалось, я могу предварительно разложить матрицу A
первый, чтобы ускорить код, но это не работает, если потенциал изменяется во времени, так как мне нужно изменить диагональный член обеих матриц A
а также B
на каждом временном шаге. Для этой конкретной проблемы, есть ли способы ускорить код, используя метод, аналогичный предварительной факторизации или другие способы?
- Я установил
umfpack
, Я пробовал настройкуuse_umfpack
Правда и Ложь, чтобы проверить это, но на удивлениеuse_umfpack=True
занимает почти вдвое больше времени, чемuse_umfpack=False
, Я думал, что наличие этого пакета должно ускорить. Есть идеи, в чем дело? (PS: я использую Anaconda Spyder для запуска кода, если это имеет значение)
У меня есть использование cProfile
профилировать мои коды, и почти все время тратится на эту линию с spsolve
, Поэтому я думаю, что оставшаяся часть кода (инициализация матрицы / задачи) в значительной степени оптимизирована, и это часть для решения матриц, которую необходимо улучшить.
Благодарю.