On a Monotone Scheme for Nonconvex Nonsmooth Optimization with Applications to Fracture Mechanics

A general class of nonconvex optimization problems is considered, where the penalty is the composition of a linear operator with a nonsmooth nonconvex mapping, which is concave on the positive real line. The necessary optimality condition of a regularized version of the original problem is solved by means of a monotonically convergent scheme. Such problems arise in continuum mechanics, as for instance cohesive fractures, where singular behaviour is usually modelled by nonsmooth nonconvex energies. The proposed algorithm is successfully tested for fracture mechanics problems. Its performance is also compared to two alternative algorithms for nonsmooth nonconvex optimization arising in optimal control and mathematical imaging.


Introduction
In this paper, we investigate a class of nonconvex and nonsmooth optimization problems, where the penalty is the composition of a nonsmooth nonconvex mapping with a linear operator and the smooth part is a least square-type term. Similar optimization problems, in the case where the operator inside the penalty coincides with the identity matrix, have attracted increasingly attention due to their applications to sparsity of solutions, feature selection, and many other related fields as, e.g. compressed sensing, signal processing, and machine learning (see, e.g. [1,2]). The convex nonsmooth case of the 1 norm has gained large popularity and has been thoroughly studied. The convexity allows to formulate efficient and globally convergent algorithms to find a numerical solution. Here, we mention [3,4] where the basis pursuit and the Lasso problems were introduced to solve 1 minimization problems.
Recently, increased interest has arisen towards nonconvex and nonsmooth penalties, such as the τ quasi-norm, with τ larger or equal to zero and less than 1 (see, e.g. [5][6][7][8][9][10]), the smoothly clipped absolute deviation (SCAD) [11,12], and the minimax concave penalty (MCP) [12,13]. The nonconvexity has been shown to provide some advantages with respect to the convex models. For example, it allows to require less data in order to recover exactly the solution (see, e.g. [14][15][16]) and it tends to produce unbiased estimates for large coefficients [11,17,18]. Note that all the previously mentioned works deal with the particular case where the operator coincides with the identity.
Nonconvex optimization problems as we consider, where the operator inside the penalty is different form the identity, arise also in the modelling of cohesive fractures in continuum mechanics, where the concavity of the penalty is crucial to model the evolution of the fracture energy released within the growth of the crack opening. Here, the operator is of importance to model the jump of the displacement between the two lips of the fractures. We refer to [19][20][21][22] and Sect. 3.1 for more details.
The study of these problems for nonconvex penalties, including τ , with τ strictly positive and less than 1, the SCAD and the MCP functionals, and for linear operators not necessarily coinciding with the identity, is also motivated by applications different from those arising in fracture mechanics. For example, in imaging the τ quasi-norm, with τ strictly positive and less than 1, of the numerical gradient of the solution has been proposed as a nonconvex extension of the total variation (like TV) regularizer (see, e.g. [6,10]) in order to reconstruct piecewise smooth solutions. The SCAD and the MCP penalties have been used for high-dimensional regression and variable selection methods in high-throughput biomedical studies [23]. We mention also that the SCAD has been proposed as a nonconvex penalty in the network estimation to attenuate the bias problem [24].
The main difficulties in the analysis of these problems come from the interplay between the nonsmoothness, the nonconvexity, and the coupling between coordinates which is described by the operator inside the penalty. Since standard algorithms are not readily available, the resolution of these problems requires the development of new analytical and numerical techniques.
In the present paper, we propose a monotonically convergent algorithm to solve this kind of problems. This is an iterative procedure which solves the necessary optimality condition of a regularized version of the original problem. A remarkable property of our scheme is the strict monotonicity of the functional along the sequence of iterates. The convergence of the iteration procedure is proved under the same assumptions that guarantee the existence of solutions.
The performance of the scheme is successfully tested to simulate the evolution of cohesive fractures for several different test configurations. Then, we turn to an issue of high relevance, namely the comparison between two alternative algorithms, the GIST "General Iterative Shrinkage and Thresholding" algorithm for τ minimization, with τ strictly positive and less than 1, and the FISTA "Fast Iterative Shrinkage-Thresholding Algorithm" for 1 minimization. The comparison is carried out with respect to the infimal value reached by the iteration procedure and with respect to computing time. Our results show that the monotone algorithm is able to reach a smaller value of the objective functional that we consider when compared to the one of GIST. Note that differently from GIST, the monotone scheme solves a system of nonlinear equations at each iteration level. We remark that in [25], GIST was compared with the IRLS "iterative reweighted least squares" algorithm, which is another popular scheme for τ minimization, with τ strictly positive and less than 1. The results of [25] show that GIST and IRSL have nearly the same performance, with only one difference which is speed, where GIST appears to be the faster one. An analogous procedure to the one proposed in the present paper was developed in [20] to solve similar problems where the nonconvex penalty coincides with τ quasinorm, with τ strictly positive and less than or equal to 1. With respect to [20], in the present paper, we deal with more general concave penalties. Moreover, we carry out several numerical experiments for diverse situations in cohesive fracture mechanics, comparing the behaviours for different concave penalties such as the SCAD, the MCP, and the τ penalty, with τ strictly positive and less than 1. Finally, in the present paper, we compare the performance of the scheme with that of GIST.
Let us recall some further literature concerning nonconvex nonsmooth optimization of the type investigated in the present paper. In [12,26], a primal-dual active set-type algorithm has been developed; in the case, the operator inside the penalty coincides with the identity. For more references, in this case, we refer to [20]. Concerning τ minimization, with τ larger than or equal to zero and less than or equal to 1 when the operator is not the identity, other techniques have recently been investigated. Here, we mention iteratively reweighted convex majorization algorithms [10], alternating direction method of multiplier (ADMM) [9] and finally a Newton-type solution algorithm for a regularized version of the original problem [6]. Finally we recall the paper [21], where a novel algorithm for nonsmooth nonconvex optimization with linear constraints is proposed, consisting of a generalization of the well-known nonstationary-augmented Lagrangian method for convex optimization. The convergence to critical points is proved and several tests were made for free-discontinuity variational models, such as the Mumford-Shah functional. The nonsmoothness considered in [21] does not allow singular behaviour of the type that the τ term, with τ larger than or equal to zero and strictly less than 1 does.
The paper is structured as follows. In Sect. 2, Sect. 2.1, we state the precise assumptions, in Sect. 2.2, we prove existence for the problem in consideration, in Sect. 2.3, we propose the monotone scheme to solve a regularized version of the original problem and we prove its convergence, and finally in Sect. 2.4, we study the asymptotic behaviour as the concavity and regularization parameters go to zero. In Sect. 3, we present the precise form of our scheme. In Sect. 3.1, we discuss our numerical experience for cohesive evolution of fracture mechanics and in Sect. 3.2, we compare the performance of our scheme to that of GIST for three different test cases, the academical M-matrix example, an optimal control problem, and a microscopy imaging example.

