Éphane St

Zheng Gaubert
email: zhengqu@hku.hk

AND Srinivas Qu

Sridharan

Stéphane Gaubert
email: stephane.gaubert@inria.fr

Srinivas Sridharan

MAXIMIZING CONCAVE PIECEWISE AFFINE FUNCTIONS ON THE UNITARY GROUP

We show that a convex relaxation, introduced by Sridharan, McEneaney, Gu and James to approximate the value function of an optimal control problem arising from quantum gate synthesis, is exact. This relaxation applies to the maximization of

(2)

Here, U (n) denotes the unitary group of degree n, i.e., the subset of unitary matrices of the space M n (C) of n × n complex matrices, B(n) denotes the nth unitary ball, i.e., the set of matrices of M n (C) with singular values at most 1; •, • denotes the canonical inner product of M n (C), thought of as a real Hilbert space, so that X 1 , X 2 = Re(trace(X * 1 X 2 )), where X * denotes the conjugate transpose of a matrix X, and finally P 0 , . . . , P m ∈ U (n), c 0 , . . . , c 1 ∈ R , are given. Since B(n) is the convex hull of U (n), Problem 2 is a convex optimization problem. Moreover, since the convex set B(n) can be represented in terms of elementary matrix inequalities, Problem 2, and so, Problem 1, can be solved efficiently by convex programming techniques. We refer the reader to [START_REF] Ben-Tal | Lectures on modern convex optimization[END_REF] for more background on complexity results in convex programming, including general polynomial time complexity bounds which apply to the formulation (2).

1.2. Motivation from Optimal Control. The present contribution arose from the pruning step of the curse of dimensionality-free numerical method, introduced by McEneaney [START_REF] Mceneaney | A curse-of-dimensionality-free numerical method for solution of certain HJB PDEs[END_REF], for solving first order Hamilton-Jacobi-Bellman Partial Differential Equations (HJB PDE). This method applies to a special class of optimal control problems in which the Hamiltonian is given or approximated by a pointwise supremum of a finite number of elementary Hamiltonians. Such Hamiltonians arise naturally when the control space comprises a discrete component, for example in switched systems. The value function V of the optimal control problem is approximated by the infimum (or supremum if the objective is maximized) of a finite number of elementary basis functions {φ 0 , . . . , φ m }, which are typically linear or quadratic forms:

V min 0 i m φ i . (3)
Theoretical estimates show that, to reach an approximation of a prescribed accuracy, the computational complexity grows polynomially with the space dimension [START_REF] Mceneaney | Convergence rate for a curse-of-dimensionalityfree method for a class of HJB PDEs[END_REF][START_REF] Qu | Contraction of Riccati flows applied to the convergence analysis of a max-plus curse-of-dimensionality-free method[END_REF]. However, the number of basis functions is multiplied by the number of elementary Hamiltonians at each propagation step, the so called curse of complexity. Therefore in order to attenuate this curse of complexity in practice, a pruning operation is introduced. It achieves this improvement in efficiency by removing, at each step, a certain number of basis functions less useful than the others.

The pruning problem has received much attention [MDG08, McE09, SGJM10, GMQ11, SMGJ14, GQS14]. Actually, the experiments in [START_REF] Gaubert | Curse of dimensionality reduction in maxplus based approximation methods: Theoretical estimates and improved pruning algorithms[END_REF] indicate that the main part of the execution time of the curse of dimensionality-free method is spent in pruning. Therefore, pruning is both a critical step and a bottleneck stagenecessitating the development of fast exact or approximate algorithms. The pruning problem appears as well in other approximate dynamic programming methods in which the value function is approximated by a maximum of affine functions, like the dual dynamic programming method of Pereira and Pinto [START_REF] Pereira | Multi-stage stochastic optimization applied to energy planning[END_REF] and its extensions, note in particular the discussion by Shapiro [START_REF] Shapiro | Analysis of stochastic dual dynamic programming method[END_REF].

Pruning was shown in [START_REF] Gaubert | Curse of dimensionality reduction in maxplus based approximation methods: Theoretical estimates and improved pruning algorithms[END_REF] to be a continuous version of the facility location problem. Most pruning methods used in practice rely on calculating the importance metric δ j of each basis function φ j , with respect to a bounded state space X . This metric is defined by

