There are a few things in the question that I don't entire understand and seem contradictory, but I think I have two candidate solutions for you. If I missed any key components you were looking for, please feel free to update the question. Here are the constraints I followed:

• Input is "a 2D array (size N x M), let's call it U, where each cell contains a non-negative value K ≥ 0"
• "An algorithm or method that returns a set of non-overlapping rectangular boxes"
• Where "every cell in the original array U will have a corresponding number of "boxes" assigned to it"
• "Balance each box's size inversely with the average density inside it"

Here I have understood "box's size" to mean number of boxes assigned to that cell.

The two candidates I have for you are proc_array_unweighted and proc_array_weighted. show_plot is just a testing function to make some images so that you can visually assess the assignments to see if they meet your expectations.

The main bit of logic is to take the density array input, invert all the values so that little numbers are big and big numbers are little, scale it so that the greatest input cells get one box, then find a square number to chop up the smaller input cells into. Because this direct calculation makes some cells have a huge number of boxes, I also propose a weighted variant which further scales against the square root of the inverted cell values, which narrows the overall range of box counts.

import matplotlib.pyplot as plt
import numpy as np

def _get_nearest_square(num: int) -> int:
    # https://stackoverflow.com/a/49875384
    return np.pow(round(np.sqrt(num)), 2)

def proc_array_unweighted(arr: np.ndarray):
    scaled_arr = arr.copy()
    # Override any zeros so that we can invert the array
    scaled_arr[arr == 0] = 1
    # Invert the array
    scaled_arr = 1 / scaled_arr
    # Scale it so that the highest density cell always gets 1
    scaled_arr /= np.min(scaled_arr)
    # Find a square value to apply to each cell
    #   This guarantees that the area can be perfectly divided
    scaled_arr = np.vectorize(_get_nearest_square)(scaled_arr)

    return scaled_arr

def proc_array_weighted(arr: np.ndarray):
    scaled_arr = arr.copy()
    # Override any zeros so that we can invert the array
    scaled_arr[arr == 0] = 1
    # Invert the array, weighted against the square root
    #   This reduces the total range of output values
    scaled_arr = 1 / scaled_arr ** 0.5
    # Scale it so that the highest density cell always gets 1
    scaled_arr /= np.min(scaled_arr)
    # Find a square value to apply to each cell
    #   This guarantees that the area can be perfectly divided
    scaled_arr = np.vectorize(_get_nearest_square)(scaled_arr)

    return scaled_arr

def show_plot(arr: np.ndarray, other_arr1: np.ndarray, other_arr2: np.ndarray):
    fig, (ax1, ax2, ax3) = plt.subplots(1, 3)
    ax1.set_axis_off(); ax1.set_aspect(arr.shape[0] / arr.shape[1])
    ax2.set_axis_off(); ax2.set_aspect(arr.shape[0] / arr.shape[1])
    ax3.set_axis_off(); ax3.set_aspect(arr.shape[0] / arr.shape[1])
    for x_pos in range(arr.shape[1]):
        for y_pos in range(arr.shape[0]):
            ax1.text(
                (x_pos - 0.5) / arr.shape[1],
                (arr.shape[0] - y_pos - 0.5) / arr.shape[0],
                f'{arr[y_pos, x_pos]}',
                horizontalalignment='center',
                verticalalignment='center',
                transform=ax1.transAxes
            )
            for ax, arrsub in (
                (ax2, other_arr1),
                (ax3, other_arr2)
            ):
                ax.add_patch(plt.Rectangle(
                    (x_pos / arr.shape[1], y_pos / arr.shape[0]),
                    1 / arr.shape[1],
                    1 / arr.shape[0],
                    lw=2,
                    fill=False
                ))
                arr_dim = round(np.sqrt(arrsub[y_pos, x_pos]))
                for x_sub in range(arr_dim):
                    for y_sub in range(arr_dim):
                        # Draw sub-divides
                        top_leftx = x_pos / arr.shape[1] + x_sub / arr.shape[1] / arr_dim
                        top_lefty = y_pos / arr.shape[0] + (y_sub + 1) / arr.shape[0] / arr_dim
                        ax.add_patch(plt.Rectangle(
                            (top_leftx, 1 - top_lefty),
                            1 / arr.shape[1] / arr_dim,
                            1 / arr.shape[0] / arr_dim,
                            lw=1,
                            fill=False
                        ))
    plt.show()

