A Blind Permutation Similarity Algorithm
Eric Barszcz
cs.DS math.CO 68R10, 68Q25, 05C50, 05C85 F.2.1; G.2.2





This paper introduces a polynomial blind algorithm that determines
when two square matrices, A and B, are permutation similar.
The shifted and translated matrices (A+β I+γ J)
and (B+β I+γ J) are used to color the vertices
of two square, edge weighted, rook's graphs. Then the orbits are found
by repeated symbolic squaring of the vertex colored and edge weighted
adjacency matrices. Multisets of the diagonal symbols from non-permutation
similar matrices are distinct within a few iterations, typically four
or less.


A Blind Permutation Similarity Algorithm
    Eric Barszcz

NASA Ames Research Center

Moffett Field, CA 94035


    December 30, 2023
======================================================================





§ INTRODUCTION


This paper introduces a polynomial blind algorithm for permutation
similarity. Specht's Theorem gives an infinite set of necessary and
sufficient conditions for two square matrices to be unitarily similar
<cit.>. Specht's Theorem compares the
traces of certain matrix products. If the traces ever fail to match,
then the matrices are not unitarily similar. Subsequent work reduced
the number of traces required to a finite number <cit.>.
The permutation similarity algorithm is somewhat analogous, an infinite
set of matrix products are checked to see if the multisets (not traces)
of the diagonal elements match. If they ever fail to match, the two
matrices are not permutation similar.

Permutation similarity plays a large role in the graph isomorphism
problem where two graphs are isomorphic iff their adjacency matrices
are permutation similar. The best theoretical results for graph isomorphism
are by Babai, showing graph isomorphism is quasipolynomial <cit.>.
The typical graph isomorphism algorithm uses equitable partitions
and vertex individualization with judicious pruning of the search
tree to generate canonical orderings that are compared to determine
if two graphs are isomorphic <cit.>. However, Neuen and
Schweitzer suggest that all vertex individualization based refinement
algorithms have an exponential lower bound <cit.>. The
permutation similarity algorithm does not perform vertex individualization
nor does it try to construct a canonical ordering.

The overall process is as follows. Square (real or complex) matrices
A and B are converted into positive integer matrices whose diagonal
entries that are distinct from the off-diagonal entries, (A_Σ+β I+γ J)
and (B_Σ+β I+γ J). The shifted and translated
matrices are used to color the vertices of edge weighted rook's graphs.
The edge weights are `1' for a column edge and `2' for a row edge.
The resultant graphs are called permutation constraint graphs (PCGs),
see Figure <ref>. The purpose of a PCG is to add symmetric
permutation constraints to the original matrix. Then the vertex colored
and edged weighted adjacency matrices of the PCGs are constructed.
Such an adjacency matrix is called a permutation constraint matrix
(PCM). It will be shown that A and B are permutation similar
iff their associated PCMs are permutation similar. Next, the PCMs
are repeatedly squared using symbolic matrix multiplication. Symbolic
squaring generates a canonical string for each inner product. Symbols
are substituted for the strings in a consistent fashion. After the
symbol substitution, the multisets of diagonal symbols are compared.
If they every differ, the matrices are not permutation similar. Symbolic
squaring monotonically refines the partitions until a stable partition
is reached. 

Experimentally, a stable partition is reached in six iterations or
less and non permutation similar matrices are separated in four iterations
or less. The canonical inner product strings generated by symbolic
squaring separate the orbits imposed by the automorphism group acting
on the set of (i,j) locations.

The necessity of the blind algorithm is proven in Section <ref>.
Sufficiency is argued in Section <ref>. Section <ref>
provides a polynomial algorithm that uses the blind algorithm to find
a permutation between permutation similar matrices. Section <ref>
provides some additional material that readers may find interesting.
Section <ref> is a summary.



§ BLIND P-SIMILARITY ALGORITHM


Algorithm <ref> contains pseudocode
for the blind permutation similarity algorithm. This section is organized
so the background for each function is discussed in order. The discussion
will make it clear that the algorithm tests necessary conditions for
two matrices to be permutation similar. A conservative complexity
bound for Algorithm <ref> is presented
in Section <ref>. Sufficiency arguments
are given in Section <ref>.






 §.§ Permutation Similarity


Real or complex square m× m matrices A and B are permutation
similar iff there exists a permutation matrix P such that PAP^T=B
<cit.>. The notation AP∼B
is used to indicate that A and B are permutation similar. Permutation
similarity will be abbreviated as p-similarity in the text.



 §.§ Symmetric Permutation and Mixes


The process of applying a permutation matrix P from the left and
right as PAP^T is called a symmetric permutation <cit.>.
Symmetric permutations have four main characteristics:


  * Elements of a matrix are moved around by symmetric permutation but
not changed, so the multiset of elements in a matrix before and after
a symmetric permutation are identical. Let the multiset of all elements
in a matrix be called the mix, i.e., mix(A)={{ A_i,j}}.
Then mix(A)=mix(PAP^T).

  * Elements of a row are moved together, so the multiset of elements
from a row before and after being moved by a symmetric permutation
are identical. Let the multiset of row multisets be called the row
mix, i.e., rowMix(A)={{[ {{ A_1,:}} ,           ⋯,   {{ A_m,:}} ]}}. Then rowMix(A)=rowMix(PAP^T).

  * Elements of a column are moved together, so the multiset of elements
from a column before and after being moved by a symmetric permutation
are identical. Let the multiset of column multisets be called the
column mix, i.e., colMix(A)={{[ {{ A_:,1}} ,           ⋯,   {{ A_:,m}} ]}}. Then colMix(A)=colMix(PAP^T).

  * Elements on the diagonal remain on the diagonal, so the multiset of
elements on the diagonal is the same before and after a symmetric
permutation. Let the multiset of diagonal elements be called the diagonal
mix, i.e., diagMix(A)={{ A_j,j}}.
Then diagMix(A)=diagMix(PAP^T).




 §.§ Triplet Notation for a Matrix


Because symmetric permutations move values around without changing
them, the values can be replaced with a different set of symbols without
changing the action. It is useful to view a m× m matrix as
a triplet consisting of i) Π: a partition of the (i,j) locations,
ii) Σ: a set of distinct symbols, and iii) g: a bijective
mapping between symbols and cells of the partition. Each cell of the
partition is associated with a single symbol/value. 

For example, let M be a m× m matrix. Define L as the
set of all (i,j) locations, 

    L={(i,j) | 1≤ i≤ m, 1≤ j≤ m} .

 Then the triplet notation for M is 

    M=(Π,Σ,g)

 where Π is a partition of L where each cell contains all
of the locations associated with the same symbol, Σ is the
set of distinct symbols in M, and g is a bijective mapping of
symbols in Σ to cells in Π.

Often Π and Σ are assumed to be ordered. Π is ordered
by first applying lexicographic ordering within each cell to order
the locations in the cell. The first location in each ordered cell
is designated the representative location for that cell. Next,
the representative locations are collected and lexicographically ordered
to order the cells themselves. The set of symbols Σ is ordered
such that the first symbol is associated with the first cell of Π
and the second symbol with the second cell, and so on. So the mapping
g from Σ to Π is simply the identity mapping. If the
notation M=(Π,Σ) appears, it means that Π
and Σ are ordered and g is the identity mapping.

To compare two partitions, Π1 and Π2, the representative
location for each cell is assigned to all locations within the cell,
essentially using the representative locations as the symbol set.
Then the partition difference, Π1-Π2, is defined as

    (Π1-Π2)_i,j=
    0    (i,j)
    
    1    (i,j)


The partition Π will be referred to as a pattern to evoke
a two dimensional image. Π(A) is a refinement of Π(B),
Π(A)≼Π(B), if each cell of Π(A)
is a subset of a cell of Π(B) <cit.>.



 §.§ Symbol Substitution: Lines 11
& 17


A consequence of symmetric permutations moving values/symbols around
without changing them is that the symbols can be replaced without
changing the p-similarity relation.

Given square matrices A
and B. AP∼B A_ΣP∼B_Σ
where A=(Π_A,Σ,g_A), B=(Π_B,Σ,g_B),
A_Σ=(Π_A,Σ,g_A),
B_Σ̅=(Π_B,Σ,g_B) and
|Σ|=|Σ|.



Assume AP∼B so there exists a permutation matrix
P such that PAP^T=B. Further assume that the symmetric permutation
moves the symbol at location (i,j) of A to location (r,s)
in B. So A_i,j=B_r,s=σ_k where σ_k∈Σ
is the k^ th symbol of Σ. Everywhere σ_k
appears in A and B will be replaced by σ_k,
the k^ th symbol of Σ in A_Σ
and B_Σ. Similarly for each of the other symbols
in Σ. Therefore, PA_ΣP^T=B_Σ
so A_ΣP∼B_Σ.

For the other direction, reverse the roles of Σ and Σ.

A consequence of Theorem <ref> is
that the p-similarity problem for real or complex matrices reduces
to solving p-similarity for positive integer matrices. 

The symbol substitution function, SymbolSubstitution( ),
appears on lines 11 and 17 in Algorithm <ref>.
The function takes two square arrays of values (real or complex numbers
in line 11, and strings in line 17) and outputs two square, positive
integer matrices. The substitutions are performed so that distinct
input values are assigned distinct output symbols and identical input
values are assigned identical output symbols.

Given m× m arrays as input, a complexity bound on symbol substitution
is based on assumption that it is performed similar to the Sieve of
Eratosthenes. Arrays are searched in column major order. When a previously
unseen input value is encountered, a new output symbol is assigned
to that input value and O(m^2) comparisons are performed
on each input array to find all locations containing that input value.
Once all matching input values have been replaced with the output
symbol, the next unseen input value is located and the process repeats.
Since there are at most 2m^2 distinct input values, there are
O(m^4) comparisons. So a symbol substitution using
this approach requires O(m^4) comparison operations
times the cost to do a comparison operation. 

The sufficiency argument in Section <ref> uses a second
symbol substitution function SymSub( ) which takes
a single square array of values and output a positive integer matrix.
For reproducibility, symbol assignments are performed in a permutation
independent fashion. Assume there a n distinct symbols where there
are n_1 distinct off-diagonal symbols and n_2 distinct diagonal
symbols, with n=n_1+n_2. The set of n_1 distinct off-diagonal
strings are lexicographically ordered and then 1 is assigned to the
first string in the ordered list and so on until n_1 is assigned
to the last string in the ordered list. Similarly, the list of n_2
distinct diagonal strings are ordered prior to assigning (n_1+1)
to n to them. This guarantees that any symmetric permutation applied
to the array of canonical inner product strings will assign the same
integer to the same string.



 §.§ Shifting and Translating: Lines 13 & 14


The next theorem shows up in spectral graph literature with integer
coefficients for β and γ. It is presented here in the
more general form where the coefficients are real numbers.

Given square matrices A and B.
AP∼B(α A+β I+γ J)P∼(α B+β I+γ J)
where I is the identity matrix, J is the matrix of all ones
and α, β, γ∈ℝ with α≠0.



By noting that PIP^T=I and PJP^T=J for all permutation matrices
P so 

    PAP^T    =     B
    α PAP^T+β I+γ J     =    α B+β I+γ J
    
    P(α A+β I+γ J)P^T    =    α B+β I+γ J.


There is overlap in the capability to replace symbols between Theorem
<ref> and Theorem <ref>.
However, if a symbol appears as a diagonal entry and as an off-diagonal
entry, Theorem <ref> can make the diagonal
entry distinct from the off-diagonal entry by shifting the spectrum,
changing the number of occurrences of that symbol. Theorem <ref>
cannot change the number occurrences of a symbol, it can only replace
a symbol with a new symbol without changing the the number of occurrences. 

The shift and translate function, ShiftAndTranslate( ) appears
in lines 13 and 14 of Algorithm <ref>.
The function takes a square, positive integer matrix M and positive
integer values for β and γ as input and generates the
square, positive integer matrix (M+β I+γ J)
as output.



 §.§ Color Matrix: Lines 11, 13 & 14


Given m× m matrices A and B, define their color
matrices as the result of performing a consistent symbol substitution
followed by shifting and translating so that all symbols are positive
integers greater than or equal to three and the symbols on the diagonal
are distinct from all off-diagonal symbols. In particular, assume
A=(Π_A,Σ,g_A) and B=(Π_B,Σ,g_B)
where |Σ|=k is less than or equal to m^2.
Perform a consistent symbol substitution using Σ={ 1,2,…,k}
to get A_Σ and B_Σ followed
by shifting and translating using β=m^2 and γ=2 to
get color matrices A_C=(A_Σ+m^2I+2J)
and B_C=(B_Σ+m^2I+2J).

The smallest integer in Σ is one and the largest
is k, so shifting the spectrum of A_Σ and
B_Σ by m^2I guarantees that the diagonal
symbols are distinct from off-diagonal symbols. Then adding 2J
guarantees the smallest symbol is greater than or equal to three.

Given square matrices A and B and their
associated color matrices A_C and B_C. AP∼B A_CP∼B_C
.



Result follows from applying Theorem <ref>
to the symbol substitution followed by Theorem <ref>
for the shifting and translating.

Color matrices are formed in Algorithm <ref>
as a two step process. Given square input matrices A and B,
the symbol substitutions are performed on line 11, followed by the
shift and translations on lines 13 and 14 for A_C and B_C
respectively.



 §.§ Permutation Constraint Graph


Given a m× m matrix M. The permutation constraint
graph (PCG) associated with M is the m× m rook's graph
with distinct column and row edge weights and vertex colors matching
M_C, M's color matrix. It is called a permutation constraint
graph because the weighted edges of the rook's graph along with the
vertex coloring add symmetric permutation constraints to M. The
rationale behind PCGs is discussed in Section <ref>.

It is well known that rows and columns of rook's graphs can be permuted
independently and that the rows and columns of a square rook's graph
can be exchanged, so the total number of automorphisms is 2(m!)^2
<cit.>. The automorphism group for a PCG is a subgroup
of the rook's graph automorphism group since all of the edges are
rook's graph edges. PCGs break the symmetry that allows rows to be
exchanged with columns by adding edge weights. A PCG has a column
edge weight of `1' and a row edge weight of `2' as shown in Figure
<ref>. 


