# Speed up summation of many inverse matrices

Physical background

I am trying to write a code that computes the Green function of a system, given the Self energy $$\hat{\Sigma}_n$$, a list of Matsubara frequencies $$z_n$$ ($$n=1,...,N_{Mats}$$) and the Hamiltonian $$\hat{H}_k$$, where $$k$$ is a momentum index that runs from $$k=1,...,LE$$. The equation is the following: $$\hat{G}_n = \sum_{k=1}^{LE} W_k \left[ z_n \hat{\mathbb{I}} - \hat{H}_k - \hat{\Sigma}_n \right]^{-1},$$ where $$W_k$$ is just a list of weights (real positive numbers), and the matrices are obviously square matrices and have an arbitrary size $$N_{orb}$$.

Implementation

I have tried to write a compiled function that does the job: i.e. given the $$(LE, N_{orb}, N_{orb})$$-dimensional tensor $$\hat{H}$$, the $$(N_{Matsubara},N_{orb},N_{orb})$$-dimensional tensor $$\hat{\Sigma}$$, and two lists $$z$$ and $$W$$ of length $$N_{Matsubara}$$ and $$LE$$ respectively, I want to return a $$(N_{Matsubara}, N_{orb}, N_{orb})$$-dimensional tensor $$\hat{G}$$. Here's my code:

G = Compile[{
{H, _Real, 3}, {W, _Real, 1}, {Sigma, _Complex, 3}, {z, _Complex, 1}
},
With[
{LE = Length[H], Norb = Length[H[[1]]], NMatsubara = Length[z]},
Total@Table[
W[[i]] * Inverse[#] & /@ (
(IdentityMatrix[Norb]*#) & /@ z -
ConstantArray[H[[i]], NMatsubara] - Sigma
)
, {i, LE}]
],
RuntimeAttributes -> {Listable}, Parallelization -> True
];


My problem

I have tested this code with randomly generated input tensors chosing $$N_{Matsubara}=5000$$, $$LE=1000$$, $$N_{orb}=2$$ and the code takes 23 seconds to run on my machine, which is not terribly bad, but I wonder if there are some tricks to speed it up, since performance really matters for what I'm doing. For example, in real calculations $$N_{orb}$$ can be as large as $$5$$, and this function is just a small step of an iterative procedure... So my question is: can you suggest a more efficient way to write this function? Thanks in andvance for your help!

• Do you need the matrix itself or is it sufficient to have an operator that can perform the operation of G on any arbitrary vector? Because in that case, you might want to use LinearSolve instead of Inverse (which is really inefficient for large matrices). Commented Nov 10, 2022 at 12:00
• Actually I need the full matrix since I never apply $\hat{G}$ to vectors Commented Nov 10, 2022 at 12:38
• Hi. I think you may have provided a little too much info here. Maybe the background is relevant, but one cannot tell from your question. The amount of notation and jargon may put people off. The idea is to ask questions in a simple, appropriately abstract form that others might find useful. Commented Nov 11, 2022 at 6:47

Note: This answer is only for the $$2\times 2$$ case, and it uses the explicit formula for the matrix inverse. The same idea may also work for other small dimensions, which is the case that OP seems to be interested in.

I will consider a simplified problem: Given an $$n \times 2 \times 2$$ array of real or complex numbers, where $$n$$ is large, we want to add the inverses of the $$n$$ matrices of size $$2 \times 2$$. Here are two implementations:

sumInverses[A_]:=Total[Map[Inverse,A]];

sumInversesNew[A_]:=With[{
a=A[[;;,1,1]],
b=A[[;;,1,2]],
c=A[[;;,2,1]],
d=A[[;;,2,2]]},
ArrayFlatten[{{d,-b},{-c,a}}].(1/(a*d-b*c))];


Example:

SeedRandom[1];
A=RandomComplex[{-1-I,1+I},{5000,2,2}];

sumInverses[A]==sumInversesNew[A]
(* true *)

RepeatedTiming[sumInverses[A];]
(* {0.024691,Null} *)

RepeatedTiming[sumInversesNew[A];]
(* {0.00051788,Null} *)


So the second implementation is quite a bit faster.

Arrays should always satisfy DeveloperPackedArrayQ`.