Oh duh. For some reason it looked weird for me to write a Lagrangian with \mathrm{Re}(\langle \Lambda, X - X^\dagger\rangle) in it, but of course if I had written it as \langle \Lambda, X - X^\dagger\rangle + \langle X - X^\dagger, \Lambda\rangle, I would have never batted an eye.
1 Like
I did a quick (single-threaded BLAS) benchmark on my M3 laptop, and the iterative m \times n solve is about 2x faster than the dense lyapc solve for 1000 \times 123, about 10x faster for 2000 \times 123, and about 700x faster for 2000 \times 12. (And, of course, both methods are vastly better than forming explicit Kronecker-product matrices, which quickly runs out of memory as m grows.)
The iterative/Sylvester-transpose form also saves memory (O(mn) vs O(m^2)) if you leave the solution X in the low-rank-factor form X = A \Lambda^\dagger + \Lambda A^\dagger, which can also be multiplied more quickly by vectors than the explicit X.
Not sure what your typical matrix sizes are, though, @araujoms?
1 Like
Right now they are tiny, like 5 \times 3, so efficiency is not a concern. I was interested in learning how to solve the problem properly, and I’m very pleased with the solutions you (and @Mason) found. They might be needed in practice, as matrix sizes do tend to get very large in non-commutative polynomial optimization.
1 Like
Thanks for posting this, it was a fun rabbit hole for me. This is actually something I wondered before.
1 Like
The reason why lyapc works (i.e., it provides almost always a solution), is that the eigenvalues of A*A’ result nonzero, because of rundoff errors. I recently updated MatrixEquations.jl, just to detect the “exact” singularity and issue an error message accordingly.
The computed solutions are in general not of minimum norm. However, a minimum norm solution can be computed using the iterative solver lyapci:
julia> m, n = 1000, 123
(1000, 123)
julia> A = randn(ComplexF64, m,n); B = hermitianpart!(randn(ComplexF64,m,m)) * A;
julia> X = lyapc(A*A', -(A*B'+B*A'));
julia> norm(X*A-B)/norm(X)
1.199116863635282e-14
julia> Xi, = lyapci(A*A', -(A*B'+B*A'),reltol=1.e-10);
julia> norm(Xi*A-B)/norm(Xi)
3.6391141711284985e-9
julia> norm(X)
3.877131194912115e7
julia> norm(Xi)
339.6758498801833
4 Likes
This equation can be also solved using the iterative CG-solver ghsylvi from MatrixEquations.jl
@time Λ, info = ghsylvi([I],[A'*A],[A],[A],B,reltol=1.e-14)
2 Likes