21

I'm looking for a fast implementation for the following, I'll call it Position2D for lack of a better term:

Position2D[ matrix, sub_matrix ]

which finds the locations of sub_matrix inside matrix and returns the upper left and lower right row/column of a match.

For example, this:

Position2D[{
 {0, 1, 2, 3},
 {1, 2, 3, 4},
 {2, 3, 4, 5},
 {3, 4, 5, 6}
}, {
 {2, 3},
 {3, 4}
}]

should return this:

{
 {{1, 3}, {2, 4}},
 {{2, 2}, {3, 3}},
 {{3, 1}, {4, 2}} 
}

It should be fast enough to work quickly on 3000x2000 matrices with 100x100 sub-matrices. For simplicity, it is enough to only consider integer matrices.

2
  • 1
    Position[] gives you the top left corners. The bottom right corners are easily computed from that by just adding the dimensions of the searched for matrix to the top left corner.
    – Brian
    Dec 3, 2011 at 2:11
  • 1
    @brian please specify how you would use Position directly to find a specific submatrix. I don't think you can. Dec 3, 2011 at 15:21

4 Answers 4

29

Algorithm

The following code is based on an efficient custom position function to find positions of (possibly overlapping) integer sequences in a large integer list. The main idea is that we can first try to eficiently find the positions where the first row of the sub-matrix is in the Flatten-ed large matrix, and then filter those, extracting full sub-matrices and comparing to the sub-matrix of interest. This will be efficient for most cases except very pathological ones (those, for which this procedure would generate a huge number of potential position candidates, while the true number of entries of the sub-matrix would be much smaller. But such cases seem rather unlikely generally, and then, further improvements to this simple scheme can be made).

For large matrices, the proposed solution will be about 15-25 times faster than the solution of @Szabolcs when a compiled version of sequence positions function is used, and 3-5 times faster for the top-level implementation of sequence positions - finding function. The actual speedup depends on matrix sizes, it is more for larger matrices. The code and benchmarks are below.

Code

A generally efficient function for finding positions of a sub-list (sequence)

These helper functions are due to Norbert Pozar and taken from this Mathgroup thread. They are used to efficiently find starting positions of an integer sequence in a larger list (see the mentioned post for details).