To see that a PCG implements all symmetric permutation constraints,
note that rook's graph automorphisms keep rows together and columns
together, similar to symmetric permutations. Symmetric permutation
also requires diagonal elements to remain on the diagonal. PCGs accomplish
this by having the diagonal colors be distinct from off-diagonal colors
as shown in Figure <ref>. A PCG automorphism must permute
a diagonal entry to another diagonal entry as required by symmetric
permutation. This forces the row permutation and the column permutation
to match, leading to the following theorem.

Given square matrices A and B. Their respective
PCGs Γ_A and Γ_B are isomorphic iff AP∼B.



Theorem <ref> establishes that AP∼B
iff A_CP∼B_C so if AP∼B
then Γ_A and Γ_B are isomorphic since their vertex
colors are A_C and B_C are p-similar. In the other direction,
assume PCGs Γ_A and Γ_B are isomorphic. The first
step is to show that the re-labeling of vertices takes the form of
a symmetric permutation. Γ_A and Γ_B have identical
off-diagonal structure. So any re-labeling must act as an automorphism
for that structure, implying the set of possible re-labelings is a
subset of the rook's graph automorphisms that do not exchange the
rows and columns. Further, since the diagonal colors of Γ_A
and Γ_B are distinct from off-diagonal colors, the set of
possible re-labelings is further restricted to only include re-labeling
that apply the same permutation to rows and columns. Therefore, any
re-labeling between Γ_A and Γ_B is applied as
a symmetric permutation. Otherwise the off-diagonal structure is not
maintained or diagonal vertices of Γ_A are not mapped to
diagonal vertices of Γ_B. By hypothesis, Γ_A and
Γ_B are isomorphic, so there exists a symmetric permutation
mapping vertices of , Γ_A to vertices of Γ_B implying
the color matrices A_C and B_C are p-similar. Then by Theorem
<ref> A and B are p-similar.

Permutation constraint graphs are not explicitly formed in Algorithm
<ref>. However, they provide the
link to the permutation constraint matrices described in the next
section.



 §.§ Permutation Constraint Matrix: Lines 13 & 14


The vertex colored adjacency matrix associated with a PCG is called
a permutation constraint matrix (PCM). Column major ordering
is used to bijectively map the vertices of the PCG to the diagonal
of the PCM. Let M_C be the m× m color matrix of the PCG.
Then the diagonal of the PCM is given by 

    D=diag(reshape(M_C,m^2,1))

 where reshape(M_C,m^2,1) uses column major ordering
to reshape the m× m color matrix into a m^2×1 vector
and diag( ) converts the m^2×1 vector into
a m^2× m^2 diagonal matrix. Applying column major ordering
to the vertices of the PCG creates a regular off-diagonal structure
in the PCM. Let R be the m^2× m^2 matrix representing
the off-diagonal structure of a PCM. R is given by 

    R=I_m⊗(1*(J_m× m-I_m))+(J_m× m-I_m)⊗(2*I_m)

 where I⊗(1*(J-I)) are the
column edges, (J-I)⊗(2*I)
are the row edges, and ⊗ is the Kronecker product. So the
PCM of a m× m color matrix M_C is written as 

    PCM(M_C)=D+R.

 Note that all m^2× m^2 PCMs have identical off-diagonal
structure. The only difference between two m^2× m^2 PCMs
is the diagonal.

For example, let M=J_3×3, the 3×3 matrix of all
ones. The associated color matrix is M_C=(J+9I+2J)
and the PCG looks exactly like Figure <ref>, where the
white diagonal vertices are associated with `12' and the off-diagonal
vertices with `3'. The PCM generated using column major ordering is

    PCM(J+9I+2J)=[[ 12  1  1  2        2;  1  3  1     2        2;  1  1  3        2        2;  2        3  1  1  2;     2     1 12  1     2;        2  1  1  3        2;  2        2        3  1  1;     2        2     1  3  1;        2        2  1  1 12 ]]

where column edges have a weight of `1', row edges a weight of `2',
and the blank areas are filled with zeros. The horizontal and vertical
lines in (<ref>) are there to emphasize the block structure.

Given m× m matrices A and B and their
associated m^2× m^2 PCMs, S=PCM(A_C) and T=PCM(B_C),
then AP∼B SP∼T.



Using Theorem <ref> and the fact that column major
ordering is a bijective mapping from the m× m arrays of vertices
of the PCGs to the diagonals of the m^2× m^2 PCMs.

If PCMs S and T are p-similar, then the permutation symmetrically
permuting S to T has the form P⊗ P. To see this, note
that the adjacency matrix of an edge weighted rook's graph with a
column edge weight of `1' and a row edge weight of `2' is given by
R in (<ref>). Since rows and columns of R
can be permuted independently, but rows cannot be exchanged with columns,
the automorphism group of R is 

    Aut(R)={ P_c⊗ P_r}

 where P_c and P_r are m× m permutation matrices
applied to the columns and rows respectively. Since the off diagonal
structure of every m^2× m^2 PCM is identical to R,
any permutation between S and T must have the form P_c⊗ P_r.
However, the vertex coloring on the diagonal of the PCG, when mapped
to the diagonal of the PCM, restrict the possible permutations to
those of the form P⊗ P.

Construction of the PCMs from the color matrices occurs in lines 13
and 14 of Algorithm <ref>. The function
PCM( ) takes a m× m color matrix M_C as
input and returns the associated m^2× m^2 PCM defined
by equations (<ref>), (<ref>), and (<ref>)
as output.

To recap, given (real or complex) square matrices A and B, a
symbol substitution is performed to convert all values/symbols to
be positive integers. Next the color matrices A_C and B_C
are constructed by shifting and translating the positive integer matrices
so the diagonal symbols differ from the off-diagonal symbols and the
smallest value is greater than or equal to three. The color matrices
are used to color the vertices of edge weighted rook's graphs where
column edges have weight one and row edges weight two creating PCGs.
Then PCMs S and T are constructed by using column major ordering
to map the PCG vertices to the diagonals. Theorem <ref>
shows that S and T are p-similar iff A and B are p-similar.




 §.§ Symbolic Matrix Multiplication




  §.§.§ Eigenspace Projector Patterns


PCMs are real symmetric matrices. Therefore they have a unique spectral
decomposition in terms of the distinct eigenvalues and their associated
spectral projectors 

    A=∑λ_iE_iI=∑ E_i

 where λ_i is an eigenvalue and E_i is its associated
spectral projector <cit.>. Assume there are k
distinct eigenvalues. Rearrange (<ref>) as

    A_1× m× m=[λ_1,⋯,λ_k]×[[ E_1;   ⋮; E_k ]]_k× m× m

 where the spectral projectors are stacked like a pages in a book
lying on a desk and A is a sheet of paper lying on the desk. Any
symmetric permutation applied to A on the lhs is also applied to
each of the spectral projectors on the rhs. Therefore we can construct
a string from each (:,i,j) column by concatenating the values.
Then substituting distinct symbols for distinct strings yields a pattern.
We call this pattern the eigenspace projector pattern.

To find the eigenspace projector pattern you need to compute the eigenvalues
and eigenvectors. This is an iterative process and may be susceptible
to floating point arithmetic errors. Alternatively, the eigenspace
projector pattern can be generated from a stack of powers of A
by taking advantage of Wilkinson's observation <cit.>
that the transpose of the Vandermonde matrix 

    V^T=[[       1       ⋯       1;     λ_1       ⋯     λ_k;       ⋮       ⋯       ⋮; λ_1^k-1       ⋯ λ_k^k-1 ]],

 is non-singular. Powers of A can be stacked on the lhs as

    [[     I;     A;     ⋮; A^k-1 ]]_k× m× m=[[       1       ⋯       1;     λ_1       ⋯     λ_k;       ⋮       ⋯       ⋮; λ_1^k-1       ⋯ λ_k^k-1 ]]×[[ E_1;   ⋮; E_k ]]_k× m× m.

Examining (<ref>), one sees that strings constructed
on the lhs of (<ref>) yield a pattern identical
to the eigenspace projector pattern since V^T is non-singular.
PCMs are symmetric integer matrices so integer arithmetic can be used
to construct the eigenspace projector pattern. However, overflow errors
may occur.

The eigenspace projector pattern is a refinement of the individual
spectral projectors. So given a real symmetric matrix, the eigenspace
projector pattern is the most refined pattern one would expect for
that matrix. In graph isomorphism terms, the eigenspace projector
pattern is like the coarsest equitable partition <cit.> achievable
using naive color refinement. For two real symmetric matrices to be
p-similar, their eigenspace projector patterns must be p-similar.



  §.§.§ Computing Eigenspace Projector Patterns for SPD Matrices


For symmetric positive definite (SPD) matrices, one doesn't need strings
constructed from the stack of powers of A, lhs of (<ref>),
to determine the eigenspace projector pattern. For a high enough power,
say n, the pattern Π(A^n) is identical to the
eigenspace projector pattern.

For real, symmetric positive definite matrix A,
there exists a finite integer n such that the pattern in A^n
is identical to the eigenspace projector pattern, and all higher powers,
A^n+l, l=1,… have the same pattern.



Assume symmetric positive definite matrix A has k distinct eigenvalues
and they are ordered 0<λ_1<λ_2<⋯<λ_k.
We only need to focus on a subset of locations from the lhs and rhs
of (<ref>), one for each distinct string. Let (i,j)
and (r,s) be two off-diagonal locations with distinct strings.
Let eij=E(:,i,j) and ers=E(:,r,s) be the
column vectors associated with the strings on the rhs of (<ref>)
and aij and ars be the column vectors associated with the strings
on the lhs of (<ref>). Note that 

    [[ 1 ⋯ 1 ]]× eij=[[ 1 ⋯ 1 ]]× ers=0

 from (<ref>) since (i,j) and (r,s) are off-diagonal
locations.

Let δ=eij-ers, then δ≠0 since the strings differ
but 

    [[ 1 ⋯ 1 ]]×δ=0.

 Equation <ref> also applies to differences
between distinct strings on the diagonal since 

    [[ 1 ⋯ 1 ]]× eii=[[ 1 ⋯ 1 ]]× ejj=1.

 Comparisons between diagonal and off-diagonals strings are ignored
since they can never match.

Given that eij≠ ers, then aij≠ ars and there is at least
one power t, 1≤ t≤ k-1, such that A_i,j^t≠ A_r,s^t.
The question remains, is there a finite integer n such that A_i,j^n≠ A_r,s^n
and A_i,j^n+l≠ A_r,s^n+l for all l=1,....

Assume for every sequence of consecutive powers t where A_i,j^t-A_r,s^t≠0,
there is a maximum integer n such that A_i,j^n-A_r,s^n≠0
but A_i,j^n+1-A_r,s^n+1=0. Now write A_i,j^n+1-A_r,s^n+1=0
as 

    A_i,j^n+1-A_r,s^n+1=[λ_1^n+1,⋯,λ_k^n+1]×δ=[[ 1 ⋯ 1 ]]×[[     λ_1^n+1δ_1;              ⋮; λ_k-1^n+1δ_k-1;     λ_k^n+1δ_k ]]=0.

 Without loss of generality, assume that δ_k≠0, and make
the following substitution 

    [[     λ_1^n+1δ_1;              ⋮; λ_k-1^n+1δ_k-1;     λ_k^n+1δ_k ]]=[[                     λ_1^n+1δ_1;                              ⋮;                 λ_k-1^n+1δ_k-1; -(λ_1^n+1δ_1+⋯+λ_k-1^n+1δ_k-1) ]]

 where the right hand vector uses (<ref>) to
substitute for λ_k^n+1δ_k. Equating the last terms
of (<ref>) yields

    λ_k^n+1δ_k=-(λ_1^n+1δ_1+⋯+λ_k-1^n+1δ_k-1)

or 

    λ_k^n+1=-(λ_1^n+1δ_1+⋯+λ_k-1^n+1δ_k-1)/δ_k.

 Without loss of generality, assume δ_k-1≠0, and factor
λ_k-1^n+1 out of the numerator 

    λ_k^n+1=λ_k-1^n+1(-((λ_1/λ_k-1)^n+1δ_1+⋯+1δ_k-1)/δ_k),

and rearrange to get 

    (λ_k/λ_k-1)^n+1    =    -((λ_1/λ_k-1)^n+1δ_1+⋯+1δ_k-1)/δ_k.

 Taking the limit of both sides of (<ref>) results
in 

    lim_n→∞(λ_k/λ_k-1)^n+1=+∞lim_n→∞(-((λ_1/λ_k-1)^n+1δ_1+⋯+1δ_k-1)/δ_k)=-δ_k-1/δ_k

 creating a contradiction. Therefore, there must be some finite integer
n such that A_i,j^n≠ A_r,s^n and A_i,j^n+l≠ A_r,s^n+l
for all l=1,….

Applying the same argument pairwise to all distinct off-diagonal strings
shows that a n exists for each distinct pair of strings. Taking
the maximum over all of the n's yields the desired result.

Theorem <ref> does not guarantee that n is less than
or equal to k-1, or that the value for k is known. It states
that mathematically the pattern in powers of a real SPD matrix converges
to the eigenspace projector pattern. In practice, computer arithmetic
is finite precision. So if n is large, the emergent pattern will
only reflect the eigenspace projector(s) associated with the dominant
eigenvalue(s).



  §.§.§ Why Symbolic Matrix Multiplication


Theorem <ref> says the eigenspace projector patterns of
integer SPD matrices can be computed using integer arithmetic. However,
it is not known when the eigenspace projector pattern will be reached,
and the computation is susceptible to overflow. When the author attempted
to find eigenspace projector patterns using recursive squaring, overflow
would occur after a few iterations. To avoid overflow, a symbol substitution
is performed to reduce the magnitude of the integers. The process
repeated alternating between recursive squaring and symbol substitution
until either the pattern converged or the number of symbols was so
large that symbol substitution could not reduce the magnitudes enough
to prevent overflow.