δ j := max x∈X min i =j φ i (x) -φ j (x) . (4)
The importance metric δ j thus measures the maximal loss caused by removing the basis function φ j . When the state space X is not bounded, δ j is replaced by a variant, involving a normalization, see [START_REF] Mceneaney | Curse-of-complexity attenuation in the curse-of-dimensionality-free method for HJB PDEs[END_REF]. In some sense the higher δ j is, the more useful the function φ j is for the approximation. In particular, if δ j < 0, then, the basis function j can be deleted without changing the approximation of V in (3). Once we get all the values {δ 0 , . . . , δ m }, the simplest algorithm consists of removing those basis functions with lowest importance metrics, as suggested in [START_REF] Mceneaney | Curse-of-complexity attenuation in the curse-of-dimensionality-free method for HJB PDEs[END_REF]. One may also apply an iterative greedy algorithm, in which the importance metric is recomputed dynamically after each removal of basis function, or use the importance metric to construct a discretized problem to which combinatorial facility location algorithms can be applied [START_REF] Gaubert | Curse of dimensionality reduction in maxplus based approximation methods: Theoretical estimates and improved pruning algorithms[END_REF].

As mentioned before, the curse of dimensionality-free method can be used to solve switched optimal control problems. In particular, it has been recently applied to the optimal quantum circuit synthesis, following a work by Gu, James, McEneaney and Sridharan, see [START_REF] Sridharan | Reduced-complexity numerical method for optimal gate synthesis[END_REF][START_REF] Sridharan | A reduced complexity minplus solution method to the optimal control of closed quantum systems[END_REF] and also [START_REF] Gaubert | Bundle-based pruning in the max-plus curse of dimensionality free method[END_REF]. The related optimal control problem is to find a least path-length trajectory on the special unitary group SU (n). For this particular application, the basis functions are chosen to be affine functions defined on the space M n (C):

φ i (X) = P i , X + c i , ∀X ∈ M n (C), i = 0, . . . , m,
with P 0 , . . . , P m ∈ U (n) and c 0 , . . . , c m ∈ R. Then calculating the importance metric of φ 0 reduces to solving the following optimization problem:

δ 0 := max X∈SU (n) min 1 i m P i -P 0 , X + c i -c 0 . (5)
The importance metric of other functions {δ 1 , . . . , δ m } can be obtained by solving similar programs (in which the index 0 is replaced by an index j). A key difficulty here is that the state space (the special unitary group SU (n)) is non-convex. Hence, in order to compute efficiently an importance metric, we need to perform a relaxation. In [START_REF] Sridharan | A reduced complexity minplus solution method to the optimal control of closed quantum systems[END_REF], the special unitary group SU (n) is first relaxed into the unitary group U (n), then replaced by the unitary ball B(n). Then, the value δ 0 is approximated by the optimal value δ0 of the convex optimization problem (2). Since δ0 δ 0 , the condition δ0 < 0 guarantees that the function φ 0 is not active in the infimum in (3) and can therefore be pruned. More generally, any relaxation leading to an importance metric δ 0 δ 0 yields a sound pruning method. The closer the relaxed importance metric is from the true one, the more efficient is the pruning. 1.3. Discussion of the results. One question which arises from [START_REF] Sridharan | A reduced complexity minplus solution method to the optimal control of closed quantum systems[END_REF] is whether the relaxation of (1) into (2) is exact. Our main result, Theorem 2.1, shows that the optimal solution of (2) contains always a unitary matrix. Thus, Program (1) is equivalent to its convex relaxation (2). This property is somehow surprising, as the maximum of a concave function over a compact convex set is generally not attained at an extreme point of a set. What makes this property valid for the special non-convex optimization problem (1) is that the gradients of the affine constituents of the objective function are differences P i -P 0 , 1 i m, where each P i is unitary, and P 0 is a fixed unitary matrix.

