Solving Algebraic Riccati Equation by using Newton-Raphson's method? 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$$ ?
 A: $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$).
A: 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.
