Accelerating the DC algorithm for smooth functions

We introduce two new algorithms to minimise smooth difference of convex (DC) functions that accelerate the convergence of the classical DC algorithm (DCA). We prove that the point computed by DCA can be used to define a descent direction for the objective function evaluated at this point. Our algorithms are based on a combination of DCA together with a line search step that uses this descent direction. Convergence of the algorithms is proved and the rate of convergence is analysed under the {\L}ojasiewicz property of the objective function. We apply our algorithms to a class of smooth DC programs arising in the study of biochemical reaction networks, where the objective function is real analytic and thus satisfies the {\L}ojasiewicz property. Numerical tests on various biochemical models clearly show that our algorithms outperforms DCA, being on average more than four times faster in both computational time and the number of iterations. The algorithms are globally convergent to a non-equilibrium steady state of a biochemical network, with only chemically consistent restrictions on the network topology.


Introduction
Many problems arising in science and engineering applications require the development of algorithms to minimise a nonconvex function. If a nonconvex function admits a decomposition, this may be exploited to tailor specialised optimisation algorithms. Our main focus is the following optimisation problem where f 1 , f 2 : R m → R are continuously differentiable convex functions and In our case, as we shall see in Section 4, this problem arises in the study of biochemical reaction networks. In general, φ is a nonconvex function. The function in problem (1) belongs to two important classes of functions: the class of functions that can be decomposed as a sum of a convex function and a differentiable function (composite functions) and the class of functions that are representable as difference of convex functions (DC functions).
In 1981, M. Fukushima and H. Mine [6,16] introduced two algorithms to minimise a composite function. In both algorithms, the main idea is to linearly approximate the differentiable part of the composite function at the current point and then minimise the resulting convex function to find a new point. The difference between the new and current points provides a descent direction with respect to the composite function, when it is evaluated at the current point. The next iteration is then obtained through a line search procedure along this descent direction. Algorithms for minimising composite functions have been extensively investigated and found applications to many problems such as: inverse covariance estimate, logistic regression, sparse least squares and feasibility problems, see e.g. [8,13,14,18] and the references quoted therein.
In 1986, T. Pham Dinh and S. El Bernoussi [23] introduced an algorithm to minimise DC functions. In its simplified form, the Difference of Convex functions Algorithm (DCA) linearly approximates the concave part of the objective function (− f 2 in (1)) about the current point and then minimises the resulting convex approximation to the DC function to find the next iteration, without recourse to a line search. The main idea is similar to Fukushima-Mine approach but was extended to the non-differentiable case. This algorithm has been extensively studied by H.A. Le Thi, T. Pham Dinh and their collaborators, see e.g. [11,12,22,9]. DCA has been successfully applied in many fields, such as machine learning, financial optimisation, supply chain management and telecommunication [5,21,9]. Nowadays, DC programming plays an important role in nonconvex programming and DCA is commonly used because of its key advantages: simplicity, inexpensiveness and efficiency [9]. Some results related to the convergence rate for special classes of DC programs have been also established [10,12].
In this paper we introduce two new algorithms to find stationary points of DC programs, called Boosted Difference of Convex function Algorithms (BDCA), which accelerate DCA with a line search using an Armijo type rule. The first algorithm directly uses a backtracking technique, while the second uses a quadratic interpolation of the objective function together with backtracking. Our algorithms are based on both DCA and the proximal point algorithm approach of Fukushima-Mine. First, we compute the point generated by DCA. Then, we use this point to define the search direction. This search direction coincides with the one employed by Fukushima-Mine in [6]. The key difference between their method and ours is the starting point used for the line search: in our algorithms we use the point generated by DCA, instead of using the previous iteration. This scheme works thanks to the fact that the defined search direction is not only a descent direction for the objective function at the previous iteration, as observed by Fukushima-Mine, but is also a descent direction at the point generated by DCA.
Moreover, it is important to notice that the iterations of Fukushima-Mine and BDCA never coincide, as the largest step size taken in their algorithm is equal to one (which gives the DCA iteration). In fact, for smooth functions, the iterations of Fukushima-Mine usually coincide with the ones generated by DCA, as the step size equal to one is normally accepted by their Armijo rule.
We should point out that DCA is a descent method without line search. This is something that is usually claimed to be advantageous in the large-scale setting. Our purpose here is the opposite: we show that a line search can increase the performance even for high-dimensional problems.
Further, we analyse the rate of convergence under the Łojasiewicz property [15] of the objective function. It should be mentioned that the Łojasiewicz property is recently playing an important role for proving the convergence of optimisation algorithms for analytic cost functions, see e.g. [1,2,3,12].
We have performed numerical experiments in functions arising in the study of biochemical reaction networks. We show that the problem of finding a steady state of these networks, which plays a crucial role in the modelling of biochemical reaction systems, can be reformulated as a minimisation problem involving DC functions. In fact, this is the main motivation and starting point of our work: when one applies DCA to find a steady state of these systems, the rate of convergence is usually quite slow. As these problems commonly involve hundreds of variables (even thousands in the most complex systems, as Recon 2 1 ), the speed of convergence becomes crucial. In our numerical tests we have compared BDCA and DCA for finding a steady state in various biochemical network models of different size. On average, DCA needed five times more iterations than BDCA to achieve the same accuracy, and what is more relevant, our implementation of BDCA was more than four times faster than DCA to achieve the same accuracy. Thus, we prove both theoretically and numerically that BDCA results more advantageous than DCA. Luckily, the objective function arising in these biochemical reaction networks is real analytic, a class of functions which is known to satisfy the Łojasiewicz property [15]. Therefore, the above mentioned convergence analysis results can be applied in this setting.
The rest of this paper is organised as follows. In Section 2, we recall some preliminary facts used throughout the paper and we present the main optimisation problem. Section 3 describes our main results, where the new algorithms (BDCA) and their convergence analysis for solving DC programs are established. A DC program arising in biochemical reaction network problems is introduced in Section 4. Numerical results comparing BDCA and DCA on various biochemical network models are reported in Section 5. Finally, conclusions are stated in the last section.

