Updated explanation from PR

examples/ex16p.cpp

@@ -65,7 +65,7 @@ protected:
    HypreSmoother M_prec; // Preconditioner for the mass matrix M
 
    CGSolver T_solver;    // Implicit solver for T = M + dt K
-   HypreSmoother T_prec; // Preconditioner for the implicit solver
+   HypreBoomerAMG T_prec; // Preconditioner for the implicit solver
 
    double alpha, kappa;
 
@@ -353,7 +353,6 @@ int main(int argc, char *argv[])
          }
 #endif
       }
-      oper.SetParameters(u);
    }
 
 #ifdef MFEM_USE_ADIOS2
@@ -414,6 +413,7 @@ ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
    T_solver.SetPrintLevel(0);
    T_solver.SetPreconditioner(T_prec);
 
+   T_prec.SetPrintLevel(0);
    SetParameters(u);
 }
 
@@ -430,19 +430,53 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
 void ConductionOperator::ImplicitSolve(const double dt,
                                        const Vector &u, Vector &du_dt)
 {
-   // Solve the equation:
-   //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
-   // for du_dt
-   if (!T)
-   {
-      T = Add(1.0, Mmat, dt, Kmat);
-      current_dt = dt;
+   // Here we use Picard iterations to solve a nonlinear equation
+   // for the Runge-Kutta stage vector k,
+   //
+   //    M*k = N(u+dt*k)                     (1)
+   //
+   // for nonlinear operator N. We assume N can be written as
+   //
+   //    N(u+dt*k) := L[u+dt*k](u+dt*k) + f(t)
+   //
+   // where L is a matrix-valued operator evaluated at u+dt*k and f(t)
+   // a (potentially zero) time-dependent forcing vector. (1) can be
+   // rewritten as a fixed-point equation
+   //
+   //    x = (M - dt*L[x])^{-1} (Mu + f)     (2)
+   //
+   // where x := u + dt*k, which can be solved using a Picard iteration,
+   // where a function G(x) = x is solved via iteraitons x_{k+1} = G(x_k).
+
+   double tol = 1e-6;
+   int maxiter = 100;
+
+   // Right-hand side for nonlinear iteration
+   Mmat.Mult(u, z);  // Add forcing function if one exists
+   du_dt = u;        // Set u as initial guess for x (2)
+   Vector temp(u);   // Vector to measure error
+   temp = u;
+   double error = 1;
+   int iter = 0;
+   while (error > tol) {
+      iter ++;
+      this->SetParameters(du_dt);   // Update nonlinear operator L[x]
+      T = Add(1.0, Mmat, dt, Kmat); // Form matrix (M - dt*L[x])
       T_solver.SetOperator(*T);
+      T_solver.Mult(z, du_dt);      // Apply (M - dt*L[x])^{-1}
+      temp -= du_dt;                // Measure error
+      error = std::sqrt(InnerProduct(MPI_COMM_WORLD, temp, temp));
+      temp = du_dt;
+      if (iter >= maxiter) {
+         mfem_warning("Nonlinear iteration did not converge!");
+         break;
+      }
    }
-   MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
-   Kmat.Mult(u, z);
-   z.Neg();
-   T_solver.Mult(z, du_dt);
+
+   // Above we solved for x = u + dt*k, where k is the desired update
+   // Map du_dt -> k.
+   du_dt -= u;
+   du_dt /= dt;
 }
 
 void ConductionOperator::SetParameters(const Vector &u)
@@ -483,4 +517,4 @@ double InitialTemperature(const Vector &x)
    {
       return 1.0;
    }
-}
+}

examples/ex1p.cpp

@@ -231,10 +231,17 @@ int main(int argc, char *argv[])
    else
    {
       prec = new HypreBoomerAMG;
+      prec->SetOperator(*A);
    }
-   CGSolver cg(MPI_COMM_WORLD);
+   // CGSolver cg(MPI_COMM_WORLD);
+   AndersonAcceleration cg(MPI_COMM_WORLD);
+   cg.SetKDim(10);
+   cg.SetRestart(true);      // WORKS
+   cg.SetAAStart(0);          // WORKS
+   cg.SetWeight(1.0);         // Not robust but seems to work
+
    cg.SetRelTol(1e-12);
-   cg.SetMaxIter(2000);
+   cg.SetMaxIter(50);
    cg.SetPrintLevel(1);
    if (prec) { cg.SetPreconditioner(*prec); }
    cg.SetOperator(*A);

linalg/operator.hpp

@@ -651,7 +651,6 @@ public:
    virtual void SetOperator(const Operator &op) = 0;
 };
 
