\section{Introduction}
\label{sec:introduction}

The international community is committed to the development of fusion energy in order to
meet growing world energy needs.  This is a complex, international, and costly endeavor;
the next experimental magnetic-fusion plasma confinement device,
ITER~\cite{iter:homepage}, is estimated to cost in excess of \$10B.  The handling of the
plasma exhaust is a challenge for the operation of ITER and for the design of the future
demonstration power plant DEMO. To achieve its design goals, ITER must operate in a
high-performance regime with steep pressure gradients close to the plasma edge
\cite{keilhacker1984,doyle2008}. While this capability improves the plasma performance,
the steep pressure gradients and associated bootstrap current can destabilize
peeling-ballooning modes \cite{connor-1998,snyder-2002,wilson-2002}. These result in
eruptions from the plasma edge, known as edge localized modes (ELMs) (see, e.g.,
\cite{maggi2010}). Scaling from existing experiments indicates that these violent events,
if uncontrolled, could produce peak power loads of around 1 GW/m$^2$ onto material
surfaces in ITER \cite{leonard1999,federici2003} and potentially more in a DEMO power
plant.  Such loads would result in unacceptable erosion of plasma facing components. Thus,
researchers are extremely interested in understanding and controlling ELMs, and work is
ongoing to use several nonlinear simulation codes in order to simulate ELMs, predict their
size in ITER, and find ways to mitigate them.

Nonlinear simulations of tokamak plasmas are challenging because of the wide range of
length scales and timescales present: The characteristic ion cyclotron motion has scales
around $1$ mm and $10^{-7}$ seconds, while the machine has length scales of $\mathord{\sim} 10$ m
and transport timescales in seconds. If electron physics is important, then the
characteristic scales become smaller by the square root of the mass ratio, a factor of
$\mathord{\sim} 60$.  These small cyclotron scales are due to the strong magnetic field (typically $1
- 5$ Teslas), which is used to confine the plasma in a tokamak. This produces strong
anisotropy in the system, as particles are confined perpendicular to the magnetic field
but are essentially free to move along it, resulting in a difference of $10^8$ in heat
conduction rates perpendicular and parallel to the field. The topology of the magnetic
field poses challenges for simulation, as the core fieldlines form helices with a pitch
angle that varies across the device. At the boundary of the core plasma, in the region we
are interested in simulating, a separatrix is formed, outside of which the fieldlines
intersect material boundaries rather than forming closed toroidal ``flux'' surfaces.
Furthermore, the edge region of tokamaks is complicated by the interaction with neutral
gas and material surfaces, bringing in atomic physics, sputtering, and deposition
processes.

Analytic reduction of the full problem, employing asymptotic expansions in small
parameters such as the ratio of the Larmor radius to machine size, has been successfully
used for many years (see, e.g., \cite{krommes2012} and references therein). In this paper
we consider fluid equations, which evolve moments of the particle velocity distribution,
and remove the smallest scales associated with cyclotron motion \cite{braginskii1965}.
Further simplifications include reduced magnetohydrodynamics (MHD), exploiting the
anisotropy of the system in which length scales parallel to the magnetic field are much
larger than those perpendicular to the field (see, e.g., \cite{strauss1976,
  hazeltinemeiss}).  By assuming force balance, this approach allows us to analytically
remove the fast magnetosonic wave, which would be the fastest wave in our system. This
wave is not important for the relatively slow MHD instabilities we wish to study here, but
it is responsible for establishing force balance along the magnetic field. After this
reduction, the system of equations still contains two fast timescales that cannot be
removed without also removing physics of interest: the shear Alfv\'en wave and the
parallel heat transport. These can severely limit the timestep in numerical simulations,
so implicit and semi-implicit methods are widely used in plasma simulation
\cite{jardin2012}.  In BOUT++ edge simulations of ELMs \cite{dudson2009bout,xu2010prl,
  dudson2011ppcf,xu2011nf}, such techniques help address challenges in the broad range of
temporal scales, complicated geometry with a range of spatial scales, and competing
physical effects of the various degrees of freedom.

This paper focuses on efficient and scalable nonlinear solvers for BOUT++ ELM simulations.
We construct a \emph{physics-based} preconditioner that reduces the time to solve
nonlinear systems with Jacobian-free Newton-Krylov (JFNK) techniques~\cite{knollkeyes:04}.  The
term \emph{physics-based preconditioner} has been used to describe various approaches to
application-specific preconditioners that incorporate insight into the physics behind the
equations being solved \cite{MousseauKnoll03,ChaconKnollEtAl02,SantilliScotti11,
  PremnathPattisonEtAl09,HirshmanSanchezEtAl11,jamroz2010,McCourtRognlienEtAl12}.
Inspired by the approach of Chac\'on~\cite{chacon2008}, we use the term here to describe
our development of a custom preconditioning algorithm based on analytic reduction and Schur
factorization of the system Jacobian. Identifying the key physics terms allows a
simplification of the resulting equations and enables efficient numerical solution.  A
similar approach has been previously attempted within the NIMROD code \cite{jamroz2010},
where a Crank-Nicholson (CN) Jacobian-free Newton-Krylov scheme was used along with
physics-based preconditioning. In that work numerical instability was encountered and
attributed to the CN scheme. In this paper we exploit the field-aligned coordinate system
used in BOUT++~\cite{dudson2009bout} to construct a preconditioner that reduces a large 3D
problem to the solution of many one-dimensional systems, each of which can be solved
independently. This approach greatly improves convergence of Newton-Krylov methods in
implicit timestepping schemes, without causing numerical instability.

We also introduce the complementary approach of nonlinear GMRES \cite{Oosterlee:2000tg},
which offers the advantage of requiring only nonlinear function evaluations from the
application, rather than also needing Jacobian-vector products and preconditioning for
Newton-based methods.  We demonstrate that both of these approaches---physics-based
preconditioned Newton-Krylov and nonlinear GMRES---significantly reduce the time for
solving the nonlinear ELM systems at each timestep compared with a Newton-Krylov approach
with no preconditioning, and we demonstrate good strong scaling on up to 512 processor
cores.

The remainder of this paper explains the details of our approach.  Section
\ref{sec:BOUT++} introduces the BOUT++ application and its approach to modeling ELMs.
Section \ref{sec:nonlinear_solvers} describes the algorithms and software that we use to
solve the resulting nonlinear systems that arise at each timestep, including both
Newton-Krylov methods with physics-based preconditioning and nonlinear GMRES.  Section
\ref{sec:results} presents experimental results, while Section \ref{sec:conclusions}
discusses conclusions and directions for future work.

\section{BOUT++ Simulation Overview}
\label{sec:BOUT++}

The numerical methods studied in this paper are quite general, but are applied here to the
solution of reduced MHD for fusion applications within the BOUT++ framework.

\subsection{\texttt{BOUT++ } Capabilities}
\label{subsec:boutsimcode}

BOUT++, a parallel MPI + OpenMP finite
difference code written in C++ \cite{dudson2009bout}, simulates a range of plasma
fluid equations in curvilinear coordinate systems---in particular, reduced MHD models in
coordinates aligned with a strong background magnetic field. Current applications of
BOUT++ include the study of ELMs \cite{xu2010prl,dudson2011ppcf,xu2011nf}, plasma
turbulence \cite{friedman2012arxiv}, and blob dynamics \cite{angus2012prl}. In this paper
BOUT++ is extended to more efficiently simulate heat conduction parallel to magnetic
fields, an important process in all these areas. Parallel heat conduction modifies the
linear structure and growth rate of plasma instabilities and plays a crucial role in the
loss of energy to material surfaces in tokamak edge simulations of turbulence and ELMs.

\subsection{Reduced MHD Equations}
\label{subsec:equations}

Since the linear stability threshold is well described by ideal MHD~\cite{snyder-2002}, it
is reasonable to start studying ELMs by using fluid models based on MHD. Several initial-value codes
including BOUT++~\cite{dudson2011ppcf}, M3D-C$^1$~\cite{ferraro2010} and NIMROD~\cite{burke2010},
have now been benchmarked and shown good agreement with linear ideal MHD codes ELITE~\cite{snyder-2002,wilson-2002} and
GATO~\cite{bernard1981}.
Nonlinearly, ideal MHD theory predicts
narrowing structures erupting outwards, leading to a finite time singularity
\cite{wilson-2004}.  Long before this happens, additional effects must become important in
order to allow plasma to decouple from magnetic fields and relax to a new
equilibrium. Nonlinear modeling of ELMs has therefore concentrated on fluid models
incorporating a range of different nonideal effects \cite{burke2010,pamela2010,
  hoelzl2012,sugiyama2009,jardin2007,ferraro2010}.

Since each simulation code currently uses different fluid models, nonlinear benchmarking
of codes is difficult. In order to enable future comparisons between codes, a set of
equations close to those used in JOREK~\cite{huysmans2009} has been implemented in BOUT++,
which we will refer to in this paper as the six-field model.  Six scalar fields are
evolved: mass density $\rho$; electron and ion temperatures $T_e$ and $T_i$, respectively;
vorticity $U=\vec{b}_0\cdot\nabla\times\vec{v} \simeq \frac{1}{B_0}\nabla_\perp^2\phi$;
parallel velocity $v_{||} = \vec{b}_0\cdot\vec{v}$; and parallel vector potential $A_{||} =
\vec{b}_0\cdot\vec{A}$. From these auxiliary quantities are calculated the total plasma
pressure $P = \rho (T_e + T_i)$ and parallel current density $J_{||} =
\vec{b}_0\cdot\vec{J} = -\nabla_\perp^2 A_{||}$.  The equilibrium magnetic field unit
vector is $\vec{b}_0=\vec{B}_0/B_0$. Here all quantities with subscript ``0'' are
calculated from the starting equilibrium and are not evolved. The evolving quantities are
deviations from equilibrium; for example, the total mass density is $\rho_0 + \rho$. The
equations for plasma density and temperatures are
\begin{align}
  \pd{\rho}{t} &= -\vec{v} \cdot \nabla \left(\rho + \rho_0\right)
                + \left(\nabla \cdot \vec{v}\right) \left(\rho + \rho_0\right)
                + D_\perp\nabla_\perp^2\rho \cr
  \pd{T_s}{t} &= -\vec{v} \cdot \nabla \left(T_s + T_{s0}\right)
               - \frac{2}{3} \left(\nabla \cdot \vec{v}\right) \left(T_s + T_{s0}\right) \cr
              &+ \frac{1}{\rho + \rho_0} \left[
                   \nabla_{||} \cdot \left(\chi_{s||}\partial_{||}T_s\right)
                   + \chi_{S\perp}\nabla_\perp^2T_s
                 \right] + \frac{2}{3\left(\rho+\rho_0\right)} W_s ,
  \label{eq:6field1}
