<p>This was meant as a comment, but too long.</p>
<p>If we fix matrix norm to be 1, we can ask - which algorithms terminate faster when smallest eigenvalue is far from zero?</p>
<p>We can rule out algebraic algorithms like in Ryan Williams <a href="https://mathoverflow.net/questions/24287/what-is-the-best-algorithm-to-find-the-smallest-nonzero-eigenvalue-of-a-symmetric/24294?_gl=1*tmjz6g*_ga*MjU2MDkyNzE2LjE2ODcwNzA4Mjc.*_ga_S812YQPLT2*MTY5MTE2MTA5NC45MC4xLjE2OTExNjI3OTUuMC4wLjA.#24294">answer</a> -- those perform the same number of steps regardless of numeric values. Instead, only consider iterative algorithms which continue until certain numeric error criterion is met.</p>
<p>My feeling is that many such algorithms terminate faster for well-conditioned matrices.</p>
<p>Here's a sketch of how to show this for linear equation solving. The trick is to reformulate gradient descent as operating on the space of eigenvalues rather than dimensions.</p>
<p>The the problem of solving <span class="math-container">$Ax=b$</span> using gradient descent with learning rate <span class="math-container">$1$</span> on least squares objective. Since gradient descent is rotationally symmetric we can take <span class="math-container">$H=A^T A$</span> to be diagonal WLOG, with diagonal entries <span class="math-container">$\{h_1,\ldots,h_d\}$</span>. Suppose we care about relative loss rather than absolute, hence can take <span class="math-container">$\|A\|=1$</span>. If the initial loss is bounded, we have a bound on <span class="math-container">$(x^*_i-x^0_i)^2$</span> for solution <span class="math-container">$x^*$</span> and initial guess <span class="math-container">$x^0$</span>. For simplicity, assume this quantity is exactly <span class="math-container">$1$</span>.</p>
<p>After <span class="math-container">$t$</span> steps, our least squares loss is the following</p>
<p><span class="math-container">$$f(t)=\sum_{i=1}^d h_i (1-h_i)^{2t}$$</span></p>
<p>We can write this sum down as integral over <span class="math-container">$\mathrm{d}h$</span>, the empirical density of <span class="math-container">$h_i$</span></p>
<p><span class="math-container">$$f(t)=\int \mathrm{d}h\ h (1-h)^{2t}$$</span></p>
<p>You can view this as a convolution with a kernel that is moving to the left as <span class="math-container">$t$</span> increases. Fixing range of the spectrum fixes the range of <span class="math-container">$t$</span> for which the loss is changing. For instance, visualize this for <span class="math-container">$h=\{1/200,1/20,1/2\}$</span>
<a href="https://i.sstatic.net/leoDw.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/leoDw.png" alt="enter image description here" /></a></p>
<p>Since we have norm <span class="math-container">$1$</span>, upper bound on condition gives a lower bound on  smallest eigenvalue <span class="math-container">$\lambda_d$</span>. The &quot;gradient descent kernel&quot; is centered at <span class="math-container">$\frac{1}{2t}$</span>, hence one would expect for the error to decay  exponentially with <span class="math-container">$t$</span> after <span class="math-container">$t&gt;\frac{1}{2\lambda_d}$</span> steps, hence needing <span class="math-container">$O\left(\frac{\log \epsilon }{2\lambda_d}\right)$</span> steps for accurate approximation</p>
<hr />
<p>Regarding Lyapunov's equation <span class="math-container">$AX+XA'=B$</span> with Hermitian <span class="math-container">$A$</span> associated with a matrix of eigenvectors <span class="math-container">$U$</span> and column vector of eigenvalues <span class="math-container">$s$</span>, you can write solution as follows (<a href="https://math.stackexchange.com/questions/3354913/solution-of-axxa-b-through-eigenvectors-of-a/3355045#3355045">proof</a>):</p>
<p><span class="math-container">$$X=U \left( \frac{U' BU}{s + s'} \right) U'$$</span></p>
<p>If our application doesn't care about choice of basis, we would just need to compute the central part</p>
<p><span class="math-container">$$\frac{U' BU}{s + s'} $$</span></p>
<p>Now the question would be -- is this cheap to approximate for well-conditioned <span class="math-container">$A$</span>?</p>

