Re^3: Fast sliding submatrix sums with PDL (inspired by PWC 248 task 2)

[wlmb]

For large submatrices one can use a Fourier transforms to perform the convolution, as in

use v5.36;
use PDL;
use PDL::FFT;
my $matrix=pdl(shift); # Format "[[m11,m21...],[m21...]...]"
my $width=my $w=shift;
my $height=my $h=shift;
my $small=ones($w%2?$w:$w+1, $h%2?$h:$h+1);
$small->slice(-1).=0, ++$w unless $w%2; # zero row and/or column for e
+ven kernels
$small->slice([],-1).=0, ++$h unless $h%2;
my $kernel=kernctr($matrix, $small); #full kernel
my $result=$matrix->copy;
$result->fftconvolve($kernel);
say "$matrix $width $height -> ",
    $result->slice([floor(($width-1)/2),floor(-($width+1)/2)],
                   [floor(($height-1)/2),floor(-($height+1)/2)]);

It would be interesting to also compare it.

[Anonymous Monk]

Interesting. FFT for this task, who would have thought a few days ago. Looks like it runs in almost constant time, regardless of kernel size; and of course slower than special solutions. I had to add a couple of checks to accommodate for (unimportant) corner test case of submatrix with dims (4,3) and example matrix.

sub sms_WxH_PDL_fftconvolve ( $m, $w, $h ) {
    my ( $W, $H ) = $m-> dims;

    my $small = ones( $w%2?$w:$w+1, $h%2?$h:$h+1 );
    $small-> slice( -1 )    .= 0 unless $w%2; # zero row and/or column
+ for even kernels
    $small-> slice( [],-1 ) .= 0 unless $h%2;

    my $result = $m-> copy;
    $result = append( $result, 0 ) if $W == $w and not $W % 2;
    $result = append( $result-> transpose, 0 )-> transpose 
                                   if $H == $h and not $H % 2;

    my $kernel = kernctr( $result, $small ); #full kernel
    $result-> fftconvolve( $kernel );
    
    $result = $result-> slice( '0:-2','' ) if $W == $w and not $W % 2;
    $result = $result-> slice( '','0:-2' ) if $H == $h and not $H % 2;

    return $result-> slice(
        [floor(($w-1)/2),floor(-($w+1)/2)],
        [floor(($h-1)/2),floor(-($h+1)/2)]
    )
}

__END__

    Time (s) vs. N (NxN submatrix, PDL: Double D [1500,1500] matrix)

      +-----------------------------------------------------------+
      |    +      +      +      +       +      +      +      +    |
      |     A                                                     |
      |                                                           |
      |                                                           |
    2 |-+                                                       +-|
      |    B   B                        B  B      B               |
      |    BB         B      B      B          B         B   B    |
      |     B     B      B      B                     B           |
      |                         C   C                             |
      |                      C          C                         |
  1.5 |-+   A                              C                    +-|
      |                  C                                        |
      |                                        C                  |
      |                                                           |
      |               C                                           |
    1 |-+                                         C             +-|
      |     A                                                     |
      |           C                                   C           |
      |                                                           |
      |                                                           |
      |                                                  C        |
  0.5 |-+  A                                                    +-|
      |        C                                                  |
      |                                                      C    |
      |    A                                                      |
      |     D  D                                                  |
      |    CC     D   D  D   D  D   D   D  D   D  D   D  D        |
    0 |-+                                                    D  +-|
      |    +      +      +      +       +      +      +      +    |
      +-----------------------------------------------------------+
           0     200    400    600     800    1000   1200   1400
                         sms_WxH_PDL_conv2d    A
                    sms_WxH_PDL_fftconvolve    B
                           sms_WxH_PDL_lags    C
                        sms_WxH_PDL_sliding    D


  +------+-------+-------+-------+-------+
  | N    | A     | B     | C     | D     |
  +------+-------+-------+-------+-------+
  | 2    | 0.023 | 1.792 | 0.021 |       |
  | 4    | 0.065 | 1.854 | 0.036 |       |
  | 8    | 0.237 | 1.896 | 0.031 |       |
  | 12   | 0.534 | 1.865 | 0.049 |       |
  | 16   | 0.930 | 1.740 | 0.039 |       |
  | 20   | 1.432 | 1.812 | 0.057 |       |
  | 25   | 2.224 | 1.716 | 0.081 | 0.099 |
  | 100  |       | 1.899 | 0.380 | 0.107 |
  | 200  |       | 1.779 | 0.826 | 0.075 |
  | 300  |       | 1.836 | 1.143 | 0.089 |
  | 400  |       | 1.760 | 1.393 | 0.065 |
  | 500  |       | 1.826 | 1.560 | 0.047 |
  | 600  |       | 1.740 | 1.659 | 0.036 |
  | 700  |       | 1.818 | 1.620 | 0.049 |
  | 800  |       | 1.919 | 1.589 | 0.031 |
  | 900  |       | 1.917 | 1.484 | 0.026 |
  | 1000 |       | 1.789 | 1.260 | 0.008 |
  | 1100 |       | 1.880 | 1.068 | 0.010 |
  | 1200 |       | 1.745 | 0.823 | 0.010 |
  | 1300 |       | 1.841 | 0.562 | 0.008 |
  | 1400 |       | 1.841 | 0.273 | 0.003 |
  +------+-------+-------+-------+-------+