-
 /// Identity Operator I: x -> x.
 class IdentityOperator : public Operator
 {

linalg/solvers.cpp

@@ -1733,6 +1733,323 @@ void NewtonSolver::AdaptiveLinRtolPostSolve(const Vector &x,
    }
 }
 
+void AndersonAcceleration::QRdelete(std::deque<Vector *> &Q, DenseMatrix &R) const
+{
+   Vector temp(Q[0]->Size());
+   for (int i=0; i<(maxVecs-1); i++) {
+      double d = sqrt( R(i,i+1)*R(i,i+1) + R(i+1,i+1)*R(i+1,i+1) );
+      double c = R(i,i+1) / d;
+      double s = R(i+1,i+1) / d;
+      R(i,i+1) = d;
+      R(i+1,i+1) = 0;
+
+      if (i < (maxVecs-2)) {
+         for (int j=(i+2); j<maxVecs; j++) {
+            d = c*R(i,j) + s*R(i+1,j);
+            R(i+1,j) = -s*R(i,j) + c*R(i+1,j);
+            R(i,j) = d;
+         }
+      }
+      // temp = c*Q[i] + s*Q[i+1];
+      add(c, *(Q[i]), s, *(Q[i+1]), temp);
+      // Q[i+1] = -s*Q[i] + c*Q[i+1];
+      *(Q[i+1]) *= c;
+      Q[i+1] -> Add(-s, *(Q[i]));
+      *(Q[i]) = temp;
+   }
+   
+   // Shift Q <- Q[:,0:(m-2)], i.e., delete last column of Q
+   delete Q.back();
+   Q.pop_back();
+
+   // Shift columns of R to the left by one, R = R[0:(m−2), 1:(m-1)]
+   for (int j=1; j<maxVecs; j++) {
+      for (int i=0; i<maxVecs; i++) {
+         R(i,j-1) = R(i,j);
+      }
+   }
+}
+
+void AndersonAcceleration::SetOperator(const Operator &op)
+{
+   oper = &op;
+   height = op.Height();
+   width = op.Width();
+   MFEM_ASSERT(height == width, "square Operator is required.");
+}
+
+void AndersonAcceleration::FixedPointMult(const Vector &b,
+   const Vector &x, Vector &y, Vector &r) const
+{
+   // Assume Operator represents a fixed-point operator *with
+   // right-hand side*, so we iterate x_{k+1} = M^{-1}G(x_k)
+   if (isFixedPointOp) {
+      oper->Mult(x, y);
+      if (prec) {
+         prec->Mult(y,r);
+         y = r;
+         r -= x;
+      }
+      else {
+         r = y;
+         r -= x;
+      }
+   }
+   // Otherwise add x to y and not to residual r = y - x
+   else {
+      oper->Mult(x, r);
+      r *= -1;
+      r += b;
+      if (prec) {
+         prec->Mult(r,y);
+         r = y;
+         y += x;
+      }
+      else {
+         y = r;
+         y += x;         
+      }
+   }
+}
+
+void AndersonAcceleration::Mult(const Vector &b, Vector &x) const
+{
+   MFEM_ASSERT(oper != NULL, "the Operator is not set (use SetOperator).");
+  
+   int n = width;
+   int numVecs = 0;
+   double resid, norm_df;
+   double min_diag = 1e-13;
+   final_norm = -1;
+
+   // Check vector is initialized, set to zero for
+   // iterative_mode = false
+   if (x.Size() != n) {
+      x.SetSize(n);
+      x = 0.0;
+   }
+   else if (!iterative_mode) {
+      x = 0.0;
+   }
+
+   // Storage containers for acceleration
+   std::deque<Vector *> G;
+   std::deque<Vector *> Q;
+   DenseMatrix R(maxVecs);
+   R = 0.0;
+   Vector g_old(n);
+   Vector g_current(n);
+   Vector f_old(n);
+   Vector f_current(n);
+   Vector gamma(maxVecs);
+   Vector rhs(maxVecs);
+   Vector correction;
+   if (omega > 0 && std::abs(omega -  1) > 1e-14) {
+      correction.SetSize(n);
+   }
+   Vector *dg;
+   Vector *df;
+
+   // Loop over AA iterations
+   int k;
+   for (k=0; k<max_iter; k++) {
+
+      // Compute g_current = G(x), f_current = G(x) - x
+      this->FixedPointMult(b, x, g_current, f_current); 
+
+      // Check norm of current approximation to fixed point G(u) = u
+      resid = Norm(f_current);
+      MFEM_ASSERT(IsFinite(resid), "||G(u) - u|| = " << resid);
+      if (print_level == 1)
+      {
+         mfem::out << "  Iteration : " << setw(3) << k
+                 << "  ||G(u) - u|| = " << resid << endl;
+      }
+      
+      // Set stopping tolerance on first iteration.
+      if (final_norm < 0) {
+         final_norm = std::max(rel_tol*resid, abs_tol);
+      }
+
+      // Check for convergence
+      if (resid <= final_norm)
+      {
+         final_norm = resid;
+         final_iter = k;
+         converged = 1;
+         goto finish;
+      }
+
+      // Start Anderson Acceleration after AAstart FP iterations
+      if (k > AAstart) {
+         // df = f_current - f_old; 
+         df = new Vector(n);
+         add(1.0, f_current, -1.0, f_old, *df);
+         // dg = g_current - g_old;
+         dg = new Vector(n);
+         add(1.0, g_current, -1.0, g_old, *dg);
+      
+         if (numVecs < maxVecs) {
+            G.push_back(dg);
+         }
+         else {
+            delete G[0];
+            G.pop_front();
+            G.push_back(dg);
+         }
+         numVecs++;
+         dg = NULL;
+      }
+      
+      f_old = f_current;
+      g_old = g_current;
+      
+      // First iteration or initial fixed-point iterations
+      if (numVecs == 0) {
+         x = g_current;
+         continue;
+      }
+
+      // All later iterations: orthogonalize and find best approximation
+      if (numVecs == 1) {
+         norm_df = Norm(*df);
+         MFEM_ASSERT(IsFinite(norm_df), "norm_df = " << norm_df);
+         (*df) /= norm_df;
+         Q.push_back(df);
+         R(0,0) = norm_df;
+         df = NULL;
+      }
+      else {
+         // Remove first column in basis F and R, reorthogonalize
+         if (numVecs > maxVecs) {
+            this->QRdelete(Q, R);
+            numVecs--;
+         }
+         // Compute last column of R
+         for (int i=0; i<(numVecs-1); i++) {
+            R(i,numVecs-1) = Dot(*(Q[i]), *df);
+            // df -= R(i,numVecs-1) * Q[i]
+            df -> Add(-R(i,numVecs-1), *(Q[i]));
+         }
+         norm_df = Norm(*df);
+         MFEM_ASSERT(IsFinite(norm_df), "norm_df = " << norm_df);
+         (*df) /= norm_df;
+         Q.push_back(df);
+         R(numVecs-1, numVecs-1) = norm_df;
+         df = NULL;
+      }
+
+      // Back solve for new weights, R\gamma = Q^T * f_current
+      rhs = 0.0;
+      gamma = 0.0;
+      for (int i=0; i<numVecs; i++) {
+         rhs(i) = Dot(*(Q[i]), f_current);   // Form right hand side
+      }
+      for (int i=(numVecs-1); i>=0; i--) {
+         double temp = rhs(i);
+         for (int j=(i+1); j<numVecs; j++) {
+            temp -= R(i,j)*gamma(j);
+         }
+         if (std::abs(R(i,i)) < min_diag) {
+            gamma(i) = 0.0;
+            std::cout << "Diagonal of R -- " << R(i,i) << " ~ 0.\n";
+         }
+         else {
+            gamma(i) = temp / R(i,i);            
+         }
+      }
+
+      /// DEBUG --> test backsolve
+      Vector test(numVecs);
+      for (int i=0; i<numVecs; i++) {
+         test(i) = 0;
+         for (int j=i; j<numVecs; j++) {
+            test(i) += R(i,j) * gamma(j);
+         }
+         if (std::abs(rhs(i) - test(i)) > 1e-10) {
+            std::cout << "Bad solve! Err = " << rhs(i) - test(i) << "\n";
+         }
+      }
+
+      // Compute updated solution x = g_current − G*\gamma 
+      x = g_current;
+      for (int i=0; i<numVecs; i++) {
+         // x -= gamma(i)*G[i]
+         x.Add(-gamma(i), *(G[i]));
+      }
+
+      // Apply damped iteration for \omega \in (0,1),
+      //   x -= (1−omega) * (f_current − Q*R*gamma);
+      if (omega > 0 && std::abs(omega -  1) > 1e-14) {
+         // Redefine rhs = R*gamma
+         for(int i=0; i<numVecs; i++) {
+            rhs(i) = 0;
+            for (int j=i; j<numVecs; j++) {
+               rhs(i) += R(i,j)*gamma(j);
+            }
+         }
+         correction = f_current;
+         for (int i=0; i<numVecs; i++) {
+            // correction -= rhs(i)*Q[i]
+            correction.Add(-rhs(i), *(Q[i]));
+         }
+         // x -= (1 - omega) * correction;
+         x.Add( -(1 - omega), correction);
+      }
+
+      // Restart AA minimization by eliminating all vectors but the most recent
+      if (restart && (numVecs == maxVecs)) {
+         for (int i=0; i<(maxVecs-1); i++) {
+            delete G[0];
+            G.pop_front();          
+            delete Q[0];
+            Q.pop_front();
+         }
+         R = 0.0;
+         R(0,0) = norm_df;
+         numVecs = 1;
+         if (print_level == 1)
+         {
+            mfem::out << "Restarting..." << '\n';
+         }
+      }
+   }
+
+   // Compute final residual, save counts for solve
+   this->FixedPointMult(b, x, g_current, f_current); 
+   resid = Norm(f_current);
+   MFEM_ASSERT(IsFinite(resid), "||G(u) - u|| = " << resid);
+   final_norm = resid;
+   final_iter = max_iter;
+   if (resid <= final_norm) converged = 1;
+   else converged = 0;
+
+finish:
+   if (print_level == 3)
+   {
+      mfem::out << "  Iteration : " << setw(3) << k
+              << "  ||G(u) - u|| = " << resid << endl;
+   }
+   else if (print_level == 2)
+   {
+      mfem::out << "Anderson Acceleration: Number of iterations: " << final_iter << '\n';
+   }
+   if (print_level >= 0 && !converged)
+   {
+      mfem::out << "Anderson Acceleration: No convergence!\n";
+   }
+
+   // Cleanup pointers
+   for (int i=0; i<numVecs; i++) {
+      delete G[0];
+      G.pop_front();          
+      delete Q[0];
+      Q.pop_front();
+   }
+   delete dg;
+   delete df;
+}
+
 void LBFGSSolver::Mult(const Vector &b, Vector &x) const
 {
    MFEM_VERIFY(oper != NULL, "the Operator is not set (use SetOperator).");

linalg/solvers.hpp

@@ -14,6 +14,7 @@
 
 #include "../config/config.hpp"
 #include "densemat.hpp"
+#include <deque>
 
 #ifdef MFEM_USE_MPI
 #include <mpi.h>
@@ -489,6 +490,64 @@ public:
                            const double gamma = 1.0);
 };
 
