Processing a 2D matrix - need to speed up my O(n^4) algorithm

I have an `n x n` matrix which I need to convert into a list sorted by value. Starting with the maximum value cell at (row `x1`, col `y1`), I must immediately exclude all cells where `(x >= x1, y <= y1)` or `(x <= x1, y >= y1)`.

There are `n x n` cells, and each of them is getting compared to `(n - x1) * y1 + x * (n - y1)` cells. This gives an overall time complexity of `O(n^4)` and is making my algorithm very slow. Can you improve the scalability of this algorithm?

Here's a visualisation:

Here's some python/pseudocode:

``````INVALID = -1

class Pair:
def __init__(self, x, y, value):
self.x = x
self.y = y
self.value = value

def get_ordered_pairs(matrix):
pairs = []
for i in range(len(matrix) * len(matrix[0])):
pairs.append(get_max_pair(matrix))
return pairs

def get_max_pair(matrix):
max_pair = Pair(INVALID, INVALID, INVALID)
for x in range(len(matrix)):
for y in range(len(matrix[x])):
if matrix[x][y] > max_pair.value:
max_pair = Pair(x, y, matrix[x][y])

# invalidate entire row
for y in range(len(matrix[0])):
matrix[max_pair.x][y] = INVALID

# invalidate entire column
for x in range(len(matrix)):
matrix[x][max_pair.y] = INVALID

# invalidate top right
for x in range(max_pair.x + 1, len(matrix)):
for y in range(max_pair.y):
matrix[x][y] = INVALID

# invalidate bottom left
for y in range(max_pair.y + 1, len(matrix[0])):
for x in range(max_pair.x):
matrix[x][y] = INVALID

return max_pair
``````
• What is the typical / worst cast sparsity of the matrix? In other words, typically / in the worst case, what percentage of the matrix elements are non-zero (or, of interest) ? Dec 18, 2020 at 10:08
• @rwong About 1% of cells have non-zero values when n = 1000. Dec 18, 2020 at 11:05
• You don’t take any advantage of the sparsity. Record the indices of non-zero entries and use that and your code will run at least 100 times faster. Dec 19, 2020 at 19:21

I apologize for initially missing the additional requirement aside from just sorting values. I think the requirements are still a bit unclear but based on the code and visual plus some assumptions, think I mostly follow.

Essentially each cell has two areas that it can 'dominate': the submatrix where it is the bottom-left cell and the one where it is the top-right cell. We can define this behavior in a class `Cell` (which I think is better name for this than `Pair`) like so:

``````class Cell:
def __init__(self, x, y, value):
self.x = x
self.y = y
self.value = value

def in_lower_left(self, other):
return self.x >= other.x and self.y <= other.y

def in_upper_right(self, other):
return self.x <= other.x and self.y >= other.y

def dominates(self, other):
if self.value < other.value:
return False

return self.in_lower_left(other) or self.in_upper_right(other)

def __eq__(self, other):
return self.value == other.value

def __lt__(self, other):
if self.value == other.value:
return other.dominates(self)
else:
return self.value < other.value

def __repr__(self):
return f"({self.x},{self.y}): {self.value}"
``````

Now that we have that defined, the algorithm follows easily:

``````import bisect

data = [
[1, 0, 0, 0, 0],
[0, 3, 6, 0, 0],
[0, 0, 0, 0, 0],
[0, 0, 7, 0, 0],
[0, 0, 0, 0, 2]
]

""" only square matrices supported """
size = len(data)

def is_dominated(list, cell):
for item in list:
if item.dominates(cell):
return True

return False

def sort_cells_by_value(matrix):
sorted = []

for row in range(size):
for col in range(size):
cell = Cell(col, row, matrix[row][col])
bisect.insort(sorted, cell)

sorted.reverse()
return sorted

def filter_dominated(sorted_cells):
results = []

for cell in sorted_cells:
if not is_dominated(results, cell):
results.append(cell)

return results

sorted_cells = sort_cells_by_value(data)
results = filter_dominated(sorted_cells)

print(results)
``````

Now, it's been a while but the complexity for the two steps (m=n*n)

Sort: O(m log m) + Filter: O(m * n)

