Speeding up filling a matrix when approximating a double integral

Question

I want to approximate a double integral in the process of solving a partial-integro differential equation, to then put it, next to the differential part. The integral is (after some simplifications): \begin{equation}\label{int104} \mathcal{L}_Iu(s_i,v_j,\tau)=\int_{0}^{s_{\max}}\int_{v_j}^{v_{\max}} \frac{1}{s_i}u\left(z_1,z_2,\tau\right) \mathfrak{p}\left(z_1/s_i,z_2-v_j\right)\ dz_2\ dz_1 =\sum_{p=1}^{p=m-1}\sum_{q=1}^{q=n-1} M_{p,q}. \end{equation}

I extract the function $u$ by its average at the four closest (un-uniform) points in each cells of a 2D mesh, and write

\begin{equation}\label{int107} u(z_1,z_2,t)\simeq\frac{1}{4}\left(u(s_{p},v_{q},\tau)+u(s_{p},v_{q+1},\tau)+u(s_{p+1},v_{q},\tau)+u(s_{p+1},v_{q+1},\tau)\right), \end{equation} where $(z_1,z_2) \in [s_p,s_{p+1}]\times[v_q,v_{q+1}]$. Thus, we have \begin{align}\label{int108} M_{p,q}=&\frac{1}{4}\left(u(s_{p},v_{q},\tau)+u(s_{p},v_{q+1},\tau)+u(s_{p+1},v_{q},\tau)+u(s_{p+1},v_{q+1},\tau)\right)\\ &\times\mathfrak{P}\left(\frac{s_{p}}{s_i},\frac{s_{p+1}}{s_i},v_q-v_j,v_{q+1}-v_j\right), \end{align}

where I define the following function (from the 2D probability distribution function on one cell of the mesh)

\begin{align}\label{int109} \mathfrak{P}(A,B,C,D)=\int_{A}^{B}\int_{C}^{D} \mathfrak{p}(z_1,z_2)dz_2 dz_1. \end{align} To not distract from the main problem, I now give the following implementation:

ClearAll["Global`*"];

m = 32; n = 32; size = m*n;
r = 0.03; TT = 0.5; e = 100.;
xsmin = 0.; smax = xsmax = 4 e; ysmin = 0.; vmax = ysmax = 3.0;
{d1 = e/4., d2 = ysmax/200.};
r1 = 0.0025;
sleft = Max[0.5, Exp[-r1*TT]]*e; sright = e;
ksiMin = ArcSinh[(xsmin - sleft)/d1]; kint = (
 sright - sleft)/d1; ksiMax = kint + ArcSinh[(xsmax - sright)/d1];
ksi = Range[ksiMin, ksiMax, (ksiMax - ksiMin)/(m - 1)];
fun[ks_] := 
  Which[ksiMin <= ks < 0, sleft + d1*Sinh[ks], 0 <= ks <= kint, 
   sleft + d1*ks, kint < ks <= ksiMax, sright + d1*Sinh[ks - kint]];
nx = xgrid1 = Chop@Table[fun[ksi[[i]]], {i, 1, m}];
del2 = 1/(n - 1) (ArcSinh[ysmax/d2]); yg1 = 
 Table[(j - 1)*del2, {j, 1, n}]; ny = ygrid1 = Chop[d2 *Sinh[yg1]];
U[t_] = Flatten@Table[Subscript[u, i, j][t], {i, 1, m}, {j, 1, n}];
U1 = U[t];

And now I try to fill a matrix which provides the approximation of the double integral on our mesh

cf[a_, b_, c1_, 
   d_] = {(1/
     2) E^(-5 (c1 + d)) (E^(
       5 d) (-Erf[(5 (1 + c1 + 2 Log[a]))/(4 Sqrt[2])]
         + Erf[(5 (1 + c1 + 2 Log[b]))/(4 Sqrt[2])]) + 
      E^(5 c1) (-1 + Erf[(5 (1 + d + 2 Log[a]))/(4 Sqrt[2])] + 
         E^(13 + 5 d) (a^10 (Erf[(21 + 5 c1 + 10 Log[a])/(4 Sqrt[2])]
               - Erf[(21 + 5 d + 10 Log[a])/(4 Sqrt[2])]) + 
            b^10 (-Erf[(21 + 5 c1 + 10 Log[b])/(4 Sqrt[2])] + 
               Erf[(21 + 5 d + 10 Log[b])/(4 Sqrt[2])])) + 
         Erfc[(5 (1 + d + 2 Log[b]))/(4 Sqrt[2])]))};

Table[
   po1[o] = 
    Flatten@Table[
      Table[Mean@
        Flatten@Table[
          Subscript[u, i, j][t], {i, l, l + 1}, {j, k, k + 1}], {k, o,
         n - 1}], {l, 1, m - 1}]
   , {o, 1, n - 1}]; // AbsoluteTiming

mat02 = PadRight[
    ParallelTable[
     Table[{{sp = nx[[h1]], vq = ny[[h2]]};
        pt0 = 
         Flatten@Table[
           cf[(Rationalize@(nx[[i]]/sp)), (nx[[i + 1]]/sp), (ny[[j]] -
               vq), (ny[[j + 1]] - vq)], {i, m - 1}, {j, h2, n - 1}];
        Total@(pt0*po1[h2])}[[1]]
      , {h2, n - 1}]
     , {h1, 2, m - 1}]
    , {(m - 2), n}]; // AbsoluteTiming

