2
votes

Say I have a NxN numpy matrix. I am looking for the fastest way of extracting all square chunks (sub-matrices) from this matrix. Meaning all CxC parts of the original matrix for 0 < C < N+1. The sub-matrices should correspond to contiguous rows/columns indexes of the original matrix. I want to achieve this in as little time as possible.

2
Should the submatrices correspond to contiguous rows/columns indexes of the original matrix (e.g. rows 0, 1, 2 and columns 0, 1, 2) or could they be in arbitrary order (e.g: rows 0, 4, 8 and rows 0, 2, 1)?rth
I have updated my post to answer your question.TweedBeetle

2 Answers

0
votes

You could use Numpy slicing,

import numpy as np

n = 20
x = np.random.rand(n, n)

slice_list = [slice(k, l) for k in range(0, n) for l in range(k, n)]

results = [x[sl,sl] for sl in slice_list]

avoiding loops in Numpy, is not a goal by itself. As long as you are being mindful about it, there shouldn't be much overhead.

0
votes

Tricky enough, but here is an example of extracting all MxM submatrices in a NxN matrix.

import numpy as NP
import numpy.random as RNG

P = N - M + 1
x = NP.arange(P).repeat(M)
y = NP.tile(NP.arange(M), P) + x

a = RNG.randn(N, N)
b = a[NP.newaxis].repeat(P, axis=0)
c = b[x, y]
d = c.reshape(P, M, N)
e = d[:, NP.newaxis].repeat(P, axis=1)
f = e[:, x, :, y]
g = f.reshape(P, M, P, M)
h = g.transpose(2, 0, 3, 1)

for i in range(0, P):
    for j in range(0, P):
        assert NP.equal(h[i, j], a[i:i+M, j:j+M]).all()