function residual_matrix_method
% residual_matrix_method.m (RMM)
%
% Determination of the eigenvectors for an NxN matrix
% when the N eigenvalues are known. Only eigenvectors for
% unique eigenvalues are determined, but not those for
% repeated eigenvalues. It can give both the column and
% row eigenvectors.
%
% This uses a method known as the 'residual matrix method' (RMM)
% in Unity Root Matrix Theory (URMT)
%
% email: richard@urmt.org
%
% The author published this method (*) in 2010 in the
% free, online PDF paper
%  'Pythagorean Triples as Eigenvectors and Related Invariants'
% available at:
%
%  http://www.urmt.org/pythag_eigenvectors_invariants.pdf
%
% The script is run within Matlab (or the Matlab-compatible 'Octave')
% as per a standard 'm-file' script. When successfully completed,
% it displays the RMM-calculated eigenvectors held in array 'Vmstore',
% followed by those calculated by Matlab, held in array 'eigvecs'.
%
% Note that eigenvectors are only unique to within a scale factor
% and complex eigenvectors may well not appear identical
% between those calculated by RMM and the Matlab-normalised
% calculated eigenvectors (for the same, unique eigenvalue),
% whereby the ratio will invariably be complex, of
% unit norm, given both RMM and Matlab eigenvectors are normalised.
%
% RMM cannot determine the eigenvectors for repeated
% eigenvalues, but only those with unique eigenvalues.
%
% The script is currently setup
% to generate a random 4x4 matrix 'A', hence generally possessing
% four unique eigenvalues from which the eigenvectors are then
% determined using RMM.
%
% There is extra test code for fixed, not random, matrices. This has been
% commented-out, but can, of course, be re-enabled to illustrate other
% cases - some problematic by design, i.e. repeated eigenvalues.
%
% See also Section (4), page 6, Generation of Eigenvectors, which
% solves for three eigenvectors X+, X0 and X- from the
% known eigenvalues C,0,-C respectively for a 3x3 Unity
% Root Matrix A, as detailed in Sections (1) and (2).
%
% (*) Note that the author made no claim of originality
% at the time but neither could an exact reference be found,
% barring a similarity to the numerical Purification
% method, see the aforementioned paper for details. albeit this appears to
% be a rather loose connection.
%
% The author welcomes comments on this and/or any other
% aspects of URMT, email richard@urmt.org
%
% There are no restrictions that the matrix be symmetric
% or Hermitian which is only required for orthogonal
% eigenvectors (to unique eigenvalues).
%
% 24/01/2025 V1.3 R J Miller minor additions and commentary
% 01/10/2019 V1.2 R J Miller sanity check on eigenvalues
% 30/09/2019 V1.1 R J Miller cleanup
% 28/09/2019 V1.0 R J Miller original
%
% Matlab/Octave compatibility
%
% Whilst this script should be compatible with both, some versions
% of Matlab may warn of unmatched 'end' statements, and you may see a few
% of the following warnings:
%
% Warning: File: C:\URMT_documents\archive\urmt_matlab\residual_matrix_method.m Line: 603 Column: 1
% Unmatched "end".
% (Type "warning off MATLAB:m_warning_end_without_block" to suppress this warning.)
%
% The warnings can be safely ignored. The offending 'end' statements are
% for 'end of function' and it would appear there are some compatibility
% issues over the decades of Matlab/Octave when these should be given.
% The warnings come after the main function code and, in fact,
% after the first return statement around line 594 the code could
% safely be deleted if not of interest.
% *****************************************************
clear; % clear workspce variables

errtol = 1e-5; % error tolerance

% Firstly, demonstrate the Residual Matrix Method for an
% order N, random matrix.

N=4; % arbitrary order N matrix, N>=2

% Construct an arbitrary, real-valued random matrix
A = rand(N,N);

% Alternative constructions of A follow, uncomment to suit
% Construct an arbitrary, complex-valued random matrix
% A = rand(N,N) + 1i*rand(N,N);