The proof of Theorem 2.1 exploits the strict convexity of the Frobenius norm. The proof also uses a regularization argument (replacing the piecewise affine objective function by a smooth log-exp type function), together with algebraic properties of the tangent cone at a given point of B(n). The analogous result holds, with obvious changes, when SU (n) is replaced by the nth real special orthogonal group SO(n). The latter also arises in optimal estimation and filtering, in particular in optimal attitude estimation [START_REF] Sridharan | Deterministic filtering for optimal attitude estimation on SO(3) using max-plus methods[END_REF], and so, the present results apply as well to this case.

To get more insight on the assumptions which govern Theorem 2.1, we finally present a variation of this theorem, showing that a similar result of independent interest holds (Proposition 2.2) for the maximization over a sphere of a concave piecewise affine functions the affine constituents of which have unitary gradients.

We finally note that Theorem 2.1 has been announced, without proof, in our recent conference article [START_REF] Gaubert | Bundle-based pruning in the max-plus curse of dimensionality free method[END_REF]. We showed there that this result can be combined with scalable bundle-type non-differentiable optimization methods to speed up the pruning procedure of [START_REF] Sridharan | A reduced complexity minplus solution method to the optimal control of closed quantum systems[END_REF].

1.4. Summary of notation. It is convenient to list here the main notation used throughout the paper. As mentioned above, M n (C) denotes the space of n × n complex matrices. For X ∈ M n (C), X * denotes its conjugate transpose. The space M n (C) is considered as a real Hilbert space endowed with the inner product •, • given by X 1 , X 2 = Re(trace(X * 1 X 2 )), ∀X 1 , X 2 ∈ M n (C). We denote by I n the n × n identity matrix, so that the nth unitary group is given by U

(n) = {X ∈ M n (C) | XX * = I n }. The space of n × n positive semidefinite (resp. positive definite) matrices is denoted by S + n (resp. Ŝ+ n ). For two Hermitian matrices A, B ∈ M n (C), we write A B (resp. A B) if A -B ∈ S + n (resp. A -B ∈ Ŝ+ n ). Hence, the unitary ball B(n) is given by B(n) = {X ∈ M n (C) | XX * I n }.

Maximizing a class of concave functions over the unitary group

In the following we show that (2) is an exact relaxation of (1). Our main result is Theorem 2.1. We begin with a few lemmas.

Lemma 2.1. Let P 0 , . . . , P m be n × n unitary matrices and let X ∈ M n (C). If XP 0 is a point in the relative interior of the convex hull of {XP i } i=1,...,m , then XP i = XP 0 for all i = 1, . . . , m.

Proof. If XP 0 belongs to the relative interior of the convex hull of {XP i } i=1,...,m , then there are α 1 , . . . , α m > 0 such that

m i=1 α i = 1, m i=1 α i XP i = XP 0 .
The Frobenius norm of XP i equals to that of X for all i = 0, . . . , m. By the strict convexity of the Frobenius norm, we deduce that XP i = XP 0 holds for all i = 1, . . . , m.

The tangent cone to B(n) at X ∈ B(n) is defined by

T B(n) (X) = cl{λ(Z -X) : λ 0, Z ∈ B(n)} ,
where cl denotes the closure of a set, see [START_REF] Rockafellar | Variational analysis[END_REF]p. 204].

Lemma 2.2. Let 1 k n and Σ be a diagonal matrix with non-negative real diagonal entries (λ 1 , • • • , λ n ) such that λ j = 1 for all j = 1, . . . , k -1 and λ j < 1 for all j = k,

• • • , n. Then X 1 X 2 X 3 X 4 ∈ M n (C) : -X 1 -X * 1 ∈ Ŝ+ k-1 ⊂ T B(n) (Σ) .
Proof. If k = 1 then it is clear that Σ is an interior point of B(n) and the tangent cone at Σ is the whole space M n (C). We next consider the case when 1 < k n.

For ease of proof, we write Σ as a block matrix Σ = I k-1 0 0 Σ4 ,where Σ 4 is a diagonal matrix such that