\end{align}
where the total plasma velocity $\vec{v} = \vec{v}_{E\times B} + \vec{b} v_{||}$ is given
by an $E\times B$ drift perpendicular to magnetic fieldlines $\vec{v}_{E\times B} =
\frac{1}{B_0}\vec{b}_0\times\nabla\phi$ and parallel flow $v_{||}$ along them. The above
equations describe the advection ($\vec{v}\cdot \nabla$) and compression ($\nabla \cdot
\vec{v}$) of density and temperatures ($s = {e,i}$), collisional diffusion perpendicular
to the magnetic field with coefficients $D_\perp$ and $\chi_{s\perp}$, and parallel
conduction of heat along magnetic fields with coefficient $\chi_{s||}$. Derivatives are
split into those parallel to the magnetic field $\nabla_{||} = \vec{b}_0\cdot\nabla$ and
those perpendicular $\nabla_\perp = \nabla - \vec{b}_0\cdot\nabla$.

The velocity of the plasma is described by the vorticity and parallel flow:
\begin{align}
  \pd{U}{t} &= -\vec{v}\cdot\nabla U + \frac{1}{\rho+\rho_0} \left[B_0^2\nabla_{||}
                \left(\frac{J_{||} + J_{||0}}{B_0}\right) \right. \cr
            &+ \left. 2\vec{b}_0 \times \vec{k}_0 \cdot \nabla P
             + \nu_{||}\partial_{||}^2U + \nu_\perp\nabla_\perp^2U \right] \cr
  \pd{v_{||}}{t} &= -\vec{v}\cdot\nabla v_{||} - \frac{1}{\rho+\rho_0}\nabla_{||} P .
  \label{eq:6field2}
\end{align}
The vorticity equation (evolving $U$) contains a term responsible for shear Alfv\'en wave
propagation ($B_0^2\nabla_{||}\frac{J_{||}}{B_0}$), instability drive due to equilibrium
magnetic field curvature $\kappa$ and pressure gradient $\nabla P$, as well as parallel
and perpendicular viscosity terms $\nu_{||}$ and $\nu_\perp$.

The evolution of the magnetic field is described by Ohm's law for the electric field
parallel to the magnetic field $\vec{b}\cdot\vec{E} = \vec{b}\cdot\eta\vec{J}$,
\begin{equation}
  \pd{A_{||}}{t} = -\nabla_{||}\phi - \eta J_{||} ,
  \label{eq:6field3}
\end{equation}
where $\eta$ is the parallel resistivity.