Preliminaries
Throughout this paper, the inner product of two vectors x, y ∈ R m is denoted by x, y , while · denotes the induced norm, defined by x = x, x . The nonnegative orthant in R m is denoted by R m + = [0, ∞) m and B(x, r) denotes the closed ball of center x and radius r > 0. The gradient of a differentiable function f : Recall that a function f : R m → R is said to be convex if for all x, y ∈ R m and λ ∈ (0, 1).
x−y 2 for all x, y ∈ R m and λ ∈ (0, 1), On the other hand, a function F : R m → R m is said to be monotone when Further, F is called strongly monotone with modulus σ > 0 when F is called locally Lipschitz continuous if for every x in R m , there exists a neighbourhood U of x such that F restricted to U is Lipschitz continuous.
We have the following well-known result. To establish our convergence results, we will make use of the Łojasiewicz property, defined next.
where we adopt the convention 0 0 = 1. The constant θ is called Łojasiewicz exponent of f atx.
(ii) The function f is said to be real analytic if for every x ∈ R n , f may be represented by a convergent power series in some neighbourhood of x.
Problem (1) can be easily transformed into an equivalent problem involving strongly convex functions. Indeed, choose any ρ > 0 and consider the functions g(x) Then g and h are strongly convex functions with modulus ρ and g(x) − h(x) = φ (x), for all x ∈ R m . In this way, we obtain the equivalent problem The key step to solve (P) with DCA is to approximate the concave part −h of the objective function φ by its affine majorisation and then minimise the resulting convex function. The algorithm proceeds as follows. 2. Solve the strongly convex optimisation problem to obtain the unique solution y k .
3. If y k = x k then STOP and RETURN x k , otherwise set x k+1 := y k , set k := k + 1, and go to Step 2.
In [6] Fukushima and Mine adapted their original algorithm reported in [16] by adding a proximal term ρ 2 x − x k 2 to the objective of the convex optimisation subproblem. As a result they obtain an optimisation subproblem that is identical to the one in Step 2 of DCA, when one transforms (1) into (4) by adding ρ 2 x 2 to each convex function. In contrast to DCA, Fukushima-Mine algorithm [6] also includes a line search along the direction d k := y k − x k to find the smallest nonnegative integer l k such that the Armijo type rule one has x k+1 = y k and the iterations of both algorithms coincide. As we shall see in Proposition 3.1, this is guaranteed to happen if α ≤ ρ.

Boosted DC Algorithms
Let us introduce our first algorithm to solve (P), which we call a Boosted DC Algorithm with Backtracking. The algorithm is a combination of Algorithm 1 and the algorithm of Fukushima-Mine [6].