I n-k+1 -Σ 4 Σ * 4 0. Let X = X1 X2 X3 X4 ∈ M n (C), such that -X 1 -X * 1 ∈ Ŝ+ k-1 . For the simplicity of notation, let H 1 ∈ M k-1 (C), H 4 ∈ M n-k+1 (C) and H 2 ∈ C (k-1)×(n-k+1) such that XX * = H1 H2 H * 2 H4 .Take any scalar ∆ > 0. Then, (Σ + ∆X)(Σ + ∆X) * -I n = ΣΣ * + ∆(XΣ * + ΣX * ) -I n + ∆ 2 XX * = ∆(X 1 + X * 1 ) ∆(X 2 Σ * 4 + X * 3 ) ∆(Σ 4 X * 2 + X 3 ) Σ 4 Σ * 4 -I n-k+1 + ∆(Σ 4 X * 4 + X 4 Σ * 4 ) + ∆ 2 H 1 ∆ 2 H 2 ∆ 2 H * 2 ∆ 2 H 4
By Schur's complement lemma, we know that (Σ + ∆X)(Σ + ∆X) * ≺ I n if and only if X 1 + X * 1 + ∆H 1 ≺ 0 and

I n-k+1 -Σ 4 Σ * 4 -∆(Σ 4 X * 4 + X 4 Σ * 4 ) -∆ 2 H 4 + ∆(Σ 4 X * 2 + X 3 + ∆H * 2 )(X 1 + X * 1 + ∆H 1 ) -1 (X 2 Σ * 4 + X * 3 + ∆H 2 ) 0 . Since X 1 + X * 1 ≺ 0 and I n-k+1 -Σ 4 Σ * 4
0, there is ∆ > 0 such that the latter two inequalities hold thus (Σ + ∆X)(Σ + ∆X) * ≺ I n . Hence, X is in the tangent cone of B(n).

In the sequel, let P 0 , . . . , P m be n × n unitary matrices, and c 0 , . . . , c m be real numbers. For all β > 0, define φ β : M n (C) → R by:

φ β (X) = -β -1 log m i=1 e -β( Pi-P0,X +ci-c0) , ∀X ∈ M n (C) . (6) Lemma 2.3. There is a constant K > 0 independent of β > 0 such that the map φ β : M n (C) → R is K-Lipschitz continuous.
Proof. Let β > 0. For all X, Y ∈ M n (C), the differential map of φ β at point X, Dφ β (X) evaluated at Y , satisfies:

Dφ β (X) • Y = m i=1 α i P i -P 0 , Y ,
where

α i = ( m j=1
e -β( Pj -P0,X +cj -c0) ) -1 e -β( Pi-P0,X +ci-c0) , i = 1, . . . , m.

Thus α 1 , . . . , α m > 0 and m i=1 α i = 1. Hence, for all X, Y ∈ M n (C), by the Cauchy-Schwarz inequality,

|Dφ β (X) • Y | m i=1 α i P i -P 0 Y max 1 i m P i -P 0 Y ,
where • denotes the Euclidean norm associated to the Frobenius scalar product •, • . It follows that the function φ β is K-Lipschitz with respect to the norm • for all β > 0, with K = max 1 i m P i -P 0 .

Proposition 2.1. For every β > 0, the set of optimal solutions of the optimization problem

max X∈B(n) φ β (X) (7)
contains a unitary matrix.

Proof. Let U 0 ∈ B(n) be an optimal solution of (7). Suppose that U 0 is not unitary. Consider the SVD decomposition of U 0 given by U 0 = V 1 ΣV 2 , where Σ is a diagonal matrix with non-negative real diagonal entries (λ 1 , . . . , λ n ), listed in non-increasing order. Let k ∈ {1, . . . , n} such that λ i = 1 for all i = 1, . . . , k -1 and λ j < 1 for all j = k, . . . , n. Then Σ is an optimal solution of the following optimization problem: max

X∈B(n) φ β (V 1 XV 2 ) . (8)
The first-order optimality condition [RW98, p.207] implies that

Dφ β (V 1 ΣV 2 ) • (V 1 Y V 2 ) 0, ∀Y ∈ T B(n) (Σ) . (9)
We have:

Dφ β (V 1 ΣV 2 ) • (V 1 Y V 2 ) = m i=1 α i P i -P 0 , V 1 Y V 2 ,
where α 1 , . . . , α m > 0 and

