(8) 最大 determinant 問題: Algo 4. combination

Algorithm 4: combination

Determinant は row の入れかえでは符号しか変化しない．したがって，順列は必要ではなく，組合せのみが必要である．nxn の matrix の場合，row は n の要素を持つ．これの {-1,1}の順列は 2^6 = 64 である．したがって，この組合せの数は，_{2^6}C_{6} = 74974368 である．これは 2^{25} の倍程度であるから，全ての場合を計算しても5時間程度で済むのではないかと考えた．その実装がprogram 4 である．

Program 4
function MaxDeterminant = algo_04(matrix_rank)
% Introduction to linear algebra Chapter 5. problem 33
% Algorithm 4: combination method
% @author Hitoshi
if nargin ~= 1;
error('Usage: la_chapt5_33_comb_row(matrix_rank).')
end

MatrixRank = matrix_rank;
% generate all the row combination (simple permutation)
CombMat = gen_combinatorial_matrix(matrix_rank);
comb_mat_size = size(CombMat);
CombRowCount = comb_mat_size(1);
curChoise = 1:MatrixRank;

global MaxDet MaxDetMat
MaxDet = 0;
MaxDetMat = [];

tic
while (1)
mat = CombMat(curChoise,:);
d = det(mat);
if d > MaxDet
MaxDet = d;
MaxDetMat = mat;
end

find_idx = 0;
for i = MatrixRank:-1:1
if curChoise(i) < CombRowCount - (MatrixRank - i)
find_idx = i;
break
end
end

if find_idx == 0
break % done
else
start_val = curChoise(find_idx) + 1;
curChoise(:,find_idx:MatrixRank) = start_val:(start_val + MatrixRank - find_idx);
end
end
toc

MaxDet
MaxDeterminant = MaxDetMat;
end


実は combination を使うアイデアは既にあったので，algorithmn 3.5 のpermutation の手法を実装するつもりはなかったのだが，combination のつもりで実装したものが，実は permultation であった．とんでもない間違いをしてしまった．permutation と combination の違いはここでは 6! 通り，つまり，720倍違う．計算すると permutation 法では 60 日ほどかかることがわかったが，720倍高速な方法というのは 2 時間しかかからない．

答えはあっている．実は web では得らなかった副作用として，matrix も得られた．以下が 6x6 の {-1,1} matrixの最大の determinant を持つ matrix の一つである．

gen_combinatorial_matrix()は row の permutation を作るものだが，これは本質ではないし，実装も簡単なので省略する．

(7) 最大 determinant 問題: Algo 3.5. another perm.

Algorithm 3.5: permutation

私の会社では同僚と Lunch に行く．そこではこういう話をすることがある．Leoと Joachim がこの話に興味を持った．minimal dot product 法は私は良いアイデアだと思って得意になって同僚に話をし，同僚も悪いアイデアではないと思ったようだ．私はこれは高速だと思っておおはしゃぎしたが，結果は前回の通り，失敗だった．

Marc は geometry 的なアプローチを示唆した．これは Hadamard matrix を考える時には良い clue になる．Hadamard の上界は次式のようになっている．

これは geometry で考えると簡単である．まず，n 次元の原点から (1,1,1,..., 1) 座標までの距離は，n 次元のピタゴラスの定理から

である．これが全て直交した図形の体積は，そのn 乗であるから，上式が上界を示していることは簡単だ．問題は直交軸を形成できない次元の場合が問題である．例えば Strang の質問がそうである．Strang の本の問題は簡単そうに見えて，常にこういうことを考えて作った問題なのだ．

Joachim は matrix の構造をもっと考えれば組合せを絞れるという話をしていたのだが，minimal dot product 法はそういうことをあまりせずとも候補を絞れるので単純で賢いと思って考えなかったのだが，失敗したために考えることにした．

