Rust

A port of @qwr's Python answer in Rust.

// Cargo.toml
// [dependencies]
// ndarray = "0.15.3"

// src/main.rs
use ndarray::{Array2, Axis, arr2};

fn main() {
    let a = arr2(&[
        [-3, -2,  3,  2, -2],
        [ 3, -2, -2, -5, -3],
        [ 0,  3,  1,  1, -3],
        [ 1,  1,  2,  1, -1],
        [ 1,  1,  1, -3,  3]
    ]);

    let n = a.nrows();
    let mut d = Array2::<i32>::zeros((n + 1, n));
    let mut c = Array2::<i32>::zeros((n + 1, n));
    let mut m = Array2::<i32>::zeros((n + 1, n));

    for j in 0..n {
        for i in 0..=n {
            c[[i, j]] = if i > 0 { c[[i - 1, j]] + a[[i - 1, j]] } else { 0 };
            d[[i, j]] = c[[i, j]] + if j > 0 { m[[i, j - 1]] } else { 0 };
        }

        for i in (0..=n).rev() {
            m[[i, j]] = if i == n { d[[i, j]] } else { d[[i, j]].max(m[[i + 1, j]]) };
        }
    }

    let last_column = d.index_axis(Axis(1), n - 1);
    let max_value = last_column.iter().max().unwrap(); // This is how you find the max value correctly

    println!("{:?}", d);
    println!("Max value of the last column: {:?}", max_value);
}

$ cargo build --release
$ cargo run --release
[[0, 2, 3, 8, 10],
 [-3, 0, 6, 10, 8],
 [0, -2, 4, 5, 2],
 [0, 1, 5, 6, -1],
 [1, 2, 7, 7, -2],
 [2, 3, 8, 4, -2]], shape=[6, 5], strides=[5, 1], layout=Cc (0x5), const ndim=2
Max value of the last column: 10

Edit:

Also considering the added condition \$K\$, the code is:

// src/main.rs

fn main() {
    let n = 5;
    let a: [[i32; 5]; 5] = [
        [-3, -2, 3, 2, -2],
        [3, -2, -2, -5, -3],
        [0, 3, 1, 1, -3],
        [1, 1, 2, 1, -1],
        [1, 1, 1, -3, 3],
    ];
    let k = 1;
    let mut d = vec![vec![vec![0; k + 1]; n]; n + 1];
    let mut c = vec![vec![0; n]; n + 1];
    let mut m = vec![vec![vec![0; k + 1]; n]; n];

    for j in 0..n {
        for i in 0..=n {
            c[i][j] = if i > 0 {
                c[i - 1][j] + a[i - 1][j]
            } else {
                0
            };

            for k in 0..=k {
                d[i][j][k] = c[i][j];
                if j > 0 {
                    d[i][j][k] += std::cmp::max(
                        d[i][j - 1][k],
                        std::cmp::max(
                            d[n][j - 1][k],
                            if i + 1 < n && k > 0 {
                                m[i + 1][j - 1][k - 1]
                            } else {
                                i32::MIN
                            },
                        ),
                    );
                }
            }
        }

        for i in (0..n).rev() {
            for k in 0..=k {
                m[i][j][k] = std::cmp::max(d[i][j][k], if i + 1 < n { m[i + 1][j][k] } else { i32::MIN });
            }
        }
    }

    for k in 0..=k {
        for i in 0..=n {
            for j in 0..n {
                print!("{} ", d[i][j][k]);
            }
            println!();
        }
        println!();
    }

    let mut max_val = i32::MIN;
    for i in 1..n {
        max_val = std::cmp::max(max_val, d[i][n - 1][k - 1]);
    }
    max_val = std::cmp::max(max_val, d[0][n - 1][k]);

    println!("{}", max_val);
}

Without the \$k\$ constraint, \$O(n^2)\$ is possible using a columnwise DP.

Define \$d[i,j]\$ as the maximum sum of columns \$0\$ to \$j\$, given some boundary line with a horizontal boundary above cell \$(i,j)\$, and \$c[i,j]\$ to be the column prefix sum up to row \$i\$. Then $$d[i,j] = \max_{i \le i' \le n} d[i',j-1] + c[i,j]$$