m i=1 α i = 1. Therefore, Dφ β (V 1 ΣV 2 ) • (V 1 Y V 2 ) = m i=1 α i (V * 1 P i V * 2 -V * 1 P 0 V * 2 ), Y
By the first-order optimality condition (9) and Lemma 2.2, we deduce that

m i=1 α i (V * 1 P i V * 2 -V * 1 P 0 V * 2 ), X 0 , for all X = X1 X2 X3 X4 ∈ M n (C) such that -X 1 -X * 1 ∈ Ŝ+ k-1 . Hence, m i=1 α i (V * 1 P i V * 2 -V * 1 P 0 V * 2 ) = Z 0 0 0 for some Z ∈ M k-1 (C). Therefore, (I n -Σ) m i=1 α i (V * 1 P i V * 2 -V * 1 P 0 V * 2 ) = 0 .
By Lemma 2.1, we know that

(I n -Σ)V * 1 P i V * 2 = (I n -Σ)V * 1 P 0 V * 2 , i = 1, . . . , m. This implies that P i -P 0 , V 1 V 2 = P i -P 0 , V 1 ΣV 2 , i = 1, . . . , m. Therefore, φ β (V 1 V 2 ) = φ β (V 1 ΣV 2 ) = max X∈B(n) φ β (X).
Theorem 2.1. Let P 0 , . . . , P m be n × n unitary matrices, and let c 0 , . . . , c m be real numbers. The set of optimal solutions of the following optimization problem:

max X∈B(n) min 1 i m P i -P 0 , X + c i -c 0 (10)
contains a unitary matrix.

Proof. Denote φ(X) = min 1 i m P i -P 0 , X + c i -c 0 , for X ∈ M n (C). By Lemma 2.3 and the Arzelà-Ascoli theorem, the function φ β defined in (6), which converges pointwise to φ as β goes to +∞, converges uniformly to the same function. For each β, by Proposition 2.1 the intersection

O β := U (n) ∩ arg max X∈B(n)
φ β is not empty. Since the convergence of φ β to φ is uniform, each cluster point of a sequence {U β } β 0 with U β ∈ O β for all β > 0 is an optimal solution of the problem (10), see [START_REF] Rockafellar | Variational analysis[END_REF]p.266]. The cluster point is unitary because U (n) is closed. Thus the optimization problem (10) must have a unitary optimal solution. By Theorem 2.1, solving (2) is equivalent to solving (1). Theorem 2.1, which concern matrices with complex entries, has a real analogue, in which SU (n) is replaced by the real special orthogonal group SO(n). The latter group appeared for instance in the application of max-plus methods to a deterministic filtering problem [START_REF] Sridharan | Deterministic filtering for optimal attitude estimation on SO(3) using max-plus methods[END_REF]. Then, U (n) is replaced by the nth real orthogonal group O(n), and B(n) is replaced by B R (n), the set of n × n matrices X with real entries such that XX I n . We next state the real analogue of Theorem 2.1.

Theorem 2.2. Let P 0 , . . . , P m be n × n orthogonal matrices, and let c 0 , . . . , c m be real numbers. Then, the set of optimal solutions of the following optimization problem:

max X∈B R (n) min 1 i m P i -P 0 , X + c i -c 0 (11)
contains an orthogonal matrix.

Proof. The proof is identical to the one of Theorem 2.1, up to trivial changes.

Example 2.1. The elementary situation in which n = 1, so that

O(n) = {1, -1}, B R (n) = [-1, 1], U (
n) is the unit circle, and B(n) is the unit disk, will allow us to see that our restrictive assumptions in Theorem 2.1 and 2.2 are useful. The latter result states that as soon as P 0 , . . . , P m ∈ {-1, 1}, the maximum of the function φ on the interval [-1, 1] is always attained at the boundary of this interval. Indeed, if P 0 = -1, the gradients of the affine constituents of φ belong to {0, 2}, and the maximum of φ is attained at point 1, whereas if P 0 = 1, the same gradients belong to {-2, 0}, and the maximum of φ is now attained at point -1. A consideration of the n = 1 case shows that Theorem 2.2 does not carry over to piecewise linear concave functions of the form φ(X) = min

