Why do siamese magic squares have real eigenvalues, symmetric around zero?

Question

There is a standard method to construct magic squares of odd size, known as the Siamese construction. I'll write $S_m$ for the $m \times m$ Siamese square. For example, here is $S_5$.

17    24     1     8    15
23     5     7    14    16
 4     6    13    20    22
10    12    19    21     3
11    18    25     2     9

As you can see, this matrix certainly isn't symmetric! It has one obvious eigenvector, the all ones vector, whose corresponding eigenvalue is $\frac{m^3+m}{2}$, and there is no obvious reason its other eigenvalues should have nice structure. Nonetheless, experimentation suggests:

(1) All eigenvalues of $S_m$ are real.

(2) If $\lambda$ is an eigenvalue of $S_m$ other than $\frac{m^3+m}{2}$, then $-\lambda$ is as well.

Can anyone explain why?

Remark In case anyone is curious how I found this, I am teaching myself MATLAB and came up with the task of listing characteristic polynomials of magic squares as a more or less random task which would involve learning some basic control structures and some mathematical tools. This piece of MATLAB documentation describes some alternate ways of thinking about the Siamese construction, which may be relevant.

Answer 1

Partial answer (explaining but not proving why the non-trivial eigenvalues come in pairs with opposite signs).

It should not be too hard to prove that construction yields a matrix with the property that rotating it by 180 degrees gives "an entrywise complementary magic square", i.e. the sum of $A=S_m$ and its rotated copy $\tilde{A}$ is the matrix with all entries equal to $m^2+1$.

In other words $$ A_{i,j}=m^2+1-A_{m+1-i,m+1-j} $$ for all $i,j$. Let $B=(A-\tilde{A})/2$. It follows that the row and column sums of $B$ are all zero, and that $$B_{i,j}=-B_{m+1-i,m+1-j}\qquad(*)$$ for all $i,j$. When $m=5$ we have $$ B=\left( \begin{array}{ccccc} 4 & 11 & -12 & -5 & 2 \\ 10 & -8 & -6 & 1 & 3 \\ -9 & -7 & 0 & 7 & 9 \\ -3 & -1 & 6 & 8 & -10 \\ -2 & 5 & 12 & -11 & -4 \\ \end{array} \right). $$

If $u=(x_1,\ldots,x_m)$ is an eigenvector of $B$ belonging to eigenvalue $\lambda$, then it is easy to show that $\tilde{u}=(x_m,x_{m-1},\ldots,x_1)$ is an eigenvector of $B$ belonging to eigenvalue $-\lambda$. Namely the assumption is that for all $i$ we have $$ \sum_j B_{i,j}x_j=\lambda x_i. $$ Taking $(*)$ into account this implies that $$ \sum_j B_{i,j}x_{m+1-j}=-\sum_j B_{m+1-i,m+1-j}x_{m+1-j}=-\lambda x_{m+1-i} $$ as required.

Clearly $B$ commutes with the all ones matrix $\mathbf{1}$ (the products are zero in either direction, because row and column sums of $B$ vanish). So if $B$ is diagonalizable, then it is simultaneously diagonalizable with $\mathbf{1}$. We can take advantage of the known eigenvalues of $\mathbf{1}$: the all ones vector spans the 1-dimensional eigenspace with $\lambda=m$, and its orthogonal complement $V$ (the zero sum subspace of $\Bbb{R}^m$) is the $(m-1)$-dimensional eigenspace belonging to $\lambda=0$.

The all ones vector belongs to eigenvalue $0$ of $B$, so all the eigenvectors of $B$ belonging to non-zero eigenvalues are in $V$. This is important because $$ A=\frac{m^2+1}2\mathbf{1}+B. $$ If $u$ belongs to eigenvalue $\lambda\neq0$ of $B$, then $\mathbf{1}u=0$ and thus $Au=\lambda u$. So the non-zero eigenvalues of $B$ are also eigenvalues of $A$ and the claim follows.

Missing:

1. Why are the eigenvalues of $B$ all real and distinct (so that $0$ has multiplicity one, and that $B$ is diagonalizable)?
2. Prove that the Siamese construction implies $(*)$.

Answer 2

Here are some partial results.

