how to make a sum of submatrices

Question

i need to make for every sub-matrix the sum of the values. for example if i have [[1,1,2],[2,3,4]] the resulting matrix will be:

 M[0][0] = 1 M[0][1] = 1+1 = 2 M[0][2] = 1+1+2 = 4
 M[1][0] = 1+2 = 3 M[1][1] = 1+1+2+3 = 7 M[1][2] = 1+1+2+2+3+4 = 13

or

M = [[1,2,4],[3,7,13]]

and i made this code

`N = []
 M = []
 summ = 0
 n= list(map(int, input().split()))
 while n != []:
     N.append(n)
     n = list(map(int, input().split()))
 for i in range(len(N)):
     M.append([0 for i in range(len(N[0]))])
     summ = 0
     for j in range(len(N[0])):
         summ += N[i][j]
         M[i][j] = M[i-1][j] + summ ` 

the problem is that when the matrix is big becomes reeeeally slow. i need to solve a 100x100 matrix at max in 0.5 sec

can anybody help me? WITHOUT IMPORTING PACKAGES!! `

Answer 1

For speed, you really want to be using NumPy which will be significantly faster than base Python in addition to giving cleaner code for matrices. From your small example, you can use numpy.cumsum() twice across different axes:

import numpy as np

arr = np.array([[1,1,2],[2,3,4]])
out = arr.cumsum(axis=1).cumsum(axis=0)

print(out)

Gives:

array([[ 1,  2,  4],
       [ 3,  7, 13]], dtype=int32)

Side note: on Windows, the default int type is 32 Bit and cumsum() is liable to silently overflow on large matrices/large numbers so you'll probably want arr = np.array([[1,1,2],[2,3,4]]).astype(np.int64) if on Windows.

Timings:

arr = np.arange(10000).reshape(100, 100)
%timeit out = arr.cumsum(axis=1).cumsum(axis=0)
56.3 µs ± 4.96 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

So, thousands of times faster than your requirement.

Answer 2

We can compute the sum of all sub-matrix by using a dp matrix. In this we compute the prefix sum of the given matrix.We create a new matrix and store the prefix sum and then use the formula

tlx: Top left x coordinate
tly: Top Left y coordinate
brx: Bottom right x coordinate
bry: Bottom Right y coordinate

Sum between (tlx, tly) and (brx, bry) :

sum += dp[brx][bry] - dp[tlx - 1][bry] - dp[brx][tly - 1] + dp[tlx - 1][tly - 1];

Let me explain: We first have a given matrix suppose Arr :

#define n 3
int row(n), col(n), sum(0);
int arr[row][col] = {{1, 2, 3},
                     {4, 5, 6},
                     {7, 8, 9}};

Then we will create a prefix matrix with one size greater and we will fill all values as 0 initially.

int dp[row + 1][col + 1];
memset(dp, 0, sizeof(dp));

Then will create our prefix matrix by:

1. Copy first row from our arr to dp

   for (int j = 1; j < col + 1; j++){dp[1][j] = arr[0][j - 1];}
   

2. Run a loop and add values of first row of dp to second col of arr and store it in second col of dp

   for (int i = 2; i < row + 1; i++)
    {
        for (int j = 1; j < col + 1; j++)
        {
            dp[i][j] += dp[i - 1][j] + arr[i - 1][j - 1];
        }
    }
   

3. Repeat this process till you fill the complete dp.

4. Then run a loop through col in which we will add each col to their next col and save it just like we create prefix array.

   for (int i = 0; i < row + 1; i++) {
   for (int j = 1; j < row + 1; j++)
        {
            dp[i][j] += dp[i][j - 1];
        }
    }
   

5. Now your dp (Prefix Matrix is ready).

    arr : {{1, 2, 3},
          {4, 5, 6},
          {7, 8, 9}};
    dp : 0 0 0 0 
         0 1 3 6 
         0 5 12 21 
         0 12 27 45 
   

We took one extra size for our dp array and put value 0 in the first row and col because with this it will be easier for us to use the formula which I mentioned at above. We divide matrices into rectangles and with the help of top left position and bottom right position we calculate the occurrences in the given matrix. Now we just have to iterate through each position from top left to bottom right and continuously add the sum and then print it.

for (int tlx = 1; tlx < row + 1; tlx++){
  for (int tly = 1; tly < col + 1; tly++){
    for (int brx = tlx; brx < row + 1; brx++){
     for (int bry = tly; bry < col + 1; bry++){
       sum += (dp[brx][bry]) - (dp[tlx - 1][bry]) - (dp[brx][tly - 1]) +   (dp[tlx - 1][tly - 1]);
    }
  }
}
} 

OUTPUT: 500