numerics - Solving underdetermined Lyapunov equations?

I'm wondering if there's an efficient way to get a solution (ie, LeastSquaressolution) for Lyapunov equation $AX+XA=C$ with symmetric positive definite $ A $ and $ C $.

I want something that would work like LyapunovSolve, but would work for underconstrained problems, i.e., LyapunovSolve[A, A] should give me something whose spectrum looks like $ I $.

I tried a naive approach which is to do Kronecker expansion followed by LeastSquares, which gives the desired result

kronExpand[x_] := Module[{ii},
   ii = IdentityMatrix[First[Dimensions[x]]];
   ii\[TensorProduct]x + Transpose[x]\[TensorProduct]ii
   ];

lyapLeastSquares[A_, B_] := Module[{d, X},
  X = LeastSquares[kronExpand[A], vec[B]];
  X = unvec[X, d];

(check this notebook for an end-to-end example)

However, this expansion is too large to be practical. IE, my matrices are on the order of 1000 which is fast using LyapunovSolve, but doing Kronecker expansion means that I have matrices on the order of 1M-by-1M. Any suggestions on how to make this feasible?