Delete rows or columns of matrix containing invalid elements, such that a maximum number of valid elements is kept

Question

Originally posted in Stack Overflow but was told to post here.

Context: I am doing a PCA on a $M \times N \, (N \gg M)$ matrix with some invalid values located in the matrix. I cannot infer these values, so I need to remove all of them, which means I need to delete the whole corresponding row or column. Of course I want to keep the maximum amount of data. The invalid entries represent ~30% of data, but most of it is completely filled in a few lines, few of it is scattered in the rest of the matrix.

Some possible approaches:

• Similar to this problem, where I format my matrix such that valid data entries are equal to 1 and invalid entries to a huge negative number. However, all proposed solutions are of exponential complexity and my problem is simpler.

• Computing the ratio (invalid data / valid data) for each row or column, and deleting the highest ratio(s). Recompute the ratios for the sub-matrix and remove the highest(s) ratios. (not sure how many lines or columns we can remove safely in one step), and so on until there is no invalid data left. It seems like an okay solution, but I am unsure it always gives the optimal solution.

My guess is that it is a standard data analysis problem, but surprisingly I could not find a solution online.

Answer 1

Here is a solution based on back-tracking. You probably need to double check but it seemed ok on my side. Should take 1h on a 30x900 matrix, but you can use multi-threading in the bottom for loop to speed it up.

import numpy as np
from tqdm import tqdm


def hash_list(l1, l2):
    return str(sorted(l1)) + str(sorted(l2))  # could be faster


def max_matrix(a):
    # the process:
    # we construct list of rows and columns to select, which keep only valid values (a[rows, columns] = 1 everywhere)
    # at each step, for each other row/column, if it can be added in a valid way, we carry on deeper
    # when it's not possible to extend anymore
    # we backtrack to where we were and try to add the next possible row/column.
    # we store in res the visited valid rows/columns, which prevents most redundant computations
    res = {}
    M, N = a.shape

    def make_grow(a, sel_columns, sel_rows, other_columns, other_rows, d='right'):
        assert len(sel_columns) + len(other_columns) == M
        assert len(sel_rows) + len(other_rows) == N
        h = hash_list(sel_columns, sel_rows) + d
        if h in res: # we went through this already
            return

        # attempting to add a column
        if d == 'right':
            for i in other_columns:
                if 0 not in a[i, sel_rows]:
                    # we can append it: let's try to go left or right from there
                    other_columns.remove(i)
                    sel_columns.append(i)
                    make_grow(a, sel_columns, sel_rows, other_columns, other_rows, d='left')
                    make_grow(a, sel_columns, sel_rows, other_columns, other_rows, d='right')
                    
        # attempting to add a row
        else:
            for j in other_rows:
                if 0 not in a[sel_columns, j]:
                    # we can append it: let's try to go left or right from there
                    other_rows.remove(j)
                    sel_rows.append(j)
                    make_grow(a, sel_columns, sel_rows, other_columns, other_rows, d='right')
                    make_grow(a, sel_columns, sel_rows, other_columns, other_rows, d='left')

        # this position has been visited.
        res[h] = (sel_columns, sel_rows)

    # now launching from all possible first positions.
    for i in tqdm(range(M)):
        for j in range(N):
            if a[i, j] == 1:
                make_grow(a, [i], [j], [elt for elt in range(M) if elt != i], [elt for elt in range(N) if elt != j], d='right')
                make_grow(a, [i], [j], [elt for elt in range(M) if elt != i], [elt for elt in range(N) if elt != j], d='left')
            
    return res


if __name__ == '__main__':
    a = np.random.uniform(size=(30, 900))
    a[a > 0.3] = 1
    a[a <= 0.3] = 0
    a = a.astype(np.uint8)
    print(a)
    res = max_matrix(a)

    # just looking through the results to see the best configuration. (all valid configurations are stored btw)
    res_max = 0
    l1_max, l2_max = [], []

    for _, (l1, l2) in res.items():
        if len(l1) * len(l2) > res_max:
            res_max = len(l1) * len(l2)
            l1_max = l1
            l2_max = l2
    print(res_max, l1_max, l2_max)