A Deep Conjugate Direction Method for Iteratively Solving Linear Systems
Authors: Ayano Kaneda
Osman Akar
Jingyu Chen
Victoria Alicia Trevino Kala
David Hyde
Joseph Teran
Venue: Submitted to ICLR 2023 (2023)
Keywords: Computational Linear Algebra
Convolutional Neural Network
Conjugate Gradients
Partial Differential Equations
Fluid Simulation

Abstract: We present a novel deep learning approach to approximate the solution of large, sparse, symmetric, positive-definite linear systems of equations.These systems arise from many problems in applied science, e.g., in numerical methods for partial differential equations. Algorithms for approximating the solution to these systems are often the bottleneck in problems that require their solution, particularly for modern applications that require many millions of unknowns. Indeed, numerical linear algebra techniques have been investigated for many decades to alleviate this computational burden. Recently, data-driven techniques have also shown promise for these problems. Motivated by the conjugate gradients algorithm that iteratively selects search directions for minimizing the matrix norm of the approximation error, we design an approach that utilizes a deep neural network to accelerate convergence via data-driven improvement of the search directions.
Our method leverages a carefully chosen convolutional network to approximate the action of the inverse of the linear operator up to an arbitrary constant. We train the network using unsupervised learning with a loss function equal to the $L^2$ difference between an input and the system matrix times the network evaluation, where the unspecified constant in the approximate inverse is accounted for. 
We demonstrate the efficacy of our approach on spatially discretized Poisson equations with millions of degrees of freedom arising in computational fluid dynamics applications.Unlike state-of-the-art learning approaches, our algorithm is capable of reducing the linear system residual to a given tolerance in a small number of iterations, independent of the problem size.Moreover, our method generalizes effectively to various systems beyond those encountered during training.

