Handling under-constrained Lyapunov equations

AX + XA = C

Background

Datasets X and Y overlap in d dimensions, and are different in 1 dimension. Let A, B be their covariance matrices. How to find d?

Eq 3 of arXiv:1610.03774 suggests that we can recover d from spectrum of X in AX+XA=B

However, this equation is underconstrained so default LyapunovSolve fails

Util

In[]:=

randomRotation[n_]:=Module[{},z=RandomVariate[NormalDistribution[0,1],{n,n}];{q,r}=QRDecomposition[z];d=Diagonal[r];ph=d/Abs[d];q*ph];vec[W_]:=Transpose@{Flatten@Transpose[W]};unvec[Wf_, rows_]:=Transpose[Flatten/@Partition[Wf,rows]];unvec[Wf_]:=Module[{d},d=Sqrt[Dimensions[Wf]//First];Assert[dFloor[d]];Assert[d*dFirst@Dimensions[Wf]];unvec[Wf,d]];CircleTimes=KroneckerProduct;kronExpand[x_]:=Module[{ii},ii=IdentityMatrix[First[Dimensions[x]]];x⊗ii+ii⊗x];randomCov[rank_,d_]:=Module[{},mean=Table[0,{rank}];cov=IdentityMatrix[rank];normal=MultinormalDistribution[mean,cov];X=RandomVariate[normal,{n}]//Transpose;X = X~Join~Array[0&,{d-rank,n}];X=randomRotation[d].X;X.X/n];covPair[sharedRank_,excessRank_,d_,strength_]:=Module[{A,B},shared=randomCov[sharedRank,d];A=shared+strength*randomCov[excessRank,d];B=shared+strength*randomCov[excessRank,d]; {A,B}];erank[A_]:=Tr[A]/Norm[A];

Implementations

(*normalizedrighteigenvectors.*)normalize[vec_]:=vec/Norm[vec];evecs[mat_]:=Transpose[normalize/@Eigenvectors[mat]];norm2[mat_]:=Total[Flatten[mat*mat]];kronExpand[x_]:=Module[{ii},ii=IdentityMatrix[First[Dimensions[x]]];ii⊗x+x⊗ii];lyapLeastSquares[A_,B_]:=Module[{d,X},X=LeastSquares[kronExpand[A],vec[B]];X=unvec[X];Print["errorLsqr=",norm2[A.X+X.A-B]];X];lyapSpectral[A_,B_]:=Module[{U,s,C,cutoff,sdiv,Y},(*AmustbePSD*)U=evecs[A];(*A=U.D.U'*)s=Eigenvalues[A];C=U.B.U;cutoff=Max[s]*

6

10

*$MachineEpsilon;sdiv=Map[If[#>cutoff,1/#,#]&,Outer[Plus,s,s],{2}];Y=C*sdiv;X=U.Y.Transpose[U];Print["errorSpectral=",N@norm2[A.X+X.A-B]];X];lyapSvd[A_,B_]:=Module[{U,s,C,cutoff,sdiv,Y},{U,S,V}=SingularValueDecomposition[A];s=Diagonal[S];C=U.B.V;smat=Outer[Plus,s,s];sinv=Map[If[#>0,1/#,#]&,smat,{2}];Y=C*sinv;X=U.Y.Transpose[V];Print["errorSvd=",N@norm2[A.X+X.A-B]];X];Needs["Notation`"]distinguishMatrixPowerAndPower[A_,d_]:=If[SquareMatrixQ[A],MatrixPower[A,d],Power[A,d]];Notation

d_

A_

⟺

distinguishMatrixPowerAndPower[A_,d_]

lyapGeneric[A_,B_,C_]:=Module{n,ones,Avals,Avecs,Bvals,Bvecs,T,U,DA,DB},n=First@Dimensions[A];Assert[Dimensions[A]=={n,n}];Assert[Dimensions[B]=={n,n}];normalize[x_]:=xNorm[x];ones:=Array[1&,{n,n}];{Avals,Avecs}=Eigensystem[A];{Bvals,Bvecs}=Eigensystem[B];T=Transpose[Avecs];U=Transpose[Bvecs];DA=DiagonalMatrix[Avals];DB=DiagonalMatrix[Bvals];(*Assert[T.DA.Inverse[T]≐A]Assert[U.DB.Inverse[U]≐B]*)X=T.

-1

T

.C.U

DA.ones+ones.DB

.

-1

U

;Print["errorGeneric=",N@norm2[A.X+X.B-C]];(*Symmetricsquareroot*)isymsqrt[m_]:=Module[{U,S,W,IS,svals,isvals},{U,S,W}=SingularValueDecomposition[m];svals=Diagonal[S];isvals=If[#>0,1/#,#]&/@svals;IS=DiagonalMatrix[isvals];U.Sqrt[IS].Transpose[W]];(*SimultaneouslydiagonalizeinvertiblePSDHandsymmetricSusingcongruencetransformation*)simDiagoalize[H_,S_]:=Module[{temp,o,s,is},temp=isymsqrt[H].S.isymsqrt[H];o=Transpose[Eigenvectors[temp]];(*o.DiagonalMatrix[Eigenvalues[temp]].o≐temp;*)s=isymsqrt[H].o;d=s.S.s;(*disdiagonal,traceofd=trace(

-1

H

S)*)Assert[DiagonalMatrixQ[d]];Assert[Tr[S.Inverse[H]]≐Tr[d]];is=Inverse[s];Print["errorH=",Norm[is.is-H]];Print["errorS=",Norm[is.d.is-S]];];(*computess,suchthats.a.sands.b.sarediagonal*)simDiagoalize[a_,b_]:=Module[{s,d},o=Eigenvectors[

-1/2

b

.a.

-1/2

b

]//Transpose;d=o.(

-1/2

b

.a.

-1/2

b

).o;s=

-1/2

b

.o];On[Assert];A={{10,17},{17,29}};B={{26,2},{2,4}};lyapLeastSquares[A,B];lyapSpectral[A,B];lyapSvd[A,B];A={{-1,4},{3,0}};B={{-4,3},{2,6}};C0={{-5,0},{-2,-6}};lyapGeneric[A,B,C0];

Compare spectra of three solutions

all 3 methods seem to produce the same solution

Non - uniqueness of solutions

Larger eval

Lyapunov general

Matlab solution

Transposed Lyapunov

Hadamard equation

Solve AXB + X⊙C = D

Conclusion

All 3 methods seem to give same result for large symmetric semi - definite matrices, but differ when SPD is violated.