matrix - Solving "Lyapunov-like" equation AX+X'B=C

Is there some way I can solve the following equation with $d-by-d$ matrices in Mathematica in reasonable time?

$$AX+X'B=C$$

My solution below calls linsolve on $d^2,d^2$ matrix, which is too expensive for my case (my d is 1000)

kmat[n_] := Module[{mat1, mat2},
   mat1 = Array[{#1, #2} &, {n, n}];
   mat2 = Transpose[mat1];
   pos[{row_, col_}] := row + (col - 1)*n;
   poses = Flatten[MapIndexed[{pos[#1], pos[#2]} &, mat2, {2}], 1];
   Normal[SparseArray[# -> 1 & /@ poses]]
   ];
unvec[Wf_, rows_] := Transpose[Flatten /@ Partition[Wf, rows]];
vec[x_] := Flatten[Transpose[x]];


solveLyapunov2[a_, b_, c_] := Module[{},
  dims = Length[a];
  ii = IdentityMatrix[dims];
  x0 = LinearSolve[
    KroneckerProduct[ii, a] + 
     KroneckerProduct[Transpose[b], ii].kmat[dims], vec[c]];
  X = unvec[x0, dims];
  Print["error is ", Norm[a.X + Transpose[X].b - c]];
  X

  ]

a = RandomReal[{-3, 3}, {3, 3}];
b = RandomReal[{-3, 3}, {3, 3}];
c = RandomReal[{-3, 3}, {3, 3}];
X = solveLyapunov2[a, b, c]

Edit Sep 30: An approximate solution would be useful as well. In my application $C$ is the gradient, and $X$ is the preconditioned gradient, so I'm looking for something that's much better than a "default" solution of $X_0=C$

Answer

General matrices

For the desired matrix sizes I have doubts that a numerical solution would be feasible. Here is a simplified code using sparse matrices.

tmSylvester[n_]:=Module[{a,b,c,sA,sB,sC,sAB},
a=RandomReal[{-3,3},{n,n}];
b=RandomReal[{-3,3},{n,n}];
c=RandomReal[{-3,3},{n,n}];
sA=SparseArray[Table[{(i-1)n+l,(k-1)n+l}->a[[i,k]],{i,n},{k,n},{l,n}]//Flatten];
sB=SparseArray[Table[{(l-1)n+j,(k-1)n+l}->b[[k,j]],{k,n},{j,n},{l,n}]//Flatten];
sAB=sA+sB;
sC=SparseArray[Table[{(i-1)n+j}->c[[i,j]],{i,n},{j,n}]//Flatten];
First[Timing[LinearSolve[sAB,sC];]]]

Now, let us plot the timing

ListLogPlot[Table[{n,tmSylvester[n]},{n,10,120,10}],Joined->True,PlotTheme->{"Frame","Monochrome"}, FrameLabel->{"Matrix Size","Time(s)"}]

Even at a very optimistic extrapolation it is unlikely that the n=1000 calculation would be routinely possible. There are, however, experts here that might be able to further tune up the linear solver.

Nonsingular matrices

According to F. M. Dopico, J. González, D. Kressner, and V. Simoncini. Projection methods for large-scale T-Sylvester equations, in Mathematics of Computation (2015), under the usual conditions of existence the following equations have equal unique solutions

$$􏰁B^{−T} A􏰂 X − X 􏰁A^{−T} B􏰂 = B^{−T} C − B^{−T} C^{T} A^{−T} B;$$ $$AX + X^T B = C, $$ where $A^{-T}\equiv(A^{-1})^T$.

Therefore, we can use the Lyapunov solver

tmDopico[n_]:=Module[{a,b,c},
a=RandomReal[{-3,3},{n,n}];
b=RandomReal[{-3,3},{n,n}];
c=RandomReal[{-3,3},{n,n}];
First[Timing[LyapunovSolve[Transpose[Inverse[b]].a,-Transpose[Inverse[a]].b,Transpose[Inverse[b]].c-Transpose[Inverse[b]].Transpose[c].Transpose[Inverse[a]].b];]]]

Let us check the timing:

ListLogPlot[Table[{n,tmDopico[n]},{n,50,1000,50}],Joined->True,PlotTheme->{"Frame","Monochrome"}, FrameLabel->{"Matrix size","Time(s)"}]

The method should therefore have $\mathcal{O}(n^3)$ scaling under favorite conditions.

Blog

Search This Blog