Copyright (C) 2022 Andreas Kloeckner
import numpy as np
import numpy.linalg as la
np.set_printoptions(linewidth=150, precision=3, suppress=True)
Consider a random diagonal matrix $A$ with entries sorted by magnitude:
n = 10
eigvals = np.array(sorted(np.random.randn(n), key=abs, reverse=True))
D = np.diag(eigvals)
D
array([[ 1.885, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 1.251, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 0. , -1.215, 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 0. , 0. , -1.19 , 0. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 0. , 0. , 0. , 1.044, 0. , 0. , 0. , 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. , -0.723, 0. , 0. , 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. , 0. , -0.716, 0. , 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , -0.538, 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.314, 0. ], [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , -0.062]])
D @ D @ D
array([[ 6.699, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 1.956, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 0. , -1.792, 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 0. , 0. , -1.687, 0. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 0. , 0. , 0. , 1.139, 0. , 0. , 0. , 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. , -0.378, 0. , 0. , 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. , 0. , -0.367, 0. , 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , -0.156, 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.031, 0. ], [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , -0. ]])
tmp = D @ D @ D
tmp / tmp[0,0]
array([[ 1. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 0.292, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 0. , -0.267, 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 0. , 0. , -0.252, 0. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 0. , 0. , 0. , 0.17 , 0. , 0. , 0. , 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. , -0.056, 0. , 0. , 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. , 0. , -0.055, 0. , 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , -0.023, 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.005, 0. ], [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , -0. ]])
tmp = D @ D @ D @ D @ D @ D @ D @ D @ D @ D @ D @ D
tmp / tmp[0,0]
array([[1. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [0. , 0.007, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [0. , 0. , 0.005, 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [0. , 0. , 0. , 0.004, 0. , 0. , 0. , 0. , 0. , 0. ], [0. , 0. , 0. , 0. , 0.001, 0. , 0. , 0. , 0. , 0. ], [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ]])
This works just as well if the matrix is not diagonalized:
X = np.random.randn(n, n)
Xinv = la.inv(X)
A = X @ D @ Xinv
A
array([[ 1.33 , 0.656, 1.094, 0.629, 1.106, -1.3 , -0.694, 1.579, 1.381, 0.522], [-0.214, 1.982, 1.776, 0.622, 0.847, -2.11 , -0.108, -0.326, 0.916, 1.782], [-0.955, -2.707, -2.992, -1.069, -1.885, 2.917, -0.027, -1.266, -1.929, -1.922], [-0.712, 0.177, 1.006, 0.605, -0.16 , -1.893, -0.055, -2.249, -0.044, 2.541], [ 0.21 , 1.607, 2.268, 1.413, 0.735, -2.988, 0.239, -0.364, 0.901, 2.205], [-0.344, -0.097, -0.005, 0.547, -0.583, -1.129, -0.408, -1.355, 0.138, 1.416], [ 1.047, 0.844, 0.494, 0.546, 0.927, -1.258, -1.565, 0.987, 1.51 , 0.856], [-0.352, -0.499, -0.523, -0.142, -0.536, 0.489, 0.264, -1.036, -0.355, -0.068], [ 0.168, -0.536, -1.716, -0.638, -0.179, 1.161, -0.455, 0.932, 0.048, -0.889], [ 0.595, 1.385, 1.777, 1.45 , 0.671, -2.42 , -0.419, -0.339, 1.312, 2.072]])
A @ A @ A
array([[ 1.557, 2.417, 2.95 , 1.724, 2.136, -3.487, -1.057, 1.349, 2.389, 2.634], [-1.002, 5.675, 6.169, 3.656, 2.382, -7.841, -0.288, -3.128, 2.131, 7.586], [-0.262, -8.049, -9.248, -5.048, -4.445, 10.674, -0.037, 1.463, -3.924, -9.556], [-1.519, 1.341, 2.351, 2.241, -0.068, -4.451, 0.013, -4.634, 0.063, 5.279], [-0.542, 4.897, 6.301, 3.962, 2.238, -7.78 , 0.344, -2.794, 1.922, 7.372], [-0.641, 0.636, 0.824, 1.309, -0.203, -2.521, -0.484, -2.502, 0.337, 2.828], [ 1.087, 3.56 , 3.438, 2.638, 2.337, -5.443, -2.596, -0.385, 3.111, 4.855], [-0.502, -1.019, -1.112, -0.498, -0.836, 1.056, 0.775, -0.723, -0.8 , -0.742], [ 0.454, 0.32 , -0.999, 0.089, 0.306, -0.517, -1.432, 0.261, 1.092, 0.449], [ 0.114, 5.324, 6.541, 4.296, 2.788, -8.307, -0.38 , -2.238, 2.76 , 7.73 ]])
At first, it doesn't look like there is much happening, however:
Apower = A @ A @ A @ A @ A @ A @ A @ A @ A @ A @ A @ A @ A
tmp = Xinv @ Apower @ X
tmp / tmp[0,0]
array([[ 1. , 0. , 0. , -0. , 0. , 0. , 0. , -0. , -0. , 0. ], [ 0. , 0.005, -0. , -0. , -0. , -0. , 0. , 0. , -0. , -0. ], [ 0. , -0. , -0.003, 0. , 0. , 0. , -0. , 0. , -0. , 0. ], [ 0. , -0. , 0. , -0.003, 0. , 0. , 0. , -0. , -0. , -0. ], [ 0. , -0. , -0. , -0. , 0. , 0. , 0. , 0. , -0. , -0. ], [ 0. , -0. , 0. , 0. , 0. , -0. , -0. , -0. , 0. , 0. ], [ 0. , -0. , 0. , 0. , 0. , 0. , -0. , -0. , 0. , 0. ], [ 0. , -0. , -0. , 0. , -0. , 0. , 0. , -0. , 0. , 0. ], [-0. , 0. , -0. , 0. , -0. , -0. , -0. , 0. , 0. , 0. ], [-0. , 0. , -0. , 0. , -0. , -0. , 0. , 0. , 0. , 0. ]])
Let $\boldsymbol{y} = A^{13} \boldsymbol{x}$?
x = np.random.randn(n)
y = Apower @ x
Anything special about $\boldsymbol{y}$?
la.norm(A @ y - D[0,0] * y, 2)/la.norm(y, 2)
0.005644766209795598
Is there a better way to compute $\boldsymbol y$?
y2 = x
for i in range(13):
y2 = A @ y2
la.norm(y2-y)
5.779049855454484e-12