CC BY-SA 4.0
2023-08-04T16:55:54.503Z
2023-08-04T17:02:33.690Z
2023-08-04T17:02:33.690Z
Yaroslav Bulatov
<p>Using the <a href="https://en.wikipedia.org/wiki/Conjugate_gradient_method" rel="noreferrer">Conjugate gradient method</a> we can solve a linear system <span class="math-container">$Ax=b$</span>, where <span class="math-container">$A\in\mathbb R^{n\times n}$</span> in time <span class="math-container">$O(n^2 \sqrt{\kappa})$</span>, where <span class="math-container">$\kappa=\frac{\sigma_\mathrm{max}(A)}{\sigma_\mathrm{min}(A)}$</span> is the condition number.
For well conditioned matrices this is better than inverting <span class="math-container">$A$</span> in <span class="math-container">$O(n^\omega)$</span> time (as well as various matrix factorization algorithms.)</p>
<p>My intuition is that well conditioned matrices are in some sense trivial, and this should be a benefit, not just for linear solving, but also other matrix related problems.
In particular I'm interested in whether there are condition number dependent algorithms (like <span class="math-container">$O(n^{2} f(\kappa))$</span> time) for any of the following problems:</p>
<ul>
<li>Matrix Inverse.</li>
<li>Solving the <a href="https://en.wikipedia.org/wiki/Lyapunov_equation" rel="noreferrer">Lyapunov equation</a>.</li>
<li>Matrix determinant.</li>
</ul>
<p>For the first, of course we can solve <span class="math-container">$n$</span> linear systems, but that's <span class="math-container">$n^3\sqrt{\kappa}$</span> time.</p>
<p>The Lyapunov matrix equation (or <a href="https://en.wikipedia.org/wiki/Sylvester_equation" rel="noreferrer">Sylvester Equation</a>) <span class="math-container">$A\Sigma + \Sigma A^T = R$</span> can be solved by the alternating-direction implicit (ADI) method.
But I'm not sure how fast that actually is. In <a href="http://www.dm.unibo.it/%7Esimoncin/sirev58-3_377.pdf" rel="noreferrer">this paper</a> (Theorem 6 and 7) they seem to suggest a condition number related bound, but unless I misunderstand something, they need roughly <span class="math-container">$n-\varepsilon/\kappa$</span> iterations (each <span class="math-container">$n^2$</span> time). Not really better than <span class="math-container">$n^3$</span>.</p>
<p>For matrix determinant it seems <a href="https://mathoverflow.net/questions/24287/what-is-the-best-algorithm-to-find-the-smallest-nonzero-eigenvalue-of-a-symmetric/24294#24294">from this mathoverflow answer</a> that <span class="math-container">$n^\omega$</span> is not even possible. So maybe it's unrealistic to beat <span class="math-container">$n^3$</span> here, even though some structured matrices can be solved faster.</p>

