(5) 最大 determinant 問題: Algo 02, orthogonality

Algorithm 2: using orthogonality


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 的に考えた手法に関して紹介しよう．

(4) 最大 determinant 問題: Algo 01, Permutation

では問題を示し，どうやって解いたかについて具体的に順を追って説明していく．まずは正しいが，実用的でない方法から説明していく．

Introduction to Linear Algebra の第5章の問題33で，Gilbert strang は{-1,1}のみを使った 6x6 の matrix の deterinant の最大値を matlab を使って計算してみることを求めている．私が最初に考えたのは単純に全ての組合せを使うものであった．それが Algorithm 1: permutation 法である．

Algorithm 1: permutation method

Matrix の要素は -1 と 1 に限られているのであるから，全ての順列を生成してその determinant を計算し，最大値を求めることができる．それが Program 1である．

ところで Amoeba ではプログラムを整形する方法がわからないので，
http://shitohichiumaya.blogspot.com/2011/04/4-max-determinant-problem-algorithm-1.html
を参照して欲しい．

Program 1
function algo_01(mat_rank)
% Introduction to linear algebra Chapter 5. problem 33
% Algorithm 1: generate all the combination and test method.
% @author Hitoshi
if nargin ~= 1;
error('Usage: algo_01(mat_rank)')
end

% possible element set A_i = {-1, 1}
SetA = [-1 1];
cand_mat = zeros(1, mat_rank * mat_rank);

global MAXDET
global MAXDET_MAT

% We know det I = 1, therefore, it is a good initial value.
MAXDET = 1;
MAXDET_MAT = zeros(1, mat_rank * mat_rank);

gen_comb_set(SetA, cand_mat, mat_rank, 1);

MAXDET_MAT
fprintf(1, 'max detderminant = %d.\n', MAXDET);
end

%%----------------------------------------------------------------------
% Execute depth num_set loop by recursion to generate matrix candidates.
%
% \param SetA element candidate set
% \param cand_mat current candidates
% \param mat_rank rank of matrix (not exactly the rank, size of n)
% \param loopdepth parameter to simulate for-loop depth by recursion.
function gen_comb_set(SetA, cand_mat, mat_rank, loopdepth)

global MAXDET;
global MAXDET_MAT;

num_set = mat_rank * 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 loopdepth > num_set
MA = reshape(cand_mat, mat_rank, mat_rank);
det_a = det(MA);
if det_a > MAXDET
MAXDET = det_a;
MAXDET_MAT = MA;
end
else
for j = 1:szSetA(2)
cand_mat(loopdepth) = SetA(j);
gen_comb_set(SetA, cand_mat, mat_rank, loopdepth + 1);
end
end
end

このプログラムを Intel Core2 2.26 GHz, octave で走らせてみたら1時間たっても終わらない．そこで 5x5 の場合を計算してみたら，2時間10分ほどかかった．まず最初に考えるべきだったが，いくつの場合を試すのか考えてみると，5x5のmatrix では 2^{25} 通りあるので， 33554432, 3.3 * 10^{7} である．6x6の場合は 2^{36} 通りあるので， 68719476736, 6.8 * 10^{10} である．順当に行けば 2200 時間ほどあれば解ける，つまり 92 日，3ヶ月あればよいのだが，それが練習問題になっているというのはどうもおかしい．もっと賢い方法があるようである．

(3) 最大 determinant 問題 補足

ここは余談であるので，飛ばしてもらってもかまわない． Fredholm equation of the second kind がどう行列式とかかわってくるかの基本的なアイデアは以下である．

- Fredholm の第二種積分方程式を離散値で考えると，行列の形になること．
- その極限をとることは，行列の次元 n を無限に持っていくこと．
- Cramer の解法に関連した手法で方程式を解くため，determinant が必要であったこと．
- この場合，determinant の絶対値の最大値が収束性に関わるので，最大 determinant が問題になった．

Fredholm Equation of the second kind は以下の形である．

離散値で考える．区間a,bをn等分すると以下のようになる．

λ (b-a)/n = h とすると，この式の係数は以下になる．

分割数 n を無限大に持っていくと，行列が無限次元になることがおわかりだろう．これを Cramer の解法で解こうとすると，determinant がでてくる．ただし，Fredholm は Cramer の解法を使ったわけではなく，もう少し詳細が必要になる(Fredholm の小行列式など)が，ここでは determinant の最大値を求めたいという動機が感じられれば良いと思うのでこれ以上は立ち入らない．(正直なところ，私が理解できていないので．)