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?
3 ビュー (過去 30 日間)
古いコメントを表示
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
2 件のコメント
John D'Errico
2023 年 6 月 21 日
編集済み: John D'Errico
2023 年 6 月 21 日
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.
採用された回答
John D'Errico
2023 年 6 月 22 日
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
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)
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)
Again, these are integers. The solve introduces crap down at the level of the least significnt bits.
S0 = round(S0)
reshape(round(S0),[3,3])
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))
This tells us that the set of all solutions will be of the general form
u = sym('u',[2 1]);
allsquares = S0 + Mnull*u
Try it out. For example, when u1=u2=-3, we get this magic square:
reshape(S0 + Mnull*[-3;-3],[3,3])
And indeed that is the same magic square as generated by magic.
magic(3)
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)
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))
reshape(Sq,[6,6])
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.
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Logical についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!