Assumptions
We consider where (iv) there exists a neighbourhood of zero where the function t → φ (t) t is monotone; Above monotonically increasing or decreasing are admitted. Throughout the rest of the paper, we will use the notation Under assumption (H), the following two cases are analysed: (a) (i) φ(t) is a constant, when |t| ≥ t 0 for some t 0 > 0; (ii) A is coercive, i.e. rank(A) = n. (b) (i) for some γ > 0, it holds φ(at) = a γ φ(t) for all t ∈ R and a ∈ R + ; (ii) Ker(A) ∩ Ker(Λ) = {0}.
Three popular examples of nonconvex penalties which satisfy (H) and the assumptions on φ in (a) or (b) are the following: satisfying (a)(i).

Remark 2.1
The singularity at the origin of the three penalties leads to sparsity of the solution. In the SCAD and the MCP, the derivative vanishes for large values to ensure unbiasedness.
The SCAD (smoothly clipped absolute deviation) ( [11,18]) has raised interest in relation to variable selection consistency and asymptotic estimation efficiency (see [18]). It can be obtained upon integration of the following formula for τ > 2 The MCP (minimax concave penalty) [13] can be recovered from the following formula It minimizes the maximum concavity sup 0<t 1 subject to the constraints φ (t) = 0 for any |t| ≥ λτ (unbiasedness) and φ (0 ± ) = ±λ (feature selection). The condition τ > 1 ensures the wellposedness of the thresholding operator.

Existence
First, we prove coercivity of the functional J in (1) under assumptions (a) or (b).

Lemma 2.1 Let assumptions (H) and either (a) or (b) hold. Then, the functional J in
Proof Under assumption (a), the coercivity of J follows trivially. Suppose now that (b) holds. Then, the result follows by similar arguments to that used in [20], Theorem 1 (where φ is the τ quasi-norm). We proceed by contradiction and we suppose that |x k | 2 → +∞ and J (x k ) is bounded. For each k, let x k = t k z k be such that t k ≥ 0, x k ∈ R n and |z k | 2 = 1. By (b) (i), we have Φ(Λz k ) = 1 t γ k Φ(Λx k ) and then since t k → +∞ and J (x k ) is bounded, we have for k → +∞ and hence lim k→+∞ |Az k | 2 2 + Φ(Λz k ) = 0. By compactness, the sequence {z k } has an accumulation pointz such that |z| = 1 andz ∈ Ker(A)∩Ker(Λ), which contradicts (b) (ii).
In the following theorem, we state the existence of at least a minimizer to (1) under either (a) or (b). We omit the proof since it follows directly by the continuity and coercivity of the functional in (1).

