Какой алгоритм использует функция scipy.sparse.linalg.spilu (A) .solve ()?

Question

Рассмотрим следующий код:

import numpy as np
import scipy.sparse.linalg

# Setup
A = scipy.sparse.csc_matrix([[1, -3], [-1, 4]])
b = np.array([1, 0])

spilu = scipy.sparse.linalg.spilu(A)  # Find ILU decomposition
x = spilu.solve(b)  # Use an iterative method to solve Ax = b ?

В результате x является решением Ax = b.

Насколько я понимаю, scipy.sparse.linalg.spilu(A) вычисляет неполное LU-разложение A.

Мой вопрос: какой именно алгоритм spilu.solve(b) использует для решения для Ax = b?

Я ожидал бы, что он будет использовать некоторый итерационный метод, такой как метод сопряженного градиента, потому что это, кажется, нормальный способ использовать неполное разложение LU. Однако я не смог найти документацию, подтверждающую или не согласную с этим. Кроме того, я запутался, потому что вижу, что некоторые люди используют LinearOperator и scipy.sparse.linalg.cg / scipy.sparse.linalg.cg в сочетании с scipy.sparse.linalg.spilu(A), что может показаться глупым sh, если моя гипотеза верна ( например ).

Answer 1

Документы говорят:

Эта функция использует библиотеку SuperLU.

Код проходит через:

return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
                      csc_construct_func=csc_construct_func,
                      ilu=True, options=_options)

, то есть просто оболочка SuperLU:

static char gstrf_doc[] = "gstrf(A, ...)\n\
\n\
performs a factorization of the sparse matrix A=*(N,nnz,nzvals,rowind,colptr) and \n\
returns a factored_lu object.\n\
\n\
arguments\n\
---------\n\
\n\
Matrix to be factorized is represented as N,nnz,nzvals,rowind,colptr\n\
  as separate arguments.  This is compressed sparse column representation.\n\
\n\
N         number of rows and columns \n\
nnz       number of non-zero elements\n\
nzvals    non-zero values \n\
rowind    row-index for this column (same size as nzvals)\n\
colptr    index into rowind for first non-zero value in this column\n\
          size is (N+1).  Last value should be nnz. \n\
\n\
additional keyword arguments:\n\
-----------------------------\n\
options             specifies additional options for SuperLU\n\
                    (same keys and values as in superlu_options_t C structure,\n\
                    and additionally 'Relax' and 'PanelSize')\n\
\n\
ilu                 whether to perform an incomplete LU decomposition\n\
                    (default: false)\n\
";

посмотрите на последний аргумент -> ilu!

splu выглядит следующим образом:

return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
                      csc_construct_func=csc_construct_func,
                      ilu=False, options=_options)

, что указывает на то, что весь полный и неполный логи c передается в SuperLU.

Теперь давайте рассмотрим SuperLU's manual :

2.7 Прекондиционер неполной LU-факторизации (ILU)

Начиная с версии SuperLU 4.0, мы предоставляем подпрограммы ILU, которые будут использоваться в качестве предобусловливателей для итерационных решателей. Наш метод ILU можно рассматривать как вариант метода ILUTP, первоначально предложенного Саадом [31], который сочетает в себе стратегию двойного отбрасывания с числовым поворотом («T» обозначает порог, а «P» обозначает поворот).

Ссылка [31]:

Саад, Юсеф. «ILUT: двухпороговая неполная факторизация LU». Числовая линейная алгебра с приложениями 1.4 (1994): 387-402.

Для меня это выглядит ПРЯМЫМ методом ( Эйген тоже использует его ) , но я думаю, что вы сможете найти то, что вам нужно знать.