This paper states that every regular magic square $\mathbf A\in\mathbb N^{n\times n}$ (that is, its entries satisfy $a_{i,j}+a_{n-i+1,n-j+1}=n^2+1$) is representable as a rank-$1$ correction of a skew-centrosymmetric matrix $\mathbf Z$ (that is, $\mathbf Z=-\mathbf J\mathbf Z\mathbf J$, where $\mathbf J$ is the exchange matrix, or the identity matrix with its columns reversed). In other words, the matrix you get when you subtract $\frac{n^2+1}{2}$ from all of the entries of $\mathbf A$ is skew-centrosymmetric. This is important since the eigenvalues of $\mathbf A$, apart from the magic eigenvalue $\frac{n^3+n}{2}$, are the same as the eigenvalues of $\mathbf Z$, which now has $0$ as an eigenvalue.

We can now use the results of this paper, which says that if $\lambda$ is an eigenvalue of $\mathbf Z$, then $-\bar{\lambda}$ must also be an eigenvalue, with the same algebraic multiplicity.

Thus, the remaining task is to prove that all the eigenvalues of a Siamese magic square are necessarily real.

---

A tantalizing result: Pressman, in this paper, shows the stronger result that if $\lambda$ is an eigenvalue of $\mathbf Z$, then $-\lambda$ is also an eigenvalue. There seems to be a criterion in the paper for saying when a skew-centrosymmetric matrix has all its eigenvalues real, but I have not yet figured how to apply it to this case.

Answer 3

This solution stems from understanding the constructions of the magic square algorithm, I will discuss only for odd case, it is a heuristic approach, consider:

n = 7;    % Must be odd
[I,J] = ndgrid(1:n);
A = mod(I+J+(n-3)/2,n);
B = mod(I+2*J-2,n);
M = n*A + B + 1
magic(n) 

If you see the matrices that are being constructed , $A$ and $B$ which inturn generate the magic square $M$m they have a very special structure. For instance when $n= 7$, this is their structure,

A =

 4     5     6     0     1     2     3
 5     6     0     1     2     3     4
 6     0     1     2     3     4     5
 0     1     2     3     4     5     6
 1     2     3     4     5     6     0
 2     3     4     5     6     0     1
 3     4     5     6     0     1     2

B =

 1     3     5     0     2     4     6
 2     4     6     1     3     5     0
 3     5     0     2     4     6     1
 4     6     1     3     5     0     2
 5     0     2     4     6     1     3
 6     1     3     5     0     2     4
 0     2     4     6     1     3     5

One can check that A is symmetric indefinite matrix, a matrix is symmetric indefinite if it is symmetric and has both positive and negative eigenvalues. Symmetric indefinite matrices are an important class of matrices arising in many applications. Observe also that this is due to the fact that A is being constructed by addition of replicated rows and columns with numbers 1 to n ( $I$ and $J$ ) plus some offset (see. symmertic indefinite solutions).

Now consider the sum, $M = nA + (B + 1)$, while $(B+1)$ is not symmetric indefinite, the first term begins to dominate as A is added successively e.g.

eig(3*A + B + 1) =  
  91.0000 + 0.0000i
  24.2980 + 0.0000i
  12.1490 + 1.3399i
  12.1490 - 1.3399i
 -24.2980 + 0.0000i
 -12.1490 + 1.3399i
 -12.1490 - 1.3399i

eig(4*A + B + 1) = 
   112.0000
   32.3220
   16.8439
   15.4781
  -32.3220
  -16.8439
  -15.4781

eig(M) =
  175.0000
  -56.4848
  -31.0882
  -25.3967
   56.4848
   25.3967
   31.0882

making $M \approx nA$ and its eigenvalue structure begins to resemble that of $A$ which is a symmetric indefinite matrix with non-dominant diagonal and hence M will have both positive and negative real eigenvalues.

We can also try to use other approaches, to prove that the eigenvalues occur in pair of positive and negative values, one can proceed as follow: using a power iteration method or something similar we can find the largest eigenvalue say $\lambda_1$, using trace of the matrix, $M$ it can then be shown that Trace($M$) = $\lambda_1$, hence owing to the symmetry, the other eigenvalues must occur in pairs of positive and negative values cancelling all but one which is the largest. To gain more insight you can also use LDL decomposition which are used for solving equations of indefinite linear systems. Hope this will help to develop a rigorous mathematical solution.

Answer 4

just doing the matrices as in the writeup at David's second link by one Cleve Moler. It would not let me do $n=105,$ but I did get $n=35.$

They give some matrices, $A, B, 1, M.$ Here the numeral $1$ refers to the matrix with all entries $1,$ I called that $C$ below. $$ 0 \leq a_{ij} \leq n-1, \; \; \; a_{ij} \equiv i + j + \frac{n-3}{2} \pmod n, $$ $$ 0 \leq b_{ij} \leq n-1, \; \; \; b_{ij} \equiv i + 2j -2 \pmod n, $$ $$ c_{ij} = 1, $$ $$ M = nA + B + C. $$

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

