matrix convolution

Community Support
1

Hi - I’m trying to figure out how to do a convolution on a matrix applied to successive overlapping sub-windows.

Say I have a larger matrix

 

q)m4:zpad 4 4#“f”$til 16


where

zpad:{0,/:flip 0,/:(flip x,:0),:0} / pad w/0s all around


and a smaller convolution matrix to apply to all 3x3 sub-windows within m4:

q)m0:3 3#-1 -1 -1 -1 8 -1 -1 -1 -1

So I want to convolve m0 against each of these sub-window in m4:

±--------±----------±----------±--------+
|0 0 0 |0 0 0 |0 0 0 |0 0 0 |
|0 0 1f |0 1f 2f |1f 2f 3f |2f 3f 0 |
|0 4f 5f |4f 5f 6f |5f 6f 7f |6f 7f 0 |
±--------±----------±----------±--------+
|0 0 1f |0 1f 2f |1f 2f 3f | 2f 3f 0|
|0 4f 5f |4f 5f 6f |5f 6f 7f | 6f 7f 0|
|0 8f 9f |8f 9f 10f |9f 10f 11f |10f 11f 0|
±--------±----------±----------±--------+
|0 4f 5f| 4f 5f 6f| 5f 6f 7f| 6f 7f 0|
|0 8f 9f| 8f 9f 10f| 9f 10f 11f|10f 11f 0|
|0 12f 13f|12f 13f 14f|13f 14f 15f|14f 15f 0|
±--------±----------±----------±--------+
|0 8f 9f| 8f 9f 10f| 9f 10f 11f|10f 11f 0|
|0 12f 13f|12f 13f 14f|13f 14f 15f|14f 15f 0|
|0 0 0 | 0 0 0 | 0 0 0 | 0 0 0|
±--------±----------±----------±--------+

I’m thinking that since m4 is a list of lists, maybe I need to generate the set of indexes into m4, like this:

q)m4[0 1 2;0 1 2]
0 0 0
0 0f 1f
0 4f 5f
q)m4[0 1 2;1 2 3]
0 0 0
0 1 2
4 5 6


and so on for 16 sets, or use an iterator but have not come up with anything non-loopy.

Does anyone have any suggestions?

Thanks,

Devon

–

Devon McCormick, CFA

Quantitative Consultant

2

Here’s what I came up with.  It’s not too elegant and it assumes both matrixes are square and that the convolution matrix is 3x3 (because we pad the target matrix with one row and column of zeros on each side):

q)m0:3 3#-1 -1 -1 -1 8 -1 -1 -1 -1 / Edge detector convolver

q)zpad:{0,/:flip 0,/:(flip x,:0),:0} / Pad mat all around w/1 row & col of 0s on each side
q)m4:zpad 4 4#“f”$til 16 / Mat to convolute
q)shape:{(count x),count x[0]} / Shape of regular rectangular array
q)convo:{[x;y] sum raze xy} / 1D convolution
q)n0:count m0
q)n1:count m4
q)n2:(count m4)-n0-1
q)sd1:raze (n1til n0)+:til n0
q)idxm4:((n2n2),n0n0)#raze (raze flip sd1+:n1*til n2)+/:til n2
q)flip (n2,n2)#((raze m4) idxm4)convo: raze m0

The last line performs the convolution, giving the expected answer of 

-10 -9 -6 9
9 0 0 24
21 0 0 36
66 51 54 85

This works by building the indexes into m4 of each sub-window into a 16x9 matrix and convoluting each row against raze m0.  Any suggestions for improvement are welcome.

Thanks,

Devon

3

I suspect there isn’t anything wrong with your approach - you flatten the matrices and deal with the flattened indices which is a good idea.

Here’s an approach that doesn’t flatten but indexes at depth, it should also handle non-square matrices. 

q)f:{til[1+count-c]+:til c:count y};

q){count[a 0]cut(sum raze y*)@/:x ./:raze a:fx;y/::f[x 0;y 0]}[m4;m0]

-10 -9 -6 9

9 0 0 24

21 0 0 36

66 51 54 85

Terry

4

Convolutions are expensive operations, O(n^2),  it’s going to be hard to be much faster

If you are working with larger convolutions I would consider instead implementing the convolution as a multiplication in the Fourier domain (i.e., FFT each matrix, multiple and inverse back) which can drop complexity to O(n log n). Implementing this in kdb borrowing from the Kx whitepaper on Signal Processing, this would require padding your matrices to a power of 2  so you need to consider the benefits. 

5

Thanks for the help.  I may look at the Fourier transform if what I have is too slow but I want to complete the process all the way to the end before I worry about optimizing it.  Does anyone have something for reading JPGs or PNGs into q?  I have a kluge for now but it would be more seamless if I could read images directly.

6

Hi Devon,

.pgm and .ppm files have a very easy to parse format. You can make a system call to imagemagik which will convert your jpg/png:

pixels:`short$@[;4] system"convert file.jpg pgm:-"


https://imagemagick.org/script/command-line-processing.php#output


This code has many flaws but may still be of interest: https://github.com/rianoc/qCam

Rian