Convolution

Introduction

In J-forum Miller pointed to a problem in RosettaCode concerning de-convolution, the inverse of convolution.
Convolution is defined by


Let $F,G,H:\mathbb{Z}^d\rightarrow\mathbb{R}$ take $d$-tuples of integers to real numbers,
satisfying:

 $H(n_1, \dots, n_d)=\sum\limits_{m_1=-\infty}^{\infty}\dots\sum\limits_{m_d=-\infty}^{\infty}F(m_1, \dots, m_d)G(n_1-m_1, \dots, n_d-m_d)$

 then it is said that $H$ is the convolution of $F$ and $G$.

Coordinates of nouns

Let's take f and g as examples.

   [f=:2 5?@$10
6 9 0 6 3
2 9 1 7 3

   [g=:3 4?@$10
5 2 9 5
4 0 9 0
0 7 9 8
)

The coordinates of a noun is the set of coordinates of its atoms. The coordinate of an atom a in a noun N is the noun c such that a -: (< c){N .
The index of a noun is the order of its atoms: So I define

   indx=:i.@$
   crds=:(#: i.)@$

   (<"1@crds ,&< indx)g
+-----------------+---------+
|+---+---+---+---+|0 1  2  3|
||0 0|0 1|0 2|0 3||4 5  6  7|
|+---+---+---+---+|8 9 10 11|
||1 0|1 1|1 2|1 3||         |
|+---+---+---+---+|         |
||2 0|2 1|2 2|2 3||         |
|+---+---+---+---+|         |
+-----------------+---------+

   g -: (<"1 crds g) { g
1

As can be seen, the coordinates of a noun are completely determined by its shape, $ .
On the other hand, the shape of a noun is completely determined by its coordinates: it's equal to 1 less the maximum of coordinates:

$ -: <:@(>./)@(#: ,@i.)@$

Notice that there is a close connection between the coordinates and the index of a noun:

   (crds -: $ #: indx) g
1
   (indx -: $ #. crds) g
1

so one can be used in a definition for the other.

Coordinates of convolution

As can be seen from the definition, an atom of the convolution of 2 nouns is equal to the product of all atoms which coordinates sum up to the coordinate of the atom of the convolution.
Let's have a look how that works out for f and g . Boxing is done for overview reasons.

   <"2 f (<@(+"1)/&crds) g              NB. plane table of sum of coordinates
+-----------------+-----------------+-----------------+-----------------+-----------------+
|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|
||0 0|0 1|0 2|0 3|||0 1|0 2|0 3|0 4|||0 2|0 3|0 4|0 5|||0 3|0 4|0 5|0 6|||0 4|0 5|0 6|0 7||
|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|
||1 0|1 1|1 2|1 3|||1 1|1 2|1 3|1 4|||1 2|1 3|1 4|1 5|||1 3|1 4|1 5|1 6|||1 4|1 5|1 6|1 7||
|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|
||2 0|2 1|2 2|2 3|||2 1|2 2|2 3|2 4|||2 2|2 3|2 4|2 5|||2 3|2 4|2 5|2 6|||2 4|2 5|2 6|2 7||
|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|
+-----------------+-----------------+-----------------+-----------------+-----------------+
|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|
||1 0|1 1|1 2|1 3|||1 1|1 2|1 3|1 4|||1 2|1 3|1 4|1 5|||1 3|1 4|1 5|1 6|||1 4|1 5|1 6|1 7||
|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|
||2 0|2 1|2 2|2 3|||2 1|2 2|2 3|2 4|||2 2|2 3|2 4|2 5|||2 3|2 4|2 5|2 6|||2 4|2 5|2 6|2 7||
|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|
||3 0|3 1|3 2|3 3|||3 1|3 2|3 3|3 4|||3 2|3 3|3 4|3 5|||3 3|3 4|3 5|3 6|||3 4|3 5|3 6|3 7||
|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|+---+---+---+---+|
+-----------------+-----------------+-----------------+-----------------+-----------------+

   (<@,./.&, indx) f (<@(+"1)/&crds) g          NB. grouping of ravel (,) of indx of table according to sum 
                                                NB. stitching (,.) is done for overview reasons
+-+--+--+--+--+--+--+--+--+--+--+---+--+---+---+--+---+---+--+---+---+--+---+---+--+--+--+---+---+---+---+---+
|0| 1| 2| 3| 4| 5| 6| 7| 8| 9|10| 11|15| 19| 23|27| 31| 35|39| 43| 47|51| 55| 59|68|69|70| 71| 83| 95|107|119|
| |12|13|14|60|16|17|18|64|20|21| 22|26| 30| 34|38| 42| 46|50| 54| 58|  |111|115|  |80|81| 82| 94|106|118|   |
| |  |24|25|  |61|28|29|  |65|32| 33|37| 41| 45|49| 53| 57|  | 99|103|  |   |   |  |  |92| 93|105|117|   |   |
| |  |  |36|  |72|62|40|  |76|66| 44|48| 52| 56|  | 87| 91|  |110|114|  |   |   |  |  |  |104|116|   |   |   |
| |  |  |  |  |  |73|63|  |  |77| 67|  | 75| 79|  | 98|102|  |   |   |  |   |   |  |  |  |   |   |   |   |   |
| |  |  |  |  |  |84|74|  |  |88| 78|  | 86| 90|  |109|113|  |   |   |  |   |   |  |  |  |   |   |   |   |   |
| |  |  |  |  |  |  |85|  |  |  | 89|  | 97|101|  |   |   |  |   |   |  |   |   |  |  |  |   |   |   |   |   |
| |  |  |  |  |  |  |96|  |  |  |100|  |108|112|  |   |   |  |   |   |  |   |   |  |  |  |   |   |   |   |   |
+-+--+--+--+--+--+--+--+--+--+--+---+--+---+---+--+---+---+--+---+---+--+---+---+--+--+--+---+---+---+---+---+

From the definition of convolution, it is rather obvious that the shape of the convolution is 1 less than the sum of the shapes, since this is the maximum of the sum of the coordinates:

   [shps=: f ,:&$ g
2 5
3 4
   [s=: <:@+/ shps
4 8

However, if I reshape the coordinate sums accordingly, we get

   4 8 $ <@:>/.~  , f (<@(+"1)/&crds) g
+---+---+---+---+---+---+---+---+
|0 0|0 1|0 2|0 3|1 0|1 1|1 2|1 3|
|   |0 1|0 2|0 3|1 0|1 1|1 2|1 3|
|   |   |0 2|0 3|   |1 1|1 2|1 3|
|   |   |   |0 3|   |1 1|1 2|1 3|
|   |   |   |   |   |   |1 2|1 3|
|   |   |   |   |   |   |1 2|1 3|
|   |   |   |   |   |   |   |1 3|
|   |   |   |   |   |   |   |1 3|
+---+---+---+---+---+---+---+---+
|2 0|2 1|2 2|2 3|0 4|1 4|2 4|0 5|
|2 0|2 1|2 2|2 3|0 4|1 4|2 4|0 5|
|   |2 1|2 2|2 3|0 4|1 4|2 4|0 5|
|   |2 1|2 2|2 3|0 4|1 4|2 4|   |
|   |   |2 2|2 3|   |1 4|2 4|   |
|   |   |2 2|2 3|   |1 4|2 4|   |
|   |   |   |2 3|   |1 4|2 4|   |
|   |   |   |2 3|   |1 4|2 4|   |
+---+---+---+---+---+---+---+---+
|1 5|2 5|0 6|1 6|2 6|0 7|1 7|2 7|
|1 5|2 5|0 6|1 6|2 6|   |1 7|2 7|
|1 5|2 5|   |1 6|2 6|   |   |   |
|1 5|2 5|   |1 6|2 6|   |   |   |
|1 5|2 5|   |   |   |   |   |   |
|1 5|2 5|   |   |   |   |   |   |
+---+---+---+---+---+---+---+---+
|3 0|3 1|3 2|3 3|3 4|3 5|3 6|3 7|
|   |3 1|3 2|3 3|3 4|3 5|3 6|   |
|   |   |3 2|3 3|3 4|3 5|   |   |
|   |   |   |3 3|3 4|   |   |   |
+---+---+---+---+---+---+---+---+

And this obviously is not the order in which I want them. The first row should have the ones starting with 0 , the second row with 1 , etc. So I have to sort the (ravel of the) table first.

   4 8 $ <@:>/.~ /:~, f (<@(+"1)/&crds) g
+---+---+---+---+---+---+---+---+
|0 0|0 1|0 2|0 3|0 4|0 5|0 6|0 7|
|   |0 1|0 2|0 3|0 4|0 5|0 6|   |
|   |   |0 2|0 3|0 4|0 5|   |   |
|   |   |   |0 3|0 4|   |   |   |
+---+---+---+---+---+---+---+---+
|1 0|1 1|1 2|1 3|1 4|1 5|1 6|1 7|
|1 0|1 1|1 2|1 3|1 4|1 5|1 6|1 7|
|   |1 1|1 2|1 3|1 4|1 5|1 6|   |
|   |1 1|1 2|1 3|1 4|1 5|1 6|   |
|   |   |1 2|1 3|1 4|1 5|   |   |
|   |   |1 2|1 3|1 4|1 5|   |   |
|   |   |   |1 3|1 4|   |   |   |
|   |   |   |1 3|1 4|   |   |   |
+---+---+---+---+---+---+---+---+
|2 0|2 1|2 2|2 3|2 4|2 5|2 6|2 7|
|2 0|2 1|2 2|2 3|2 4|2 5|2 6|2 7|
|   |2 1|2 2|2 3|2 4|2 5|2 6|   |
|   |2 1|2 2|2 3|2 4|2 5|2 6|   |
|   |   |2 2|2 3|2 4|2 5|   |   |
|   |   |2 2|2 3|2 4|2 5|   |   |
|   |   |   |2 3|2 4|   |   |   |
|   |   |   |2 3|2 4|   |   |   |
+---+---+---+---+---+---+---+---+
|3 0|3 1|3 2|3 3|3 4|3 5|3 6|3 7|
|   |3 1|3 2|3 3|3 4|3 5|3 6|   |
|   |   |3 2|3 3|3 4|3 5|   |   |
|   |   |   |3 3|3 4|   |   |   |
+---+---+---+---+---+---+---+---+

Verb conv

Now everything is ready for the final verb conv .
But since I have to ravel the output anyway to sort it, I introduce the index-raveled and the coordinates-raveled:

   indxr=:,@i.@$
   crdsr=: (#: ,@i.)@$

   (<"1@crdsr ,:&< indxr)g
+-------------------------------------------------+
|+---+---+---+---+---+---+---+---+---+---+---+---+|
||0 0|0 1|0 2|0 3|1 0|1 1|1 2|1 3|2 0|2 1|2 2|2 3||
|+---+---+---+---+---+---+---+---+---+---+---+---+|
+-------------------------------------------------+
|0 1 2 3 4 5 6 7 8 9 10 11                        |
+-------------------------------------------------+

NB. we still have:

   (crdsr -: $ #: indxr) g
1
   (indxr -: $ #. crdsr) g
1

This leads to

conv =: 4 : 0
 t=. /: sc=. , x <@(+"1)/&crdsr y
 (x <:@+&$ y) $ sc +//.&(t&{) , x */ y
)

Let's check this from one of the examples from RosettaCode:

h2=:".;._2]0 :0
 _8  1 _7 _2 _9  4
  4  5 _5  2  7 _1
 _6 _3 _3 _6  9  5
)
f2=:".;._2]0 :0
 _5  2 _2 _6 _7
  9  7 _6  5 _7
  1 _1  9  2 _7
  5  9 _9  2 _5
 _8  5 _2  8  5
)
g2=:".;._2]0 :0
  40  _21  53   42  105    1  87   60   39 _28
 _92  _64  19 _167  _71  _47 128 _109   40 _21
  58   85 _93   37  101  _14   5   37  _76 _56
 _90 _135  60 _125   68   53 223    4  _36 _48
  78   16   7 _199  156 _162  29   28 _103 _10
 _62  _89  69  _61   66  193 _61   71   _8 _30
  48   _6  21   _9 _150  _22 _56   32   85  25
)

   g2 -: f2 conv h2
1

Multinomials

Unnoticed so far is that a convolution also gives the product of two multinomials.
A multinomial is for variables $x_1,x_2,\dots,x_n$ what a polynomial is for $x$:
$F(x_1, \dots, x_n)=\sum\limits_{m_1}\dots\sum\limits_{m_n} f(m_1 \dots, m_n)x_1^{m_1} \dots, x_n^{m_n}$ and like wise for $G$.
Obviously the product of 2 such multinomials is F conv G .

   T1 =:  (f2 (,&(, <@:,. crdsr)) h2) p."0"0 _ T0 =: <"(1) 20 2 ?.@$ 10
                                                        NB. T1 is equal to f2(T0) ,: g2(T0) (f2 and g2 applied to T0)

   (*/T1) -: ((, <@:,. crdsr) f2 conv3 h2) p. T0        NB. */T1 is equal to (f2 conv g2)(T0)
1

Deconvolution

Setting up the equations

Let's define deconvolution by y -: z deconv x if and only if z -: x conv y .
So suppose f and h are given and g =: h deconv f is to be determined.

   [f=: 2 5 $ 6 9 0 6 3 2 9 1 7 3
6 9 0 6 3
2 9 1 7 3

   [h=: 4 8 $ 30 57 72 141 72 60 57 15 34 85 95 233 95 128 89 15 8 78 139 238 135 138 102 24 0 14 81 104 130 92 83 24
30 57  72 141  72  60  57 15
34 85  95 233  95 128  89 15
 8 78 139 238 135 138 102 24
 0 14  81 104 130  92  83 24

We know from the definition of convolution the relation between the shapes of f, g and h. So I can immediately say:

   [sg =: h >:@-&$ f            NB. shape of g
3 4

and from that, with the knowledge of conv the the pairs of indexes of f and g can be determined

   $poi=:  ,<"1 ($f) ,"0/&(,@i.) sg             NB. pairing of indexes
120

Now each box in poi  represents a known atom in f and an unknown atom in g .
I also need the sum of the coordinates, and its grade up.

t=: /: sc=: , <@(+"1)/&(#: ,@i.)/ ($f),:sg

This I use to group the items of poi (displaying only the first 12 items).

   12{. sc <@:>/.&(t&{)poi
+---+---+---+---+---+---+---+---+---+---+---+---+
|0 0|0 1|0 2|0 3|1 3|2 3|3 3|4 3|0 4|0 5|0 6|0 7|
|   |1 0|1 1|1 2|2 2|3 2|4 2|   |5 0|1 4|1 5|1 6|
|   |   |2 0|2 1|3 1|4 1|   |   |   |5 1|2 4|2 5|
|   |   |   |3 0|4 0|   |   |   |   |6 0|5 2|3 4|
|   |   |   |   |   |   |   |   |   |   |6 1|5 3|
|   |   |   |   |   |   |   |   |   |   |7 0|6 2|
|   |   |   |   |   |   |   |   |   |   |   |7 1|
|   |   |   |   |   |   |   |   |   |   |   |8 0|
+---+---+---+---+---+---+---+---+---+---+---+---+

The next phrase needs some explanation. Let's present it first:

   12 {."1 T0 =: (<"0,h) ,:~ (]/:"1 {.)&.> (< , f) ({:@] ,: ({"1~ {.))&.>  sc <@|:@:>/.&(t&{) poi
+--+---+-----+-------+-------+-----+---+--+---+-------+-----------+---------------+
|0 |1 0|2 1 0|3 2 1 0|3 2 1 0|3 2 1|3 2|3 |4 0|5 4 1 0|6 5 4 2 1 0|7 6 5 4 3 2 1 0|
|6 |6 9|6 9 0|6 9 0 6|9 0 6 3|0 6 3|6 3|3 |6 2|6 9 2 9|6 9 0 2 9 1|6 9 0 6 2 9 1 7|
+--+---+-----+-------+-------+-----+---+--+---+-------+-----------+---------------+
|30|57 |72   |141    |72     |60   |57 |15|34 |85     |95         |233            |
+--+---+-----+-------+-------+-----+---+--+---+-------+-----------+---------------+

First of all, the former output is transposed in each box and then the first row in each box is the index of the unknown atom in g and the second row are the corresponding values in f by which these unknowns have to be multiplied and summed to give the result in the lower box, which is the value of h .
Notice the sorting by t which is done to have the upper boxes in the same order as the lower ones.
Also the sorting within the boxes to have the unknowns in the same order

So this leaves me with a set of (in this case 32) equations of (here 12) unknowns. There are several ways to solve these. I will present three, but first I will collect what I got so far for the verb deconv .

deconv =: 4 : 0
 sz  =. x >:@-&$ y                              NB. shape of z
 poi =.  ,<"1 ($y) ,"0/&(,@i.) sz               NB. pair of indexes
 t=. /: sc=. , <@(+"1)/&(#: ,@i.)/ ($y),:sz     NB. order of ,y
 (<"0,x) ,:~ (]/:"1 {.)&.> (<, y) ({:@] ,: ({"1~ {.))&.>  sc <@|:@:>/.&(t&{) poi        NB. set of boxed equations
)

Solving the equations 1

This uses the boxed result to solve the equations right away, using the (unproven) fact that there is always at least one equation with exactly one unknown.
The deconv-verb will be extended to

deconv1 =: 4 : 0
 sz  =. x >:@-&$ y                              NB. shape of z
 poi =.  ,<"1 ($y) ,"0/&(,@i.) sz               NB. pair of indexes
 t=. /: sc=. , <@(+"1)/&(#: ,@i.)/ ($y),:sz     NB. order of ,y
 T0=. (<"0,x) ,:~ (]/:"1 {.)&.> (<, y) ({:@] ,: ({"1~ {.))&.>  sc <@|:@:>/.&(t&{) poi   NB. set of boxed equations
 sz $ ({:/:{.) >{: solveq^:_ T0;a: 
)

So the set of boxed equations will be extended with a box which will hold the results.
Here is the commented solution.

solveq =: 3 : 0                                         NB. last box holds known unknowns, first one the boxed equations
 if. 0 = {:$>{.y do. y return. end.                             NB. done if no equations are left
 Y=. >{.y                                                       NB. Y is set of boxed equations
 a=. ,.&.>/ ~.({.@[,(%{:)~)&.>/ Y #"1~ t=. 1 = #@{.@>@{. Y      NB. select boxes with 1 unknown, leading to known unknowns
 b=. (-.t) #"1 Y                                                NB. remaining equations with at least 2 unknowns
 (a,"1~&.>{:y) ;~ (#"1~ *@#@{.@>@{.) |:(>a) slv"_ 1 |:b         NB. join a to last box, select non empty reduced equations
)

slv=: 4 : 0                                                     NB. x is set of known unknowns, y is set of boxed equations
 'x0 x1'=. x                                                    NB. x0 is index of unknown, x1 is value
 'y0 y1'=. y0 ['y0 y2'=: y              NB. y0 are indexes of unknown in equation, y1 multipliers, y2 result in equation
 c=. x0 ([ #~ e.) y0
 ixy=. (<c) e.~&.> x0;y0
 (y2 - +/*/ ixy #&> x1 ; y1) ;~ (y0 ; y1) (#~ -.)&> {:ixy       NB. reducing an equation by all the known unknowns
)

   g ; h deconv1 f
+-------+-------+
|5 2 9 5|5 2 9 5|
|4 0 9 0|4 0 9 0|
|0 7 9 8|0 7 9 8|
+-------+-------+

And for the examples from RosettaCode:

   h2 ; g2 deconv1 f2
+-----------------+-----------------+
|_8  1 _7 _2 _9  4|_8  1 _7 _2 _9  4|
| 4  5 _5  2  7 _1| 4  5 _5  2  7 _1|
|_6 _3 _3 _6  9  5|_6 _3 _3 _6  9  5|
+-----------------+-----------------+

   h3 -: g3 deconv1 f3
1

Solving the equations 2

The second way of solving the equations is done in a more natural way. Given the deconv1-verb T0 is transformed in a normal set of linear equations:

   [T1=: (,h) ,.~ (<12#0) (({:@])`({.@])`[})&> /:~ {.T0         NB. 12 is */ sz where sz is shape of result h deconv f
6 0 0 0 0 0 0 0 0 0 0 0  30
9 6 0 0 0 0 0 0 0 0 0 0  57
0 9 6 0 0 0 0 0 0 0 0 0  72
3 6 0 9 0 0 0 0 0 0 0 0 141
6 0 9 6 0 0 0 0 0 0 0 0  72
3 7 1 9 3 6 0 9 0 0 0 0  60
7 1 9 2 6 0 9 6 0 0 0 0  57
1 9 2 0 0 9 6 0 0 0 0 0  15
9 2 0 0 9 6 0 0 0 0 0 0  34
2 0 0 0 6 0 0 0 0 0 0 0  85
0 3 6 0 0 0 0 0 0 0 0 0  95
0 3 7 1 0 3 6 0 0 0 0 0 233
0 0 3 6 0 0 0 0 0 0 0 0  95
0 0 3 7 0 0 3 6 0 0 0 0 128
0 0 0 3 0 0 0 0 0 0 0 0  89
0 0 0 3 0 0 0 3 0 0 0 0  15
0 0 0 0 3 7 1 9 3 6 0 9   8
0 0 0 0 7 1 9 2 6 0 9 6  78
0 0 0 0 1 9 2 0 0 9 6 0 139
0 0 0 0 9 2 0 0 9 6 0 0 238
0 0 0 0 2 0 0 0 6 0 0 0 135
0 0 0 0 0 3 7 1 0 3 6 0 138
0 0 0 0 0 0 3 7 0 0 3 6 102
0 0 0 0 0 0 0 3 0 0 0 3  24
0 0 0 0 0 0 0 0 2 0 0 0   0
0 0 0 0 0 0 0 0 9 2 0 0  14
0 0 0 0 0 0 0 0 1 9 2 0  81
0 0 0 0 0 0 0 0 3 7 1 9 104
0 0 0 0 0 0 0 0 7 1 9 2 130
0 0 0 0 0 0 0 0 0 3 7 1  92
0 0 0 0 0 0 0 0 0 0 3 7  83
0 0 0 0 0 0 0 0 0 0 0 3  24

Since the number of rows is (far) more than the number of unknowns, these equations are mutually dependent.
To my surprise I did not find any standard solution in J to solve such equations. So I constructed one myself.

NB. y is a matrix representing a set of equations given by the multipliers (}:"1 y) and the results ({:"1 y).
NB. the result is a vector, or a matrix, the minimal set of independent equations

sle=: 3 : 0                                             NB. solve linear equations
 Y=. (#"1~ +./@:*)&.|: y                                NB. eliminate dependent equations
 Y=. \:~ rdeq^:(<:@{:@$)^:_ Y                           NB. reduce the equations and sort them in descending order
 if. 1 = -~/ $Y do. {:"1 Y return. end.                 NB. give the solution in the natural order
 Y                                                      NB. give the minimal set of independent equations
)

rdeq=: 3 : 0                                            NB. reduce linear equations
 y0=. {.y
 assert. (0 ~: +./ }:y0) +. 0 = {:y0                    NB. conflicting equations
 t=. ((%~<:)@{. , }. % {.) y {"1~ 1 i.~ * | y0          NB. coefficients to be used
 Y=. ({. ,~ }.) t (] - (*/ {.)) y                       NB. cleaning the column
 (#"1~ +./@:*)&.|: Y                                    NB. eliminate dependent equations
)

So now we can construct the verb deconv2 based on deconv .

deconv2 =: 4 : 0
 sz  =. x >:@-&$ y                                      NB. shape of z
 poi =.  ,<"1 ($y) ,"0/&(,@i.) sz                       NB. pair of indexes
 t=. /: sc=: , <@(+"1)/&(#: ,@i.)/ ($y),:sz             NB. order of ,y
 T0=. (<"0,x) ,:~ (]/:"1 {.)&.> (<, y) ({:@] ,: ({"1~ {.))&.>  sc <@|:@:>/.&(t&{) poi   NB. set of boxed equations
 T1=. (,x),.~(<0 #~ */sz) (({:@])`({.@])`[})&> {.T0     NB. set of linear equations
 sz $ sle T1
)

   h deconv2 f
5 2 9 5
4 0 9 0
0 7 9 8

   h (deconv2 -: deconv1) f
1

   h3 -: g3 deconv2 f3
1   

Solving the equations 3

As Hui pointed out, I could have used the primitive dyad %. for solving the equations. Then I get

deconv3 =: 4 : 0
 sz  =. x >:@-&$ y                                      NB. shape of z
 poi =.  ,<"1 ($y) ,"0/&(,@i.) sz                       NB. pair of indexes
 t=. /: sc=: , <@(+"1)/&(#: ,@i.)/ ($y),:sz             NB. order of ,y
 T0=. (<"0,x) ,:~ (]/:"1 {.)&.> (<, y) ({:@] ,: ({"1~ {.))&.>  sc <@|:@:>/.&(t&{) poi   NB. set of boxed equations
 T1=. (,x),.~(<0 #~ */sz) (({:@])`({.@])`[})&> {.T0     NB. set of linear equations
 sz $ 1e_8 round ({:"1 %. }:"1) T1
)
   h3 -: g3 deconv3 f3
1

Collected Definitions

indx=:i.@$
crds=:(#: i.)@$

indxr=:,@i.@$
crdsr=: (#: ,@i.)@$

conv =: 4 : 0
 t=. /: sc=. , x <@(+"1)/&crdsr y
 (x <:@+&$ y) $ sc +//.&(t&{) , x */ y
)

==============
deconv1 =: 4 : 0
 sz  =. x >:@-&$ y                              NB. shape of z
 poi =.  ,<"1 ($y) ,"0/&(,@i.) sz               NB. pair of indexes
 t=. /: sc=. , <@(+"1)/&(#: ,@i.)/ ($y),:sz     NB. order of ,y
 T0=. (<"0,x) ,:~ (]/:"1 {.)&.> (<, y) ({:@] ,: ({"1~ {.))&.>  sc <@|:@:>/.&(t&{) poi   NB. set of boxed equations
 sz $ ({:/:{.) >{: solveq^:_ T0;a: 
)

solveq =: 3 : 0                                         NB. last box holds known unknowns, first one the boxed equations
 if. 0 = {:$>{.y do. y return. end.                             NB. done if no equations are left
 Y=. >{.y                                                       NB. Y is set of boxed equations
 a=. ,.&.>/ ~.({.@[,(%{:)~)&.>/ Y #"1~ t=. 1 = #@{.@>@{. Y      NB. select boxes with 1 unknown, leading to known unknowns
 b=. (-.t) #"1 Y                                                NB. remaining equations with at least 2 unknowns
 (a,"1~&.>{:y) ;~ (#"1~ *@#@{.@>@{.) |:(>a) slv"_ 1 |:b         NB. join a to last box, select non empty reduced equations
)

slv=: 4 : 0                                                     NB. x is set of known unknowns, y is set of boxed equations
 'x0 x1'=. x                                                    NB. x0 is index of unknown, x1 is value
 'y0 y1'=. y0 ['y0 y2'=: y              NB. y0 are indexes of unknown in equation, y1 multipliers, y2 result in equation
 c=. x0 ([ #~ e.) y0
 ixy=. (<c) e.~&.> x0;y0
 (y2 - +/*/ ixy #&> x1 ; y1) ;~ (y0 ; y1) (#~ -.)&> {:ixy       NB. reducing an equation by all the known unknowns
)

==============
deconv2 =: 4 : 0
 sz  =. x >:@-&$ y                                      NB. shape of z
 poi =.  ,<"1 ($y) ,"0/&(,@i.) sz                       NB. pair of indexes
 t=. /: sc=: , <@(+"1)/&(#: ,@i.)/ ($y),:sz             NB. order of ,y
 T0=. (<"0,x) ,:~ (]/:"1 {.)&.> (<, y) ({:@] ,: ({"1~ {.))&.>  sc <@|:@:>/.&(t&{) poi   NB. set of boxed equations
 T1=. (,x),.~(<0 #~ */sz) (({:@])`({.@])`[})&> {.T0     NB. set of linear equations
 sz $ sle T1
)

sle=: 3 : 0                                             NB. solve linear equations
 Y=. (#"1~ +./@:*)&.|: y                                NB. eliminate dependent equations
 Y=. \:~ rdeq^:(<:@{:@$)^:_ Y                           NB. reduce the equations and sort them in descending order
 if. 1 = -~/ $Y do. {:"1 Y return. end.                 NB. give the solution in the natural order
 Y                                                      NB. give the minimal set of independent equations
)

rdeq=: 3 : 0                                            NB. reduce linear equations
 y0=. {.y
 assert. (0 ~: +./ }:y0) +. 0 = {:y0                    NB. conflicting equations
 t=. ((%~<:)@{. , }. % {.) y {"1~ 1 i.~ * | y0          NB. coefficients to be used
 Y=. ({. ,~ }.) t (] - (*/ {.)) y                       NB. cleaning the column
 (#"1~ +./@:*)&.|: Y                                    NB. eliminate dependent equations
)

==============
deconv3 =: 4 : 0
 sz  =. x >:@-&$ y                                      NB. shape of z
 poi =.  ,<"1 ($y) ,"0/&(,@i.) sz                       NB. pair of indexes
 t=. /: sc=: , <@(+"1)/&(#: ,@i.)/ ($y),:sz             NB. order of ,y
 T0=. (<"0,x) ,:~ (]/:"1 {.)&.> (<, y) ({:@] ,: ({"1~ {.))&.>  sc <@|:@:>/.&(t&{) poi   NB. set of boxed equations
 T1=. (,x),.~(<0 #~ */sz) (({:@])`({.@])`[})&> {.T0     NB. set of linear equations
 sz $ 1e_8 round ({:"1 %. }:"1) T1
)

RE Boss/J-blog/(De)Convolution (last edited 2010-03-14 15:08:39 by RE Boss)