# Implementing Gelfand’s formula for the spectral radius in Python - lack of convergence

For context: Gelfand's formula for the spectral radius is $$\lim_{k\rightarrow \infty}|A^k|^{1/k}$$ where $$|\cdot|$$ is any well-defined operator norm.

I naively coded a function to calculate the $$k$$th term in Python as

from numpy import linalg as la

def kth_term(matrix, k):
matrix = la.matrix_power(matrix, k)
f_norm = la.norm(matrix, 'fro')
a = f_norm**(1.0/k)

return a


It appears to converge up to a point, and then veers off the track.

For example, using the matrix

$$\left[ {\begin{array}{ccc} 9 & -1 & 2 \\ -2 & 8 & 4 \\ 1 & 1 & 8 \\ \end{array} } \right]$$

that appears in the wikipedia article on the spectral radius :

k=1: 15.362291495737216
k=2: 12.328294348193777
k=3: 11.532450663575863
k=4: 11.151002985846981
k=5: 10.921242234560514
k=6: 10.766714723560009
k=7: 10.655756642574673
k=8: 10.572406231885966
k=9: 10.507628501663131
k=10: 10.455910429510872
k=11: 10.41370221340334
k=12: 10.378620929581876
k=13: 10.349011593187207
k=14: 10.323691295053795
k=15: 10.301793374505857
k=16: 10.282669031598717
k=17: 10.265823265946006
k=18: 10.250872030293372
k=19: 10.237512882117132
k=20: 8.999022729694703
k=21: 8.267420074010055
k=22: 7.497283519619893
k=23: 6.853431593167068
k=24: 5.895962973131352
k=25: 5.867018843533908


To the eye everything is fine up to $$k=20$$, at which point it stops converging to the correct answer (which is $$10$$). What is likely to be the cause of this? There are only three function calls, so obviously the error is creeping in either when I calculate the matrix power, the Frobenius norm, or its root (or some combination of these things).

The Frobenius norm is not an operator norm, it is a norm on the vector space of linear operators/matrices, which is not the same thing. Just change it to any other preset norm and it should work.

It is also the case that your method of computing matrix powers is not stable. The algorithm used in Numpy is basic repeated squaring, which has no normalization or regularization steps and can blow up for large $$k$$. A more stable way to compute high matrix powers would be to use an eigendecomposition, but you can compute the spectral radius directly from that so it seems moot. Can you more accurately describe your objectives and limitations? Are you doing this for school? Are your matrices prohibitively large? etc.

• Thank you for your answer. I have the habit of loosely conflating operators and their matrix representations. I think your comment that numpy’s matrix power function is unstable for large k is what I was looking for. – Martin C. Dec 29 '19 at 6:36
• Do you know of any python libraries that would be suitable? Alternatively, the name of a suitable algorithm would let me implement it myself. – Martin C. Dec 29 '19 at 6:41
• It turns out the issue was trivial (see my answer). You might want to edit your answer. I still don't know of an appropriate way of doing this for large k... – Martin C. Dec 29 '19 at 18:21

It turns out the issue was very simple (and not related to a limitation of numpy's matrix power function as such). I initially thought there was some numerical floating point error being propagated - but the example matrix I was testing on contains only integers. The problem was that at $$k=20$$ the value of the matrix entries exceeded numpy's maximum possible int64 value - they become of the order of $$10^{20}$$ there.

If I change

a = np.array([[9, -1, 2], [-2, 8, 4], [1, 1, 8]])


to

a = np.array([[9.0, -1, 2], [-2, 8, 4], [1, 1, 8]])


in my code, Gelfand's formula converges as expected as we increase the value of $$k$$.

Some plots to demonstrate the issue. The first (incorrect) one corresponds to the numpy array of type integer above, while the second (correct) one corresponds to the one with a float type.

This resolution is hardly amazing, but maybe posting it here will help someone else who runs into a similar problem. You will still run into numerical issues with this as $$k$$ increases (I start getting inf values at $$k=154$$), so I suppose more sophisticated methods will be necessary if you want to do much better than that.

Code:

import matplotlib.pyplot as plt
import numpy as np
from numpy import linalg as la

def kth_term(matrix, k=20):
"""Gelfand's formula"""
matrix = la.matrix_power(matrix, k)
f_norm = la.norm(matrix, 'fro')
term = f_norm**(1.0/k)

return term

def test_convergence_of_gelfands_formula():
a = np.array([[9.0, -1, 2], [-2, 8, 4], [1, 1, 8]])

for k in range(1, 100):
term = kth_term(a, k)
print("k={}: {}".format(k, term))
plt.plot(k, term, 'rx')

plt.plot(k, term, 'rx', label="kth term")
plt.title("Convergence of Gelfand's formula for the spectral radius")
plt.xlabel("k")
plt.axhline(y=10, color='b', linestyle='-', label='Spectral radius')
plt.legend()
plt.show()

def main():
test_convergence_of_gelfands_formula()

if __name__ == '__main__':
main()

• You could also explicitly set the dtype of the matrix a, then there is no automatic construction as integer matrix. – Lutz Lehmann Dec 31 '19 at 15:09