Key observation: There exist cases where the symbol substitution
resulted in a matrix whose eigenspace projector pattern is a strict
refinement of the original eigenspace projector pattern. These
“new” matrices have the same p-similarity relationship as the
original matrices but a more refined eigenspace projector pattern.

An example is the first graph of case 10 from the website http://funkybee.narod.ru/graphs.htm.
The graph has 11 vertices and its adjacency matrix is given in Figure
<ref>.

 Its eigenspace projector pattern has 29 cells but recursive squaring
with symbol substitution resulted in a pattern with 30 cells. An exhaustive
search revealed four automorphisms. Applying the automorphisms yielded
30 orbits. The orbits are identical to the pattern generated by the
recursive squaring with symbol substitution. Comparing the orbits
with the eigenspace projector pattern shows that the cell associated
with symbol `6' in the eigenspace projector pattern contains two orbits,
see Figure <ref>.


The question is how to find these more refined eigenspace projector
patterns without resorting to random symbol substitutions. This led
to examining symbolic matrix multiplication as a possible method for
simultaneously checking all possible symbol substitutions while avoiding
overflow.



  §.§.§ Canonical Inner Product Strings


It turns out symbolic matrix multiplication also has an advantage
over regular matrix multiplication in separating orbits. Assume AP∼B
then there exists a P such that PAP^T=B. If P:i→ rj→ s
then A_i,j=B_r,s. Also, 

    PAP^T× PAP^T    =     B× B
    
    P(AP^T× PA)P^T    =     B× B

 implying 

    [AP^T× PA]_i,j=[B× B]_r,s

where 

    A_i,:P^T× PA_:,j=B_r,:× B_:,s.

 Since symmetric application of P moves row i of A to become
row r of B (with the values permuted) and column j of A
becomes column s of B (with the values permuted in the same
way) the multiset of terms in the inner products are identical. Further,
the order of the factors within a term are also identical. To see
why this is interesting, assume the inner products from two locations
are given by 

    [[ α β ]]×[[ β; α ]]=αβ+βα=γ

and 

    [[ α α ]]×[[ β; β ]]=αβ+αβ=γ.

 Using regular matrix multiplication, both inner products appear to
be identical, evaluating to γ. However, symbolically the strings
formed from the vector of terms 

    [[ αβ; βα ]]

are different, implying the two locations are not in the same orbit.
So symbolic inner products can separate orbits that appear to be the
same to regular matrix multiplication.

Define the canonical inner product string of u^T× v
as the ordered multiset of inner product terms where the first factor
in each term is from the row vector u^T and the second from the
column vector v 

    canonicalString(u^T× v)=order({{(u_k,v_k)}}).

To create the canonical inner product string, the multiset of terms
is split into two parts: terms involving diagonal symbols and terms
involving two off-diagonal symbols. For terms involving diagonal symbols,
they are ordered such that the row vector diagonal symbol comes first,
followed by the term involving the column diagonal symbol. For terms
involving two off-diagonal symbols, they are lexicographically ordered.
Then the ordered parts are concatenated, diagonal terms followed by
off-diagonal terms, to construct the canonical inner product string.





  §.§.§  𝐒𝐲𝐦𝐒𝐪𝐫: Line 17


A PCM is a symmetric positive integer matrix whose diagonal symbols
are distinct from off-diagonal symbols. Therefore the output of symbolically
squaring a PCM has a pattern that is symmetric and the diagonal canonical
inner product strings are distinct from off-diagonal canonical inner
product strings. The canonical inner product strings are constructed
as described in Section <ref>. 

SymSqr() is a function that takes a symmetric matrix
whose diagonal symbols are distinct from off-diagonal symbols as input
and generates an array of canonical inner product strings as output.
Since the canonical inner product strings at (i,j) and (j,i)
likely differ, the output array is made symmetric by choosing the
lessor of the canonical inner product strings at (i,j) and (j,i)
to represent both locations. 


Using a Jordan product to construct matching canonical strings for
(i,j) and (j,i) was considered. If a single string is constructed
from the combined terms of the inner products at (i,j) and (j,i),
it is possible that other locations (r,s) and (s,r), with distinct
canonical inner product strings, will generate the same string if
their terms are mixed together. This would hide the fact that the
locations are actually in different orbits. However, if the Jordan
product is interpreted as lexicographically ordering the pair of canonical
strings from (i,j) and (j,i) and concatenating them into a single
string, the final strings will be distinct and equivalent to choosing
the lessor string.

Locations in the same orbit have identical
canonical inner product strings.



By the arguments used in Section <ref> to derive and
construct the canonical inner product strings.

The contrapositive of Theorem <ref> causes patterns
to be refined. If two locations in the same cell of Π(M)
have distinct canonical strings when M is symbolically squared,
they cannot be in the same orbit and will be assigned to different
cells in the new pattern Π(SymSqr(M)).

Locations whose canonical inner product
strings differ are not in the same orbit.



Contrapositive of Theorem <ref>.


If M is a square matrix whose diagonal
symbols are distinct from off-diagonal symbols, then the cells of
Π(SymSqr(M)) represent disjoint sets of
orbits.



This is a direct consequence of Theorem <ref>.
Canonical inner product strings for locations in the same orbit are
identical. Therefore locations in the same orbit are assigned to the
same cell.

Theorem <ref> does not guarantee there is one
orbit per cell. Only that locations in an orbit will not be split
across multiple cells.

The symbolic squaring function, SymSqr(), appears twice
on line 17 of Algorithm <ref>. Each
input is a m^2× m^2 symmetric, positive integer matrix
whose diagonal symbols are distinct from the off-diagonal symbols.
The outputs are m^2× m^2 square arrays of canonical inner
product strings, where each string is composed of m^2 terms.



  §.§.§ Monotonic Refinement


Let M be a square matrix whose diagonal symbols are distinct from
the off-diagonal symbols. The next theorem addresses repeated symbolic
squaring of M. Theorem <ref> shows that
symbolic squaring monotonically refines the pattern since elements
of two different cells will never be assigned to the same cell.

Repeated symbolic squaring of a square
matrix whose diagonal symbols are distinct from off-diagonal symbols
monotonically refines the pattern until a stable pattern is reached.



Let M be a square matrix whose diagonal symbols are distinct from
the off-diagonal symbols. To show that refinement is monotonic, it
is enough to show that two locations with distinct symbols can never
have identical canonical inner product strings. To see this, let (i,j)
and (r,s) be two off-diagonal locations where M_i,j≠ M_r,s.
Without loss of generality assume that i<j and r<s. Let D
be the diagonal of M. Then the inner products [M× M]_i,j
and [M× M]_r,s look like 

    [M× M]_i,j=[[     ⋯ D_i,i     ⋯ M_i,j     ⋯ ]]×[[     ⋮; M_i,j;     ⋮; D_j,j;     ⋮ ]]=[[ 1 ⋯ 1 ]]×[[          ⋮; D_i,iM_i,j;          ⋮; M_i,jD_j,j;          ⋮ ]]

 and 

    [M× M]_r,s=[[     ⋯ D_r,r     ⋯ M_r,s     ⋯ ]]×[[     ⋮; M_r,s;     ⋮; D_s,s;     ⋮ ]]=[[ 1 ⋯ 1 ]]×[[          ⋮; D_r,rM_r,s;          ⋮; M_r,sD_s,s;          ⋮ ]].

 The canonical inner product string for [M× M]_i,j
includes terms { (D_i,i,M_i,j), (M_i,j,D_j,j)}
and the canonical string for [M× M]_r,s includes
terms { (D_r,r,M_r,s), (M_r,s,D_s,s)}. Since
M_i,j≠ M_r,s the strings do not match, and they will be
assigned distinct symbols.

A similar argument holds for diagonal inner product strings where
M_i,i≠ M_r,r. Their canonical inner product strings have
terms (D_i,i,D_i,i) and (D_r,r,D_r,r)
respectively that do not match. Diagonal locations and off-diagonal
locations have a different structure to their canonical inner product
strings. Diagonal locations have a single term involving diagonal
symbols, whereas off-diagonal locations have two terms involving diagonal
symbols. Therefore, symbolic squaring either strictly refines the
pattern or has reached a stable pattern.



If M is a symmetric m× m matrix
whose diagonal symbols are distinct from off-diagonal symbols, then
symbolic squaring followed by symbol substitution converges to a stable
pattern within m iterations.



Given M, a symmetric m× m matrix whose diagonal symbols
are distinct from off-diagonal symbols. Recall that Theorems <ref>
and <ref> imply that any pattern that is symmetric
and has diagonal values distinct from off-diagonal values can be replaced
by a SPD matrix with positive integer values. So, M and the output
of symbolically squaring M can always be made to be SPD matrices
with positive integer values.

Now consider equation (<ref>). The eigenspace projector
pattern on the rhs of (<ref>) can be generated from
the stack of k-1 powers of M on the lhs, where k is the number
of distinct eigenvalues. A m× m matrix has at most m distinct
eigenvalues. So at most m-1 powers of M are needed on the lhs
of (<ref>) to construct the eigenspace projector
pattern.

Next consider replacing successive powers of M by symbolic squaring
with symbol substitution. At most m-2 symbolic squarings are needed
to fill the stack since I and M are the first two layers. Theorem
<ref> guarantees that repeated symbolic
squaring monotonically refinements the pattern. Therefore, each successive
layer is a refinement of all prior layers. 

Theorem <ref> states that the powers of a SPD matrix converge
to the eigenspace projector pattern. The bottom layer of the stack
constructed from repeated symbolic squaring is a refinement of all
prior layers. Therefore the bottom layer is at least as refined as
the eigenspace projector pattern of M.

So symbolic squaring converges to a stable pattern within m iterations.

Symbolic squaring followed by symbol substitution occurs on line 17
of Algorithm <ref>.

Let M be a square matrix whose diagonal symbols differ from off-diagonal
symbols. Define Π(M^*) as the stable pattern reached
by repeated symbolic squaring and symbol substitution.


One might wonder why symbolic squaring does not converge in a single
iteration. The answer appears to be related to the number of initial
symbols and the off-diagonal structure of a PCM. For the first iteration,
a PCM has at least five distinct symbols and the off-diagonal is highly
structured. For the five symbol case, see (<ref>), it can
be shown that after the first iteration, there are 9 distinct symbols.
A second symbolic squaring results in 11 distinct symbols and is the
stable pattern. For another example, the PCM of a Petersen graph has
six symbols. The stable pattern has 65 symbols and is reached on the
third symbolic squaring. So the structure and limited number of initial
symbols typically causes symbolic squaring to take several iterations
to converge. In the authors experience, no case has taken more than
six iterations to converge and the number of distinct symbols can
be in the 10's of millions.



 §.§ Comparing Multisets: Line 18


Theorem <ref> shows that matching canonical inner
product strings is a necessary condition for two locations to be in
the same orbit. This implies the multiset of canonical inner product
strings for SymSqr(A) and SymSqr(B) must
match for A and B to be p-similar, mix(SymSqr(A))=mix(SymSqr(B))[See Section <ref> for the definitions
of the various mixes.]. Conversely, if the mixes do not match, then A and B cannot
be p-similar. However it is possible to detect when the mixes do not
match by looking at a smaller multiset. If the mixes do not match,
mix(SymSqr(A))≠ mix(SymSqr(B)),
then the column mixes do not match, colMix(SymSqr(A))≠ colMix(SymSqr(B)),
which in turn implies the diagonal mixes will not match, diagMix(SymSqr(A))≠ diagMix(SymSqr(B)).
Column mixes may differ one iteration before the diagonal mixes if
the differences are in off-diagonal symbols. However, after the next
iteration, the off-diagonal differences are reflected in the diagonal
symbols causing the diagonal mixes to differ. Comparing diagonal mixes
is performed on line 18 in Algorithm <ref>.

Comparing diagonal mixes to determine when two matrices are not p-similar
is consistent with Specht's Theorem <cit.>.
Specht's Theorem gives necessary and sufficient conditions for showing
when two matrices are unitarily similar. Specht's Theorem works by
comparing the traces of sequences of matrix products of various lengths.
For p-similarity it is not enough to show that the traces match, the
diagonal mixes must match. So, matching diagonal mixes, after every
round of symbolic squaring, is a necessary condition for A and
B to be p-similar. This implies differing diagonal mixes is a sufficient
condition for declaring two matrices non p-similar. Note that symbolic
squaring with symbol substitution means the diagonal mixes of an uncountably
infinite number of matrix pairs are compared each iteration.



 §.§ Comparing Patterns, Line 20


The last function in Algorithm <ref>
compares two patterns and returns TRUE if the patterns are identical.
If the pattern before, Π(M), and after symbolic squaring,
Π(SymSqr(M)), match then the stable pattern
Π(M^*) has been reached for that matrix. Necessary
conditions for A and B to be p-similar are that the diagonal
mixes match and they reach stable patterns on the same iteration.
The first condition is checked on line 18 and the second on line 20
of Algorithm <ref>.

Note that the case where diagonal mixes match and one PCM has reached
a stable pattern but the other has not, will be detected as non p-similar.
The non stable pattern is refined by the symbolic squaring in the
next iteration, changing the column mixes. This in turn implies the
diagonal symbols will change, causing the diagonal mixes to differ.
Ultimately resulting in the PCMs being declared non p-similar.



 §.§ Algorithm <ref>
Complexity


The overall complexity of Algorithm <ref>
for comparing two m× m matrices is O(m^12).
The bound is conservative and computed using simple, unoptimized approaches.
In particular, symbol substitutions are assumed to use the method
described in Section <ref>.

Below are complexity bounds on individual functions of Algorithm <ref>
given two m× m matrices as input:


  Line 11, SymbolSubstitution O(m^4) since there
are 2m^2 locations containing values and O(m^2)
comparison operations per substitution with a constant cost per comparison
operation.

  Lines 13&14, ShiftAndTranslate O(m^2) since
there are 2m^2 locations to update.

  Lines 13&14, PCM O(m^4) since there are 2m^4
locations to initialize.

  Line 17, SymSqr O(m^8) reasoned as follows:


  * Construct terms from factors: O(m^6) since there are
2m^4 canonical strings and each string has m^2 terms and
a constant cost to construct a term from its factors.

  * Construct canonical inner product strings from the terms: O(m^8)
since there are 2m^4 strings and it takes a maximum of O(m^4)
comparison operations to order the m^2 terms of a string at a
constant cost per comparison operation.

  * Select lessor string for (i,j) and (j,i) locations: O(m^6)
since there are 2m^4 locations and comparing two strings takes
O(m^2) comparison operations at a constant cost per
comparison operation.


  Line 17, SymbolSubstitution O(m^10) since there
are 2m^4 locations and O(m^4) string comparison
operations per substitution at a cost of O(m^2) to
compare two strings.

  Line 18, CompareDiagMultisets O(m^4) since
it takes a maximum of O(m^4) comparison operations
to sort a multiset with m^2 values at a constant cost per comparison
operation. Comparing the sorted multisets takes at most O(m^2)
comparison operations.

  Line 20, ComparePatterns O(m^8) since the comparison
is performed by replacing values with their representative locations
and then comparing the results. To replace values with representative
locations involves 2m^4 locations and O(m^4)
comparison operations per substitution at a constant cost per comparison
operation. This is followed by at most O(m^4) comparison
operations to compare the two patterns. 

  Lines 16-24, Repeat Until Iterations O(m^2)
by Theorem <ref> applied to a m^2× m^2
SPD matrix.

  Overall Complexity O(m^12) derived from iterating
the repeat until loop a maximum of O(m^2) times, and
a O(m^10) cost to substitute symbols for canonical
inner product strings each iteration.

As mentioned previously, in the author's experience, no case has taken
more the six iterations to converge. The greatest reduction in overall
complexity will come from optimizing the symbolic substitution applied
to canonical inner product strings. Even better would be to find a
way to avoid symbolic matrix multiplication all together, see Section
<ref>. In terms of space, Algorithm <ref>
requires O(m^6) locations to hold the 2m^4 canonical
inner product strings, where a location is large enough to hold a
term.



§ SUFFICIENCY ARGUMENT


In Section <ref>, Algorithm <ref>
is shown to test necessary conditions for two PCMs to be p-similar
and that differing diagonal mixes are sufficient to determine that
two matrices are not p-similar. The last theoretical hole is to determine
whether the patterns of non p-similar PCMs can stabilize on the same
iteration with matching diagonal mixes. The purpose of this section
is to argue that it is not possible, implying Algorithm <ref>
is necessary and sufficient for determining when two matrices are
p-similar. 

An outline of the argument goes as follows: 


  * Given a m^2× m^2 PCM and its automorphism group, the set
of all symmetric positive definite matrices with the same automorphism
group is constructed.

  * The set is shown to be convex. 

  * Repeated symbolic squaring with symbol substitution is shown to act
as a descent method on the convex set, implying the orbits are found
within m^2 steps. 

  * Having the orbits does not guarantee the symbol to orbit assignments
are consistent across matrices. 

  * Using properties of PCMs, it is shown that an Oracle either matches
all of the columns S^* and T^*, or none of the columns.
This differs from cases where an Oracle is able to partially match
non p-similar matrices, Section <ref>.

  * Next, the PCM of the color matrix for the direct sum of the initial
pair of matrices, M=PCM((A⊕ B)_C),
is constructed and symbolically squared until its orbits, Π(M^*),
are found. PCMs S and T are embedded within PCM M.

  * The patterns at locations associated with S and T in Π(M^*)
must match Π(S^*)and Π(T^*) since
they are the orbits imposed by Aut(S) and Aut(T),
respectively. 

  * Since S and T are not p-similar, the Oracle guaranteed no column
of S^* can be matched with any column of T^*. So orbits
associated with locations of S^* in M^* are distinct from
orbits associated with locations of T^* in M^*. 

  * So the diagonal symbols in M^* associated with locations of S^*
and T^* are distinct. 

  * Therefore, the diagonal symbols of S^* and T^* are distinct
as long as consistent symbol substitution is performed each symbolic
squaring.

  * This implies comparing diagonal mixes of S^* and T^* is
sufficient for detecting non p-similar matrices.




 §.§ Crux of the Argument


Showing that repeated symbolic squaring of a PCM converges to the
orbits is the crux of the argument. Given that it is true, a simplified
argument constructs M=PCM((A⊕ B)_C)
and then either the multisets of diagonal locations associated with
A_C and B_C in M^* match or they differ, implying one
or more of the orbits are distinct. This would work as a blind permutation
similarity algorithm. This follows the graph isomorphism argument
that two graphs are isomorphic iff there exists an automorphism of
their direct sum that exchanges the graphs.

Repeated symbolic squaring of the direct sum (A⊕ B),
without constructing the PCM of the direct sum, is not guaranteed
to separate non p-similar matrices. For example, the non p-similar
adjacency matrices `had-sw-32-1' and `had-sw-32-2' from the Bliss
graph collection <cit.> do not have their orbits separated
by repeated squaring of the direct sum. However, their associated
PCM's S and T have distinct diagonal mixes after the third symbolic
squaring.

