Suppose I have semi-positive definite matrices $A$ and $C$, is there an efficient approach to get top singular values of X entering the following expression?

$$ AX+XA=C $$

My matrices are 4k-by-4k and are known to be low rank, rank somewhere between 50 and 1000.

I've been using SVD on top of scipy.linalg.svd(scipy.linalg.solve_lyapunov(A, C)), but it's too slow for my application.

The second issue is the issue of stability. Using scipy.linalg to solve $AX+XA=2A$ gives me solution with norm much higher than 1 which I expected since $I=X$ is a valid solution. Since I'm interested in the spectrum, I may need some extra constraints to make this problem well defined. Also, there seems to be significant error introduced (plugging X into the equation and looking at norm gives me relative errors around 2)


enter image description here


These singular values come up in the problem of determining length of step size for linear least-squares estimation, described here: https://arxiv.org/abs/1610.03774

Suppose we are solving linear least squares problem in $d$ dimensions by minimizing the following objective:

$$L(w) = \frac{1}{2}E_{xy}[y-\langle w, x\rangle^2]$$ Residual is defined as $$\epsilon_{x,y}=y-\langle w, x\rangle$$

Hessian of this loss is $$H=E[xx']$$

While covariance matrix of gradients $$\Sigma = E[\epsilon_{xy} xx']$$

The problem is to minimize this loss using stochastic gradient descent. How big can we make the step size? We can show that when errors are uncorrelated with observations, following step size can be taken while still maintaining convergence

$$\gamma = \frac{2}{R^2}$$

where $R^2$ is an upper bound on $\|x\|^2$

However, when errors are correlated (heteroscedactic/misspecified case), this rate must be reduced to the following

$$\gamma = \frac{2}{\rho R^2}$$

Where $1\le \rho\le d$ is a measure of misspecification, and is computed as follows

Let $X$ be the solution of the following $$HX+XH=\Sigma$$

Then $$\rho=d \frac{\|X\|_2}{\text{tr}(X)}$$

The peculiarity of data generating process (these are features generated by neural network) means that $x$ lie in an unknown low-dimensional subspace of $\mathbb{R}^d$ . $d\in[4000,10000]$, but true dimensionality is 50-1000

The task is:

  1. Compute $\rho$ efficiently in order to obtain step size to use for SGD
  2. Compute full spectrum of $X$ efficiently for visualization purposes
  • $\begingroup$ Which part is too slow- the solution of the Lyapunov equation or the SVD that follows? Or, are they taking about equal amounts of time? $\endgroup$ Sep 11, 2019 at 16:33
  • $\begingroup$ What are your time constraints? Are you using an optimized BLAS/LAPACK library? $\endgroup$ Sep 11, 2019 at 16:34
  • $\begingroup$ The lyapunov equation solver is the slow part, several times slower than SVD. It prints warnings about numerical stability and perturbing coefficients, probably not expecting A,C to be low-rank $\endgroup$ Sep 11, 2019 at 16:36
  • 1
    $\begingroup$ You understand that there are nonuniqueness problems when A and C aren't positive definite, right? You may need to regularize the Lyapunov equation. $\endgroup$ Sep 11, 2019 at 16:40
  • 1
    $\begingroup$ I'm afraid you need to back up and ask a different question about how to deal with the Lyapunov equation when $A$ is singular. Start by explaining where you Lyapunov equation came from. $\endgroup$ Sep 11, 2019 at 18:21

2 Answers 2


The speed problem can be fixed: there is literature on methods for large and sparse Lyapunov equations that return $X \approx ZZ^T$ already in factored form (hence sparing you most of the work also in the SVD part). See for instance https://www.mpi-magdeburg.mpg.de/projects/mess and http://www.dm.unibo.it/~simoncin/software.html for the two leading methods.

But, first of all, you need to understand the issue with multiple solutions.

  • $\begingroup$ Thanks for the reference! Is there a standard approach to get the least-squares solution to the Lyapunov equation? (besides using kronecker expansion -> least squares which is expensive) $\endgroup$ Sep 12, 2019 at 18:05
  • $\begingroup$ @YaroslavBulatov Not that I know of. The only thing that comes to mind is that you can use LSQR (a large sparse least-squares algorithm) on the Kronecker form implicitly via its matvec operator. But I doubt that it's going to be competitive; the storage needed seem still quite high. $\endgroup$ Sep 12, 2019 at 18:17
  • $\begingroup$ BTW, Matlab's fast eigenvector based lyap2.m is giving me same result as explicit least squares (had to modify it to skip dividing by 0 eigenvalues) -- github.com/msubhransu/matrix-sqrt/blob/master/lyap2.m (evals wolframcloud.com/obj/yaroslavvb/newton/lyapunov.nb) $\endgroup$ Sep 14, 2019 at 23:49
  • 1
    $\begingroup$ @YaroslavBulatov That's not surprising, because essentially with that method (since $A$ is symmetric) you are computing an SVD of the Kronecker product $I\otimes A + A^T \otimes I$, and then skipping the divisions by 0 gives you the minimum-norm solution, by the usual least squares theory. $\endgroup$ Sep 15, 2019 at 6:35

To address the singularity issue: have you thought about projecting the equation on $\operatorname{range}(A)$? Write $A= Q\hat{A}Q^T$, with $Q\in \mathbb{R}^{d\times k}$ a tall thin matrix with orthogonal columns ($Q^TQ=I$) and $\hat{A}\in\mathbb{R}^{k\times k}$ nonsingular, and accordingly $X=Q\hat{X}Q^T$. Then you can multiply the equation on the left and right by $Q^T$ and $Q$, to get the small-scale equation $$\hat{A}\hat{X}+\hat{X}\hat{A}=\hat{C},$$ with $\hat{C} = Q^TCQ$. This will give you a solution with $\operatorname{range}(X) \subseteq \operatorname{range}(A)$. Not sure if that makes sense in your application.

I did not check, but I would not be surprised if it turned out that this is the same solution as your rank-deficient least-squares approach.

Also, I guess you already have all the pieces available to compute the coefficients of the small-scale equation without even forming a single $d\times d$ matrix, just with a QR of the data matrix containing the vectors $x$.

  • $\begingroup$ Thanks, that approach makes sense. Is there a name for this decomposition of A? For my application, I actually just need the trace(X)/norm(X) for some solution X, which I suspect is the same for all X satisfying the constraints $\endgroup$ Sep 15, 2019 at 23:29
  • $\begingroup$ @YaroslavBulatov I would call it a "low-rank decomposition", generically. It doesn't have a name, as far as I know, also because it's not uniquely defined: you can choose to take $\hat{A}$ diagonal with nonnegative elements, and then it becomes a thin/economy-sized SVD, for instance. I think it's a sort of "folklore idea" in linear algebra. $\endgroup$ Sep 16, 2019 at 6:11
  • $\begingroup$ BTW, I found something relevant in Golub/Loan 5.4.7: Complete Orthogonal Decomposition -- you can get $\hat{A}$ by applying QR with pivoting followed by regular QR $\endgroup$ Sep 16, 2019 at 21:18
  • 1
    $\begingroup$ This is not the least-squares solution. Let $P$ be the orthogonal complement of $Q$. Observe that $Q^TCP \ne 0$ in general, but we can always choose $Q^TXP$ such that $Q^TAXP = Q^TCP$ because $Q^TAXP=Q^TAQQ^TXP$ and $Q^TAQ$ is nonsingular by hypothesis. $\endgroup$ Sep 18, 2019 at 18:28
  • 1
    $\begingroup$ The unique choice of $X$ satisying $Q^TXQ = \hat{X}$ as above and $Q^TXP = (P^TXQ)^T = (Q^TAQ)^{-1}(Q^TCP)$, however, is indeed the least-squares solution, which you can implicitly compute in time linear to dimension $d$ but cubic to rank $k$ $\endgroup$ Sep 18, 2019 at 18:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.