The filter complexity is based on the fact that you cannot have more than `N` results in the output (which I'm fairly certain is true.) If I'm not mistaken about that, you can exit the filter loop once the results are of length `N`. You might also be able to get fancy around keeping track of non-dominated cells. Once there are no more unoccupied, non-dominated cells in the results, you are done. These don't necessarily change the O-time but can have impact on the actual performance of a real program.

If you think I've made an error here or have questions, let me know.

To address Helena's comment about keeping only two dominating cells, here's an example matrix which I think shows this doesn't work in general:

The above shows the state of the filter after the first two dominating cells are found. We can see there are 3 areas still free, including a 4X4 middle area. By induction if we can create empty sub-matrices, then these can be further subdivided as well in the same manner which can be seen once the 7 cell is handled. There's some upper limit to the number to this for any given N that, if generalized, might be potentially useful for optimization.

• Thanks for pointing out I cannot have more than `N` results; this one-line change has brought it down to `O(n^3)`, which has reduced the execution time from 32 hours to 2 minutes; a thousand times faster for `n = 1000`. Dec 18, 2020 at 5:13
• @Irfan434 OK, I see. One thing to keep in mind is that for real execution times, things that we ignore for O times can be significant. For example you have 4 O(N) loops in `get_max_pair`. The factor of 4 doesn't count in O times but it does have an impact on real processing times. Additionally, the O times I give are worst case. The best case and average case are significantly better. In your approach, you always get the worst case O(n^3) time. Dec 18, 2020 at 14:14
• @rwong Wouldn't the best case be zero free zones when the largest value is in position (0, N-1) or (N-1, 0)? This is what I alluded to about keeping track of the covered cells. The question is whether you can do that without adding more work than you save (in general.) Dec 18, 2020 at 14:26
• @JimmyJames Turns out my solution was wrong after all. For some reason I assumed that there would always ever be two allowed subrectangles. Dec 18, 2020 at 17:30
• @Helena Yes but I think the approach of looking at remaining space as opposed to covered space has a lot of merit. Depending on the data it will outperform the approach here. So with some heuristics and a strategy pattern, it might close in on optimal. Too much for a SO question though. Thanks for helping me see my initial erorr. Dec 18, 2020 at 17:35

The following algorithm recursively computes the result by finding the maximum element in the current matrix, then subdividing in two sub-matrices (just like in Helena's answer) and repeating the same calculation on both sub-matrices. The results are collected in a list, which is sorted at the end.

The complexity of finding the maximum value in a matrix is O(n^2). We must repeat this a total of n times, so in order to find the unsorted result we have a complexity of O(n^3). Sorting at the end has a complexity of O(n * log n), so the total complexity is still O(n^3).

The following snippet is pseudo-code. It resembles Python, but you'll still have to adapt it to make it work.

``````def recursive_maximum(matrix):
if matrix.empty():
return []

# find the maximum inside this matrix -> O(n^2)
max_element_x, max_element_y = argmax(matrix)

top_left_submatrix = matrix[:max_element_x, :max_element_y] # note: max_element_x, max_element_y should not be included
bottom_right_submatrix = matrix[max_element_x+1:, max_element_y+1:]

# return a list containing this matrix's maximum and the results for both sub-matrices
# In the end, we'll have to call recursive_maximum at most n times on a non-empty matrix
return [matrix[max_element_x, max_element_y]] +
recursive_maximum(top_left_submatrix) +
recursive_maximum(bottom_right_submatrix)

unsorted_result = recursive_maximum(matrix)
output = sort(recursive_maximum(matrix)) # -> O(n * log n)
``````

First, whenever you pick one maximum, you remove the whole row containing it, so you cannot ever pick two maxima in the same row, so you can pick at most N values. This means the execution time would be at most O(N^3). Unfortunately, this just improves your estimate, and doesn't do anything to the actual time.

Here's what I do: For each row, we start by finding the maximum, and recording the maximum, its index, which range was searched (1 .. N), and which range should have been searched (also 1.. N initially). We put these six values into a struct, and put it into a priority queue. The "range that should be searched" gets changed every time we pick a maximum value.

As long as the queue is not empty: Extract the row with the largest maximum from the priority queue. If the range that should have been searched is empty, try again (the row has been removed from the priority queue, so it won't be searched again). Check if the maximum happened in the range that should have been searched. If not, search for the maximum in the range that should have been searched, and update the maximum value, its position, and the range that was searched, then put the row back into the priority queue and try again.

If the maximum happened within the range that should have been searched, then it is the maximum that you were looking for. Record the maximum. Go through all the rows in the priority queue, and update the range that should be searched.

How does this help? When you find a new maximum value, you do very little work, just updating at most N ranges. There is a good chance that the next maximum is still not removed from the search, so you can find it quickly. You only need to do extra work for rows that contained the largest remaining value, which is now excluded. I don't think it reduces the worst case, but it will reduce the time for typical cases a lot.

PS. If only one percent of the entries are non-zero, then you obviously store the indices of the non-zero values in each row first, which will reduce the time looking for the largest value in a row by a factor 100.

Algorithm, assuming we are only interested in maximum values > 0:

``````Define a structure "row" consisting of the following:
The index of the row
The "searched" range
The "shouldSearch" range
An array of all column indices in the "searched" range with non-zero values
The maximum value in the "searched" range
The column index of the maximum value

If there are any non-zero items:
Set "nonZero" to the range from the first to the last non-zero item.
Set "searched" and "shouldSearch" to "nonZero"
Find the maximum value in the shouldSearch range and its column index.
Add the "row" structure to the priority queue.

Create a priority queue containing rows. The order is:
Rows with higher maximum value have higher priority.
If maximum values are equal, rows where "searched" = "shouldSearch" have higher priority.

For each row:
Create a "row" structure.
Find the indices of non-zero items and store them

As long as the priority queue is not empty:
Get the "row" with highest priority and remove from the queue.
If the "shouldSearch" range is empty then loop again.
If the maximum index is not in the "shouldSearch" range:
Remove the indices of non-zero values outside "shouldSearch"
Loop again.
If the maximum index is in the "shouldSearch" range:
Record the maximum value.
For each "row" in the priority queue:
Update the "shouldSearch" range.
``````

How long does this all take, assuming each entry is nonzero with the same probability p? The initial setup takes O(n^2) steps. Every time we process a maximum we change up to n rows, total O(n^2).

Re-calculating a row maximum reduces the number of entries to be considered in the row by 1, therefore this happens at most pn^2 times and takes at most pn steps, for a total of p^2 n^3 steps. In practice this happens much less often because we only do this if the previous maximum of the row was found in an excluded range, and is larger than the currently highest value, and because finding a new maximum can substantially reduce the number of points that need checking, often by half.

PS. I implemented this in C++ and tested it on a full 10,000 x 10,000 matrix. I counted how many comparisons are made to calculate the maximum of a subrange of a row. The number was never more than 1.01 n^2. n^2 comparisons are needed to find the global maximum, and less than a million comparisons in total on top of that. Runtime < 1 second. Setting 100 million values using arc4random() took significantly more time than that.

• I upvoted this answer on the basis that it mentions an optimization insight that I didn't find in the other answers: for each row (and likewise for each column), at most one element can become the dominating cell, namely, the maximum value for that row (or column). The filtering can be done with `O(n)` space where `n` is the side length of the matrix, and each test/update op is `O(1)` (per cell). This puts an upper limit on `K` (the number of candidates needed to be sorted), which is `K <= n`. Dec 18, 2020 at 14:22

How do you make something more efficient?

• Do less work
• Organize better

Regarding doing less work: instead of 'x'ing out boxes and rescanning the whole matrix, handle only what's left after finding each max.

Each found max subdivides the remaining space into two smaller spaces, which can be handled independently of its larger containing matrix.  (And what's more, you're already basically computing the corners of these two smaller spaces.)

Regarding better organization: apply @JimmyJames suggestion to scan the matrix only once gathering all values into sorted order; then while traversing that sorted series, track the subdivisions of the space, which are then used to tell you whether or not to allow/output the next element in the series.

Coding left as exercises for the reader.

• Unfortunately this is still O(n^4) Dec 17, 2020 at 16:48