CC BY-SA 4.0
2023-11-03T13:35:27.790Z
BTW, answer [here](https://scicomp.stackexchange.com/a/43427/18786) gives a way to bound convergence of CG linear solver when $A$ matrix eigenvalues are contained in two clusters of width $\epsilon$. High condition number is a special case of this analysis with eigenvalues contained in just one cluster of small width
Yaroslav Bulatov
<cc.complexity-theory><ds.algorithms><reference-request><approximation-algorithms><matrices>
Condition Number dependent algorithms for matrix operations
Thomas Ahle
# Condition Number dependent algorithms for matrix operations

(Asked by: username_0 on 2023-08-02)

Using the [Conjugate gradient method](https://en.wikipedia.org/wiki/Conjugate_gradient_method) we can solve a linear system $Ax=b$, where $A\in\mathbb R^{n\times n}$ in time $O(n^2 \sqrt{\kappa})$, where $\kappa=\frac{\sigma\_\mathrm{max}(A)}{\sigma\_\mathrm{min}(A)}$ is the condition number.
For well conditioned matrices this is better than inverting $A$ in $O(n^\omega)$ time (as well as various matrix factorization algorithms.)

My intuition is that well conditioned matrices are in some sense trivial, and this should be a benefit, not just for linear solving, but also other matrix related problems.
In particular I'm interested in whether there are condition number dependent algorithms (like $O(n^{2} f(\kappa))$ time) for any of the following problems:

* Matrix Inverse.
* Solving the [Lyapunov equation](https://en.wikipedia.org/wiki/Lyapunov_equation).
* Matrix determinant.

For the first, of course we can solve $n$ linear systems, but that's $n^3\sqrt{\kappa}$ time.

The Lyapunov matrix equation (or [Sylvester Equation](https://en.wikipedia.org/wiki/Sylvester_equation)) $A\Sigma + \Sigma A^T = R$ can be solved by the alternating-direction implicit (ADI) method.
But I'm not sure how fast that actually is. In [this paper](http://www.dm.unibo.it/%7Esimoncin/sirev58-3_377.pdf) (Theorem 6 and 7) they seem to suggest a condition number related bound, but unless I misunderstand something, they need roughly $n-\varepsilon/\kappa$ iterations (each $n^2$ time). Not really better than $n^3$.

For matrix determinant it seems [from this mathoverflow answer](https://mathoverflow.net/questions/24287/what-is-the-best-algorithm-to-find-the-smallest-nonzero-eigenvalue-of-a-symmetric/24294#24294) that $n^\omega$ is not even possible. So maybe it's unrealistic to beat $n^3$ here, even though some structured matrices can be solved faster.

**Tags:** cc.complexity-theory, ds.algorithms, reference-request, approximation-algorithms, matrices
> **Comments:**
> **username_1**: BTW, answer [here](https://scicomp.stackexchange.com/a/43427/18786) gives a way to bound convergence of CG linear solver when $A$ matrix eigenvalues are contained in two clusters of width $\epsilon$. High condition number is a special case of this analysis with eigenvalues contained in just one cluster of small width  

---

## 1 answers

### Answer by username_1 on 2023-08-04. Score: 1

This was meant as a comment, but too long.

If we fix matrix norm to be 1, we can ask - which algorithms terminate faster when smallest eigenvalue is far from zero?

We can rule out algebraic algorithms like in Ryan Williams [answer](https://mathoverflow.net/questions/24287/what-is-the-best-algorithm-to-find-the-smallest-nonzero-eigenvalue-of-a-symmetric/24294?_gl=1*tmjz6g*_ga*MjU2MDkyNzE2LjE2ODcwNzA4Mjc.*_ga_S812YQPLT2*MTY5MTE2MTA5NC45MC4xLjE2OTExNjI3OTUuMC4wLjA.#24294) -- those perform the same number of steps regardless of numeric values. Instead, only consider iterative algorithms which continue until certain numeric error criterion is met.

My feeling is that many such algorithms terminate faster for well-conditioned matrices.

Here's a sketch of how to show this for linear equation solving. The trick is to reformulate gradient descent as operating on the space of eigenvalues rather than dimensions.

The the problem of solving $Ax=b$ using gradient descent with learning rate $1$ on least squares objective. Since gradient descent is rotationally symmetric we can take $H=A^T A$ to be diagonal WLOG, with diagonal entries $\{h\_1,\ldots,h\_d\}$. Suppose we care about relative loss rather than absolute, hence can take $\|A\|=1$. If the initial loss is bounded, we have a bound on $(x^\*\_i-x^0\_i)^2$ for solution $x^\*$ and initial guess $x^0$. For simplicity, assume this quantity is exactly $1$.

After $t$ steps, our least squares loss is the following

$$f(t)=\sum\_{i=1}^d h\_i (1-h\_i)^{2t}$$

We can write this sum down as integral over $\mathrm{d}h$, the empirical density of $h\_i$

$$f(t)=\int \mathrm{d}h\ h (1-h)^{2t}$$

You can view this as a convolution with a kernel that is moving to the left as $t$ increases. Fixing range of the spectrum fixes the range of $t$ for which the loss is changing. For instance, visualize this for $h=\{1/200,1/20,1/2\}$
[![enter image description here](https://i.sstatic.net/leoDw.png)](https://i.sstatic.net/leoDw.png)

Since we have norm $1$, upper bound on condition gives a lower bound on smallest eigenvalue $\lambda\_d$. The "gradient descent kernel" is centered at $\frac{1}{2t}$, hence one would expect for the error to decay exponentially with $t$ after $t>\frac{1}{2\lambda\_d}$ steps, hence needing $O\left(\frac{\log \epsilon }{2\lambda\_d}\right)$ steps for accurate approximation

---

Regarding Lyapunov's equation $AX+XA'=B$ with Hermitian $A$ associated with a matrix of eigenvectors $U$ and column vector of eigenvalues $s$, you can write solution as follows ([proof](https://math.stackexchange.com/questions/3354913/solution-of-axxa-b-through-eigenvectors-of-a/3355045#3355045)):

$$X=U \left( \frac{U' BU}{s + s'} \right) U'$$

If our application doesn't care about choice of basis, we would just need to compute the central part

$$\frac{U' BU}{s + s'} $$

Now the question would be -- is this cheap to approximate for well-conditioned $A$?


---