\subsection{Timescales for Shear Alfv\'en Wave and Parallel Heat Conduction}
\label{subsec:alfvenwave}

The timescales of interest in equations~\eqref{eq:6field1}-\eqref{eq:6field3} are those of
the linear growth of peeling-ballooning modes and nonlinear eruption of filaments,
typically of the order of $\mathord{\sim} 100\mu$s up to milliseconds. The above set of equations
contains timescales much faster than this, which must be preconditioned in an implicit
scheme.  The vorticity and $A_{||}$ equations contain the terms
\[
  \pd{U}{t} =  \frac{B_0^2}{\rho_0}\nabla_{||} \left(\frac{J_{||}}{B_0}\right)
  \qquad
  \pd{A_{||}}{t} = -\nabla_{||}\phi ,
\]
which reduce to
\[
  \pdd{\phi}{t} = \nabla_\perp^{-2} \left[
                    \frac{B_0^3}{\rho_0}\nabla_{||} \left(\nabla_\perp^2\nabla_{||}\phi / B_0\right)
                  \right]
           \simeq \frac{B_0^2}{\rho_0}\nabla_{||}^2 \phi ,
\]
corresponding to the shear Alfv\'en wave traveling along magnetic fieldlines. Given the
strong magnetic field $B_0$ and low mass density $\rho_0$ in tokamaks, wave speeds of $v_A
\mathord{\sim} 10^7$m/s are typical, which for parallel grid spacing of $\mathord{\sim} 10$ cm $ - 1$ m give
CFL timesteps $\mathord{\sim} 10^{-8}$s, 3 to 4 orders of magnitude shorter than timescales of
interest.

The second fast timescale is introduced through the ion and electron temperature equation
\[
  \pd{T_s}{t} = \frac{1}{\rho_0}\nabla_{||}\left(\chi_{s||}\partial_{||}T_s\right) ,
\]
which is a diffusion equation along magnetic fieldlines. At the high temperatures relevant
to fusion, electron thermal speeds along fieldlines of $10^6-10^7$m/s are typical, and
hence parallel thermal diffusion is also fast compared with nonlinear evolution
timescales.

The above equations contain one other wave along magnetic fieldlines: the sound wave
described by
\[
  \pd{v_{||}}{t} = -\frac{1}{\rho_0}\nabla_{||} P \qquad
  \pd{\rho}{t} = -\rho_0\nabla_{||}v_{||} \qquad
  \pd{T_s}{t} = -\frac{2}{3}T_{s0}\nabla_{||}v_{||} .
\]
The speed of this wave is typically around $10^5$ m/s for the problems in which we are
interested, several orders of magnitude slower than the Alfv\'en wave and parallel
electron heat conduction. It is trivial to extend the preconditioner described below to
include this wave, but this extension was found to make little difference to the
convergence rate and so is not included here.

In the following sections we outline the implicit timestepping method and construct a
preconditioning scheme that can be applied to address both the Alfv\'en wave and
temperature conduction timescales.

\section{Scalable Nonlinear Solvers}
\label{sec:nonlinear_solvers}

The coupled set of reduced MHD equations~\eqref{eq:6field1}-\eqref{eq:6field3} are
discretized by the method of lines, using 4th-order central differencing for diffusive terms,
and 3rd order WENO \cite{jiang-1996} for upwinding. Because this work focuses on nonlinear
solvers, we chose the timestepping algorithm to be the same for each test case, namely, a
variable order backward differential formulation (BDF)~(see, e.g., \cite{Ascher_B1998}).
This family comprises methods that are implicit and linear multistep, and that allow
variable timestepping.  At each timestep, by writing the system state as a vector, we
obtain from equations~\eqref{eq:6field1}-\eqref{eq:6field3} a nonlinear system of the form
\begin{equation}
f(u) = f \left( \begin{array}{c}
  \rho \\
  T_i \\
  T_e \\
  U \\
  v_{||} \\
  A_{||}\end{array}\right) = 0 ,
\label{eq:f=0}
\end{equation}
where $f: \, \bb{R}^n \to \bb{R}^n$.
This requires the solution of a nonlinear system, here using Newton-Krylov and NGMRES
methods.

\subsection{Newton-Krylov Methods}
\label{subsec:newtonprecon}

In a Newton-Krylov method, we solve the nonlinear system~\eqref{eq:f=0} at each timestep
using a Newton iteration (see, e.g., \cite{NW99}),
\[
  u_{k+1} = u_k - [ f'(u_k)]^{-1} f(u_k), \;\; k=0,1, \ldots.
\]
Here $u_0$ is an initial approximation to the solution, and $f'(u_k)$, the Jacobian, is
nonsingular at each iteration. A Krylov iterative method (here GMRES) is then used to
solve this linear system at each Newton iteration.  In practice, the Newton iteration is
implemented by the following two steps:
\begin{enumerate}
  \item (Approximately) solve \;\; $f'(u_k) \Delta u_k = -f(u_k)$.
  \item Update \;\; $u_{k+1} = u_k + \alpha \Delta u_k$.
\end{enumerate}
Here $0 < \alpha \leq 1$ is a scalar. In these experiments, we use a JFNK variant, where
Jacobian-vector products are computed matrix-free~\cite{knollkeyes:04} by an approximation
\[
  f'(v) a \approx \frac{f(v + h \cdot a) - f(v)}{h},
\]
for a differencing parameter $h$.  By computing finite-difference Jacobian-vector products
for the iterative solve, assembling the actual Jacobian is never required.

If the nonlinear system being evolved contains a wide range of timescales, as in most
plasma simulations, then it is said to be ``stiff,'' and the Jacobian can be
ill-conditioned. In this case the Krylov method will require an unacceptably large number
of iterations.  Preconditioning the linearized Newton systems, which have the form
$\mathbb{J} x = b$, where $\mathbb{J}$ is the nonsingular $n \times n$ Jacobian matrix and
$x$ and $b$ are $n$-dimensional vectors, transforms these to equivalent systems
$\mathbb{P} \mathbb{J} x = \mathbb{P} b$ through the action of a preconditioner,
$\mathbb{P}$, whose inverse action approximates that of the Jacobian but at smaller cost.
Choosing an effective approximation is one of the goals of this work, since the
preconditioning phase is crucial to achieving low computational cost and scalable
parallelism.

\subsection{Physics-Based Preconditioning}
\label{subsec:preconditioning}

Physics-based preconditioning \cite{chacon2008,chacon02,ChaconKnollEtAl02} is a method by
which insight into the physics behind the equations being solved is used to precondition
the Newton-Krylov method.  After first identifying the stiff terms in the equations, an
approximate inverse is constructed with block Gaussian elimination (Schur factorization) of
the simplified analytic Jacobian.  The resulting matrices have a parabolic form, even if
the starting system is hyperbolic, as a result of the implicit timestepping
procedure. Multigrid methods are often used to solve these systems (see, e.g.,
\cite{chacon2008}), but these are difficult to implement in complex geometries and an
existing codebase that was not designed with them in mind. Here we describe a simpler
scheme that can be used to efficiently solve the factorized Jacobian in field-aligned
coordinates without the complexity of a multigrid method.

To derive a preconditioner, we begin by simplifying the nonlinear six-field system,
equations~\eqref{eq:6field1}-\eqref{eq:6field3}, retaining the linear vorticity terms
driving perpendicular plasma $E\times B$ motion, the fast parallel Alfv\'en wave that
propagates the motion along fieldlines, the linear terms describing perpendicular
advection of density and temperature, and the parallel diffusion of temperature
perturbations. This procedure reduces the equations to
\begin{align}
  \pd{\rho}{t}   &\simeq -\vec{v}_{E\times B} \cdot \nabla_\perp \rho_0 \cr
  \pd{T_s}{t}    &\simeq -\vec{v}_{E\times B} \cdot \nabla_\perp T_{s0}
                 + \frac{1}{\rho_0} \nabla_{||} \cdot \left(\chi_{s||}\partial_{||}T_s\right) \cr
  \pd{U}{t}      &\simeq \frac{1}{\rho_0} \left[B_0^2\nabla_{||} \left(\frac{J_{||}}{B_0}\right)
                 + 2\vec{b}_0 \times \vec{k}_0 \cdot \nabla P\right] \cr
  \pd{v_{||}}{t} &\simeq 0 \cr
  \pd{A_{||}}{t} &\simeq -\nabla_{||}\phi .
  \label{eq:simplified}
\end{align}

We have therefore neglected all nonlinear terms (those involving two or more evolving
quantities), all perpendicular diffusive terms, viscosity, and resistivity. The neglect of
nonlinear terms is justified because we wish to follow the nonlinear dynamics and thus
should not step over these timescales; instead, we focus on preconditioning the fast
linear dynamics. Perpendicular diffusion and resistive effects are also slow relative to
the parallel dynamics and so can be dropped.  We stress that these simplifications made in
the preconditioner do not affect the solution to the full six-field model, only the
convergence rate toward it.  When the assumptions made here are violated, the
preconditioner becomes less effective, but the solution is not compromised (provided that
the Krylov method still converges and does not stall).

In an implicit timestepping scheme (Section~\ref{subsec:newtonprecon}) the operator to be
inverted is $\mathbb{A} = \mathbb{I} - \gamma\mathbb{J}$, where $\mathbb{J}$ is the system
Jacobian; $\gamma$ is usually proportional to the timestep and is determined by the choice
of integration scheme. For preconditioning, a simplified form of the Jacobian is derived
from Equations~\eqref{eq:simplified} and is given in Equation \eqref{eq:jacobian}.

{\footnotesize
\begin{equation}
  \bb{J} = \left[
    \begin{array}{c|c|c|c|c|c}
      0 & 0 & 0 & 0 & 0 & \frac{1}{B}\vec{b}\times\nabla\rho_0\cdot\nabla\nabla_\perp^{-2}[B_0 \\

      0 & \frac{\chi_{i,||}}{\rho_0}\nabla_{||}^2 & 0 & 0 & 0 &
      \frac{1}{B}\vec{b}\times\nabla T_{i0}\cdot\nabla\nabla_\perp^{-2}[B_0 \\

      0 & 0 & \frac{\chi_{e,||}}{\rho_0}\nabla_{||}^2 & 0 & 0 &
      \frac{1}{B}\vec{b}\times\nabla T_{e0}\cdot\nabla\nabla_\perp^{-2}[B_0 \\

      0 & 0 & 0 & 0 & 0 & 0 \\

      0 & 0 & 0 & 0 & 0 & -\nabla_{||}\nabla_\perp^{-2}[B_0 \\

      \frac{2}{\rho_0}\vec{b}\times\kappa\cdot\nabla\left[\left(T_{e0}+T_{i0}\right)\right. &
      \frac{2}{\rho_0}\vec{b}\times\kappa\cdot\nabla\left[\rho_0\right. &
      \frac{2}{\rho_0}\vec{b}\times\kappa\cdot\nabla\left[\rho_0\right. &
      0 &
      -\frac{1}{\rho_0}B_0\nabla_{||}\nabla_\perp^2 &
      0
    \end{array}
  \right]
\label{eq:jacobian}
\end{equation}
}

Following the methodology of \cite{chacon2008}, we first factorize the matrix
$\bb{A}=\bb{I} - \gamma\bb{J}$ using Schur factorization:
\begin{equation}
  \bb{A}^{-1} = \left( \bb{I} - \gamma\bb{J} \right)^{-1} =
  \boutmatrix
    {\bb{E}} {\bb{U}}
    {\bb{L}} {\bb{D}}^{-1} =
  \boutmatrix
    {\bb{I}} {\bb{E}^{-1}\bb{U}}
    {0} {\bb{I}}
  \boutmatrix
    {\bb{E}^{-1}} {0}
    {0} {\bb{P}_{Schur}^{-1}}
  \boutmatrix
    {\bb{I}} {0}
    {-\bb{L}\bb{E}^{-1}} {\bb{I}},
\label{eq:factorisation}
\end{equation}
where $\bb{L}$ is a row vector containing the vorticity drive terms due to curvature and
parallel current, $\bb{U}$ is a column vector containing the linear $E \times B$ advection
terms, and $\bb{E}$ contains the parallel heat conduction terms, which have the form
$\bb{I} - \gamma \frac{\chi_{||e,i}}{\rho_0}\nabla^2_{||}$.

The matrix $\bb{P}_{Schur} = \bb{D} - \bb{L}\bb{E}^{-1}\bb{U}$ is the Schur complement;
and since nonlinear terms are neglected, $\bb{D}$ is the identity matrix.  To simplify the
form of $\bb{P}_{Schur}$, we drop terms involving only perpendicular derivatives, leaving
only the shear Alfv\'en wave:
\[
  \bb{P}_{Schur} \simeq \bb{I} - \gamma^2\frac{B_0^2}{\rho_0}
  \nabla_{||}\nabla_\perp^2\nabla_{||}\nabla_\perp^{-2} .
\]
In order to remove the perpendicular derivatives, the following approximation is used:
\[
  \nabla_{||}\nabla_\perp^2\nabla_{||}\nabla_\perp^{-2}
  \simeq
  \nabla_{||}^2 -\frac{2\nabla_{||}\left(RB_\theta\right)}{RB_\theta} ,
\]
where $R$ and $B_\theta$ are the equilibrium major radius and poloidal magnetic field,
respectively.

The above factorization contains three applications of $\bb{E}^{-1}$; however, it was
found that only the second had a significant impact on the convergence rate, and so the
first and last were omitted, significantly speeding the calculation. The resulting
preconditioning operator can be divided into three stages, $\bb{P} = \bb{P}_3 \bb{P}_2
\bb{P}_1$, and written as an operation on a unpreconditioned state vector $v^{(U)}$, which
produces a state vector $v^{(P)} = \bb{P} v^{(U)}$. In the following description,
superscripts on state vectors or components of state vectors will indicate the stage of
the preconditioning operation.

Note that in the simplest backwards Euler implicit scheme when $\bb{P} = \bb{A}^{-1}$, the
vector $v^{(U)}$ is the solution at the last timestep, and $v^{(P)}$ the solution at the
next.  The preconditioner can therefore be interpreted as a predictor step, calculating an
approximation of the state vector at the next timestep from the previous one. The
application of $\bb{P}$, Equation \eqref{eq:factorisation}, is performed as follows.
\begin{enumerate}
  \item Given an unpreconditioned vector $v^{(U)} = \left(\rho^{(U)}, T_i^{(U)},T_e^{(U)},
    v_{||}^{(U)}, A_{||}^{(U)}, U^{(U)}\right)$, apply the first matrix in
    Equation~\eqref{eq:factorisation}:
    \[
      v^{(P1)} = \boutmatrix
                  {\bb{I}} {0}
                  {-\bb{L}} {\bb{I}}
                v^{(U)} ,
    \]
    where the $\bb{E}^{-1}$ operator has been dropped, as discussed above. When this is
    expanded, only the vorticity field is modified:
    \[
      U^{(P1)} = U^{(U)} + \gamma \frac{B_0}{\rho_0}\nabla_{||}J_{||}^{(U)}
              + \gamma \frac{1}{\rho_0}2\vec{b}\times\kappa\cdot\nabla P^{(U)} ,
    \]
    where $J_{||}^{(U)} = -\nabla_\perp^2A_{||}^{(U)}$ and $P^{(U)} =
    \rho_0\left(T_i^{(U)} + T_e^{(U)}\right) + \rho^{(U)}\left(T_{i0} +
      T_{e0}\right)$. This operation represents a forward Euler step to predict the
    vorticity after a timestep $\gamma$.

  \item Apply the second matrix, solving for parallel temperature conduction and shear
    Alfv\'en wave dynamics,
    \[
    v^{(P2)} = \boutmatrix
                {\bb{E}^{-1}} {0}
                {0} {\bb{P}_{Schur}^{-1}}
              v^{(P1)} ,
    \]
    which can be written as
    \begin{align*}
      T_{s}^{(P2)} &= \left(
                       \bb{I} - \gamma \frac{\chi_{||s}}{\rho_0}\nabla^2_{||}
                     \right)^{-1}T_{s}^{(U)} \cr
      U^{(P)} &= \left(
                  \bb{I} + \frac{2\gamma^2\nabla_{||}(RB_\theta)}{RB_\theta B_0^2}
                  - \gamma^2\frac{B_0^2}{\rho_0}\nabla_{||}^2
                \right)^{-1}U^{(P1)} .
    \end{align*}
    Since all these operators are of the form $\left(a + b\nabla_{||}^2\right)^{-1}$ with
    different scalars $a$ and $b$, the same parallel solver can be used. Since the
    equations for temperature and vorticity are decoupled, a further optimization, not
    used in this current work, would be to also solve the three equations at the same
    time.
    This operation represents a backward Euler step for the propagation of vorticity and
    temperature along magnetic fieldlines over a timestep $\gamma$ and is the most
    computationally demanding.  Efficiently solving these systems is critical to the
    performance of the preconditioner, so the method used is described in detail below.

  \item Apply the third matrix, updating the other fields based on the new vorticity
    \[
      v^{(P)} = \boutmatrix
                 {\bb{I}} {-\bb{U}}
                 {0} {\bb{I}}
               v^{(P2)} ,
    \]
    where again the $\bb{E}^{-1}$ operator has been dropped. This can be written as
    \begin{align*}
      \rho^{(P)} &= \rho^{(U)} - \gamma\vec{v}_{E\times B}^{(P)}\cdot\nabla\rho_0 \cr
      T_s^{(P)} &= T_s^{(P2)} - \gamma \vec{v}_{E\times B}^{(P)}\cdot\nabla T_{s0} \cr
      A_{||}^{(P)} &= A_{||}^{(U)} - \gamma\nabla_{||}\phi^{(P)} ,
    \end{align*}

    where $\vec{v}_{E\times B}^{(P)} = \frac{1}{B_0} \vec{b}_0\times\nabla\phi^{(P)}$ is
    the perpendicular $E\times B$ velocity calculated by using $\phi^{(P)} =
    \nabla_\perp^{-2}\left(B_0 U^{(P)}\right)$, the electrostatic potential calculated
    from the updated vorticity.  This represents an implicit update of $\rho,T_s,A_{||}$,
    advecting plasma quantities based on an approximate $E\times B$ velocity at the next
    timestep.

    The only variable that is not modified in any of these operations is the parallel
    velocity, so $v_{||}^{(P)} = v_{||}^{(U)}$.
\end{enumerate}

The preconditioner described above contains the key physics responsible for fast
timescales in the system, and it has greatly reduced the complexity and computational cost
of calculating $\bb{P}$.  These features can be seen by considering the sizes of the
matrices involved.  Since 6 variables are evolved on $N$ mesh points, the matrix $\bb{A}$
to be inverted has $(6N)^2$ elements, most of which are however zero. A typical simulation
will use on the order of 100 points in each dimension (between 64 and 512 in most cases,
though the radial mesh size $N_\psi$ can be several thousand for ITER simulations) and so
have $\mathord{\sim} 10^6$ mesh points, requiring the inversion of a $\left(6\times 10^6\right)^2$
matrix at each timestep. By block factorizing and keeping only the fast components, we
have reduced this to solving three separate matrices (in $\bb{E}$ and $\bb{P}_{Schur}$),
which can be performed in parallel, each of which is of size $(N)^2 \simeq (10^6)^2$. This
approach in itself does not represent a large savings; but by exploiting the structure of
these blocks in the field-aligned coordinate system employed by BOUT++, further
simplifications can be made.

Since each block contains only derivatives along magnetic fieldlines after we remove
perpendicular derivatives, and equilibrium magnetic fieldlines form helices around a torus
covering 2D ``flux'' surfaces, we can independently solve for each flux surface,
immediately reducing the $(10^6)^2$ inversion problem to $\mathord{\sim} 100$ separate $(10^4)^2$
problems. This assumption of equilibrium flux surfaces will break down if the magnetic
field becomes significantly perturbed in the later stages of an ELM crash, but this is not
a critical problem because at these times nonlinear dynamics are fast and a smaller
timestep is unavoidable.

To further simplify the problem, we exploit the toroidal symmetry of the equilibrium and
decompose into toroidal Fourier harmonics. This approach decouples the toroidal and
poloidal directions, splitting the 2D problem into a set of one-dimensional equations in
poloidal angle, with a complex phase shift representing the pitch of the fieldlines. If we
started with $N_\psi$ flux surfaces, each containing $N_\theta$ poloidal points and
$N_\zeta$ toroidal points (so $N = N_\psi\times N_\theta\times N_\zeta$), then we are left
with $N_\psi\times \left(N_\zeta/2\right) \sim 10^4$ complex cyclic tridiagonal systems,
each of length $N_\theta \sim 10^2$. These are then solved efficiently in parallel by
using a variant of the cyclic reduction algorithm with interface equations
\cite{austin2004}. Since the three equations being solved are also independent, all
$3\times N_\psi\times \left(N_\zeta/2\right)$ tridiagonal systems can be solved
simultaneously, resulting in a good overlap of communication and computation on a large
number of processors.

\subsection{Nonlinear GMRES}
\label{subsec:ngmres}

As an alternative to preconditioned Newton-Krylov methods, we consider nonlinear GMRES
(NGMRES) to solve the nonlinear system~\eqref{eq:f=0}.  NGMRES and other similar methods
(quasi-Newton, Anderson mixing) of solution for nonlinear systems are analogous to
multistep Krylov methods used for linear systems \cite{Oosterlee:2000tg,Anderson:1965wp}.
At each iteration, an updated solution is constructed as a linear combination of several
stored previous solutions and a new candidate solution.  A minimum-residual linear
combination is found by solving a small minimization problem.

The new candidate solution is constructed by perturbing the current solution by some
negative damping of the current residual, typically found through a line search.  One may
also construct the candidate using a small number of iterations of some other nonlinear
solution technique. For the present problem we use a damping factor of $-0.001$ on the
residual to construct the candidate and store 50 previous solutions. By using these Krylov
subspace accelerators, we aim to reduce the overall CPU time for each nonlinear solve and
also reduce memory usage.

\section{Experimental Results}
\label{sec:results}

We use the nonlinear solvers in the PETSc library \cite{petsc-user-ref,efficient} to solve
these nonlinear systems over varying numbers of processors.  The library provides a simple
interface for nonlinear solvers that enables selection of algorithms and parameters at
runtime, with no changes to the application code.

Our experiments were run on the Fusion cluster at Argonne National Laboratory.  This
cluster has 320 compute nodes, each with an Intel Nehalem dual-socket, quad-core 2.53 GHz
processor (or 2,560 compute processors), connected by Infiniband.

We specified a relative nonlinear convergence tolerance of $10^{-3}$ for all methods,
which is consistent to what is typically employed for implicit timestepping.  For each
Newton-Krylov run, we used GMRES with a restart of 30 and a relative linear convergence
tolerance of $10^{-2}$.  The method NGMRES, by definition, has no linear solver associated
with it.  Consequently, each nonlinear solve using NGMRES has many more iterations, in
contrast with the Newton-Krylov approaches, where an average of two or three nonlinear
iterations achieved convergence.

\subsection{Convergence of Nonlinear Solvers}
\label{subsec:convergencenonlinear}

To test the feasibility of our approach to physics-based preconditioning, we conducted
initial experiments with a simpler problem: a two-field subset of the six-field model that
contains only the Alfv\'en wave,
\begin{align}
  \pd{A_{||}}{t} &= -\nabla_{||}\phi - \eta J_{||} \cr
  \pd{U}{t} &= -\vec{v}_{E\times B}\cdot\nabla U
             +\frac{1}{\rho_0} B_0^2\nabla_{||} \left(\frac{J_{||} + J_{||0}}{B_0}\right) ,
  \label{eq:2field}
\end{align}
with auxiliary equations
\[
  U = \frac{1}{B}\nabla_\perp^2\phi \qquad J_{||} = -\nabla_\perp^2 A_{||} .
\]
Though greatly simplified, this model can describe the essential physics of magnetic
reconnection and island formation in tokamak plasmas~\cite{fitzpatrick1998}. Thus, the
efficient solution of this system of equations is of interest in itself. We used a
rectangular (slab) grid of dimension $68\times 32 \times 16$.  As seen in
Figure~\ref{fig:intro-convergence-16}, physics-based preconditioned JFNK showed a significant
reduction in overall runtime in comparison with the default (JFNK with no
preconditioning). NGMRES also significantly decreased overall runtime.

\begin{figure}[ht]
  \centering
  \begin{tikzpicture}
    \begin{semilogyaxis}[
      ylabel={Function Norm},
      xlabel={Wall-clock Time (sec)},
      cycle list name=bwcolors,
      legend style={
        legend columns=1,
        align=left,
        nodes=right,
      },]
      \addplot table {reconnect-2field-none-16.snes};
      \addplot table {reconnect-2field-ngmres-16.snes};
      \addplot table {reconnect-2field-precon-16.snes};

      \addlegendentry{JFNK with no\\ preconditioning};
      \addlegendentry{Nonlinear GMRES};
      \addlegendentry{JFNK with physics-\\based preconditioning};
    \end{semilogyaxis}
  \end{tikzpicture}
  \caption{Function norm ($l_2$) versus wall-clock time (sec) on 16 cores for two-field
    test case, Equation~\ref{eq:2field}. Labels on the figure indicate cumulative number
    of linear solver iterations for Newton-Krylov methods.}
  \label{fig:intro-convergence-16}
\end{figure}

For the remaining experiments, we used the six-field model, Equations
\eqref{eq:6field1}-\eqref{eq:6field3}, on a grid size of $N_\psi=64$, $N_\theta=64$, and
$N_\zeta=128$, partitioned only in the $\psi\theta$-plane because of using the discrete
Fourier transform in the $\zeta-$direction. Each example has six coupled independent
variables $(\rho, T_e, T_i, U, v_{||}, A_{||})$, as given by Equation~\eqref{eq:f=0}. One
reason for using a relatively small grid is to test the scaling of these algorithms when
the number of grid cells per processor becomes small. Since domain decomposition is in
$\psi$ and $\theta$, 512 processors (the maximum used here) will divide the mesh into $4\times 2\times 128$
blocks. This is a more useful test than simply using a large mesh to demonstrate scaling.

\begin{figure}[ht!]
  \centering
  \subfloat[$t=0$]{
    \begin{tikzpicture}
      \begin{semilogyaxis}[
        ylabel={Function Norm},

        xlabel={Wall-clock Time (sec)},
        legend to name=convergence-legend,
        cycle list name=bwcolors,
        ]
        \addplot table {t0-none-512.snes};
        \addplot table[each nth point={10}] {t0-ngmres-512.snes};
        \addplot table {t0-precon-512.snes};

        \addlegendentry{JFNK with no\\ preconditioning};
        \addlegendentry{Nonlinear GMRES};
        \addlegendentry{JFNK with physics-\\based preconditioning};
      \end{semilogyaxis}
    \end{tikzpicture}
    \label{fig:convergence0}
  }
  \subfloat[$t=40$]{
    \begin{tikzpicture}
      \begin{semilogyaxis}[
        xlabel={Wall-clock Time (sec)},
        cycle list name=bwcolors,
        ]
        \addplot table {t40-none-512.snes};
        \addplot table[each nth point={40}] {t40-ngmres-512.snes};
        \addplot table {t40-precon-512.snes};
      \end{semilogyaxis}
    \end{tikzpicture}
    \label{fig:convergence40}
  }\\
  \vspace{3ex}
  \hspace{3ex}
  \pgfplotslegendfromname{convergence-legend}
  \caption{Function norm ($l_2$) versus wall-clock time (sec) on 512 cores for the full
    six-field model (Section~\ref{subsec:equations}) at fixed problem size.  {\em Left:}
    Results for the initial timestep $t=0$. {\em Right:} Results for timestep $t=40$.
    Labels indicating the cumulative number of linear solver iterations for Newton-Krylov
    methods show that the linearized Newton system is relatively easy to solve initially,
    though it becomes more difficult at timestep $t=40$.}
\end{figure}

As shown in Figures \ref{fig:convergence0} and \ref{fig:convergence40}, physics-based
preconditioning achieves an order of magnitude decrease in overall time for the
Newton-Krylov nonlinear solve and nearly two orders of magnitude reduction in the number
of Krylov iterations in comparison with the unpreconditioned variant.  NGMRES also
significantly reduces overall time for the nonlinear solve compared with the BOUT++
default of JFNK with no preconditioning.  Note that while the nonlinear relative
convergence tolerance for each run is $10^{-3}$, the final value of the function norm for
the Newton-Krylov runs is much smaller compared with the NGMRES run because of the
substantial decrease in function norm at each Newton step.  While the physics-based
preconditioner code is hand-written for each type of physics model, NGMRES requires no
application-specific code and can be activated at runtime, making it an attractive
alternative for efficient and scalable nonlinear solves.

At early simulation times ($t=0$, Figure~\ref{fig:convergence0}), nonlinear terms in the
evolution equations are negligible, and so the solution is essentially exponential
growth. At later times ($t=40 / \omega_A$, Figure~\ref{fig:convergence40}), nonlinear
terms become important, filaments begin to erupt outwards from the plasma, and the
solution becomes more complicated, even turbulent.  As a result the number of iterations
and wall clock time increase significantly for later timesteps.  The physics-based
preconditioner is marginally less effective at later times, reducing the iteration count
over unpreconditioned JFNK by a factor of $67$ at $t=40$, compared with a factor of $86$
at $t=0$. This result is expected because nonlinear terms were not included in
construction of the preconditioner.  Nevertheless, the linear preconditioner remains
highly effective, indicating that linear dynamics are still limiting the timestep even
when the plasma evolution is nonlinear.

\subsection{Scalability of Nonlinear Solvers}
\label{subsec:scalabilitynonlinear}

\pgfplotstableread{
procs, none, ngmres, precon
16, 6.079748e+03, 1.349796e+03, 1.431308e+02
32, 3.624823e+03, 8.953517e+02, 9.952766e+01
64, 1.685646e+03, 4.164732e+02, 4.088985e+01
128, 9.180763e+02, 2.658314e+02, 2.383816e+01
256, 5.977239e+02, 1.971505e+02, 1.552403e+01
512, 4.410280e+02, 1.653697e+02, 1.161303e+01
}\scalability

\begin{figure}[ht!]
  \centering
  \begin{tikzpicture}
    \begin{groupplot}[
      group style=
      {
        group size=2 by 2,
        xlabels at=edge bottom,
        ylabels at=edge left,
        horizontal sep=0.075\textwidth,
        group name=scalability
      },
      scalabilitybar/.style={
        ybar,
        ymin=0,
        cycle list name=bwcolors-bar,
      },
      xtick=data,
      enlarge x limits=0.2,
      xticklabels from table={\scalability}{procs},
      xmode=log,
      ylabel={Wall-clock Time (sec)},
      xlabel={Number of Processors},
      table/x=procs,
      legend style={
        draw=none,
        nodes=right,
        align=left,
        /tikz/every even column/.append style={column sep=5pt},
      },
      ]

      \nextgroupplot[scalabilitybar]
      \pgfplotsset{cycle list shift=3}
      \addplot table[y=none] {\scalability};
      \pgfplotsset{cycle list shift=-3}
      \addplot table[y=ngmres] {\scalability};
      \pgfplotsset{cycle list shift=-2}
      \addplot table[y=precon] {\scalability};
      \addlegendentry {JFNK with no\\ preconditioning};
      \addlegendentry {Nonlinear GMRES};
      \addlegendentry {JFNK with physics-\\based preconditioning};

      \nextgroupplot[scalabilitybar]
      \pgfplotsset{cycle list shift=2}
      \addplot table[y=ngmres] {\scalability};
      \addlegendentry {Nonlinear GMRES};

      \nextgroupplot[scalabilitybar]
      \addplot table[y=precon] {\scalability};
      \addlegendentry {JFNK with physics-\\based preconditioning};

      \nextgroupplot[
        ymode=log,
        log base 10 number format code/.code={
          $\pgfmathparse{10^(#1)}\pgfmathprintnumber{\pgfmathresult}$
        },
        cycle list name=bwcolors,
        legend pos=south west,
        legend style={
          nodes=right,
          align=left,
          fill=none,
          font=\footnotesize,
        }
      ]

      \addscalingplot [table/y=none] {\scalability};
      \addscalingplot [table/y=ngmres] {\scalability};
      \addscalingplot [table/y=precon] {\scalability};
      \addlegendentry {JFNK with no\\ preconditioning};
      \addlegendentry {Nonlinear GMRES};
      \addlegendentry {JFNK with physics-\\based preconditioning};
    \end{groupplot}
  \end{tikzpicture}
  \caption{The upper left plot compares the scalability of each algorithm for the full
    six-field model.  The upper right and lower left plots are close-ups of NGMRES and
    physics-based preconditioning, respectively. The lower right plot compares each
    algorithm with ideal scaling. We can see that NGMRES and especially physics-based
    preconditioning significantly reduced overall runtime.}
  \label{fig:scalability}
\end{figure}

\pgfplotstableread{
procs, rhs, VecMDot, VecMAXPY
16, 5.9893e+03, 1.0253e+03, 3.9946e+01
32, 3.5834e+03, 5.1720e+02, 1.7514e+01
64, 1.6591e+03, 2.3304e+02, 9.9039e+00
128,9.0413e+02, 1.3117e+02, 4.8912e+00
256,5.8619e+02, 8.6498e+01, 2.9258e+00
512,4.3136e+02, 8.1301e+01, 1.5468e+00
}\logsummarytimenone

\pgfplotstableread{
procs, rhs, VecNorm, VecDot, VecMDot
16, 1.1648e+03, 1.9524e+02, 3.8848e+01, 2.3147e+01
32, 7.7791e+02, 1.4300e+02, 2.8514e+01, 1.0978e+01
64, 3.2883e+02, 7.6011e+01, 2.6913e+01, 5.9854e+00
128, 1.8896e+02, 5.6367e+01, 2.7379e+01, 3.4278e+00
256, 1.1926e+02, 5.0114e+01, 3.0705e+01, 2.1295e+00
512, 8.4222e+01, 5.1892e+01, 3.4524e+01, 1.3097e+00
}\logsummarytimengmres

\pgfplotstableread{
procs, rhs, PCApply, VecMDot, VecNorm
16, 1.2007e+02, 2.1316e+01, 1.6914e+01, 1.3464e+00
32, 7.9991e+01, 1.8606e+01, 1.1552e+01, 9.6238e-01
64, 3.3514e+01, 6.7648e+00, 4.7187e+00, 4.3807e-01
128,1.9236e+01, 4.0946e+00, 2.8131e+00, 3.1593e-01
256,1.2155e+01, 2.8135e+00, 1.9681e+00, 2.8729e-01
512,8.5291e+00, 2.4913e+00, 1.7678e+00, 3.2046e-01
}\logsummarytimeprecon

\begin{figure}[ht!]
  \centering
  \includegraphics{paper-figure4.pdf}
  \caption{Results showing that function evaluation (RHS) dominates the time for the
    six-field model}
  \label{fig:stacked}
\end{figure}

Figure~\ref{fig:scalability} shows the scaling of the JOREK test case with the Jacobian-free
Newton-Krylov (no preconditioning), nonlinear GMRES, and Jacobian-free Newton-Krylov (with
physics-based preconditioning) algorithms in the first subfigure with lighter bars
denoting the comparison with NGMRES and physics-based preconditioning.  Notice that the
physics-based preconditioning is barely legible in the first subfigure because of the
overall runtime being an order of magnitude faster than the default solver.  For this
reason the individual plots of NGMRES and physics-based preconditioning are also displayed
as individual subfigures.  All three algorithms then are plotted in a log-log graph in
conjunction with ideal scaling (the dashed lines). One can see that all three algorithms
have almost identical scaling.  This relatively small test case becomes too small for core
counts larger than 128, as evidenced by the tapering of scalability.

Clearly, using a physics-based preconditioner achieves the fastest runtime.  Physics-based
approaches, however, require detailed knowledge of the modeling equations being used
\cite{chacon2008,chacon02,ChaconKnollEtAl02}, enforcing a tight coupling between the
programmer and physics-based preconditioning analysis. Even if such a coupling exists,
potential challenges arise with implementation because of complexities of physics and
preconditioning algorithms.  In contrast, NGMRES requires no application-specific
customization, since the application provides code only for the nonlinear function
evaluation.

The bar graphs in Figure~\ref{fig:stacked} show the proportion of time spent in various
phases of the nonlinear solve.  For each run, the evaluation of the function (or right-hand
side, RHS) dominates the runtime. For Jacobian-free variants, a function evaluation is
required for each Jacobian-vector product (one per Krylov iteration).  The preconditioned
variant of JFNK requires fewer function evaluations because of better-conditioned linear
systems.  NGMRES, alternatively, needs roughly ten times more function evaluations than
does preconditioned JFNK. Also, we can see that the time needed for computing the search
direction in NGMRES (VecNorm) is nontrivial. This result further indicates that the
parallel scaling shown in Figure~\ref{fig:scalability} is due primarily to the scaling of the
nonlinear function evaluation; improving this scaling is the subject of ongoing work.

\section{Conclusions}
\label{sec:conclusions}

We have introduced a physics-based preconditioner that leverages insight into the physics
of the equations being solved in order to improve convergence of matrix-free Newton-Krylov
methods in BOUT++ ELM simulations.  We also introduced the use of a complementary
approach, nonlinear GMRES, which requires no Jacobian-vector products or preconditioning.
We demonstrated that both approaches reduce time for solving the nonlinear systems that
arise at each timestep of an implicit ELM simulation.

The physics-based preconditioning scheme described in this paper is straightforward to
construct once the analytic Jacobian is known, can be efficiently solved in magnetic
field-aligned coordinate systems commonly used in plasma simulations, and is highly
effective for both hyperbolic (Alfv\'en wave) and parabolic (heat conduction) problems.
The preconditioning steps have intuitive interpretations, making it straightforward to
experiment with adding and removing terms in the equations in order to optimize the
preconditioner.  By applying this method to a six-field reduced MHD model,
order-of-magnitude reductions in wall clock time have been achieved, even in regimes where
nonlinear terms in the evolution equations are important.  The application of this
improved code to ELM simulations, and predictions for ITER, will be reported in a separate
paper.

Using a relatively small problem size, we have demonstrated good scaling up to 512
processors. Since the costly matrix inversion step in the preconditioner is trivially
parallelized in two dimensions ($\psi$ and $\zeta$), increasing the problem size should
allow scaling to a much larger number of processors. Testing this will be the subject of
future work.

Other directions of future work include the development of additional physics-based
preconditioning variants, such as fieldsplits for fast Alfv\'en waves and thermal
conductivity along the fieldlines.  We are also exploring broader issues in
time-dependent simulations, including the use of implicit-explicit (IMEX) approaches for
flexible timestepping. Work is under way to incorporate BOUT++ into the FACETS
framework~\cite{facets-scidac09,facets-fec2010} in order to facilitate research on
core-edge-wall coupling.

\section*{Acknowledgments}

We gratefully acknowledge use of the Fusion cluster in the Laboratory Computing Resource
Center at Argonne National Laboratory.  S. Farley and L. McInnes were supported by the
Office of Advanced Scientific Computing Research, Office of Science, U.S. Department of
Energy, under Contract DE-AC02-06CH11357. B. Dudson was supported by EFDA WP11 and HPC-FF
computing resources. This work was carried out largely as part of the FACETS SciDAC
project supported by DOE. We also thank L.Chacon for his presentation and insightful
discussions at the 2011 BOUT++ workshop, which provided the impetus for this work.

\bibliographystyle{elsarticle-num}


\section{Introduction}

The problem of estimating a signal from the magnitude of its short-time Fourier transform (STFT) has vexed audio signal processing researchers for decades \cite{HayesLimOpp80,NawabQuatieriLim83,Slaney94}. The fundamental challenge can be summarized as:

\begin{quote}
\emph{Fourier coefficients are complex numbers, but in many applications we only know their magnitudes and not their phases.}
\end{quote}

For example, in many frequency-domain audio processing techniques such as time-scale modification and speech enhancement, it is typically straightforward to modify the magnitudes of the Fourier coefficients but not their phases. In the case of time stretching, it is clear that the energy (i.e., magnitudes) should be spread over more time bins, but unclear how to specify the phases accordingly.

The absence of phase information is problematic when we want to reconstruct the modified signal, for which it is necessary to invert the STFT---an operation which requires the full complex coefficients (i.e., both magnitudes and phases). Hence, the magnitude-only reconstruction problem is also termed \emph{phase retrieval} in the literature. Using incompatible phase information in the reconstruction can result in amplitude-modulation distortion \cite{LarocheDolson97}.

Recently, interest in the magnitude-only reconstruction problem has resurged in the source separation community, where matrix factorization has become a standard workhorse \cite{SmaragdisBrown03}. Since these approaches operate entirely on the magnitude spectrogram, the same difficulties arise in reconstructing the time-domain signals. A heuristic that is commonly used in this literature is to set the phases of each of the recovered signals equal to the phases of the original mixture signal \cite{Smaragdis10}. Some justification has been offered for this approach; under certain modeling assumptions, it is a consequence of the Wiener filter reconstruction \cite{Fevotte09}.

Griffin and Lim \cite{GL84} proposed a more general solution to the magnitude-only reconstruction problem. As their approach is now the standard, we briefly review their algorithm in the next section so that we may later contrast it with our approach.

\section{Existing Approaches}

Let $|Y_w(mR,\omega)|$ denote the desired magnitude STFT, where $R$ is the hop size, $w$ denotes the window, and $m$ and $\omega$  index the time frame and frequency, respectively. Our goal is to find a signal $x$, whose magnitude STFT $|X_w(mR,\omega)|$ is as close to the desired as possible, in a least squares sense:
\begin{equation}
\sum_{m=-\infty}^\infty \int_{\omega=-\pi}^\pi \left( |X_w(mR,\omega)| - |Y_w(mR,\omega)| \right)^2\,d\omega
\label{obj}
\end{equation}

The Griffin-Lim algorithm attempts to minimize this objective by starting with an initial guess for the signal $x^{(0)}$ and iterating the following steps until convergence:

\begin{enumerate}

\item Compute the STFT of the current signal estimate $x^{(i)}$:
\begin{equation}
X_w^{(i)}(mR,\omega) := \sum_{n=-\infty}^\infty w(mR-n) x^{(i)}(n) e^{-j\omega n}.
\end{equation}
Retain the phases of $X_w^{(i)}$, but replace its magnitudes by the desired magnitudes:
\begin{equation}
X_w^{(i+1)}(mR,\omega) := |Y_w(mR,\omega)|\frac{X_w^{(i)}(mR,\omega)}{|X_w^{(i)}(mR,\omega)|}.
\end{equation}

\item Compute the IDTFT of $X_w^{(i+1)}(mR,\cdot)$, denoted $\hat x_{w, mR}^{(i+1)}$, for each $m$.\footnote{Ideally we would like for $\hat x_{w, mR}^{(i+1)}$ to be the windowed version of some signal $x$, i.e.
\begin{equation}
\hat x_{w, mR}^{(i+1)}(n) = w(mR-n)x(n).
\label{gl2a}
\end{equation}
To see that this is not always possible, consider the case where $\hat x_{w,mR}(n) \neq 0$ but $w(mR-n) = 0$. Furthermore, we would need a single signal $x$ such that (\ref{gl2a}) holds for all $m$.} Then, find a signal $x$ such that the windowed version of $x$, $x_{w,mR}$, is close to $\hat x_{w,mR}$ for each $m$ in a least-squares sense. The solution turns out to be:
\begin{equation}
x^{(i+1)} := \frac{\sum_{m=-\infty}^\infty w(mR-n)\hat x_{w,mR}^{(i+1)}(n)}{\sum_{m=-\infty}^\infty w^2(mR-n)}.
\end{equation}

\end{enumerate}

Griffin and Lim showed that this algorithm decreases the objective (\ref{obj}) on each iteration and converges to a stationary point. However, the objective is nonconvex, so different initializations $x^{(0)}$ can lead to different solutions, and there are no guarantees on convergence to a global optimum. In fact, as we will show, the algorithm often fails to recover a signal from its true magnitude spectrogram, even when several initializations are used.

It is worth noting that the magnitude-only reconstruction problem has been formulated in other ways: as a probabilistic model \cite{AchanRoweis03} and as a root-finding problem \cite{BouvrieEzzat06}. However, these methods suffer from the same drawbacks as the Griffin-Lim algorithm---namely, that they try to solve a nonconvex problem and hence are prone to local optima. To our knowledge, the algorithm that we describe in the next section is the first convex formulation for this problem.

\section{Convex Formulation}

The Griffin-Lim algorithm resembles a number of earlier alternating projections algorithms \cite{GS72} for solving the basic phase retrieval problem (i.e., where the magnitudes are of only a single Fourier transform). Recently it was shown in \cite{Candes11a} that this basic problem could be recast in terms of estimating a rank-1 matrix $X$ which is the outer product of the signal with itself, i.e., $X = xx^*$. Once $X$ is obtained, the putative signal can be recovered by factorizing $X$, provided that it is rank one.

This technique of \emph{lifting} the problem of estimating a vector $x$ to estimating a matrix $X = xx^*$ is standard in the optimization literature \cite{AGJL05}. Although it drastically increases the dimensions of the problem, it leads naturally to a tractable convex program, as we shall see. In this section, we formulate this problem in the context of the STFT. In the spirit of \cite{Candes11a}, who term their algorithm PhaseLift, we call ours STliFT.

First, we consider the problem of recovering a signal from its true magnitude spectrogram (i.e., when a signal with that exact magnitude spectrogram exists). Suppose we use an FFT of length $N$ for each frame. Letting $s_k(n) = e^{j2\pi n k/N}$, $k=0,...,N-1$ denote the sampled complex sinusoids corresponding to this FFT length, the $k^\t{th}$ Fourier coefficient for the $m^\t{th}$ time frame is $\langle x_{w,mR}, s_k\rangle$, where $x_{w,mR}$ denotes the windowed signal $x_{w,mR}(n) = w(mR-n)x(n)$, zero-padded to the FFT length $N$. Since we want the magnitude of each of these coefficients to be $|Y_w(mR,\omega_k)|$, the problem can be recast as:
\begin{equation}
\begin{aligned}
& \text{find} & & x \\
& \text{subject to} & & \left|\langle x_{w,mR}, s_k \rangle\right| = |Y_w(mR,\omega_k)| \\
& & & \indent 0 \leq k < N,\,\, 0 \leq m < T
\end{aligned}
\end{equation}
Next, we note that we can write the windowed signal $x_{w,mR}$ as a matrix product $W_{mR} x$, where $W_{mR}$ is of the form:
\[ W_{mR} = \begin{pmatrix}
& & & w(-\frac{N}{2}+1) & & & \\
 & 0 & & & \ddots & & & 0 &\\
& & & & & w(\frac{N}{2}) &
\end{pmatrix}.
 \]
Now we square both sides of the equality constraint. Using the fact that $x_{w,mR} = W_{mR}x$, we can write the left-hand side as
\begin{align*}
\left|\langle x_{w,mR}, s_k \rangle\right|^2 &=  \t{tr}((W_{mR} x)^T s_k s_k^* (W_{mR} x)) \\
&= \t{tr}( \underbrace{(W_{mR}^T s_k) (W_{mR}^T s_k)^*}_{S_{k,m}} \underbrace{x x^T}_X)
\end{align*}
where in the second step we have used the trace identity $\t{tr}(AB) = \t{tr}(BA)$ to group $x$ and $x^T$ into a single rank-1 matrix $X$. Hence the problem becomes:
\begin{equation}
\begin{aligned}
& \text{find} & & X \\
& \text{subject to} & & \t{rank}(X) = 1,\, X \succeq 0 \\
& & & \t{tr}( S_{k,m} X ) = |Y_w(mR,\omega_k)|^2 \\
& & & \indent 0 \leq k < N,\,\, 0 \leq m < T
\end{aligned}
\label{rank}
\end{equation}
This is equivalent to the original problem and thus also nonconvex, except now the nonconvexity plainly arises from a rank constraint. We can relax (\ref{rank}) to a convex problem by replacing the rank constraint with its convex surrogate, a trace minimization:
\begin{equation}
\begin{aligned}
& \text{minimize} & & \t{tr}(X) \\
& \text{subject to} & & X \succeq 0 \\
& & & \t{tr}( S_{k,m} X ) = |Y_w(mR,\omega_k)|^2 \\
& & & \indent 0 \leq k < N,\,\, 0 \leq m < T
\end{aligned}
\label{noiseless}
\end{equation}
This problem is not only convex; it is in fact the dual semidefinite program (SDP) described in \cite{BL96}.

Now we consider the case where there is not necessarily any signal with the specified magnitude spectrogram, and we wish to find a signal whose magnitude spectrogram comes as close to the specification as possible. Because the above problem is formulated for the \emph{power} spectrogram rather than the \emph{magnitude} spectrogram, it cannot readily be adapted to minimize the distance between the magnitude spectrograms as in (\ref{obj}). However, it is straightforward to extend the above framework to minimize the distance between the power spectrograms, by solving
\begin{equation}
\begin{aligned}
& \text{minimize} & & \sum_{k,m} (\t{tr}(S_{k,m}X)-|Y_w(mR,\omega_k)|^2)^2 \\
& & & \indent+ \lambda\cdot\t{tr}(X) \\
& \text{subject to} & & X \succeq 0 \\
\end{aligned}
\label{noisy}
\end{equation}
where $\lambda$ controls the tradeoff between minimizing the distance between the spectrograms and the low-rank constraint on $X$. In practice, we initialize $\lambda$ large enough so that the solution is $X = 0$ and solve (\ref{noisy}) repeatedly for decreasing values of $\lambda$ until we obtain an (approximately) rank-1 solution.

Before concluding this discussion of the model, we comment that Euclidean distance may be inappropriate as a measure of divergence between the estimated and desired spectrogram. Euclidean distance is symmetric on a linear scale, so if the desired power in a time-frequency bin is $0.1$, then Euclidean distance would favor a reconstruction which assigns a power of $0.01$ to that bin, over one which assigns $0.2$. However, the human ear perceives energy on a logarithmic scale; $0.01$ represents a discrepancy of $20$ dB, whereas $0.2$ represents one of only $6$ dB. For these reasons, logarithmic divergences such as Kullback-Leibler or Itakura-Saito may be preferable to Euclidean distance for audio \cite{Fevotte09}.

For purposes of comparison with Griffin-Lim, we only deal with the Euclidean case in this work. However, our formulation above can be readily extended to these other divergences by simply replacing $(\t{tr}(S_{k,m}X)-|Y_w(mR,\omega_k)|^2)^2$ by the appropriate divergence $D(\t{tr}(S_{k,m}X),|Y_w(mR,\omega_k)|^2)$. Such problems can be solved using the algorithm described in the next section, by modifying the gradient $\nabla L$ appropriately.

\section{Practical Implementation}

In this section, we describe a practical algorithm for solving (\ref{noisy}). Although the problem can in principle be solved using standard solvers, the scale of the problem ($X$ is an $n\times n$ matrix, where $n$ is the length of the signal) makes this all but impossible for any real-world problem.

A projected gradient descent algorithm was proposed to solve the PhaseLift problem \cite{Candes11a}. Letting $L(X)$ denote the objective function in (\ref{noisy}), the algorithm starts with an initial guess $X_0$ and iterates:
\[ X_{k+1} = \mathcal{P}(X_k - t_k\nabla L(X_k)) \]
where $\mathcal{P}$ denotes projection onto the feasible set $X \succeq 0$, $\nabla L$ the gradient, and $t_k$ a ``suitable'' step size. $\mathcal{P}$ performs an eigendecomposition on the matrix $X$ and thresholds all negative eigenvalues to 0.

For Euclidean distance as in (\ref{noisy}), the gradient is:
\[ \nabla L(X) =  2\sum_{k,m} (\t{tr}(S_{k,m}X) - |Y_w(mR,\omega_k)|^2)S_{k,m} + \lambda I \]
We emphasize that this is a conceptual, rather than a practical, formula. In practice, the structure of each $S_{k,m}$ (its nonzero entries form a block within the $n\times n$ matrix) should be exploited for fast computation. Such optimizations are implemented in the Matlab code for this algorithm, which will be made available at the first author's webpage.\footnote{{\tt http://www-stat.stanford.edu/$\sim$dlsun}}

Although gradient descent scales well to large problems, it can converge slowly for ill-conditioned problems, where the direction of steepest descent can be nearly orthogonal to the direction of the minimum. This results in the characteristic ``zig-zag'' trajectory of gradient descent \cite{BV}. One way to accelerate the convergence is to evaluate the gradient not at $X_k$ but at an auxiliary point $Y_k$ which is a function of the past values of $X$, in order to capture ``momentum'' of the trajectory. \cite{Candes11a} suggests a variant of the FISTA acceleration scheme \cite{FISTA}:
\begin{align}
X_{k+1} &= \mathcal{P}(Y_k - t_k \nabla L(Y_k)) \nonumber\\
\theta_{k+1} &= 2 \left( 1 + \sqrt{1 + 4/\theta_k^2} \right)^{-1} \nonumber \\
\beta_{k+1} &= \theta_{k+1}(\theta^{-1}_k - 1) \nonumber\\
Y_{k+1} &= X_{k+1} + \beta_k (X_{k+1} - X_k)
\label{nesterov}
\end{align}
We have found this algorithm to scale to reasonably large problems and have confirmed that the accelerations defined by (\ref{nesterov}) drastically improve the performance of projected gradient descent in our setting.

The final piece of the algorithm is the step size $t_k$. Although $t_k$ can be chosen using a line search, evaluating the objective is prohibitively expensive for large problems since it involves an eigendecomposition. We have found choosing a fixed step size by trial and error to work well in practice. For step sizes that are too large, the objective values diverge within a few iterations, so we use the largest step size for which the objective values still converge.

\section{Results}

First, we consider the ``noiseless'' problem of recovering a signal from its true magnitude STFT. We first generated 20 random signals each of length $n=16$ and $n= 32$. We then determined the magnitude STFT of each signal using Hann windows of length $M=5,7,9,11$ for both the $n=16$ and $n=32$ signals, and also windows of length $M=13,15,17,19,21$ for the $n=32$ signals. We considered all hop sizes that resulted in constant overlap-add for each window length. The FFT size was taken to be the same as the signal length: $N=n$.

We then applied the Griffin-Lim algorithm and the program (\ref{noiseless}) to each magnitude spectrogram to obtain a signal estimate. For Griffin-Lim, we used 10 random initializations and retained only the ``best'' (see below) of the 10 initializations. To evaluate performance, we tabulated the percentage of the 20 signals whose estimates achieved a value of (\ref{obj}) less than $10^{-3}$. Results are presented in Table \ref{comp16}.

\begin{table}

{\small\begin{tabular}{cccc}
\hline
\multicolumn{4}{c}{$n=16$} \\
\hline
\parbox{.67in}{Window Length ($M$)} & \parbox{.47in}{Hop Size ($R$)} & \parbox{.6in}{Accuracy (\%) of G-L} & \parbox{.6in}{Accuracy (\%) of (\ref{noiseless})} \\
\hline
5 & 2 & 0 & 65 \\
5 & 1 & 30 & 100 \\
\hline
7 & 3 & 0 & 50 \\
7 & 2 & 10 & 100 \\
7 & 1 & 60 & 100 \\
\hline
9 & 4 & 5 & 60 \\
9 & 2 & 70 & 100 \\
9 & 1 & 80 & 100 \\
\hline
11 & 5& 5 & 60 \\
11 & 2 & 95 & 100 \\
11 & 1 & 95 & 100 \\
\hline
\hline
\multicolumn{4}{c}{$n=32$} \\
\hline
\parbox{.67in}{Window Length ($M$)} & \parbox{.47in}{Hop Size ($R$)} & \parbox{.6in}{Accuracy (\%) of G-L} & \parbox{.6in}{Accuracy (\%) of (\ref{noiseless})} \\
\hline
5 & 2 & 0 & 20 \\
5 & 1 & 0 & 100 \\
\hline
7 & 3 & 0 & 15 \\
7 & 2 & 0 & 100 \\
7 & 1 & 15 & 100 \\
\hline
9 & 4 & 0 & 15 \\
9 & 2 & 10 & 100 \\
9 & 1 & 20 & 100 \\
\hline
11 & 5& 0 & 20 \\
11 & 2 & 55 & 100 \\
11 & 1 & 55 & 100 \\
\hline
13 & 6 & 0 & 55 \\
13 & 4 & 0 & 100 \\
13 & 3 & 45 & 100 \\
13 & 2 & 80 & 100 \\
13 & 1 & 65 & 100 \\
\hline
15 & 7 & 0 & 55 \\
15 & 2 & 85 & 100 \\
15 & 1 & 80 & 100 \\
\hline
17 & 8 & 0 & 60 \\
17 & 4 & 60 & 100 \\
17 & 2 & 80 & 100 \\
17 & 1 & 75 & 100 \\
\hline
19 & 9 & 0 & 60 \\
19 & 6 & 0 & 100 \\
19 & 3 & 85 & 100 \\
19 & 2 & 85 & 100 \\
19 & 1 & 90 & 100 \\
\hline
21 & 10 & 0 & 20 \\
21 & 5 & 65 & 100\\
21 & 4 & 90 & 100 \\
21 & 2 & 85 & 100 \\
21 & 1 & 100 & 100 \\
\hline
\end{tabular}}

\caption{\small Comparison of Griffin-Lim and program (\ref{noiseless}) for different configurations in the ``noiseless'' setup. Accuracy is the percentage of the 20 test signals for which the objective (\ref{obj}) of the solution was $< 10^{-3}$.}
\label{comp16}
\end{table}

\begin{figure}

\includegraphics[width=.5\textwidth]{gl.png}

\includegraphics[width=.5\textwidth]{us.png}

\caption{We applied the Griffin-Lim algorithm and program (\ref{noiseless}) to a random signal of length $n=32$ with a window of length $M=15$ and a hop size of $R=7$. Shown at top is an instance where Griffin-Lim fails to converge to the global optimum of (\ref{obj}) and returns an unreasonable solution. Shown at bottom: our algorithm basically recovers the signal from the magnitude spectrogram, up to the numerical precision of the algorithm.}

\label{signals}

\end{figure}

We found that our algorithm recovers the signal more consistently than the Griffin-Lim algorithm. Furthermore, the objective values attained by the latter when it did not converge to the optimum were much larger than $1$, which empirically supports the claim that alternating projections algorithms have a tendency to get stuck in local optima for nonconvex problems. Recovery is a problem for both algorithms when the hop size is large, but especially for Griffin-Lim, most likely because there are fewer constraints with a large hop size, so the landscape is dotted with more local optima.

Next, we consider a more realistic ``noisy'' setup in which there may not be any signal with the given magnitude STFT and we wish only to find a signal minimizing (\ref{obj}). We take the same 20 signals as before, calculate their power spectrogram, and add i.i.d. Gaussian noise with mean 0 and standard deviation 0.2. (We set to zero any entries which end up negative as a result.) We then apply the Griffin-Lim algorithm and the program (\ref{noisy}) to this modified magnitude STFT, and compare the values of (\ref{obj}) for the two estimates. Results are presented in Table \ref{comp32}.

Here the advantages of our algorithm are not so clear cut. Griffin-Lim tends to perform better for large hop sizes, whereas our algorithm does better for smaller hop sizes. However, the two algorithms appear to become increasingly similar as window length increases.

\begin{table}

{\small\begin{tabular}{cccc}
\hline
\multicolumn{4}{c}{$n=16$} \\
\hline
\parbox{.67in}{Window Length ($M$)} & \parbox{.47in}{Hop Size ($R$)} & \parbox{.74in}{Median Rel.\\\% Err.\,of G-L} & \parbox{.7in}{Median Rel.\\\% Err.\,of (\ref{noisy})} \\
\hline
5 & 2 & 1.18 & 5.66 \\
5 & 1 & 0.82 & 0.86\\
\hline
7 & 3 & 1.25 & 4.27 \\
7 & 2 & 1.85 & 1.00 \\
7 & 1 & 0.44 & 0.35 \\
\hline
9 & 4 & 1.03 & 4.32 \\
9 & 2 & 0.66 & 0.41\\
9 & 1 & 0.31 & 0.27 \\
\hline
11 & 5& 0.85 & 5.60 \\
11 & 2 & 0.35 & 0.31 \\
11 & 1 & 0.20 & 0.21\\
\hline
\hline
\multicolumn{4}{c}{$n=32$} \\
\hline
\parbox{.67in}{Window Length ($M$)} & \parbox{.47in}{Hop Size ($R$)} & \parbox{.74in}{Median Rel.\\\% Err.\,of G-L} & \parbox{.7in}{Median Rel.\\\% Err.\,of (\ref{noisy})} \\
\hline
5 & 2 & 0.89 & 16.7 \\
5 & 1 & 1.72 & 0.22 \\
\hline
7 & 3 & 1.00 & 3.91 \\
7 & 2 & 2.10 & 0.24 \\
7 & 1 & 0.77 & 0.12 \\
\hline
9 & 4 & 0.46 & 7.80 \\
9 & 2 & 1.43 & 0.11 \\
9 & 1 & 1.07 & 0.06 \\
\hline
11 & 5& 0.52 & 8.65 \\
11 & 2 & 0.74 & 0.09 \\
11 & 1 & 0.67 & 0.05 \\
\hline
13 & 6 & 0.40 & 9.82 \\
13 & 4 & 1.06 & 0.15 \\
13 & 3 & 1.15 & 0.12 \\
13 & 2 & 0.39 & 0.08 \\
13 & 1 & 0.05 & 0.04 \\
\hline
15 & 7 & 0.40 & 4.62 \\
15 & 2 & 0.07 & 0.07 \\
15 & 1 & 0.04 & 0.04 \\
\hline
17 & 8 & 0.39 & 5.77 \\
17 & 4 & 0.12 & 0.11 \\
17 & 2 & 0.06 & 0.06 \\
17 & 1 & 0.04 & 0.04 \\
\hline
19 & 9 & 0.31 & 3.08 \\
19 & 6 & 0.93 & 0.14 \\
19 & 3 & 0.09 & 0.08 \\
19 & 2 & 0.05 & 0.05 \\
19 & 1 & 0.04 & 0.03 \\
\hline
21 & 10 & 0.51 & 4.83 \\
21 & 5 & 0.15 & 0.14\\
21 & 4 & 0.09 & 0.09 \\
21 & 2 & 0.05 & 0.05 \\
21 & 1 & 0.03 & 0.03 \\
\hline
\end{tabular}}

\caption{\small Comparison of Griffin-Lim and program (\ref{noisy}) for different configurations in the ``noisy'' setup. The median relative error (defined as (\ref{obj}) divided by the total power $\sum_{m,k} |Y_w(mR,\omega_k)|^2$) over 10 test signals is shown.}
\label{comp32}

\end{table}

\section{Discussion}

Although the algorithm is promising on small examples, it remains to scale it to real audio data. To appreciate the difficulties of scaling the algorithm, note that the algorithm estimates an $n \times n$ matrix, where $n$ is the length of the audio clip. For a 3-second clip sampled at 44.1 kHz, this amounts to estimating a $10^5 \times 10^5$ matrix.

We make several observations which may be useful for scaling the algorithm. First, we note that the algorithm is entirely parallelizable in the sense that we can divide the signal into smaller segments and run the algorithm separately on each chunk. The only thing that we lose is the information that would have been provided by the time frames that stretch across two segments. Furthermore, it may not be necessary to compute the phase over the entire signal. Since phase is not audible over stationary stretches, one could apply the algorithm only on the transient time frames and allow phase to ``free run'' elsewhere. Since transients are short in duration and often involve smaller windows, these problems will be much smaller in scale.

\section{Conclusion}

This paper has introduced a novel method for estimating a signal from its magnitude STFT. Although it remains to scale the algorithm to real-life audio applications, our approach is highly competitive with existing algorithms in stylized simulations.

\section{Acknowledgements}

We would like to thank Rahul Mazumder for helpful discussions.