This task is a straightforward generalization of Deconvolution/1D to higher dimensions. For example, the one dimensional case would be applicable to audio signals, whereas two dimensions would pertain to images. Define the discrete convolution in 





d




{\displaystyle {\mathit {d}}}

 dimensions of two functions





H
,
F
:


Z


d


→

R



{\displaystyle H,F:\mathbb {Z} ^{d}\rightarrow \mathbb {R} }

 
taking 





d




{\displaystyle {\mathit {d}}}

-tuples of integers to real numbers as the function





G
:


Z


d


→

R



{\displaystyle G:\mathbb {Z} ^{d}\rightarrow \mathbb {R} }


also taking 





d




{\displaystyle {\mathit {d}}}

-tuples of integers to reals and satisfying





G
(

n

0


,
…
,

n

d
−
1


)
=

∑


m

0


=
−
∞


∞


…

∑


m

d
−
1


=
−
∞


∞


F
(

m

0


,
…
,

m

d
−
1


)
H
(

n

0


−

m

0


,
…
,

n

d
−
1


−

m

d
−
1


)


{\displaystyle G(n_{0},\dots ,n_{d-1})=\sum _{m_{0}=-\infty }^{\infty }\dots \sum _{m_{d-1}=-\infty }^{\infty }F(m_{0},\dots ,m_{d-1})H(n_{0}-m_{0},\dots ,n_{d-1}-m_{d-1})}


for all 





d




{\displaystyle {\mathit {d}}}

-tuples of integers 



(

n

0


,
…
,

n

d
−
1


)
∈


Z


d




{\displaystyle (n_{0},\dots ,n_{d-1})\in \mathbb {Z} ^{d}}

. Assume






F




{\displaystyle {\mathit {F}}}

 and 





H




{\displaystyle {\mathit {H}}}

 (and therefore 





G




{\displaystyle {\mathit {G}}}

) are non-zero over only a finite domain bounded by the origin, hence possible to represent as finite multi-dimensional arrays or nested lists 





f




{\displaystyle {\mathit {f}}}

, 





h




{\displaystyle {\mathit {h}}}

, and 





g




{\displaystyle {\mathit {g}}}

.

For this task, implement a function (or method, procedure, subroutine, etc.) deconv to perform deconvolution (i.e., the inverse of convolution) by solving for 





h




{\displaystyle {\mathit {h}}}

 given 





f




{\displaystyle {\mathit {f}}}

 and 





g




{\displaystyle {\mathit {g}}}

. (See Deconvolution/1D for details.)

 The function should work for 





g




{\displaystyle {\mathit {g}}}

 of arbitrary length in each dimension (i.e., not hard coded or constant) and 





f




{\displaystyle {\mathit {f}}}

 of any length up to that of 





g




{\displaystyle {\mathit {g}}}

 in the corresponding dimension.
 The deconv function will need to be parameterized by the dimension 





d




{\displaystyle {\mathit {d}}}

 unless the dimension can be inferred from the data structures representing 





g




{\displaystyle {\mathit {g}}}

 and 





f




{\displaystyle {\mathit {f}}}

.
 There may be more equations than unknowns. If convenient, use a function from a library that finds the best fitting solution to an overdetermined system of linear equations (as in the Multiple regression task). Otherwise, prune the set of equations as needed and solve as in the Reduced row echelon form task.
 Debug your solution using this test data, of which a portion is shown below. Be sure to verify both that the deconvolution of 





g




{\displaystyle {\mathit {g}}}

 with 





f




{\displaystyle {\mathit {f}}}

 is 





h




{\displaystyle {\mathit {h}}}

 and that the deconvolution of 





g




{\displaystyle {\mathit {g}}}

 with 





h




{\displaystyle {\mathit {h}}}

 is 





f




{\displaystyle {\mathit {f}}}

. Display the results in a human readable form for the three dimensional case only.
dimension 1:

h: [-8, 2, -9, -2, 9, -8, -2]
f: [ 6, -9, -7, -5]
g: [-48, 84, -16, 95, 125, -70, 7, 29, 54, 10]