3:

parisize = 4000000, primelimit = 500509
?  a = [   2,    0,    1;    0,    1,    2;    1,    2,    0]
%1 = 
[2 0 1]

[0 1 2]

[1 2 0]

? b = [    1,    0,    2;    2,    1,    0;    0,    2,    1]
%2 = 
[1 0 2]

[2 1 0]

[0 2 1]

? c = [    1,    1,    1;    1,    1,    1;    1,    1,    1]
%3 = 
[1 1 1]

[1 1 1]

[1 1 1]

? m = 3 * a + b + c
%4 = 
[8 1 6]

[3 5 7]

[4 9 2]

? factor(charpoly(a))
%5 = 
[x - 3 1]

[x^2 - 3 1]

? factor(charpoly(b))
%6 = 
[x - 3 1]

[x^2 + 3 1]

? factor(charpoly(m))
%7 = 
[x - 15 1]

[x^2 - 24 1]

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

    5:
    a:
[3 4 0 1 2]

[4 0 1 2 3]

[0 1 2 3 4]

[1 2 3 4 0]

[2 3 4 0 1]


b: 
[1 3 0 2 4]

[2 4 1 3 0]

[3 0 2 4 1]

[4 1 3 0 2]

[0 2 4 1 3]

     m = 5 * a + b + c
    %4 = 
    [17 24 1 8 15]

    [23 5 7 14 16]

    [4 6 13 20 22]

    [10 12 19 21 3]

    [11 18 25 2 9]

    ? factor(charpoly(a))
    %5 = 
    [x - 10 1]

    [x^4 - 25*x^2 + 125 1]

    ? factor(charpoly(b))
    %6 = 
    [x - 10 1]

    [x^4 - 125 1]

    ? factor(charpoly(m))
    %7 = 
    [x - 65 1]

    [x^4 - 625*x^2 + 78000 1]

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

7:
? factor(charpoly(a))
%5 = 
[x - 21 1]

[x^6 - 98*x^4 + 2401*x^2 - 16807 1]

? factor(charpoly(b))
%6 = 
[x - 21 1]

[x^6 - 16807 1]

? factor(charpoly(m))
%7 = 
[x - 175 1]

[x^6 - 4802*x^4 + 5764801*x^2 - 1988873152 1]

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

9:

? factor(charpoly(a) ) 
%7 = 
[x - 36 1]

[x^2 - 27 1]

[x^6 - 243*x^4 + 13122*x^2 - 177147 1]

? factor(charpoly(b) ) 
%8 = 
[x - 36 1]

[x^2 + 27 1]

[x^6 + 177147 1]

? factor(charpoly(m) ) 
%9 = 
[x - 369 1]

[x^2 - 2160 1]

[x^6 - 19683*x^4 + 86093442*x^2 - 94143001680 1]

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

11:

? factor(charpoly(a))
%5 = 
[x - 55 1]

[x^10 - 605*x^8 + 102487*x^6 - 7086244*x^4 + 214358881*x^2 - 2357947691 1]

? factor(charpoly(b))
%6 = 
[x - 55 1]

[x^10 + 2357947691 1]

? factor(charpoly(m))
%7 = 
[x - 671 1]

[x^10 - 73205*x^8 + 1500512167*x^6 - 12553713506884*x^4 + 45949729863572161*x^2 - 61159090446056598600 1]

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

13:

? 
? factor(charpoly(a))
%5 = 
[x - 78 1]

[x^12 - 1183*x^10 + 399854*x^8 - 57921708*x^6 + 4078653605*x^4 - 137858491849*x^2 + 1792160394037 1]

? factor(charpoly(b))
%6 = 
[x - 78 1]

[x^12 - 1792160394037 1]

? factor(charpoly(m))
%7 = 
[x - 1105 1]

[x^12 - 199927*x^10 + 11420230094*x^8 - 279577021469772*x^6 + 3327083045915899205*x^4 - 19004963774880799438801*x^2 + 41753905413411324206651760 1]

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

15:

? factor(charpoly(a))
%5 = 
[x - 105 1]

[x^2 - 75 1]

[x^4 - 225*x^2 + 10125 1]

[x^8 - 1800*x^6 + 708750*x^4 - 79734375*x^2 + 2562890625 1]

? factor(charpoly(b))
%6 = 
[x - 105 1]

