Tadao Takaoka algorithm for finding the maximum submatrix or Maximum Subarray Problem
Not so long ago passed the competition of parallel programming Acceler8 2011 . The essence of the task was to find the maximum submatrix in this matrix (the sum of the elements of the found submatrix must be maximum). After a short googling, it was found that a certain algorithm of Tadao Takaoka solves this problem faster than others.
“The call is accepted!”, And I began to look for this algorithm wherever possible, setting a goal to implement it. Despite the fact that it is parallelized badly and in its complexity it contains a rather big constant.
However, all that was found was articles in English by this very Tadao Takaoki (here is one of these articles ). I had to translate.
')
The very idea of the algorithm at first seemed to be ugly simple:
Here is what is written about finding this solution in the original article:
" Now the column-centered problem can be solved as follows.
A_column = max (k = 1: i-1, l = 0: n / 2-1, i = 1: m, j = n / 2 + 1: n) {s [i, j] - s [i, l] - s [k, j] + s [k, l]}.
We fix i and k, and maximize the remaining expression by changing l and j. Then it is equivalent to the following expression:
i = 1, ..., m and k = 1, ..., i - 1.
A_column [i, k] = max (l = 0: n / 2-1, j = n / 2 + 1: n) {−s [i, l] + s [k, l] + s [i, j ] - s [k, j]}
Let s ∗ [i, j] = −s [j, i]. Then this expression can be reduced to:
A_column [i, k] = −min (l = 0: n / 2-1) {s [i, l] + s ∗ [l, k]} + max (j = n / 2 + 1: n) { s [i, j] + s ∗ [j, k]} "
The following is a reference to some other problem of Tadao Takaoka: Distance Matrix Multiplication (I don’t undertake to translate this term).
The bottom line is:
The goal of DMM is to calculate the product of distances (distance product) C = AB for two n-dimensional matrices A = a [i, j] and B = b [i, j], whose elements are real numbers.
c [i, j] = min (k = 1: n) {a [i, k] + b [k, j]}
The essence of c [i, j] is the smallest distance from vertex i from the first layer to vertex j in the third layer in an acyclic oriented graph consisting of three layers of vertices. These vertices are numbered 1, ..., n in each layer, and the distance from i in the first layer to j in the second layer is a [i, j] and from i in the second layer to j in the third layer is b [i, j ]. Obviously, this expression can be easily given to calculate the maximum instead of the minimum - this will be the problem of finding the largest paths.
So using this definition, you can get the following:
In the formula A_column [i, k] = −min (l = 0: n / 2-1) {s [i, l] + s ∗ [l, k]} + max (j = n / 2 + 1: n ) {s [i, j] + s ∗ [j, k]} it is clear that the first part of the formula is DMM for a minimum, the second is DMM for a maximum. Let S1 and S2 be matrices whose elements (i, j) are s [i, j - 1] s [i, j + n / 2], respectively. For any matrix T, let T ∗ be the matrix obtained from T by transposition and negation. Then the upper expression can be calculated by “multiplying” S1 and S1 ∗ for a minimum and obtaining a lower triangular matrix, “multiplying” S2 and S2 ∗ for a maximum and obtaining a lower triangular matrix, and finally subtracting from the last previous matrix. Calculating the maximum in it, we get the answer:
A_column = max (DMM_max (S2 * S2 ∗) - DMM_min (S1 * S1 ∗)).
Next, pseudo-code is translated to solve the DMM problem for matrices A and B of size nxn:
The article recommends using a binary heap for queues. In later versions it was stated that the use of the fibonacci heap could speed up the process even more, however, to my shame, I could not fully implement it. Therefore, I used binary.
However, even by implementing this pseudocode, I found that the solution was not working. The obvious error was found at the stage of calculating S1 and S2 - the described rule for constructing these matrices was not enough for the algorithm to work properly. And I have not found a more detailed description in any article. Thus it was necessary to build the matrix itself.
So, after long hours spent searching and checking, the following rules were obtained.
Let S be the prefix matrix for matrix A.
Then the matrix S1 of size m on n / 2 contains the first column of 0, and the remaining n / 2-1 columns contain the first columns of the prefix matrix S:
The matrix S1 ∗ with dimensions n / 2 on m contains the first row and a column of 0, the remaining submatrix is the negation and transposition of the first rows and columns of the prefix matrix S:
(k = n div 2 + n mod 2)
The matrix S2 of size m on k contains the columns of the prefix matrix S, starting with the k-th:
Matrix S2 ∗ with dimensions of k on m: the first column of 0, then the transposition and negation of the columns of the prefix matrix S, starting with the k-th column:
In general, the algorithm turned out to be very difficult to understand and implement and did not meet its expectations for exceptional speed (no additional optimizations were carried out).
The same article claims that this algorithm works in O (n ^ 3 * log (log (n)) / log (n)) time, and also mentions that the algorithm runs faster than other algorithms on large matrices ( apparently very large).
If someone is interested or vital, here are the sources of my implementation (it is written on the pluses, so I do not guarantee “intelligibility”). I would be very happy if this article will help someone (no wonder I spent so much time on the proceedings with this algorithm!).
“The call is accepted!”, And I began to look for this algorithm wherever possible, setting a goal to implement it. Despite the fact that it is parallelized badly and in its complexity it contains a rather big constant.
However, all that was found was articles in English by this very Tadao Takaoki (here is one of these articles ). I had to translate.
')
The very idea of the algorithm at first seemed to be ugly simple:
First, we need to calculate all prefix sums s [i, j] for such matrices as a [1..i, 1..j], for all i and j with the constraint s [i, 0] = s [0, j] = 0. Obviously, this is done in O (mn) time.If A_left and A_right are through recursion, then you must manually find A_column.
The main algorithm is:
If the matrix is one element, we return its value.
Otherwise:Thus, m ≤ n.
- if m> n, rotate the matrix by 90 degrees.
- Let A_left be the solution for the left half.
- A_right is the solution for the right half.
- A_column - solution for column-centered tasks (finding such a maximum submatrix that would cross the center line)
- The overall solution is a maximum of these three.
Here is what is written about finding this solution in the original article:
" Now the column-centered problem can be solved as follows.
A_column = max (k = 1: i-1, l = 0: n / 2-1, i = 1: m, j = n / 2 + 1: n) {s [i, j] - s [i, l] - s [k, j] + s [k, l]}.
We fix i and k, and maximize the remaining expression by changing l and j. Then it is equivalent to the following expression:
i = 1, ..., m and k = 1, ..., i - 1.
A_column [i, k] = max (l = 0: n / 2-1, j = n / 2 + 1: n) {−s [i, l] + s [k, l] + s [i, j ] - s [k, j]}
Let s ∗ [i, j] = −s [j, i]. Then this expression can be reduced to:
A_column [i, k] = −min (l = 0: n / 2-1) {s [i, l] + s ∗ [l, k]} + max (j = n / 2 + 1: n) { s [i, j] + s ∗ [j, k]} "
The following is a reference to some other problem of Tadao Takaoka: Distance Matrix Multiplication (I don’t undertake to translate this term).
The bottom line is:
The goal of DMM is to calculate the product of distances (distance product) C = AB for two n-dimensional matrices A = a [i, j] and B = b [i, j], whose elements are real numbers.
c [i, j] = min (k = 1: n) {a [i, k] + b [k, j]}
The essence of c [i, j] is the smallest distance from vertex i from the first layer to vertex j in the third layer in an acyclic oriented graph consisting of three layers of vertices. These vertices are numbered 1, ..., n in each layer, and the distance from i in the first layer to j in the second layer is a [i, j] and from i in the second layer to j in the third layer is b [i, j ]. Obviously, this expression can be easily given to calculate the maximum instead of the minimum - this will be the problem of finding the largest paths.
So using this definition, you can get the following:
In the formula A_column [i, k] = −min (l = 0: n / 2-1) {s [i, l] + s ∗ [l, k]} + max (j = n / 2 + 1: n ) {s [i, j] + s ∗ [j, k]} it is clear that the first part of the formula is DMM for a minimum, the second is DMM for a maximum. Let S1 and S2 be matrices whose elements (i, j) are s [i, j - 1] s [i, j + n / 2], respectively. For any matrix T, let T ∗ be the matrix obtained from T by transposition and negation. Then the upper expression can be calculated by “multiplying” S1 and S1 ∗ for a minimum and obtaining a lower triangular matrix, “multiplying” S2 and S2 ∗ for a maximum and obtaining a lower triangular matrix, and finally subtracting from the last previous matrix. Calculating the maximum in it, we get the answer:
A_column = max (DMM_max (S2 * S2 ∗) - DMM_min (S1 * S1 ∗)).
Next, pseudo-code is translated to solve the DMM problem for matrices A and B of size nxn:
Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.
Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.Copy Source | Copy HTML / n B list[ i ], .. list[ i ] i ; / V = { 1 , ..., n}; for i := 1 to n do begin for k := 1 to n do begin cand[k]:=first of list[k]; d[k] := a [ i , k] + b [k, cand[k]]; end; / V d[ 1 ], ..., d[n]; / S - ; /* 1 : */ while |S| ≤ n − n log n do begin / v ; / cand[v] S; c[ i , cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; for w in W do while cand[w] is in S do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; U := S; /* 2 : */ while |S| < n do begin / v ; if cand[v] is not in S then begin / cand[v] S; c[v, cand[v]] := d[v]; / W = {w|cand[w] = cand[v]}; end else W = {v}; for w in W do cand[w]:=next of list[w]; while cand[w] is in U do cand[w]:= next of list[w]; / W d[w] = a [ i , w]+ b [w, cand[w]]; end; end.
The article recommends using a binary heap for queues. In later versions it was stated that the use of the fibonacci heap could speed up the process even more, however, to my shame, I could not fully implement it. Therefore, I used binary.
However, even by implementing this pseudocode, I found that the solution was not working. The obvious error was found at the stage of calculating S1 and S2 - the described rule for constructing these matrices was not enough for the algorithm to work properly. And I have not found a more detailed description in any article. Thus it was necessary to build the matrix itself.
So, after long hours spent searching and checking, the following rules were obtained.
Let S be the prefix matrix for matrix A.
Then the matrix S1 of size m on n / 2 contains the first column of 0, and the remaining n / 2-1 columns contain the first columns of the prefix matrix S:
The matrix S1 ∗ with dimensions n / 2 on m contains the first row and a column of 0, the remaining submatrix is the negation and transposition of the first rows and columns of the prefix matrix S:
(k = n div 2 + n mod 2)
The matrix S2 of size m on k contains the columns of the prefix matrix S, starting with the k-th:
Matrix S2 ∗ with dimensions of k on m: the first column of 0, then the transposition and negation of the columns of the prefix matrix S, starting with the k-th column:
In general, the algorithm turned out to be very difficult to understand and implement and did not meet its expectations for exceptional speed (no additional optimizations were carried out).
The same article claims that this algorithm works in O (n ^ 3 * log (log (n)) / log (n)) time, and also mentions that the algorithm runs faster than other algorithms on large matrices ( apparently very large).
If someone is interested or vital, here are the sources of my implementation (it is written on the pluses, so I do not guarantee “intelligibility”). I would be very happy if this article will help someone (no wonder I spent so much time on the proceedings with this algorithm!).
Source: https://habr.com/ru/post/136188/
All Articles