# A trace-constrained maximization problem in the cone of positive definite matrices

Let $A\in\mathbb{R}^{n\times n}$ be a matrix having eigenvalues with strictly negative real part (in other words, $A$ is supposed to be Hurwitz stable). Let $\mathrm{tr}(\cdot)$ denote the trace operator and $I$ the identity matrix. Let $>$ be the standard partial order in the cone of positive definite matrices.

Consider the following optimization problem $$\tag{\star}\label{prob} m(A) :=\sup_{\substack{X\in\mathbb{R}^{n\times n}\\ X>0,\ \mathrm{tr}(X)=1}} \mathrm{tr}\left(X(I+P)^{-1}\right),$$ where $P$ is the (unique) positive definite solution of the following Lyapunov equation $AP+PA^\top=-X$.

My question. Suppose that $A$ is upper triangular, does the following inequality $$m(A)\ge \frac{2\,\mathrm{tr}\,A}{2\,\mathrm{tr}\,A-1}$$ always hold true?

A special case. Notice that if $A+A^\top<0$, then the answer is in the affirmative. To see this, just pick $P^\star= -\frac{1}{2\mathrm{tr}(A)}I$ and observe that $X^\star=-(AP^\star+P^\star A^\top)$ is positive definite and satisfies $\mathrm{tr}(X^\star)=1$.

Numerical evidences. Quite surprisingly, after runnning an extensive number of numerical simulations, it seems that the answer is in the affirmative for any (upper triangular Hurwitz stable) $A$. More precisely it seems that (modulo numerical errors of magnitude $\sim 10^{-6}$) the conjectured inequality actually holds with equality, that is $$m(A) = \frac{2\,\mathrm{tr}\,A}{2\,\mathrm{tr}\,A-1}.$$

