函数scipy.sparse.linalg.spilu(A).solve()使用什么算法?

iibxawm4  于 2022-11-10  发布在  其他
关注(0)|答案(3)|浏览(258)

请考虑以下代码:

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 ?

这导致xAx = b的解。
我的理解是scipy.sparse.linalg.spilu(A)计算的是A的不完全LU分解。
我的问题是:spilu.solve(b)使用什么算法来求解Ax = b
我希望它使用一些迭代方法,如共轭梯度法,因为这似乎是使用不完全LU分解的正常方式。然而,我一直无法找到支持或不同意这一点的文档。此外,我感到困惑,因为我看到一些人将LinearOperatorscipy.sparse.linalg.cg/scipy.sparse.linalg.cgscipy.sparse.linalg.spilu(A)结合使用。如果我的假设是正确的(for example),那么这看起来是愚蠢的。

wecizke3

wecizke31#

医生说:
此函数使用SuperLU库。
代码通过:

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

它只是一个SuperLU Package 器:

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)

其指示整个完整的逻辑与不完整的逻辑被传递到SuperLU。
现在,我们来看看SuperLU's manual
2.7不完全LU分解(ILU)预条件子
从SuperLU版本4.0开始,我们提供了ILU例程,用作迭代求解器的预处理器。我们的ILU方法可以被认为是最初由Saad [31]提出的ILUTP方法的变体,它结合了双重丢弃策略和数值主元(“T”代表阈值,“P”代表主元)。
参考文献[31]是:
萨阿德,优素福。“ILUT:双阈值不完全LU分解〉,《数值线性代数与应用》1.4(1994):387-402.

在我看来这是一个DIRECT方法Eigen does use it too),但我认为您将能够查找您需要了解的内容。

lqfhib0f

lqfhib0f2#

此求解方法 * 通过使用A的不完全因式分解近似求解 * Ax=b。这就是为什么您会看到它在那些迭代方法中用作预处理器的原因。此求解方法是使用不完全因式的经典稀疏前向/后向替换算法,但由于它们是不完全的,因此不太可能给予准确的结果。

70gysomp

70gysomp3#

萨莎补充道:
_superlu.gstrf返回的对象的solve函数可以在scipy/sparse/linalg/_dsolve/_superluobject. c中找到。主要部分是对gstrs的调用(请注意,它是gstr**s**而不是gstr**f**):

static PyObject *SuperLU_solve(SuperLUObject * self, PyObject * args,
                               PyObject * kwds)
{
...
    gstrs(self->type,
          trans, &self->L, &self->U, self->perm_c, self->perm_r,
          (SuperMatrix *)&B, (SuperLUStat_t *)&stat, (int *)&info);
...
}

查看documentation of SuperLU

DGSTRS solves a system of linear equations A*X=B or A'*X=B
with A sparse and B dense, using the LU factorization computed by
DGSTRF.

所以Python代码基本上是调用*gstrf来计算LU分解,然后把它给*gstrs来计算一个系统的解。
看一下dgstrs的源代码,看起来它只是使用了LU分解,不管它是否是不完整的(也没有参数允许指定),所以它只是用不完整的分解给予你一个近似的解。

相关问题