Iterating the eigenvalue function (eig) produces patterns --- is this a well known result?

Suppose M is a square matrix and [N lambda] = eig(M), so that the columns of N contain the eigenvectors of M and the diagonal matrix elements of lambda are the eigenvalues. If M is unitary, then N is also unitary and the eigvals are phase factors e^(i theta) that lie on the unit circle. << Theoretically each eigvec could have its own norm, but Matlab eig normalizes each of them to 1 which assures that N is unitary. >>

So a unitary matrix M produces another unitary matrix N. Now suppose that this process is iterated, with (going to lower case for the matrices)

[u2 lam1] = eig(u1)
[u3 lam2] = eig(u2)
...
[u(k+1) lam(k)] = eig(u(k))

Since the eigvals lie on the unit circle their only parameter is angle. It turns out that plotting those angles by iteration number produces some interesting repeating patterns.

oo one question is, is this well known behavior?

The function eig( u(k) ) produces u(k+1); the process does not depend explicitly on any eigvals. The eigvals can be seen as just a compact indicator of the states themselves and are much easier to plot since an mxm matrix has m^2 complex elements but only m real eigval angles. << Other plotting possibilities include Re&Im uk(:), Re&Im (trace(uk)) = Re&Im (trace(lamk)), angle (det(uk)) = angle(det(lamk)). >>

The plot below shows an example of eigvals vs iteration when u is 2x2 and the starting matrix u1 is random complex unitary. Since each u has two eigvals, there are two horizontally spaced dots on the plot for each iteration. On the upper plot, eig spends some time making up its mind before settling into a pattern of period 14. The lower plot shows a basic pattern of period 22.

I don’t know how much the details of settling into a pattern depend on the level of precision, double in this case.

At least for matrix sizes 2x2 through 4x4, every initial condition leads to a repeating pattern, although it appears that the number of different-looking patterns is limited. Here are periods for repeating patterns, in decreasing probability for random starting matrices:

2x2 14, 22, 4, 2, 1, 6 3x3 106, 4 ,3 4x4 1966, 898, 256, 154, 84 5x5 622956, 62783

I certainly don’t claim to have found them all. Running a search on a given size of matrix for more than a few minutes takes more patience than I have.

The next plot has patterns for a 3x3 and a 4x4. Most random 3x3 starting matrices produce the interesting pattern with period 106. As the size of the matrix increases, the plots are visually denser. Since an mxm matrix with its m eigvals has m dots per iteration, the actual density per iteration goes up only in proportion to m. Most of the visual density increase is because as m increases, the periods become much longer and many more iterations must be fitted into the plot to show the pattern. As the size proceeds from 2x2 to 4x4, the eigvals concentrate more and more into the small angle regions, producing a dark band. They also concentrate around +-180 degrees, although to a lesser extent. I don’t know the reasons. This is a bit reminiscent of the Wigner semicircle distribution for random matrices, although angles are being plotted and not real values.

For all sizes of matrix, small periods tend to produce stripes and zippers such as on the previous plot.

The next plot shows an iteration period of 1966. To reduce the visual density, every other iteration is plotted. The human eye is pretty good at discerning patterns and one can see repeated dark striations and white blobs, as with the 4x4 in the previous plot.

The 5x5 periods are based on just a few examples. The patterns are too large to plot, and it takes trials with 5e6 iterations to have a good chance of finding one. Even with that many iterations, many random initial matrices don’t find a period and I don’t know if every initial matrix leads to a pattern.

A period of 622956 might seem unlikely but it was verified numerically. This was done simply by, for a given starting matrix, taking the stored matrix A of eigvals for every iteration (double precision, not approximated with bins) and comparing two submatrices: the first consisting of the last 622956 iterations (rows) and the second whose row indices are translated back by the assumed 622956. The largest difference in the 622956 x 5 = 3114780 elements of each submatrix was 4.5e-9. << In order that the period be fundamental, one also has to show that possible shorter periods of commensurate length are not periods. >>

