We introduce reduced order model (ROM) based Hessian approximations for use in inexact Newton methods for the solution of large-scale smooth optimization problems min u∈R m J (u) = J (y(u), u), (1) where for given u ∈ R m the vector y(u) ∈ R n is the solution of c(y(u), u) = 0.
(
The optimization problem (1,2) often comes from a discretization of optimal control problems and in this case u ∈ R m is the discretized control, y ∈ R n is the discretized state, and (2) is the discretized state equation. Our approach can be easily extended to a setting where the control space U and the state space Y are Hilbert spaces, but because of space restrictions we limit ourselves to the finite dimensional case. In our applications, the state Eq. ( 2) is a discretized parabolic partial differential equation (PDE). In this case, each objective function (1) evaluation requires the solution of a discretized PDE, the evaluation of the objective function gradient requires the solution of the (linear) adjoint PDE, and the evaluation of a Hessian-times-vector multiplication requires the additional solution of two (linear) PDE. This is computationally expensive. In addition, the additional PDEs that have to be solved for gradient and Hessian-times-vector computations depend on the solution of y(u) ∈ R n of the state equation, which for time dependent PDEs is also memory intensive. ROMs can be used to reduce the computation time and memory requirements.
Previous approaches of using ROMs in optimization have approximated the mapping u → y(u) ∈ R n using a ROM applied to the state Eq. (2). Typically, the state Eq. ( 2) is approximated by a projection based ROM
where V ∈ R n×r is a matrix with rank r n. The ROM state solution is then used to approximate (1) by min u∈R m J (V y(u), u).
For many problems, including example problems in Sect. 4 of this paper, the ROM approximation (4,3) of the states implies via the optimality conditions that the controls u are also contained in a low dimensional subspace related to V.
In special cases, one can extend ROM approaches for dynamical systems [5] to find one ROM such that the solution of (4, 3) is a good approximation of the solution of (1,2). See, e.g., [2,3], and the survey [7]. For general problems, however, the ROM is only valid in a potentially small neighborhood of the current control u c and corresponding state y c = y(u c ). In this case trust-region based model management approaches have been proposed. See [1,10,13,16,19,20], and the surveys [7,17].
While in principle, one can use trust-region based model management approaches that generate a new ROM at each new iterate, in practice their use is limited by the computational cost of ROM generation versus the computational savings resulting from that ROM.
Therefore, instead of generating ROMs for the optimization problem, we generate ROM approximations of the Hessians. These ROMs are computationally cheaper than a ROM that also has to approximate the objective function (1) and its gradient over a range of controls. Ideally the ROM Hessians still generate optimization steps that are close to the Newton steps, and therefore lead to fast local convergence of the optimization algorithm at a fraction of the cost of a corresponding Newton-type methods applied to the full order model (FOM) problem (1,2).
In Sect. 2 we will outline our new approach and in Sect. 3 we will specify it further for a discretized parabolic optimal control problem. We demonstrate our approach on semilinear parabolic optimal control problems in Sect. 4.
This section provides a general outline of our proposed approach. We begin with a review of gradient and Hessian computation and the line-search Newton-Conjugate-Gradient (Newton-CG) method. Then we will present the general structure of our ROM Hessian approximation, and a heuristic for choosing the ROM subspace within this approximation. Inexact Newton Method and ROM Hessian Approximation. To make the definition of the objective function J and its gradient and Hessian calculation rigorous, we make the following assumptions throughout. Here we use ∇ y J (y, u) ∈ R n , ∇ u J (y, u) ∈ R m to denote the partial gradients and c y (y, u) ∈ R n×n , c u (y, u) ∈ R n×m to denote the partial Jacobians. Similar notation is used for second derivatives.
The gradient of the objective J in (1, 2) can be computed using the so-called adjoint equation approach. See, e.g., [12,Sect. 1.6]. Define the Lagrangian L(y, u, λ) := J (y, u) + λ T c(y, u) corresponding to (1,2). Given a control u and corresponding state y(u) the gradient of the objective J is computed by first solving the adjoint equation
for λ = λ(u) and then setting
Given a control u, the corresponding state y(u), and the corresponding adjoint λ(u), the Hessian-times-vector product ∇ 2 J (u) d is computed by solving the linearized state equation c y (y(u), u) w = c u (y(u), u) d (6a)
for w, then the second order adjoint equation
for p, and then computing
For optimal control problems, the computation of y(u) requires the solution of a nonlinear discretized PDE. The gradient computation requires the solution of a linear discretized adjoint PDE to compute λ(u). Each Hessian-times-vector operation ∇ 2 J (u) d requires the solution of two linear discretized PDEs to compute w and p, respectively.
The problem (1) is solved using a line-search Newton-CG Method. However, a trust-region method could be used as well, and the proposed ROM Hessian approximation provides the same computational benefits in a trust-region setting.
Given a current iterate u c the line-search Newton-CG Method [15, Algorithm 7.1] applies the truncated CG method to approximately solve the Newton subproblem
Then it computes a suitable step-size α c ∈ (0, 1] such that the sufficient decrease condition
is satisfied. The new iterate is given by u
The ROM Hessian approximation is computed by applying a projection ROM to the Eq. (6a, b). Let V ∈ R n×r , r n, with linearly independent columns such that V T c y (y(u), u)V is invertible.
The ROM Hessian-times-vector product is computed by first solving the projected linearized state equation
for w, then solving the projected second order adjoint equation
) for p, and then computing
The ROM Hessian approximation ( 9) is used to compute the direction in the inexact Newton method. Instead of solving (7), given a current iterate u c the direction d is computed by using the truncated CG method to approximately solve
The step size α c is computed as before using ( 8) with d replaced by d. The new iterate is given by u
The global convergence of the original line-search Newton method and the linesearch Newton method with ROM Hessian approximation can be proven using standard arguments. See, e.g., [15,Theorem 3.2]. While first order global convergence can be ensured under standard conditions if the ROM Hessian approximation is used, the speed with which the inexact Newton method converges in this case depends on the quality of the ROM Hessian approximation, and in particular on V. Next we motivate our choice of V.
Basic Properties of the ROM Hessian. Using (6) it can be seen that the Hessian of J is given by
where we have set y = y(u) and λ = λ(u) to simplify notation. We will use this simplified notation frequently during the remainder of this section.
Applying the Implicit Function Theorem to the state Eq. ( 2) shows that the sensitivity of u → y(u) is given by
Similarly, applying the Implicit Function Theorem to the adjoint Eq. (5a) shows that the derivative of u → λ(u) is given by
Using these derivatives, the Hessian can be written as
Using ( 9) it can be seen that the ROM Hessian approximation is given by
If we define
and
) then the ROM Hessian approximation can be written as
The next result ensures the positive definiteness of the Hessian and its ROM approximation in situations that are often encountered in classes of applications, such as those in Sect. 4.
and ∇ yu L(y, u, λ) = 0, then ∇ 2 J (u) and ∇ 2 J (u) are positive definite with smallest eigenvalue bounded from below by the smallest eigenvalue of ∇ uu L(y, u, λ).
The statement follows immediately from (11) and (15).
Comparing (14) with (18) shows that ∇ 2 J (u) ≈ ∇ 2 J (u) if y u (u) ≈ y u (u) and λ u (u) ≈ λ u (u). The latter two approximation conditions motivate our choice for the ROM matrix V. Proposition 2 Let V ∈ R n×r have rank r < n and satisfy that V T c y (y(u), u)V is invertible. If y(u + d) ∈ R(V) for all sufficiently small d, then
If, in addition, λ(u + d) ∈ R(V) for all sufficiently small d, then Applying the Implicit Function Theorem to this equation shows that
and ( 16) can be written as y u (u) = V y u (u). Since y(u + d) = V y(u + d) for all sufficiently small d,
which is (19). If, in addition, λ(u + d) ∈ R(V) for all sufficiently small d, then we can use similar argument to show the first equality in (20).
The second identity in (20) follows immediately from ( 14), ( 18), (19), and first identity (20).
In general there is no matrix V ∈ R n×r with small rank r n, such that the conditions in Proposition 2 hold. However, if there exists V ∈ R n×r , r n, with orthogonal columns such that the projections of y(u + d) and λ(u + d) onto R(V) are close to y(u + d) and λ(u + d), respectively, or equivalently such that
for a smal tolerance tol 1, then we expect y u (u) ≈ y u (u), λ u (u) ≈ λ u (u), and consequently that the ROM Hessian ∇ 2 J (u) is a good approximation of the original FOM Hessian ∇ 2 J (u).
However, finding V ∈ R n×r such that (21) is guaranteed to hold for all sufficiently small d would be computationally expensive. In practice we compute V ∈ R n×r using proper orthogonal decomposition (POD) applied to the state y(u) and the adjoint λ(u) at the current control.
ROM Hessians versus ROM Optimization Proposition 2 and even (21) describe optimistic scenarios. If we were able to find a V such that (21) holds at the current control u = u c , then it could be used to approximate the original problem (1, 2) in a trust-region based model management framework and approximately solve (4, 3) over all u in a neighborhood (trust-region) around the current control u c . Rather than using the ROM model u → J (V y(u), u) in a neighborhood of the current control u c , our approach uses the approximate quadratic Taylor expansion d →
The true function J (u) and this Taylor model have the same function and gradient value at u c , no matter how good V is. If ∇ 2 J (u c ) well approximates ∇ 2 J (u c ), we will essentially compute a Newton step. If the resulting V is less good, our approach can still generate a descent direction that is much better than a simple gradient step. Computationally inexpensive, but possibly less accurate ROMs can still be used to accelerate the overall optimization.
In this section we apply the ROM Hessian approach to a minimization problem that arises, e.g., from a discretized parabolic optimal control problem. Model Problem. Given functions : R n y → R, σ : R n u → R, F : R n y → R n y , and G : R n u → R n y , and given y 0 ∈ R n y consider the following minimization problem in the variables u = (u T 1 , . . . , u T n t ) T ∈ R n u n t , y = (y T 1 , . . . , y T n t ) T ∈ R n y n t .
Minimize
where y and u satisfy an implicit constraint,
It is easily possible to extent our ROM Hessian approach to generalizations of ( 22), e.g., replace (y k ) and G(u k ) by (y k , u k ) and G(y k , u k ). We consider (22) to simplify notation.
Throughout this section we make the following assumptions, which are adaptations of the assumptions (A1)-(A3) to the problem (22). (A1') There exists an open set D u ⊂ R n u n t such that for every u = (u T 1 , . . . , u T n t ) T ∈ D the state Eq. (22b) has a unique solution y(u) = (y 1 (u) T , . . . , y n t (u) T ) T . (A2') There exists an open set D y ⊂ R n y n t such that {y(u) : u ∈ D u } ⊂ D y , and the functions : D y → R, σ : D u → R, F : D y → R n y , and G : D u → R n y are twice continuously differentiable. (A3') For every y = (y T 1 , . . . , y T n t ) T ∈ {y(u) : u ∈ D u } and every r 1 , . . . , r n t ∈ R n y there exists a unique solution w 1 , . . . , w n t ∈ R n y of the equations
Under the assumptions (A1')-(A3') we can define the function J : 22), is a special case of (1, 2). Gradient and Hessian Computation. Under the assumptions (A1')-(A3') the function J : D u → R is twice continuously differentiable. Its gradient can be computed using the adjoint equation approach in Algorithm 1. Hessian-times-vector operations are computed using Algorithm 2.
T solve the state Eq. (22b). 2: Solve the adjoint equations for λ nt , . . . , λ 1 :
Note that since the objective function and the implicit constraints in the model problem ( 22) are separable, ∇ uy L(y, u, λ) = 0.
Hessian Approximation by Model Order Reduction. As mentioned towards the end of Sect. 2, we compute a ROM from state y and adjoint λ information. In this example the state y and the adjoint λ are vectors of vectors y k ∈ R n y and λ k ∈ R n y respectively. We compute a matrix V ∈ R n y ×r such that the components y k and λ k are approximately contained in the range of V (this V is a component of the projection matrix in Sect. 2) The projection matrix for y and λ is the n t n y × n t r block diagonal matrix with identical diagonal blocks given by V. Algorithm 2: Hessian-Times-Vector Computation -∇ 2 J (u)v 1: Given u = (u T 1 , . . . , u T nt ) T let y = (y T 1 , . . . , y T nt ) T solve the state Eq. ( 22b) and let λ = (λ T 1 , . . . , λ T nt ) T solve the adjoint Eq. ( 23). 2: Set w 0 = 0 and solve the linearized state equation w 1 , . . . , w nt :
3: Solve the second order adjoint equations for p nt , . . . , p 1 .
4: Compute the action of the Hessian
We apply POD to the computed state and adjoint to compute V ∈ R n y ×r . Since states and adjoint typically have very different scales, we do not apply POD to the combined snapshots, but individually, as stated in Algorithm 3. 2: Construct by truncated SVD matrices V y ∈ R ny ×ry and V λ ∈ R ny ×rλ so that
The ROM Hessian-time-vector computation is specified in Algorithm 4. The subspace matrix V ∈ R n y ×n r used in steps 2 and 3 is computed via Algorithm 3.
Computational Efficiency of ROM Hessian. While matrices like V T F y (y k )V in the ROM Hessian computation are smaller than F y (y k ), the dependence of F y on y k , which changes with time step k makes the computation of V T F y (y k )V expensive. This well-known issue [4,6,9,11,18] is addressed via hypereduction. Specifically, we use a so-called unassembled form of the Discrete Empirical Interpolation Method (DEIM), which originates from the (D)EIM [6,9] and is described in more detail Algorithm 4: ROM Hessian-Times-Vector Computation-∇ 2 J (u)v 1: Given u = (u T 1 , . . . , u T nt ) T let y = (y T 1 , . . . , y T nt ) T be the solution of the state Eq. (22b) and let λ = (λ T 1 , . . . , λ T nt ) T be the solution of the adjoint Eq. ( 23). 2: Set w 0 = 0 and solve the POD linearized state equation for w k ,
3: Solve the POD second order adjoint equation for p k ,
4: Compute the approximation to the action of the Hessian,
in [4,18]. This leads to a further approximation of the ROM Hessian computed by Algorithm 4. The error in this additional approximation is controlled by the tolerance tol DEIM ) used in DEIM. We refer to [14, Sects. 5.4, 6.2] for further details. We will use ∇ 2 J (u) to denote the final ROM Hessian approximation. We can replace y k ≈ V y k , where y k = V T y k ∈ R r , and λ k ≈ V λ k , where λ k = V T λ k ∈ R r to reduce storage requirements.
We apply the ROM Hessian approximation to semi-linear parabolic optimal control problems of the form min u∈L 2 ( ×(0,1))
where for given function u ∈ L 2 ( × (0, 1)) the function y(•, •; u) is the solution of and = (0, 1) 2 . The problem is discretized in space using piecewise linear finite elements on a triangulation obtained by dividing = (0, 1) 2 into 40 × 40 squares and then dividing each square into two triangles. This results in a semi-discretization with n u = 1681 and n y = 1521 degrees of freedom in the control and state respectively. The resulting semi-discretization is then discretized in time using the backward Euler method using n t = 100 times steps. This results in a problem of the type (22) with G(u k ) = M u u k and σ (u k ) = (α/2)u T k M u u k , where M u is the mass matrix, and
), where M is the mass matrix that also appears in (22b) (M u and M differ in size because of bounday conditions for y) and y k,d comes from a discretization of the desired state y d in (25a).
All optimization runs are initialized at u = 0. The truncated CG is stopped when negative curvature is detected or when the CG residual satisfies
depending on whether a FOM or ROM Hessian is used. Example 1: Cubic Reaction. In the first example the nonlinearity in (25b) is cubic, f (y) = y 3 . The desired state is y d (x, t) = 2e t + 2x 1 (x 1 -1) + 2x 2 (x 2 -1), the initial state is y 0 (x, t) = sin(2π x 1 ), and the control penalty is α = 10 -4 .
Both optimization with FOM Hessians and with ROM Hessians converged to virtually the same solution. The computed optimal control u at t = 0.1, 0.5, 1 is shown in Fig. 2. Table 1 shows that the optimization histories for the two approaches are nearly the same.
While the optimization histories for the two approaches are nearly the same, the ROM Hessian approach is much faster as shown in Table 2. The timing reported Fig. 1 Convergence history for Newton-CG with and without POD+DEIM approximations in the Hessian-vector computation in Table 2 is for computations using MATLAB running on a MacBook Pro (13inch, Retina, Mid 2014). The difference in computing time is due to the cost of Hessian-times-vector multiplications. In this example, using FOM Hessian requires 98 Hessian-times-vector multiplications and therefore 196 linear discretized PDE solves. In contrast, using ROM Hessian requires 102 Hessian-times-vector multiplications and therefore 204 linear ROM PDE solves. Figure 1 is another presentation of the convergence results in Tables 1 and2. The left plot shows that both approaches have nearly the same iteration history. However, when convergence history is plotted against computational work performed, measured in terms of PDE solves, the ROM Hessian approach is much faster. Note that here we do not differentiate between the nonlinear state PDE solve and the linear PDE solves needed for gradient, adjoint and FOM Hessian-times-vector computations. While the discretized state Eq. ( 22b) is nonlinear and computing y k+1 is performed via Newton's method, only 1-2 Newton iterations per time step in the state computation are needed in this example. Therefore the difference between a nonlinear state PDE solve and linear PDE solve is small. Example 2: Solid Fuel Ignition Model. This example is modeled after [8]. The nonlinearity in (25b) is f (y) = -δe y with δ = 5. The desired state is y d (x, t) = π -2 sin(π x 1 ) sin(π x 2 ), the initial state is y 0 ≡ 0, and the penalty is α = 5 • 10 -3 .
The numerical results for this example mirror those of the previous example. Optimization with FOM Hessians and with ROM Hessians converged to virtually the same solution, and the optimization histories for the two approaches are nearly the same, as shown in Table 3. The computed optimal control u at t = 0.1, 0.5, 1 is shown in Fig. 4.
The ROM Hessian approach is much faster as shown in Table 4. Again, the timing reported in Table 4 is for computations using MATLAB running on a MacBook Pro (13-inch, Retina, Mid 2014) and the reason for the difference in computing time is the cost of Hessian-times-vector multiplications. The ROM Hessian requires 96 ROM PDE solves for ROM Hessian-vector operations in contrast to the 96 FOM PDE solves in the FOM Hessian approach. Figure 3 shows that when convergence history is plotted not against iteration count, but against computational work performed, measured in terms of PDE solves, the ROM Hessian approach is much faster.
We have introduced ROM Hessian approximations for use in inexact Newton methods for large-scale smooth optimization problems obtained from discretizations of optimal control problems governed by parabolic PDEs. In contrast to other ROM approaches, which approximate the optimization problem, our approach retains the original FOM objective function and gradient. The Hessian ROM is computed by applying POD to state and adjoint snapshots that have to be computed anyway, and therefore ROM computation is relatively inexpensive. Since original objective functions and gradients are used, the qualitative global convergence properties of the original line-search Newton method and the line-search Newton method with ROM Hessian approximation are the same. However, since Hessian approximations are significantly cheaper, the ROM Hessian approach can lead to substantial savings. This is confirmed by numerical experiments with two semilinear parabolic optimal control problems. The two optimization approaches had essentially the same convergence behavior, but the ROM Hessian approach had a factor 5-8 speedup because of computational savings due to the Hessian approximations. However, it is possible to construct examples, where the ROM Hessian approach does not lead to computational savings. When the ROM Hessian too poorly approximates the true one, the ROM Hessian approach can require substantially more optimization iterations, which erase savings enjoyed within an optimization iteration. Additional analysis and numerical tests are part of future work.