Algorithm <ref> uses S and T
instead of the PCM of the direct sum M because S and T are
four times smaller. This makes working with S and T more practical
than working with M. 



 §.§ WLOG Symbolic Squaring Generates
SPD Matrices


As pointed out in Section <ref>, inputs to the symbolic
squaring function SymSqr() are symmetric, positive integer matrices
whose diagonal symbols are distinct from off-diagonal symbols. The
output is a symmetric array of canonical inner product strings whose
diagonal strings are distinct from off-diagonal strings. This is immediately
followed by a symbol substitution where the canonical inner product
strings are replaced by positive integers resulting in an output that
is a symmetric, positive integer matrix whose diagonal symbols are
distinct from off-diagonal symbols. So, one is free to choose the
symbol set, as long as they are positive integers. 

In particular, assume there are n distinct canonical inner product
strings, of which there are n_1 distinct off-diagonal strings
and n_2 distinct diagonal strings, with n_1+n_2=n. Choose
the consecutive integers 1 to n_1 as the set of symbols to
replace off-diagonal strings and consecutive integers (m^2n_1+1)
to (m^2n_1+n_2) as the set of symbols to replace
diagonal strings. Then the resulting matrix will be symmetric and
diagonally dominant by construction, and so symmetric positive definite
(SPD). Therefore in the sufficiency arguments below, symbolic squaring
followed by symbol substitution generates SPD matrices unless otherwise
noted.



 §.§ SymSub


Algorithm <ref> uses SymbolSubstitution
to perform a consistent symbol substitution on two arrays of canonical
inner product strings. In this section, symbol substitution is often
performed on a single array of canonical inner product strings using
a function called SymSub. The substitutions are performed in a
permutation independent fashion. The set of n_1 distinct off-diagonal
strings are lexicographically ordered and then 1 is assigned to the
first string in the ordered list and so on until n_1 is assigned
to the last string in the ordered list. Similarly, the list of n_2
distinct diagonal strings are ordered prior to assigning (m^2n_1+1)
to (m^2n_1+n_2) to them. This guarantees that
any symmetric permutation applied to the array of canonical inner
product strings will assign the same integer to the same string and
that the result is SPD.



 §.§ A Convex Set


Theorem <ref> shows that determining p-similarity between
two square matrices is equivalent to determining the p-similarity
between two PCMs. PCMs are symmetric matrices whose diagonal symbols
differ from off-diagonal symbols. Since the spectrum can be shifted
without affecting the p-similarity relationship, there exist SPD PCMs
with the same pattern. One way to generate a SPD PCM is to re-define
the color matrix to use γ=2m^2 instead of γ=2. Using
γ=2m^2 guarantees the PCM is SPD since it is symmetric and
diagonally dominant by construction. 

Given a m^2× m^2 SPD PCM M, let 𝕄={ N∈ SPD | Aut(N)=Aut(M)}
be the set of all real SPD matrices whose automorphism group is Aut(M).



 𝕄 is a convex set. That is for
all M1, M2∈𝕄, the linear combination α M1+(1-α)M2
is also in 𝕄 for 0≤α≤1.



Let M1 and M2 be two matrices in 𝕄. It is well known
that the set of real SPD matrices is convex, so the focus of the proof
is on showing that Aut(α M1+(1-α)M2)=Aut(M)
for 0≤α≤1.

⇒ For all P∈ Aut(M), PM1P^T=M1 and
PM2P^T=M2 since M1 and M2 are in 𝕄. So we have

    α PM1P^T+(1-α)PM2P^T    =    α M1+(1-α)M2
    
    P(α M1+(1-α)M2)P^T    =    α M1+(1-α)M2

 showing that Aut(M)≤ Aut(α M1+(1-α)M2).

⇐ For all P∈ Aut(α M1+(1-α)M2),
we have 

    P(α M1+(1-α)M2)P^T=α M1+(1-α)M20≤α≤1.



  * For α=0, (<ref>) equals PM2P^T=M2
so Aut(α M1+(1-α)M2)≤ Aut(M)
for α=0.

  * For 0<α≤1, (<ref>) can be rewritten
as 

    PM1P^T+γ PM2P^T=M1+γ M2

 where γ=(1-α)/α≥0. Rearrange
to get 

    PM1P^T-M1=γ(M2-PM2P^T).

 Now the matrix differences on the lhs and rhs of (<ref>)
are fixed for a given P and yet (<ref>)
holds for all γ≥0. This creates a contradiction unless
(PM1P^T-M1) and (M2-PM2P^T) are zero
matrices, implying PM1P^T=M1 and PM2P^T=M2. So Aut(α M1+(1-α)M2)≤ Aut(M)
for 0<α≤1.

Therefore Aut(α M1+(1-α)M2)=Aut(M)
for 0≤α≤1 so each α M1+(1-α)M2
is in 𝕄 implying 𝕄 is a convex set.

𝕄 is not closed under matrix multiplication. For any N∈𝕄,
its inverse N^-1 is also in 𝕄, but the automorphism
group of N× N^-1=I is S_m, the symmetric group, which
is not a subgroup of { P⊗ P}. However, the
powers of a matrix in 𝕄 are in 𝕄.




 §.§ 𝐌_𝐀𝐮𝐭 and 𝐌^* are in 𝕄


Matrices in 𝕄 can be separated into equivalence classes
by their patterns (partitions). Elements of a class differ from each
other by a symbol substitution. 

Let M be a SPD PCM. Define M_Aut as a SPD
matrix representing the orbits imposed by Aut(M) on
L, the set of (i,j) locations, see Section <ref>.
Locations in the same orbit are assigned the same symbol and locations
in distinct orbits have distinct symbols. Symbols on the diagonal
of M_Aut are distinct from its off-diagonal symbols since the
orbits of diagonal elements are distinct from the orbits of off-diagonal
elements. M_Aut is symmetric since M is symmetric and the
diagonal symbols are chosen so that M_Aut is SPD.


M_Aut is in 𝕄 since it is SPD and Aut(M_Aut)=Aut(M)
by the definition of M_Aut.

Let M be a SPD PCM. Let M^(i)=SymSub(SymSqr(M^(i-1)))
be the result of the i^ th symbolic squaring and symbol
substitution, where M^(0)=M. Let s be the iteration where
the stable pattern is reached, then Π(M^(s+j))=Π(M^(s))
for j=1,2,…. Define M^* as a SPD matrix whose pattern
matches M^(s), Π(M^*)=Π(M^(s)).
Each M^(i), i=1,…, is SPD by Section <ref>.


The next theorem shows that the M^(i) and M^* are in 𝕄.

Given SPD PCM M and 𝕄.
If M^(i)=SymSub(SymSqr(M^(i-1))) where
M^(0)=M and M^* is a SPD matrix representing the stable
pattern resulting from repeated symbolic squaring. Then the M^(i),
for i=1,…, and M^* are in 𝕄.



Given PCM SPD M and 𝕄, the set of all SPD matrices
whose automorphism group is Aut(M), to show that M^(i)
and M^* are in 𝕄 we need to show that they are SPD
and their automorphism groups match Aut(M). Since M^(i)
and M^* are SPD by Definition <ref>, the focus is on
showing that their automorphism groups match Aut(M).

Let M^(i)=SymSub(SymSqr(M^(i-1))) where
M^(0)=M. 

For i=1: M^(1)=SymSub(SymSqr(M^(0))).
To see that Aut(M^(1))=Aut(M), first note
that Theorem <ref> guarantees that 

    Π(M^(1))=Π(SymSub(SymSqr(M^(0))))≼Π(M^(0))=Π(M).


⇒ For each P∈ Aut(M^(1)), if P symmetrically
permutes (i,j) to (r,s) then M_i,j^(1)=M_r,s^(1)
and M_i,j=M_r,s since Π(M^(1))≼Π(M).
So PMP^T=M implying P∈ Aut(M) and Aut(M^(1))≤ Aut(M). 

⇐ For each P∈ Aut(M), if P symmetrically
permutes (i,j) to (r,s) then (i,j) and (r,s) are in the
same orbit and Theorem <ref> guarantees both
locations have identical canonical inner product strings. Symbol substitution
assigns the same symbol to identical canonical inner product strings
resulting in M_i,j^(1)=M_r,s^(1), so PM^(1)P^T=M^(1)
and P∈ Aut(M^(1)) implying Aut(M)≤ Aut(M^(1)). 

Therefore Aut(M^(1))=Aut(M).

By induction, the automorphism group of M^(i) is Aut(M)
for i=2,…. 

For M^*, let M^*=M^(s) where s is the iteration a stable
pattern is reached, i.e., Π(M^(s+1))=Π(M^(s)).
Then M^* is in 𝕄 since M^(s) is in 𝕄. 