dimension 2:

h: [
      [-8, 1, -7, -2, -9, 4], 
      [4, 5, -5, 2, 7, -1], 
      [-6, -3, -3, -6, 9, 5]]
f: [
      [-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]]
g: [
      [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]]

dimension 3:

h: [
      [[-6, -8, -5, 9], [-7, 9, -6, -8], [2, -7, 9, 8]], 
      [[7, 4, 4, -6], [9, 9, 4, -4], [-3, 7, -2, -3]]]
f: [
      [[-9, 5, -8], [3, 5, 1]], 
      [[-1, -7, 2], [-5, -6, 6]], 
      [[8, 5, 8],[-2, -6, -4]]]
g: [
      [
         [54, 42, 53, -42, 85, -72], 
         [45, -170, 94, -36, 48, 73], 
         [-39, 65, -112, -16, -78, -72], 
         [6, -11, -6, 62, 49, 8]], 
      [
         [-57, 49, -23, 52, -135, 66], 
         [-23, 127, -58, -5, -118, 64], 
         [87, -16, 121, 23, -41, -12], 
         [-19, 29, 35, -148, -11, 45]], 
      [
         [-55, -147, -146, -31, 55, 60], 
         [-88, -45, -28, 46, -26, -144], 
         [-12, -107, -34, 150, 249, 66], 
         [11, -15, -34, 27, -78, -50]], 
      [
         [56, 67, 108, 4, 2, -48], 
         [58, 67, 89, 32, 32, -8], 
         [-42, -31, -103, -30, -23, -8],
         [6, 4, -26, -10, 26, 12]]]

import std.stdio, std.conv, std.algorithm, std.numeric, std.range;
 
class M(T) {
    private size_t[] dim;
    private size_t[] subsize;
    private T[] d;
 
    this(size_t[] dimension...) pure nothrow {
        setDimension(dimension);
        d[] = 0; // init each  entry to zero;
    }
 
    M!T dup() {
        auto m = new M!T(dim);
        return m.set1DArray(d);
    }
 
    M!T setDimension(size_t[] dimension ...) pure nothrow {
        foreach (const e; dimension)
            assert(e > 0, "no zero dimension");
        dim = dimension.dup;
        subsize = dim.dup;
        foreach (immutable i; 0 .. dim.length)
            subsize[i] = reduce!q{a * b}(1, dim[i + 1 .. $]);
        immutable dlength = dim[0] * subsize[0];
        if (d.length != dlength)
            d = new T[dlength];
        return this;
    }
 
    M!T set1DArray(in T[] t ...) pure nothrow @nogc {
        auto minLen = min(t.length, d.length);
        d[] = 0;
        d[0 .. minLen] = t[0 .. minLen];
        return this;
    }
 
    size_t[] seq2idx(in size_t seq) const pure nothrow {
        size_t acc = seq, tmp;
        size_t[] idx;
        foreach (immutable e; subsize) {
            idx ~= tmp = acc / e;
            acc = acc - tmp * e; // same as % (mod) e.
        }
        return idx;
    }
 
    size_t size() const pure nothrow @nogc @property {
        return d.length;
    }
 
    size_t rank() const pure nothrow @nogc @property {
        return dim.length;
    }
 
    size_t[] shape() const pure nothrow @property { return dim.dup; }
 
    T[] raw() const pure nothrow @property { return d.dup; }
 
    bool checkBound(size_t[] idx ...) const pure nothrow @nogc {
        if (idx.length > dim.length)
            return false;
        foreach (immutable i, immutable dm; idx)
            if (dm >= dim[i])
                return false;
        return true;
    }
 
    T opIndex(size_t[] idx ...) const pure nothrow @nogc {
        assert(checkBound(idx), "OOPS");
        return d[dotProduct(idx, subsize)];
    }
 
    T opIndexAssign(T v, size_t[] idx ...) pure nothrow @nogc {
        assert(checkBound(idx), "OOPS");
        d[dotProduct(idx, subsize)] = v;
        return v;
    }
 
