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 *)
```