# Solve linear system under Hermitian constraint

**URL:** <https://discourse.julialang.org/t/solve-linear-system-under-hermitian-constraint/130158>\
**Category:** Numerics\
**Created:** [June 24, 2025, 9:23am UTC](https://discourse.julialang.org/t/solve-linear-system-under-hermitian-constraint/130158 "2025-06-24T09:23:20Z")\
**Posts on this page:** 6\
**Page:** 2

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [June 24, 2025, 4:25pm UTC](https://discourse.julialang.org/t/solve-linear-system-under-hermitian-constraint/130158/21 "2025-06-24T16:25:34Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 24, 2025, 4:36pm UTC](https://discourse.julialang.org/t/solve-linear-system-under-hermitian-constraint/130158/22 "2025-06-24T16:36:14Z")

</div>

> [@stevengj](#):
>
> However, note that this method is O(m^3), whereas the minimum-norm iterative solution from above is O(mn^2) (modulo convergence rate). This might matter for large m.

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?

---

<div class="post-metadata">

**Author:** ![araujoms](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/araujoms/32/217734_2.png) [@araujoms](https://discourse.julialang.org/u/araujoms)\
**Post date:** [June 24, 2025, 4:57pm UTC](https://discourse.julialang.org/t/solve-linear-system-under-hermitian-constraint/130158/23 "2025-06-24T16:57:35Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [June 24, 2025, 5:05pm UTC](https://discourse.julialang.org/t/solve-linear-system-under-hermitian-constraint/130158/24 "2025-06-24T17:05:11Z")

</div>

Thanks for posting this, it was a fun rabbit hole for me. This is actually something I wondered before.

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [June 26, 2025, 12:53pm UTC](https://discourse.julialang.org/t/solve-linear-system-under-hermitian-constraint/130158/25 "2025-06-26T12:53:02Z")

</div>

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`](https://github.com/andreasvarga/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
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

```

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [June 30, 2025, 1:26pm UTC](https://discourse.julialang.org/t/solve-linear-system-under-hermitian-constraint/130158/26 "2025-06-30T13:26:05Z")

</div>

> [@stevengj](#):
>
> AΛ†A+ΛA†A=B

This equation can be also solved using the iterative CG-solver `ghsylvi` from [MatrixEquations.jl](https://github.com/andreasvarga/MatrixEquations.jl)

`@time Λ, info = ghsylvi([I],[A'*A],[A],[A],B,reltol=1.e-14)`

[Previous page](https://discourse.julialang.org/t/solve-linear-system-under-hermitian-constraint/130158.md?page=1)