This would take linear time, but since max is associative, we can do it with prefix sums too: $$ m[i,j] = \max (d[i,j], m[i+1,j]) = \max_{i \le i' \le n} d[i',j] \\ d[i,j] = c[i,j] + m[i,j-1] \\ $$

As requested, here is python code:

import numpy as np
a = np.array(
    [[-3, -2,  3,  2, -2],
     [ 3, -2, -2, -5, -3],
     [ 0,  3,  1,  1, -3],
     [ 1,  1,  2,  1, -1],
     [ 1,  1,  1, -3,  3]])
n = a.shape[0]

d = np.zeros((n+1, n))
c = np.zeros((n+1, n))
m = np.zeros((n+1, n))

for j in range(n):
    for i in range(n+1):
        c[i,j] = c[i-1,j] + a[i-1,j] if i > 0 else 0
        d[i,j] = c[i,j] + (m[i,j-1] if j > 0 else 0)

    for i in range(n,-1,-1):
        m[i,j] = d[i,j] if i == n else max(d[i,j], m[i+1,j])


print(d, max(d[:,-1]))

Finally, the added condition \$K\$ is a third parameter of the DP counting down how many new horizontal lines are left to use. This requires keeping track if the horizontal boundary row has changed in the max operation. In your example, you don't count the very bottom \$i=n\$ horizontal boundary, so this requires a special case to handle. Since you don't count the very top \$i=0\$ boundary either, your final result could look at \$d[0,n-1,K]\$ for the top row and \$d[1:n,n-1,K-1]\$ for the rest below. You should verify that this implements your problem correctly.

(This may also work if you modify the problem to allow the boundary line to go up and down.)

$$d[i,j,k] = c[i,j] + \max(d[i,j-1,k], \max_{i < i' < n} d[i', j-1, k-1], d[n,j-1,k])$$

With max prefix sum: $$ m[i,j,k] = \max(d[i,j,k], m[i+1,j,k]) \\ d[i,j,k] = c[i,j] + \max(d[i,j-1,k], m[i+1,j-1,k-1], d[n,j-1,k]) $$

K = 1
d = np.zeros((n+1, n, K+1))
c = np.zeros((n+1, n))
m = np.zeros((n, n, K+1))

for j in range(n):
    for i in range(n+1):
        c[i,j] = c[i-1,j] + a[i-1,j] if i > 0 else 0

        for k in range(K+1):
            d[i,j,k] = c[i,j]
            if j > 0:
                d[i,j,k] += max(d[i,j-1,k],
                                d[n,j-1,k],
                                m[i+1,j-1,k-1] if i+1 < n and k > 0
                                else -np.inf)

    for i in range(n-1,-1,-1):
        for k in range(K+1):
            m[i,j,k] = max(d[i,j,k], m[i+1,j,k] if i+1 < n else -np.inf)


for k in range(K+1):
    print(d[:,:,k])

print(max(d[0,-1,K], max(d[1:n,-1,K-1])))

The final time complexity is \$O(n^2 k)\$.

Pure JavaScript:

const _staircase = (rowIndex, colIndex, prevRowIndex, maxRowIndex, k, sum, rowIndices) => {
  result.numIterations += 1;
  if (maxRowIndex === 0 || k === kMax) {
    rowIndices += (',' + prevRowIndex).repeat(width - colIndex);
    for (let j = colIndex; j < width; j++) sum += c[prevRowIndex][j];
    colIndex = width; // skip further iteration
  }
  if (colIndex === width) { // extract result
    if (sum > result.sum) [result.sum, result.rowIndices] = [sum, rowIndices];
    return;
  }
  const kInc = colIndex === 0 && rowIndex === height ? 0 : rowIndex !== prevRowIndex;
  _staircase(0, colIndex + 1, rowIndex, rowIndex, k + kInc, sum + c[rowIndex][colIndex], rowIndices + ',' + rowIndex);
  if (rowIndex < maxRowIndex) _staircase(rowIndex + 1, colIndex, prevRowIndex, maxRowIndex, k, sum, rowIndices);
};

Try it online!

The function is a closure that iterates the pre-computed static matrix c recursively. The pre-computed matrix holds cumulative column sums and has a row of zeros prepended.

The paths that are followed are stored by concatenating 1-indexed row numbers in rowIndices which is a simple delimited text string like ,5,5,5,1,0.

The function keeps track of how deep to go in each column by observing maxRowIndex. The first if() is a simple optimization that skips iteration of further columns when we're on the first row (that only contains zeros in c) and when we're out of horizontal lines to insert.

The harness is at the TIO link above. It constructs c, calls the recursive function, and shows the result.

Given the constraint that the partition line may only move up or right, this corresponds to the number of combination words defined by north-east lattice-paths: https://en.wikipedia.org/wiki/Lattice_path

For example, for n=3, we would write:

NNN
NNE
NEE
ENN
EEN
EEE

The total number of paths is: C(n+n,n) = C(2n,n)

This is the central binomial coefficient: https://en.wikipedia.org/wiki/Central_binomial_coefficient

This next section may be considered as speculation on a technicality.

By this constraint you can define that the left-side contains no entries while the right-side contains all entries:

NNN

Or that the left-side contains all entries and the right-side contains no entries.

EEE

These considers the following sub-cases:

• The maximal left-side partition is the entire matrix, because every entry is positive.
• The maximal left-side partition is empty, because the right-side partition contains negative integers for every entry. Subsequently, the entire matrix contains negative entries.

Calculate the total value of the matrix based on the left-hand side of the partition created by the NE path. Words of 'N' and 'E' string movements corresponds from the bottom-left to the top-right corner.

def calculate_partition_value(matrix, path):
    n = matrix.shape[0]
    row, col = n - 1, 0  # Start from the bottom-left corner
    total_value = 0

    # Track boundary created by the path
    boundary = [(row, col)]

    for move in path:
        if move == 'N':
            row -= 1  # Move up
        elif move == 'E':
            col += 1  # Move right
        boundary.append((row, col))

    for i in range(n):
        for j in range(n):
            if is_left_of_path((i, j), boundary):
                total_value += matrix[i, j]

    return total_value

Determine if a cell is on the left-hand side of the path based on the boundary coordinates by checking if the cell is below any segment, then return True. Otherwise, the cell is either right of the path or above, return False.

def is_left_of_path(cell, boundary):
    for k in range(1, len(boundary)):
        if cell[1] < boundary[k][1] and cell[0] > boundary[k][0]:
            return False
    return True

Sort matrices based on the partition value:

def sort_matrices_by_value(matrices_with_values):
    return sorted(matrices_with_values, key=lambda x: x[1])

And to provide a simpler example, create an nxn matrix with random integers between -1 and 1.

def randomize_matrix_entries(n):
    return np.random.randint(-1, 1+1, size=(n, n))

Ungolfed code:

import numpy as np
from itertools import combinations

def generate_ne_paths(n):
    steps = 'N' * n + 'E' * n
    return set([''.join(p) for p in combinations(steps, n)])

def randomize_matrix_entries(n):
    return np.random.randint(-1, 1+1, size=(n, n))

def calculate_partition_value(matrix, path):
    n = matrix.shape[0]
    row, col = n - 1, 0
    total_value = 0
    boundary = [(row, col)]

    for move in path:
        if move == 'N':
            row -= 1
        elif move == 'E':
            col += 1
        boundary.append((row, col))

    for i in range(n):
        for j in range(n):
            if is_left_of_path((i, j), boundary):
                total_value += matrix[i, j]

    return total_value

def is_left_of_path(cell, boundary):
    for k in range(1, len(boundary)):
        if cell[1] < boundary[k][1] and cell[0] > boundary[k][0]:
            return False
    return True

def sort_matrices_by_value(matrices_with_values):
    return sorted(matrices_with_values, key=lambda x: x[1])

n = 3
ne_paths = generate_ne_paths(n)
matrices_with_values = []

for path in ne_paths:
    matrix = randomize_matrix_entries(n)
    value = calculate_partition_value(matrix, path)
    matrices_with_values.append((matrix, value, path))

sorted_matrices = sort_matrices_by_value(matrices_with_values)

print("First sorted matrix (example):")
print("Matrix:\n", sorted_matrices[0][0])
print("Total Value:", sorted_matrices[0][1])
print("Path:", sorted_matrices[0][2])