mat02 = Last@CoefficientArrays[Flatten@mat02, U1]; // AbsoluteTiming
vec0 = SparseArray[{i_} -> 0, size];
mat01 = SparseArray@Table[vec0, {i, 1, n}];
integral = SparseArray@ArrayFlatten[{{mat01}, {mat02}, {mat01}}];

The above implementation is correct and give the correct matrix. But filling the matrix mat02 is very time consuming when $m$ and $n$ are high.

For example, when $m=100$, and $n=80$, it takes too much time.

Can anyone help me how to speed up filling such a matrix. Maybe, there is a way to Compile filling the matrix by Table!

Answer 1

You drive with handbreak on if you use symbolic computations for such a task. While this use of CoefficientArrays is very idiomatic and may greatly help to produce a correctly working prototype, it has to be optimized away for "production quality" code.

Compilation

In the following I compiled the four loops needed to assemble the system matrix into a CompiledFunction cL. The basic strategy is to create a matrix full L of zeros in the beginning and adding contributions of each cell into L from within loop.

cP = Block[{a, b, c, d},
   With[{code = 
      N[(1/2) E^(-5 (c + 
             d)) (E^(5 d) (-Erf[(5 (1 + c + 2 Log[a]))/(4 Sqrt[2])] + 
             Erf[(5 (1 + c + 2 Log[b]))/(4 Sqrt[2])]) + 
          E^(5 c) (-1 + Erf[(5 (1 + d + 2 Log[a]))/(4 Sqrt[2])] + 
             E^(13 + 
                 5 d) (a^10 (Erf[(21 + 5 c + 10 Log[a])/(4 Sqrt[2])] -
                    Erf[(21 + 5 d + 10 Log[a])/(4 Sqrt[2])]) + 
                b^10 (-Erf[(21 + 5 c + 10 Log[b])/(4 Sqrt[2])] + 
                   Erf[(21 + 5 d + 10 Log[b])/(4 Sqrt[2])])) + 
             Erfc[(5 (1 + d + 2 Log[b]))/(4 Sqrt[2])]))]},
    Compile[{{a, _Real}, {b, _Real}, {c, _Real}, {d, _Real}}, code]
    ]
   ];

cL = With[{cP = cP},
   Compile[{{s, _Real, 1}, {v, _Real, 1}},
    Block[{m, n, L},
     m = Length[s];
     n = Length[v];
     L = Table[0., {m n}, { m n}];
     Do[
      L[[m (i - 1) + j, m (p - 1) + q]] += 0.25 Times[
         cP[
          Compile`GetElement[s, p]/Compile`GetElement[s, i],
          Compile`GetElement[s, p + 1]/Compile`GetElement[s, i + 1],
          Compile`GetElement[v, q] - Compile`GetElement[v, j],
          Compile`GetElement[v, q + 1] - Compile`GetElement[v, j + 1]
          ]
         ]
      , {j, 1, n - 1}, {i, 2, m - 1}, {p, 2, m - 1}, {q, j, n - 1}];
     L
     ],
    CompilationTarget -> "C",
    CompilationOptions -> {"InlineCompiledFunctions" -> True},
    RuntimeOptions -> "Speed"
    ]
   ];

Let's make a test run with you given dimenension (m = 32 and n = 32):

nx = Threshold[Table[fun[ksi[[i]]], {i, 1, m}]];
ny = Threshold[d2*Sinh[yg1]]; 

L = cL[nx, ny]; // AbsoluteTiming // First

0.193937

Actually, I get a different matrix than yours. Maybe it is due to some differing numbering convention or maybe I got something else wrong. I mostly implemented cL from your $\LaTeX$ formulas. Still, this may show you how to proceed.

Parallelization

A parallelized version of this code could look like this

cLij = With[{cP = cP},
   Compile[{{s, _Real, 1}, {v, _Real, 
      1}, {i, _Integer}, {j, _Integer}},
    Block[{m, n, Lij},
     m = Length[s];
     n = Length[v];
     Lij = Table[0., {m n}];
     Do[
      Lij[[m (p - 1) + q]] += 0.25 cP[
         Compile`GetElement[s, p]/Compile`GetElement[s, i],
         Compile`GetElement[s, p + 1]/Compile`GetElement[s, i + 1],
         Compile`GetElement[v, q] - Compile`GetElement[v, j],
         Compile`GetElement[v, q + 1] - Compile`GetElement[v, j + 1]
         ],
      {p, 2, m - 1}, {q, j, n - 1}];
     Lij
     ],
    CompilationTarget -> "C",
    RuntimeAttributes -> {Listable},
    Parallelization -> True,
    CompilationOptions -> {"InlineCompiledFunctions" -> True},
    RuntimeOptions -> "Speed"
    ]
   ];

Application and timing example (quad core CPU):

{i, j} = Transpose@Tuples[{Range[2, m - 1], Range[1, n - 1]}];
Lij = cLij[nx, ny, i, j]; // AbsoluteTiming // First

0.05457

This appears to scale quite well with the numbers of processors. Note that this produces only the part of L corresponding to mat02 and not the full matrix L.

Further remarks

There is no point in using SparseArray, here; the final matrix will have about 0.5 density which is way too much for an efficient use of SparseArray. (That's typical for integral operators by the way.)

Something else that a noticed: You use Chop on list. Since it will replace values close to 0. by 0 (and not by 0.!), it will unpack arrays which can slow down later processing of these lists.