def _main():
    test_points = [
        np.array([
            [1, 9, 1],
        ]),
        np.array([
            [0],
            [4],
            [1],
        ]),
        np.array([
            [1, 1, 1],
            [1, 1, 1],
            [1, 1, 1]
        ]),
        np.array([
            [1, 1, 1],
            [1, 8, 1],
            [1, 1, 1]
        ]),
        np.array([
            [1, 2, 1],
            [4, 8, 4],
            [1, 2, 1]
        ]),
        np.array([
            [ 1,   2,   4],
            [ 8,  16,  32],
            [64, 128, 256]
        ]),
        np.array([
            [1,  1, 1],
            [1, 72, 1],
            [1,  1, 1]
        ]),
        np.array([
            [1,  1,  1,  1, 1],
            [1, 72, 72, 72, 1],
            [1, 72, 72, 72, 1],
            [1, 72, 72, 72, 1],
            [1,  1,  1,  1, 1]
        ])
    ]

    for i, tp in enumerate(test_points):
        sol_unweighted = proc_array_unweighted(tp)
        sol_weighted = proc_array_weighted(tp)
        print('Array U:')
        print(tp)
        print('Array W (unweighted):')
        print(sol_unweighted)
        print('Array W (weighted):')
        print(sol_weighted)
        print('\n')
        show_plot(tp, sol_unweighted, sol_weighted)

if __name__ == '__main__':
    _main()

Here is the console print:

Array U:
[[1 9 1]]
Array W (unweighted):
[[9 1 9]]
Array W (weighted):
[[4 1 4]]


Array U:
[[0]
 [4]
 [1]]
Array W (unweighted):
[[4]
 [1]
 [4]]
Array W (weighted):
[[1]
 [1]
 [1]]


Array U:
[[1 1 1]
 [1 1 1]
 [1 1 1]]
Array W (unweighted):
[[1 1 1]
 [1 1 1]
 [1 1 1]]
Array W (weighted):
[[1 1 1]
 [1 1 1]
 [1 1 1]]


Array U:
[[1 1 1]
 [1 8 1]
 [1 1 1]]
Array W (unweighted):
[[9 9 9]
 [9 1 9]
 [9 9 9]]
Array W (weighted):
[[4 4 4]
 [4 1 4]
 [4 4 4]]


Array U:
[[1 2 1]
 [4 8 4]
 [1 2 1]]
Array W (unweighted):
[[9 4 9]
 [1 1 1]
 [9 4 9]]
Array W (weighted):
[[4 1 4]
 [1 1 1]
 [4 1 4]]


Array U:
[[  1   2   4]
 [  8  16  32]
 [ 64 128 256]]
Array W (unweighted):
[[256 121  64]
 [ 36  16   9]
 [  4   1   1]]
Array W (weighted):
[[16  9  9]
 [ 4  4  4]
 [ 1  1  1]]


Array U:
[[ 1  1  1]
 [ 1 72  1]
 [ 1  1  1]]
Array W (unweighted):
[[64 64 64]
 [64  1 64]
 [64 64 64]]
Array W (weighted):
[[9 9 9]
 [9 1 9]
 [9 9 9]]


Array U:
[[ 1  1  1  1  1]
 [ 1 72 72 72  1]
 [ 1 72 72 72  1]
 [ 1 72 72 72  1]
 [ 1  1  1  1  1]]
Array W (unweighted):
[[64 64 64 64 64]
 [64  1  1  1 64]
 [64  1  1  1 64]
 [64  1  1  1 64]
 [64 64 64 64 64]]
Array W (weighted):
[[9 9 9 9 9]
 [9 1 1 1 9]
 [9 1 1 1 9]
 [9 1 1 1 9]
 [9 9 9 9 9]]

And here are the plots:

Let me know if you have any questions, or if there is some feature you were hoping to see which is not presented.