решение большой разреженной линейной системы, хуже с переупорядочением и предварительным кондиционированием? - PullRequest
Новая
       88

решение большой разреженной линейной системы, хуже с переупорядочением и предварительным кондиционированием?

2 голосов
/ 04 мая 2020

У меня есть линейная система с матрицей 60000x60000, которую я буду решать sh, с около 6 000 000 ненулевых записей в ней.

Мой текущий подход заключается в переупорядочении матрицы с помощью обратного среза mckee, разложить матрица, а затем решить ее с помощью предопределенного сопряженного градиента, но я не получаю очень хорошие результаты, и я не понимаю, почему. Изменение порядка выглядит разумным.

Ниже я привел простой пример, в котором я использую только подсистему матрицы, которую пытаюсь решить. 

import matplotlib
matplotlib.use('TkAgg') #TkAgg for vizual
from matplotlib import pyplot as plt
import time
import numpy as np
import scipy
from scipy import sparse
from scipy.sparse.linalg import LinearOperator, spilu, cg
from numpy.linalg import norm
L = sparse.load_npz("L_Matrix.npz")
n = 20000
b = np.random.randn((n))
L2 = L[0:n,0:n].copy()

t00 = time.time()
perm = scipy.sparse.csgraph.reverse_cuthill_mckee(L2, symmetric_mode=True)
I,J = np.ix_(perm,perm)
bp = b[perm]
L2p = L2[I, J]
t01 = time.time()

fig = plt.figure(0, figsize=[20, 10])
plt.subplot(1, 2, 1)
plt.spy(L2)
plt.subplot(1, 2, 2)
plt.spy(L2p)
plt.pause(1)
# plt.pause(1)

t0 = time.time()
print("reordering took {}".format(t0-t00))
ilu = spilu(L2p)
t1 = time.time()
print("Factorization took {}".format(t1-t0))

Mx = lambda x: ilu.solve(x)
M = LinearOperator((n, n), Mx)
x,stat = cg(L2p, bp, tol=1e-12, maxiter=500, M=M)
t2 = time.time()
print("pcg took {} s, and had status {}".format(t2-t1,stat))
print("reorder+pcg+factor = {} s".format(t2-t00))
bsol = L2p @ x
R = norm(bsol - bp)
print("pcg residual = {}".format(R))

x,stat = cg(L2, b, tol=1e-12, maxiter=500)
t3 = time.time()
print("cg took {} s, and had status {}".format(t3-t2,stat))
bsol = L2 @ x
R = norm(bsol - b)
print("pcg residual = {}".format(R))

Результаты, которые я получаю из этого являются: 

reordering took 66.32699060440063
Factorization took 64.96741151809692
pcg took 12.732918739318848 s, and had status 500
reorder+pcg+factor = 144.0273208618164 s
pcg residual = 29.10655954230801
cg took 1.2132720947265625 s, and had status 500
pcg residual = 2.5236861383747353

Таким образом, не только переупорядочение и факторизация занимают очень много времени, но и решение с помощью cg не go быстрее и не дает правильного решения, в На самом деле остаток намного хуже!

Может кто-нибудь сказать мне, что именно я здесь делаю неправильно?

1 Ответ

2 голосов
/ 04 мая 2020

Существует высокая вероятность того, что вы не делаете ничего плохого в рамках вашего текущего подхода (по крайней мере, я не смог обнаружить очевидную ошибку).

Пара замечаний:

  1. Остаток 29.10655954230801 и 2.5236861383747353 после 500 итераций фактически одинаковы: ваше итеративное решение не сходилось.
  2. Похоже, вы запрашиваете очень высокий допуск итерационного решателя 1E-12. Это не имеет значения, так как у вас есть проблема, которая не сходится вообще.
  3. Факторизация (ILU) должна занять примерно это время. Я не удивлен, увидев этот номер для такой системы. Не очень знаком с этой реализацией Cuthill-McKee.

Не зная, откуда взялась ваша система, было бы очень сложно что-либо сказать. Однако:

  1. Проверьте номер условия для вашей уменьшенной версии матрицы (если она в некоторой степени отражает вашу первоначальную проблему). Большое число условий указывает на проблему с подготовкой матрицы; таким образом, потенциальная плохая сходимость или плохая сходимость итеративного решения (или любого типа решения в крайнем случае).
  2. Градиент сопряженности предназначен для систем, которые симметричны c и положительно определенные . Может сходиться для других случаев; однако, это может потерпеть неудачу для хорошо обусловленных проблем, которые не являются положительно определенными.
...