2

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.

TweedBeetle
  • 13
  • 4
  • 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 Jun 18 '15 at 19:15
  • I have updated my post to answer your question. – TweedBeetle Jun 18 '15 at 19:34

2 Answers2

0

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.

rth
  • 10,680
  • 7
  • 53
  • 77
0

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()
Quan Gan
  • 69
  • 1
  • 1