Solving Algebraic Riccati Equation by using Newton-Raphson's method?

Question

Assume that we have a function called $$F(X) = 0$$

That function can be written as the Algebraic Riccati Equation:

$$F(X) = A^T X + XA - XBR^{-1}B^TX + Q = 0$$

Where the solution to the equation is matrix $X$.

To do a Newton-Raphson solveing method. I need to use this equation:

$$X_{i+1} = X_{i} - F'(X_{i})^{-1}F(X_{i})$$

Where $F'(X)$ is the jacobian function of $F(X)$ and $X_{i+1}$ is going the be the solution to the Algebraic Riccati Equation.

My question is:

How can I find the derivative of $F(X)$ ? Can I just assume that the derivative is:

$$F'(X) = A^T + A - XBR^{-1}B^T$$ ?

Answer 1

$F'(X):H\in M_n\rightarrow A^TH+HA-HUX-XUH$ where $U=BR^{-1}B^T$.

Let $K=F'(X_i)^{-1}F(X_i)$.

At each step, you must solve in $K$ this linear equation:

$A^TK+KA-KUX_i-X_iUK=A^TX_i+X_iA-X_iUX_i+Q$, or

$(A^T-X_iU)K+K(A-UX_i)=A^TX_i+X_iA-X_iUX_i+Q$, which is a Sylvester equation.

EDIT 1. Answers to the first three comments.

I don't need any dot product or "fourth order tensor"!!! $F'(X)$ is a derivative and not a gradient. Stop to read the matrix cookbook; you don't understand that you do.

$F$ is a function from $M_n$ to $M_n$. Then, its derivative, in $X\in M_n$, $F'(X)$ is a LINEAR application from $M_n$ to $M_n$; we assume that it is a bijection in $X_i$ and, therefore $F'(X_i)(K)=F(X_i)$. One cannot do simpler!

Last point. A Sylvester equation in the unknown matrix $K$, is a linear equation in the form $AK+KB=C$.

cf. https://en.wikipedia.org/wiki/Sylvester_equation

Matlab has a solver dedicated to this type of equation.

EDIT 2. Answer to Daniel Mårtensson . Your below comment: "but where do I update $K$ then?" shows that you did not understand one word of my post...

I rewrite. One iteration consists in that follows.

i) input. $X_i$

ii) Solve in $K$ (using Matlab): $(A^T-X_iU)K+K(A-UX_i)=A^TX_i+X_iA-X_iUX_i+Q$

The equation depends on $X_i$; then $K$ depends on $X_i$.

iii) output. $X_{i+1}=X_i-K$ (and not $X_i+K$).

Answer 2

Let a lowercase letter stand for the corresponding vectorized matrix, e.g. $$f={\rm vec}(F), \,\,\, x={\rm vec}(X)$$ Write the function and its differential $$\eqalign{ F &= A^TX + XA - XPX + Q \cr dF &= A^T\,dX + dX\,A - dX\,PX -XP\,dX \cr }$$ and vectorize $$\eqalign{ df &= (I\otimes A^T+A^T\otimes I-I\otimes XP-X^TP^T\otimes I)\,dx \cr J=\frac{\partial f}{\partial x} &= M - (I\otimes XP) - (PX\otimes I)^T \cr }$$ Note that $$M=(I\otimes A^T+A^T\otimes I)$$ is constant between iterations and only needs to be calculated once.

Then we have the standard Newton iteration $$\eqalign{ Js &= f \cr x_+ &= x - s \cr }$$ All of the Kronecker products are sparse, so by utilizing the sparse solvers from your matrix library, you can solve moderately-sized systems with this simple method.