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)

https://colab.research.google.com/drive/113XQ88puFNNSQHPnO-zoAgm-8w01ANuh#scrollTo=n6OQHQ5jYw53

Background

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
2. Compute full spectrum of $X$ efficiently for visualization purposes

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)

https://colab.research.google.com/drive/113XQ88puFNNSQHPnO-zoAgm-8w01ANuh#scrollTo=n6OQHQ5jYw53

Background

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

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)

https://colab.research.google.com/drive/113XQ88puFNNSQHPnO-zoAgm-8w01ANuh#scrollTo=n6OQHQ5jYw53

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)

https://colab.research.google.com/drive/113XQ88puFNNSQHPnO-zoAgm-8w01ANuh#scrollTo=n6OQHQ5jYw53

Background

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
2. Compute full spectrum of $X$ efficiently for visualization purposes

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. An example loglog (base10) plot of the spectrum of X is below -- singular values cluster bottom out around 1/2, with a cliff-like fall-off in the tail. I'm curious if this fall-off is an artifact of A, C being low rank.

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)

https://colab.research.google.com/drive/113XQ88puFNNSQHPnO-zoAgm-8w01ANuh#scrollTo=n6OQHQ5jYw53