However, this fact does not seem trivial to prove (I spent quite some time thinking about this, but I didn't manage to prove it); so any help in clarifying this conjecture is greatly appreciated. Thanks!

An (perhaps useful?) equivalent formulation. By plugging $X=-AP-PA^\top$ into the trace functional in \eqref{prob}, the latter can be rewritten as \begin{align}\tag{$\star\star$} m(A) &=\sup_{\substack{X\in\mathbb{R}^{n\times n}\\ X>0,\ \mathrm{tr}(X)=1}} -\mathrm{tr}\left((AP+PA^\top)(I+P)^{-1}\right)\notag \\ &=\sup_{\substack{X\in\mathbb{R}^{n\times n}\\ X>0,\ \mathrm{tr}(X)=1}} -2\mathrm{tr}\left(A(I+P^{-1})^{-1}\right)\\ &=-2\inf_{\substack{X\in\mathbb{R}^{n\times n}\\ X>0,\ \mathrm{tr}(X)=1}} \mathrm{tr}\left(A(I+P^{-1})^{-1}\right)\notag\\ &=-2\inf_{\substack{P\in\mathbb{R}^{n\times n}\\ P>0,\ \mathrm{tr}(AP)=-\frac{1}{2}\\ AP+PA^\top<0}} \mathrm{tr}\left(A(I+P^{-1})^{-1}\right)\\ &\overset{(\#)}{=}-2\inf_{\substack{P\in\mathbb{R}^{n\times n}\\ P>0,\ \mathrm{tr}(AP)=-\frac{1}{2}\\ AP+PA^\top<0}} \mathrm{tr}\left(A-A(I+P)^{-1}\right) \\ &=-2\,\mathrm{tr}\,A + 2\sup_{\substack{P\in\mathbb{R}^{n\times n}\\ P>0,\ \mathrm{tr}(AP)=-\frac{1}{2}\\ AP+PA^\top<0}} \mathrm{tr}\left(A(I+P)^{-1}\right), \\\label{prob-eq} \end{align} where in (#) I used the Woodbury matrix identity.

• @Mahdi: I've fixed this, thanks! – Ludwig Sep 4 at 14:33
• Notice that in the last line of your equivalent formulation, the infimum can be taken over $P > 0, \text{tr}(AP) = -\frac{1}{2}$. This is because whenever $X>0$ and $A$ is Hurwitz, it is guaranteed that $P>0$. Not sure if this is easier to deal with ... – Abhishek Halder Sep 13 at 18:50
• @AbhishekHalder: Nice observation, thanks! I've included it in the equivalent formulation of the problem. – Ludwig Sep 14 at 0:34
• @AbhishekHalder: On second thought, I think that the infimum can be taken over $P>0$, $\mathrm{tr}(AP)=-1/2$, only if we add the further constraint $AP+PA^\top<0$. – Ludwig Sep 15 at 22:14

In fact, the solution seems to be invariant under similarity transformations of $A$, i.e. $A\mapsto M A M^{T}$, where $M\in O(n)$.

Edit: I think I need to add the assumption of $X$ being symmetric to apply the law of inertia. But since $X$ appears in the Lyapunov equation, it is standard to assume that $X$ is symmetric.

Indeed, let $m(A)$ be the maximum as defined by $$m(A)=\max_{\substack{X\in\text{Sym}(n)\\ X>0,\ \mathrm{tr}(X)=1}} \mathrm{tr}\left(X(I+P(X,A))^{-1}\right),$$

where $P(X,A)$ is the solution to the Lyapunov equation with parameters $X,A$. It can be shown that $P(MXM^T,MAM^{T})=MP(X,A)M^T$ for all orthogonal $M$. Using invariance of the trace under cyclical permutations, as well as the fact that the constraint space is preserved under $X\mapsto MXM^T$ (Sylvester's law of inertia plus invariance of the trace), we obtain $$m(MAM^{T})=m(A).$$

• Yes, $X$ is supposed to be symmetric and positive definite. I'm wondering whether your comment can be used to show that $m(A)$ depends only on the eigenvalues of $A$. – Ludwig Sep 4 at 16:02
• This follows from the last equation. The eigenvalues of $A$ and $MAM^T$ are the same. – S.Surace Sep 4 at 18:21
• Or maybe I misunderstood your comment. Are you saying that you found two matrices $A,A'$ that are not related by an orthogonal transformation, but have the same eigenvalues, and $m(A)=m(A')$? Then I do not know how to show that this holds in general. – S.Surace Sep 4 at 19:01
• Yes, this is what I observed in my simulations and I would like to prove (or disprove). – Ludwig Sep 4 at 19:13
• If $m$ indeed depends only on the spectrum of $A$ this should reflect a symmetry of your stochastic system. But you are probably right that we don't get much more insight from that; the claim should be simply provable from the explicit form of $m$. – S.Surace Sep 7 at 19:46

Not a full answer, but some ideas. From the last line of your equivalent formulation, we need to solve

$$J^{*} = \sup_{\substack{P\in\mathbb{R}^{n\times n}\\ P>0,\ \mathrm{tr}(AP)=-\frac{1}{2}}} \mathrm{tr}\left(A(I+P)^{-1}\right),$$

which, after letting $Q:=I + P$, and $a:=\text{tr}(A)$, becomes

$$J^{*} = \sup_{\substack{Q\in\mathbb{R}^{n\times n}\\ Q>I,\ \mathrm{tr}(AQ)=a-\frac{1}{2}}} \mathrm{tr}\left(AQ^{-1}\right).$$

If $A$ were negative diagonal then surely $\mathrm{tr}\left(AQ^{-1}\right)$ would have been concave, making the above a convex optimization problem. (Question: does concavity still hold for Hurwitz $A$?)

In any case, the first order optimality conditions give: $\lambda Q A Q = A$, where $\lambda >0$ is the Lagrange multiplier, and $\mathrm{tr}(AQ)=a-\frac{1}{2}$. Taking trace of $\lambda Q A Q = A$ gives $J^{*}=\lambda(a-1/2)$, where $\lambda>0$ remains to be determined (not sure yet how).

• Thanks for sharing your ideas! Shouldn't be $Q> I$ instead of $Q> 0$? – Ludwig Sep 15 at 18:45
• You are right, corrected. – Abhishek Halder Sep 15 at 21:49
• Actually, solving $Q>I$ from $\lambda Q A Q = A$ is same as solving for $P>0$ from the algebraic Riccati equation: $PAP + PA + AP + (1 - 1/\lambda)A = 0$. – Abhishek Halder Sep 15 at 22:07