Solving matrix quadratic equation

Question

Let $\mathbf{A}, \mathbf{C}\in \mathbb{R}^{n\times n}$. Given a skew-symmetric matrix $\mathbf{G}$ I am looking for any numerical procedure so solve the following quadratic equation. $$ \mathbf{B}^\top\left[\mathbf{A}^\top \mathbf{C}^\top -\mathbf{C}\mathbf{A}\right] \mathbf{B} + \mathbf{B}^\top \mathbf{A}^\top - \mathbf{A}\mathbf{B} = \mathbf{G} $$ for $\mathbf{B}\in\mathbb{R}^{n\times n}$. Does anyone have suggestions?

Answer 1

Your equation, in the unknown $X=[x_{i,j}]\in M_n$, has the form

(1) $X^TLX+X^TA^T-AX-G=0$ where $L\in M_n$ is skew-symmetric.

It's an algebraic system of $n(n-1)/2$ equations in the $n^2$ unknowns $x_{i,j}$. We can show that these equations are algebraically independent when $L,A,G$ are generic matrices (for example, randomly choose them). Thus the general solution of (1) depends on $n^2-n(n-1)/2=n(n+1)/2$ complex parameters. Anyway, there are too many parameters to be able to calculate them all. It is easier to find particular solutions, for example symmetric matrices. In this case, we are brought back to solve an equation of the form

(2) $XLX+XA^T-AX-G=0$, that is, a Riccati equation. Since $X$ is symmetric, one has $n(n-1)/2$ supplementary relations linking the $(x_{i,j})$; again, these relations are algebraically independent and the general solution of (2) depends on $n(n+1)/2-n(n-1)/2=n$ complex parameters.

There are many solvers for this type of equation. Example for $n=3$

A symmetric solution is

EDIT. Answer to the OP

here is a Maple procedure which provides a particular symmetric solution

restart; with(LinearAlgebra); n := 5; L := RandomMatrix(n, n, shape = skewsymmetric); G := RandomMatrix(n, n, shape = skewsymmetric); A := RandomMatrix(n, n);

M := Matrix(2*n, 2*n, [[-Transpose(A), -L], [-G, -A]]);

roll := rand(-10 .. 10); u := Vector(2*n, proc (i) options operator, arrow; roll() end proc); Y := Matrix([[u]]); v := u; for i to n-1 do w := M.v; Y := Matrix([[Y, w]]); v := w end do; U := SubMatrix(Y, [1 .. n], [1 .. n]); V := SubMatrix(Y, [n+1 .. 2*n], [1 .. n]); X := V.(1/U);