Details on eig conventions

I used straight Matlab eig, but many variations on the output of eig can still provide a valid eig decomposition. Here

[u2 lam1] = eig(u1)

will be used as an example. From the equivalent expression

u1*u2 = u2*lam1

it’s easy to show that for u2 --> u2x and with lam1 unchanged, eig is free to choose one or both of two options:

[1] Multiply each eigvec by a separate arbitrary phase factor. << equivalent to u2x = u2*g where g is a diagonal matrix of phase factors. >> This means that for each eigvec, any nonzero element can be chosen to be positive real.

[2] Reorder the eigvals and permute the eigvecs (columns of u2) to match. << equivalent to u2x = u2*P and lam1x = P’*lam1*P, where P is a permutation matrix. >> For this Question, reordering is considered to be the same solution, not a different one.

So u2 --> u2x preserves the eig decomposition. But now consider the next iteration

[u3 lam2] = eig(u2)           vs.          [u3x lam2x] = eig(u2x)

In general, u3 and u3x are totally different from each other with no obvious relationship, and the same for lam2 and lam2x. Since the eigvals differ, there is no process similar to [1] and [2] that will make the two cases equivalent. Consequently the iteration path depends very much on how eig makes its choices, but not every property is tied down, exactly.

For complex matrices Matlab eig does in fact specify [1]. For each column, eig sensibly sets the element of max abs value to be positive real. << One also could have, for example, first row pos real; last row pos real; diagonal elements pos real; etc. >>

For [2], Matlab eig does not sort the columns in an obvious way. Without knowing the rule one can’t pin down the iteration process completely.

<< Some possible sorts are: sort by the largest abs value in each column; sort columns so that (by abs value) each diagonal element is larger than any element in the row to its right; or just sort by the angles of the eigvals. But since the eigvals lie on a circle, picking a starting point for the sort is pretty arbitrary. >>

oo second question: How does Matlab eig sort its eigvec columns, or does it depend on the input matrix (as with, e.g. Gaussian elimination with pivoting) so that there is no easy way to say?

Trying a few different choices for [1] and [2] did not seem to make a major difference in the allowable patterns. In the end I just used Matlab eig as is.

% script called eigit
tic
% -------- set parameters --------
rng(1234)
m = 3                    % matrix is mxm
nit = 1000;              % number of iterations
ntrial = 1000;           % number of trials, each w/ random start matrix
nd = 360;                % number of discrete digitization bins
% nd needs to be even.
u0 = zeros(m,m,ntrial);         % initial matrix for each trial
% U = zeros(nit,m^2,ntrial);      % eigvec matrices, saved as u(:)
A = zeros(nit,m,ntrial);        % eigval angles
Frs = zeros(nit,ntrial);
per = zeros(1,ntrial);          % period for each trial
perchk = zeros(1,ntrial);       % checks the period
% ----- u0 starting matrices -----
for q = 1:ntrial
    r = randn(m,m) + i*randn(m,m);
    [u s v] = svd(r);
    u0(:,:,q) = u;
end
% ------- results for each trial ----------
h = waitbar(0,'trials');
for q = 1:ntrial  
    u = u0(:,:,q);
    C = zeros(nit,nd);  % C = 1 if one or more values fall into a given bin
    % ------ iterate for qth trial -------
    for k = 1:nit
        [u lam] = eig(u,'vector');
        anorm = (angle(lam)).'/pi;                        % -1 <= anorm <= 1 
        A(k,:,q) = 180*anorm;
        % see binning note
        % U(k,:,q) = u(:);
        ind = floor(rem((nd/2)*anorm +(nd+1)/2,nd))+1;    % 1 <= ind <= nd
        C(k,ind) = 1;
    end
    % ----- find the fundamental period for qth trial --------
    ind = sum(C)>2;
    C = C(:,ind);
    temp = fft(C);
    F = (real(ifft(temp.*conj(temp))));  % autocorrelation of C, down cols
    temp = sum(F,2);                     % row sums
    Frs(:,q) = temp;                     % can plot row sums to troubleshoot
    temp  = temp(2:round(end/2));
    [~, indP] = max(temp);               % period; method not foolproof, not
                                         % great at finding periods 1 and 2
    chk = eigitchk(A(:,:,q),indP);       % verify
    if chk == 4
        perchk(q) = 4;
        per(q) = indP;
    elseif eigitchk(A(:,:,q),2) == 4;    % check for period 2
        perchk(q) = 4;
        per(q) = 2;
    elseif eigitchk(A(:,:,q),1) == 4;    % check for period 1
        perchk(q) = 4;
        per(q) = 1;
    else
        perchk(q) = chk;
        per(q) = indP;
    end
    waitbar(q/ntrial)
