2011-04-28

(6) Max determinant problem: Algorithm 3.5, another permutation

Algorithm 3.5: another permutation

I often go to lunch with my colleagues. At this point, I started to talk about this problem. It seems, Leo, Joachim, and Marc are interested in this story. I thought minimal dot product method is a good idea, so, I was kind of proud of that I found this simple method. My friends also agreed that this might work. But, the result is complete failure as shown in the last article.

Marc suggested me a geometrical approach. The max determinant of Hadamard matrix is
If you think about this is a geometrical problem, it is simple. The distance from origin to (1,1,1, ..., 1) coordinate in n-dimensional space is

(n-dimensional Pythagoras theorem). This is one of the longest distance edge. If these length edges are all perpendicular, then the volume of such object has \sqrt{n}^n. This is exactly the Hadamard's bound. The problem arises when these vectors can not be perpendicular. For example, this Strang's question. The problems in Strang's book are always like this. It seems simple problems, but, they are deep.

Joachim also suggested me that exploiting the matrix structure. But when he told me that, I thought minimal dot product is the answer. Therefore, I stupidly disregarded it. But my approach was completely failed, I needed to think about that.

First of all, if the matrix has the same row, the determinant is always zero. (Think about geometry, two axes are the same, then, no volume. Like two edges of a triangle are the same, the triangle is degenerated.) This we can omit it because any determinant value if exists, the max determinant is always more than zero. Even we have a minus value, just exchange two rows, then you got a plus determinant. (This means that if I have the max determinant, then I immediately have the min determinant by exchanging any two rows. Also think about geometry, flip two axes changes the direction of the volume, means the sign is changed.)

Based on this observation, we have the following algorithm: generate all the row candidate, then permute the rows. This narrows down 2^{36} candidates to _{2^6}P_{6} candidates in n = 6 case. But this is 53981544960 cases. This is less than 2^{36}, but the almost the same order.

This algorithm can only improve 27% of the time. It is not so much gain. Let's see the next idea.

2011-04-21

(6) Max determinant problem: Algorithm 3, min dot product

Algorithm 3: minimal dot product

I extended the geometry idea: if a pair of axes has minimal dot product, it could be a max determinant matrix. Minimal dot product means as perpendicular as possible. I implemented this idea in 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 generates a matrix that has relative large determinant value. But, they are not always the largest. For example, this program generates 32 when n = 5. But, the correct answer is 48 when n = 5. I could not proof why this doesn't work, if anybody know it, please let me know.

I needed to continue to seek for a better idea.

2011-04-20

(5) Max determinant problem: Algorithm 2, Orthogonality

Algorithm 2: using orthogonality

First I looked into the matrix pattern in 2x2 and 4x4. I saw the rows are orthogonal. I thought, ``Aha, because the determinant is volume and when a simplex has the maximal volume when the edge vector length is fixed? Orthogonal vectors!'' This is quite intuitive for me.

Therefore, I implemented a method that looked up the orthogonal vectors. This is 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

But, this program gives me the max determinant value is zero when 3x3, 5x5, and 6x6 matrix. This is strange. For instance, I can easily find a non zero determinant matrix, for instance, [1 1 1; 1 -1 -1 ; 1 1 -1] for 3x3. The determinant is 4. Also I realized I can not make a orthogonal rows in 3x3 case as the following.

When I think about the geometry, it is also easy to see it is not possible. Figure 1 shows we can not generates orthogonal vectors in 3D case when the coordinates value are only allowed {-1, 1}.
Figure 1: 3D, can not make orthogonal vectors by using {-1,1} coordinates
Figure 2 shows that this method works in 2D case. There are cases that we could make orthogonal vectors even the coordinate values are limited. Figure 2 also shows that the volume (= area, in 2D) that represents the determinant. This is (√2)^2 = 2, and the max determinant of 2x2 matrix is also 2.
Figure 2: Orthogonal vectors by using {-1,1} coordinates in 2D case.
This method can be applied to only 1,2,4n (n >=1) cases. At this point, I found these kind of matrices are called Hadamard's matrix. This problem is called Hadamard's Maximum Determinant Problem. On the Web, there is even the number (max determinant value) for 6x6 case. I am surprised that a lot of cases are known. In the case of 1,2,4n, there is a construction method to generate a Hadamard matrix. The number 1,2,4n is called Hadamard number.

Moreover, matlab/octave has function hadamard(), this generates a Hadamard matrix.

But, I didn't know how to compute the max determinant value of non-Hadamard number matrix.  According to http://mathworld.wolfram.com/HadamardsMaximumDeterminantProblem.html, the max determinant value sequence of Hadamard matrix is known in 1962. There should be a clever method.

2011-04-19

(4) Max determinant problem: Algorithm 1, Permutation method

Now, I introduce the problem and explain how I have developed my last answer. My first solution is correct, but, it is practically useless.  In chapter five of Introduction to Linear Algebra, Strang asked us a question, what is max determinant value of 6x6 {-1, 1} matrix? (Problem 33) This is a matlab question, so we can use matlab/octave. My first answer is generate all the combination of {-1,1}. This is Algorithm 1: permutation method.

Algorithm 1: permutation method

Since the component of matrix is limited to -1 or 1, we can generate all the permutation of 6x6 matrix and compute their determinant, then find the max determinant value. Program 1 shows the implementation.

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