Solve the strongly convex minimisation problem
to obtain the unique solution y k .
3. Set d k := y k − x k . If d k = 0, STOP and RETURN x k . Otherwise, go to Step 4.
5. Set x k+1 := y k +λ k d k . If x k+1 = x k then STOP and RETURN x k , otherwise set k := k +1, and go to Step 2.
The next proposition shows that the solution of (P k ), which coincides with the DCA subproblem in Algorithm 1, provides a decrease in the value of the objective function. For the sake of completeness, we include its short proof.
Proof. Since y k is the unique solution of the strongly convex problem (P k ), we have which implies On the other hand, the strong convexity of h implies Adding the two previous inequalities, we have which implies (6).
If λ k = 0, the iterations of BDCA-Backtracking coincide with those of DCA, since the latter sets x k+1 := y k . Next we show that d k = y k − x k is a descent direction for φ at y k . Thus, one can achieve a larger decrease in the value of φ by moving along this direction. This simple fact, which permits an improvement in the performance of DCA, constitutes the key idea of our algorithms.
that is, d k is a descent direction for φ at y k .
Proof. The function h is strongly convex with constant ρ. This implies that ∇h is strongly monotone with constant ρ; whence, Further, since y k is the unique solution of the strongly convex problem (P k ), we have and completes the proof.
As a corollary, we deduce that the backtracking Step 4 of Algorithm 2 terminates finitely when ρ > α.
Proof. If d k = 0 there is nothing to prove. Otherwise, by the mean value theorem, there is somet ∈ (0, 1) such that As ∇φ is continuous at y k , there is some δ > 0 such that and the proof is complete.
Therefore, Algorithm 2 uses the same direction as the Fukushima-Mine algorithm [6], where x k+1 = x k + β l d k = β l y k + 1 − β l x k for some 0 < β < 1 and some nonnegative integer l. The iterations would be the same if β l = λ + 1. Nevertheless, as 0 < β < 1, the step size λ = β l − 1 chosen in the Fukushima-Mine algorithm [6] is always less than or equal to zero, while in Algorithm 2, only step sizes λ ∈ ]0,λ ] are explored. Moreover, observe that the Armijo type rule (5), as used in [6], searches for an l k such that φ (x k + β l k d k ) < φ (x k ), whereas Algorithm 2 searches for a λ k such that φ (y k + λ k d k ) < φ (y k ). We know from (6) and (9) that thus, Algorithm 2 results in a larger decrease in the value of φ at each iteration than DCA, which sets λ := 0 and x k+1 := y k . Therefore, a faster convergence of Algorithm 2 compared with DCA is expected, see Figure 1 and Figure 3.
The following convergence results were inspired by [2], which in turn were adapted from the original ideas of Łojasiewicz; see also [4, Section 3.2].