end
delete(h)
toc
% ----------- summarize results for a given m ---------
fper = per(perchk == 4);           %  fundamental periods   
fpers = unique(fper);              %  rather use uniq
nfpers = histcounts(fper,[fpers inf]);
eigitsummary = [fpers;nfpers];       
disp(['top row is fundamental periods, bot row is number of occurrences'])
disp(eigitsummary)
% ========= end of calcs  =========

% -------- sample plot -------
if m ==3
    trialq = min(find(fper==106));
    figure(1)
    eigitplot(A(:,:,trialq))
    title('3x3    period = 106')
    xlabel('eigenvalue angle')
    ylabel('iteration number')
    grid on
end
% typical prameters:  m = (2,3,4): rng(1234) (all); nd = 360 (all);
% nit = (1e3,1e3,12e3); ntrial = (1e3,4e3,12e3); // m = 5: rng(11112);
% nd = 100; nit = 5e6; ntrial = 2; produces periods 62783 and 622956
% binning note: for nd = 360 bins, code puts -179 degrees at center of
% first bin, +-180 degrees at center of last bin.  This avoids deceptive
% binning when the angle bounces from 180-const*eps to -180+const*eps
% due to numerical imprecision.  Works since angle is cyclical mod 2 pi.

% ---------
function result = eigitchk(X,nn)
% check to see if matrix X of size (nit x m) i.e. (iterations x values)
% has a presumed fundamental period nn. result: 
% 4: nn is a fundamental period    3: nn is a period, not fundamental
% 2: nn is not a period            1: not possible

a = max(max(abs(X(end-nn+1:end,:) - X(end-2*nn+1:end-nn,:)))); 
nn1 = nn./unique(factor(nn));
b = [];
for k = nn1;
  mx = max(max(abs(X(end-k+1:end,:) - X(end-2*k+1:end-k,:)))); 
  b = [b mx];
  minb = min(b);
end
del = [a;minb];
if nn==1 & all(del<1e-5)    % 1e-5 is boundary between tiny and mid
    result = 4;
else
    result =  2*(a<=1e-5) + (minb>=1e-5) + 1;
end
% Subtracts a sample of nn rows from another sample whose row index
% is translated down by nn.  If a pattern does (does not) repeat
% with period nn, difference is tiny (mid).
% First element of del checks for patterns of period nn.  Second element
% checks for periods smaller than, and commensurate with, nn.
% del: [tiny;mid]-->4  [tiny;tiny]-->3  [mid;mid]-->2  [mid;tiny]-->1
end

%----------
function eigitplot(X,cin)
% plot a matrix X of size (nit x m) i.e. (iterations x values).
% the function plots every cth iteration (reduced dot density).
% input c is optional, default is 1
%
% eigitplot(X,c)

if nargin==2
  c = cin;
else
    c = 1;
end
[nrow ncol] = size(X);
k = (1:nrow)';
kc = k(1:c:end);
Xc = X(1:c:end,:);
q = size(Xc);
% creating columns with iteration number unsorted and
% dependent variable to match is all right for a 'dot' plot
Xplot = Xc(:);
kplot = repmat(kc,ncol,1);
plot(Xplot,kplot,'.','markersize',3)
end