Theorem 2.1 Let assumptions (H)
and either (a) or (b) hold. Then, there exists at least one minimizer to problem (1).

Remark 2.2
We remark that when assumption (a) (i) holds but A is not coercive, existence can still be proven in case Λ ∈ R n×n is invertible. Indeed, by the invertibility of Λ, one can defineȳ = Λ −1x , wherex is a minimizer ofJ (x) = 1 2 |(AΛ −1 )x − b| 2 2 + Φ(x) and prove thatȳ is a minimizer of (1). The existence of a minimizer for the functionalJ was proven in [26], Theorem 2.1.
However, in our analysis, we cover the two cases (a) and (b) since when (a) (ii) is replaced by the invertibility of Λ, we cannot prove the coercivity of J , which is a key element for the convergence of the algorithm that we analyse (see the following section).

A Monotone Convergent Algorithm
Following [7], in order to overcome the singularity of the function φ(t) near t = 0, we consider for ε > 0 the following regularized version of (1) where for t ≥ 0 and hence Ψ ε is C 1 and by assumption (H) (iii) is concave on [0, ∞[. Moreover, we have also that the map t → Ψ ε (t 2 ) ∈ C 1 (] − ∞, ∞[). Concerning the denominator in (7), let us point out that due to concavity of φ and (H)(i), the derivative φ is non-

The necessary optimality condition for (5) is given by
the second addend is short for the vector with l-component For convenience of exposition in the following, we write (8) in the more compact notation where the l-component of the second addend is given by

This can equivalently be expressed as
In order to solve (9), the following iterative procedure is considered: Existence of a unique solution to (10) follows from (a) (ii), respectively, (b) (ii) and the fact that is considered as derivative from the right. We have the following convergence result.

Theorem 2.2 Assume (H) and either (a) or (b).
For ε > 0, let {x k } be generated by (10). Then, J ε (x k ) is strictly monotonically decreasing, unless there exists some k such that x k = x k+1 , and x k satisfies the necessary optimality condition (9). Moreover, every cluster point of x k , of which there exists at least one, is a solution of (9).
Proof The proof strongly depends on the coercivity of the functional J and it follows arguments similar to those of [7, Theorem 4.1].
Multiplying (10) by x k+1 − x k , we get Note that for each i = 1, . . . , n, we have By assumption (H) (iii), the function t → Ψ ε (t) is concave on [0, ∞), and thus Using (12) and (13), we obtain the estimate Then, using (11), (14) and the definition of J ε , we get From (15) and the coercivity of J ε , it follows that {x k } ∞ k=1 and thus {y k } ∞ k=1 are bounded. Consequently, from (15) and (7), there exists a constant κ > 0 such that In the latter case, from (10), we infer that x k solves (9), from which we conclude the first part of the theorem. From (16), we conclude that Since {x k } ∞ k=1 is bounded, there exists a subsequence andx ∈ R n such that x k l →x. By (17), we get Then, by using the coercivity of A under assumption (a) and the fact that We can now pass to the limit with respect to k in (10), to obtain thatx is a solution to (9).
In the following proposition, we establish the convergence of (5) to (1) as ε goes to zero. (5). Then any cluster point of {x ε } ε>0 , of which there exists at least one, is a solution of (1).

Proposition 2.1 Assume (H) and either (a) or (b). Denote by {x ε } ε>0 a solution to
Proof From the coercivity of J ε , we have that {x ε } ε is bounded for ε small. Hence, there exists a subsequence andx ∈ R n such that x ε l →x as l → ∞. By By the concavity of the function φ, , and, by choosing s = ε and t = 0 and by (18), we get for ε small enough By the definition of Ψ ε , (18) and (19), we obtain that Ψ ε (t) converges uniformly to Since x ε l solves (5) for ε = ε l , by letting l → ∞ and using (20), we easily get thatx is a solution of (1).

Asymptotic Behaviour as 0 and 0 for the Power Law
We discuss the asymptotics as λ and τ go to zero in (1) for φ(t) = |t| τ , τ ∈]0, 1[, which we repeat for convenience min x∈R n where A, b, Λ are as in (1) First, we analyse the convergence as λ → 0 for any fixed τ > 0. We denote by P the orthogonal projection of R m onto Ker(A * ) and setb = ( Theorem 2.3 Assume that ker (A) ∩ ker(Λ) = {0} and let τ > 0 be fixed. For each λ > 0, let x λ be a minimizer of (21). Then, every cluster point of {x λ }, of which there exists at least one, is a solution to (22).
Proof Letx ∈ R n be an arbitrary solution to Ax −b = 0. By optimality of x λ , we have We conclude that lim λ→0 |Ax λ −b| 2 2 = 0, and |Λx λ | τ τ ≤ |Λx| τ τ for allx satisfying Ax =b. (24) In particular, the families {Ax λ } and {Λx λ } are bounded in λ. Since by assumption we have ker(A) ∩ ker(Λ) = {0}, it follows that {x λ } is bounded. Hence, there exists a convergent subsequence x λ with some limitx. From (24), it follows thatx is a solution to (22). Now, we prove the convergence as τ → 0 for any fixed λ > 0 of (21) to the related 0 -problem min x∈R n where for any x ∈ R n |x| 0 = n k=1 |x k | 0 = number of nonzero elements of x. The precise statement is given in the following theorem.

Theorem 2.4
Assume that rank(A) = n, and that Λ ∈ R n×n is regular, and let λ > 0 be fixed. Then, every cluster point (of which there exists at least one) of solutions {x τ } to (21) converges as τ 0 to a solution of (25).

Remark 2.3
The assumption on ker(A) is caused by the fact that the | · | 0 functional is not radially unbounded on R n . Since Theorem 2.4 provides an existence result for (25), such an assumption is natural. To the best of our knowledge, existence of minimizers of (25) has only been addressed under assumption (a) (ii) (or assumptions implying it). In the context of algorithm development further, 1 -or 2 -regularization is often added, we refer to [12,29,30].

Algorithm and Numerical Results
For convenience, we recall the algorithm in the following form.

Remark 3.1
Note that an ε-continuation strategy is performed, that is, the procedure is performed for an initial value ε 0 and then ε is decreased up to a certain value. More specifically, in all our experiments, ε is initialized with 10 −1 and decreased up to 10 −12 .

Remark 3.2
The stopping criterion is based on the l ∞ -norm of Eq. (9) and the tolerance is set to 10 −3 in all the following examples, except for the fracture problem where it is of the order of 10 −15 .
In the following subsection, we present our numerical results in cohesive fracture mechanics. Then, in Sect. 3.2, the performance of our algorithm is compared to two other schemes for nonconvex and nonsmooth optimization problems.

Application to Quasi-Static Evolution of Cohesive Fracture Models
In this section, we focus on the numerical realization of quasi-static evolutions of cohesive fractures. These kinds of problems require the minimization of an energy functional, which has two components: the elastic energy and the cohesive fracture energy. The underlying idea is that the fracture energy is released gradually with the growth of the crack opening. The cohesive energy, denoted by θ , is assumed to be a monotonic non-decreasing function of the jump amplitude of the displacement, denoted by u . Cohesive energies were introduced independently by Dugdale [31] and Barenblatt [32]; we refer to [19] for more details on the models. Among the vast existing literature on fracture mechanics, we also point out [33][34][35], as some of the most significant references to us. Let us just remark that the two models differ mainly in the evolution of the derivative θ ( u ), that is, the bridging force, across a crack amplitude u . In Dugdale's model, this force keeps a constant value up to a critical value of the crack opening and then drops to zero. In Barenblatt's model, the dependence of the force on u is continuous and decreasing.
In this section, we test the τ -term 0 < τ < 1 as a model for the cohesive energy. In particular, the cohesive energy is not differentiable in zero and the bridging force goes to infinity when the jump amplitude goes to zero. Note also that the bridging force goes to zero when the jump amplitude goes to infinity.
We denote by u : Ω → R the displacement function. The deformation of the domain is given by an external force which we express in terms of an external displacement function g : Ω ×[0, T ] → R. We require that the displacement u coincides with the external deformation, that is, u| ∂Ω = g| ∂Ω . We denote by Γ the point of the (potential) crack, and by θ( u ) Γ the value of the cohesive energy θ on the crack amplitude of the displacement u on Γ . Since we are in a quasi-static setting, we introduce the time discretization 0 = t 0 < t 1 < · · · < t T = T and look for the equilibrium configurations which are minimizers of the energy of the system. This means that for each i ∈ {0, . . . , T } we need to minimize the energy of the system with respect to a given boundary datum g: u * ∈ argmin u=g(t i ) on ∂Ω J (u). The function a(·) measures the degree of homogeneity of the material, e.g. a(x) ≡ 1 means that the material is homogeneous.
In Sects. 3.1.1 and 3.1.2, we show our results for one-dimensional and twodimensional experiments, respectively.

One-Dimensional Experiments
We consider the one-dimensional domain Ω = [0, 1] and we chose the point of crack as the midpoint Γ = 0.5. We divide Ω into 2N intervals and approximate the displacement function with a function u h that is piecewise linear on Ω\Γ and has two degrees of freedom on Γ to represent correctly the two lips of the fracture, denoting with u − N the one on [0, 0.5] and u + N the one on [0.5, 1]. We discretize the problem in the following way where, if i ≤ N , we identify u N = u − N , while for i > N , u N = u + N and a i denotes the piecewise linear approximation of the material inhomogeneity function. We remark that the jump of the displacement is not taken into account in the sum, and the gradient of u is approximated with finite difference of first order. The Dirichlet condition is applied on ∂Ω = {0, 2l} and the external displacement is chosen as u(0, t) = 0, u(2l, t) = 2lt. To enforce the boundary condition in the minimization process, we add it to the energy functional as a penalization term. Hence, we solve the following unconstrained minimization problem where the operator A ∈ R (2N +1)×(2N +1) is given by We consider the three different penalizations given by the τ , τ ∈]0, 1[, the SCAD, and the MCP penalties. Note that KerA = 0, hence assumptions (a)(ii) and (c)(ii) are satisfied and existence of a minimizer for (29) is guaranteed.
Our numerical experiments were conducted with a discretization in 2N intervals, N = 100. The time step, in the time discretization of [0, T ], with T = 3, is set to dt = 0.01. The parameters of the energy functional J h (u h ) are set to λ = 1, γ = 50.
We remark that in the following experiments, the material function a(x) was always chosen as the identity. For tests with more general a(x), we refer to the two-dimensional experiments reported in the following subsection. In Figs. 1 and 2, we report our results obtained by Algorithm 1, respectively, for the models p and SCAD. In each figure, We observe the three phases that we expect from a cohesive fracture model: -Pure elastic deformation: in this case, the jump amplitude is zero and the gradient of the displacement is constant in Ω\Γ ; -Prefracture: the two lips of the fracture do not touch each other, but they are not free to move. The elastic energy is still present. -Fracture: the two parts are free to move. In this final phase, the gradient of the displacement (and then the elastic energy) is zero.
We point out that the model we consider describes a fracture process in which the crack amplitude u is allowed to differ from zero. In particular, there is no possibility of breaking of the material with the crack surfaces remaining in contact. We remark that the formation of the crack is anticipated for smaller values of τ . As we see in Fig. 1, for τ = .01, prefracture and fracture are reached at t = .3 and t = 1.5, respectively. As τ is increased to τ = .1, prefracture and fracture occur at t = 1 and t = 3, respectively. We observe the same phenomenon for the SCAD (see Fig. 2).
We tested our algorithm also for the MCP model, where no prefracture phase can be observed, that is, the displacement breaks almost instantaneously to reach the complete fracture.
Finally, we remark that in our experiments, the residue is O(10 −16 ) and the number of iterations is small, e.g. 12, 15 for τ = .01, .1, respectively.

Two-Dimensional Experiments
We consider the two-dimensional domain Ω =]0, 1[×]0, 1] and we chose the onedimensional subdomain 0.5×]0, 1[ as the line of crack. We proceed in the discretization of the problem analogously as in Sect. 3.1.1, that is, we divide [0, 1] into 2N intervals and approximate the displacement function with a function u h , that is, piecewise linear in Ω\Γ and has two degrees of freedom on Γ to represent correctly the two lips of the fracture. Define the operator A ∈ R 3m(m−1)×m 2 where m = 2N + 2 and R 1 ∈ R (m(m−1))×m(m−1) and R 2 ∈ R m(m−2)×m(m−2) are two diagonal operators approximating the degree of homogeneity of the material, D 2 ∈ R m(m−2)×m 2 is defined as where G 1 , G 2 ∈ R m(m−1)×m 2 are defined as follows G 1 = kron(I m , D), G 2 = kron(D, I m ) and D =: diag(−ones(m, 1)) + diag(ones(m − 1, 1), 1) ∈ R (m−1)×m without the last row. Again, we enforce the boundary condition by adding it to the energy functional as a penalization term. Hence, we solve the following unconstrained minimization problem min |Au h − b| 2 2 + θ(D f u), (30) where b ∈ R 3m(m−1) in (30)  We perform two different series of experiments with boundary data, respectively, resulting from evaluating g 1 , x 2 ) ∈ Ω and the other one with boundary datum g 2 (t)(x) = 2t cos(4(x 2 − 0.5))( Figs. 3, 4, 5, and 6, we show the results obtained with boundary datum g 1 for each of the considered models, that is, τ , SCAD, and MCP and in Fig. 7, the ones with boundary datum g 2 for the τ model. In the case of boundary datum g 2 , we tested our algorithm also on the SCAD and the MCP models, obtaining similar results to the ones shown in Fig. 7. In these first experiments, the diagonal operators R 1 , R 2 are taken as the identity, that is, we suppose to have an homogeneous material. As expected from a cohesive fracture model, we observe the three phases of pure elastic deformation, prefracture, and fracture (see Sect. 3.1.1 for an explanation of the model and the three phases). When the boundary datum is g 2 , that is, not constant in the y direction, we note that the fracture is reached before in the part of the fracture line corresponding to the part of the boundary where the datum is bigger.
In Figs. 8, 9, and 10, we tested the algorithm in case of a non-homogeneous material. In Fig. 8, we show the result for a two-material specimen, that is, we took Note that for the above values of R 1 , R 2 , the slides of the specimen show an asymmetric behaviour, namely the displacement is flatter where the material function is bigger (that is, when R ii (x) = 600). In Figs. 9 and 10, we report the results when R 1 , R 2 are the discretization of the following function r (x, y) = 400ex p(y), for x ≤ N r (x, y) = 400y otherwise (33) Note that in Fig. 10, the boundary datum is chosen as g 3 (t) = 1 100 cos(2(y − 0.5))(x − 0.5). As expected due to the choice of R 1 , R 2 , we remark an asymmetric behaviour of the fracture in the y direction, namely the specimen brakes before where the material function is higher.

Comparison with GIST
In this section, we present the result of experiments to compare the performance of Algorithm 1 with the following two other algorithms for nonconvex and nonsmooth minimization. We first compare with the GIST "General Iterative Shrinkage and Thresholding" algorithm for τ , τ < 1 minimization. We took advantage of the fact that for GIST 1 an open source toolbox is available, which facilitated an unbiased comparison. Moreover, in [25], several tests were made to compare GIST and IRLS "Iteratively reweighted least squares", showing that the two algorithms have nearly the same performance, with only significant difference in speed, where GIST appears to be the faster one.
We remark that the results of [25] show no particular differences in the performance of the algorithm for different values of τ , except that the speed becomes much worse  for p near to 1, say τ = 0.9. Motivated also by this observations, the comparisons explained in the following were made for one fixed value of τ .
The comparison is carried out through the following three examples, the academical M-matrix problem, an optimal control problem, and a microscopy imaging reconstruction example.
The monotone algorithm is stopped when the ∞ -residue of the optimality condition (9) is of the order of 10 −3 in the M-matrix and optimal control problems and of the order of 10 −8 in the imaging example. GIST is terminated if the relative change of the two consecutive objective function values is less than 10 −5 or the number of iterations exceeds 1000. We remark that no significant changes were remarked by setting a lower tolerance than 10 −5 or a bigger number of maximal iteration for GIST.
Since both GIST and the FISTA solve the problem (1) when the operator Λ coincides with the identity, we also make this choice in the following subsections. Finally, we remark that the three examples were analysed already in [20] with different aims.

M-Matrix Example
We consider min x∈R n×n A is the forward finite difference gradient A = G 1 , G 2 , with G 1 := I ⊗ D ∈ R n(n+1)×n 2 and G 2 := D ⊗ I ∈ R n(n+1)×n 2 , I is the n ×n identity matrix, ⊗ the tensor product, D = (n + 1)D, andD := diag(ones(n + 1, 1)) + diag(−ones(n, 1) − 1) ∈ R (n+1)×n , without the last column. Then, A T A is an M matrix coinciding with the 5point star discretization on a uniform mesh on a square of the Laplacian with Dirichlet boundary conditions. Moreover, (34) can be equivalently expressed as min x∈R n×n where f = A T b. If λ = 0, this is the discretized variational form of the elliptic equation For λ > 0, the variational problem (35) gives a sparsity-enhancing solution for the elliptic equation (36), that is, the displacement y will be 0 when the forcing f is small. Our tests are conducted with f chosen as discretization of f = 10x 1 sin(5x 2 )cos(7x 1 ). The initialization is chosen as the solution of the corresponding non-sparse optimization problem.
We remark that in [37] and [20], the algorithm was also tested in the same situation for different values of τ and λ, showing, in particular and consistent with our expectations that the sparsity of the solution increases with λ.
Here, we focus on the comparison between the performances of Algorithm 1 and GIST. In order to compare the two schemes, we focus on the value of the unregularized functional J in (34) reached by both algorithms, the time to acquire it, and the number of iterations. Our tests were conducted for τ = 0.5, and λ incrementally increasing from 10 −3 to 0.3. The parameter ε was decreased from 10 −1 to 10 −6 . We report the values in Table 1 for λ = 0.05, since for the other values of λ, the results we obtained are comparable. We observe that Algorithm 1 achieves always lower values of the functional J, but in a longer time. The number of iterations needed by Algorithm 1 is smaller than the number of iterations of GIST for small values of λ, more precisely for λ < 0.1. Note that for smaller λ the number of iterations of Algorithm 1 is smaller than the one of GIST. This suggests, consistent with our expectation, that the monotone scheme is slower than GIST mainly because it solves a nonlinear equation at each iteration.
We carried out a further test in order to measure the timing performance of Algorithm 1, that is, the algorithm is stopped as soon as the value of J achieved by GIST is reached. In Table 1, we report the time, the number of iterations, the values of J, and the value of ε reached. We observe that the time is almost always smaller than the one of GIST.

Optimal Control Problem
We consider the linear control system d dt y(t) = Ay(t) + Bu(t), y(0) = 0, that is, where the linear closed operator A generates a C 0 -semigroup e At , t ≥ 0 on the Hilbert space X . More specifically, we consider the one-dimensional controlled heat equation for y = y(t, x): with homogeneous boundary conditions y(t, 0) = y(t, 1) = 0 and thus X = L 2 (]0, 1[). The differential operator Ay = y x x is discretized in space by the secondorder finite difference approximation with n = 49 interior spatial nodes (Δx = 1 50 ). We use two time-dependent controls − → u = (u 1 , u 2 ) with corresponding spatial control distributions b i chosen as step functions: 6,.7[ . The control problem consists in finding the control function − → u that steers the state y(0) = 0 to a neighbourhood of the desired state y d at the terminal time T = 1. We discretize the problem in time by the mid-point rule, i.e.
where − → u = (u 1 1 , . . . , u m 1 , u 1 2 , . . . u m 2 ) is a discretized control vector whose coordinates represent the values at the mid-point of the intervals (t k , t k+1 ). Note that in (39), we denote by B a suitable rearrangement of the matrix B in (37) with some abuse of notation. A uniform step-size Δt = 1 50 (m = 50) is utilized. The solution of the control problem is based on the sparsity formulation (1), where Λ = I and φ λ,τ (x) = λ|x| τ and b in (1) is the discretized target function chosen as the Gaussian distribution y d (x) = 0.4 exp(−70(x − .7) 2 )), centred at x = .7. That is, we apply our algorithm for the discretized optimal control problem in time and space where x from (1) is the discretized control vector u ∈ R 2m , which is mapped by A to the discretized output y at time 1 by means of (39). Moreover, b from (1) is the discretized state y d with respect to the spatial grid Δx. The parameter ε was initialized with 10 −3 and decreased down to 10 −8 .
Similarly as in the previous subsection, we compare the values of the functional, the time and the number of iterations. The experiments are carried out for τ = 0.5 and λ in the interval 10 −3 -0.2. We report only the values for the second control u 2 since the first control u 1 is always zero (as expected).
As can be seen from Table 2, the same kind of remarks as in the previous subsection apply. In particular, GIST is faster but less precise than Algorithm 1, but Algorithm 1 overcomes the value reached by GIST more rapidly. Note that we reported again only the results we obtained for the two values of λ = 0.001 and λ = 0.01 since for the other values of λ tested, we got comparable results.

Compressed Sensing Approach for Microscopy Image Reconstruction
We compare Algorithm 1 and GIST in a microscopy imaging problem, in particular we focus on the STORM (stochastic optical reconstruction microscopy) method, based on stochastic switching and high-precision detection of single molecules to achieve an image resolution beyond the diffraction limit. The literature on the STORM has been intensively increasing, see e.g. [38][39][40][41]. We refer in particular to [20] for a detailed description of the method and for more references.
Our approach is based on the following constrained-minimization problem: where τ ∈]0, 1], x is the up-sampled, reconstructed image, b is the experimentally observed image, and A is the impulse response (of size m × n, where m and n are the numbers of pixels in b and x, respectively). A is usually called the point spread First, we tested the procedure for same resolution images, in particular the conventional and the true images are both 128 × 128 pixel images. Then, the algorithm was tested in the case of a 16 × 16 pixel conventional image and a 128 × 128 true image. The values for the impulse response A and the measured data b were chosen according to the literature, in particular A was taken as the Gaussian PSF matrix with variance σ = 8 and size 3 × σ = 24, and b was simulated by convolving the impulse response A with a random 0-1 mask over the image adding a white random noise so that the signal to noise ratio is .01.
We carried out several tests with the same data for different values of τ, λ. We report only our results for τ = .1 and λ = 10 −6 , λ = 10 −9 for the same and the different resolution case, respectively, since for these values the best reconstructions were achieved. We focus on two different types of images, a sparse 0-1 cross-like image and the standard phantom image. In order to compare the performance of Algorithm 1 and the GIST algorithm, we focus on the number of surplus emitters (Error+) and missed emitters (Error−) recovered in the case of the cross image and different resolution. The errors are computed on an average over six recoveries for different values of the noise. The graphics of the errors against the noise are reported in Figs. 11 and 12 for Algorithm 1 and GIST, respectively. We remark that these quantities are typically used as a measure of the efficacy of the reconstruction method, see for example [42] (where, under certain conditions, a linear decay with respect to the noise is proven) and [43]. The results shows that by GIST the Error− is always 197, whereas by Algorithm 1 is always under 53 and even smaller for small values of the noise. On the other hand, the Error+ by GIST is always 0 and by Algorithm 1 is zero for small values of the noise and then monotonically increasing until it reaches 175 when the noise is equal to 0.1. Consistently with what expected, by Algorithm 1, the graphics show a linear decay w.r.t. the noise, differently from the behaviour showed by GIST. Moreover, the results found by Algorithm 1 lead to more accuracy in the recovery, in the sense that the quantity of missed emitters is smaller, whereas on the other hand, GIST seems to lead to a more sparser solutions (since the Error+ is 0 by GIST).
Finally, we remark that in the case of the cross image, GIST is faster than our algorithm, consistently with the result presented in the previous subsection and as expected, since our algorithm solves a nonlinear equation for each minimization problem. On the other hand, in the case of the standard phantom image, GIST results to be far slower than Algorithm 1.
In Fig. 12, we report the results obtained in the same situation by the FISTA "Fast Iterative Shrinkage-Thresholding Algorithm" for 1 minimization. We remark that by the FISTA, the Error+ is always above 400, whereas by Algorithm 1 is zero for small value of the noise. This shows that Algorithm 1 leads to more sparsity with respect to the FISTA, consistently with our expectation since the FISTA is based on 1 minimization.

