# Speeding up filling a matrix when approximating a double integral

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): $$\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}.$$

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

$$\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),$$ 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

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!

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[
CompileGetElement[s, p]/CompileGetElement[s, i],
CompileGetElement[s, p + 1]/CompileGetElement[s, i + 1],
CompileGetElement[v, q] - CompileGetElement[v, j],
CompileGetElement[v, q + 1] - CompileGetElement[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[
CompileGetElement[s, p]/CompileGetElement[s, i],
CompileGetElement[s, p + 1]/CompileGetElement[s, i + 1],
CompileGetElement[v, q] - CompileGetElement[v, j],
CompileGetElement[v, q + 1] - CompileGetElement[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.

• Thanks for the very detailed response and comments. I have some inquiries. 1. It seems I do not have a C comipler in my machine, and thus I faced with the following error:CCompilerDriverCreateLibrary::nocomp: A C compiler cannot be found on your system. Please consult the documentation to learn how to set up suitable compilers. 2. After I remove the "CompilationTarget -> "C"", I still face with the following error: CmCompiledFunction::cfn: Numerical error encountered at instruction 1; proceeding with uncompiled evaluation. How did you run your implementation without any problems? Jun 8, 2018 at 12:34
• I think $Log[0.]$ could provide "Indeterminate". Because of this I used Chop[] in lieu of Threshold to print $-Infinity$ and could proceed with the Erf function. May you please revise your implementation or write some further tips to get rid of these errors? Please describe how we could do on Windows. Jun 8, 2018 at 12:39
• I have a C-compiler installed. I suppose that you don't use any Linux distribution as operating system; those usually come along with gcc included. On macos, make sure to install XCode and its command line tools. On Windows, you can install MinGW (which provides you also with gcc) or Microsoft Visual Studio. Jun 8, 2018 at 12:45
• You don't want a -\[Infinity] pop up in any numerical computations since this will slow down everything (you could not use compiled code). With a C-compiler installed, all error messages will disappear. In particular, there are no numerical errors related to any infinite expressions. Jun 8, 2018 at 12:47
• Yes, exactly. I just need an improvement of my written version. Your version is superb but it is not doable in my machine. May you revise it, in order to be useful in ordinary (Windows-based without C comilper) machines as well? Jun 8, 2018 at 12:51