Answer 1

[Christine Tobler]

Hi David,

This is very interesting! To answer your questions:

oo one question is, is this well known behavior?

I'm not aware of this behavior, or of this iteration itself.

oo second question: How does Matlab eig sort its eigvec columns, or does it depend on the input matrix (as with, e.g. Gaussian elimination with pivoting) so that there is no easy way to say?

MATLAB returns the eigenvalues and eigenvectors in arbitrary order, defined by the implementation. That is, different devices or releases can return different order. This is because every possible choice for how to sort eigenvalues on the complex plane will have cases where a small perturbation will have the order change drastically.

This makes it very interesting that you're finding a period. Possibly with small matrices like 2x2 (where there are just two options for ordering the columns), the iteration is jumping back and forth between both options for reordering?

Since these are unitary matrices, you could order the columns coming out of eig based on the angle of the eigenvectors, and pass this reordered matrix back to eig. It could be interesting to see what periodicities this results in.

To understand more of what's going on, one could try if this also works for just a real 2x2 rotation matrix as initial input. One could try to interpret the next matrices as rotation matrices (possibly would need some reflections in there?) and to come up with a formula for the angle of rotation. I would guess that this iteration must then have some attractive periodic points.

There's no clear pattern in the angles from the first plot (at least to me), so this might be quite a complicated formula.

[David Goodmanson]

Hi Christine, thanks for taking a look. I'll try and investigate some more possibilities.

At least for 2x2 I do know of two types of exact solution, both of which are found by iteration if you take enough random trials. One is simply

u1 = u2 = lam1 = lam2 = eye(2) 

which produces a vertical line of (double) dots with iteration spacing 1 at eigval angle = 0. Not very interesting, but it is found.

A few other plots show two vertical lines of dots with iteration spacing 1, at +-45 degrees. In this case the eigvals swap places with each iteration as you supposed, with eigval period 1 since swapped eigvals are considered to be the same solution. On a dot plot the swapped and not-swapped eigvals look the same.

In this case

[u2 lam1] = eig(u1) 
[u1 lam2] = eig(u2) 

with
u1 = [1   g;   -g'  1]/sqrt(2)
u2 = [1 -i*g; -i*g' 1]/sqrt(2)
lam1 = [a' 0; 0 a]     
lam2 = [a 0; 0 a']            % swapped
a =  e^( i*pi/4)
a' = e^(-i*pi/4)
g = e^(i*theta)  

Here g is an arbitrary phase factor, and for the trials that settle into this solution, each one comes up with a different g. I don't know if a solution of this kind exists for eigval angles other than +-45.

Answer 2

[John D'Errico]

I'd have a hard time believing this is well known. @Christine Tobler would be the one to know though, and she agrees.

It just feels like some strange interaction in the eig code, maybe having to do with double precision. So it might be interesting to see if a higher precision eig (so the syms version) saw the same behavior. Periodicity tells me there is something finite happening, but my expectation would be there is no reason for choosing specific sets of eigenvalues on the unit circle. Another possibility is if there is some sort of high order interaction with rand. We do know that in the past random number generators did generate periodicities when things were viewed in higher dimensions. As such, looking at other random generators might prove interesting.

If I were to look more deeply, I would consider that the eigenvalues are the roots of a low order polynomial, just quadratic for the 2x2. So again, if you were to do the same effective computations using syms in high precision, BUT using this time the roots of the characteristic polynomial, would the same behavior emerge? If it does, then you know this is not an artifact that arises from the eig code itself and double precision arithmetic, but is more fundamental.

In the end, my goal would be to isolate where the behavior arises. That might make it easier to understand.

Edit: If you can show this behavior is not eig dependent, nor is it just an artifact of precision in some strange way, then I would want to look at how the periodicity arises. You might be able to do something for the 2x2 case, where you can show just what is "stable" about those periodic sets. But until you manage to show this is not just some artifact of the eig code and its implementation, or something to do with precision, you can do nothing.