改良手法として，まったく同じ row を持つものは意味がない．その場合にはdeterminant は 0 になるからどんな determinant でも 0 以外ならば最大値になる候補になる．(0になる理由は，Geometry から考えると簡単である．ある体積を張るようなベクトルの二つが同一になる場合，体積は0である．たとえば，三角形の2つの辺が重なってしまう場合にはその三角形の面積は0になるのと同じ理由で
ある．) row 交換は determinant の符号を変化させることを考慮すると，負の determinant でも 0 でなければただちに正の determinant を求められる．これは最大 determinant が求まれば，2つの row を交換することで，直ちに最小も求まることを意味している．(また Geometry を考えると簡単である．ある2つの辺を入れ替えると，その体積は符号が変化するだけで体積の絶対値は変化しない．)


したがって，row の候補を組合せで生成し，その permutation を考えることにする．これは 2^{36} の候補ではなく，_{2^6}P_{6} を考える方法である．しかし，これはまだ，53981544960 通りある．確かに 2^{36} よりは小さいのだが，ほとんど同じ orderである．計算してみると 27% 程度の改善しかみこめない．

この方法は速度の改善をあまりしないので次の方法に行くことにしよう．

(6) 最大 determinant 問題: Algo 03, min dot product

Algorithm 3: minimal dot product

私は geometry の手法を押しすすめて各軸が最小の dot product を持つものがいいのではと考えた．それが program 3 である．

Program 3

function algo_03(matrix_rank)
% Introduction to linear algebra Chapter 5. problem 33
% Algorithm 3: find minimal dot product vector
% @author Hitoshi
if nargin ~= 1;
error('Usage: algor_03(matrix_rank).')
end
global MatrixRank;
MatrixRank = matrix_rank;
global CombMat;
CombMat = gen_combinatorial_matrix(matrix_rank);
% initialize candidate matrix
cand_mat = zeros(matrix_rank, matrix_rank);
cand_mat(1,:) = CombMat(1,:);
comb_mat_size = size(CombMat);
global CombRowCount;
CombRowCount = comb_mat_size(1);

for i = 2:matrix_rank
min_dotprod_row_idx = get_min_dotprod_row(cand_mat, i);
cand_mat(i,:) = CombMat(min_dotprod_row_idx, :);
end

cand_mat
det(cand_mat)
end
%%----------------------------------------------------------------------
% get minimal dotproduct row
% \param[in] cand_mat candidate matrix
% \param[in] cur_row_idx current last row index
function min_dp_row_idx = get_min_dotprod_row(cand_mat, cur_row_idx)
global MatrixRank;
global CombMat;
global CombRowCount;
% init dot prod with the large one
min_dotprod = dot(cand_mat(1,:), cand_mat(1,:)) * MatrixRank;
min_dotprod_rowidx = 1;

for i = 2:CombRowCount % 1 has been used
check_row = CombMat(i,:);
% skip the same row if exists
if check_same_row(cand_mat, check_row, cur_row_idx) == 1
continue;
end
% check min dot product row
sum_dotprod = 0;
for j = 1:cur_row_idx
sum_dotprod = sum_dotprod + abs(dot(cand_mat(j,:), check_row));
end
if min_dotprod > sum_dotprod
min_dotprod = sum_dotprod;
min_dotprod_rowidx = i;
end
end
min_dp_row_idx = min_dotprod_rowidx;
end
%%----------------------------------------------------------------------
% check the same row exists
% \param[in] cand_mat candidate matrix
% \param[in] check_row checking row entry
% \param[in] cur_row_idx current row index
% \return 1 when found the row
function ret = check_same_row(cand_mat, check_row, cur_row_idx)
is_found = 0;
% check the same entry in the candidate?
for j = 1:cur_row_idx
if all(check_row == cand_mat(j,:)) == 1
is_found = 1;
break; % the same entry found
end
end
ret = is_found;
end


program 3 はある程度の大きさの determinant を持つ matrix を生成するが，知られている数列の答えと合わない．5 の場合には 32 を答えとして返すが，この場合は48が最大値であり，最小の dot product を選択するという方法では真の最大の determinant を求められない場合がある．(これが何故か証明できる方は御一報下さい．)

そこでもっと賢い方法はないか考えることにした．