Proposition 3.3.
For any x 0 ∈ R m , either Algorithm 2 returns a stationary point of (P) or it generates an infinite sequence such that the following holds.
(i) φ (x k ) is monotonically decreasing and convergent to some φ * .
(ii) Any limit point of {x k } is a stationary point of (P). If in addition, φ is coercive then there exits a subsequence of {x k } which converges to a stationary point of (P).
Proof. Because of (7), if Algorithm 2 stops at Step 3 and returns x k , then x k must be a stationary point of (P). Otherwise, by Proposition 3.1 and Step 4 of Algorithm 2, we have Hence, as the sequence {φ (x k )} is monotonically decreasing and bounded from below by (2), it converges to some φ * , which proves (i). Consequently, we have Thus, by (10), one has d k Taking the limit as i → ∞ in (7), as ∇h and ∇g are continuous, we have ∇h(x) = ∇g(x). If φ is coercive, since the sequence {φ (x k )} is convergent, then the sequence {x k } is bounded. This implies that there exits a subsequence of {x k } converging tox, a stationary point of (P), which proves (ii).
To prove (iii), observe that (10) implies that Summing this inequality from 0 to N, we obtain whence, taking the limit when N → ∞, and the proof is complete.
We will employ the following useful lemma to obtain bounds on the rate of convergence of the sequences generated by Algorithm 2. This result appears within the proof of [ Then (i) if α = 0, the sequence {s k } converges to 0 in a finite number of steps; (ii) if α ∈ (0, 1], the sequence {s k } converges linearly to 0 with rate 1 − 1 β ; (iii) if α > 1, there exists η > 0 such that Proof. If α = 0, then (13) implies and (i) follows. Assume that α ∈ (0, 1]. Since s k → 0, we have that s k < 1 for all k large enough. Thus, by (13), we have Therefore, s k+1 ≤ 1 − 1 β s k ; i.e., {s k } converges linearly to 0 with rate 1 − 1 β . Suppose now that α > 1. If s k = 0 for some k, then (13) implies s k+1 = 0. Then the sequence converges to zero in a finite number of steps, and thus (iii) trivially holds. Hence, we will assume that s k > 0 and that (13) holds for all k ≥ N, for some positive integer N. Consider the decreasing function ϕ : (0, +∞) → R defined by ϕ(s) := s −α . By (13), for k ≥ N, we have As α − 1 > 0, this implies that for all k ≥ N. Thus, summing for k from N to j − 1 ≥ N, we have which gives, for all j ≥ N + 1, Therefore, there is some η > 0 such that which completes the proof. (ii) if θ ∈ 0, 1 2 then the sequences {x k } and {φ (x k )} converge linearly to x * and φ * , respectively; (iii) if θ ∈ 1 2 , 1 then there exist some positive constants η 1 and η 2 such that for all large k.
Proof. By Proposition 3.3, we have lim k→∞ φ (x k ) = φ * . If x * is a cluster point of {x k }, then there exists a subsequence {x k i } of {x k } that converges to x * . By continuity of φ , we have that Hence, φ is finite and has the same value φ * at every cluster point of Therefore, x k = x k+p for all p ≥ 0 and Algorithm 2 terminates after a finite number of steps. From now on, we assume that φ (x k ) > φ * for all k. As φ satisfies the Łojasiewicz property, there exist M > 0, ε 1 > 0 and θ ∈ [0, 1) such that Further, as ∇g is locally Lipschitz around x * , there are some constants L ≥ 0 and ε 2 > 0 such that ∇g(x) − ∇g(y) ≤ L x − y , ∀x, y ∈ B(x * , ε 2 ).
From (17), as λ k ∈ (0,λ ], we deduce for all k ≥ N such that x k ∈ B(x * , ε). We prove by induction that x k ∈ B(x * , ε) for all k ≥ N. Indeed, from (16) the claim holds for k = N. We suppose that it also holds for k = N, N + 1, . . ., N + p − 1, with p ≥ 1. Then (20) is valid for k = N, N + 1, . . ., N + p − 1. Therefore where the last inequality follows from (16). Adding (20) from k = N to P one has Taking the limit as P → ∞, we can conclude that This means that {x k } is a Cauchy sequence. Therefore, since x * is a cluster point of {x k }, the whole sequence {x k } converges to x * . By Proposition 3.3, x * must be a stationary point of (P).
We know that s i := ∑ ∞ k=i x k+1 − x k is finite by (22). Notice that x i − x * ≤ s i by the triangle inequality. Therefore, the rate of convergence of x i to x * can be deduced from the convergence rate of s i to 0. Adding (20) from i to P with N ≤ i ≤ P, we have Then by (14) and (15), we get Hence, taking K 2 := MLK By applying Lemma 3.1 with α := θ 1−θ and β := K 2 , we see that the statements in (i)-(iii) regarding the sequence {x k } hold.
The iteration given by DCA (Algorithm 1) satisfies On the other hand, the iteration defined by Algorithm 2 is , we have φ ( x 1 ) < φ (x 1 ). The optimal step size is attained at λ opt = 25 24 with x 1 = 1, which is the global minimiser of φ .
As shown in Proposition 3.2, the function is decreasing at 0. The value λ = 0 corresponds to the next iteration chosen by DCA, while the next iteration chosen by algorithm [6] sets λ ∈ ] − 1, 0]. Algorithm 2 chooses the next iteration taking λ ∈ ]0,λ ], which permits to achieve an additional decrease in the value of φ . Here, the optimal value is attained at λ opt = 25 24 ≈ 1.04.
Observe in Figure 1 that the function behaves as a quadratic function nearby 0. Then, a quadratic interpolation of this function should give us a good candidate for choosing a step size close to the optimal one. Whenever ∇φ is not too expensive to compute, it makes sense to construct a quadratic approximation of φ with an interpolation using three pieces of information: φ k (0) = φ (y k ), φ ′ k (0) = ∇φ (y k ) T d k and φ k (λ ). This gives us the quadratic function see e.g. [19,Section 3.5]. When φ k (λ ) > φ k (0) +λ φ ′ k (0), the function ϕ k has a global minimiser at This suggests the following modification of Algorithm 2.
2. Solve the strongly convex minimisation problem to obtain the unique solution y k .
3. Set d k := y k − x k . If d k = 0 STOP and RETURN x k . Otherwise, go to Step 4.
5. Set x k+1 := y k +λ k d k . If x k+1 = x k then STOP and RETURN x k , otherwise set k := k +1, and go to Step 2.
Corollary 3.2. The statements in Theorem 3.1 also apply to Algorithm 3.
Proof. Just observe that the proof of Theorem 3.1 remains valid as long as the step sizes are bounded above by some constant and below by zero. Algorithm 3 uses the same directions than Algorithm 2, and the step sizes chosen by Algorithm 3 are bounded above by λ max and below by zero.
Another option here would be to construct a quadratic approximation ψ k using φ k (−1) = φ (x k ) instead of φ k (λ ). This interpolation is computationally less expensive, as it does not require the computation of φ k (λ ). Nevertheless, our numerical tests for the functions in Section 4 show that this approximation usually fits the function φ k more poorly. In particular, this situation occurs in Example 3.1, as shown in Figure 2.  Figure 2: Plots of the quadratic interpolations ϕ 0 and ψ 0 for the function φ 0 (λ ) = φ (y 0 + λ d 0 ) from Example 3.1, withλ = 2. Note that ψ 0 (λ ) fits poorly φ 0 (λ ) for λ > 0.
One could also construct a cubic function that interpolates φ k (−1), φ k (0), φ ′ k (0) and φ k (λ ), see [19,Section 3.5]. However, for the functions in Section 4, we have observed that this cubic function usually fits the function φ k worse than the quadratic function ϕ k in (24).

