# Possible speed improvement for numerical linear algebra?

I have a question about further optimising a numerical linear algebra calculation. I've trawled through many useful posts on this site and have ended up with something quite efficient - I was wondering if there is any further speed to be gained?

The problem: I'm computing an operator that takes a constant hermitian matrix and outputs another hermitian matrix. The new matrix is determined by computing the value of some functions at a point, forming a kind of inner product on the functions, and then summing this contribution for each point with some weight.

More precisely, I calculate

$$T(h)_{\alpha\bar\beta}=\sum_{i=1}^{N_\text{p}}\frac{s_\alpha(p_i)\,\overline{s_\beta (p_i)}\,w_i}{\sum_{\gamma,\bar{\delta}=1}^{N_\text{f}}h^{\gamma\bar{\delta}}s_\gamma(p_i)\,\overline{s_\delta (p_i)}},\qquad \alpha,\bar\beta=1,\ldots,N_\text{f}$$

where $$N_\text{p}$$ is the number of points, $$N_{\text{f}}$$ is the number of functions, $$p_i$$ is a point, $$h^{\gamma\bar{\delta}}$$ is the hermitian matrix, $$s_\gamma(p_i)$$ and $$w_i$$ are the values of the functions and the weight for the point $$p_i$$. The new $$h$$ matrix is then given by taking an inverse transpose of this. One then iterates this procedure 10 times.

The number of points ranges from 200,000 to 1,000,000 and the number of functions from 5 to 2,000. I already have the value of each function at each point and the weights in a packed array.

Following what I've seen on this site, I've tried to 'vectorise' the calculation as much as possible. I'm also using what seems to be a faster version of MapThread and have swapped transposes of the hermitian matrix for conjugates.

I was wondering if there is any further speed up to be found? I feel like I'm computing the products of things that stay the same with every iteration when I don't need to. Is it worth compiling this or are Mathematica's inbuilt functions as fast as it gets? I have access to a 16 core machine, but I haven't found a way to parallelize the calculation.

Thanks in advance for any suggestions!

(* number of points and functions *)
Np = 200000;
Nf = 200;

(* value of each function at each point - Np x Nf array; save conjugate etc for speed *)
s = RandomComplex[{-0.1, 0.1}, {Np, Nf}];
s$$conj = Conjugate[s]; s$$trans = Transpose[s];

(* weight for each point - Np array *)
weights = RandomReal[{0, 0.001}, Np];

(* pre-compute suggested by eyorble *)
tf = s$conj*weights; (* operator to find new h matrix *) T[hmatrix_] := s$$trans.(tf*(1/ NDSolveFEMMapThreadDot[s, s$$conj.Conjugate[hmatrix]])); (* iterate kk times, ensure hermitian after numerical errors *) find$$h[hmat_, iterations_] := Block[{htemp = hmat, kk = iterations}, Do[h$$temp = 1/2 (# + ConjugateTranspose[#]) &[Transpose[Inverse[T[h$$temp]]]], kk]; h$$temp] (* define non-singular initial guess - Nf x Nf array *) h = DiagonalMatrix[ConstantArray[1.0, Nf]]; (* iterate 10 times *) h = find$h[h, 10]; // AbsoluteTiming

(* 16.1577s *)
$$$$

• Not really my area of expertise, but have you taken a look at numerical BLAS routines which are exposed in Mma since version 11.2 (and as an undocumented feature earlier)? See reference.wolfram.com/language/LowLevelLinearAlgebra/guide/… May 29, 2019 at 4:23
• Also of potential interest: mathematica.stackexchange.com/questions/149597/… May 29, 2019 at 4:48
• Thanks for the suggestions. I'm looking into BLAS at the moment though, from what I've seen, Mathematica's implementation of many of the functions is of the same performance (if not slightly faster) than BLAS itself. I'll keep reading! May 30, 2019 at 10:14

I found some modest improvements with only a marginal difference in the final result by using Compile and one more stored calculation.

Firstly, s$conj * weights is a constant: tf = s$conj * weights;


This is the source of the difference in the calculation, so if that's an issue this would be where to look. It's on the order of $$10^{-13}$$, however, so it's probably negligible.

Let Tc be the compiled version of T:

Tc = Compile[{{hmatrix, _Complex, 2}, {s, _Complex, 2},
{s$$trans, _Complex, 2}, {tf, _Complex, 2}, {sconj, _Complex, 2}}, s$$trans.(tf*(1/
NDSolveFEMMapThreadDot[s, s$conj.Conjugate[hmatrix]])), CompilationTarget -> "C", Parallelization -> True];  Note that it takes 5 arguments, rather than just 1. Note that more improvement comes from precalculating Tf than from compiling Tc, but it's still a slight improvement. I suspect NDSolveFEMMapThreadDot is simply passed through in the compilation, but it works and doesn't cause any issues. Compiling the iterating operation also helps a bit: htipci = Compile[{{x, _Complex, 2}}, 1/2 (Transpose[Inverse[x]] + Conjugate[Inverse[x]]), CompilationTarget -> "C", Parallelization -> True];  Note that this arrangement results in one less Transpose being performed, and I do hope that Compile manages to reuse the Inverse[x]. find$$h2[hmat_, iterations_] := Nest[htipci[Tc[#, s, s$$trans, tf, s$conj]] &, hmat, iterations];


I replaced Do with Nest. I will admit, I'm not entirely sure why, but it does result in a slight speedup. I would expect it's because Nest knows that it'll have to reallocate new memory for its next result, but I have no way of verifying that.

Compared to the original implementation of find$h, the speed is about 18% faster (15% less time) on my machine. h0 = DiagonalMatrix[ConstantArray[1.0, Nf]]; h = find$h[h, 10]; // AbsoluteTiming


{17.8535, Null}

h2 = find\$h2[h0, 10]; // AbsoluteTiming


{15.0781, Null}

Max[Abs[h2 - h]]


4.54747*10^-13

Some notes: light testing suggests that the Parallelization -> True option in both of the above compilations is only a very slight improvement over omitting it. However, it is still an improvement. Most of the parallelization being used is baked into the matrix operations to begin with. You should notice a significant speed increase on a machine with more cores if you haven't already tested it, even without explicitly calling the Parallel family of functions.

• Thanks for the suggestions! As you say, the greatest speed increase seems to come from pre-computing tf. The compilation gives an extra couple of percent increase. I was wondering if there is any way to take advantage of {Listable} in Compile, since I'm evaluating the T operator for each point? I can't quite think of how to do this though... May 30, 2019 at 10:13
• @ala10 I definitely tried to look for a way to convert that to Listable, but I can't say I see it. I don't think they'll help with the number of places Transpose and Inverse are used. I'm going to try and look at the BLAS subroutines that kirma mentioned too, but I'm not sure there's much room for improvement on this one. May 30, 2019 at 16:45
• I had another think about this and managed to find a 10% speed up by using a compiled version of MapThread[Dot...] that uses Listable to run over each point: MapThreadDot11 = Compile[{{m, _Complex, 1}, {x, _Complex, 1}}, Dot[m, x], RuntimeOptions -> "Speed", Parallelization -> True, RuntimeAttributes -> {Listable}, CompilationTarget -> "C"] and then replacing NDSolveFEMMapThreadDot with MapThreadDot11 everywhere. For the MapThread bit, Np=200000 and Nf=500`, the timings are 38s for MapThread[Dot...], 25s for NDSolve... and 17s for the above. Jun 23, 2019 at 11:36