Fast solution of quadratic matrix equation

Hi everyone,

I have a quadratic matrix equation of the form XCX+X+D=0. Currently I am solving it using the code kindly provided at https://www.mathworks.com/matlabcentral/fileexchange/26956-solve-bilateral-matrix-quadratic-equation . Since I have to perform this operation repetitively, my code (which involves multiple calls to this function) takes around an hour to run. I'm running a parallel pool but I'm limited by the number of cores.

I was wondering if there is any faster way of solving the equation?

Thanks

[John D'Errico]

That code does seem a bit inefficient, using a rather simple iteration scheme. I'd take a read through this to see if they offer a better solution.

https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=3&ved=0ahUKEwjl9-6iua_YAhUB74MKHUXPANwQFgg9MAI&url=http%253A%252F%252Fwww.maths.manchester.ac.uk%252F~higham%252Fnarep%252Fnarep347&usg=AOvVaw23sbdmt8fEHi5g4o7CN3bd

Then I'd look at this alternative code from the FEX:

https://www.mathworks.com/matlabcentral/fileexchange/36263-algebraic-riccati-equation-solver

[Jan]

The code in the posted link does not run: The function sylv() is missing. What are the sizes of your inputs? It matters if you run millions of iterations on a 100x100 system, or 100 iterations on a 1e5x1e5 system.

[Umer Abdullah]

Thanks guys.

Jan, my apologies I forgot to mention that I had to change the sylv() function to lyap() in order to make it run. The sizes are variable, but maximum 50x50. But I have to call the function around 5 million times.

[Jan]

@Umer: Can you provide a working version with some test input? I tried it with sylvester, but it did not converge.

[Umer Abdullah]

@Jan Thanks for the help in this. I'm attaching the code, a sample input set is as follows

Amat=[0.0100,0;0,0.0100];

Bmat=[0,0;0,0];

Cmat=[0.0312,0.0129;0.0129,0.0054];

Dmat=[-0.0015,-0.0036;0.0036,-0.0090];

result=Solve_BQMEv4(Amat,Bmat,Cmat,Dmat);

Answer 1

[David Goodmanson]

Hi Umer,

If matrix C is well behaved enough, meaning not too large a condition number, then the following might work

[v lambda] = eig(D*C,'vector');
r = (-1 + sqrt(1-4*lambda))/2;
x = (v*diag(r))/(C*v);

For 5 million calls, all 50x50, this would take about 4 hours on my pc, but I only have the one processor.

For an NxN matrix there are lots of solutions, since for the + sign in front of the sqrt you have an independent choice of +- for each element of lambda, so 2^N in all.

[Umer Abdullah]

Hi David,

Thanks a lot for this. Indeed a closed form solution seems more elegant.