Making a magic square matrix singular

We know that any magic square matrix of odd order is not singular. When each element of the matrix is subtracted by the sum-average of the total elements, then this perturbed matrix becomes singular, and the determinant of the resulted matrix is zero. That is,

     det(magic(n)-ones(n)*((1+n*n)/2)) = 0,    for any odd n.

Can anyone help me the proof or find literture in this subject?

[Walter Roberson]

For, say, n = 9, then (1+n*n)/2 is 49, and ones(49) is a 49 x 49 matrix... which you try to subtract from the 9 x 9 matrix magic(9).

What is the "sum-average of the total elements" ? Is it (1 + 2 + 3 + ... n)/2 ? Is it (1 + 2 + 3 + ...) / n ? The first of those formulas is n * (n+1) / 2 / 2 or n * (n+1) / 4, and the second of those formulas is n * (n+1) / 2 / n or (n+1)/2 -- neither of which is (1+n*n)/2 . Perhaps you meant sum(sum(magic(n)) / n, which would be (n*n)*(n*n+1)/2/n^2 which would be (n*n+1)/2, the formula you gave? If so then test with det(magic(9)-41) and you will find the result is not 0.

[Feng Cheng Chang]

Dear Walter:

Thank you for your comment.

1. For n=9, (1+n*n)/2 is not 49, but 41.

2. "sum-average of the total elements" means sum(sum(magic(n)))/n, which is (1+n*n)/2.

3. Test with det(magic(n)-ones(n)*(1+n*n)/2) for n=9 will find the result equal to zero.

4. Please try for n = 3,5,7,.... . If it works, I would like to know why.

F C Chang 3/30/2013

[Walter Roberson]

Odd oh odd!!!

R2013a on OS-X:

>> n=9;det(magic(n)-ones(n)*(1+n*n)/2)

ans =

            -0.4816269069579

R2012a on OS-X gives -0.3210846046386

I constructed the matrix and det() it on both sides to be sure that the same matrix was being used, but the det() differ!!! I checked with "which" to be sure I was only using the official MATLAB det.

det([6 17 28 39 -40 -29 -18 -7 4;16 27 38 -32 -30 -19 -8 3 5;26 37 -33 -31 -20 -9 2 13 15;36 -34 -23 -21 -10 1 12 14 25;-35 -24 -22 -11 0 11 22 24 35;-25 -14 -12 -1 10 21 23 34 -36;-15 -13 -2 9 20 31 33 -37 -26;-5 -3 8 19 30 32 -38 -27 -16;-4 7 18 29 40 -39 -28 -17 -6])

[Matt J]

2. "sum-average of the total elements" means sum(sum(magic(n)))/n, which is (1+n*n)/2.

I think you really mean

 sum(sum(magic(n)))/n^2 %divide by n^2

[Feng Cheng Chang]

Dear Matt J:

Thank you again for your reply.

I have rechecked again for n=3,5,7,9,11, and get the results det( )=0 all the time. I use OCTAVE instead of MATLAB. However, I cannot figure out how you get 2 different non-zero results for computing det( ).

The matrix in det([... ]) is what I expected for the 9x9 perturbed matrix. Its determinant shold be equal to zero.

Thabk you for your correction for "sum-average of the total elements".

F C Chang 04/01/2013

[Matt J]

The matrix in det([... ]) is what I expected for the 9x9 perturbed matrix. Its determinant shold be equal to zero.

Feng, the determinant of the perturbed matrix is zero. Nobody disputes that and in fact, you've been given analytical proofs of why the matrix is singular.

The point, though, is that determinant computations are sensitive to floating point errors. I don't know how OCTAVE contends with that. What happens in OCTAVE when you compute the determinant as prod(eig(P))? In MATLAB, you get a very bad result,

>> P=magic(9)-41; prod(eig(P))

ans =

   -2.2875

However, MATLAB's rank() command clearly knows that P is singular

>> rank(P)

ans =

       8

[Feng Cheng Chang]

Dear Matt J,

I checked the results using OCTAVE, and the results give prod(eig(P))= -3.7867 (but not -2.2875) and rank(P)=8. But what for?

F C 04/01/2013

[Matt J]

But what for?

Well, it shows you that there are bad ways to compute determinants, even for integer matrices. The formula det(P)=prod(eig(P)) clearly doesn't work here very well, again because of precision issues.

Answer 1

[Ahmed A. Selman]

I don't think details are required since

A=magic(n)-ones(n)*((1+n*n)/2)

is changed into an antisymmetric matrix, any such A matrix must satisfy (basic math.. etc)

det(A) = -1^n * det(A)

since n is odd, det(A) must be zero (thus, A is singular). Changing A from magic(n) to (magic(n)-ones(n)*((1+n*n)/2) ) as mentioned in the question is enough to destroy the symmetry of A.

Yet, since this is too basic, and it works the same for magic(n) with n is odd or even, (also, produces antisymmetric), I'm afraid you already know this. I tried (quickly, to be honest) other means like the nice arguments above, but didn't got anything useful so I thought to share, it might help. Regards.

[Feng Cheng Chang]

Dear Ahmed:

Thank you very much for your reply.

Any n order magic square matrix perturbed by -ones(n)*(1+n*n)/2 is always an antisymmetric matrix. and it also becomes singular with zero determinant (within machine precession). The order n can be either odd or even.

This is the perfect and correct answer what I want. Thank you again for your answer. I would be very much pleased to accept your answer.

By the way, any non-singular matrix perturbed by the matrix D will become singular. For example,

 n=5; A=randn(n), D=ones(n)/sum(sum(inv(A)), P=A-D, ch=det(P)

It will give the result ch=0 (within machine precession).

Please try some matrices, such as pascal(4), magic(7), etc. and any type-in square matrices.

F C Chang 04/01/2013

[Matt J]

I'm doubtful about the claim that the perturbed matrix is anti-symmetric. An anti-symmetric matrix should satisfy A+A.'=0, which examples will show is not true,

>> A=magic(3)-5; A+A.'

ans =

       6    -6     0
      -6     0     6
       0     6    -6

[Feng Cheng Chang]

I think that the anti-symmetry means A(i,j)=-A(n+1-i,n+1-j), and the skew-symmetry means A(i,j)=-A(j,i). FC 4/1/13

[Matt J]

If you say so. I can't find that distinction either at Wolfram or Wikipedia. It's also not clear to me how

 det(A) = -1^n * det(A)

still applies when A(i,j)=-A(n+1-i,n+1-j). It clearly applies to skew-symmetric matrices.

[Ahmed A. Selman]

If (A) was a square, i.e., (n by n), and antisymmetric matrix, with (n) is odd, the relation is valid. How is that basic math? Please be patient to read the following:

The definition of a determinant of any square matrix A (size is n by n, for any integer n > 0 ) is

   det(A)= S( (i=1 to n) a(i,j) F(i,j) );

Here, (S) refers to summation from (i=1 to n); a(i,j) is the matrix element of (A), and F(i,j) is the co-factor of (A), found from

   F(i,j)=(-1)^(i+j)*Mij

and (Mij) being the minor of (A).

Now, let (b) be any real constant. By definition,

   det(bA)=b^n *det(A).

This comes true if we substituted (b) in the basic, first definition of determinant of (A).

Now, a matrix (A) is called (symmetric) iff (if and only if) the condition was satisfied:

   transpose (A) = A,

or, by definition of transpose (I don't believe proof is needed), then:

   Aij=Aji

But if

   Aij= - Aji

the matrix (A) is said to be (antisymmetric, or skew-symmetric, the same are equal last time I checked). Now, for any odd (n) value, an (n by n) square, antisymmetric matrix, we can and will use

b=-1 , so, from

   det(bA)= b^n det(A) 

   det(-A)=(-1)^n det(A),

the L.H.S. for such equality is the same as det(A) - again from the definition of determinant and transpose...Thus

    LHS= det(-A)=det(A)

but the R.H.S. of the last equality was

   RHS= (-1)^n det(A);

equaling RHS with LHS,

   det(A)=(-1)^n det(A)

which can be true iff

   det(A)=0.

ANd for any matrix (A), if the determinant was (ZERO) then (A) is said to be (singular) matrix.

Thank you, Feng and Matt, for all the earlier, useful comments.

[Ahmed A. Selman]

And this basic, primitive derivation, is found (must be found) in any textbook dealing with matrices and determinant properties. I did find it on Wolfram search that:

   det(-A)=(-1)^n det(A)

from: <http://mathworld.wolfram.com/Determinant.html >

and found that for antisymmetric matrix A then

   Aij= - Aji

from

< http://mathworld.wolfram.com/AntisymmetricMatrix.html>

The rest is, however, a plane and direct substitution.

[Matt J]

Ahmed,

We established several Comments ago that Aij=-Aji is not satisfied for the modified magic square matrix.

There may be a way to extend the determinant equation to the weirder kind of asymmetry that this matrix exhibits, but it looks like it would take some work. Showing that ones(n,1) is a null-vector of the matrix seems to me like the quicker proof, not to mention that it also covers even-valued n.

Answer 2

[Jonathan Epperl]

The row-sum, column-sum and diag-sum of a magic square are all the same, and the magic square uses all the integeres 1:n^2. Thus, the sum of all elements must be n^2*(n^2+1)/2, and each row, column, diag sum must be n*(n^2+1)/2.

Now look at what you wrote, multiply it from the right by ones(n,1), and you'll see that you will get zero. Voila, thus the matrix is singular.

[Feng Cheng Chang]

Dear Jonathan Epperl:

Thank you for your reply.

The first paragragh is the well-known definition or description of magic square matrix.

The second paragragh is unclear to me how to get the matrix singular. The test statement is listed here again for your reference

   det(magic(n)-ones(n)*(1+n*n)/2), for any odd n.

So I cannot accept your answer at this time.

F C Chang 3/30/2013

[Jonathan Epperl]

Just to fill that hole in your knowledge: A square Matrix A is singular if and only if

• inv(A) does NOT exist
• det(A)==0
• The range of A is not all of R^n
• THE KERNEL OF A IS NONTRIVIAL
• ...

That last point there means that if you can find a nonzero vector v such that A*v==0, then you have proven that A is singular, and ones(n,1) is such a vector.

Answer 3

[Matt J]

Let x=ones(n,1)/n and P be the perturbed matrix. You can verify that

P*x=0

proving that P is singular.

[Feng Cheng Chang]

Dear Matt J:

Thank you for your answer.

Let n=5, x=ones(n,1)/n, and P=magic(n)-ones(n)*(1+n*n)/2, then we get det(P)=0. Thus P is singular. But we can not get P*x=0 proving P is singular as you stated in your answer.

Sorry I cannot accept your answer.

F C Chang 3/30/2012

[Matt J]

Feng, bear in mind that finite machine precision is in play here. You will not get P*x=0 exactly, but you will get very close

    n=5;
    x=ones(n,1)/n;
    P=magic(n)-ones(n)*(1+n*n)/2;

    >> norm(P*x)

    ans =

       4.5776e-16

Incidentally, for the same reasons also you will NOT get det(P)=0 exactly

    >> isequal(det(P),0)

    ans =

         0

And, as Walter has shown you, the calculation of det(P) is numerically unstable. You can get highly nonzero-looking results for larger n on some machines. I also get similar results to his

 >> det(magic(9)-41)

ans =

   -0.3211

Here's more of an analytical proof. You know that the row-sums of all magic squares are the same. If S is the "sum-average", then the row-sums in terms of S are all S/n. Therefore

 P*x= mean(magic(n) - ones(n)*S/n , 2)
    = sum(magic(n),2)/n - sum(ones(n,1)*S/n , 2)/n
    = ones(n,1)*S/n - ones(n,1)*S/n
    = zeros(n,1)

Voila! And by the way, this is clearly true for all n, both even and odd.

[Feng Cheng Chang]

Dear Matt J,

Thank you again for your reply.

I do not think the calculation of det(P) is numerically unstable since elements of P in the case are all integers.

By the way I got the result det(magic(9)-41)=0, instead of =-0.321.

F C Chang 4/1/2013

[Matt J]

I do not think the calculation of det(P) is numerically unstable since elements of P in the case are all integers.

You don't know that MATLAB computes determinants in steps that always result in integers. In fact, I've just learned that MATLAB uses LU factorization, which is consistent with the results Walter and I have been getting,

>> P=magic(9)-41; [L,U]=lu(P); det(L)*det(U)

ans =

     -0.3211

[Feng Cheng Chang]

Dear Matt J,

The result I got (by using OCTAVE within machine precession) is 0, but not -0.3211. I hope someone can verify it.

F C 04/01/2013

Answer 4

[Feng Cheng Chang]

Dear Readers:

I have already accepted the expected simple answer from Ahmed A Selman, so I would like to conclude this "Question-Answer" contribution.

Also I would like to thank Walter Roberson, Matt J and Jonathan Epperl for the further extensive discussion related to the subject.

It gives me the great opportunity to find the interesting result unexpectedly. Hope anyone can see it.

For any n x n non-singular square matrix A, I can make it singular by perturbing the same amount to each of the matrix elements. That is

    det(A-ones(n)/sum(sum(inv(A))) == 0

I would like to extend "the making a matrix singular" further to include the perturbation of any single element, any row or coulumn elements, or even the entire matrix elements in the prescribed distribution. Please watch my contribution in MATLAB Files Exchange in the near future.

For now, please try the followings:

    for n=1:9, A=rand(n); P=A-1/sum(sum(inv(A))); detP=det(P); n,A,P,detP, end;

Best wishes,

Feng Cheng Chang 04/06/2013