Perspectives and Open Problems
An open problem of interest to us is the study of problems like (1) for the case where the linear mapping A is replaced by a nonlinear, smooth operator f : R n → R m . One of the motivations arises from control of nonlinear dynamical system. We could proceed by iteratively applying the monotone scheme in Sect. 2.3 to auxiliary problems arising from linearization of f at the current iterate and by updating the sequence such obtained in an outer loop. We are particularly interested in the case in which the nonlinear operator f is nonconvex. From a fracture mechanics point of view, this would mean considering not only small-strain energy as in the current paper, but possibly polyconvex strain energy as in [44]. Considering a nonconvex energy would be more consistent from a mechanical point of view, and in particular in line with Coleman-Noll's theorem.

Conclusions
We have developed a monotone convergent algorithm for a class of nonconvex nonsmooth optimization problems arising in the modelling of fracture mechanics and in imaging reconstruction, including the τ , τ ∈]0, 1], the smoothly clipped absolute deviation and the minimax concave penalty. Theoretically, we established the existence of a minimizer of the original problem under assumptions implying coercivity of the functional. Then, we derived necessary optimality conditions for a regularized version of the original problem. The optimality conditions for the regularized problem were solved through a monotonically convergent scheme based on an iterative procedure. We proved the convergence of the iteration procedure under the same assumptions that guarantee existence. A remarkable result is the strict monotonicity of the functional along the sequence of iterates generated by the scheme. Moreover, we proved the convergence of the regularized problem to the original one, as the regularization parameter goes to zero.
The procedure is very efficient and accurate. The efficiency and accuracy of the procedure was verified by numerical tests simulating the evolution of cohesive fractures and microscopy imaging. An issue of high relevance to us was the comparison of the scheme to two alternative algorithms, the GIST "General Iterative Shrinkage and Thresholding" algorithm for τ minimization, with τ strictly positive and less than 1 and the FISTA "Fast Iterative Shrinkage-Thresholding Algorithm" for 1 minimization. We first compared with GIST by focusing on the infimal value reached by the iteration procedure and on the computing time. Our results showed that the monotone algorithm is able to reach a smaller value of the objective functional when compared to GIST's, therefore leading to a better accuracy. Finally we compared our scheme with FISTA in sparse recovery related to microscopy imaging. The results showed that the monotone scheme leads to more sparsity with respect to FISTA, as expected since FISTA concerns 1 minimization.