时间戳

首页 > 解决方案 > 如何在 scipy.optimize.newton_krylov 中更新预处理器

python - 如何在 scipy.optimize.newton_krylov 中更新预处理器

问题描述

我正在尝试解决一个示例非线性问题，并且我想为 scipy.optimize 的 newton_krylov 求解器设置一个预处理器。

我使用 spilu 由 LinearOperator 设置预调节器。矩阵随着自变量而变化，即矩阵M(x)，并且M(x)需要像残差F(x)一样在每次非线性迭代中更新。

似乎 newton_krylov 中的“inner_M”选项仅在流程开始时调用一次 LinearOperator。如何设置“inner_M”预处理器选项以在每次牛顿迭代中执行 M(x) 更新？

具体来说，我无法理解 “inner_M”设置的https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.newton_krylov.html中的描述：

“...如果预条件器有一个名为 'update' 的方法，它将在每个非线性步骤之后被称为 update(x, f)，其中 x 给出当前点，f 给出当前函数值。”

非常感谢。

标签： pythonoptimizationscipy

解决方案

Following the example under doc/scipy/ ... still-too-slow-precondition, set M.update to a function which copies x to a global, here xglo. Or, make a class with self.xupdate.

    xglo = None

    def get_preconditioner():
        ...
        J1_ilu = spilu(J1)  # better spsolve ?

        M = LinearOperator(shape=(nx*ny, nx*ny), matvec=J1_ilu.solve)

        #.......................................................................
        def xglobal( x, F ):
            """ newton_krylov calls this at each iteration: xglo = x """
            global xglo
            xglo = x.copy()
            print( "update: x %s  F %s " % (nu.asum(x), nu.asum(F)) )  # test
                # nu.asum: array summary, size min av max

        M.update = xglobal

        return M

(I don't know if this is really a good idea -- depends on your preconditioner.
scicomp.stack might be a better place to ask, see e.g. when-is-newton-krylov-not-an-appropriate-solver .)

推荐阅读