aa r X i v : . [ m a t h . A P ] A ug PDES WITH COMPRESSED SOLUTIONS
RUSSEL E. CAFLISCH, STANLEY J. OSHER, HAYDEN SCHAEFFER, AND GIANG TRAN
Abstract.
Sparsity plays a central role in recent developments in signal processing, linear algebra, statistics,optimization, and other fields. In these developments, sparsity is promoted through the addition of an L norm (or related quantity) as a constraint or penalty in a variational principle. We apply this approachto partial differential equations that come from a variational quantity, either by minimization (to obtainan elliptic PDE) or by gradient flow (to obtain a parabolic PDE). Also, we show that some PDEs can berewritten in an L form, such as the divisible sandpile problem and signum-Gordon. Addition of an L term in the variational principle leads to a modified PDE where a subgradient term appears. It is knownthat modified PDEs of this form will often have solutions with compact support, which corresponds to thediscrete solution being sparse. We show that this is advantageous numerically through the use of efficientalgorithms for solving L based problems. Introduction
Sparsity has played a central role in recent developments in fields such as signal processing, linear algebra,statistics and optimization. Examples include compressed sensing [12, 17], matrix rank minimization [31],phase retrieval [10] and robust principal component analysis [11, 16, 30], as well as many others. A keystep in these examples is the use of an L norm (or related quantity) as a constraint or penalty term ina variational formulation. In all of these examples, sparsity is for the coefficients ( i.e. , only a small set ofcoefficients are nonzero) in a well-chosen set of modes for representation of the corresponding vectors orfunctions.The use of sparse techniques in physical sciences and partial differential equations (PDEs) has beenlimited, but recent results have included numerical solutions of PDEs with multiscale oscillatory solutions[32], efficient material models derived from quantum mechanics calculations [26], “compressed modes” forvariational problems in mathematics and physics [27], and “compressed plane waves” [28]. In the latter twoexamples, sparsity is used in a new way, in that the solutions are sparse and localized in space (as opposed tosparsity of the coefficients in some modal representation). Sparse solutions with respect to low-rank librariesare used in modeling and approximating dynamical systems, see for example [9].Motivated by these works and by the early theoretical framework established in [5, 6, 7, 8], we investigatePDEs with L subdifferential terms. The PDE is either an elliptic PDE coming from a variational principleor a parabolic PDE coming from a gradient flow of a convex functional. In either case, the L term in theconvex functional leads to a subgradient term in the PDE. Fortunately, the subgradient term has a simpleexplicit form, so that the PDEs are amenable to analysis and computation.The goal of this work is to present fast computational schemes for these modified PDEs, provide someadditional theoretical insights, and show some connections to known physical equations. Our starting pointis the convex functional:(1.1) E ( u ) = Z
12 ( ∇ u ) · M ( ∇ u ) − uf + γ | u | dx, where γ ≥ M = M ( x ) is a symmetric, positive definite matrix as a function of x , and f = f ( x ) or f = f ( x, t ) will be a specified function depending on x or ( x, t ). Define the partial differential operator Au = −∇ · ( M ∇ u ). Minimization of E ( u ) for f = f ( x ) leads to the following elliptic PDE(1.2) Au = f − γp ( u ) , and gradient descent ∂ t u = − ∂ u E ( u ), starting from initial data g ( x ), leads to the following parabolic PDE(1.3) u t + Au = f − γp ( u ) u ( x,
0) = g ( x ) , n which p ( u ) is a subgradient of k u k L , i.e. , k v k ≥ k u k + h v − u, p ( u ) i , for any u and v , where h , i denotesthe L inner product.The paper is divided as follows: in Section 2, we provide the general formulation of the problem. InSection 3, we review known results and present various properties of solutions to the modified PDEs. Thenumerical implementation and simulations are presented in Sections 4 and 5, and we conclude in Section 6.2. Problem Formulation
The problem we consider in this work is to numerically solve the following PDE(2.1) u t + Au = f − γp ( u ) u ( x,
0) = g ( x ) , and to verify theoretical results. The difficulty with such equations is the multivalued nature of the subgra-dient term. Fortunately for this type of equation, we can explicitly identified the subgradient as(2.2) p ( u ) = sign( u ) if | u | > | q |≤ | f − γq | if u = 0.Note that if u = 0 and | f ( x ) | ≤ γ , then p = f ( x ) /γ . This specification for u was proved in general in [4, 5].It can be shown directly from Equations (1.2) and (1.3), as follows. For u = 0 in an open set, the left sideof the equations is 0 so that f ( x ) − γp ( u ) = 0, which is only possible if f ( x ) ≤ γ and p ( u ) = f ( x ) /γ . Thevalue of p ( u ) on a lower dimensional set does not matter, since the value of the forcing terms on a lowerdimensional set does not affect the solution u of the differential equations. For the elliptic equation (1.2)one can also show directly that this identification of p ( u ) gives u = 0 as the unique minimizer of E ( u ) (seeAppendix). 3. Various Properties
In this section we recall the established existence theory for the elliptic equation (1.2) and the parabolicequation (1.3), and provide some further insights to the behavior of solutions.3.1.
Review of Theoretical Results.
Equation (1.2) is related to the general class of elliptic equation: − ∆ u = F ( u ) , where F contains a discontinuous component. The existence and uniqueness of the solution u are studied in[20, 19, 15]. Solutions also satisfy the standard maximum and comparison principles given the correct signof F . The solutions are compactly supported in both the elliptic and parabolic case, under some additionalconditions [7, 8]. For the parabolic equation, the solutions are Lipschitz continuous and right differentiable intime. Furthermore, solutions exhibit finite speed of propagation [8]. More precisely, let S ( t ) be the supportset of u ( x, t ), then for small times t : • if u ( x,
0) does not vanish on ∂S (0) , then S ( t ) ⊂ S (0) + B ( c p t log( t )) , • if u ( x,
0) and ∇ u ( x,
0) vanishes on ∂S (0), then S ( t ) ⊂ S (0) + B ( c √ t ) , where B ( r ) is the ball of radius r centered at the origin. In a simple case, we can construct the exact boundsin order to verify the convergence of the method to a known solution.At a number of places in the manuscript, we will simplify the presentation by assuming that x ∈ R andthat M = 1, so that the elliptic PDE (1.2) becomes Laplace’s equation with nonlinear forcing:(3.1) u xx = − f + γp ( u ) , and the parabolic PDE (1.3) becomes the heat equation with nonlinear forcing:(3.2) u t − u xx = f − γp ( u ) , .2. A Free Boundary Formula.
In 1D, consider the following equation(3.3) u t − u xx = ( f ( x ) − γ, | x | < a ( t )0 , | x | > a ( t ) u ( x,
0) = 0 . For simplicity assume that f ( x ) = f ( | x | ) and f is a decreasing function with f ( | x | ) → | x | → ∞ . Denote a ≥ f ( a ) = γ and assume that f x ( a ) = 0. Then, the free boundary’s endpoint is governed by(for small time t ):(3.4) a ( t ) = a + a √ t, for some a ≥ Support Size.
Since it is known that the support is compact, we would like to estimate its size. Infact, by integrating Equation (1.2) (see Appendix), the support of u satisfies(3.5) | supp( u ) | ≤ γ − Z supp( u ) | f | dx. A slight modification of (3.5) shows that for any nonnegative α and β with α + β = 1, we have(3.6) | supp( u ) | ≤ ( αγ ) − Z ( | f | − βγ ) + dx. In this inequality, the superscript + denotes the positive part; i.e. , ( x ) + = max( x, | supp ( x,t ) u ( x, t ) | ≤ ( αγ ) − (cid:18)Z | g | dx + Z Z ( | f | − βγ ) + dx dt (cid:19) , (3.7)for any nonnegative α and β with α + β = 1.3.4. L Contraction and Total Variation Diminishing.
Let u and v be solutions of Equation (3.2)with initial data g ( x ) and h ( x ), respectively. First, note that for any subgradient p of a convex functional,we have(3.8) sign( u − v )( p ( u ) − p ( v )) ≥ . We wish to show that the solutions are L contractive and TVD by computing the following: ddt || u − v || L = ddt Z | u − v | > | u − v | dx = Z | u − v | > sign( u − v )( u t − v t ) dx = Z | u − v | > sign( u − v )( u − v ) xx − γ sign( u − v )( p ( u ) − p ( v )) dx. The first term is zero by the divergence theorem and the second term is negative by Equation (3.8), so wehave ddt || u − v || L ≤
0, and thus the modified PDE is an L contraction. Moreover, if we take h ( x ) = g ( x + δ )for any δ > ddt k u ( x, t ) − u ( x + δ, t ) k L ≤ . Dividing the equation above by δ and taking the supremum over all δ , the following inequality holds: ddt || u || T V ≤ . Therefore, Equation (3.2) is TVD. .5. Entropy Condition.
The L contraction and TVD results are directly analogous to those that areobtained by solving the viscosity regularized nonlinear conservation laws: w ǫt = ǫ w ǫxx − f ( w ǫ ) x , for ǫ >
0. Then by letting ǫ →
0, one recovers the unique inviscid limit, see [21].We can also easily obtain an “entropy inequality” in the same spirit. Consider the scaled modified heatequation:(3.9) u t = ǫ u xx − γp ( u ) . We deliberately put an ǫ in front of the diffusion term to emphasize the similarities to the theory of scalarconservation laws. The following argument holds in more general cases.Let K ( u ) be a convex function of u with subgradient q ( u ). Multiplying Equation (3.9) by the subgradient(as in [21]) yields:(3.10) ddt K ( u ) ≤ ǫ d dx K ( u ) − γq ( u ) p ( u ) . For example, if K ( u ) = | u | , then whenever u = 0, we have(3.11) | u | t ≤ ǫ | u | xx − γ. We integrate Equation (3.10) over the region S ( t ), the support set of u ( x, t ) defined in Section 3.3, to get(3.12) ddt Z S ( t ) K ( u ) dx ≤ − γ Z S ( t ) q ( u ) p ( u ) dx, since the spatial gradient is zero along the boundary. By choosing K ( u ) = a | u | a for a ≥
1, Equation (3.12)provides L a estimates of the solutions. Furthermore, if K ( u ) = ( u − c ) + for c >
0, then(3.13) ddt Z S + c ( t ) ( u − c ) + dx ≤ − γ |S + c ( t ) | , where S + c ( t ) is the set of x for which u ( x ) > c .3.6. Regularity.
We can show that the solutions of the Laplace’s equation (3.1) and of the heat equation(3.2) are smooth. Let Ω + , Ω − , and Ω denote the sets { u > } , { u < } and { u = 0 } , respectively. Thenthe solution u of the Laplace’s equation (3.1) can be represented by(3.14) u ( x ) = Z Ω + G ( x − y )( f ( y ) − γ ) dy + Z Ω − G ( x − y )( f ( y ) + γ ) dy, and the solution of the heat equation (3.2) can be written as(3.15) u ( x, t ) = Z G ( x − y, t ) g ( y ) dy + Z t Z Ω + ( s ) G ( x − y, t − s )( f ( y ) − γ ) dyds + Z t Z Ω − ( s ) G ( x − y, t − s )( f ( y ) + γ ) dyds, in which the Green’s function G ( x, t ) for the heat equation and the Green’s function G ( x ) for the Laplace’sequation are given by(3.16) G ( x ) = | x | / ,G ( x, t ) = (4 πt ) − / exp( − x / t ) . From these formulas, if f is continuous, then one can see that u is C ( x ) and C ( t ) away from u = 0 andthat u is C ( x ) everywhere. .7. Traveling Wave.
To demonstrate finite speed of propagation, consider the 1D-traveling wave solution u ( x, t ) = v ( s ) for s = x − σt , of the Equation (3.2) with no forcing term. To be specific, we will assume that v ( s ) ≥ s ≥ v ( s ) = 0 for s ≤
0. We see that v must satisfy the ODE(3.17) v ss + σv s − γ = 0 , subject to the conditions(3.18) v (0) = v ′ (0) = 0 . The general solution of Equation (3.17) is v ( s ) = ( γσ s + c e − σs + c , s ≥ , otherwise . (3.19)The boundary conditions imply c = − c = γσ , so that the traveling wave solution of Equation (3.2) is u ( x, t ) = ( γσ ( x − σt ) + γσ (cid:0) e − σ ( x − σt ) − (cid:1) , x ≥ σt , otherwise . We see that in this case we have one sided support.
Remark . This traveling wave solution is used as a reference solution to compute the error for our numericalscheme (see Section 5.1). Also, the simple analytic form shows that solutions with non-trivial support setsare easy to find in the modified PDE.3.8.
An Exact Solution.
We construct the exact solution of Equation (3.1) with nonnegative force f =(1 + x ) − / and γ ∈ [0 , u = − (1 + x ) / + 12 γx + c, | x | ≤ a , | x | > a. where, c = γ + γ − , a = p γ − − . The boundary value a and constant c are determined so that u ( ± a ) = u x ( ± a ) = 0. At the boundary of thesupport, f ( ± a ) = γ < γ . The results show that the solution is nonnegative for nonnegative f , and thathaving | f ( x ) | ≤ γ does not imply p ( u ( x )) = f ( x ) γ .4. Numerical Implementation
Given an elliptic operator A , we would like to solve problems of the form: Au + ∂ k u k L ∋ f (4.1)or u t + Au + ∂ k u k L ∋ f (4.2)which corresponds to the elliptic or parabolic equations, respectively. We will present two methods to doso. The first scheme is semi-implicit (also known as implicit-explicit or proximal gradient method), wherethe subgradient term is discretized forward in time and the diffusion term is lagged. We apply this methodto solve the time dependent equations. The second scheme is the Douglas-Rachford method, which we useto solve both the elliptic problem and the parabolic problem. Both methods can handle the multivaluednature of the subgradient ∂ k u k L . In this section, we denote h and τ the space and time steps of the finitedifference schemes. .1. Implicit-Explicit Scheme (Proximal Gradient Method).
From the numerical perspective, themultivalued term ∂ k u k L is the main source of difficulties, since the value is ambiguous. However, anoperator of the form I + σ ∂F ( where F is convex) has an easy-to-compute inverse. The inverse operator( I + σ ∂F ) − , also known as the resolvent or proximal operator, prox σF ( · ), can be found by solving thefollowing optimization: ( I + σ ∂F ) − ( z ) = argmin v || v − z || L + σF ( v ) . (4.3)For example, if F ( u ) = || u || L and thus ∂F ( u ) = ∂ k u k L , we have:( I + σ ∂ k · k L ) − ( z ) = argmin v || v − z || L + σ || v || L = S ( v, σ ) , where the shrink operator, S , is defined point-wise as S ( v, σ ) := max( | v | − σ, v | v | .Using the proximal operator, we will write the discretization of Equation (1.3) in a semi-implicit form.We first discretize Equation (1.3) in time: u n +1 − u n + τ Au + τ ∂ k u k L ∋ τ f. (4.4)Then to apply the proximal gradient method, the last two terms on the left are evaluated at n and n + 1 asfollows: u n +1 − u n + τ Au n + τ ∂ k u n +1 k L ∋ τ f. (4.5)The resulting iterative scheme is: u n +1 = S ( u n − τ Au n + τ f, τ ) . (4.6)For example, for the heat equation, where A = − ∆, the iterative scheme is: u n +1 = S ( u n + τ ∆ u n + τ f, τ ) , (4.7)and is convergent given τ ≤ h . This scheme has the same complexity as the corresponding standard explicitmethods for PDE.4.2. Alternating Direction Implicit (Douglas-Rachford) Method.
The Douglas-Rachford algorithmfor nonlinear multivalued evolution equation was studied in [25]. Denote Bu := ∂ k u k L , the iterative schemefor Equation (1.3) is(4.8) u n +1 = ( I + τ B ) − (cid:2) ( I + τ A ) − ( I − τ B ) + τ B (cid:3) u n , which can be rewritten as:(4.9) u n +1 = ( I + τ B ) − ˜ u n ˜ u n +1 = ˜ u n + ( I + τ A ) − (2 u n +1 − ˜ u n ) − u n +1 . It was shown that the method is unconditionally stable and convergent for all τ > u n converges to a solution of the stationary equation (1.2). For the sandpile problem [24],the operators A and B are chosen specifically as follows:(4.10) Au = − ∆ u − f, Bu = ∂ k u k L , so that the operation for u n +1 in the iterative process, Equation (4.9), is a shrink. The correspondingproximal operators are prox τF ( z ) = ( I + τ A ) − ( z ) = ( I − τ ∆) − ( z + τ f )prox τG ( z ) = ( I + τ B ) − ( z ) = S ( z, τ ) , where F ( z ) = k∇ z k L − h f, z i and G ( z ) = k z k L . To compute ( I − τ ∆) − numerically, we use the FFT,where the discrete Laplacian ∆ h u is viewed as the convolution of u with the finite difference stencil. Remark . Since the shrink operator is the last step of the iterative process, this method provides anumerically well-defined support set for u , making it easier to locate the free boundary. a) t = 0 (b) t = 0 . (c) t = 1 . (d) t = 2 . Figure 1.
Numerical solution starting with an initial traveling wave profile with σ = 2 and γ = 0 .
05 computed using 500 grid points.5.
Computational Simulations
In this section we show convergence of our numerical scheme to known solutions, approximations to thesupport set evolution, and numerical solutions for higher dimension.5.1.
Numerical Convergence.
In Figure 1, we solve Equation (3.2) (with γ = 0 .
05) using the implicit-explicit scheme (Equation (4.7)). The initial data is taken to be the traveling wave profile (Equation (3.19))with speed σ = 2. The numerical solution has the correct support set and speed of propagation, validatingthe traveling wave solution as well as the numerical method.This is further confirmed in Figure 2, where the numerical solution is compared to the exact solution. Tocompute the error, we use the following norms:Error q ( h ) = max n || u nh − u exact || q , where q = 1 , , ∞ and u nh is the solution at t n with space resolution h . The errors in these three norms areplotted along side the line representing the second order (dashed line) convergence.To test the stability of these traveling wave solutions, we initialize our numerical scheme with the travelingwave profile perturbed by uniformly random noise sampled from [0 , . One Dimensional Heat Equation.
In Figure 4, the plot shows the modified heat equation (Equa-tion (3.2)) with zero initial data and force f ( x ) = 2 e − x . The solutions evolves upward in time with theirsupport sets marked by red circles. We see that the computed solutions are indeed compactly supported inspace, as the theory states. The corresponding table provides a least squares fit to estimate the coefficient a from Equation (3.4) under grid refinement. We see that the coefficient a approaches the value 1 quicklywithin some small approximation error, which is used to verify that our numerical approximation is valid. Figure 2.
Convergence analysis using the L (dotted line), L (dashed line), and L ∞ (solidline) norms in space and L ∞ norm in time. The x -axis is the log of the grid resolution h andthe y -axis is the log of the Error. The blue dashed line represents second order convergence.Grid Size 128 256 512 1024 2048 4096 8192 L -Error 0.4601 0.2319 0.1133 0.0569 0.0284 0.0143 0.0072 Table 1.
Error between our numerical solution and the analytic solution of the Signum-Gordon Equation.5.3.
Two Dimensional Heat Equation.
In Figure 5, we compute the solution of Equation (3.2) with γ = 2 and f = 0. In this case, we apply the parabolic Douglas-Rachford algorithm, which allows for largertime-steps. The initial data is a smoothed indicator function on the star shaped domain. In Figure 6, thecorresponding support set of Figure 5 is shown. The support set grows outward to a maximum size andretracts inward as the solution decays to zero. The solution is identically zero at time t = 0 . Graph Diffusion.
In higher dimensions, we can consider the standard normalized diffusion equation:(5.1) u t = L g u := − (cid:16) I − D − / AD − / (cid:17) uu ( x,
0) = g ( x ) , where L g is the graph Laplacian, A is the adjacency matrix, and D is the degree matrix. For more on thegraph Laplacian, see [13, 34].In Figure 7, the points represent the projection of vectors from R and each point is connected to manyothers in a non-local fashion. For the initial data, we concentrate the mass on one point in the far left,specifically, the u ( x j ,
0) = δ j, where δ j,k is the Kronecker delta function. As the system evolves governedby Equation (5.1), the solution becomes strictly positive quickly.The modified equation is:(5.2) u t = − (cid:16) I − D − / AD − / (cid:17) u − γp ( u ) ,u ( x,
0) = g ( x ) . In Figure 8, we begin with the same initial condition and see that over time the support set does notgrow past a bounded region if u evolves as in (5.2). Therefore, numerically we show that the support is offinite size for the case of graph diffusion. In Figure 8(d), the solution begins to decay to zero which causesits support set to retract towards the initial support before vanishing. a) t = 0 (b) t = 0 . (c) t = 0 . (d) t = 0 . (e) t = 1 . Figure 3.
Numerical solution starting with an initial traveling wave profile perturbed byuniformly random noise sampled from [0 , .
05] with σ = 2 and γ = 0 .
05. This is solved ona grid of 500 points.5.5.
Signum-Gordon Equation.
The signum-Gordon equation has an interpretation as an approximationto certain physical models [2, 1, 3]. The equation takes the form of a second order nonlinear hyperbolicequation:(5.3) u tt − ∆ u = − sign( u ) u ( x,
0) = g ( x ) u t ( x,
0) = g ( x ) , and exhibits both compactly supported traveling waves and oscillatory (stationary) soliton-like structures.This equation can be derived from the Lagrangian with the following L potential: L = Kinetic − Potential = 12 | u t | − |∇ u | − | u | . umber of grid points Estimate of a
256 0.948512 0.9791024 0.9852048 0.9914096 0.9958192 0.99716384 0.997
Figure 4.
The graph is a 1D simulation of the heat equation with the subgradient term,zero initial data, and a Gaussian forcing function centered around zero, f ( x ) = 2 e − x . Thesolutions are growing upward in time and their support sets are marked by red circles. Thetable shows the estimate of the coefficient a from Equation (3.4) under grid refinement.The equation of motion can be derived from the Lagrangian: u tt − ∆ u = − p ( u ) u ( x,
0) = g ( x ) u t ( x,
0) = g ( x ) , which is the same as Equation (5.3) by replacing the sign( u ) term with the subgradient p ( u ).To discretize the problem, we apply the ideas from the proximal gradient method, by placing p ( u ) in thefuture: u n +1 − u n + u n − − τ ∆ u n = − τ p ( u n +1 ) , and thus, u n +1 = S (2 u n − u n − + τ ∆ u n , τ ) . In Figure 9, we plot our numerical approximation to the traveling wave solution found in [1]. Since thetraveling wave profile is also known analytically, we show numerical convergence of our scheme as h → + (see Table 1). Also, in Figure 10, we show the time evolution of an oscillatory compact soliton-like structurewhich appears in [2, 3]. These examples show the range of behaviors that appear via the addition of an L subgradient term.5.6. Divisible Sandpile.
As a model for self-assembly and internal diffusion limited aggregation, the sand-pile problem has received attention recently [29, 24, 18, 22, 23]. The problem is posed discretely, but hasthe following continuous formulation for the divisible sandpile problem [24, 18]:∆ u = 1 − f, if u ≥ , (5.4)where f is some non-negative external force. By multiplying Equation (5.4) with u and integrating over R ,the associated variational energy is: min u Z u ≥ |∇ u | + u − uf dx . (5.5)There are several choices for relaxing the constraint u ≥
0, in particular, we use the following:min u Z |∇ u | + | u | − uf dx . (5.6)It can be shown (via maximum principle) that for f ≥ L sandpile problem is:∆ u = p ( u ) − f, (5.7) a) t = 0 (b) t = 3 . × − (c) t = 3 . × − (d) t = 1 . × − (e) t = 4 . × − (f) t = 0 . Figure 5.
Solutions of the initial value problem (no forcing term) computed on a 500 by500 grid with γ = 2 at times indicated. The solution smoothes out and decays to zero.and is solved numerically via the Douglas-Rachford algorithm (see Equation (4.9)). Note that if the externalforce is a finite sum of characteristic functions f = P Nj =1 α j χ S j where S j are compact sets and α j ≥
0, thenby integrating Equation (5.7) over R we get:(5.8) | supp( u ) | = N X j =0 α j | S j | , since u ≥ u ) is compact. This refers to preservation of mass. a) t = 0 (b) t = 3 . × − (c) t = 3 . × − (d) t = 1 . × − (e) t = 4 . × − (f) t = 0 . Figure 6.
Support set of the initial value problem in Figure 5. The support set growsoutward to a maximum size and retracts inward as the solution decays to zero.In Figure 11, we take f = χ S + χ S , where S and S are the two overlapping square domains (on theleft). The support set of u , given in Figure 11 (right), agrees with direct numerical simulation of the discretesandpile problem. The direct simulation follows a topping rule described in [24].In Figures 12-14, we take f = αχ S where S is the shape given in Figures 12-14 (the top left), and α = 2 . , . .
5, respectively. The support set of u is given in Figures 12-14 (the bottom right) withintermediate calculation shown in Figures 12-14 (the remaining plots). To verify that the solutions fromour algorithm correspond to the correct solutions for the sandpile problem, we use the mass conservationproperty, Equation (5.8). Unlike direct simulation, our method also calculates the function u as shown in a) t = 0 (b) t = 47 . (c) t = 475 (d) t = 1425 Figure 7.
Solution of the initial value problem diffusing standard normalized graph Laplacian. (a) t = 0 (b) t = 2 . (c) t = 6 . (d) t = 11 . Figure 8.
Solution of the initial value problem with the subgradient term, γ = 5 × − . igure 9. Our numerical approximation to a compact traveling wave solution to theSignum-Gordon equation.Figure 15. One of the benefits of our approach is that the solutions can be computed quickly, for example,our method is at least 8 times faster than direct simulation (76 seconds vs. 652 seconds) at approximatingthe solution found in Figure 15. 6.
Conclusion
By adding the subdifferential of L to certain PDEs, we have shown (numerically and theoretically)various properties of the solutions. These problems arise from physical models as well as exact relaxationof other PDEs, and could provide useful tools in computing fast approximations to nonlinear problems witha compactly supported free boundary. This is all in the spirit of borrowing the key idea from compressedsensing, that L regularization implies sparsity of discrete systems [17], and transferring it to classicalproblems in PDE. See [32, 27] for earlier work in this direction. Acknowledgments
The authors would like to thank Farzin Barekat, Jerome Darbon, William Feldman, Inwon C. Kim andJames H. von Brecht for their helpful discussions and comments.R. Caflisch was supported by ONR N00014-14-1-0444. S. Osher was supported by ONR N00014-14-1-0444and N000141110719. H. Schaeffer was supported by NSF 1303892 and University of California PresidentsPostdoctoral Fellowship Program. G. Tran was supported by ONR N00014-14-1-0444 and N000141110719.
References [1] H. Arod´z, P. Klimas, and T. Tyranowski. Scaling, self-similar solutions and shock waves for V-shaped field potentials.
Physical Review E , 73(4):046609, 2006.[2] H. Arod´z, P. Klimas, and T. Tyranowski. Compact oscillons in the signum-Gordon model.
Physical Review D , 77(4):047701,2008.[3] H. Arod´z and Z. ´Swierczy´nski. Swaying oscillons in the signum-Gordon model.
Physical Review D , 84(6):067701, 2011.[4] J. P. Aubin and A. Cellina.
Differential inclusions: set-valued maps and viability theory . Springer-Verlag New York, Inc.,1984.[5] H. Br´ezis.
Operateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert , volume 5. Elsevier,1973.[6] H. Br´ezis.
Monotone operators non linear semi-groups and applications . Universit´e Pierre et Marie Curie, Laboratoired’Analyse Num´erique, 1974.[7] H. Br´ezis. Solutions with compact support of variational inequalities.
Russian Mathematical Surveys , 29(2):103–108, 1974.[8] H. Br´ezis and A. Friedman. Estimates on the support of solutions of parabolic variational inequalities.
Illinois Journal ofMathematics , 20(1):82–97, 1976.[9] S. L. Brunton, J. H. Tu, I. Bright, and J. N. Kutz. Compressive sensing and low-rank libraries for classification of bifurcationregimes in nonlinear dynamical systems. arXiv preprint arXiv:1312.4221 , 2013. igure 10. The dynamics of an oscilatory compact solution of the Signum-Gordon equation. [10] E. J. Cand`es, Y. C. Eldar, T. Strohmer, and V. Voroninski. Phase retrieval via matrix completion.
SIAM Journal onImaging Sciences , 6(1):199–225, 2013.[11] E. J. Cand`es, X. Li, Y. Ma, and J. Wright. Robust principal component analysis?
Journal of the ACM (JACM) , 58(3):11,2011.[12] E. J. Cand`es, J. Romberg, and T. Tao. Robust uncertainty principles: Exact signal reconstruction from highly incompletefrequency information.
Information Theory, IEEE Transactions on , 52(2):489–509, 2006.[13] F. R. Chung.
Spectral graph theory , volume 92. American Mathematical Soc., 1997.[14] P. L. Combettes and J.-C. Pesquet. Proximal splitting methods in signal processing. In
Fixed-point algorithms for inverseproblems in science and engineering , pages 185–212. Springer, 2011.[15] M. G. Crandall and A. Pazy. Semi-groups of nonlinear contractions and dissipative sets.
Journal of functional analysis ,3(3):376–418, 1969.[16] A. d’Aspremont, L. El Ghaoui, M. I. Jordan, and G. R. Lanckriet. A direct formulation for sparse PCA using semidefiniteprogramming.
SIAM review , 49(3):434–448, 2007.[17] D. L. Donoho. Compressed sensing.
Information Theory, IEEE Transactions on , 52(4):1289–1306, 2006.[18] A. Fey, L. Levine, and Y. Peres. Growth rates and explosions in sandpiles.
Journal of Statistical Physics , 138(1-3):143–159,2010.[19] T. Kato. Accretive operators and nonlinear evolution equations in Banach spaces. In
Proc. Symp. in Pure Math , volume 18,pages 138–161, 1970.[20] Y. Komura. Differentiability of nonlinear semigroups.
Journal of the Mathematical Society of Japan , 21(3):375–402, 1969. igure 11. A two-region sandpile problem, where each of the larger squares defines theset S j , for j = 1 ,
2. The darker blue region has no mass, the lighter blue and yellow regionhas a mass density of 1, and the red (overlap) region has a density of 2. On the left, theregion of positive mass is displayed.
Figure 12.
The iterative evolution of our sandpile problem algorithm applied to a flower-shaped region S on the top left with f = 2 χ S . The final state appears in the bottom rightcorner. igure 13. The iterative evolution of our sandpile problem algorithm applied to the fractalregion S on the top left with f = 1 . χ S . The final state appears in the bottom right corner. [21] P. D. Lax. Hyperbolic systems of conservation laws and the mathematical theory of shock waves , volume 11. SIAM, 1973.[22] L. Levine, W. Pegden, and C. K. Smart. Apollonian structure in the abelian sandpile. arXiv preprint arXiv:1208.4839 ,2012.[23] L. Levine, W. Pegden, and C. K. Smart. The apollonian structure of integer superharmonic matrices. arXiv preprintarXiv:1309.3267 , 2013.[24] L. Levine and Y. Peres. Strong spherical asymptotics for rotor-router aggregation and the divisible sandpile.
PotentialAnalysis , 30(1):1–27, 2009.[25] P.-L. Lions and B. Mercier. Splitting algorithms for the sum of two nonlinear operators.
SIAM Journal on NumericalAnalysis , 16(6):964–979, 1979.[26] L. J. Nelson, G. L. Hart, F. Zhou, and V. Ozoli¸nˇs. Compressive sensing as a paradigm for building physics models.
PhysicalReview B , 87(3):035125, 2013.[27] V. Ozoli¸nˇs, R. Lai, R. Caflisch, and S. Osher. Compressed modes for variational problems in mathematics and physics.
Proceedings of the National Academy of Sciences , 110(46):18368–18373, 2013.[28] V. Ozoli¸nˇs, R. Lai, R. Caflisch, and S. Osher. Compressed plane waves yield a compactly supported multiresolution basisfor the Laplace operator.
Proceedings of the National Academy of Sciences , 111(5):1691–1696, 2014.[29] W. Pegden and C. K. Smart. Convergence of the abelian sandpile.
Duke Mathematical Journal , 162(4):627–642, 2013.[30] X. Qi, R. Luo, and H. Zhao. Sparse principal component analysis by choice of norm.
Journal of multivariate analysis ,114:127–160, 2013.[31] B. Recht, M. Fazel, and P. A. Parrilo. Guaranteed minimum-rank solutions of linear matrix equations via nuclear normminimization.
SIAM review , 52(3):471–501, 2010.[32] H. Schaeffer, R. Caflisch, C. D. Hauck, and S. Osher. Sparse dynamics for partial differential equations.
Proceedings of theNational Academy of Sciences , 110(17):6634–6639, 2013.[33] S. Setzer. Split Bregman algorithm, Douglas-Rachford splitting and frame shrinkage. In
Scale space and variational methodsin computer vision , pages 464–476. Springer, 2009.[34] U. Von Luxburg. A tutorial on spectral clustering.
Statistics and computing , 17(4):395–416, 2007. igure 14. The iterative evolution of our sandpile problem algorithm applied to the fractalregion S on the top left with f = 1 . χ S . The final state appears in the bottom right corner. Figure 15.
The solution u from Figure 14 (bottom right). Appendix A. appendix A.1.
Proof of a free boundary formula in Section 3.2.
We derive the short time asymptotic equationfor the support set Equation (3.2). First, we provide a natural boundary condition for the problem.
Flux Condition.
Let u ( x, t ) ∈ C ( C ( R ); (0 , T )) and u t ∈ L ∞ ( C ( R ); (0 , T )) be a solution to (A.1) u t − u xx = h ( x, t, γ ) . ssume that there exists a positive valued function a ∈ C (0 , T ) such that h = 0 for | x | > a ( t ) , g = 0 for | x | > a (0) , and the exterior mass, m ( t ) = ∞ Z a ( t ) u ( x, t ) dx, is conserved, then u ( a ( t ) , t ) = 0 and u x ( a ( t ) , t ) = 0 . To derive this condition, consider the heat equation (A.1). Differentiate the one sided mass in time yields: dmdt = − u ( a ( t ) , t ) a ′ ( t ) + ∞ Z a ( t ) u t ( x, t ) dx = − u ( a ( t ) , t ) a ′ ( t ) + ∞ Z a ( t ) u xx ( x, t ) dx = − u ( a ( t ) , t ) a ′ ( t ) − u x ( a ( t ) , t )= − F ( t ) , in which F is the flux across the moving boundary x = a ( t ).We now can see that if the flux across a moving boundary x = a ( t ) is zero ( i.e. the mass is conserved),we have(A.2) F ( t ) = u ( a ( t ) , t ) a ′ ( t ) + u x ( a ( t ) , t ) = 0 . This is the natural boundary condition for this problem. In the time-dependent region F = { ( x, t ) : x > a ( t ) } ,the initial data g , force h and and incoming flux F are all zero, so that the solution is identically zero. Inparticular, u = u x = 0 on x = ± a ( t ) . Next, consider the following equation: u t − u xx = ( f ( x ) − γ, | x | < a ( t )0 , | x | > a ( t ) u ( x,
0) = 0 . For simplicity assume that f ( x ) = f ( | x | ) and f is a decreasing function with f ( | x | ) →
0. Denote a ≥ f ( a ) = γ and w.l.o.g. f x ( a ) = 0. By studying the exterior mass of Equation (A.1), we want toshow that in small time: a ( t ) = a + a √ t, for some a ≥ a ( t ) such that the exterior mass of Equation (A.1) is zero: m ( t ) = ∞ Z a ( t ) dx t Z ds a ( s ) Z − a ( s ) G ( x − y, t − s )( f ( y ) − γ ) dy. where we use the Greens formula to represent u . Since a ( t ) is an increasing function, we have y ≤ a ( s ) ≤ a ( t ) ≤ x. Therefore, for t small, the Green’s function G ( x − y, t − s ) is sharply peaked near the point y = a ( t ) , s = t, x = a ( t ) . So we can replace ( f ( y ) − γ ) by the first few terms in its Taylor expansion f ( y ) − γ = ( y − a ) f + O (( y − a ) ) , n which f = f x ( a ). Also, since G ( x − y, t − s ) decays exponentially as y → −∞ , we replace the lower limit y = − a ( s ) by −∞ . Now the mass can be approximated by m ( t ) = f ∞ Z a ( t ) dx t Z ds a ( s ) Z −∞ ( y − a ) G ( x − y, t − s ) dy Next we show the existence of a satisfying the following approximations a ( t ) = a + a √ t, and m ( t ) = 0 . We change the variables to x = x √ t + a , x ∈ [ a , ∞ ) ,y = y √ t + a , y ∈ ( −∞ , a √ s ] ,s = s t, s ∈ [0 , , and x = x a , x ∈ [1 , ∞ ) ,y = y a , y ∈ ( −∞ , √ s ] , and note that G ( x − y, t − s ) = t − / G ( x − y , − s ) = t − / G ( a ( x − y ) , − s ). Then m ( t ) = f t ∞ Z a dx Z ds a √ s Z −∞ y G ( x − y , − s ) dy = a f t ∞ Z dx Z ds √ s Z −∞ y G ( a ( x − y ) , − s ) dy . Consider the rescaled masses f m ( a ) = m ( t ) / ( f t ) and f m ( a ) = m ( t ) / ( a f t ); i.e. , f m ( a ) = ∞ Z a dx Z ds a √ s Z −∞ y G ( x − y , − s ) dy , f m ( a ) = a ∞ Z dx Z ds √ s Z −∞ y G ( a ( x − y ) , − s ) dy . As a → f m ( a ) goes to f m (0) = ∞ Z dx Z ds Z −∞ y G ( x − y , − s ) dy with f m (0) <
0. This shows that m ( t ) < a = 0.On the other hand, for a ≫ a G ( a ( x − y ) , − s ) is approximately the Dirac delta function at x = y , s = 1. At this point, we have y >
0, therefore f m ( a ) >
0. This shows that m ( t ) > a . Thus there exists a positive value a so that m ( t ) = 0.A.2. Proof of support size estimate in Section 3.3.
Proof.
First, observe that if γ ≥ max | f | , then the unique solution of Equation (1.2). is u ≡ . Indeed, if u = 0, since fγ ∈ [ − , p ( u ) = fγ and Equation (1.2) is satisfied.Now, take S = supp( u ) and integrating both sides of Equation (1.2) gives us Z ∂ S M ∇ u · N ds = − Z S f dx + γ sign( u ) |S| . ince the left hand side is nonpositive, we have | supp( u ) | ≤ γ − Z supp( u ) | f | dx. For the parabolic case, define the time dependent support set S ( t ) := supp( u ( x, t )). Differentiating theintegral of u over S ( t ) and using the boundary conditions ( i.e. , u = 0 on ∂ S ( t )) yields: ddt Z S ( t ) u ( x, t ) dx = Z S ( t ) u t dx = Z S ( t ) ∇ · M ∇ u + f − γp ( u ) dx. Because of the divergence theorem and the fact that M is positive definite, we have ddt Z S ( t ) | u ( x, t ) | dx ≤ Z S ( t ) | f | dx − γ |S ( t ) | . Integrating the expression in time yields the following bound on the support size: | supp ( x,t ) u ( x, t ) | ≤ Z S ( t ) | g | dx + Z Z S ( t ) | f | dx dt. (cid:3)(cid:3)