Therefore all of the intermediate M^(i),i=1,…, and M^*
are in 𝕄.




 §.§ Lattice and Chain


The patterns (partitions) of the matrices in 𝕄 form equivalence
classes that can be partially ordered by refinement. Let M_sup
be a PCM in 𝕄 with the fewest distinct symbols/cells.
Then the partial ordering forms a complete lattice with Π(M_sup)
is the supremum and Π(M_Aut) as the infimum <cit.>. 

The subset of patterns {Π(M), Π(M^(1)), ⋯, Π(M^*)},
where

    Π(M^*)≺⋯≺Π(M^(1))≺Π(M)

 is a totally ordered set, called a chain, between Π(M)
and Π(M^*) <cit.>. 

Both Π(M^*) and Π(M_Aut) are fixed
point patterns (stable partitions) of symbolic squaring, 

    Π(SymSqr(M^*))     =    Π(M^*)

 and 

    Π(SymSqr(M_Aut))     =    Π(M_Aut).

Tarski's Lattice-Theoretical Fixed Point Theorem <cit.>
guarantees that the set of fixed points of symbolic squaring is not
empty, and that the fixed points form a complete lattice partially
ordered by refinement. Unfortunately it does not say that Π(M^*)
equals Π(M_Aut). If they are identical, then the
process of symbolic squaring PCMs finds the orbits.




 §.§ Symbolic Squaring as a Descent
Method


Symbolic squaring with symbol substitution monotonically refines a
PCM M until a stable pattern, Π(M^*), is reached.
Cells of Π(M^(i)) represent disjoint sets of orbits,
Theorem <ref>. 𝕄 is a convex set
composed of all SPD matrices whose automorphism group matches Aut(M).
Symbolic squaring creates a chain of strictly refined patterns from
Π(M) to Π(M^*). The question is
whether Π(M^*) is equal to Π(M_Aut).

The process of refinement can be viewed as a descent method. Given
M and 𝕄, for all P∈ Aut(M) one can
write 

    M× M-P(MP^T× PM)P^T=0.

 Equation <ref> can be converted to a pattern
difference that is symbol independent as 

    ‖Π(M× M)-Π(P(MP^T× PM)P^T)‖ _F^2=0

 where ‖‖ _F^2 is the square of the Frobenius
norm and the pattern difference assigns zero to locations with matching
representative locations and one to locations where the representative
locations do not match, see Section <ref>.
Convert (<ref>) to using symbolic squaring with symbol substitution
as 

    ‖Π(SymSub(SymSqr(M)))-Π(P(SymSub(SymSqr(M)))P^T)‖ _F^2=0

where the function SymSub() performs symbol substitution on a single
matrix using a permutation independent mapping from symbol to canonical
inner product string, as described in Section <ref>.

Let M_Aut be a SPD matrix representing the orbits imposed by
Aut(M), as defined in Definition <ref>. Let
M^(i)=SymSub(SymSqr(M^(i-1))) where Π(M^(0))=Π(M)
and M^(i)∈𝕄. Then the objective is to minimize the
cost function 

    ‖Π(SymSub(SymSqr(M_Aut)))-Π(P(SymSub(SymSqr(M^(i-1))))P^T)‖ _F^2=SSE

for i=1,⋯,* where SSE is the sum square error. Since Π(M_Aut)
is a stable pattern, and each M^(i) is in 𝕄, (<ref>)
can be rewritten as 

    ‖Π(M_Aut)-Π(M^(i))‖ _F^2    =     SSE.

 The SSE in (<ref>) decreases monotonically.
Each cell of Π(M^(i)) that contains multiple orbits,
has only one orbit with the correct representative location, the remaining
orbits in the cell have the wrong representative location. As symbolic
squaring refines a cell, at least one additional orbit will have the
correct representative location. So the SSE of the pattern difference
decreases every iteration until the stable pattern Π(M^*)
is reached.

Although we don't know Π(M_Aut) the monotonic refinement
of symbolic squaring along with the convexity of 𝕄 drives
the iterations towards Π(M_Aut). So an iterative
method given by 

    M^(i)=SymSub(SymSqr(M^(i-1)))

 that halts when 

    ‖Π(M^(i))-Π(M^(i-1))‖ _F^2=0

causes the successive patterns Π(M^(i)) to approach
Π(M_Aut). Algorithm <ref>
uses a variant of (<ref>), namely testing Π(M^(i))
and Π(M^(i-1)) for equality, to detect when the
process has reached a stable pattern, line 20 of Algorithm <ref>. 

In summary, repeated symbolic squaring with symbol substitution finds
the most refined pattern. The refinement is monotonic and converges
to a stable pattern where locations in an orbit are in the same cell.
By the convexity of 𝕄, Π(M^*) equals
Π(M_Aut). 

Given that symbolic squaring finds the orbits, one can construct PCM
M=PCM((A⊕ B)_C) from the color matrix
for the direct sum of the input matrices and check whether locations
on the diagonal of M^* associated with A and B are in the
same orbits. If not, A and B are not p-similar. The sufficiency
argument would end here. However, the author chooses to continue the
sufficiency argument because a) it will be shown that not only do
the orbits differ for non p-similar matrices but no orbits match,
and b) comparing multisets of diagonal elements to determine p-similarity
fits nicely with Specht's Theorem comparing traces to determine unitary
similarity.

Although the process of symbolic squaring a PCM finds its orbits,
it says nothing about the association of symbol to orbit. Two PCMs,
A and B, if processed separately can have the same automorphism
group Aut(A)=Aut(B) and partitions Π(A)=Π(B)
so the orbits match M_Aut(A)=M_Aut(B).
But if they have different symbol sets, Σ_A≠Σ_B
or the same symbol set with different mappings of symbol to cell g_A≠ g_B,
then the matrices are not p-similar. This issue is addressed in the
following sections.



 §.§ SymMult


Before addressing the issue of symbols to orbits, there is one more
symbolic matrix multiplication function that is useful. It is called
SymMult() and takes two diagonal matrices and a square
matrix as input. The first diagonal matrix is applied to the square
matrix from the left and the second from the right, SymMult(D1, M, D2).
For example, if 

    D1=[[ δ_1   0;   0 δ_2 ]]M=[[ M_1,1 M_1,2; M_2,1 M_2,2 ]]D2=[[ δ_1   0;   0 δ_2 ]]

 then SymMult(D1,M,D2) is the canonical inner product
string array 

    SymMult(D,M,D)=[[ δ_1M_1,1δ_1 δ_1M_1,2δ_2; δ_2M_2,1δ_1 δ_2M_2,2δ_2 ]].



Given square matrices A and B
whose diagonal symbols are distinct from off-diagonal symbol. If there
exists a permutation matrix P such that PAP^T=B and pairs
of column vectors (u,v) and (x,y) such
that Pu=v and Px=y. Then 

    P(SymMult(D_u,A,D_x))P^T=SymMult(D_v,B,D_y)

 where D_u=diag(u), D_v=diag(v), D_x=diag(x), and D_y=diag(y)
are diagonal matrices. Further, 

    Π(SymMult(D_u,A,D_x))≼Π(A)Π(SymMult(D_v,B,D_y))≼Π(B).




Assume PAP^T=B and (u,v) and (x,y)
are column vector pairs such that Pu=v and Px=y. Then PD_uP^T=D_v
and PD_xP^T=D_y where D_u=diag(u), D_v=diag(v),
D_x=diag(x), and D_y=diag(y) are diagonal matrices. 

Using regular matrix multiplication, multiplying PAP^T=B on the
left and right yields

    PAP^T    =     B
    (PD_uP^T)PAP^T(PD_xP^T)     =     D_vBD_y
    
    P(D_uAD_x)P^T    =     D_vBD_y

so P(D_uAD_x)P^T=D_vBD_y. 

From (<ref>), the output arrays of canonical inner product
strings of 

    SymMult((PD_uP^T),(PAP^T),(PD_xP^T))SymMult(D_v,B,D_y)

 are identical. To go from (<ref>) to (<ref>)
requires 


    P(SymMult(D_u,A,D_x))P^T=SymMult((PD_uP^T),(PAP^T),(PD_xP^T)).

 This is true since symmetric permutations move values around without
changing them. So a symmetric permutation can move strings as easily
as real or complex numbers. Therefore, 


    P(SymMult(D_u,A,D_x))P^T=SymMult(D_v,B,D_y).


It is easy to see from (<ref>) that Π(SymMult(D_u,A,D_x))≼Π(A).
Let α be the symbol associated with one cell of A. Locations
within that cell are the only locations that will have α as
the middle symbol on the lhs of (<ref>). Therefore,
Π(SymMult(D_u,A,D_x)) is a refinement
of Π(A). A similar argument holds for Π(SymMult(D_v,B,D_y))
and Π(B). Therefore Π(SymMult(D_u,A,D_x))≼Π(A)
and Π(SymMult(D_v,B,D_y))≼Π(B)
as desired.

For PCMs it can be shown that the pattern from the first symbolic
squaring is the same as symbolically applying the diagonal from the
left and right.

Given m^2× m^2 PCM M=D+R where
D is the diagonal and R=I⊗(1*(J-I))+(J-I)⊗(2*I)
are the off-diagonal values. Then 

    Π(SymMult(D,M,D))=Π(SymSqr(M)).




The result is by careful consideration of the canonical inner product
strings for each (i,j) location. 

