5

For a given $A=\begin{bmatrix}a_1&&&\\&a_2&1&\\&&a_2&\\&&&a_2\end{bmatrix}$ and $B=\begin{bmatrix}b_{11}&b_{12}\\b_{21}&b_{22}\\b_{31}&b_{32}\\0&b_{42}\end{bmatrix}$, where $a_i>1$, $a_1\neq a_2$, $a_i$ and $b_{ij}$ are real,

\begin{array}{ll} \underset{X \in \mathbb{R}^{4\times 4}}{\text{minimize}} & \mathrm{tr} \left( B^T X B \right)\\ \text{subject to} & X=A^TX(I+BB^TX)^{-1}A.\end{array}

Here $X$ is the unique solution to DARE (discrete-time algebraic Riccati equation). For a specific $A$ and $B$, I can use Matlab to solve the DARE and then insert the resulting $X$ to find $\mathrm{tr}(B^TXB)$. Is it possible to get the general answer in terms of $A$ and $B$?

I have asked this question for general $A$ and $B$ here, but now trying to solve it for simpler case.

Currently I have the following bounds:

$$a_1^2a_2^4+a_2^2-2\geq \mathrm{tr}(B^TXB) \geq 2a_1a_2^3-2.$$

Matlab code to calculate the $\mathrm{tr}(B^TXB)$ when $A$ is fixed for various $B'$s:

