Here we mean by submatrix any set of columns (not necessarily contiguous). I can model this as a linear MIP:
The results look like:
Notice I don’t do a max of a max here. I just let the model pick the elements that form the maximum objective. Automatically the solver will pick the maximum in each row.
Note this model only works when all entries are non-negative. If there are negative numbers things may not work anymore. But a simple preprocessing step can help. Just take the most negative value. Multiply by minus one and add to all entries of the matrix. Then apply the model.