First convert the off-diagonal values of R into symbols. The off-diagonal
values of R are `1', `2', and `0'. The locations associated with
`0' are not explicitly called out in the equation for R. Rewrite
R to include those locations, 

    R=I⊗(1*(J-I))+(J-I)⊗(2*I+0*(J-I)).

Now substitute symbols in place of `1', `2', and `0' to get 

    R=I⊗(σ_1*(J-I))+(J-I)⊗(σ_2*I+σ_3*(J-I)).

 R has numeric zeros on its diagonal because the splitting
M=D+R places the diagonal terms of M in D. 

Let D be a diagonal matrix where symbols have been substituted
for values on the diagonal. All off-diagonal locations of D
are a numeric zeros, due to the splitting. The i^ th diagonal
symbols is denoted as δ_i, i.e., D_i,i=δ_i.

Let M=D+R. Symbolically squaring
M can be written as 

    SymSqr(M)     =    M×M
        =    (D+R)×(D+R)
        =    DD+DR+RD+RR

where the rhs abuses notation by using regular matrix notation applied
to arrays of symbols. Terms in the canonical inner product strings
involving one or more diagonal symbols come from DD+DR+RD.
Terms in the canonical inner product strings involving two off-diagonal
symbols come from RR.

Table <ref> shows all of the terms
making up the canonical inner product strings from symbolically squaring
M. The first column in Table <ref>
shows the symbol associated with an (i,j) location in M.
The second column shows the term(s) involving a diagonal symbol, coming
from DD+DR+RD.
The third columns shows all of the terms involving two off-diagonal
symbols, coming from RR. 

To use Table <ref> to construct
the canonical inner product string at (i,j), first find the row
matching the symbol at M_i,j in the first column.
Then using that row, concatenate the term(s) in the second column
followed those in the third column after lexicographically ordering
the terms in the third column. Note that the term(s) in the second
column are already in the order specified for terms involving diagonal
symbols. 


Examining the third column of Table <ref>,
one notes that the terms come in matching sets. For example, using
the row with σ_2 in the first column, the third column has
(m-1) copies of σ_1σ_3 and (m-1) copies of
σ_3σ_1. In all cases if a row in column three contains
k copies of αβ then it also contains k copies of
βα. So the contribution to the canonical inner product
strings at (i,j) and (j,i) from (RR)_i,j
is identical to the contribution from (RR)_j,i
implying RR is symmetric. Further, one sees
from Table <ref> that each off-diagonal
symbol of R (the σ_k in column one) is associated
with a distinct multiset of terms in RR (column
three of Table <ref>). The diagonal
of RR is associated with a fourth multiset
of terms that is distinct from the three off-diagonal multisets. So,
(σ_4I+R) and RR
have identical patterns,

    Π(σ_4I+R)=Π(RR),

 where the symbol σ_4 is chosen to be distinct from all
other symbols in M.

Now focus on the second column of Table <ref>.
Assume that M_i,j=σ_1 and M_r,s=σ_1
where (i,j)≠(r,s). Then the only difference in the canonical
inner product strings are the terms involving diagonal symbols, δ_iσ_1+σ_1δ_j
and δ_rσ_1+σ_1δ_s respectively. Further,
note that the row associated with an off-diagonal symbol use the same
off-diagonal symbol in column two. So if column one has off-diagonal
symbol σ_k, k∈{ 1,2,3}, then column two
has the terms δ_iσ_k+σ_kδ_j. Each two
factor pair of terms δ_iσ_k+σ_kδ_j
could be replaced by a single three factor term δ_iσ_kδ_j
without changing its ability to distinguish canonical inner product
strings. Further, δ_iδ_i can be replaced by δ_iσ_4δ_i
without changing its ability to distinguish canonical inner product
strings. Performing these substitutions allows column two to be expressed
as D(σ_4I+R)D.
Now comparing the pattern of D(σ_4I+R)D
to the pattern of (DD+DR+RD),
one sees that they are equal,

    Π(D(σ_4I+R)D)=Π(DD+DR+RD).


Note that Π(D(σ_4I+R)D)≼Π(σ_4I+R)
and 

    Π(D(σ_4I+R)D)=Π(D(D+R)D̅)=Π(DMD)

since Π(DDD)=Π(D(σ_4I)D).
Therefore, 

    Π(SymMult(D,M,D))=Π(SymSqr(M))

 for PCM M as was to be shown.

Theorem <ref> only applies to the first symbolic squaring
of a PCM. Attempting subsequent symbolic squaring by symbolically
multiplying on the left and right by the new diagonal D^3 does
not refine the pattern since it has the same pattern as D, Π(D^3)=Π(D).



 §.§ PCMs and Oracles


In Algorithm <ref>, a PCM is constructed
for each input matrix, A and B. Let S=PCM(A_C) and T=PCM(B_C).
Then symbolic squaring is applied to S and T simultaneously
and the canonical inner product strings compared across products.
This section examines how many columns an Oracle can match between
S^* and T^* when restricted to permutations of the form
P⊗ P. 

Given m^2× m^2 PCMs S and T
and their refined patterns Π(S^*) and Π(T^*)
respectively. An Oracle asked to match columns of S^*and T^*
using permutations of the form P⊗ P will either match all
of the columns or none of the columns.



Proof by contradiction. Assume SPT and the Oracle
matches k, 1≤ k<m^2, columns of S^* and T^* using
a permutation P_O∈{ P⊗ P}. Then k columns
of P_OS^*P_O^T and T^* match. Symmetrically permute
the matched columns to the lhs using P_lhs to get 

    P_lhsP_OS^*P_O^TP_lhs^T=[[   C1 C2^T;   C2 S2^* ]]P_lhsT^*P_lhs^T=[[   C1 C2^T;   C2 T2^* ]]

 where 

    [[ C1; C2 ]]

 are the matched columns. Both P_lhsP_OS^*P_O^TP_lhs^T
and P_lhsT^*P_lhs^T are stable patterns, so symbolic squaring
does not change the pattern.

Now apply the permutations to S and T to get 

    S^(0)=P_lhsP_OSP_O^TP_lhs^T

 and 

    T^(0)=P_lhsTP_lhs^T

 which have a conforming block structure 

    S^(0)=[[     C1^(0) (C2^(0))^T;     C2^(0)     S2^(0) ]]T^(0)=[[     C1^(0) (C2^(0))^T;     C2^(0)     T2^(0) ]].

S^(0) and T^(0) have identical off-diagonal structure. The
application of P_O does not change the off-diagonal structure
because P_O∈{ P⊗ P} is an automorphism
of R. Then P_lhs is applied to P_OSP_O^T and T,
so the off-diagonal structure changes the same way for both. The diagonal
of C1^(0) for S^(0) and T^(0) must match. If the diagonals
in C1^(0) differ, they will also differ in C1 of P_lhsP_OS^*P_O^TP_lhs^T
and P_lhsT^*P_lhs^T creating a contradiction.

Now consider the first symbolic squaring. From Theorem <ref>
the pattern from the first symbolic squaring matches the pattern from
SymMult(D,M,D). This still holds after symmetrically
permuting columns to the lhs. So SymSqr(S^(0)) and
SymSqr(T^(0)) have patterns matching the patterns
from SymMult(D_S^(0),S^(0),D_S^(0)) and SymMult(D_T^(0),T^(0),D_T^(0))
respectively. In block form that is 

    D_S^(0)S^(0)D_S^(0)=[[ D_C1^(0);          D_S2^(0) ]]×[[     C1^(0) (C2^(0))^T;     C2^(0)     S2^(0) ]]×[[ D_C1^(0);          D_S2^(0) ]]

 and 

    D_T^(0)T^(0)D_T^(0)=[[ D_C1^(0);          D_T2^(0) ]]×[[     C1^(0) (C2^(0))^T;     C2^(0)     T2^(0) ]]×[[ D_C1^(0);          D_T2^(0) ]].


Focus on what happens to C2^(0). We see that the patterns resulting
from symbolic multiplication of D_S2^(0)C2^(0)D_C1^(0)
and D_T2^(0)C2^(0)D_C1^(0) do not match since D_S2^(0)≠ D_T2^(0).
Recall that S and T are not p-similar by hypothesis. Therefore
the C2 from P_lhsP_OS^*P_O^TP_lhs^T and P_lhsT^*P_lhs^T
cannot match, implying the Oracle failed to match any columns of S^*
and T^* using a permutation of the form P⊗ P.

So the Oracle either matches all columns of S^* and T^*
or none of the columns using amenable permutations.

There is no partial matching of S^* and T^* using a permutation
of the form P⊗ P. Either S^* and T^* are p-similar
and the columns can be matched or S^* and T^* are not p-similar
and no columns can be matched using a permutation of the form P⊗ P.
This is independent of comparing multisets of diagonal symbols for
S^* and T^*.



 §.§ Symbolic Squaring of 𝐏𝐂𝐌((𝐀⊕ 𝐁)_𝐂)


Theorem <ref> showed that given two matrices A and
B and their PCMs S and T, where S=PCM(A_C)
and T=PCM(B_C), that either S^* is p-similar
to T^* or none of the columns can be matched using an admissible
permutation of the form P⊗ P. However, it does not say anything
about diagonal symbols.

Given non p-similar m× m matrices A and B, construct
the PCM of the direct sum, M=PCM((A⊕ B)_C).
Also construct the PCMs S=PCM(A_C) and T=PCM(B_C).
The color matrix of the direct sum (A⊕ B) is given
by 

    (A⊕ B)_C=[[ A_C  3J;  3J B_C ]].

 Applying column major ordering to construct the diagonal of M
results in 2m×2m diagonal blocks where the k^ th
diagonal block in the upper left quadrant of M is given by 

    [[ (diag(A_C(:,k))+(J-I))                      J;                      J                 (J+2I) ]]_2m×2m

 and the (m+k)^ th diagonal block in the lower
right quadrant of M is given by 

    [[                 (J+2I)                      J;                      J (diag(B_C(:,k))+(J-I)) ]]_2m×2m

 where J is the m× m matrix of all ones and I is the
conforming identity matrix. Off-diagonal blocks of M are given
by 2I_2m×2m. 

So, all of the entries in S and T are embedded in M. The
k^ th m× m diagonal block of S, (diag(A_C(:,k))+(J-I)),
is the upper left portion of each 2m×2m diagonal block in
the upper left quadrant of M, see (<ref>). The m× m
off-diagonal blocks of S are the upper left portion of each 2m×2m
off-diagonal blocks in the upper left quadrant of M. A similar
case holds for T except the blocks occupy the lower right portion
of each 2m×2m block in the lower right quadrant of M.

Perform symbolic squaring on M until the stable pattern Π(M^*)
is reached. By the argument in Section <ref>,
Π(M^*)=Π(M_Aut) so each cell of M^*
represents an individual orbit. The patterns associated with locations
of S and T in M^* are identical to Π(S^*)
and Π(T^*) since Π(S^*) and Π(T^*)
are the orbits of S and T respectively. Therefore, Theorem <ref>
applies to the columns of M^* associated with S and T. 

Since A and B are non p-similar, none of the columns in M^*
associated with S and T can be matched. Therefore the multisets
of diagonal symbols corresponding to locations associated with S
and T in M^* are distinct. Then since the patterns in M^*
associated with S and T are identical to Π(S^*)
and Π(T^*), the diagonal mixes from S^* and
T^* are different as long as consistent symbol substitution is
performed each symbolic squaring. 

This is the final piece showing that comparing the multisets of diagonal
symbols from S^* and T^* is necessary and sufficient to
determine whether A and B are p-similar.



§ FINDING 𝐏 USING THE BLIND ALGORITHM


Given that Algorithm <ref> is a
polynomial blind permutation similarity algorithm, it can be used
as part of a second polynomial algorithm to find a permutation between
p-similar matrices. Such an algorithm is presented in Algorithm <ref>.
The remainder of this section describes the concepts behind Algorithm
<ref>.





 §.§ Overall Concept


Algorithm <ref> uses the blind permutation similarity
algorithm, Algorithm <ref>, to find
a permutation vector p, such that M1(p,p)=M2 when M1 and
M2 are p-similar. If M1 and M2 are not p-similar, then psim=FALSE
and p=[1:m]^T, the identity permutation vector, are
returned as output. 

The main loop in Algorithm <ref>, lines 23-43, looks
for the column to use as p(c), where c is the iteration variable.
The preconditions for each iteration are that i) A and B are
p-similar, ii) their diagonal symbols are distinct from the off-diagonal
symbols, and iii) the first (c-1) locations of p,
p(1:(c-1)), are correct. 

Upon entering the loop for the first time with c=1, the preconditions
are satisfied since:


  * A and B are p-similar by the call to Algorithm <ref>
on line 13 of Algorithm <ref>. 

  * The transformations applied in lines 16-18 are similar to those constructing
color matrices from the input matrices. So the diagonal symbols of
A and B are distinct from off-diagonal symbols. 

  * (c-1)=0 columns have been matched correctly, satisfying
the third precondition. 

For each iteration, A and B are the remaining trailing principal
submatrices to be matched. Only columns of A with a diagonal symbol
A_j,j=B_1,1 are potential matches for the first column of B.
The nested loop, lines 25-41 searches down the diagonal of A for
columns where A_j,j=B_1,1. When it finds such a column, it
is symmetrically permuted to become the first column of AJ1. In
Algorithm <ref> the exchange happens on line 28. The
exchange can be expressed as 

    AJ1=P_1jAP_1j^T

 where permutation matrix P_1j exchanges rows and columns 1
and j. Assuming the first column of AJ1 is a correct match for
the first column of B. Then there exists a complementary permutation
P_c such that 

    P_c(AJ1)P_c^T=B

 where the structure of P_c is given by 

    P_c=[[   1;     P22 ]].

 Then 

    [[   1;     P22 ]]×[[     AJ1_1,1   AJ1_1,2:n;   AJ1_2:n,1 AJ1_2:n,2:n ]]×[[   1;     P22 ]]^T=[[     B_1,1   B_1,2:n;   B_2:n,1 B_2:n,2:n ]]

 where 

    P22× AJ1_2:n,1    =     B_2:n,1
    
    AJ1_1,2:n× P22^T    =     B_1,2:n

 and 

    P22× AJ1_2:n,2:n× P22^T    =     B_2:n,2:n.


Although P22 is unknown, Algorithm <ref>
can be used to determine if AJ1_2:n,2:n and B_2:n,2:n in
(<ref>) are p-similar. However, just testing the
(2:n,2:n) trailing principal submatrices for p-similarity
may give a false result. It may be that the only permutation(s), that
match the trailing principal submatrices, violate (<ref>)
and (<ref>). To eliminate this possibility, the
first column and first row of AJ1 and B are used to refine the
(2:n,2:n) trailing principal submatrices prior to checking
if they are p-similar. For AJ1, this is done by constructing diagonal
matrices from the first column and first row of AJ1, diag(AJ1_:,1)
and diag(AJ1_1,:) respectively, and symbolically applying
them from the left and the right to AJ1, line 29. The (2:n,2:n)
trailing principal submatrix of B is similarly refined using the
first column and first row of B, line 24. Theorem <ref>
guarantees that the symbolic multiplication will refine the (2:n,2:n)
trailing principal submatrices while preserving the p-similarity relationship
between AJ1 and B. The left and right diagonal matrices are
constructed from the first column an first row respectively because
A and B might not be symmetric. 

Algorithm <ref> is applied to S22
and T22, the refined (2:n,2:n) trailing principal submatrices
of DAJ1D and DBD respectively, to test for p-similarity without
the possibility of a false response, line 35. If S22 and T22
are p-similar, then the amenable permutation also satisfies (<ref>)
and (<ref>). 

After finding a S22 and T22 that are p-similar, the global index
of the column permuted into the first column of AJ1 is saved as
p(c), line 37, and A and B are set to S22 and T22 respectively,
line 38. The size of the trailing principal submatrix is decremented
in preparation for searching for the next column, line 42. 

P-similar S22 and T22 are guaranteed to exist because the precondition
for entering the loop is that A and B are p-similar. 

For the next iteration c=2, again the preconditions are met. A
and B are p-similar, tested on line 35, their diagonal symbols
are distinct from off-diagonal symbols since the inputs to the symbol
substitution, line 30, have distinct diagonal symbols, and the first
(c-1) locations in p, p(1), are correct. 

Therefore, for p-similar input matrices M1 and M2, Algorithm
<ref> finds a permutation vector p such that M1(p,p)=M2
in O(m^2) iterations.



 §.§ Algorithm <ref> Complexity


Given m× m input matrices, the nested loops in Algorithm <ref>
make O(m^2) calls to BPSAY, Algorithm <ref>.
So the overall complexity of Algorithm <ref> is O(m^14)
using symbolic squaring. If widely space primes matrices are used,
Section <ref>, the overall complexity drops
to O(m^12). Therefore finding the permutation between
two p-similar matrices is polynomial.



§ DISCUSSION


In practice, it does not appear to be necessary to add edge weights
to the PCG. The edge weights are added to eliminate the possibility
of exchanging the rows with the columns. Edge weights simplified the
arguments by avoiding the need to carry the transpose based permutations
through the proofs. If rook's graphs are compact then color refinement
works on PCMs with equal edge weights <cit.>.
Although rook's graphs have many of the characteristics of compact
graphs, the author has not been able to show that rook's graphs are
compact.



 §.§ Origin of PCGs


For a period of time the author tried to develop a graph isomorphism
algorithm. For every algorithm, non isomorphic graphs were found that
were not discriminated. Eventually the investigation turned to looking
at why color refinement fails.

Given two non-isomorphic graphs, iterative color refinement can be
applied to the adjacency matrices until they converged to a stable
partition <cit.>. It is well known that strongly regular
graphs cannot be separated using naive color refinement. However,
one could ask an Oracle to maximally column match the matrices. The
result would look like 

    A=[[   C1 C2^T;   C2   A2 ]]B=[[   C1 C2^T;   C2   B2 ]]

 after symmetrically permuting the matched columns to the lhs. Differences
are confined to the lower right blocks, where A2≠ B2. We can
assume the column multisets match (otherwise A and B are known
to be non-isomorphic). Next consider ordering C2 into groups of
identical rows. The permutation to do this can be applied symmetrically
to both A and B. Now the matrices look something like 

    A=[[    C1 r_1^T r_2^T r_3^T r_4^T;   r_1     δ     γ     α     β;   r_2     γ     δ     β     α;   r_3     α     β     δ     γ;   r_4     β     α     γ     δ ]]B=[[    C1 r_1^T r_2^T r_3^T r_4^T;   r_1     δ     γ     β     α;   r_2     γ     δ     α     β;   r_3     β     α     δ     γ;   r_4     α     β     γ     δ ]]

 where α≠β and each r_j is a pairwise distinct
row vector. Consider the first column of the unmatched region. Notice
that A has α at the intersection of r_3 and r_1^T
whereas B has a β at that location and the reverse at the
intersection of r_4 and r_1^T. By hypothesis, r_1≠ r_3,
r_3≠ r_4, and r_4≠ r_1, so the author looked for
ways to take advantage of this situation by constructing invariants
such as all of the four-cycles. Four-cycles are permutation independent.
The process always broke down when trying to compare results using
multisets. The multisets either maintained row constraints and relaxed
column constraints or maintained column constraints and relaxed row
constraints. This happens for any comparison method based on multisets
such as those described in <cit.>.

This led to asking the question: “Can symmetric permutation constraints
be added to the original graphs?” If so, the symmetric permutation
constraints would be part of the original problem and using multisets
for comparison should not experience the same difficulty. PCGs were
the result. PCGs also changed the focus from graph isomorphism to
permutation similarity.



 §.§ Refinement using Regular Matrix
Multiplication


It is not clear that the sequence of refinements generated by symbolic
squaring can be replicated using regular matrix multiplication. If
the sequence cannot be replicated, it makes arguments using the existing
machinery for SPD matrices more difficult to accept. It turns out
there are two different ways the sequence of refinements can be replicated
using regular matrix multiplication. They both use matrices constructed
from widely-spaced prime numbers. 

A widely-spaced primes matrix
W, where W=wspm(M) is constructed from a symmetric
matrix M, whose diagonal symbols are distinct from off-diagonal
symbols. The construction uses the following method:


  * Let M=(Π, Σ, g) be a m× m symmetric
matrix whose diagonal symbols are distinct from the off-diagonal symbols.

  * Assume M has n=|Σ| distinct symbols where there
are n_1 off-diagonal symbols and n_2 diagonal symbols with
n_1+n_2=n.

  * Let Π(W)=Π(M).

  * Choose n prime numbers, p_i, to be the symbols of W using
the following recurrence.


  * Let p_1>mk^2 be the first prime number, where k∈ℕ^+
(typically k=1).

  * For i=2,⋯,n, chose p_i>m*p_i-1^2.


  * Use the first n_1 primes, p_1,…,p_n_1, as off-diagonal
symbols of W and use the remaining n_2 primes, p_n_1+1,…,p_n,
as diagonal symbols of W. 




A widely-spaced primes matrix W is SPD since it is symmetric and
strictly diagonally dominant by construction. 


Consider squaring a m× m widely-spaced primes matrix W.
Let W=D+O where O is the matrix of off-diagonal symbols and
D is the diagonal. Recall that widely-spaced primes matrix has
diagonal symbols that are distinct from off-diagonal symbols. Using
regular matrix multiplication, the value at (i,j) is given by 

    [W× W]_i,j=
    D_i,iO_i,j+O_i,jD_j,j+∑_[ k≠ i; k≠ j ]O_i,kO_k,j   
    
    D_i,iD_i,i+∑_[ k≠ i ]O_i,kO_k,i

 where terms involving diagonal symbols come first, followed by terms
with two off-diagonal symbols. This is similar to the order used for
canonical inner product strings. 

Each term in the inner product [W× W]_i,j is
the product of two primes.


 

Let Σ={ p_1,⋯,p_n}
be the set of widely-spaced primes in W. Then the set of all possible
distinct terms, 

    { p_1p_1, p_1p_2, p_2p_2, p_1p_3, p_2p_3, p_3p_3,⋯,p_1p_n,⋯,p_np_n}

is totally ordered by less than 

    p_1p_1<p_1p_2<p_2p_2<p_1p_3<p_2p_3<p_3p_3<⋯<p_1p_n<⋯<p_np_n

and the ratio between two adjacent terms is greater than mp_1,

    p_i+1p_j/p_ip_j>mp_ij>ip_1p_j+1/p_jp_j>mp_1i=j.




The inner product [W× W]_i,j
is unique to a set of distinct terms and their multiplicities.



Let W be a m× m widely-spaced primes matrix. Then each
inner product [W× W]_i,j has m terms and can
be written as 

    [W× W]_i,j=∑ m_kt_k

where each t_k is a distinct term in the inner product with multiplicity
m_k and ∑ m_k=m. 

Proof is by contradiction. Assume [W× W]_i,j=[W× W]_r,s
where 

    [W× W]_r,s=∑m_lt_l

 and ∑ m_k=∑m_l=m but the sets of distinct
terms { t_k} and {t_l}
do not match, { t_k}≠{t_l}.
Subtract terms appearing in both [W× W]_i,j and
[W× W]_r,s from ∑ m_kt_k and ∑m_lt_l.
Let the resultant sums be σ_ij and σ_rs
respectively. Now σ_ij=σ_rs>0 but they
have no common terms. WLOG assume that σ_ij contains the
largest distinct term between σ_ij and σ_rs.
Let p_k1p_k2 be that largest term in σ_ij and p_l1p_l2
be the largest term in σ_rs. So p_k1p_k2>p_l1p_l2.
All of the terms in σ_rs are less than or equal
to p_l1p_l2. Therefore the largest σ_rs
can be is mp_l1p_l2. But by construction, the next largest
distinct term greater than p_l1p_l2 is strictly greater than
mp_l1p_l2, see Remark <ref>. Therefore σ_ij
cannot equal σ_rs contradicting the assumption
that two inner products of W× W that are equal can be constructed
from different sets of distinct terms with different multiplicities. 

A consequence of Lemma <ref> is that the inner
product [W× W]_i,j is decomposable using modulo
arithmetic by starting with the largest prime and working downward,
recursively applying modulo arithmetic to both the integer multiple
of the prime and the remainder. 

The multiplicity m_k of a distinct term, t_k=p_lp_r=p_rp_l,
in the inner product [W× W]_i,j is equal to the
sum of the multiplicities n_1 of (p_l,p_r)
and n_2 of (p_l,p_r) in the canonical inner
product string SymSqr(W)_i,j, m_k=n_1+n_2.



Consequence of Lemma <ref>. The inner product
is uniquely determined by the set of distinct terms in the inner product
and their multiplicities. A distinct term, t_k, is the product
of a unique pair of primes, t_k=p_lp_r=p_rp_l. Canonical
inner product strings distinguishes between the ordered pairs (p_l,p_r)
and (p_l,p_r) whereas the inner product does not.
The result follows.

The next theorem shows that the square of a widely-spaced primes matrix
is a refinement of the widely-spaced primes matrix.

Given a real m× m symmetric matrix
M whose diagonal symbols are distinct from off-diagonal symbols,
there exists a widely-spaced primes matrix W=wspm(M)
such that Π(W)=Π(M) and Π(W× W)≼Π(W).



Let M be a real m× m symmetric matrix whose diagonal symbols
are distinct from off-diagonal symbols. Let W=wspm(M)
be a widely-spaced primes matrix. Then Π(W)=Π(M)
by Definition <ref>. 

Let W=D+O where O is the matrix of off-diagonal symbols and
D is the diagonal. To see that Π(W× W)≼Π(W),
it is enough to show that entries that differ in W also differ
in W× W. Consider two off-diagonal locations (i,j) and
(r,s) such that W_i,j≠ W_r,s, i.e., O_i,j≠ O_r,s.
Without loss of generality, assume i<j and r<s. Then using the
first row of (<ref>) 

    [W× W]_i,j    =     D_i,iO_i,j+O_i,jD_j,j+∑_[ k≠ i; k≠ j ]O_i,kO_k,j

 and 

    [W× W]_r,s    =     D_r,rO_r,s+O_r,sD_s,s+∑_[ k≠ r; k≠ s ]O_r,kO_k,s

 where the terms involving two off-diagonal symbols sum to less than
the smallest diagonal symbol p_n_1+1. So if 

    O_i,j(D_i,i+D_j,j)O_r,s(D_r,r+D_s,s)

 then [W× W]_i,j is guaranteed to be distinct
from [W× W]_r,s. To see that the two sides of
(<ref>) are distinct when O_i,j≠ O_r,s, consider
the following cases:


  * (D_i,i+D_j,j)=(D_r,r+D_s,s): Then
O_i,j(D_i,i+D_j,j)≠ O_r,s(D_r,r+D_s,s)
since O_i,j≠ O_r,s .

  * (D_i,i+D_j,j)≠(D_r,r+D_s,s): Without
loss of generality, assume that (D_i,i+D_j,j)<(D_r,r+D_s,s).
Further, assume that D_s,s is maximal among the diagonal symbols
under consideration, D_s,s=max(D_i,i,D_j,j,D_r,r,D_s,s),
and that D_s,s=p_l, the l^ th prime number. 


  * (D_i,i, D_j,j<D_s,s) and (D_r,r≤ D_s,s): Substitute
into (<ref>) to get O_i,j(D_i,i+D_j,j)≤2p_n_1p_l-1<D_s,s
and O_r,s(D_r,r+D_s,s)>D_s,s. Therefore O_i,j(D_i,i+D_j,j)≠ O_r,s(D_r,r+D_s,s).

  * (D_i,i<D_r,r & D_j,j=D_s,s) and (D_r,r<D_s,s):
Substitute into (<ref>) and rearrange to get O_ijD_i,i-O_r,sD_r,r
versus (O_r,s-O_i,j)D_s,s. Now |O_ijD_i,i-O_r,sD_r,r|<2p_n_1p_l-1<D_s,s
and |O_r,s-O_i,j|D_s,s>D_s,s. Therefore O_i,j(D_i,i+D_j,j)≠ O_r,s(D_r,r+D_s,s).

  * (D_i,i<D_r,r & D_j,j=D_s,s) and (D_r,r=D_s,s):
Substitute into (<ref>) to get O_i,jD_i,i versus
(2O_r,s-O_i,j)D_s,s. Similar to the prior case,
O_i,jD_i,i<2p_n_1p_l-1<D_s,s and |2O_r,s-O_i,j|D_s,s>D_s,s.
Therefore O_i,j(D_i,i+D_j,j)≠ O_r,s(D_r,r+D_s,s).


A similar argument holds when diagonal symbols W_i,i and W_r,r
differ. In this case, only the D_i,i^2 and D_r,r^2 terms
need to be considered since the remaining terms sum to less than the
smallest diagonal symbol p_n_1+1. Inner products for diagonal
and off-diagonal locations differ in W× W. Off-diagonal locations
have two terms involving a diagonal symbol and an off-diagonal symbol
whereas diagonal locations only have one term involving the square
of a diagonal symbol, as shown in (<ref>). Again, the separation
between primes guarantees they cannot be equal. Therefore if two locations
of W differ, those locations will also differ in W× W implying
Π(W× W)≼Π(W) as desired.

Theorem <ref> does not say that Π(W× W)=Π(SymSqr(M)).
There may be some inner products that symbolic squaring distinguishes
that regular matrix multiplication does not. However, if SymSqr(W)
refines the diagonal, then W× W also refines the diagonal.

For a widely-space primes matrix W,
if SymSqr(W) strictly refines the diagonal, 

    Π(diag(SymSqr(W)))≺Π(diag(W)),

 then W× W also strictly refines the diagonal, 

    Π(diag(W× W))≺Π(diag(W))

 and Π(diag(W× W))=Π(diag(SymSqr(W))).



The inner product for a diagonal symbol has the same symbol for both
factors of a term, see the second row of (<ref>). Therefore,
the factor order doesn't matter for the canonical inner product strings
of diagonal symbols and Lemma <ref> says the
inner product is unique to the set of distinct terms and their multiplicities.
So, if SymSqr(W) refines the diagonal, Π(diag(SymSqr(W)))≺Π(diag(W)),
then W× W refines the diagonal, Π(diag(W× W))≺Π(diag(W))
and 

    Π(diag(W× W))=Π(diag(SymSqr(W)))

.

For off-diagonal locations, there may be cases where SymSqr(W)
refines an off-diagonal location but W× W does not. However,
if a widely-spaced primes matrix W is generated from PCM M,
M=PCM(A_C) and W=wspm(M), then for the
first squaring, Π(W× W)=Π(SymSqr(M)).
This can be seen by looking at Table <ref>
and noting that each canonical inner product string with (m-1)αβ's
also has (m-1)βα's, so each distinct canonical inner product
string in SymSqr(M) is associated with a distinct numeric
inner product of W× W. Recall that symbolic squaring uses
the lexicographically lessor to represent locations (i,j) and (j,i)
so switching the roles of δ_i and δ_j in column
2 of Table <ref> don't matter.

Next is to look at how close squaring a widely-space primes matrix
comes to matching a canonical inner product string for off-diagonal
locations after the first symbolic squaring. Start by looking at off-diagonal
locations that are not refined by symbolic squaring.

Let W be a m× m widely-spaced
primes matrix. If SymSqr(W)_i,j=SymSqr(W)_r,s
then (W× W)_i,j=(W× W)_r,s.



Given SymSqr(W)_i,j=SymSqr(W)_r,s, then
the canonical inner product strings at (i,j) and (r,s) match,
so the multisets of terms match. Therefore, the inner products (W× W)_i,j
and (W× W)_r,s match since they are summing over
the same multiset of terms.

Next is to look at off-diagonal locations that do get refined by symbolic
squaring. This is the more complex case. Before addressing the general
case, the instance where column mixes of the inner product vectors
match but the results from symbolic squaring differ is examined.

Let W be a m× m widely-spaced
primes matrix. Let (i,j) and (r,s) be two distinct off-diagonal
locations. If colMix([W_i,W_j])=colMix([W_r,W_s])
and SymSqr(W)_i,j≠ SymSqr(W)_r,s,
then (W× W)_i,j≠(W× W)_r,s.



Assume W is a m× m widely-spaced primes matrix constructed
as described in Definition <ref> and (i,j)
and (r,s) are two distinct off-diagonal locations such that colMix([W_i,W_j])=colMix([W_r,W_s])
and SymSqr(W)_i,j≠ SymSqr(W)_r,s.
WLOG assume that colMix(W_i)=colMix(W_r)
and colMix(W_j)=colMix(W_s). 

Proof by contradiction. Assume that (W× W)_i,j=(W× W)_r,s.
Let A=P_ijWP_ij^T where P_ij symmetrically permutes
columns i and j to be columns 1 and 2 of A respectively.
Similarly, let B=P_rsWP_rs^T where P_rssymmetrically
permutes columns r and s to be columns 1 and 2 of B respectively.

Note that A and B are symmetric and colMix(A_1)=colMix(B_1)
and colMix(A_2)=colMix(B_2). So there
exists a permutation P_A that when symmetrically applied to A,
sorts rows 3:m of A_1 into contiguous groups of distinct symbols.
Similarly there exists a permutation P_B that sorts rows 3:m
of B_1 into the same contiguous groups of distinct symbols. Let

    A← P_AAP_A^TB← P_BBP_B^T.

Now A_1=B_1 and A and B have the form shown in (<ref>)
where the row signatures are the groups of distinct symbols in rows
3:m of A_1 and B_1. Note that colMix(A_3:m,2)=colMix(B_3:m,2)
because colMix(A_2)=colMix(B_2) and A
and B are symmetric, so the first two elements of A_2 and
B_2 are identical.

By hypothesis SymSqr(W)_i,j≠ SymSqr(W)_r,s
so SymSqr(A)_1,2≠ SymSqr(B)_1,2 but
A_1^TA_2=B_1^TB_2. Consider the terms in the canonical
inner product strings of SymSqr(A)_1,2 and SymSqr(B)_1,2.
The first two terms involving diagonal symbols match, so the differences
must be in the terms from rows 3:m, 

    {{[ (A_3,1,A_3,2);             ⋮; (A_m,1,A_m,2) ]}}≠{{[ (B_3,1,B_3,2);             ⋮; (B_m,1,B_m,2) ]}}

where the symmetry of A and B is used to write the left factors
as entries from column 1. Since A_1=B_1, rewrite (<ref>)
as 

    {{[ (B_3,1,A_3,2);             ⋮; (B_m,1,A_m,2) ]}}≠{{[ (B_3,1,B_3,2);             ⋮; (B_m,1,B_m,2) ]}} .

 Now the left hand factors on both sides of (<ref>) match.
Also the column mixes match, colMix(A_3:m,2)=colMix(B_3:m,2),
and the inner products match but the multisets of terms differ. Comparing
the lhs and rhs of (<ref>), one sees that A_3:m,2,
when compared to B_3:m,2, must have distinct symbols exchanged
across the row signature boundaries established in B_1. The symbol
exchange(s) act similar to the roles of α and β in
(<ref>). This is the only way for the multisets of
terms to differ while colMix(A_2) equals colMix(B_2).
By Lemma <ref>, the inner product is unique
to the set of distinct terms and their multiplicities. The exchange
of distinct symbols across the row signature boundaries change the
multiplicities of distinct terms and possibly the set of distinct
terms. So, the inner products associated with the two sides of (<ref>)
cannot be equal, contradicting the hypothesis that (W× W)_i,j=(W× W)_r,s.

Therefore, if colMix([W_i,W_j])=colMix([W_r,W_s])
and SymSqr(W)_i,j≠ SymSqr(W)_r,s,
then (W× W)_i,j≠(W× W)_r,s
as was to be shown.

Widely-spaced primes matrices and regular matrix multiplication can
generate a sequence of refinement that is at least as refined as the
one generating Π(M^*), it just might take twice
as many iterations.

For a widely-spaced primes matrix W,
let W=wspm(W× W) be a widely-spaced
primes matrix generated from the result of squaring W. Then the
pattern from squaring W is a refinement of the pattern
from symbolically squaring W,


    Π(W×W)≼Π(SymSqr(W)).




Let W be a widely-spaced primes matrix and W=wspm(W× W)
be a widely-spaced primes matrix generated from the result of squaring
W.

From Theorem <ref> and Theorem <ref>
one gets that, 

    Π(diag(W×W))≼Π(diag(W× W))=Π(diag(SymSqr(W))).

 So diagonal cells of (W×W)
are a refinement of the diagonal cells of SymSqr(W).
The remaining question is whether the off-diagonal cells of (W×W)
are a refinement of the off-diagonal cells of SymSqr(W).

Assume W_i,j and W_r,s are in the same off-diagonal cell
so W_i,j=W_r,s with (i,j)≠(r,s). 


  * If SymSqr(W)_i,j≠ SymSqr(W)_r,s and
(W× W)_i,j≠(W× W)_r,s, then
done since (W×W)_i,j≠(W×W)_r,s
by Theorem <ref>.

  * If SymSqr(W)_i,j≠ SymSqr(W)_r,s and
(W× W)_i,j=(W× W)_r,s then
consider the column mixes of W_i, W_j, W_r, and W_s.


  * Assume colMix([W_i,W_j])≠ colMix([W_r,W_s]).
Theorem <ref> says W× W and SymSqr(W)
refine the diagonal the same way, so the multisets of diagonal symbols
are not equal {{(W× W)_i,i,(W× W)_j,j}}≠{{(W× W)_r,r,(W× W)_s,s}}.
Then the second squaring results in (W×W)_i,j≠(W×W)_r,s
since the diagonal terms {{W_i,i,W_j,j}}
and {{W_r,r,W_s,s}}
differ in (<ref>) and (<ref>).

  * Assume colMix([W_i,W_j])=colMix([W_r,W_s]).
Theorem <ref>, says the condition where colMix([W_i,W_j])=colMix([W_r,W_s])
and SymSqr(W)_i,j≠ SymSqr(W)_r,s means
(W× W)_i,j≠(W× W)_r,s. This
contradicts the hypothesis that (W× W)_i,j=(W× W)_r,s.
So this is not a valid case.


Therefore, Π(W×W)≼ SymSqr(W)
as was to be shown.

Theorem <ref> implies a sequence of constructing
and squaring widely-spaced matrices will converge to a fixed point
that is a refinement of M^* in at most twice as many iterations
as symbolic squaring. 

A second method generates a sequence of refinement identical to the
symbolic squaring sequence of refinement. However, instead of squaring
a widely-spaced primes matrix, it multiplies two different widely-space
matrices. 

Given a real m× m symmetric
matrix M whose diagonal symbols are distinct from off-diagonal
symbols, there exist widely-spaced primes matrices W1 and W2
such that Π(W1)=Π(W2)=Π(M)
and 

    Π(min((W1× W2),(W1× W2)^T))=Π(SymSqr(M))

 where min() is the element-wise matrix minimum.



Similar to the proof of Theorem <ref>, widely-spaced
primes are used as symbols. Assume M has n distinct symbols
and n_1 are off-diagonal symbols and n_2 are diagonal symbols
with n_1+n_2=n. Let W1=wspm(M) be a widely-space
primes matrix as described in Definition <ref>
where the primes are { p_1,⋯,p_n}, starting
with p_1>m. Then let W2=wspm(M) be a second widely-spaced
primes matrix where the primes are {p_1,⋯,p_n},
starting with p_1>mp_n^2.

Use the first n_1 primes of { p_1,⋯,p_n},
p_1,…,p_n_1, as off-diagonal symbols of W1 and the
next n_2 primes, p_n_1+1,…,p_n, as the diagonal
symbols of W1. Then use the first n_1 primes of {p_1,⋯,p_n},
p_1,…,p_n_1, as off-diagonal symbols
of W2 and the last n_2 primes, p_n_1+1,…,p_n,
as diagonal symbols of W2. W1 and W2 are SPD by construction.

Now consider the matrix product W1× W2. For each (i,j)
location, the inner product is constructed from a row vector of W1
and a column vector from W2. Since the symbol sets for W1 and
W2 are distinct, the factors of each term can be considered to
be ordered. To see this, consider the canonical string at (i,j)
from W1× W2. Each term in an inner product is given by a p_rp_c
where p_r∈{ p_1,…,p_n} and p_c∈{p_1,⋯,p_n}.
Without loss of generality, assume that p_rp_c is
the largest term in the inner product inner product string. Then the
inner product is less than or equal to mp_rp_c.
But the next larger distinct factor after p_rp_c
is greater than mp_r^2p_c. So the inner product
can be decomposed using modulo arithmetic starting with the largest
prime and working downward to determine the exact number of times
each term appears. This implies distinct canonical inner product strings
in SymSqr(M) have distinct inner products in W1× W2.

The remaining part is to note that SymSqr(M) selects
the lessor of the canonical strings at (i,j) and (j,i) to represent
both locations for symmetric matrices. For W1× W2, this is
accomplished by using the minimum of [W1× W2]_i,j
and [W1× W2]_j,i to represent both locations.
Recall that the product of two real symmetric matrices is not symmetric
unless they have the same eigenvectors. So 

    Π(min((W1× W2),(W1× W2)^T))=Π(SymSqr(M))

 as desired.

So the refinement sequence generated by symbolic squaring can be reproduced
using regular matrix multiplication. However, it involves multiplying
widely-spaced primes matrices. 

One consequence of replacing symbolic squaring with regular matrix
multiplication on widely-spaced primes matrices is that the bounds
for Algorithm <ref> reduces from
O(m^8) for symbolic squaring to O(m^6)
for matrix multiplication of widely-spaced primes matrices and from
O(m^10) for symbol substitution of canonical inner
product string arrays to O(m^8) for symbol substitution
of positive integer matrices. Dropping the overall complexity of the
algorithm from O(m^12) to O(m^10). 



 §.§ Symbolic “Walks” 


Two vertices that are not neighbors in a rook's graph are part of
a unique four-cycle <cit.>. This is also true of PCGs since
the connectivity, ignoring edge weights, is identical to a rook's
graph. From Theorem <ref> the first symbolic squaring
has a pattern identical to the pattern from symbolically multiplying
DMD where the m^2× m^2 matrix M=R+D is a PCM. Consider
the value of DMD at (i,j). It is given by 


    [DMD]_i,j=
    D_iiR_ijD_jj   
    
    D_ii^3

which, for i≠ j, may be interpreted as the relationship between
the pair of vertices D_ii and D_jj.

Now symbolically square DMD for the second symbolic squaring. We
have 

    (DMD)(DMD)     =    (D(R+D)D)×(D(R+D)D)
        =     DRD^2RD+DRD^4+D^4RD+D^6.

Focus on the DRD^2RD term, symbolically, 

    Π(SymMult(D,R,D^2,R,D))≼Π(DRD^4+D^4RD+D^6)

 if the zeros on the diagonal of R are treated as symbols. Implying
that 

    Π(SumMult(D,R,D^2,R,D))=Π(SymSqr(DMD)).

 One sees that its value at (i,j), i≠ j, is equal to 

    [DRD^2RD]_ij=∑_k=1^m^2D_iiR_ikD_kk^2R_kjD_jj

 where each term is a “walk” of length two, with segments from
vertex D_ii to vertex D_kk and then from D_kk to vertex
D_jj. Each “walk” is a term in the canonical inner product.
“Walk” is in quotes since a walk normally implies an edge connecting
the vertices. Here “walk” is being used more generally to represent
a string of relationships between vertices.

For the third symbolic squaring, (DMD)^4 has a term
DRD^2RD^2RD^2RD whose symbolic inner product terms for (i,j),
i≠ j, are given by 

    [DRD^2RD^2RD^2RD]_ij=∑_k_1=1^m^2∑_k_2=1^m^2∑_k_3=1^m^2D_iiR_ik_1D_k_1k_1^2R_k_1k_2D_k_2k_2^2R_k_2k_3D_k_3k_3^2R_k_3jD_jj.

Now let i equal to j, (<ref>) contains all of the
four-cycles as well as other “walks” of length four. For the n^ th
symbolic squaring, there are m^2n “walks” of length 2^(n-1)
at each (i,j) location. 

Section <ref> described an unsuccessful attempt
to just use the four-cycles to separate graphs. However, it is also
pointed out that there must be unique row and column signatures that
can be used to pin down the value at a location. As symbolic squaring
continues, the “walks” will incorporate the row and column signatures
as well as the value at the intersection.



 §.§ Practical Considerations


From a practical point of view, Algorithm <ref>
is extremely slow and requires a large amount of memory. A m× m
matrix has a m^2× m^2 PCM which is symbolically dense.
Testing used adjacency matrices of graphs from the Bliss graph collection
<cit.> and the Neuen-Schweitzer graph database <cit.>.
Only matrices up to 180×180 were tested since their PCMs are
32400×32400 and take up about 8GB of memory per PCM when using
64-bit values. To make the execution time tractable, double precision
matrix multiplication with prime numbers as the symbols is used instead
of symbolic matrix multiplication. Symbol substitutions are performed
every iteration. When there are more than 10K symbols the symbol substitution
algorithm switches to using the integers 1,…,k as the symbols
to reduce magnitudes. IEEE 754 double precision format has a 52 bit
mantissa which can represent integers up to 15 digits long. So the
process is halted if an inner product goes over 2^52. Lastly
the heuristic compares column mixes instead of diagonal mixes since
column mixes often detect non p-similar cases an iteration before
the diagonal mixes as discussed in Section <ref>.



§ SUMMARY


Permutation similarity between two (real or complex) square matrices
reduces to asking whether their associated PCGs are isomorphic. PCGs
add symmetric permutation constraints to matrices so multisets can
be used for comparisons. Usually, multisets maintain either row or
column constraints while giving up the other constraint. PCGs make
row and column constraints an integral part of the problem. Then the
associated PCMs are repeatedly symbolically squared and symbols substituted
until stable partitions are reached or the diagonal mixes differ.
If the diagonal mixes do not match, the matrices are not p-similar.
The case where stable partitions are reached with matching diagonal
mixes are p-similar. The necessity and sufficiency of the algorithm
is shown. A second algorithm using the blind algorithm to find the
permutation between two p-similar matrices is also given. Therefore
permutation similarity and graph isomorphism are in P. 

plain