A=[5 0 0 0; 0 2 1 0; 0 0 2 0; 0 0 0 2];
Q=zeros(4,4);
for i=1:100
    B=2.*rand(4,2)-1*ones(4,2);
    [X,K,L] = dare(A,B,Q);
    t(i)=trace(B'*X*B);
end
Lee
  • 1,978
  • 3
    One thing that might help simplify things: $\operatorname{tr}(B^T X B) = \operatorname{tr}(BB^TX)$, so $B$ always appears as $BB^T$. – eyeballfrog May 06 '22 at 14:29
  • Since $A$ is invertible, the constraint can also be rewritten as $XA^{-1}(I + BB^TX) = A^TX$, which makes it into a quadratic in $X$. Though, with 16 equations in 16 unknowns, doesn't this have a finite number of solutions $X$? – eyeballfrog May 06 '22 at 14:59
  • @eyeballfrog DARE has unique solution $X$ for any given $A$ and $B$. However, I want to find $\mathrm{tr}(B^TXB)$ in terms of $A$ and $B$, so that I don't have to solve DARE at all. – Lee May 06 '22 at 15:30
  • Can $X$ be assumed to be symmetric? – eyeballfrog May 06 '22 at 15:38
  • 2
    @eyeballfrog yes, it is symmetric and positive definite – Lee May 06 '22 at 16:03
  • @Lee Where is that DARE coming from? It does not seem to have the usual structure. – KBS May 08 '22 at 22:54
  • @KBS Since $(I+PR)^{-1}=I-P(I+RP)^{-1}R$, we can find that $X=A^TX(I+BB^TX)^{-1}A$ is equivalent to $X=A^TXA-A^TXB(I+B^TXB)^{-1}B^TXA$. – Lee May 09 '22 at 00:35
  • @Lee In that case, why not using the explicit solution for the DARE? – KBS May 09 '22 at 10:49
  • @KBS what is the explicit solution to DARE? I know that there exists a closed form solution only for the case when all eigenvalues of $A$ are distinct. – Lee May 10 '22 at 08:28
  • 1
    @Lee As far as I remember, the only condition I remember is that $A$ is nonsingular. This is for instance what is considered on Wikipedia and in various books on Riccati equations. – KBS May 10 '22 at 08:39
  • 1
    @KBS For fixed $A$ and $B$, I have to solve DARE using Matlab right now, which I want to avoid. And I don't know any closed form solution that I can use. If you have any information on closed form solution, it would really help me. I found only one paper, where closed form solution is provided when all eigenvalues of $A$ are distinct – Lee May 10 '22 at 08:45
  • @Lee It depends what you mean by closed-form solution but check the Wikipedia page, as already mentioned, or the book by Abou Kandil et al. – KBS May 10 '22 at 08:47

1 Answers1

2

$ \def\LR#1{\left(#1\right)} \def\BR#1{\Big(#1\Big)} \def\vc#1{\operatorname{vec}\LR{#1}} \def\Mat#1{\operatorname{Mat}\LR{#1}} \def\trace#1{\operatorname{tr}\LR{#1}} \def\qiq{\quad\implies\quad} \def\p{\partial} \def\grad#1#2{\frac{\p #1}{\p #2}} $Since $X$ is symmetric and positive definite, it has an inverse. For typing convenience, define the matrix variables $$\eqalign{ V &= A^{-1} \\ Y &= X^{-1} \,= Y^T \\ C &= BB^T = C^T \\ Q &= VC \\ M &= {A\otimes I-I\otimes V} \\ P &= \LR{I-M^+M} \,= P^T \\ }$$ and use them to transform the constraint $$\eqalign{ X &= A^TX(I+CX)^{-1}A \\ YV^TXV &= (I+CX)^{-1} \\ AYA^TX &= I + CX \\ A^TX &= XV + XQX \\ YA^T &= VY + Q \\ Q &= YA^T - VY \\ }$$ into a Sylvester equation for $Y=X^{-1}$ which can be vectorized and solved as a linear equation $$\eqalign{ \vc{Q} &= M\,\vc{Y} \\ \vc{Y} &= M^+\vc{Q} + P\,\vc{U} \\ }$$ where $M^+$ denotes the pseudoinverse of $M,\,$ and $U$ is an unconstrained matrix variable.

I'm getting tired of writing $\vc{}$ so let's define the vectors $$\eqalign{ x &= \vc{X},\quad &c = \vc{C},\quad &q = \vc{Q},\quad &u = \vc{U} \\ X &= \Mat{x},\quad &C = \Mat{c},\quad &Q = \Mat{q},\quad &U = \Mat{u} \\ y &= M^+q + Pu,\quad &dy = P\,du \\ }$$ Write the differential of $X$ in vector form $$\eqalign{ dX &= -X\,dY\,X \\ dx &= -\LR{X\otimes X}\,dy \;=\; -\LR{X\otimes X}P\,du \\ }$$

Now revisit the objective function and calculate its gradient with respect to $u$. $$\eqalign{ \phi &= \trace{CX} \\ &= c^Tx \\ d\phi &= c^Tdx \\ &= -c^T {\LR{X\otimes X}P\,du} \\ &= -du^T\;{P\LR{X\otimes X}c} \\ \grad{\phi}{u} &= -P\LR{X\otimes X}c \;\doteq\; G(u) \\ }$$ Solve for the value of $u$ which makes $G(u)=0$ using your favorite gradient descent algorithm.

Keep in mind that $$X \,=\, Y^{-1} \,=\, \BR{\Mat{M^+q+Pu}}^{-1}$$ so ultimately, everything is a function of the $u$ vector.

The gradient descent step will look like $$u_{k+1} = u_k - \lambda\,G(u_k)$$ where $k$ is the iteration counter and $\lambda$ is the step-length which will be dictated by the algorithm that you choose.

greg
  • 40,033
  • thanks for your answer. I am trying to understand it. Given $A$ and $B$, I need to first compute $M$ and $P$ to get the $G(u)$? in that case I still need to calculate $X$ – Lee May 14 '22 at 16:05
  • Look closely, $(M,P,V)$ are functions of $A$ only. Likewise $(C,c)$ are functions of $B$ only. All of these quantities can be calculated once before the iterations are started. – greg May 14 '22 at 16:17
  • $\ldots$as for $Q,$ it is (ultimately) a function of both $A$ and $B$. – greg May 14 '22 at 16:23
  • when $M$ is nonsingular, then $P$ is zero matrix, thus $G(u)$ zero vector? – Lee May 14 '22 at 17:56
  • If $P$ is equal to zero, that means that you have used up all of your degrees of freedom satisfying the constraint and there's nothing left to optimize. – greg May 14 '22 at 18:26
  • Actually, since $A$ is invertible then $M$ is, too. So I guess I don't understand your constraint. Given $A$ and $B,,$ the constraint has a unique solution, which implies that it is the solution to the overall problem $$X = {\rm Mat}\left(M^{-1}q\right)^{-1}$$ – greg May 14 '22 at 18:40
  • The trivial solution $,X=0,$ also satisfies the constraint. – greg May 14 '22 at 19:47