+
+/// Nonlinear Anderson Acceleration
+class AndersonAcceleration : public IterativeSolver
+{
+protected:
+   int maxVecs; // see SetKDim()
+   int AAstart;
+   bool restart;
+   bool isFixedPointOp;
+   double omega;
+
+   /// Apply fixed-point Mult and compute residual.
+   //  For !isFixedPointOp:
+   //    y <-- G(x) = x + A(x) - b and r <-- G(x) - x = A(x) - b
+   //  For isFixedPointOp:
+   //    y <-- G(x) = A(x) - b and r <-- G(x) - x = A(x) - b - x
+   void FixedPointMult(const Vector &b, const Vector &x, Vector &y, Vector &r) const;
+
+   // Helper function for Anderson Acceleration
+   void QRdelete(std::deque<Vector *> &Q, DenseMatrix &R) const;
+
+public:
+   AndersonAcceleration() : maxVecs(25), AAstart(0), omega(1),
+      restart(false), isFixedPointOp(false) { }
+
+#ifdef MFEM_USE_MPI
+   AndersonAcceleration(MPI_Comm _comm) : IterativeSolver(_comm),
+      maxVecs(25), AAstart(0), omega(1), restart(false),
+      isFixedPointOp(false) { }
+#endif
+
+   /// Boolean describing whether the action of the operator is such
+   //  that we want to solve G(x) = x (true) or G(x) = 0 (false) for
+   //  zero right-hand side vector b passed into Mult(). Default is
+   //  false in construction of class.
+   void IsFixedPointOperator(bool isFixedPointOp_)
+   { isFixedPointOp = isFixedPointOp_; }
+
+   /// Maximum number of vectors to store in Krylov-like space
+   void SetKDim(int dim) { maxVecs = dim; }
+
+   /// Number of fixed-point iterations to do before starting AA
+   void SetAAStart(int start_) { AAstart = start_; }
+
+   /// Boolean to restart, that is, erase entire space after maxVecs
+   //  are stored (AAstart=true) or use a sliding space (AAstart=false)
+   //  where one vector is deleted to make room for a new one.
+   void SetRestart(bool restart_) { restart = restart_; }
+
+   /// Set relaxation weight
+   void SetWeight(double omega_) { omega = omega_; }
+
+   virtual void SetOperator(const Operator &op);
+
+   virtual void Mult(const Vector &b, Vector &x) const;
+};
+
+
 /** L-BFGS method for solving F(x)=b for a given operator F, by minimizing
     the norm of F(x) - b. Requires only the action of the operator F. */
 class LBFGSSolver : public NewtonSolver