Clear[seqPos];
fdz[v_] := Rest@DeleteDuplicates@Prepend[v, 0];
seqPos[list_List, seq_List] :=
  Fold[
     fdz[#1 (1 - Unitize[list[[#1]] - #2])] + 1 &, 
     fdz[Range[Length[list] - Length[seq] + 1] *
       (1 - Unitize[list[[;; -Length[seq]]] - seq[[1]]])] + 1,
     Rest@seq
  ] - Length[seq];

Example of use:

In[71]:= seqPos[{1,2,3,2,3,2,3,4},{2,3,2}]
Out[71]= {2,4}

A faster position-finding function for integers

However fast seqPos might be, it is still the major bottleneck in my solution. Here is a compiled-to-C version of this, which gives another 5x performance boost to my code:

seqposC = 
  Compile[{{list, _Integer, 1}, {seq, _Integer, 1}},
    Module[{i = 1, j = 1, res = Table[0, {Length[list]}], ctr = 0},
       For[i = 1, i <= Length[list], i++,
          If[list[[i]] == seq[[1]],
             While[j < Length[seq] && i + j <= Length[list] && 
                     list[[i + j]] == seq[[j + 1]], 
                j++
             ];
             If[j == Length[seq], res[[++ctr]] = i];
             j = 1;
          ]
       ];
       Take[res, ctr]
    ], CompilationTarget -> "C", RuntimeOptions -> "Speed"]

Example of use:

In[72]:= seqposC[{1, 2, 3, 2, 3, 2, 3, 4}, {2, 3, 2}]
Out[72]= {2, 4}

The benchmarks below have been redone with this function (also the code for main function is slightly modified )

Main function

This is the main function. It finds positions of the first row in a matrix, and then filters them, extracting the sub-matrices at these positions and testing against the full sub-matrix of interest:

Clear[Position2D];
Position2D[m_, what_,seqposF_:Automatic] :=
  Module[{posFlat, pos2D,sp = If[seqposF === Automatic,seqposC,seqposF]},
     With[{dm = Dimensions[m], dwr = Reverse@Dimensions[what]},
         posFlat = sp[Flatten@m, First@what];
         pos2D = 
            Pick[Transpose[#], Total[Clip[Reverse@dm - # - dwr + 2, {0, 1}]],2] &@
                {Mod[posFlat, #, 1], IntegerPart[posFlat/#] + 1} &@Last[dm];
         Transpose[{#, Transpose[Transpose[#] + dwr - 1]}] &@
            Select[pos2D,
                m[[Last@# ;; Last@# + Last@dwr - 1, 
                   First@# ;; First@# + First@dwr - 1]] == what &
            ]
     ]
  ];

For integer lists, the faster compiled subsequence position-finding function seqposC can be used (this is a default). For generic lists, one can supply e.g. seqPos, as a third argument.

How it works

We will use a simple example to dissect the code and explain its inner workings. This defines our test matrix and sub-matrix:

m  = {{0, 1, 2, 3}, {1, 2, 3, 4}, {2, 3, 4, 5}};
what = {{2, 3}, {3, 4}}; 

This computes the dimensions of the above (it is more convenient to work with reversed dimensions for a sub-matrix):

In[78]:= 
dm=Dimensions[m]
dwr=Reverse@Dimensions[what]

Out[78]= {3,4}
Out[79]= {2,2}

This finds a list of starting positions of the first row ({2,3} here) in the Flattened main matrix. These positions are at the same time "flat" candidate positions of the top left corner of the sub-matrix:

In[77]:= posFlat = seqPos[Flatten@m, First@what]
Out[77]= {3, 6, 9}

This will reconstruct the 2D "candidate" positions of the top left corner of a sub-matrix in a full matrix, using the dimensions of the main matrix:

In[83]:= posInterm = Transpose@{Mod[posFlat,#,1],IntegerPart[posFlat/#]+1}&@Last[dm]
Out[83]= {{3,1},{2,2},{1,3}}

We can then try using Select to filter them out, extracting the full sub-matrix and comparing to what, but we'll run into a problem here:

In[84]:= 
Select[posInterm,
   m[[Last@#;;Last@#+Last@dwr-1,First@#;;First@#+First@dwr-1]]==what&]

During evaluation of In[84]:= Part::take: Cannot take positions 3 through 4 
in {{0,1,2,3},{1,2,3,4},{2,3,4,5}}. >>

Out[84]= {{3,1},{2,2}}

Apart from the error message, the result is correct. The error message itself is due to the fact that for the last position ({1,3}) in the list, the bottom right corner of the sub-matrix will be outside the main matrix. We could of course use Quiet to simply ignore the error messages, but that's a bad style. So, we will first filter those cases out, and this is what the line Pick[Transpose[#], Total[Clip[Reverse@dm - # - dwr + 2, {0, 1}]], 2] &@ is for. Specifically, consider

In[90]:= 
Reverse@dm - # - dwr + 2 &@{Mod[posFlat, #, 1],IntegerPart[posFlat/#] + 1} &@Last[dm]
Out[90]= {{1,2,3},{2,1,0}}

The coordinates of the top left corners should stay within a difference of dimensions of matrix and a sub-matrix. The above sub-lists were made of x and y coordiantes of top - left corners. I added 2 to make all valid results strictly positive. We have to pick only coordiantes at those positions in Transpose@{Mod[posFlat, #, 1], IntegerPart[posFlat/#] + 1} &@Last[dm] ( which is posInterm), at which both sub-lists above have strictly positive numbers. I used Total[Clip[...,{0,1}]] to recast it into picking only at those positions at which this second list has 2 (Clip converts all positive integers to 1, and Total sums numbers in 2 sublists. The only way to get 2 is when numbers in both sublists are positive).

So, we have:

In[92]:= 
pos2D=Pick[Transpose[#],Total[Clip[Reverse@dm-#-dwr+2,{0,1}]],2]&@
           {Mod[posFlat,#,1],IntegerPart[posFlat/#]+1}&@Last[dm]
Out[92]= {{3,1},{2,2}} 

After the list of 2D positions has been filtered, so that no structurally invalid positions are present, we can use Select to extract the full sub-matrices and test against the sub-matrix of interest:

In[93]:= 
finalPos = 
  Select[pos2D,m[[Last@#;;Last@#+Last@dwr-1,First@#;;First@#+First@dwr-1]]==what&]
Out[93]= {{3,1},{2,2}}

In this case, both positions are genuine. The final thing to do is to reconstruct the positions of the bottom - right corners of the submatrix and add them to the top-left corner positions. This is done by this line:

In[94]:= Transpose[{#,Transpose[Transpose[#]+dwr-1]}]&@finalPos 
Out[94]= {{{3,1},{4,2}},{{2,2},{3,3}}}

I could have used Map, but for a large list of positions, the above code would be more efficient.

Example and benchmarks

The original example:

In[216]:= Position2D[{{0,1,2,3},{1,2,3,4},{2,3,4,5},{3,4,5,6}},{{2,3},{3,4}}]
Out[216]= {{{3,1},{4,2}},{{2,2},{3,3}},{{1,3},{2,4}}}

Note that my index conventions are reversed w.r.t. @Szabolcs' solution.

Benchmarks for large matrices and sub-matrices

Here is a power test:

nmat = 1000;
(* generate a large random matrix and a sub-matrix *)
largeTestMat = RandomInteger[100, {2000, 3000}];
what = RandomInteger[10, {100, 100}];
(* generate upper left random positions where to insert the submatrix *)    
rposx = RandomInteger[{1,Last@Dimensions[largeTestMat] - Last@Dimensions[what] + 1}, nmat];
rposy = RandomInteger[{1,First@Dimensions[largeTestMat] - First@Dimensions[what] + 1},nmat];
(* insert the submatrix nmat times *)
With[{dwr = Reverse@Dimensions[what]},
    Do[largeTestMat[[Last@p ;; Last@p + Last@dwr - 1, 
                     First@p ;; First@p + First@dwr - 1]] = what, 
       {p,Transpose[{rposx, rposy}]}]]

Now, we test:

In[358]:= (ps1 = position2D[largeTestMat,what])//Short//Timing
Out[358]= {1.39,{{{1,2461},{100,2560}},<<151>>,{{1900,42},{1999,141}}}}

In[359]:= (ps2 = Position2D[largeTestMat,what])//Short//Timing
Out[359]= {0.062,{{{2461,1},{2560,100}},<<151>>,{{42,1900},{141,1999}}}}

(the actual number of sub-matrices is smaller than the number we try to generate, since many of them overlap and "destroy" the previously inserted ones - this is so because the sub-matrix size is a sizable fraction of the matrix size in our benchmark).

To compare, we should reverse the x-y indices in one of the solutions (level 3), and sort both lists, since positions may have been obtained in different order:

In[360]:= Sort@ps1===Sort[Reverse[ps2,{3}]]
Out[360]= True

I do not exclude a possibility that further optimizations are possible.

15
  • 1
    I have a oneliner using PartitionMap. Unfortunately it's about 600 times slower than your solution. Better not post it less I make a fool of myself. Dec 4, 2011 at 0:07
  • 1
    @Sjoerd If it's short and not very difficult to follow, why not posting it? The question emphasized efficiency, that's true, but elegant and short code is appreciated by lots of people. Dec 5, 2011 at 11:39
  • 1
    @berniethejet Upon a cursory look, most of the code deals with positions, and so should be valid. The part that is different, should be rather straightforward to change / adjust. May 30, 2021 at 13:47
  • 1
    I've been trying out PositionIndex and Position with Alternatives and the former is faster with Reals and the latter with Integers. Dispatch is faster than PositionIndex for Reals, but still rather slow, I will have to try out these solutions here of yours, thanks. May 31, 2021 at 12:44
  • 1
    @berniethejet when using Position with reals, one has to be careful, because exact equality for reals may fail due to machine-precision errors, so one has to account for that and define some tolerance. The two Q/A-s I was able to find quickly are 1 and 2, and there are some interesting further references therein. May 31, 2021 at 16:26
13

This is my implementation:

position2D[m_, k_] :=
 Module[{di, dj, extractSubmatrix, pos},
  {di, dj} = Dimensions[k] - 1;
  extractSubmatrix[{i_, j_}] := m[[i ;; i + di, j ;; j + dj]];
  pos = Position[ListCorrelate[k, m], ListCorrelate[k, k][[1, 1]]];
  pos = Select[pos, extractSubmatrix[#] == k &];
  {#, # + {di, dj}} & /@ pos
 ]

It uses ListCorrelate to get a list of potential positions, then filters those that actually match. It's probably faster on packed real matrices.

5
  • +1 In your version, it is rare for the Position command to find more than a few matches. So almost no time is spent on the Select command -- most time is spent in the ListCorrelate command. This makes it faster than my code in most cases. However, I did find that for large matrices, your code does chew up a lot of memory and is more likely to swamp my computer than my code. (Also, my code can handle patterns! :P )
    – Simon
    Dec 3, 2011 at 12:14
  • @belisarius, Simon There's room for improvement here, so feel free to steal ideas and improve! I'll be away for some days, probably with no internet.
    – Szabolcs
    Dec 3, 2011 at 12:59
  • @Simon In my very limited test it was only 2x faster, which is probably not worth the extra memory use. Have you tried foating point arithmetic?
    – Szabolcs
    Dec 3, 2011 at 13:05
  • No, I haven't tested on non-integer matrices... -- No internet can be both heaven and hell at the same time. Enjoy your break(?)!
    – Simon
    Dec 3, 2011 at 13:07
  • +1 - your code set a high standard for speed, making it pretty hard to come up with something worth posting. Also, have a nice break! SO can be addictive, and I personally came to enjoy my time w/o internet as much as time spent here :) Dec 3, 2011 at 13:26
8

As per Leonid's suggestion here's my solution. I know it isn't very efficient (it's about 600 times slower than Leonid's when I timed it) but it's very short, rememberable, and a nice illustration of a rarely used function, PartitionMap. It's from the Developer package, so it needs a Needs["Developer`"] call first.

Given that, Position2D can be defined as:

Position2D[m_, k_] :=  Position[PartitionMap[k == # &, m, Dimensions[k], {1, 1}], True]

This only gives the upper-left coordinates. I feel the lower-right coordinates are really redundant, since the dimensions of the sub-matrix are known, but if the need arises one can add those to the output by prepending {#, Dimensions[k] + # - {1, 1}} & /@ to the above definition.

0
6

How about something like

Position2D[bigMat_?MatrixQ, smallMat_?MatrixQ] := 
 Module[{pos, sdim = Dimensions[smallMat] - 1}, 
  pos = Position[bigMat, smallMat[[1, 1]]];
  Quiet[Select[pos, (MatchQ[
      bigMat[[Sequence@@Thread[Span[#, # + sdim]]]], smallMat] &)], 
    Part::take]]

which will return the top left-hand positions of the submatrices. Example:

Position2D[{{0, 1, 2, 3}, {1, 2, 3, 4}, {2, 3, 4, 5}, {3, 5, 5, 6}}, 
   {{2, 3}, {3, _}}]
(* Returns: {{1, 3}, {2, 2}, {3, 1}} *) 

And to search a 1000x1000 matrix, it takes about 2 seconds on my old machine

SeedRandom[1]
big = RandomInteger[{0, 10}, {1000, 1000}];

Position2D[big, {{1, 1, _}, {1, 1, 1}}] // Timing
(* {1.88012, {{155, 91}, {295, 709}, {685, 661}, 
              {818, 568}, {924, 45}, {981, 613}}} *)
3
  • Thanks. This is still very slow in this, more general, case: r = RandomInteger[{0, 10}, {3000, 2000}]; p = Part[r, 400 ;; 500, 400 ;; 500]; Position2D[r, p] // Timing Dec 3, 2011 at 2:33
  • I guess using SameQ would speed up the comparisons (if you're not using patterns or approximate numerics). Or maybe there is some way to speed up multiple applications of Part. Finally, I guess if you only want a numerical Position2D, then the whole thing could be Compliled to "C"...
    – Simon
    Dec 3, 2011 at 2:45
  • Of course, looking at the other solutions, the main problem in mine is that the first pass gets too many hits, so the second has to make too many comparisons.
    – Simon
    Dec 3, 2011 at 12:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.