Permutation 法の結果，2x2, 4x4 の場合を見ると，この matrix は orthogonal な row を持っているようである．geometric に考えてみれば，determinant は volume である．特定の座標しか使えないのであれば，軸が直交するものの体積が最大というのは直感的である．
そこで，row ごとに orthogonal な vector を見るという手法をimplement したのが program 2 である．
Program 2
function algo_02(mat_rank)
% Introduction to linear algebra Chapter 5. problem 33
% Algorithm 2: generate Hadamard matrix (each row is orthogonal), but
% this only can gives me 1,2,4k matrices
% @author Hitoshi
if nargin ~= 1;
error('Usage: algo_02(mat_rank).')
end
% possible element set A_i = {-1, 1}
SetA = [-1 1];
cand_mat = zeros(mat_rank, mat_rank);
cand_mat(1,:) = ones(1, mat_rank);
cand_row = zeros(1, mat_rank);
global MAXDET
global MAXDET_MAT
MAXDET = 0;
MAXDET_MAT = zeros(1, mat_rank * mat_rank);
cur_row_index = 2;
loopdepth = 1;
gen_comb_set(SetA, cand_mat, cand_row, mat_rank, loopdepth, cur_row_index);
MAXDET_MAT
fprintf(1, 'max detderminant = %d.\n', MAXDET);
end
%%----------------------------------------------------------------------
% Looking for the orthogonal rows and compute the determinant.
% \param SetA element candidate set
% \param cand_mat current candidate matrix
% \param cand_row current candidate row
% \param mat_rank rank of matrix (not exactly the rank, size of n)
% \param loopdepth parameter to simulate for-loop depth by recursion.
% \param cur_row current row index to look for
function gen_comb_set(SetA, cand_mat, cand_row, mat_rank, loopdepth, cur_row)
global MAXDET;
global MAXDET_MAT;
num_set = mat_rank;
num_cand = size(SetA);
szSetA = size(SetA);
% This should be assert(sum(szSetA == [1 2]) == 2)
if sum(szSetA == [1 2]) ~= 2
error('Not assumed set candidate matrix (should be 1x2)')
end
if cur_row > mat_rank;
% cand_mat;
det_a = det(cand_mat);
if det_a > MAXDET
MAXDET = det_a;
MAXDET_MAT = cand_mat;
end
elseif loopdepth > num_set
if check_orthogonal_row(cand_mat, cand_row, cur_row) == 1
cand_mat(cur_row, :) = cand_row;
cand_row = zeros(1, mat_rank);
cur_row = cur_row + 1;
gen_comb_set(SetA, cand_mat, cand_row, mat_rank, 1, cur_row);
end
else
% raw is not yet ready, generate it.
for j = 1:szSetA(2)
cand_row(loopdepth) = SetA(j);
gen_comb_set(SetA, cand_mat, cand_row, mat_rank, loopdepth + ...
1, cur_row);
end
end
end
%%----------------------------------------------------------------------
% check the rows are orthogonal with rows < cur_row
% \param cand_mat current candidate matrix
% \param cand_row current candidate row
% \param cur_row current row index to look for the orthogonal
function ret_is_all_orthogonal = check_orthogonal_row(cand_mat, cand_row, cur_row)
is_all_orthogonal = 1;
for i = 1:(cur_row - 1)
if dot(cand_mat(i,:), cand_row) ~= 0
is_all_orthogonal = 0;
break
end
end
ret_is_all_orthogonal = is_all_orthogonal;
end
しかし，これは 3, 5, あるいは問題にある 6 の場合に最大値を 0 と出力する．これはおかしい．なぜなら少なくとも determinant が 0 でないような行列を作成することは容易であり (ex. 3x3 の場合[1 1 1; 1 -1 -1 ; 1 1 -1])，これは0 より大きい determinant を持つ．しかし，ちょっと考えてみれば，3 の場合に1 と -1だけでは以下のように orthogonal な row ができないことがわかる．
geometry に戻ってみてもそれはわかる．1,-1 のみしか座標値としてとれない場合を考える．図 1に示すように直交するベクトルは 3次元の場合に既に作れない．
Figure 1: 3D, can not make orthogonal vectors by using {-1,1} coordinates
一応．2D の場合に直交法が上手くいく理由を，図2に示す．限られた座標であっても，直交するベクトルが作れることがわかる．図2にはdeterminant の性質としての面積も示しておいた．この大きさは (√2)^2 =2であり，2x2 の時の最大 determinant の 2 に等しい．
Figure 2: Orthogonal vectors by using {-1,1} coordinates in 2D case.
この方法は 1,2,4n (n >=1) の場合にしか使えない．ここで私はこのようなmatrix は Hadamard's matrix と呼ばれていることを知った．この問題は，Hadamard's Maximum Determinant Problem というようだ．Web で調べたところ，6の場合の答えはわかった．驚いたことにさらにかなりの数がわかっている．1,2,4n の場合には Hadamard matrix を構成する方法までわかっている．この数，1,2,4n を Hadamard 数と言うそうだ．4n の場合には，常に最大のdeterminant になるようだ．seehttp://mathworld.wolfram.com/HadamardsMaximumDeterminantProblem.html (これを見るともっと 4 の倍数でない場合には sharp な bound がわかるとあるので，4n の時最大でないのかと思ったが，実は，4n の時には Hadamard bound に到達できるので，これが最大であり，それ以外の場合にはもっと sharp (tight)な bound があるという意味である．)その上，matlab/octave には hadamard() という関数まであって， Hadamard 数の場合には matrix を作ってくれる．
しかし，Hadamard 数以外の大きさの最大 determinant をどうやって計算するのかがわからない．
http://mathworld.wolfram.com/HadamardsMaximumDeterminantProblem.html によると，1962 年に既にこの数列が知られていたらしい．おそらく賢い方法があるに違いない．
次回はもう少し geometry 的に考えた手法に関して紹介しよう．