I just ran this program on Intel Core2 2.26GHz machine's octave. But it doesn't stop after an hour. I switched to 5x5 matrix, it took two hours and j10 minutes. I should have first thought how large is the combination. In 5x5 matrix case, it is 2^{25}, which means, 33554432, 3.3 * 10^{7}.  In 6x6 case, it is 2^{36}, which means 68719476736, 6.8 * 10^{10}. If I run this program on my machine, it took more than 2200 hours, i.e., 92 days, almost three months. This is too much time for just a practice problem. I suspect there is a better way.

2011-04-18

(3) Max determinant problem, Appendix

This is just a side talk and you can skip this. It is a detail about relationship between Fredholm equation of the second kind and the max determinant problem. The basic idea is as following.

  • The discrete form of Fredholm equation of the second kind is a matrix form (a linear system).
  • To get the limit of the discrete form, the dimension of matrix n goes to infinity.
  • The solution of linear system is obtained by a variant of Cramer's rule. This needs the determinant.
  • The absolute value of the determinant relates with the system's convergence. Therefore, the max determinant problem was interesting.

Fredholm equation of the second kind is as following.
Assume this equation as a discrete problem, the range a,b is n-subdivided, then
Set λ (b-a)/n = h, the coefficient of this equation becomes following.
When we let the subdivision number n to infinity, the dimension of matrix becomes infinite. When we solve this linear system by Cramer's rule, we need determinant. Please note that Fredholm didn't use Cramer's rule directory. We need a bit more details (Fredholm minor, Fredholm determinant), however, I just want to know what was the motivation of the max determinant problem. (Also I don't understand whole the Fredholm equation discussion.)

(2) Max determinant problem

I would like to talk about why mathematicians are interested in max determinant problem. This is just my personal theory and I could not find an article that say this directly. So, I warn you that I might be completely wrong.

The max determinant problem is mentioned in a context of partial differential equation.  Is a partial differential equation interesting? I safely say, yes. This includes heat and wave problem. We can design buildings, computers, cars, ships, airplanes, ... and so on. There are so many applications of this in our world.

Hadamard is one of the mathematicians who contributed the max determinant problem. His one of the interests was partial differential equation. A basic partial differential equation, for instance, a wave equation is like this.
This can be re-written as
(By the way, if we wrote it as above, an operator d^2/d x^2 looks like to have an eigenvalue λ. Like Mu = -λu. This is a clue of relationship between integral equation and linear algebra.)

Fredholm wrote this kind of integral equation in a finite sum form. This is a matrix form. He had an idea to solve this equation by taking the limit of it, e.g., the dimension of matrix goes to infinite.

If we can write an integral equation in his form, each finite sum equation can be solved by solving linear equations. At that time, people solved the linear equations by Cramer's rule. Cramer's rule has 1/det(A) form. If a matrix A's size becomes larger, the solution can converge or not is depends on the (absolute) max determinant value, whether it is less than 1 or not. I imagine that this is the reason of mathematicians are interested in the max determinant problem. Though, I could not find a direct article about the motivation of why they are interested in the max determinant problem.

(Note: I explained only a little that how an integral equation has a matrix form. Because this is quite detail, I will add appendix in case one is interested in this.)

Although Fredholm didn't use Cramer's rule directly, the proof of convergence needs maximal value of determinant (see. 30 lectures of Eigenproblem, Shiga Kouji, p.121. (in Japanese)).

Hilbert removed the determinant from this problem and established eigenvalue based solution --- Hilbert space. He climbed up the view one level. I think determinant is still an important subject, though, eigenanalysis is much interesting. This is also just my impression, but, when Hilbert established the Hilbert space, people's interest moved to eigenvalues from the determinant. I just imagine this.

2011-04-16

(1) Max determinant problem

Abstract

Gilbert Strang asked us what is the maximal determinant if the matrix has only specific numbers in his book, Introduction to Linear algebra. I enjoyed this problem for almost three weeks also as a programming problem. So I would like to introduce this problem in this article.


Introduction

My primary school has words, ``Be one day as one step of your life (一日生きることが一歩生きることであれ.)'' by Yukawa Hideki. These days I finally start to understand these words. I can only do something if I could do every day. Even for five minutes, if I do something every day, I found quite difference. Recently, I joined an activity. It took some significant time from my Sunday research time, though I would like to continue both my activity and my Sunday research.

At the end of March, I learn max determinant problem that exists. I didn't have any dedicated time for this problem. But, I use my commune time and elevator waiting time, I solved this problem. (Our company's elevator gives a lot of time, I usually use it for reading a book.) Using fraction time might be a point of continue something.

I didn't know why max determinant problem caught an interest of mathematicians until I started to solve this problem. In a Gilbert Strang's class, he said ``Determinant used to be very important for linear algebra''.  It was a past tense. I didn't recall that he mentioned why it was once important and not now anymore. If we think about a matrix as a linear operator, determinant is zero or not is important since it tells the system has a solution or not. I thought maximal value is not so important comparing to this.

The determinant of a matrix is magnification factor when we think the matrix is an operator. Why this is interesting? I could imagine that the absolute maximal value is less than one or not is interesting. We usually think about multiplication of matrix. For example, M^k v. But, in this case, eigenvalue is much interesting. Since if we could know the eigenvalues, this becomes M^ k v = λ^k v. This is much simpler and easy because a matrix becomes now one scalar value.

I can also think about another property of determinant, geometrical meaning. This is a volume of limited coordinates geometry. (Marc also noted this to me.) I like geometry, so, in the following articles, I will use this approach once. But, is it really interesting? This was a question to me.

I research why mathematicians are interested in max determinant problem a bit. I could not find the direct answer, but, I have an idea about that. So, I will tell about that in the next article.