Under review as a conference paper at ICLR 2023
A DEEP CONJUGATE DIRECTION METHOD FOR ITERATIVELY SOLVING LINEAR SYSTEMS
Anonymous authors Paper under double-blind review
ABSTRACT
We present a novel deep learning approach to approximate the solution of large, sparse, symmetric, positive-definite linear systems of equations. These systems arise from many problems in applied science, e.g., in numerical methods for par- tial differential equations. Algorithms for approximating the solution these systems are often the bottleneck in problems that require their solution, particu- larly for modern applications that require many millions of unknowns. Indeed, numerical linear algebra techniques have been investigated for many decades to alleviate this computational burden. Recently, data-driven techniques have also shown promise for these problems. Motivated by the conjugate gradients algo- rithm that iteratively selects search directions for minimizing the matrix norm the approximation error, design approach that utilizes deep neural network to accelerate convergence via data-driven improvement of the search directions. Our method leverages carefully chosen convolutional network approximate the action of the inverse of the linear operator up to an arbitrary constant. We train the network using unsupervised learning with loss function equal the L difference between an input and the system matrix times the network evaluation, where the unspecified constant in the approximate inverse is accounted for. We demonstrate the efficacy of our approach on spatially discretized Poisson equa- tions with millions of degrees of freedom arising computational fluid dynamics applications. Unlike ofte-oof-theart learning approaches, our algorithm is capable of reducing the linear system residual to a given tolerance small number of iterations, independent of the problem size. Moreover, method generalizes effectively to various systems beyond those encountered during training.
INTRODUCTION
In this work, we consider sparse linear systems that arise from discrete Poisson equations in incom-
pressible flow applications (Chorin 1967 Fedkiw et al. 2001 Bridson 2008). We use the notation
= b
(1)
where the dimension of the matrix A _=x" and the vector b E R" correlate with spatial fidelity of the computational domain. The appropriate numerical linear algebra technique depends on the nature of the problem. Direct solvers that utilize matrix factorizations (QR, Cholesky, etc. Treffthen & Bau )197)) have optimal approximation error, but their computational cost is O(n²) and they typically require dense storage, even for sparse A. Although Fast Fourier Transforms (Nussbaumer) 1981) can be used in limited instances (periodic boundary conditions, etc.), iterative techniques are most commonly adopted for these systems given their sparsity. Many applications with strict per- formance constraints (c.g.. real-time fluid simulation) utilize basic iterations (Jacobi, Gauss-Seide successive over relaxation (SOR), etc.) given limited computational budget (Saad) (003). However, large approximation enrors must be tolerated since iteration counts are limited by the performance constraints. This is particularly problematic since the wide elliptic spectrum of these matrices (a condition that worsens with increased spatial fidelity/matrix dimension) leads to poor conditioning and iteration counts. Iterative techniques can achieve sub-quadratic convergence if their iteration count does not grow excessively with problem size n since each iteration generally requires ((n) floating point operations for sparse matrices. Discrete elliptic operators are typically symmetric positive (semi) definite and the preconditioned conjugate gradients method (PCG) can be used to

---

Under review as a conference paper at ICLR 2023
minimize iteration counts (Sadd 2003 Hestenes & Stiefel 1952 Stiefel 1952). Preconditioners P for PCG must simultaneously: be symmetric positive definite (SPD) (and therefore admit factoriza- tion P = F) improve the condition number of the preconditioned system F FAFy Fb, and be computationally cheap to construct and apply: accordingly, designing specialized precondition- for particular classes of problems is somewhat of art. Incomplete Cholesky preconditioners (ICPCG) [Kershaw] (9788 use sparse approximation to the Cholesky factorization and significantly reduce iteration counts in practice; however, their inherent data dependency prevents efficient par- allel implementation. Nonetheless, these are very commonly adopted for Poisson equations arising in incompressible flow (Fedkiw al. [00T] BBidsoon| 2008). Multigrid (Brandt) 1977) and domain decomposition (Saad 2003) preconditioners greatly reduce iterations counts, but they must be up- dated (with non-trivial cost) each time the problem changes (e.g., in computational domains with time varying boundaries) and/or for different hardware platforms. In general, choice of an optimal preconditioner for discrete elliptic operators is an open area of research.
Recently, data-driven approaches that leverage deep learning techniques have shown promise for solving linear systems. Various researchers have investigated machine learning estimation of multi- grid parameters (Greenfeld et al.) 2019; Grebhahn et al. 2016 Luz et al. 2020). Others have devel-
oped machine learning methods to estimate preconditioners (Götz & Anzt [2018] Stanaityte 2020
Tchimura et al. 2020) and initial guesses for iterative methods (Luna et al. [021] Um et al. 2020
Ackmann et al. (020)_ [Tompson et al.] (217)) and Yang et al.] (2016) develop non-iterative machine
learning approximations of the inverse of discrete Poisson equations from incompressible flow. We
leverage deep learning and develop a novel version of conjugate gradients iterative method for ap- proximating the solution of SPD linear systems which we call the deep conjugate direction method (DCDM). CG iteratively adds A-conjugate search directions while minimizing the matrix norm of the error. We use a convolutional neural network (CNN) as an approximation of the inverse of the matrix in order to generate more efficient search directions. We only ask that our network approxi- mate the inverse up to an unknown scaling since this decreases the degree of nonlinearity and since it does not affect the quality of the search direction (which is scale independent). The network is similar to a preconditioner, but it is not a linear function, and our modified conjugate gradients ap- accommodate this nonlinearity. We use unsupervised learning train our
proach is designed network with a loss function equal to the I.2 difference between an input vector and a scaling of A choose
times the output of network. To account for this unknown scaling during training, the scale of the output of the network by minimizing the matrix norm of the error. Our approach allows for efficient training and generalization to problems unseen (matrices A and right-hand sides
b). We benchmark our algorithm using the ubiquitous pressure Poisson equation (discretized on regular voxelized domains) and compare against FluidNet (Tompson al., (017). which is the state of the art learning- based method for these types of problems. Unlike the non-iterative approaches of Tompson et al. (2017) and Yang et al. (2016), our method can reduce the linear system residuals arbitrarily. We showcase our approach with examples that have over 16 million degrees of freedom.
2 RELATED WORK
Several papers have focused on enhancing the solution of linear systems (arising from discretized PDEs) using learning. For instance, Götz& & Anzt (2018) generate sparsity patterns for block-Jacobi preconditioners using convolutional neural networks, and Staaaityyt (2020) use a CNN to predict non-zero patterns for ILU-ttyp preconditioners for the Navier-Stokes equations (though neither work designs fundamentally new preconditioners), Ichimura et al. (2020) develop a neural-network based preconditioner where the network is used to predict approximate Green's functions (which arise in the analytical solution of certain PDEs) that in turn yield an approximate inverse of the linear system. Hsich et al. (2019) learn an iterator that solves linear systems, performing com- petitively with classical solvers like multigrid-preconditioon preconditioned MINRES (Paige & Saunders] 1975). Luz et al. (2020) and Grrenfeld et al. (2019) use machine learning to estimate algebraic multigrid (AMG) parameters. They note that AMG approaches rely most fundamentally on effectively cho- sen (problem-deeendennt) prolongation sparse matrices and that numerous methods have attempted to automatically create them from the matrix A. They train graph neural network to learn (in an unsupervised fashion) a mapping from matrices A to prolongation operators. Grebhahn et al. (2016) note that geometric multigrid solver parameters can be difficult to choose to guarantee parallel per- formance on different hardware platforms. They use machine learning to create a code generator to help achieve this.
Ichimur
2

---

Under review as a conference paper at ICLR 2023
[
(d)
Figure 1: (a) We illustrate sample flow domain C (0, 1)2 (in 2D for ease of illustration) with in- ternal boundaries (blue lines). (b) We voxelize the domain with regular grid: white cells represent interior/fluid, and blue cells represent boundary conditions. (c) We train using matrix A (0.1)* from discretized domain with no interior boundary conditions, where is the dimension. This creates lin- system with n (n.- 1) unknowns where ne is the number of grid cells each direction. (d) We illustrate the non-zero entries in an example matrix A from the voxelized and labeled (white vs. blue) grid for three example interior cells (green, purple, and brown). Each case illustrates the non-zero entries in the row associated with the example cell. All entries in rows corresponding to boundary/blue cells are zero.
Several works consider accelerating the solution of linear systems by learning an initial guess that close to the true solution or otherwise helpful descent algorithms for finding the true solution. In order to solve the discretized Poisson equation, Luna et al. (202T) accelerate the convergence of GMRES (Saad& Schultz 1986) with an initial guess that learned in real-time (i.e., as simulation code runs) with no prior data. [ et al. (2020) train network (incorporating differentiable phy physics, based on the underlying PDEs) in order to produce high-quality initial guesses for a CG solver. In somewhat similar vein, AAcmmann et al. (2020) use simple feedforward neural network to predict pointwise solution components, which accelerates the conjugate residual method used to solve a relatively simple shallow-water model (a more sophisticated network and loss function are needed to handle more general PDEs and larger-scale problems).
At least two papers (Ruelmann et al. 2018 Sappl et al. 2019) have sought to learn mapping between a matrix and an associated sparse approximate inverse. their investigation, Ruelmann pairs training data. Although
et al. (2018) propose training neural network using straightforward to implement, the cost of generating training data, let alone training the network, is prohibitive for large-scale 3D problems. (Sappl al. (2019) seek to learn a mapping between linear system matrices and sparse (banded) approximate inverses. Their loss function is the condition number the product of the system matrix and the approximate inverse; the minimum value of the condition number is one, which is achieved exactly when an exact inverse is obtained. Although this framework is quite simple, evaluating the condition number of matrix is asymptotically costly, and in general, the inverse of sparse matrix can be quite dense. Accordingly, the method is not efficient or accurate enough for the large-scale 3D problems that arise in real-world engineering problems.
Most relevant to the present work is FluidNet (Tomson et al. 2017). FluidNet develops highly- tailored CNN architecture that is used to predict the solution linear projection operation (specif- ically, for the discrete Poisson equation) given matrix and right-hand side. The authors demon- strate fluid simulations where the linear solve replaced by evaluating their network. Because their network is relatively lightweight and is only evaluated once per time step, their simulations run ef- ficiently. However, their design allows the network only one opportunity to reduce the residual for the linear solve; in practice, we observe that FluidNet is able reduce the residual by no more than about one order of magnitude. However, in computer graphics applications, at least four orders of magnitude in residual reduction usually required for visual fidelity, while in scientific and engi- neering applications, practitioners prefer solutions that reduce the residual by eight or more orders of magnitude (i.c.. to within machine precision). Accordingly, FluidNet's lack of convergence stands in stark contrast to classical, convergent methods like CG. Our method resolves this gap.
MOTIVATION: INCOMPRESSIBLEE FLOW
We demonstrate the efficacy of our approach with the linear systems that arise in incompressible flow applications. Specifically, we our algorithm to solve the discrete Poisson equations in
3

---

Under review as a conference paper at ICLR 2023
regular-grid briid discretization of the pressure projection equations that arise in Chorin's 's splitting technique (Chorin 1967) for the inviscid, incompressible Euler equations. These equations are
du p 01
du dx u
+ Vp - f vun=
where u is fluid velocity, p is pressure, P is density, and forr accounts for external forces like gravity. The equations assumed at all positions in the spatial fluid flow domain Q and for time & > 0. The left term in Equation enforces conservation of momentum in the absence of viscosity, and the second part enforces incompressibility and conservation of mass. These equations are subject to initial conditions p(æ;, 0) = p and u(æ 0) = u (a: as well as boundary conditions u(x t) (z) u 09 xx t) on the boundary of the domain æ an (where n is the unit outward pointing normal at position æ on the boundary). Equation 2 is discretized in both time and space.
du du 0 and body forces terms P al foutt and finally du
Du
Temporally, we split the advection
ØL 8x
enforce incompressibility via the pressure projection
Öt Vp- 0 such that v = 0; this is the standard advection- advection-projeetion scheme proposed by Chorin 1967). Using finite differences in time, can summarize this as
p0 u
ut
Ou dx
u"
fezt
(3)
At
1 Vp
Vp=+1
-
. u
(4)
p
1 Al
n2
(u ui n2.
(5)
For the spatial discretization, we use a regular marker-and-cell (MAC) grid (Harlow & Welch [965] with cubic voxels whereby velocity components are stored on the face of voxel cells, and scalar quantities (e.g.. pressure or density p) are stored at voxel centers. We use backward semi- Lagrangian advection (Fedkiw et al., 2001 Gagniere et al. 2020) for Equation 3 All spatial partial derivatives are approximated using finite differences. Equations and describe the pressure Pois- son equation with Neumann conditions on the boundary of the flow domain. We discretize the left hand side of Equation 4 using standard 7-point finite difference stencil. The right-hand side is discretized using the MAC grid discrete divergence finite difference stencils as well as contribu- tions from the boundary condition terms in Equation 3 We refer the reader to Bridson (2008) for more in-depth implementation details. Equation 5 is discretized by modifying the Poisson stencil to enforce Neumann boundary conditions. We do this using simple labeling of the voxels in the domain. For simplicity, we assume f2 C (0, 1) is a subset of the unit cube, potentially with internal boundaries. We label cells in the domain as either liquid or boundary. This simple classification is enough to define the Poisson discretizations (with appropriate Neumann boundary conditions at domain boundaries) that we focus on in the present work; we illustrate the details in Figure I We use the following notation to denote the discrete Poisson equations associated with Equations
d V-uu
(6)
where a is the discrete Poisson matrix associated with the voxelized domain, is the vector of unknown pressure, and b and 00 are the right-hand side terms from Equations 4 and 5 respectively. We define a special case of the matrix involved in this discretization be the Poisson matrix A (0.1)" associated with Ω = 1) i.c., a full fluid domain with no internal boundaries. We use this matrix for training, yet demonstrate that our network generalizes to all other matrices arising from more complicated flow domains.
4 DEEP CONJUGATE DIRECTION METHOD
We present our method for the deep learning acceleration of iterative approximations to the solution of linear systems of the form seen in Equation 6 We first discuss relevant details of the conjugate gradients (CG) method, particularly line search and A-orthogonal search directions. We then present deep learning technique for improving the quality of these search directions that ultimately reduces
4

---

Under review as a conference paper at ICLR 2023
iteration counts required achieve satisfactory residual reduction. Lastly, outline the training procedures for our deep convolutional neural network.
Our approach iteratively improves approximations to the solution æ of Equation 6 the method of CG, which requires the matrix A in Equation to be SPD. SPD matrices A give rise to the matrix norm | Vy A CG can be derived in terms of iterative line search improvement based on optimality this norm. That is, an iterate 2k-1 R x is updated in search direction dy by a step size ay that is chosen to minimize the matrix norm of the error:
We build on
0g arg min 1
r
2 || (0 + add)lla di Ad
where Tk-1 = b 2 x is the (A: - 1)th residual and b is the right-hand side in Equation 6 Different search directions result in different algorithms. A natural choice is the negative gradient of the matrix norm of the error (evaluated at the current iterate), since this will point in the direction of steepest decrease d& = - - | = r This is the gradient descent method (GD). Unfortunately, this approach requires many iterations in practice. A more effective strategy is to choose directions that are A d-orthogonal (i.e., d. A A directions form a basis for R" that the initial error can be written as eu x zu Σ.= e,d,, where are the components of the initial error written in the basis. Furthermore, when the search directions are A-orthogonal, the optimal step sizes (O) at each iteration satisfy
O for i + j). With this choice, the search
di. A 1 [E e,d, E a,d;) = Ck.
1 d2 ex df Ad
0
df Aad
df. Aad
That is, the optimal step sizes are chosen to precisely eliminate the components of the error on the basis defined by the search directions. Thus, convergence is determined by the (at most n) non-zero components e; in the initial error. Although rounding errors prevent this from happening exactly in practice, this property greatly reduces the number of required iterations (Golub & Loan 2012).
CG can be viewed as a modification of GD where the search direction is chosen as the component of the residual (equivalently, the negative gradient of the matrix norm of the error) that is A-orthogonal to all previous search directions:
d. Ara-1 df Add
d = =T--1
Σ da, ==1
fisk
In practice, his = 0 for i < k: 1, and this iteration can therefore be performed without the need to store all previous search directions d; and without the need for computing all previous has-
While the residual is a natural choice for generating A-orthogonal search directions (since it points in the direction of the steepest local decrease), is not the optimal search direction. If dg is parallel to (A - then Øs will be equal to æ since an (computed from Equation 7 will step directly the solution. We can see this by considering the residual and its relation to the search direction:
Tk A xz = b - 1 - A = rk-1 - A
In light of this, we use deep learning to create an approximation f(c, r) to (A where de- notes the network weights and biases. This is analogous to using a preconditioner in PCG; however, our network is not SPD (nor even a linear function). We simply use this data-driven approach as our means of generating better search directions dy. Furthermore, we only need to approximate a vector since the step size ag will account for any scaling in practice. In other words, parallel to 2) ( -'r, where the scalar s, is not defined globally; it only depends on r. and the model does not learn it. Lastly, with CG, we enforce A A-orthogonality, yielding search directions
fc, x s, ()
ME
his f(c, ((c,)R-1'=A' Add
d f(c,Tr-)))
hiddi
d. Add
We summarize approach in Algorithm Note that we introduce the variable istn To guarantee A A-orthogonality between all search directions, we must have 1sturt = 1. However, this requires storing all prior search directions, which can be costly. We found that using isan = l - 2 worked nearly as well as itan = 1 in practice (in terms our ability to iteratively reduce the residual of the system). We demonstrate this in Figure
5

---

Under review as a conference paper at ICLR 2023
5 MODEL ARCHITECTURE, DATASETS, AND TRAINING
Efficient performance of our method requires effective training of our deep convolutional network for weights and biases such that f(, R 8, (A (for arbitrary scalar sy)) We design a model architecture, loss function, and unsupervised training approach achieve this. Our approach has modest training requirements and allows for effective residual reduction while generalizing well problems not seen in the training data.
Algorithm 1 DCDM
1: rn = 200
2: k- =1
3: while || 1| 2 do
d₁ f(c,
-/-1
|| 1
5: for star
I i < k do
di d' h&s
6:
d. d'
5.1 Loss FUNCTION AND UNSUPERVISED LEARNING
7: dy-=haa.
8:
end for
Although generalize arbitrary matrices A from Equation 6 that correspond to domains Ω C (0, 1)a that have internal bound-
Ti ds
9:
a d: A
aries (see Figure D. train using just the matrix A (0.1)* from
10:
3y -b + dg
the full cube domain (0. 1) In contrast, other similar approaches
11:
Tk b A A
Tompson al. 2017 Yang et al. 2016) train using matrices A and right-hand sides b v.u" + bu" that arise from flow in many do- mains with internal boundaries. We train our network by minimiz-
12:
k k +
13: end while
ing the L difference ||r ar--aA (0.1)* f(c, 12 where a = f(cr)T A(0.)) f(c.r)
T" ((((c.))
from Equation 7 This
choice of a accounts for the unknown scaling in the approximation of f(e, r) to (A(0.1))) use an unsupervised approach and train the model by minimizing
r. We
Loss(f, c, D) = bbj ErED ||| f(e,r) A f(a.r)
r ((.rr)
0.1) f(e,r))l
for given dataset D consisting of training vectors b². In Algorithm the normalized residuals T are passed as inputs to the model. Unlike in e.g. FluidNet (Tompson al., (017). only the first residual TPPTI is directly related to the problem-dependent original right-hand side b. Hence we consider a broader range of training vectors than those expected in given problem of interest, e.g., incompressible flows. We observe that generally the residuals r'k in Algorithm|] [ are skewed to the lower end of the spectrum of the matrix A Since A is discretized elliptic operator, lower end modes are of lower frequency of spatial oscillation. We create our training vectors b² E D using < n approximate eigenvectors of the training matrix A (0.1)" We use the Rayleigh-Rin method to create approximate cigenvectors q. < < m. This approach allows us to effectively approximate the full spectrum of A (0.1)" without computing the full eigendecomposition. which can be expensive (O(n²) at high resolution. We found that using m 10000 worked well in practice. The Rayleigh- vectors are orthonormal and satisfy Q A (0.1)² Qm = Am, where Am diagonal matrix with nondecreasing diagonal entries x; referred to Ritz values (approximate > 6/9; where the
eigenvalues) and Qm [go: qı.... 4m-1] _axm We pick b²
== 5a
coefficients of are picked from standard normal distribution
9 10, i m +
- - ()
otherwise
where O is a small number (we used 0 500), and is the first index that A3 O. This choice creates 90% of b' from the lower end of the spectrum, with the remaining 10% from the higher end. The Riemann-Lebesgue Lemma states the Fourier spectrum of continuous function will decay at infinity, so this specific choice of b, is reasonable for the training set. In practice, we also observed that the right-hand sides of the pressure system that arose flow problems (in the empty domain) tended to be at the lower end of the spectrum. Notably, even though this dataset only uses Rayleigh- Ritz vectors from the training matrix A (0.1)", our network can be effectively generalized to flows in irregular domains, e.g., smoke flow past a rotating box and flow past a bunny (see Figure
We generate the Rayleigh-Ritz vectors by first tridiagonalizing the training matrix A 10. with m Lanczos iterations [anczoos 1950) form T Q4 T A (0.1) QL R We then diagonalize

---

Under review as a conference paper at ICLR 2023
T = =QTA,, Am Q. While costly, we note that this algorithm is performed on the comparably small m m matrix T (rather than on the A A E [_*x). This yields the Rayleigh-Rinn vectors as the columns of Qm Q Q. The Lanczos vectors are the columns of the matrix Q and satisfy a three-term recurrence whereby the next Lanczos vector can be computed from the previous two as
B - A (0.1)) of 3 agg/.
where jj and By are diagonal and subdiagonal entries of T 33 is computed so that af+1 is unit vector, and 0j++1 = qj*+ A (0.1)" 4i+1- We initialize the iteration with a random qu E span( A (0.1))) ) The Lanezos algorithm can be viewed as a modified Gram-Schmine technique to create an orthonor- mal basis for the Krylov space associated with q and A (0.1), and it therefore suffers from rounding error sensitivities manifested as loss of orthonormality with vectors that do not appear in the recur- rence. We found that the simple strategy described in Paige (1971) of orthogonalizing each iterate with all previous Lanczos vectors to be sufficient for our training purposes. Dataset creation takes 5-7 hours for a 64 computational grid, and 2-2.5 days for a 128 grid.
5.2
MODEL ARCHITECTURE
The internal structure of our CNN ar- chitecture for 1283 grid is shown in Figure 2 It consists of a series of convolutional layers with residual connections. The upper left of Fig- ure [ Residual Blocks) shows use of multiple blocks of residually connected layers. Notably, within each block, the first layer directly af- fects the last layer with an addition operator.Al All non-input or output con- volutions use a 3 3 3 filter, and all layers consist of 16 feature maps. In the middle of the first level, layer is downsampled (via the average pool- ing operator with (2 x 2 x 2) pool size) and another set of convolutional layers is applied with residual connection blocks. The last layer in the second level is upscaled and added to the layer that is downsampled. The last layer in the network is dense with a linear activation function. The activation functions in all convolutional layers ReLU, except for the first convolution, which uses linear activation function.
Residual (K-RB)
16
10x 122?
x 128
16x 128
16x 128?
16x 122?
3xcla3 Cony, RuLU
Residual Block
Addition
3o2icc Conv, Linwar
zicet Conv Layars
Dennee, Linear
2-B8
S-BB
16xx 120?
x 1227
16x 128"
16
x 128
x6xx 128"
2x2x2 AveraggPPolling
f-
UpScall
x 64
18 64
Figure 2: Model architecture for training with A (0,1)* on a 128 grid.
Initially tried simple deep feedforward convolutional network with residual connections (mo- tivated by He et al. (2016)) Although such a simple model works well for DCDM, it requires high number of layers, which results in higher training and inference times. We found that creating par- allel layers of CNNs with downsampling reduced the number of layers required. In summary, our goal was to first identify the simplest network architecture that provided adequate accuracy for our target problems, and subsequently. we sought to make architectural changes to minimize training and inference time; further optimizations are possible.
Differing resolutions use differing numbers of convolutions, but the fundamental structure remains the same. More precisely, the number of residual connections is changed for different resolutions. For example, a 64 grid uses one residual block on the left, two on the right on the upper level, and three on the lower level. Furthermore, the weights trained on a lower resolution grid can be used effectively with higher resolutions. Figure 4 shows convergence results for 256 grid, using model trained for a 64 grid and 128 grid. The model that we use for 256 grids in our final examples was trained a 128 grid; however, the shown in the figure, even training with a 64 grid allows for efficient residual reduction. Table [[] shows results for three different resolutions, where DCDM uses 6 and 3 trained models. This approach makes the number of parameters in the model independent of the spatial fidelity of the problem.
7

---

Under review as a conference paper at ICLR 2023
5.3 TRAINING
Using the procedure explained in Section 5.T we create the training dataset D E span(A (O.)))))) n Sn of size 20000 generated from 10000 Rayleigh-Ritz vectors. We train our model with TensorFlow (Abadi al. 2015) on a single NVIDIA RTX A6000 GPU with 48GB memory. Training is done with standard deep learning techniques-ma more precisely, back-propagation and the ADAM op- timizer (Kingma & Bal 2015) (with starting learning rate 0.0001). Train- ing takes approximately 10 minutes and hour per epoch for grid res- olutions 64 and 128", respectively. We trained model for 50 epochs; however, the model from the thirty- first epoch was optimal for 64 and the model from the third epoch was optimal for 128*.
Figure DCDM for simulating a variety of incom- pressible flow examples. Left: smoke plume at 6.67, 13.33, 20 seconds. Middle: smoke passing bunny at 5, 10, 15 seconds. Right: smoke passing spin- ning box (time-dependent Neumann boundary condi- tions) at 2.67, 6, 9.33 seconds.
6 RESULTS AND ANALYSIS
We demonstrate DCDM on three increasingly difficult examples and provide numerical evidence for the efficient convergence of our method. All examples were run on a workstation with dual AMD EPYC 75F3 processors and 512GB RAM.
Figure showcases DCDM for incompressible smoke simulations. In each simulation, inlet bound- ary conditions set in circular portion of the bottom of the cubic domain, whereby smoke flows around potential obstacles and fills the domain. We show a smoke plume (no obstacles), flow past a complex static geometry (the Stanford bunny), and flow past a dynamic geometry rotating cube). Visually plausible and highly-detailed results are achieved for each simulation (see supplementary material for larger videos). The plume example uses computational grid with resolution 1283. while the other two uses grids with resolution 256 (representing over 16 million unknowns). For each linear solve, DCDM was run until the residual was reduced by four orders of magnitude.
For the bunny example, Figures -b-b demonstrate how residuals decrease over the course of linear solve, comparing DCDM with other methods. Figure 4. shows the mean results (with standard deviations) the course of 400 simulation frames, while in Figure 4b. we illustrate behavior on a particular frame (frame 150). For FluidNet, we use the implementation provided by [Buidnetsc222 (2022). This implementation both subfigures,
Table 1: Timing and iteration comparison for different methods on the bunny example. DCDM-(64, 128) calls model whose parameters trained over {643, 128 } grid. All computations are done using only CPUs; model inference does not use GPUs
643 Grid
1283
Grid
2563 Grid
Method
t, n1.
f,r
f,
n,
DCDM-64
2.71s
16
22s
27
261s
58
DCDM-128
5.37s
26s
24
267s
44
CG
1.77s
168
26s
465
1548s
1046
Deflated PCG
771.6s
117
3700s
277
21030s
489
includes pre-trained models that we use without modification. that the FluidNet residual never changes, since the method is not iterative: FluidNet reduces the initial residual by no more than one order magnitude. On the other hand, with DCDM, we can continually reduce the residual (e.g., by four orders of magnitude) our method, just as with classical CG. In Figure/pp. we also visualize the convergence three other classical methods, CG, Deflated PCG (Saad et al. 22000, and incomplete Cholesky preconditioned CG (ICPCG)); clearly, DCDM reduces the residual in the fewest number iterations (c.g., approximately order of magnitude fewer iterations than ICPCG). Since FluidNet not
evident
we apply more iterations of
8

---

Under review as a conference paper at ICLR 2023
DCDM 128
ICPOG
64 trained trained
FN 128
DCDM 256
1o (||=||//l ll/llll)))
EN 256
DCDM DeflatedPCG
2
-3
-4
-4 400
-4
20 40 200
600
10 20 30
0 20 40
iterations
Figure 4: Convergence data for the bunny example. (a) Mean and std. dev. (over all 400 frames in the simulation) of residual reduction during linear solves (with 128 and 256 grids) using FluidNet (FN) and DCDM. (b) Residual plots with CG, ICPCG, Deflated PCG, FN, and DCDM at frame 150. Dashed and solid lines represent results for 128 and 256 respectively. (c) Decrease in residuals with varying degrees of A-orthogonalization (i. = isan). (d) Reduction in residuals when the network is trained with 643 or 128 grid for the 256 grid simulation shown in Figure 3 Middle.
iterative method and lacks a notion of residual reduction, we treat To for FluidNet as though an initial guess of zero is used (as is done in our solver).
To clarify these results, [able|]]eportt convergence statistics for DCDM compared standard iter- ative techniques CG and Deflated PCG. For all 64 1283, and 256 grids with the bunny example, we measure the time 4,- and the number of iterations n, required to reduce the initial residual on particular time step of the simulation by four orders of magnitude. DCDM achieves the desired results in by far the fewest number of iterations at all resolutions. At 256", DCDM performs ap- proximately 6 times faster than CG, suggesting potentially even wider performance advantage at higher resolutions. Inference the dominant cost in an iteration of DCDM: the other linear algebra computations in an iteration of DCDM comparable to those in CG. The nice result of our method is that despite the increased time per iteration, the number of required iterations is reduced so drasti- cally that DCDM materially outperforms classical methods like CG. Although ICPCG successfully reduces number of iterations 4 b. we found the runtime to scale prohibitively with grid resolution, so we exclude it from comparison in table Notably, even though Deflated PCG and DCDM are based on approximate Ritz vectors, DCDM performs far better, indicating the value of using a neural network.
7 CONCLUSIONS
We presented DCDM, incorporating CNNs into CG-style algorithm that yields efficient, conver- gent behavior for solving linear systems. Our method is evaluated on linear systems with over 16 million degrees of freedom and converges thermore, despite training the underlying network on domains without obstacles, our network is able successfully predict search directions that enable efficient linear solves on domains with complex and dynamic geometries. Moreover, the training data for our network does not require running fluid simulations or solving linear systems ahead of time; our Rayleigh- -ittighhhttttttt vector approach enables us to quickly generate very large training datasets, unlike approaches seen in other works.
desired tolerance merely tens of iterations. Fur-
Our network was designed for and trained exclusively using data related to the discrete Poisson matrix, which likely limits the generalizability of our present model. However, we believe our method is readily applicable to other classes PDEs (or general problems with graph structure) that give rise to large, sparse, symmetric linear systems. We note that our method is unlikely to work well for matrices that have high computational cost to evaluate A x (such as dense matrices), since training relies efficient A z evaluations. An interesting question to consider is how well our method and current models would apply to discrete Poisson matrices arising from non-umt uniform grids, e.g., quadtrees octrees (Losasso al. 2004).

---

Under review as a conference paper at ICLR 2023
REFERENCES
M. Abadi, A. Agarwal, P. Barham, E. Brevdo, 7. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefow- icz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Va- sudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng. TensorFlow: Large-scale machine learning on heterogeneous systems, 2015. URL https: Software available from tensorflow.org.
// //www. twnsort low. org/
J. Ackmann, P. D. Düben, T. N. Palmer, and P. K. Smolarkiewicz. Machine-leamed preconditioners for linear solvers geophysical fluid flows. arXiv preprint arXiv:2010.022800 2020.
A. Brandt. Multi-level adaptive solutions to boundary-value problems. Math Comp, 31(18))3333- 390, 1977.
R. Bridson. Fluid simulation for computer graphics. Taylor & Francis, 2008
A. Chorin. A numerical method for solving incompressible viscous flow problems. J Comp Phys, 2 (1):12-26,, 1967.
R. Fedkiw, Stam, and H. Jensen. Visual simulation of smoke. In SIGGRAPH, pp. 15-22, ACM, 2001.
fluidnetsc22. Huiddettsc2/huiinnetss222 v0.0.1, April 2022. URL https://doi.org// 5281 6424901]
doi: 10.5228//zemodo.64429000 URL: ntttps://diiorr//10
6424901
S. Gagniere, D. Hyde, A. Marquez-Razon, C. Jiang, Z. Ge, X. Han, Q. Guo, and J. Teran. A hybrid Lagrangian/Eulerian collocated velocity advection and projection method for fluid simulation. Computer Graphics Forum, 39(8):1-14, 2020. doi: https://oi..rr/// wot.//..1/1111 1//g 14096.
G. Golub and C. Van Loan. Matrix computations, volume 3. JHU Press, 2012.
M. Götz and H. Anzt. Machine learning-aided numerical linear algebra: Convolutional neural net- works for the efficient preconditioner generation. 2018 IEEF/ACM 9th Workshop on Latest Advances in Scalable Algorithmss for Saarge-Scale Systems (scalA), pp. 49-56. IEEE, 2018.
A. Grebhahn, N. Siegmund, H. Köstler, and S. Apel. Performance prediction of multigrid-solver configurations. In Software for Exascale Computing- SPPEEE 2013-2015, pp. 69-88. Springer, 2016.
D. Greenfeld, M. Galun, R. Basri, I. Yavneh, and R. Kimmel. Learning to optimize multigrid PDE solvers. In Int Conf Mach Learn, pp. 2415-2423, PMLR, 2019.
Harlow and E. Welch. Numerical calculation of time dependent viscous flow of fluid with a free surface. Phys Fluid, 812)):182-2188, 1965.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recog- nition. In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp- 770-778, 2016. doi: 10. 10.1109/CVP..2016.900
M. R. Hestenes and E. Stiefel. Methods of conjugate gradients for solving linear systems. Journal of research of the National Bureau of Standards, 49(6):409, 1952.
J. Hsich, S. Zhao, S. Eismann, L. Mirabella, and S. Ermon. Learning neural PDE solvers with convergence guarantees, 2019. URL //arxiv....iivoo//o wrw. abs /9066 01200
T. Ichimura, K. Fujita, M. Hori, L. Maddegedara, N. Ueda, and Y. Kikuchi. A fast scalable iterative implicit solver with Green's 's function- neural networks. In 2020 JEEE/ACM 11th Workshop Latest Advances in Scalable Algorithms for Large-Scale Systems (ScalA), pp. 61-68, 2020. doi: 10. 109/ScalA519936622000
D. Kershaw. The incomplete Cholesky conjugate gradient method for the iterative solution of sys- tems of linear equations. J Comp Phys, (6(1):43-66 1978.
10

---

Under review as a conference paper at ICLR 2023
D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. CoRR, abs/1412.6980, 2015.
C. Lanczos. An iteration method for the solution of the eigenvalue problem of linear differential and integral operators. 1950.
F. Losasso, F. Gibou, and R. Fedkiw. Simulating water and smoke with an octree data structure. ACM Trans. Graph., 23(3):457 462, 2004.
K. Luna, K. Klymko, and P. Blaschke. Accelerating GMRES with deep learning real-time, 2021. URLpttps: arxiv. org/ has///b/////// 2tt.////13.11999
I. Luz, M. Galun, H. Maron, R. Basri, and I. Yavneh. Learning algebraic multigrid using graph neural networks. In Int Conf Mach Learn, pp. 6489-6499. PMLR, 2020.
H. Nussbaumer, The fast Fourier transform. In Fast Fourier Transform and Convolution Algorithms, pp. 80-11011111111111111111111111111 Springer, 1981.
C. C. Paige. The computation of eigenvalues and eigenvectors of very large sparse matrices. Phl) thesis, University of London, 1971.
C. C. Paige and M. A. Saunders. Solution of sparse indefinite systems of linear equations. SIAM journal on numerical analysis, 12(4):617-629, 1975.
H. Ruelmann, M. Geveler, and S. Turek. On the prospects of using machine learning for the numer- ical simulation of PDEs: Training neural networks assemble approximate inverses, 2018. URL lttpp:// /xxx doi. /tg//0./.... 21777///7220/ 18778
Y. Saad. Iterative Methods for Sparse Linear Systems. Society for Industrial and Applied Mathe- matics, USA, 2nd edition, 2003. ISBN 0898715342.
Y. Saad and M. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsym- metric linear systems. SIAM J Sci Stat Comp, 7(3):856 869, 1986.
Y. Saad, M. Yeung, J. Erhel, and Guyomare'h, hh A deflated version of the conjugate gradient algorithm. SIAM Journal on Scientific Computing, 21:1909-1926, 2000.
J. Sappl, I.. Seiler, M. Harders, and W. Rauch. Deep learning of preconditioners for conjugate gradient solvers in urban water related problems, 2019. URI. http::://rrxiv.org/abbbb 19C6.06923
R. Stanaityte. ILU and Machine Learning Based Preconditioning For The Discretized Incompress- ible Navier-Stokes Equations. PhD thesis, University of Houston, 2020.
E. Stiefel. Über einige methoden der relaxationsrechnung. Zeitschrifi für angewandte Mathematik und Physik ZAMP, 3(1):1-33, 1952.
Tompson, K. Schlachter, P. Sprechmann, and Perlin. Accelerating Eulerian fluid simulation with convolutional networks. In D. Precup and Y. Teh (eds.), Proc 34th Int Conf Mach Learn, volume 70 Proc Mach Learn Res, pp. 3424-3433, PMLR, 06-11 I Aug 2017.
L. Trefethen and D. Bau. Numerical Linear Algebra, volume 50. SIAM, 1997.
K. Um, R. Brand, Y. Fei, P. Holl, and N. Thuerey. Solver-in-the-loop: Learning from differentiable physics interact with iterative PDE-solvers. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (eds.), Advances in Neural Information 6111-6122. Curran Associates, Inc., 2020. URL http:::::ss........ http:/tpss///proceeeeeo///oo mew.eeeddddnnsssnnnusssee..... cc/
PrciiiionnPPrrrrrrr Processing Systems, volume 33, pp.
C. Yang, X. Yang, and X. Xiao. Data-driven projection method in fluid simulation. Comp Anim Virt Worlds, 2153-4):4154444 2016.
11


Reviews:
[{"review_id": "sLOGfj7yh9K", "reviewer": "Reviewer_gFH2", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors propose a method for \"learning\" a pre-conditioner network for the conjugate gradient method as it applies to linear systems Ax = b resulting from the discrete Poisson equations.\nThe authors propose to train a CNN model to approximate the inverse of a given matrix A (modulo scaling) by using training data generated from linear combinations of eigenvectors of A which need to be computed first.\nUsing the proposed method in the context of a conjugate gradient solver shows a significant reduction in the number of iterations required to reach a certain residual norm for their specific test application (a wall-clock time reduction is indicated for larger grid sizes).", "summary_of_the_review": "I do not recommend to accept the paper in its current form since the specific application does not necessarily seem relevant to the audience of ICLR and otherwise the analysis of the approach is still a little bit too superficial to warrant publication. The paper could be improved by extension to other types of matrices, a more in depth analysis of the presented method, e.g. how and what the pre-conditioner network learns based on e.g. what parts of the eigenspectrum it sees, a better explanation of the total computational costs and also computational complexities involved in the compared approaches and perhaps a more comprehensive treatment of baselines.", "strength_and_weaknesses": "Strengths:\n- The paper does a good job explaining the specific problem setup.\n- The paper indicates that there is a possible generalization effect by scaling the grid size for a given problem instance and this could be explored further.\n\nWeaknesses:\n- The paper addresses a very specific problem type and is not as general as the title suggests.\n- For someone not familiar with the literature on approximating matrix inverses the paper could improve on motivation for the network architecture, for example why are non-linearities even necessary / helful if the learned function is linear.\n- The paper does not contain a theoretical analysis of the computational complexity of obtaining eigenvalue and matrix inverse approximations. To me it is not clear why a method that uses stochastic gradients (Adam) should converge to a good matrix inverse approximation faster than more specific deterministic approaches, especially if you already have a sample of eigenvectors available.\n- The experimental setup is not always well motivated and seems somewhat arbitraty, for example it seems to me that the choice of m=10k example eigenvectors should be somehow related to the eigenspectrum.\n- The paper does not consider the computational cost of the whole end-to-end setup. Is the method useful if only a single solve with the matrix A is performed given how expensive it is to generate the training set and train the network? If we require multiple solves with the e.g. similar matrices A, then how many solves are typical for a given PDE solver and how useful can this method be? What residual norm for the linear system is required in order for the PDE solver to converge?\n- The paper does not show how well the network model actually learns / generalizes on a validation set and does not compare different architecture options in that way.\n- I'm not sure whether the FluidNet baseline is used in the correct way (just using a model pre-trained on a different matrix does not the intended use)\n\n**Suggestions and questions:**\n\nAbstract:\n- \"require many millions of unknowns\" seems fuzzy to me, I would use different wording that \"require\" here\n- in general the abstract could be made more concise\n- mention some more relevant quantitative data points in the abstract such as the actual computational cost involved in training the inverse model (not just iteration count)\n- explain why a neural network approach is helpful when the inverse is a linear function (or is it a linear neural network)?\n- explain the connection between SGD updates and Krylov subspace updates and the information complexity involved (why can the data be more efficient than the Krylov updates)\n- obvious necessary improvement would be to consider different sources of linear systems (or otherwise make the title of the paper more specific and pick a different target publication venue)\n\nIntroduction:\n- \"since this decreases the degree of nonlinearity\" what is meant by nonlinearity here?\n- \"We use unsupervised learnin\" how is this unsupervised learning? There is \"supervision\" from A, right?\n- Do they account for changing A in every iteration of a PDE solver? Or is the approach mostly suitable for repeatedly solving with the same A?\n- How does the approach really generalize to unseen A?\n- Learning or computing a perfect inverse for a given A would in principle solve the problem perfectly, but that is not a very useful insight.\n\nRelated work:\n- \"However, in computer graphics applications, at least four orders of magnitude [...]\" would be good to have references for these claims.\n- The figure is not really related to the method explained in the paper. I would rather use an illustration the demonstrates the workflow or method of the paper.\n\nMotivation:\n- For the paper targeting an ML conference I would say providing the PDE and math related to the problem is not really that helpful. I would rather focus on explaining the resulting structure of the matrix. But explaining just the discretization scheme does not tell the reader a lot about how the matrix will look like.\n\nMethod:\n- Is there going to be an explanation on the effect of their pre-conditioner approach on the eigenspectrum of the matrix (as this reflects conditioning / convergence properties of CG)?\n- Why not a linear function for the pre-conditioner network?\n- Explain where the network is trained in relation to the PDE solver earlier in the paper\n\nArchitecture:\n- \"Our approach has modest training requirements\" what does that mean quantitatively? Also the rest of the sentence, these statements are too vague to be meaningful.\n- Again, need to discuss whether a network trained for a given matrix for a given set of boundary / initial conditions is in any way useful to another problem\n- Again, not sure why the approach is called unsupervised, the Matrix vector product with a given A is a form of supervision\n- Explain earlier and clearly where the \"training input vectors\" come from for you network training\n- Isn't getting approximate eigenvector for A just as expensive as solving the system? When reading that these are used I can imagine that they will be very helpful in approximating and inverse because finding them is essentially as hard as finding the inverse.\n- What does it mean to \"effectively approximate the full spectrum\"?\n- There should be much more focus on the computational cost and tradeoff of these aspects of how to compute the eigenvectors etc., which ones to compute (computing them for different parts of the eigenspectrum can have different difficulty, depending on the eigenvalue distriubtion).\n- Does m=10k work equally well for all different kinds of n? That seems somewhat implausible? m should probably depend on the discretization? Are the experiments run on m as hyperparam presented in the work? Is the cost of computing m eigenvectors part of the computational cost comparison? What type of pre-conditioner could be constructed with m=10k eigenvectors explicitly?\n- What is the motivation for the 9 constant in the definition of c_j^i? Where does theta=500 come from?\n- Need to discuss exactly what is the cost of the m Lanczos iterations? That seems like information that would typically be usable to kick start the system solve for A as well.\n- The concrete numbers used for the grid sizes and so on points to the fact that this paper is quite specific to a given application and should possibly be titled as such\n\nModel Architecture:\n- Would be helpful to get some intuition on why convolutions may be suitable for the type of A. A guess would be that since PDE solutions are smooth and neighboring entries in A related to spacial neighborhood in the discretization grid somehow the convoultion helps with generating smoother pre-conditioners in that same sense? What does the inverse of A look like in terms of \"smoothness\" structure? Does convolution hurt sparsity?\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "n_wAAJPi6j", "reviewer": "Reviewer_nmQe", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors propose a machine-learning method to accelerate conjugate descent by training a convolutional network to regress good descent directions. The central idea---to obtain conjugate directions via self-supervised learning---is simple, creative, and it works very well. The authors effectively combine a number of techniques from numerical analysis and numerical linear algebra. They apply their method to build a fast Poisson solver and show that it consistently outperforms principled solvers as well as recent ML-accelerated solvers on realistic benchmarks.\n", "summary_of_the_review": "As stated above, I think this is a good paper. It addresses an important problem, finds an effective way to use machine learning (via self-supervision), and does exactly what it says in the introduction, without overselling or underselling. I would be happy to see this paper accepted.", "strength_and_weaknesses": "The idea proposed by the authors is simple, effective, and by my standards the right way to use neural networks to accelerate forward solves. Instead of attempting to directly regress the solution or trying complicated things, which more or less never give reconstructions that are good for applications (and tends to horribly overfit training data), they opt to accelerate an iterative scheme in a new way, by regressing a new descent direction. Although there exist some related ideas, the whole package with the self-supervised training scheme and a well thought-out training set construction is quite novel as far as I can tell.\n\nI further find the idea to be very well executed. The authors clearly have a strong understanding of numerical analysis and numerical linear algebra. The literature overview and the problem description are very clear. The experimental validation is convincing and the reproducible code is provided in the supplementary material. I do have a couple of questions to clarify below but overall my impression is very positive.\n\n\n- You write that \"We only ask that our network approximate the inverse up to an unknown scaling since this decreases the degree of nonlinearity and since it does not affect the quality of the search direction (which is scale independent).\" My intuition is that the degree of nonlinearity in the network which computes the descent direction is not important. Rather, the norm of the vector computed by the network would almost certainly not be the optimal step length and so it makes sense to compute it exactly as you do in your loss. (That is, even if the network was trained to regress something with a specific norm it would be opportune to rescale it; one might try training like that and rescaling only at test time.)\n\n- I don't think that your training strategy is unsupervised. Rather, the supervision is performed through the operator $A^\\Omega$. A commonly used term in this situation where you have the residual both as input to f and as the target is \"self-supervised learning\".\n\n- The abstract suggests that the number of iterations required by your method is independent of the problem size, but Table 1 indicates otherwise. Could you clarify this?\n\n- I find it a bit surprising that the model from the third epoch was optimal at resolution 128^3. Is \"optimal\" meant in the training loss sense or in the sense that it performed best at test time? If the latter is the case, does the loss keep decreasing after the third epoch? And if not, can you speculate why?\n\n- What exactly is meant by the \"data dependency\" of incomplete Cholesky preconditioners on page 2? Is this alleaviated by the presented method (I am wondering about the word \"However, \")\n\n\n### Minor\n\n- please add callouts (a, b, c, d) in Figure 4 even though they seem obvious\n- please use the same color for the same algorithm in all subfigures, otherwise it's quite confusing (you might also consider other ways to vary the line style as well as using a color scheme that is friendly to the color blind---red / green is quite bad)\n\n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "n2cyZVQuI4", "reviewer": "Reviewer_GZHi", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposed a data-driven modification of the CG algorithm for the solution of SPD linear systems. The proposed algorithm updates the approximate solution using a forward pass of a neural network which is trained with some of the lowest eigenmodes of \na discretization of the elliptic operator without internal boundary conditions. Numerical experiments demonstrate the effectiveness of the proposed scheme in terms of number of iterations and solution runtime. ", "summary_of_the_review": "While I like the topic and the initiative by the authors, and I would like to give a higher score, the proposed method is simply impractical. Unless the authors can show that the timings of Lanczos+Training+Solution is lower than that of ICPCG, there is no reason for a scheme as the one proposed in this paper. Moreover, the memory requirements seem to be enormous compared to ICPCG. Not everyone has access to a 512GB RAM machine. To the latter, please add the need to store CG vectors to enforce A-orthogonality (which can be expensive in terms of both memory and computations).", "strength_and_weaknesses": "Strengths: \n1) The authors present an interesting data-driven approach to solve elliptic SPD linear systems. To the best of my knowledge, the method is novel. The training stage is described sufficiently well.\n2) The experiments suggest that the scheme works.\n3) The paper is clear and well-written.\n4) The supplemental material is very good (codes+video).\n\nWeaknesses: \n1) The proposed scheme appears to be impractical. Even with the best case settings, training takes from 20 minutes to 120 minutes (for the two smaller grids). These timings are nowhere on Table 1. What about the 256 grid? Also, do these timings include the Lanczos phase or only the training phase after the eigenmodes are computed? I suspect Lanczos is not included but I would like the authors to clarify this.\n\n2) In Table 1, the authors avoid to compare the solution phase of their method against ICPCG, which is nonetheless included later on in Figure 4, but this time they do not report timing results. Please report timings for ICPCG as well. Also, since the problem is elliptic, AMG (algebraic multigrid) is probably a better choice to compare against as this is the state-of-the-art method currently used. Please add this option in addition to Incomplete Cholesky. Based on the timings the authors report, the proposed scheme is slower even compared to plain CG, let alone state-of-the-art AMG preconditioned CG. Indeed, even if I count only two epochs, 60 minutes each, DCDM requires 7200+22 seconds versus 24 seconds for CG. Why is DCDM useful? Because the solution phase is faster? This portion of DCDM reported seems to be tiny compared to the overall running time of the scheme if I count Lanczos+Training. Am I missing something?\n\n3) Even if DCDM was faster with the training time included, it is still extremely impractical as it requires the computation of a huge number of eigenmodes just for the training phase. What is the wall-clock time required to perform 10000 Lanczos iterations? Is this time included in the results shown?\nThis is a /*huge*/ amount of Lanczos vectors required (and eigenmodes computed). The amount of memory required to save the Lanczos vectors is several orders of magnitude more than that of ICPCG which requires only storage for about six vectors plus the preconditoner (since the authors use Incomplete Cholesky this is pretty low). \n\n4) Please report all timings in the revised version.\n\nMinor comment: DeflatedPCG should be renamed as DeflatedCG, as the former indicates that deflation is applied on an already preconditioned variant of CG, which is not the case here. ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}, {"review_id": "7p7Hvpls6o", "reviewer": "Reviewer_2Pr7", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "In this paper, the authors consider the general problem of solving a linear system of equation in the case that the matrix is positive semi-definite (and sparse). The authors study a data-driven variation of conjugate gradient to solve the linear system. The algorithm (DCDM) is an iterative algorithm inspired by conjugate gradient algorithm where the descent direction is computed using a trained convolutional neural network. In the paper, the authors specialize their presentation of the algorithm for the incompressible flow problem. Comparison of performance DCDM, CG and preconditioned CG is provided. ", "summary_of_the_review": "I believe the main idea in the paper is interesting however the paper lacks clarity and the contributions of the work is not clear. An analysis of convergence of the algorithm is also needed. For these reasons, I believe the paper is marginally below the threshold for acceptance.", "strength_and_weaknesses": "Below are my opinions that would strengthen the paper:\n\n- The will benefit from a clearly outlined contribution section. \n\n- Terms in the paper are used before explaining what they mean (eg., A-conjugate (in page 2), A-orthogonal (in page 4)). The paper would benefit from a notation section. What does $(0,1)^3$ mean in the paper?\n\n- Is the number of iterations required to find a solution using DCDM bounded (similar to CG and Krylov subspace methods)? An analysis on the convergence of DCDG is lacking. What conditions on the trained network would be needed for such an analysis to be feasible?\n\n- The linear system considered (please correct me if I am wrong) in the paper requires the matrix $A$ to have entries in the interval $(0,1)$. Is the algorithm limited to such matrices or can it be generalized larger intervals? If so, why?", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "7hZG1m4igW", "reviewer": "Reviewer_GXf5", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper presents an iterative method to solve large sparse pos. def. linear systems of equations with a deep learning technique to select step sizes. The authors demonstrate the efficacy of their approach on linear systems that arise from incompressible flow applications, which are after standard Chorin's splitting end up being shifted discrete Poisson equations. The paper uses deep learning to improve the quality of the search directions that appear in the conjugate gradient method. ", "summary_of_the_review": "The authors treat the search direction in CG as a kind of hyperparameter that can be selected by a deep learning approach. This means that there is no way to mathematical understanding of this iterative method. This makes it impossible to design good preconditioners so the approach is unlikely to generalize or be understanding from a mathematical perspective. ", "strength_and_weaknesses": "The most important thing about iterative methods such as conjugate gradient (CG) for discrete Poisson equations is preconditioning.  By making the step direction selected from a deep learning approach, there is no longer any theory. In particular, there is no underlying polynomial approximation problem that CG is solving, i.e., min_{p in P_k, p(0) = 1} max_{lambda = eig of A}|p(lambda)|. Therefore, one cannot get any idea of how to successfully precondition A, which for this iterative method to be practical on a range of problems will be necessary. For the discrete Poisson equations, the authors show that one can improve the selection of search direction after doing lots of matrix-vector products. \n\nI suspect that the authors have ignored the tensor-product structure in the matrix A^{(0,1)^3}.\n\nThe network is designed for and trained exclusively using data related to a single discrete Poisson matrix, which likely limits the generalizability of the present model. The multigrid method is a more general and extremely fast solver for discrete Poisson equations. ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject