This problem can be represented using graph theory and then solved by using linear programming techniques.
By treating all elements of the matrix as vertices the goal is to find a set of edges that minimizes the weighted path through the matrix under certain constraints (e.g. how slices can be built).
We can use scipy.optimize.linprog (which uses the Simplex algorithm) for solving the linear programming problem c.T @ x. Each element of the solution represents a possible connection between any of the N nodes (i.e. the solution vector has size N**2). The coefficients c determine the weight of a connection: that is the weight of the node connected to (the value of the matrix element), except for the first layer where we need add the starting weights as well.
We also need to apply several constraints in order to get a valid solution:
- The result must have exactly
N - 1 edges.
- Each node must not be connected to itself.
- Each row must connect to exactly one node except for the last row.
- Each node can connect to at most one other node (except for the ones of the last layer which must have zero connections).
- Each node must respect its possible successors (which are given by the "shape" of slices).
- For each of the layers 1 through (N-1) any node connected to must connect to the next layer.
These are quite a lot of pieces we must fit together and so the resulting code is a bit lengthy and might appear overwhelming at first. However when taking a closer look it should be possible to identify the single pieces and their structure (I tried to comment as much as possible + added several prints). So here goes the example code:
import numpy as np
from scipy.optimize import linprog
M = np.arange(9).reshape(3, 3) + 1
print('Matrix:')
print(M)
print('\n\n')
N = len(M)
# Compute all possible connections between nodes (1: possible, 0: forbidden).
pc = np.zeros(shape=(N**2, N**2), dtype=int)
# Connect to nodes below (except the last layer).
i = np.arange(N**2 - N)
pc[i, i + N] = 1
# Connect to left nodes (except the last layer and leftmost column).
pc[i, i + N - 1] = 1
pc[i[::N], i[::N] + N - 1] = 0
# Connect to left nodes (except the last layer and rightmost column).
r = i + N + 1
mask = r < N**2
pc[i[mask], r[mask]] = 1
r = r[N-1::N]
mask = mask[N-1::N]
pc[i[N-1::N][mask], r[mask]] = 0
print('Possible connections:')
print(pc)
print('\n\n')
# Coefficients for linear programming problem represent the weight of connections.
c = np.zeros(shape=(N**2, N**2), dtype=int)
# Add weights for connections.
c = np.tile(M.ravel(), (N**2, 1))
# Add additional weights for first layer.
c[:N] += M[0, :][:, None]
print('Coefficient matrix:')
print(c)
print('\n\n')
# === Add constraints ===
A_eq_1 = np.concatenate((
# Exactly N-1 connections.
np.ones(N ** 4, dtype=int)[None, :],
# No node can connect to itself.
np.diag([1] * N**2).flatten()[None, :]
), axis=0)
b_eq_1 = np.asarray([N - 1, 0], dtype=int)
print('Exactly N-1 connections and no self-connecting nodes:')
print(A_eq_1)
print(b_eq_1)
print('\n\n')
# Each layer connects to exactly one other node (except the last layer).
A_eq_2 = np.zeros((N, N ** 4), dtype=int)
for j in range(N):
A_eq_2[j, j * N**3:(j + 1)*N**3] = 1
b_eq_2 = np.ones(N, dtype=int)
b_eq_2[-1] = 0
print('Each layer connects to exactly one other node (except the last layer):')
print(A_eq_2)
print(b_eq_2)
print('\n\n')
# Each node connects to at most one other node (except the ones in the last layer).
N = N ** 2
A_ub_1 = np.zeros((N, N ** 2), dtype=int)
for j in range(N):
A_ub_1[j, j * N:j * N + N] = 1
b_ub_1 = np.ones(N, dtype=int)
b_ub_1[-1] = 0
print('Each node connects to at most one other node (except the ones in the last layer):')
print(A_ub_1)
print(b_ub_1)
print('\n\n')
# Each node respects its possible succesors (i.e. forbid all other connections).
A_eq_3 = np.zeros((N, N ** 2), dtype=int)
for j in range(N):
A_eq_3[j, j * N:j * N + N] = 1 - pc[j, :]
b_eq_3 = np.zeros(N, dtype=int)
print('Each node respects its possible succesors (i.e. forbid all other connections):')
print(A_eq_3)
print(b_eq_3)
print('\n\n')
# For the layers 1 through (N-1) each node connected to must connect to the next layer.
A_eq_4 = np.zeros((N, N ** 2), dtype=int)
for j in range(len(M), N-len(M)):
A_eq_4[j, j::N] = 1
A_eq_4[j, j*N:(j+1)*N] = -1
b_eq_4 = np.zeros(N, dtype=int)
print('For the layers 1 through (N-1) each node connected to must connect to the next layer:')
print(A_eq_4)
print(b_eq_4)
print('\n\n')
# Concatenate all constraints.
A_eq = np.concatenate([A_eq_1, A_eq_2, A_eq_3, A_eq_4])
b_eq = np.concatenate([b_eq_1, b_eq_2, b_eq_3, b_eq_4])
A_ub = np.concatenate([A_ub_1])
b_ub = np.concatenate([b_ub_1])
res = linprog(c.ravel(), A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=(0, 1))
print(res.success)
print(res.x.reshape(N, N)) # Edges.
The very last output is the result and it is of the form:
[[0. 0. 0. 1. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 1. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0.]]
Which tells us that in order to get the minimum path we should connect node 0 (row index) to node 3 (column index) as well as node 3 (row index) to node 6 (column index). This represents the path (or "slice") (1, 4, 7). We can deduce the path by starting at the first row and then walking through the graph as the edges point:
edges = res.x.reshape(N, N)
for i, r in enumerate(edges):
# Numerical instabilities can cause some elements to have very small values > 0.
if r.sum() > 0.5:
print('Connect {} -> {}'.format(i, r.argmax()))
path = [edges[:len(M)].ravel().argmax() // N]
while edges[path[-1]].max() > 0.5:
path.append(edges[path[-1]].argmax())
print('Path: ', path)
print('Elements: ', M.ravel()[path])
print('Path weight: ', M.ravel()[path].sum())
Final note
The above example code leaves lots of room for performance improvements. For example it considers all possible connections between nodes as solutions which scale as M.size**2. Though we constrain the possible connections the number of computations is still much larger than if we constrained it right from the beginning, by including it in the problem architecture. That means that instead of having M.size**2 coefficients we could go with just 2*(M.shape[0] - 1) + 3*(M.shape[1] - 2)*(M.shape[0] - 1) which scales only as M.size. Plus we can use a much smaller constraint matrix since we've built these constraints already in the architecture of the problem. Taking the above code example as a basis it should be possible to adapt it accordingly. Thus I'll conclude at this point and leave the implementation of possible performance improvements to the interested reader.
(The above implementation also only works on square matrices, generalization to non-square matrices should be straightforward though.)