    override bool opEquals(Object o) const pure {
        const rhs = to!(M!T)(o);
        return dim == rhs.dim && d == rhs.d;
    }
 
    int opApply(int delegate(ref size_t[]) dg) const {
        size_t[] yieldIdx;
        foreach (immutable i; 0 .. d.length) {
            yieldIdx = seq2idx(i);
            if (dg(yieldIdx))
                break;
        }
        return 0;
    }
 
    int opApply(int delegate(ref size_t[], ref T) dg) {
        size_t idx1d = 0;
        foreach (idx; this) {
            if (dg(idx, d[idx1d++]))
                break;
        }
        return 0;
    }
 
    // _this_ is h, rhs is f, output g.
    M!T convolute(M!T rhs) const pure nothrow {
        auto dm = dim.dup;
        dm[] += rhs.dim[] - 1;
        M!T m = new M!T(dm); // dm will be reused as m's idx.
        auto bound = m.size;
        foreach (immutable i; 0 .. d.length) {
            auto thisIdx = seq2idx(i);
            foreach (immutable j; 0 .. rhs.d.length) {
                dm[] = thisIdx[] + rhs.seq2idx(j)[];
                immutable midx1d = dotProduct(dm, m.subsize);
                if (midx1d < bound)
                    m.d[midx1d] += d[i] * rhs.d[j];
                else
                    break; // Bound reach, OK to break.
            }
        }
        return m;
    }
 
    // _this_ is g, rhs is f, output is h.
    M!T deconvolute(M!T rhs) const pure nothrow {
        auto dm = dim.dup;
        foreach (i, e; dm)
            assert(e + 1 > rhs.dim[i],
                   "deconv : dimensions is zero or negative");
        dm[] -= (rhs.dim[] - 1);
        auto m = new M!T(dm); // dm will be reused as rhs' idx.
 
        foreach (immutable i; 0 .. m.size) {
            auto idx = m.seq2idx(i);
            m.d[i] = this[idx];
            foreach (immutable j; 0 .. i) {
                immutable jdx = m.seq2idx(j);
                dm[] = idx[] - jdx[];
                if (rhs.checkBound(dm))
                    m.d[i] -=  m.d[j] * rhs[dm];
            }
            m.d[i] /= rhs.d[0];
        }
        return m;
    }
 
    override string toString() const pure { return d.text; }
}
 
auto fold(T)(T[] arr, ref size_t[] d) pure {
    if (d.length == 0)
        d ~= arr.length;
 
    static if (is(T U : U[])) { // Is arr an array of arrays?
        assert(arr.length > 0, "no empty dimension");
        d ~= arr[0].length;
        foreach (e; arr)
            assert(e.length == arr[0].length, "Not rectangular");
        return fold(arr.reduce!q{a ~ b}, d);
    } else {
        assert(arr.length == d.reduce!q{a * b}, "Not same size");
        return arr;
    }
}
 
auto arr2M(T)(T a) pure {
    size_t[] dm;
    auto d = fold(a, dm);
    alias E = ElementType!(typeof(d));
    auto m = new M!E(dm);
    m.set1DArray(d);
    return m;
}
 
void main() {
    alias Mi = M!int;
    auto hh = [[[-6, -8, -5, 9], [-7, 9, -6, -8], [2, -7, 9, 8]],
               [[7, 4, 4, -6], [9, 9, 4, -4], [-3, 7, -2, -3]]];
    auto ff = [[[-9, 5, -8], [3, 5, 1]],[[-1, -7, 2], [-5, -6, 6]],
               [[8, 5, 8],[-2, -6, -4]]];
    auto h = arr2M(hh);
    auto f = arr2M(ff);
 
    const g = h.convolute(f);
 
    writeln("g == f convolute h ? ", g == f.convolute(h));
    writeln("h == g deconv f    ? ", h == g.deconvolute(f));
    writeln("f == g deconv h    ? ", f == g.deconvolute(h));
    writeln("         f = ", f);
    writeln("g deconv h = ", g.deconvolute(h));
}