% construct a random Hermitian matrix
% tmp = rand(N,N); % generate a temporary random matrix
% A = (tmp + tmp')/2; % real-valued symmetric
% tmpi = 1i*rand(N,N); % imaginary random matrix
% A = A + (tmpi + tmpi')/2; % Hemitian

% ***********************************************
% very simple 2x2 demonstrative case
% N=2;
% A = [1 3;2 2]; L1=-1, L2=4

% % info only do not enable these next four lines
% % E1 = [-3 3;2 -2]
% % E2 = [2 3;2 3]
% % V1 = [3;-2] or [2;-2], L1 = -1
% % V2 = [2;2] or [3;3], L2 = 4

% ***********************************************
% another simple 2x2 demonstrative case
% N=2;
% A = [4 1;2 3]; % L1=2, L2=5

% % info only do not enable these next four lines
% % E1 = [-1 1;2 -2]
% % E2 = [4 1;2 3]
% % V1 = [1;-2], row = [1 -1]
% % V2 = [1;1],  row = [2 1]

% ***********************************************

% ***********************************************
% A 3x3 demonstrative failure case,
% repeated eigenvalue L1=5, L2=L3=-2, RMM fails.
% for the repeated -3 eigenvalue, ok for L1=5
% N=3;
% A = [-2 2 -3;2 1 -6;-1 -2 0];

% % info only do not enable these next four lines
% % E1 = 8*[1 2 -3;2 4 -6;-1 -2 3]
% % E2 = E3 = 0 for repeated -3 eigenvalue
% % V1 = [1;2;-1], row = [1 2 -3]
% ***********************************************

% ***********************************************
% Extra matrix tests - use orthogonal eigenvectors
% and spectral decomposition to create test matrices
% for which all the eigenvectors are known (specified)
% and their eigenvector to be determined using RMM.

% N=4; % four orthogonal eigenvectors

X1=[1;1;1;1]/2;
X2=[1;-1;1;-1]/2;
X3=[-1;-1;1;1]/2;
X4=[-1;1;1;-1]/2;

% spectral decomposition - four unique eigenvalues 1,2,3,4
% A = X1*X1' + 2*X2*X2' + 3*X3*X3' + 4*X4*X4';

% two repeated zero eigenvalues for X1 and X2, RMM can
% only determine the eigenvectors X3 and X4 for non-zero
% eigenvalues 3 and 4 only
% A = 3*X3*X3' + 4*X4*X4';

% single zero eigenvalue OK for RMM
% A = 2*X2*X2' + 3*X3*X3' + 4*X4*X4';

% two repeated eigenvalues, '3' - problems with minimal
% polynomial
% A = X1*X1' + 2*X2*X2' + 3*X3*X3' + 3*X4*X4';

% three repeated eigenvalues, '3' - problems with minimal
% polynomial
% A = X1*X1' + 3*X2*X2' + 3*X3*X3' + 3*X4*X4';

% the above test matrices A are real-valued and symmetric
% by construction using X1,X2,X3,X4

% ***********************************************
% Find eigenvectors and eigenvalues independently via
% Matlab - the eigenvalues are required as inputs to
% the Residual Matrix Method

[eigvecs,eigvals] = eigs(A);% Matlab eigenvectors,values

% The eigenvalues must be unique for this method and
% we should really check at this stage that the eigenvalues
% are all unique - but left unchecked here, assume
% uniqueness for now given A is a smallish random matrix
% or purposefully constructd for uniqueness (or
% non-uniqueness to test for errors

I = eye(N,N); % NxN identity matrix

% store of col eigenvectors determined by RMM
Vmstore = zeros(N,N);

% store of col eigenvalues determined by Matlab
Evastore = zeros(1,N);

mV = zeros(1,N); % mth row eigenvector of matrix A

% store of row eigenvectors determined by RMM
mVstore = zeros(N,N);

% Find eigenvector for mth eigenvalue, m=1..N
for mm=1:N

    % Initialise residual matrix for mth eigenvector
    Em = I;
    Lm = eigvals(mm,mm); % mth eigenvalue
    Evastore(mm) = Lm; % store eigenvalue for later

    % form residual matrix from N-1 factors of
    % Cayley-Hamilton matrix polynomial
    for nn=1:N

        % exclude mth eigenvalue from residual matrix,
        % hence there are only N-1 matrix factors in
        % the residual matrix
        if nn ~= mm
            Em=Em*(A - eigvals(nn,nn)*I);
        end
    end

    % Check to see if residual matrix polynomial is not zero
    % which occurs for repeated eigenvalues
    if is_zero(Em, errtol)
        tmpstr = sprintf(...
        'Repeated eigenvalues, zero residual matrix m=%d',mm);
        warning(tmpstr);
        continue; % skip rest of loop - gives a zero eigenvector
    end

    % Each column of Em is a multiple of Vm, i.e.
    % each column of Em is an eigenvector Vm of A
    % scaled by a constant Ym where the Ym are generally
    % unique for each column, but not necessarily so,
    % and may be zero hence giving an all-zero in Em
    % and thus an invalid, all zero eigenvector.
    % However, there will always be
    % at least one, non-zero column (*) and so the mth
    % eigenvector can always be found.

    % (*) The Ym actually form a row eigenvector
    % and neither are they zero.
    % The Ym can be found by the entries in one column by
    % another. For example, providing both the first and mth
    % column of Em, i.e. Em(:,1) and Em(:,m) are both not
    % zero then Em(:,1)./Em(:,m) is a constant vector
    % where every element is the same non-zero constant 'Ym'
    % as in Em(:,1)./Em(:,m) = Ym*ones(N,1).

    for kk=1:N

        % Vm = mth eigenvector from kth column of residual matrix
        Vm = Em(:,kk);

        % Check eigenvector equation
        % This check will fail if Vm is all zero (plausibly
        % possible) or the eigenvalue is repeated and the
        % residual matrix is then zero - this latter case
        % is an error, flagged as warning further above
        tmp = A*Vm - Lm*Vm; % check == 0
        if ~is_zero( tmp, errtol )
            warning('A Eigenvector error %d',mm);
        end

        % if non-zero magnitude eigenvector, to within tolerance,
        % then store unit vector form for printing at end
        tmp = sqrt(Vm'*Vm); % eigenvector magnitude
        if tmp > errtol
            Vmstore(:,mm) = Vm/tmp;
        end;

    end;

    % At this point the mth column eigenvector has
    % been determined and can now be used with the residual
    % matrix to find the row eigenvectors.
    % The row eigenvectors are known as 'conjugate' or
    % 'reciprocal' eigenvectors in URMT - see the paper,
    % Section (3) 'Divisibilty factors' for the
    % vector components alpha,beta,gamma.
    % The residual matrix Em is an outer product of the
    % column eigenvector 'Vm' with its row eigenvector 'mV' to
    % form the matrix projection operator Vm*mV.
    % This below code uses this fact to determine the row eigenvector
    % from the normalised column eigenvector

    for jj=1:N

        for kk=1:N

            % find first non-zero component of col eigenvector jj
            if abs(Vmstore(kk,mm)) > errtol

            % Note that Vmstore(:,mm) is the mth normalised col
            % eigenvector, and Vmstore(kk,mm) is thus the kth component.
            % Em(kk,jj) is the product of Vm(kk)*mV(jj) where
            % Vmstore(kk,mm) ~ Vm to within normalisation

                mV(jj) = Em(kk,jj)./Vmstore(kk,mm);

                break; % quit loop now, only needs one non-zero component

            end; % if

        end; % kk

    end; % jj

    % normalise row eigenvector to unity
    tmp = sqrt(mV*mV'); % magnitude of row vector
    if tmp > errtol
       mVstore(mm,:) = mV/tmp;
    end

    % Check row eigenvector equations
    % This check will fail if mV is zero (plausibly
    % possible) or the eigenvalue is repeated and the
    % residual matrix is then zero - this latter case
    % is an error, flagged as warning further above
    tmp = mV*A - Lm*mV; % check == 0
    if ~is_zero( tmp, errtol )
        warning('Row eigenvector error %d',mm);
    end

end; % mth eigenvector

% *****************************************************
% print column eigenvectors determined via RMM,
fprintf('\nResidual Matrix Method eigenvectors:\n'); % 01/10/2019
fout = 1;
for mm=1:N

    tmpstr = sprintf('L%d=%0.3f + %0.3fi\nVmstore(%d)=',...
        mm,real(Evastore(mm)),imag(Evastore(mm)),mm);
    print_nx1vector_cmplx( Vmstore(:,mm), tmpstr, fout, N )

    % final check on eigenvectors.
     Vm = Vmstore(:,mm);
     Lm = Evastore(mm);
     if is_zero( Vm, errtol )
         warning('Eigenvector %d is zero, eigenvalue %0.3f'...
             ,mm,Lm); % repeated eigenvalue
     end

     tmp = A*Vm - Lm*Vm; % check == 0
%      if ~is_zero( tmp, errtol )
%          warning('B Eigenvector error %d',mm);
%      end

end;

% for reference, print out the eigenvectors determined via Matlab
fprintf('\nMatlab eigenvectors:\n'); % 01/10/2019

for mm=1:N
    tmpstr = sprintf('L%d=%0.3f + %0.3fi\neigvecs(%d)=',...
        mm,real(eigvals(mm,mm)),imag(eigvals(mm,mm)),mm);
    print_nx1vector_cmplx( eigvecs(:,mm), tmpstr, fout, N )
end;

% 01/10/2019 sanity check on eigenvalues - uncomment to use
% A*Vmstore./Vmstore
% mVstore*A./mVstore

% Symmetric or Hermitian matrices only
% The row eigenvectors are orthogonal to the col
% eigenvectors for different eigenvalues.
% Thus, to within sign,
% the matrix product of the row eigenvectors with the col
% row vectors is the identity, i.e.

% tmp = abs(mVstore*Vmstore) - I;
% if ~is_zero( tmp, errtol )
%     error('row eigenvectors not orthogonal to col eigenvectors');
% end

% check RMM eigenvectors are a multiple of Matlab eigenvectors
% 24/01/2025 check code verified for the default, random matrix A case
% Other cases (currently commented-out) may not agree with Matlab
% because the eigenvalue ordering may not be the same between RMM and Matlab,
% e.g. the first RMM eigenvalue may appear as, say, the third in Matlab and,
% in such an example case, it is the first RMM and third Matlab eigenvector
% that should be compared. So, be wary if using these test cases.
 for mm=1:N
    tmp = Vmstore(:,mm)./eigvecs(:,mm);
    tmp2 = norm(tmp)/sqrt(N) - 1; % normed ratio should be unity
    if ~is_zero( tmp2, errtol )
          warning('Eigenvector RMM/Matlab ratio error eigenvector(%d)',mm);
    end
 end

return; % skip the rest for now

%%

% ***********************************************************
% Now demonstrate the Residual Matrix Method as in the original
% paper, 'Pythagorean Triples as Eigenvectors', Appendix (C),
% see header above.
%
% This is for a 3x3 integer 'Hermitian-like' Unity Root Matrix
% with eigenvalues 1,0,-1 and two eigenvectors that are
% Pythagorean Triples (eigenvalues +/-1), and a third Hyperbolic
% Eigenvector

A = [0 -2 2;2 0 -1;2 -1 0]; % Unity Root Matrix - Pythagorean form

[eigvecs,eigvals] = eigs(A); % eigenvalues C = {-1,0,1}

% Matlab-computed eigenvectors 'eigvecs' not required further
% as they will be calculated.

I = eye(3,3);

% Calculate residual matrix Ep for eigenvalue C=1
Ep = (A - 0*I)*(A - -1*I); % written in full to show eigenvalues
% Ep = A*(A + I); % written as one might normally

C = 1; % eigenvalue = 1

Xp = Ep(:,1); % column 1 eigenvector

% check eigenvector for unity eigenvalue
if ~is_zero( A*Xp - C*Xp, errtol )
    warning('A*Xp - C*Xp eigenvector error %d',1);
end

Xp = Ep(:,2); % column 2 eigenvector

% check eigenvector for unity eigenvalue
if ~is_zero( A*Xp - C*Xp, errtol )
    warning('A*Xp - C*Xp eigenvector error %d',1);
end

Xp = Ep(:,3); % column 3 eigenvector

% check eigenvector for unity eigenvalue
if ~is_zero( A*Xp - C*Xp, errtol )
    warning('A*Xp - C*Xp eigenvector error %d',1);
end

% Calculate residual matrix E0 for eigenvalue C=0
E0 = (A - 1*I)*(A - -1*I); % written in full to show eigenvalues
% E0 = (A - I)*(A + I); % written as one might normally

C = 0; % eigenvalue = 0
X0 = E0(:,1); % column 1 eigenvector for

% check eigenvector X0
if ~is_zero( A*X0 - C*X0, errtol )
    warning('A*X0 - C*X0 eigenvector error %d',0);
end

X0 = E0(:,2); % column 2 eigenvector for eigenvalue C=0

% check eigenvector X0
if ~is_zero( A*X0 - C*X0, errtol )
    warning('A*X0 - C*X0 eigenvector error %d',0);
end

X0 = E0(:,3); % column 3 eigenvector for eigenvalue C=0

% check eigenvector X0
if ~is_zero( A*X0 - C*X0, errtol )
    warning('A*X0 - C*X0 eigenvector error %d',0);
end

% Calculate residual matrix Em for eigenvalue C=-1
Em = (A - 1*I)*(A - 0*I); % written in full to show eigenvalues
% Em = (A + I)*A; % written as one might normally

C = -1; % eigenvalue = -1

Xm = Em(:,1); % column 1 eigenvector for

% check eigenvector Xm
if ~is_zero( A*Xm - C*Xm, errtol )
    warning('A*Xm - C*Xm eigenvector error %d',-1);
end

Xm = Em(:,2); % column 2 eigenvector for eigenvalue C=0

% check eigenvector Xm
if ~is_zero( A*Xm - C*Xm, errtol )
    warning('A*Xm - C*Xm eigenvector error %d',-1);
end

Xm = Em(:,3); % column 3 eigenvector for eigenvalue C=0

% check eigenvector Xm
if ~is_zero( A*Xm - C*Xm, errtol )
    warning('A*Xm - C*Xm eigenvector error %d',-1);
end

% The eigenvectors Xp,X0,Vm are not 'normalised to the DCE' in URMT parlance
% Calculate reciprocal (row) eigenvectors using URMT 'conjugate relations'

T = [1 0 0;0 1 0;0 0 -1]; % Pythagorean T operator

X0 = X0/2; % convert to integer
Xm = Xm/-5; % convert to integer

pX = transpose(T*Xm);
mX = transpose(T*Xp);
oX = transpose(T*X0);

% check reciprocal (row) eigenvectors

% eigenvalue C=+1
if ~is_zero( pX*A - pX, errtol )
    warning('pX*A - pX row eigenvector error %d',1);
end

if ~is_zero( oX*A, errtol )
    warning('oX*A row eigenvector error %d',0);
end

% eigenvalue C=-1
if ~is_zero( mX*A + mX, errtol )
    warning('mX*A + mX row eigenvector error %d',-1);
end

% invariant eigenvector inner product
if ~is_zero( pX*Xp - 2*C^2, errtol )
    warning('pX*Xp inner product error %d',2);
end

% invariant zero eigenvector inner product - the DCE
if ~is_zero( oX*X0 - C^2, errtol )
    warning('oX*X0 inner product error %d',1);
end

% invariant eigenvector inner product
if ~is_zero( mX*Xm - 2*C^2, errtol )
    warning('mX*Xm inner product error %d',2);
end

% print the URM3 eigenvectors

N = 3;
print_nx1vector_cmplx( Xp, 'Xp=', fout, N )
print_nx1vector_cmplx( X0, 'X0=', fout, N )
print_nx1vector_cmplx( Xm, 'Xm=', fout, N )

return; % skip the rest for now

%%
% *****************************************************
% test the identity for Hermitian matrices as in the
% paper  https://arXiv.org/abs/1908.03795
% Authors Denton, Parke, Tao, Zhang
%
% The eigenvalues of M and A do not seem related - TBC?
%
N=4; % arbitrary order N matrix, N>=2
% construct a random Hermitian matrix
% tmp = rand(N,N); % real-valued only
tmp = rand(N,N) + 1i*rand(N,N); % generate a temporary random matrix
A = (tmp + tmp')/2; % real-valued symmetric
[VA, tmp] = eigs(A);
LA = diag(tmp); % eigenvalues of A

Vij = zeros(N,N); % eigenvector matrix as in paper

for jj=1:N

    % Form sub-matrix M = N-1 x N-1 sub-matrix
    if jj==1
        M = A(2:N,2:N);
    else
        % jj > 1
        M = zeros(N-1,N-1);
        M(1:jj-1,1:jj-1) = A(1:jj-1,1:jj-1);
        M(1:jj-1,jj:N-1) = A(1:jj-1,jj+1:N);
        M(jj:N-1,1:jj-1) = A(jj+1:N,1:jj-1);
        M(jj:N-1,jj:N-1) = A(jj+1:N,jj+1:N);
    end

    [VM, tmp] = eigs(M); % eigenvalues of sub-matrix
    LM = diag(tmp);

    for ii=1:N

        % calc denominator as continued product
        Ld = 1;

        for kk=1:N
            % exclude ith eigenvalue from product
            if kk  ~= ii
                Ld = Ld*(LA(ii) - LA(kk));
            end
        end

        % calc numerator as continued product
        Ln = 1;
        for kk=1:N-1
            Ln = Ln*(LA(ii) - LM(kk));
        end

        % calcate square norm eigenvector element
        Vij(ii,jj) = Ln/Ld;

    end; % ii

end; % jj

% VA.*conj(VA) % mag squared of each element of eigenvectors
% transpose(Vij)
% if identity if ok, this should be approx zero < 1e-10
transpose(Vij) - VA.*conj(VA)

end % end of function

% utility functions - not central to the method

function print_nx1vector_cmplx( V, start_msg, fout, nn )
error_tol = 1e-3;

% determine if any element is imaginary
imagf = any( (abs( imag(V(:)) ) > error_tol) > 0 );

fprintf(fout,'%s\n',start_msg);

if imagf > 0
    for ii=1:nn
        if abs( real( V(ii) ) ) > error_tol
            dp = abs( imag(V(ii)) - floor(imag(V(ii))) );
            if dp > error_tol
                fprintf(fout,'   %7.3f ',real(V(ii)) );
            else
                fprintf(fout,'   %7.3f ',real(V(ii)) );
            end
        else
            fprintf(fout,'       0 ');
        end
        if abs( imag( V(ii) ) ) > error_tol
            oneunit = abs( abs(imag(V(ii))) - 1 );
            if oneunit < error_tol
                if imag(V(ii)) > 0
                    fprintf(fout,'+      i,' );
                else
                    fprintf(fout,'+     -i,' );
                end
            else
                dp = abs( imag(V(ii)) - floor(imag(V(ii))) );
                if dp > error_tol
                    fprintf(fout,'+ %5.1fi,',imag(V(ii)) );
                else
                    fprintf(fout,'+ %5.0fi,',imag(V(ii)) );
                end
            end
        else
            fprintf(fout,'        ,');
        end
        fprintf(fout,'\n');

    end% ii

else

    % real-only, no imag component
    for ii=1:nn
        dp = abs( real(V(ii)) - floor(real(V(ii))) );
        if dp > error_tol
            fprintf(fout,'   %7.3f ',real(V(ii)) );
        else
            fprintf(fout,'   %7.3f ',real(V(ii)) );
        end
        fprintf(fout,'\n');
    end% ii

end

end

function rtn = is_zero(Xin, error_tol)
    rtn = 0;
    if norm(Xin) < error_tol
        rtn = -1;
    end
end
