Is it possible to use the Laplace Transform to calculate eigenvalues?

Question

The relationship of Eigenvalues with Gradient Descent

Let $A$ be a symmetric (and thus diagonalizable) matrix, with diagonalization

$$A=VDV^T.$$

Let us define the quadratic function

$$f(x) = x^T A x.$$

Then gradient descent on $f$ for starting point $x_0$ with learning rate $h$ can be written as

$$x_n = x_{n-1} -h\nabla f(x_{n-1}) = (I-hA)x_{n-1}.$$

Or to highlight the eigenvalues

$$V^T x_n = (I-hD)V^Tx_{n-1}.$$

Taking an eigenvector $v_i$ with eigenvalue $\lambda_i$ out of $V$ this implies

$$x^{(i)}_n := \langle v_i, x_n\rangle =(1-h\lambda_i)x_{n-1}^{(i)} = (1-h\lambda_i)^nx_0^{(i)}.$$

Relation to the Laplace Transform

Now this means that for the discrete measure

$$\mu(x) = \sum_{i=1}^d x_0^{(i)}1_{\{x=-\log(1-h\lambda_i)\}}$$

we get some evaluations of the Laplace Transform of $\mu$

$$ \|x_n\|^2 = \sum_{i=1}^d (x_0^{(i)})^2(1-h\lambda_i)^{2n} =\int_{[0,\infty)} \exp(-2n x)d\mu(x) = \mathcal{L}\{\mu\}(2n) $$

Similarly we get

$$ f(x_n)=(x_0^{(i)})^2(1-h\lambda_i)^{2n+1} = \mathcal{L}\{\mu\}(2n+1) $$

Inverting the Laplace Transform somehow?

If we could invert our Laplace transform we would of course immediately get the eigenvalues as the non-zero "frequencies" in our discrete measure $\mu$. And since the discrete FFT blurs frequencies over adjacent ones if you do not hit them exactly I would expect this to at least provide a useful insight into the upper and lower bound (i.e. condition number) with very little precision. Since the condition number can be used for calculating optimal learning rates, this would be very helpful in accelerating gradient descent (although in that case we would have to do with function evaluations only, since we generally do not know where the minimum is in contrast to this centred example where it is always in 0).

Now the issue is of course that we do not get arbitrary evaluations of our Laplace transform and the Inverse Laplace Transform demands evaluations on the complex axis.

But maybe someone else has an idea how to salvage this situation. Maybe there is a less efficient algorithm than the Z-Transform which can do with those evaluations only, maybe we can reformulate gradient descent a bit to get a representation which allows complex evaluations.

I thought this insight was neat and Laplace Transform would be considerably faster than current eigenvalue calculation algorithms if this could be resolved - so I wanted to throw this out here.

Answer 1

A very unsatisfying answer is matrix exponentiation (which is very expensive itself):

$$x_n=\exp(ihA)x_{n-1}$$

results in

$$ x_n^{(j)} = \exp(ih\lambda_jn) x_0^{(j)} $$

and thus

$$ \|x_n\|^2 = \sum_{j=1}^d (x_0^{(j)})^2 \exp(i2hn\lambda_i) $$

This process seems to be basically the idea of Quantum Phase Estimation.

Answer 2

There's indeed a connection. First you can get rid of $x_0$ as follows -- starting with random $x_0$ initialized on a sphere, expected squared error $f(t)$ at time $t$ can be shown to be

$$f(t)=\operatorname{Tr}((I-A)^{2t})$$

In high dimensions, this is not only an average value, but also a "typical value", almost all trajectories will concentrate around this average.

Now if we use $\exp(-t)\approx (I-A)^t$ approximation we can write

$$f(t)=E[e^{-2 \lambda t}]$$

Where expectation is taken over empirical eigenvalue density of $A$ eigenvalues.

This is the moment generating function, so we can invert it symbolically using inverse Laplace transform. Given eigenvalue density, we can get formula for $k$th eigenvalue from the CDF of this density.

A complete example is worked out here for the harder case when average loss rather than average error is observed.

The issue with this approach is that it requires explicit symbolic form of $f(t)$, it's not clear how to adapt it to the case when $f(t)$ is known at a number of discrete intervals.