Possible speed improvement for numerical linear algebra?

Question

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/
       NDSolve`FEM`MapThreadDot[s, s$conj.Conjugate[hmatrix]]));

(* iterate kk times, ensure hermitian after numerical errors *)
find$h[hmat_, iterations_] := Block[{h$temp = 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 *)
```

Answer 1

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}, {s$conj, _Complex, 2}},
       s$trans.(tf*(1/
         NDSolve`FEM`MapThreadDot[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 NDSolve`FEM`MapThreadDot 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.