Remark 3.2. Observe that Algorithm 2 and Algorithm 3 still work well if we replace
Step 2 by the following proximal step as in [17] (P k ) minimise for some positive constants c k .
Example 3.2. (Finding zeroes of systems of DC functions) Suppose that one wants for find a zero of a system of equations where p : R m → R m + and c : R m → R m + are twice continuously differentiable functions such that p i : R m → R + and c i : R m → R + are convex functions for all i = 1, . . ., m. Then, Observe that all the components of p(x) and c(x) are nonnegative convex functions. Hence, both f 1 (x) := 2 p(x) 2 + c(x) 2 and f 2 (x) := p(x) +c(x) 2 are continuously differentiable convex functions, because they can be expressed as a finite combination of sums and products of nonnegative convex functions. Thus, we can either apply DCA or BDCA in order to find a solution to (26) by setting φ (x) Suppose thatx is an accumulation point of the sequence {x k } generated by either Algorithm 2 or Algorithm 3, and assume that ∇ f (x) is nonsingular. Then, by Proposition 3.3, we must have ∇ f (x) f (x) = 0 m , which implies that f (x) = 0 m , as ∇ f (x) is nonsingular. Moreover, for all x close tox, we have where · also denote the induced matrix norm and M is an upper bound of 1 2 (∇ f (x)) −1 aroundx. Thus, φ has the Łojasiewicz property atx with exponent θ = 1 2 . Finally, for all ρ > 0, the function g(x) := f 1 (x) + ρ 2 x 2 is twice continuously differentiable, which in particular implies that ∇g is locally Lipschitz continuous. Therefore, either Theorem 3.1 or Corollary 3.2 guarantee the linear convergence of {x k } tox.