1 i m P i -Q i , X + d i , (12)
where P i , Q i are orthogonal, and d i ∈ R. Indeed, since P i , Q i ∈ {±1}, P i -Q i can achieve any value in {-2, 0, 2}, and so, φ(X) = min(2X, -2X) is of the form (12). The maximum of the latter function over B R = [-1, 1] is equal to 0, whereas its maximum over O(1) = {1, -1} is equal to -1. A similar counter example holds in the complex case. Consider P

1 = 1, Q 1 = -1, P 2 = -1, Q 2 = 1, P 3 = i, Q 3 = -i, P 4 = -i, Q 4 = i, we get φ(X) = min(-2| Re X|, -2| Im X|) .
The maximum of this function over B(1) is equal to 0, whereas it is readily seen that its maximum over U (1) is equal to -√ 2 (attained by the unitary vector X

= (1 + i)/ √ 2).
We next state an elementary variation of Theorem 2.1, in which a concave piecewise affine function is maximized over a sphere.

Let • denote a norm on R n , and let • * denote the dual norm, so that

p * = max y 1 p • y = max y =1 p • y, ∀p ∈ R n , (13) and
x = max

q * 1 q • x = max q * =1 q • x, ∀x ∈ R n ,
where • denotes the standard scalar product of R n , and the notation "max" means that each supremum is attained. Let S, B denote the primal unit sphere and the primal ball, so that S := {x ∈ R n | x = 1} and B := {x ∈ R n | x 1}. Similarly, let S * := {p ∈ R n | p * = 1} and B * := {p ∈ R n | p * 1} denote the dual unit sphere and ball. We refer the reader to [AB99, Ch. 6] for more background on convex duality and norms.

Proposition 2.2. Let p 0 , . . . , p m ∈ S * , c 0 , . . . , c m ∈ R, and let:

φ(x) = inf 1 i m p i -p 0 , x + c i -c 0 .
Then, max Since p 0 * = 1, by (13), we can find a vector y ∈ R n such that y = 1 and p 0 • y = 1. Moreover, p i • y p i * y = 1 holds for all i. We deduce from (14) that 1 = p 0 • y = 1 i m α i p i • y, and so, 1 = p 0 • y = p i • y, for 1 i m. It follows that for all t 0, φ β (x + ty) = φ β (x). Since x + ty = x < 1 for t = 0, and x + ty → ∞ when t → ∞, using the intermediate value theorem, we deduce that there is a parameter t such that x + ty ∈ S minimizes φ β over B. Then, we conclude as in the proof of Theorem 2.1, letting β tend to ∞.

Remark 2.2. An inspection of the n = 1 case shows that the optimization problem (1) is not equivalent to the optimization problem (5), so that the present results only yield an exact convex relaxation for the former problem. An exceptional case in which there is an exact convex relaxation of a problem of type (5) is when the state space of the control problem is the group SO(2). This yields a program of type (5) where the matrices {P 0 , . . . , P m } are all in SO(2). Since SO(2) can be identified to the unit circle of R 2 , Proposition 2.2 yields an exact convex relaxation. We are not aware of exact convex relaxations of (5) in general situations.

  Main result. The main object of this paper is to show the equivalence between the following non-convex optimization problem maxX∈U (n) min 1 i m P i -P 0 , X + c i -c 0 P i -P 0 , X + c i -c 0 .

  Let φ β denote the log-exp regularization of the function φ, defined as in (6). Then, arguing as in the proof of Proposition 2.1, we see that any maximum point x of φ β on B satisfies1 i m α i p i -p 0 , z 0, ∀z ∈ T B (x) ,where α 1 , . . . , α m are positive numbers with sum 1. If x < 1, we have T B (x) = R n , and so,1 i m α i (p i -p 0 )

The second author carried out parts of this work when she was with INRIA and CMAP, École Polytechnique, CNRS and subsequently with the School of Mathematics of University of Edinburgh. The third author carried out part of this work when he was with UMA, ENSTA, Palaiseau, France.The first two authors were partially supported by the PGMO Program of FMJH and EDF, and by the program "Ingénierie Numérique & Sécurité" of the French National Agency of Research, project "MALTHY", number ANR-13-INSE-0003.