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

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.

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

Subroutine "addRow":
    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
    Call the "addRow" subroutine.

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"
        Call the "addRow" subroutine.
        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.

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.