I try to find a general vector base for all magic squares of dimension n . Why does this program work for n = 3, but not vor n > 3?

The program below is supposed to find a general vector base for all magic squares of dimension n . Why does this program work for n = 3, but not vor n > 3? Any help greatly appreciated - I am afraid I am making a very stupid mistake here. Maybe I did not properly understand the usage of rref ?

Brief introduction: a magic square is an n x n matrix, where the sums in all lines, rows and the two diagonals are the same. Translating this into an quation system Mx = const you get 2*n+2 equations for n^2 variables. The solution of this (underdetermined) system is a vector space spanning all possible magic squares of dimension n. The solution of Mx = const should appear in Mind, the base matrices of the magic-square-vector-space in Mbase

function Mbase = MagicEq(n)

% identify the vector base of all magic squares of length n .

% The return values of this function are an array of matrices, M1, M2, .... where all linear combinations a*M1 + b*M2 + .... form magic squares

M = zeros(2*n+2,n^2); % initialize a matrix holding n^2 variables and 2*n + 2 constraints of the magic square of the form Mx = const

for k = 1:n

M(k,1+n*(k-1):n+n*(k-1)) = 1; % constraints on row sums

end

for k = 1:n

M(n+k,k+(0:n-1)*n) = 1; % constraints on column sums

end

M(2*n+1,(0:n-1)*n + (1:n)) = 1; % constraint for nw-se diagonal

M(2*n+2,(1:n)*n - (0:n-1)) = 1; % constraint for ne-sw diagonal

M = rref(M); % solve the equation system by computing the rref form of M

Mold = M; % just for debugging

r = rank(M);

Mind = [];

k = 1;

while ( k <= size(M,2)) % sort out all columns which are non-unit vectors

if sum(abs(M(:,k))) > 1

Mind = [Mind, -M(:,k)]; % and move them to the right side of Mx = const

M(:,k) = [];

else

k = k + 1;

end

end

Mind = [Mind(1:r,:);eye(n^2-r)]; % variables from row r+1:end are free. So add a unit matrix

for k = 1:(n^2-r)

Mbase(:,:,k) = reshape(Mind(:,k),n,n); % 'un-vectorize' the right side of the matrix equation to create the base matrices of the magic square vector space

end

Mbase(:,:,n^2-r+1) = ones(n,n); % and add the trivial magic square

end