[x^2 + 75 1]

[x^4 - 10125 1]

[x^4 + 50625 2]

? factor(charpoly(m))
%7 = 
[x - 1695 1]

[x^2 - 16800 1]

[x^4 - 50625*x^2 + 512568000 1]

[x^8 - 405000*x^6 + 35880570000*x^4 - 908244868359375*x^2 + 6568148865600000000 1]

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

21:
? factor( charpoly(a))
%9 = 
[x - 210 1]

[x^2 - 147 1]

[x^6 - 882*x^4 + 194481*x^2 - 12252303 1]

[x^12 - 7056*x^10 + 11668860*x^8 - 6689757438*x^6 + 1664205811884*x^4 - 183478690760211*x^2 + 7355827511386641 1]

? factor( charpoly(b))
%10 = 
[x - 210 1]

[x - 21 2]

[x + 21 2]

[x^2 - 21*x + 441 2]

[x^2 + 147 1]

[x^2 + 21*x + 441 2]

[x^6 - 12252303 1]

? factor( charpoly(m))
%11 = 
[x - 4641 1]

[x^2 - 64680 1]

[x^6 - 388962*x^4 + 37822859361*x^2 - 1051059451035132 1]

[x^12 - 3111696*x^10 + 2269371561660*x^8 - 573754546059690240*x^6 + 62945022637525450097340*x^4 - 3060402710940787304923125687*x^2 + 54108197115511006692907045070400 1]

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

25:


? factor(charpoly(a))
%7 = 
[x - 300 1]

[x^4 - 625*x^2 + 78125 1]

[x^20 - 15625*x^18 + 68359375*x^16 - 130615234375*x^14 + 131225585937500*x^12 - 76389312744140625*x^10 + 27120113372802734375*x^8 - 5960464477539062500000*x^6 + 791624188423156738281250*x^4 - 58207660913467407226562500*x^2 + 1818989403545856475830078125 1]

? factor(charpoly(b))
%8 = 
[x - 300 1]

[x^4 - 78125 1]

[x^20 - 1818989403545856475830078125 1]

? factor(charpoly(m))
%9 = 
[x - 7825 1]

[x^4 - 390625*x^2 + 30517500000 1]

[x^20 - 9765625*x^18 + 26702880859375*x^16 - 31888484954833984375*x^14 + 20023435354232788085937500*x^12 - 7285052561201155185699462890625*x^10 + 1616484723854227922856807708740234375*x^8 - 222044604925031308084726333618164062500000*x^6 + 18431436932253575378126697614789009094238281250*x^4 - 847032947254300339068322500679641962051391601562500*x^2 + 16543612251060553497428173839580267667770385742187500000 1]

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

35:

? factor( charpoly(a))
%5 = 
[x - 595 1]

[x^4 - 1225*x^2 + 300125 1]

[x^6 - 2450*x^4 + 1500625*x^2 - 262609375 1]

[x^24 - 58800*x^22 + 942392500*x^20 - 6430253156250*x^18 + 22590813918750000*x^16 - 45546375353896484375*x^14 + 56625598053505126953125*x^12 - 45323879544822393798828125*x^10 + 23792863499842512817382812500*x^8 - 8143807322924036556243896484375*x^6 + 1750204205365253470420837402343750*x^4 - 214400015157243550126552581787109375*x^2 + 11419131242070580387175083160400390625 1]

? factor( charpoly(b))
%6 = 
[x - 595 1]

[x^4 - 300125 1]

[x^4 + 1500625 2]

[x^6 - 262609375 1]

[x^8 - 1500625*x^4 + 2251875390625 2]

? factor( charpoly(m))
%7 = 
[x - 21455 1]

[x^4 - 1500625*x^2 + 450374778000 1]

[x^6 - 3001250*x^4 + 2251875390625*x^2 - 482768305881750000 1]

[x^24 - 72030000*x^22 + 1414177745312500*x^20 - 11820513337182128906250*x^18 + 50871697917821843261718750000*x^16 - 125641833194720435000514984130859375*x^14 + 191350382223376715547899626959845458984375*x^12 - 187620244496627071229425371028983150177001953125*x^10 + 120652249258767712686244839929865230231475830078125000*x^8 - 50588556607865112913542176134821560108671009540557861328125*x^6 + 13318325045557471063558895256155938430252362415194511413574218750*x^4 - 1998581152148968001166733191860870554971460885345004498958587646484375*x^2 + 130396558323632395895356353118015771568144032806708128677368164062500000000 1]

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=