# How can I find the eigenvectors of a Hamiltonian matrix to solve a Riccati equation?

Given the Algebratic Riccati Equation (ARE) $$A^T X + XA + XRX + Q = 0$$ where $$A,R,Q \in \mathbb R^{n \times n}$$, we are interested in the matrix $$X$$ that solves this equation. If we define the $$2n \times 2n$$ Hamiltonian matrix $$H = \begin{bmatrix} A & R \\ -Q & -A^T\end{bmatrix}$$ and let $$X = X_2 X_1^{-1}$$, then $$H$$ satisfies the equation $$H\begin{bmatrix} X_1 \\ X_2\end{bmatrix} = \begin{bmatrix} X_1 \\ X_2\end{bmatrix} \Lambda$$ which represents an eigenvalue problem. That is, if $$v_i \in \mathbb R^{2n}$$ is the $$i$$th column of $$\begin{bmatrix} X_1 \\ X_2\end{bmatrix}$$, and $$\lambda_i$$ is the $$i$$th corresponding eigenvalue, then $$H$$ satisfies the equation $$Hv_i = \lambda_i v_i$$ I'm aware that there are better ways of solving this ARE, but I'm interested in solving this eigenvalue problem to obtain $$X_1$$ and $$X_2$$, and then $$X$$. However, when computing the eigenvalues and eigenvectors of $$H$$ using, for example, NumPy, I obtain a $$2n \times 2n$$ matrix of eigenvectors. I'm not sure how to relate this matrix to $$\begin{bmatrix} X_1 \\ X_2\end{bmatrix}$$, which is $$2n \times n$$.

• The algebraic Riccati usually admits several solutions, and most likely you are interested in finding the so-called stabilizing solution. That is the eigenvalues of the closed-loop matrix $A-BB^T X$ all have negative real-part. If this is the case, You use the eigenvectors of the Hamiltonian corresponding to the eigenvalues with negative real-part as well. Stack all these eigenvectors together and you get the desired $2n\times n$ matrix. Nov 11, 2022 at 10:09
• @DerZwirbel ahh I see. If I’m not necessarily interested in a stabilizing solution, does this mean that any combination of $n$ eigenvectors of $H$, out of a total of $2n$, can be used to construct $\begin{bmatrix} X_1 \\ X_2\end{bmatrix}$, and then from this obtain $X$? Nov 11, 2022 at 11:24
• Roughly Speaking Yes. But It might be the case that $X_1$ is singular for some combinations. But anyway, If $X_1$ is nonsingular, then $X_2 X_1^{-1}$ solves the algebraic Riccati equation. Nov 11, 2022 at 12:34
• @DerZwirbel Please post your contribution as an answer; comments are only meant to suggest improvements to posts. Nov 11, 2022 at 15:04

For reference, here is some Python code implementing the suggestions by @DerZwirbel. The following ARE is solved in this case: $$XA + A^T X - XBR^{-1}B^T X + Q = 0$$

import numpy as np
import scipy.linalg

# generate random test matrices (Q and R need to be positive semi-definite and positive definite respectively)
A = 3 * np.random.randn(2,2)
B = 3 * np.random.randn(2,2)
Q = 3 * np.random.randn(2,2)
Q = Q.T @ Q
R = 3 * np.random.randn(2,2)
R = R.T @ R

# compute the Hamiltonian matrix, its eigenvalues, and its eigenvectors
H = np.block([[A,-B @ np.linalg.solve(R,B.T)],[-Q,-A.T]])
eigs,eigvs = np.linalg.eig(H)

# find out which eigenvalues are negative
stable_idxs = np.where(np.real(eigs)<0)[0]

# split the array of eigenvectors into two separate blocks to create the X_1 and X_2 matrices,
# such that X = X_2 * X_1^{-1}
H_eigvs_top,H_eigvs_bottom = np.split(eigvs,2,axis=0)
X1 = H_eigvs_top[:,stable_idxs]
X2 = H_eigvs_bottom[:,stable_idxs]
X_hat = X2 @ np.linalg.inv(X1)

# as a sanity check, compute the solution using SciPy's implementation and compare
# this solution to the one we obtained manually
X = scipy.linalg.solve_continuous_are(A,B,Q,R,balanced=True)
np.testing.assert_allclose(actual=X_hat,desired=X,err_msg="X and X_hat do not match.")

# evaluate the value of the ARE using X and X_hat (both arrays should be close to 0)
ARE = A.T @ X + X @ A - X @ B @ np.linalg.solve(R,B.T) @ X + Q
ARE_hat = A.T @ X_hat + X_hat @ A - X_hat @ B @ np.linalg.solve(R,B.T) @ X_hat + Q

print(f"Eigenvalues of Hamiltonian matrix:\n{eigs}\n")
print(f"X1:\n{X1}\n")
print(f"X2:\n{X2}\n")
print(f"Ground truth solution X:\n{X}\n")
print(f"Solution X using manual method:\n{X_hat}\n")
print(f"Ground truth value of ARE (should be close to 0):\n{ARE}\n")
print(f"Value of ARE using manual method (should be close to 0):\n{ARE_hat}\n")