Add initial_guess parameter to apply_inverse

src/pymor/algorithms/genericsolvers.py

@@ -97,18 +97,22 @@ def solver_options(lgmres_tol=1e-5,
 
 
 @defaults('check_finite', 'default_solver', 'default_least_squares_solver')
-def apply_inverse(op, rhs, options=None, least_squares=False, check_finite=True,
+def apply_inverse(op, V, initial_guess=None, options=None, least_squares=False, check_finite=True,
                   default_solver='generic_lgmres', default_least_squares_solver='generic_least_squares_lsmr'):
     """Solve linear equation system.
 
-    Applies the inverse of `op` to the vectors in `rhs` using a generic iterative solver.
+    Applies the inverse of `op` to the vectors in `V` using a generic iterative solver.
 
     Parameters
     ----------
     op
         The linear, non-parametric |Operator| to invert.
-    rhs
+    V
         |VectorArray| of right-hand sides for the equation system.
+    initial_guess
+        |VectorArray| with the same length as `V` containing initial guesses
+        for the solution.  Some implementations of `apply_inverse` may
+        ignore this parameter.  If `None` a solver-dependent default is used.
     options
         The |solver_options| to use (see :func:`solver_options`).
     least_squares
@@ -126,13 +130,15 @@ def apply_inverse(op, rhs, options=None, least_squares=False, check_finite=True,
     |VectorArray| of the solution vectors.
     """
 
+    assert V in op.range
+    assert initial_guess is None or initial_guess in op.source and len(initial_guess) == len(V)
     options = _parse_options(options, solver_options(), default_solver, default_least_squares_solver, least_squares)
 
-    R = op.source.empty(reserve=len(rhs))
+    R = op.source.empty(reserve=len(V))
 
     if options['type'] == 'generic_lgmres':
-        for i in range(len(rhs)):
-            r, info = lgmres(op, rhs[i],
+        for i in range(len(V)):
+            r, info = lgmres(op, V[i], initial_guess[i] if initial_guess is not None else None,
                              tol=options['tol'],
                              maxiter=options['maxiter'],
                              inner_m=options['inner_m'],
@@ -142,8 +148,8 @@ def apply_inverse(op, rhs, options=None, least_squares=False, check_finite=True,
             assert info == 0
             R.append(r)
     elif options['type'] == 'generic_least_squares_lsmr':
-        for i in range(len(rhs)):
-            r, info, itn, _, _, _, _, _ = lsmr(op, rhs[i],
+        for i in range(len(V)):
+            r, info, itn, _, _, _, _, _ = lsmr(op, V[i],
                                                damp=options['damp'],
                                                atol=options['atol'],
                                                btol=options['btol'],
@@ -156,8 +162,8 @@ def apply_inverse(op, rhs, options=None, least_squares=False, check_finite=True,
             getLogger('pymor.algorithms.genericsolvers.lsmr').info(f'Converged after {itn} iterations')
             R.append(r)
     elif options['type'] == 'generic_least_squares_lsqr':
-        for i in range(len(rhs)):
-            r, info, itn, _, _, _, _, _, _ = lsqr(op, rhs[i],
+        for i in range(len(V)):
+            r, info, itn, _, _, _, _, _, _ = lsqr(op, V[i],
                                                   damp=options['damp'],
                                                   atol=options['atol'],
                                                   btol=options['btol'],

src/pymor/algorithms/timestepping.py

@@ -171,7 +171,7 @@ def implicit_euler(A, F, M, U0, t0, t1, nt, mu=None, num_values=None, solver_opt
             dt_F = F.as_vector(mu) * dt
         if F:
             rhs += dt_F
-        U = M_dt_A.apply_inverse(rhs, mu=mu)
+        U = M_dt_A.apply_inverse(rhs, mu=mu, initial_guess=U)
         while t - t0 + (min(dt, DT) * 0.5) >= len(R) * DT:
             R.append(U)
 

src/pymor/bindings/fenics.py

@@ -185,10 +185,12 @@ def _real_apply_adjoint_one_vector(self, v, mu=None, prepare_data=None):
             self.matrix.transpmult(v.impl, r.impl)
             return r
 
-        def _real_apply_inverse_one_vector(self, v, mu=None, least_squares=False, prepare_data=None):
+        def _real_apply_inverse_one_vector(self, v, mu=None, initial_guess=None,
+                                           least_squares=False, prepare_data=None):
             if least_squares:
                 raise NotImplementedError
-            r = self.source.real_zero_vector()
+            r = (self.source.real_zero_vector() if initial_guess is None else
+                 initial_guess.copy(deep=True))
             options = self.solver_options.get('inverse') if self.solver_options else None
             _apply_inverse(self.matrix, r.impl, v.impl, options)
             return r

src/pymor/bindings/ngsolve.py

@@ -134,7 +134,8 @@ def _real_apply_adjoint_one_vector(self, v, mu=None, prepare_data=None):
             mat.Mult(v.impl.vec, u.impl.vec)
             return u
 
-        def _real_apply_inverse_one_vector(self, v, mu=None, least_squares=False, prepare_data=None):
+        def _real_apply_inverse_one_vector(self, v, mu=None, initial_guess=None,
+                                           least_squares=False, prepare_data=None):
             inv = prepare_data
             r = self.source.real_zero_vector()
             r.impl.vec.data = inv * v.impl.vec

src/pymor/bindings/pyamg.py

@@ -152,17 +152,22 @@ def solver_options(tol=1e-5,
                                 'maxiter': sa_maxiter}}
 
     @defaults('check_finite', 'default_solver')
-    def apply_inverse(op, V, options=None, least_squares=False, check_finite=True, default_solver='pyamg_solve'):
+    def apply_inverse(op, V, initial_guess=None, options=None, least_squares=False,
+                      check_finite=True, default_solver='pyamg_solve'):
         """Solve linear equation system.
 
-        Applies the inverse of `op` to the vectors in `rhs` using PyAMG.
+        Applies the inverse of `op` to the vectors in `V` using PyAMG.
 
         Parameters
         ----------
         op
             The linear, non-parametric |Operator| to invert.
-        rhs
+        V
             |VectorArray| of right-hand sides for the equation system.
+        initial_guess
+            |VectorArray| with the same length as `V` containing initial guesses
+            for the solution.  Some implementations of `apply_inverse` may
+            ignore this parameter.  If `None` a solver-dependent default is used.
         options
             The |solver_options| to use (see :func:`solver_options`).
         least_squares
@@ -178,6 +183,7 @@ def apply_inverse(op, V, options=None, least_squares=False, check_finite=True, d
         """
 
         assert V in op.range
+        assert initial_guess is None or initial_guess in op.source and len(initial_guess) == len(V)
 
         if least_squares:
             raise NotImplementedError

src/pymor/bindings/scipy.py

@@ -149,18 +149,22 @@ def solver_options(bicgstab_tol=1e-15,
 
 
 @defaults('check_finite', 'default_solver', 'default_least_squares_solver')
-def apply_inverse(op, V, options=None, least_squares=False, check_finite=True,
+def apply_inverse(op, V, initial_guess=None, options=None, least_squares=False, check_finite=True,
                   default_solver='scipy_spsolve', default_least_squares_solver='scipy_least_squares_lsmr'):
     """Solve linear equation system.
 
-    Applies the inverse of `op` to the vectors in `rhs` using SciPy.
+    Applies the inverse of `op` to the vectors in `V` using SciPy.
 
     Parameters
     ----------
     op
         The linear, non-parametric |Operator| to invert.
-    rhs
+    V
         |VectorArray| of right-hand sides for the equation system.
+    initial_guess
+        |VectorArray| with the same length as `V` containing initial guesses
+        for the solution.  Some implementations of `apply_inverse` may
+        ignore this parameter.  If `None` a solver-dependent default is used.
     options
         The |solver_options| to use (see :func:`solver_options`).
     least_squares
@@ -180,6 +184,7 @@ def apply_inverse(op, V, options=None, least_squares=False, check_finite=True,
     """
 
     assert V in op.range
+    assert initial_guess is None or initial_guess in op.source and len(initial_guess) == len(V)
 
     if isinstance(op, NumpyMatrixOperator):
         matrix = op.matrix
@@ -190,12 +195,14 @@ def apply_inverse(op, V, options=None, least_squares=False, check_finite=True,
     options = _parse_options(options, solver_options(), default_solver, default_least_squares_solver, least_squares)
 
     V = V.to_numpy()
+    initial_guess = initial_guess.to_numpy() if initial_guess is not None else None
     promoted_type = np.promote_types(matrix.dtype, V.dtype)
     R = np.empty((len(V), matrix.shape[1]), dtype=promoted_type)
 
     if options['type'] == 'scipy_bicgstab':
         for i, VV in enumerate(V):
-            R[i], info = bicgstab(matrix, VV, tol=options['tol'], maxiter=options['maxiter'])
+            R[i], info = bicgstab(matrix, VV, initial_guess[i] if initial_guess is not None else None,
+                                  tol=options['tol'], maxiter=options['maxiter'])
             if info != 0:
                 if info > 0:
                     raise InversionError(f'bicgstab failed to converge after {info} iterations')
@@ -214,7 +221,8 @@ def apply_inverse(op, V, options=None, least_squares=False, check_finite=True,
                         permc_spec=options['spilu_permc_spec'])
         precond = LinearOperator(matrix.shape, ilu.solve)
         for i, VV in enumerate(V):
-            R[i], info = bicgstab(matrix, VV, tol=options['tol'], maxiter=options['maxiter'], M=precond)
+            R[i], info = bicgstab(matrix, VV, initial_guess[i] if initial_guess is not None else None,
+                                  tol=options['tol'], maxiter=options['maxiter'], M=precond)
             if info != 0:
                 if info > 0:
                     raise InversionError(f'bicgstab failed to converge after {info} iterations')
@@ -270,7 +278,7 @@ def apply_inverse(op, V, options=None, least_squares=False, check_finite=True,
             raise InversionError(e)
     elif options['type'] == 'scipy_lgmres':
         for i, VV in enumerate(V):
-            R[i], info = lgmres(matrix, VV,
+            R[i], info = lgmres(matrix, VV, initial_guess[i] if initial_guess is not None else None,
                                 tol=options['tol'],
                                 maxiter=options['maxiter'],
                                 inner_m=options['inner_m'],
@@ -287,7 +295,8 @@ def apply_inverse(op, V, options=None, least_squares=False, check_finite=True,
                                                   btol=options['btol'],
                                                   conlim=options['conlim'],
                                                   maxiter=options['maxiter'],
-                                                  show=options['show'])
+                                                  show=options['show'],
+                                                  x0=initial_guess[i] if initial_guess is not None else None)
             assert 0 <= info <= 7
             if info == 7:
                 raise InversionError(f'lsmr failed to converge after {itn} iterations')
@@ -299,7 +308,8 @@ def apply_inverse(op, V, options=None, least_squares=False, check_finite=True,
                                                         btol=options['btol'],
                                                         conlim=options['conlim'],
                                                         iter_lim=options['iter_lim'],
-                                                        show=options['show'])
+                                                        show=options['show'],
+                                                        x0=initial_guess[i] if initial_guess is not None else None)
             assert 0 <= info <= 7
             if info == 7:
                 raise InversionError(f'lsmr failed to converge after {itn} iterations')

src/pymor/operators/block.py

@@ -201,15 +201,23 @@ def apply_adjoint(self, V, mu=None):
         U_blocks = [self.blocks[i, i].apply_adjoint(V.block(i), mu=mu) for i in range(self.num_source_blocks)]
         return self.source.make_array(U_blocks)
 
-    def apply_inverse(self, V, mu=None, least_squares=False):
+    def apply_inverse(self, V, mu=None, initial_guess=None, least_squares=False):
         assert V in self.range
-        U_blocks = [self.blocks[i, i].apply_inverse(V.block(i), mu=mu, least_squares=least_squares)
+        assert initial_guess is None or initial_guess in self.source and len(initial_guess) == len(V)
+        U_blocks = [self.blocks[i, i].apply_inverse(V.block(i), mu=mu,
+                                                    initial_guess=(initial_guess.block(i)
+                                                                   if initial_guess is not None else None),
+                                                    least_squares=least_squares)
                     for i in range(self.num_source_blocks)]
         return self.source.make_array(U_blocks)
 
-    def apply_inverse_adjoint(self, U, mu=None, least_squares=False):
+    def apply_inverse_adjoint(self, U, mu=None, initial_guess=None, least_squares=False):
         assert U in self.source
-        V_blocks = [self.blocks[i, i].apply_inverse_adjoint(U.block(i), mu=mu, least_squares=least_squares)
+        assert initial_guess is None or initial_guess in self.range and len(initial_guess) == len(U)
+        V_blocks = [self.blocks[i, i].apply_inverse_adjoint(U.block(i), mu=mu,
+                                                            initial_guess=(initial_guess.block(i)
+                                                                           if initial_guess is not None else None),
+                                                            least_squares=least_squares)
                     for i in range(self.num_source_blocks)]
         return self.range.make_array(V_blocks)
 
@@ -285,15 +293,17 @@ def apply_adjoint(self, V, mu=None):
                     V.block(0) - self.E.apply_adjoint(V.block(1), mu=mu)]
         return self.source.make_array(U_blocks)
 
-    def apply_inverse(self, V, mu=None, least_squares=False):
+    def apply_inverse(self, V, mu=None, initial_guess=None, least_squares=False):
         assert V in self.range
+        assert initial_guess is None or initial_guess in self.source and len(initial_guess) == len(V)
         U_blocks = [-self.K.apply_inverse(self.E.apply(V.block(0), mu=mu) + V.block(1), mu=mu,
                                           least_squares=least_squares),
                     V.block(0)]
         return self.source.make_array(U_blocks)
 
-    def apply_inverse_adjoint(self, U, mu=None, least_squares=False):
+    def apply_inverse_adjoint(self, U, mu=None, initial_guess=None, least_squares=False):
         assert U in self.source
+        assert initial_guess is None or initial_guess in self.range and len(initial_guess) == len(U)
         KitU0 = self.K.apply_inverse_adjoint(U.block(0), mu=mu, least_squares=least_squares)
         V_blocks = [-self.E.apply_adjoint(KitU0, mu=mu) + U.block(1),
                     -KitU0]
@@ -381,8 +391,9 @@ def apply_adjoint(self, V, mu=None):
                     - self.E.apply_adjoint(V.block(1), mu=mu) * self.b.conjugate()]
         return self.source.make_array(U_blocks)
 
-    def apply_inverse(self, V, mu=None, least_squares=False):
+    def apply_inverse(self, V, mu=None, initial_guess=None, least_squares=False):
         assert V in self.range
+        assert initial_guess is None or initial_guess in self.source and len(initial_guess) == len(V)
         aMmbEV0 = self.M.apply(V.block(0), mu=mu) * self.a - self.E.apply(V.block(0), mu=mu) * self.b
         KV0 = self.K.apply(V.block(0), mu=mu)
         a2MmabEpb2K = (self.a**2 * self.M - self.a * self.b * self.E + self.b**2 * self.K).assemble(mu=mu)
@@ -393,8 +404,9 @@ def apply_inverse(self, V, mu=None, least_squares=False):
                     + a2MmabEpb2KiV1 * self.a]
         return self.source.make_array(U_blocks)
 
-    def apply_inverse_adjoint(self, U, mu=None, least_squares=False):
+    def apply_inverse_adjoint(self, U, mu=None, initial_guess=None, least_squares=False):
         assert U in self.source
+        assert initial_guess is None or initial_guess in self.range and len(initial_guess) == len(U)
         a2MmabEpb2K = (self.a**2 * self.M - self.a * self.b * self.E + self.b**2 * self.K).assemble(mu=mu)
         a2MmabEpb2KitU0 = a2MmabEpb2K.apply_inverse_adjoint(U.block(0), mu=mu, least_squares=least_squares)
         a2MmabEpb2KitU1 = a2MmabEpb2K.apply_inverse_adjoint(U.block(1), mu=mu, least_squares=least_squares)