A DC problem in biochemistry
Consider a biochemical network with m molecular species and n reversible elementary reactions 2 . Define forward and reverse stoichiometric matrices, F, R ∈ Z m×n ≥0 , respectively, where F i j denotes the stoichiometry 3 of the i th molecular species in the j th forward reaction and R i j denotes the stoichiometry of the i th molecular species in the j th reverse reaction. We use the standard inner product in R m , i.e., x, y = x T y for all x, y ∈ R m . We assume that every reaction conserves mass, that is, there exists at least one positive vector l ∈ R m >0 satisfying (R−F) T l = 0 n [7] where R − F represents net reaction stoichiometry. We assume the cardinality 4 of each row of F and R is at least one, and the cardinality of each column of R − F is at least two, usually three. Therefore, R − F may be viewed as the incidence matrix of a directed hypergraph. The matrices F and R are sparse and the particular sparsity pattern depends on the particular biochemical network being modelled.
Let u ∈ R m >0 denote a variable vector of molecular species concentrations. Assuming constant nonnegative elementary kinetic parameters k f , k r ∈ R n ≥0 , we presume elementary reaction kinetics for forward and reverse elementary reaction rates as s(k f , u) := exp(ln(k f ) + F T ln(u)) and r(k r , u) := exp(ln(k r ) + R T ln(u)), respectively, where exp(·) and ln(·) denote the respective componentwise functions. Then, the deterministic dynamical equation for time evolution of molecular species concentration is given by Investigation of steady states plays a crucial role in the modelling of biochemical reaction systems. If one transforms (28) to logarithmic scale, by letting then, up to a sign, the right-hand side of (28) is equal to the function where [ · , · ] stands for the horizontal concatenation operator. Thus, we shall focus on finding the points x ∈ R m such that f (x) = 0, which correspond to the steady states of the dynamical equation (27). A pointx will be a zero of the function f if and only if f (x) 2 = 0. Denoting one obtains, as in Example 3.2, Again, as all the components of p(x) and c(x) are positive and convex functions 5 , both are convex functions. In addition to this, both f 1 and f 2 are smooth, having see e.g. [19, pp. 245-246], with where EXP (·) denotes the diagonal matrix whose entries are the elements in the vector exp (·).
, the problem of finding a zero of f is equivalent to the following optimisation problem: minimise We now prove that φ satisfies the Łojasiewicz property.
Since b i j are nonnegative integers for all i and j, we conclude that the function φ is real analytic (see Proposition 2.2.2 and Proposition 2.2.8 in [20]). It follows from Proposition 2.2 that the function φ satisfies the Łojasiewicz property with some exponent θ ∈ [0, 1). Finally, as in Example 3.2, for all ρ > 0, the function g(x) := f 1 (x) + ρ 2 x 2 is twice continuously differentiable, which implies that ∇g is locally Lipschitz continuous. Therefore, either Theorem 3.1 or Corollary 3.2 guarantee the convergence of the sequence generated by BDCA, as long as the sequence is bounded.
Remark 4.1. In principle, one cannot guarantee the linear convergence of BDCA applied to biochemical problems for finding steady states. Due to the mass conservation assumption, Therefore, the reasoning in Example 3.2 cannot be applied.
In Table 1 we report the numerical results comparing DCA and BDCA with quadratic interpolation (Algorithm 3) for 14 models arising from the study of systems of biochemical reactions. We only provide the numerical results for Algorithm 3 because it normally gives better results than Algorithm 2 for biochemical models, as it is shown in Figure 3. In Figure 4 we show a comparison of the rate of convergence of DCA and BDCA with quadratic interpolation for two big models. The parameters used were α = 0.4, β = 0.5, ρ = 100 andλ = 5.   Table 1: Performance comparison of BDCA and DCA for finding a steady state of various biochemical reaction network models. For each model, we selected a random kinetic parameter w ∈ [−1, 1] 2n and we randomly chose 10 initial points x 0 ∈ [−2, 2] m . For each x 0 , BDCA was run 1000 iterations, while DCA was run until it reached the same value of φ (x) as obtained with BDCA.

Concluding Remarks
In this paper, we introduce two new algorithms for minimising smooth DC functions, which we term Boosted Difference of Convex function Algorithms (BDCA). Our algorithms combine DCA together with a line search, which utilises the point generated by DCA to define a search direction. This direction is also employed by Fukushima-Mine in [6], with the difference that our algorithms start searching for the new candidate from the point generated by DCA, instead of starting from the previous iteration. Thus, our main contribution comes from the observation that this direction is not only a descent direction for the objective function at the previous iteration, as observed by Fukushima-Mine, but is also a descent direction at the point defined by DCA. Therefore, with the slight additional computational effort of a line search one can achieve a significant decrease in the value of the function. We prove that every cluster point of the algorithms are stationary points of the optimisation problem. Moreover, when the objective function satisfies the Łojasiewicz property, we prove global convergence of the algorithms and establish convergence rates. We demonstrate that the important problem of finding a steady state in the dynamical modelling of systems of biochemical reactions can be formulated as an optimisation problem involving a difference of convex functions. We have performed numerical experiments, using models of systems of biochemical reactions from various species, in order to find steady states.
The tests clearly show that our algorithm outperforms DCA, being able to achieve the same decrease in the value of the DC function while employing substantially less iterations and time. On average, DCA needed five times more iterations to achieve the same accuracy as BDCA. Furthermore, our implementation of BDCA was also more than four times faster than DCA. In fact, the slowest instance of BDCA was always at least three times faster than DCA. This substantial increase in the performance of the algorithms is especially relevant when the typical size of the problems is big, as is the case with all realistic biochemical network models.