Calculate the biggest sum of submatrix

 Let's continue the topic discussed in last article. What if we want to calculate the biggest sum of the submatrices of a matrix, which can be treated as a 2-D array? The usual idea is that we can solve this problem by "reducing" the "size" of the problem, that is, reduce 2-D to 1-D, and then, we have already the solution of linear complexity at hand. So how to reduce the problem? Here is one way: 1. Choose an upper and a lower bound of lines and treat every fragment of the columns within these bounds as a whole, which means only the sum of that fragment is taken into consideration; 2. Enumerate the upper and lower bounds and repeat the 1st step. Finally we can get the biggest sum of submatrix (as well as the starting and ending elements' line indices and column indices). Let's look closer. The 2nd step is a m*n loop, where m is the number of lines and n is the number of columns of the original matrix. To calculate biggest sum of 1-D array is of complexity O(n), where n is the number of array elements. Thus the total complexity of this process is O(m*n*n')*(the complexity of calculate the sum of column fragments) (here n' = n). In the code shown below, a small trick trading space cost off time cost is used to make the complexity of summing the column fragments be O(1) and thus the total complexity is O(m*n*n).

 

1. #include <iostream>
2. using namespace std;
4. // Calculate the sums of the column fragments.
5. int* partial_sums_of_matrix(int* A, int m, int n) {
6.     // Allocate matrix P.
7.     int new_m = m + 1, new_n = n + 1;
8.     int* P = new int[new_m * new_n];
9.     // Initialize matrix P.
10.     for(int i = 0; i < new_m; ++i) { // Zero the left-most column.
11.         *(P+i*new_n) = 0; // *(P+i*new_n+0) = 0;
12.     }
13.     for(int i = 0; i < new_n; ++i) { // Zero the up-most line.
14.         *(P+i) = 0; // *(P+0*new_n+i) = 0;
15.     }
16.     // Calculate the sums of column fragments of matrix A and Store them in matrix P.
17.     for(int i = 1; i < new_m; ++i) {
18.         for(int j = 1; j < new_n; ++j) {
19.             *(P+i*new_n+j) = *(P+(i-1)*new_n+j) + *(P+i*new_n+(j-1)) - *(P+(i-1)*new_n+(j-1)) + *(A+(i-1)*n+(j-1)); // !Notic here: NOT *(A+i*n+j)!
20.         }
21.     }
22.     return P;
23. }
25. // Get the partial sum of a particular submatrix region.
26. inline int get_a_partial_sum(int* P, int tm, int bm, int on, int cn) {
27.     return *(P+bm*on+cn) - *(P+bm*on+(cn-1)) - *(P+(tm-1)*on+cn) + *(P+(tm-1)*on+(cn-1));
28. }
30. // Calculate max sum of subarray of 1-D arrays (just the same as shown in last article).
31. int max_sum_of_subarray(int* A, int n, int& start_index, int& end_index) {
32.     int max_all, max_pre; // Store the max sum of the pre-subarray.
33.     int temp_start_index, temp_end_index; // Bookkeeping during the process.
34.     max_pre = max_all = A[0];
35.     temp_start_index = temp_end_index = start_index = end_index = 0;
36.     for(int i = 1; i < n; ++i) {
37.         if(A[i] > max_pre + A[i]) { // A new "bigger" subarray appears.
38.             max_pre = A[i];
39.           temp_start_index = temp_end_index = i;
40.         }
41.         else { // Just add the ith element into the old subarray.
42.             max_pre += A[i];
43.             temp_end_index = i;
44.         }
45.         if(max_all < max_pre) { // Adjust the final max sum and the indices.
46.             max_all = max_pre;
47.             start_index = temp_start_index;
48.             end_index = temp_end_index;
49.         }
50.     }
51.     return max_all;
52. }
54. // Calculate max sum of submatrix of 2-D matrices.
55. int max_sum_of_submatrix(int* P, int m, int n, int& start_m, int& start_n, int& end_m, int& end_n) {
56.     int max_pre, max_all = 0x80000000;
57.     int temp_start_m, temp_start_n, temp_end_m, temp_end_n;
58.     int* temp_p = new int[n];
59.     for(int t = 1; t < m; ++t) {
60.         temp_start_m = t - 1; // temp_start_m is the starting line index of the original matrix.
61.         for(int b = t; b < m; ++b) {
62.             temp_end_m = b - 1; // temp_end_m is the starting line index of the original matrix.
63.             for(int j = 1; j < n; ++j) { // Put the partial sums in a 1-D array.
64.                 temp_p[j] = get_a_partial_sum(P, t, b, n, j);
65.             }
66.             // Calculate 1-D.
67.             // temp_start_n and temp_end_n are the ending element's indices of the original matrix.
68.             max_pre = max_sum_of_subarray(temp_p+1, n-1, temp_start_n, temp_end_n);
69.             if(max_pre > max_all) {
70.                 max_all = max_pre;
71.                 start_m = temp_start_m;
72.                 start_n = temp_start_n;
73.                 end_m = temp_end_m;
74.                 end_n = temp_end_n;
75.             }
76.         }
77.     }
78.     delete[] temp_p;
79.     return max_all;
80. }
82. void print_matrix(int* A, int m, int n) {
83.     for(int i = 0; i < m; ++i) {
84.         for(int j = 0; j < n; ++j) {
85.             cout << *(A + i * n + j) << " ";
86.         }
87.         cout << endl;
88.     }
89. }
91. int main() {
92.     int A[3][5] = {
93.         {-2, 4, 1, 6, -7},
94.         {3, 5, 1, -8, 0},
95.         {10, -9, 2, 0, -3}};
96.     int* P = partial_sums_of_matrix(&A[0][0], 3, 5);
97.     int max_sum, start_m, start_n, end_m, end_n;
98.     max_sum = max_sum_of_submatrix(P, 3+1, 5+1, start_m, start_n, end_m, end_n);
99.     cout << "max sum: " << max_sum << "/trange: [" << start_m << ", " << start_n << "], [" << end_m << ", " << end_n << "]" << endl;
100.     return 0;
101. }

The output is like this:

max sum: 15        range: [0, 0], [2, 2]

 

Similarly, we can also enumerate the left and right bounds with the sums of line fragments to get the final biggest sum of submatrix, and the complexity becomes O(m*n*m). If m and n are of the same scale, there is little difference. Otherwise we need to find the smaller of the two to get the most efficient process. In a word, the complexity is O(m*n*min(m, n))).

 

Referenced to <Beauty of Programming>.