[John D'Errico]

Question. If you can use rank in your code, then why would you not just use a tool like orth? Since both rank AND orth implicitly use svd, if you can use one, then surely you can use the other?

Next, your code does NOT return an array of matrices! It will return an array, where you may choose to interpret the rows or columns as vectors.

[Thomas]

I could use orth or svd, but how could this help me?

If I stop the execution of the function above in the line "Mind = [Mind(1:r,:);eye(n^2-r)];", then

Mind =

1 0 0 1 1 1 1

1 -1 1 1 0 1 2

-1 1 -1 -1 -1 -2 -2

-1 0 0 -1 0 0 -1

-1 1 1 0 0 0 0

-1 0 -1 -1 -1 -1 -2

1 -1 0 1 1 1 2

0 -1 -1 -1 0 0 0

0 0 0 0 -1 -1 -1

0 0 0 0 0 0 0

which is already wrong at this point. The following "un-vectorization" of the first column of Mind gives me a first element of the magic-square-vector-base that would be:

Mbase1 =

1 1 -1 -1

-1 -1 1 0

0 0 NaN NaN

NaN NaN NaN NaN

which is already wrong, since the sum of line 2 is not equal the sum of line 1.

A very similar problem would occur with column 2 of Mind, etc.

Answer 1

[John D'Errico]

No resonse yet? Again, I don't see why you are using rref here. If all you want are basis vectors, then tools like null and orth will do it for you.

But looking at your code, you are implicitly unrolling the nxn array into a vector of length n^2. Then you form a matrix M that defines what you know about those elements, in terms of a magic square.

First, define M. No need for loops. Kron can be a very useful tool for such problems.

n = 3;
M = kron(ones(1,n),eye(n)); % row sums
M = [M;kron(eye(n),ones(1,n))]; % column sums
M(2*n+1,(0:n-1)*n + (1:n)) = 1;           % constraint for nw-se diagonal
M(2*n+2,(1:n)*n - (0:n-1)) = 1;           % constraint for ne-sw diagonal

M

M = 8×9

1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 1 0 0 0 1 0 0 1 0 1 0 1 0 0

So M is an 8x9 array here, for the 3x3 case. It should have rank 7 for the 3x3 magic square. Essentially, that means there are 2 elements we could arbitrarily set.

rank(M)

ans = 7

Now, We will first find a specific solution to the magic square of size 3x3. We know the sums will be 15 for the basic case. backslash will suffice. Ignore the warning, because we already know this matrix is singular. What, me worry?

S0 = M\repmat(15,8,1)

Warning: Rank deficient, rank = 7, tol = 3.996803e-15.

S0 = 9×1

5.0000 0 10.0000 10.0000 5.0000 0 -0.0000 10.0000 5.0000

Again, these are integers. The solve introduces crap down at the level of the least significnt bits.

S0 = round(S0)

S0 = 9×1

5 0 10 10 5 0 0 10 5

reshape(round(S0),[3,3])

ans = 3×3

5 10 0 0 5 10 10 0 5

And that clearly is a magic square. All row, and column sums, as well as the diagonals are 15. It is a pretty boring one though.

Next, the symbolic version of null is very nice, in that it will return integer solutions for the nullspace vectors for an integer array like this.

Mnull = null(sym(M))

Mnull = 

This tells us that the set of all solutions will be of the general form

u = sym('u',[2 1]);
allsquares = S0 + Mnull*u

allsquares = 

Try it out. For example, when u1=u2=-3, we get this magic square:

reshape(S0 + Mnull*[-3;-3],[3,3])

ans = 

And indeed that is the same magic square as generated by magic.

magic(3)

ans = 3×3

8 1 6 3 5 7 4 9 2

Other linear combinations of the basis vectors in Mnull will generate correspondingly different magic squares (when then added to S0.)

Be careful though, as there is no magic square of size 2x2.

So I think you are trying to build a relation like that which I built in allsquares, but it is not clear. Note that for larger magic squares, the problem gets more complicated, since you always have n^2 unknowns, but then you will always have only 2*n+2 linear equality constraints on the problem. And of course for larger values of n, n^2 grows much faster than 2*n+2. The disparity between the two gives the number of degrees of freedom.

n = (3:6)';
table(n,n.^2,2*n+2)

ans = 4×3 table
    n    Var2    Var3
    _    ____    ____

    3      9       8 
    4     16      10 
    5     25      12 
    6     36      14 

And that would tell you there is a large space of possible solutions for a 6x6 magic square, so a 22 dimension subspace of solutions, most of which are not integer valued. A cute trick would be to use a tool like intlinprog to find an integer solution. For example...

n = 6;
M = kron(ones(1,n),eye(n));
M = [M;kron(eye(n),ones(1,n))];
M(2*n+1,(0:n-1)*n + (1:n)) = 1;
M(2*n+2,(1:n)*n - (0:n-1)) = 1;

f = ones(n^2,1);
intcon = 1:n^2;
Sq = intlinprog(f,intcon,[],[],M,ones(2*n+2,1))

LP:                Optimal objective value is 6.000000.                                             

Heuristics:        Found 1 solution using total rounding.                                           
                   Upper bound is 6.000000.                                                         
                   Relative gap is 0.00%.                                                          


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0. The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05.

Sq = 36×1

0 2 0 0 -1 0 0 -1 0 0

reshape(Sq,[6,6])

ans = 6×6

0 0 2 -1 0 0 2 -1 0 0 0 0 0 0 0 0 1 0 0 0 -1 1 1 0 -1 2 0 0 0 0 0 0 0 1 -1 1

This is a magic square, in the sense that the sum of each row, column, and diagonal are all equal to 1, and all are integers, though they are not distinct integers. But, having found one solution, all other solutions can be found by adding some linear combination of the columns of M.

[Thomas]

Hi John,

thank you very much for your input - this is a great answer, that helps me a lot!

Thomas