# I need to speed up this summation

I tried solving this problem on my own but I cant find a proper solution.

I need to evaluate the following equation: $$\sum\limits_{m=1\\m\neq i}^n{\left(\dfrac{v_m~(\lambda_i\cdot B-A)~u_i}{(\lambda_m-\lambda_i)\cdot x_m}\right)\cdot u_m-\left(\dfrac{v_i~B~u_i}{2\cdot x_i}\right)\cdot u_i~,~i=1,...,n}$$

$$A$$ and $$B$$: real non-symmetric matricies

$$v_i$$: $$i$$th left-Eigenvector

$$u_i$$: $$i$$th right-Eigenvector

$$\lambda_i$$: $$i$$th Eigenvalue

The computation time is kinda fast for $$n\leq250$$ and it gets significantly slower for greater $$n$$. I need to evaluate this equation for at least $$n=500$$ and this takes forever. I posted my current "solution" below and would like to ask if there is a chance to speed up this computation.

n = 250;
A = RandomReal[{1*10^-7, 1*10^-6}, {n, n}];
B = RandomReal[{1*10^-12, 1*10^-11}, {n, n}];
X = RandomReal[{1*10^-12, 1*10^-11}, {n}];

{EW, EVR} = Eigensystem[{A, B}];
{EW, EVL} = Eigensystem[{Transpose[A], Transpose[B]}];

ParallelTable[
Sum[If[m == i, 0,
(EVL[[m]].(EW[[i]]*B - A).EVR[[i]])/((EW[[m]] - EW[[i]])*X[[m]])*EVR[[m]]], {m, 1, n}]
- EVR[[i]]*(EVL[[i]].B.EVR[[i]])/(2*X[[i]]), {i, 1, n}];


• No MMA on this machine to check, but from eyeballing your code the first thing I'd test is not checking If[m == I at every iteration; just iterate over the indices where m != i. Generic lesson - don't bury if statements at the bottom of deep loop nests if they can be lifted out, Dec 3 '21 at 11:45
• If $u_i$ is a right eigenvector to the eigenvalue $\lambda_i$, i.e., $(\lambda_i \, B - A) \, u_i = 0$, shouldn't the sum simply vanish? Only the summand for $m = i$ should remain. Dec 3 '21 at 12:00

Suppose that

$$A \, u_i = \lambda_i \, B \, u_i$$

and

$$v_m^T \, A = \lambda_m \, v_m^T \, B.$$

Then $$v_m^T( \lambda_i \,B - A) \, u_i = (\lambda_i \, v_m^T \, B - \lambda_m \, v_m^T \, B) \, u_i = (\lambda_i - \lambda_m) \, v_m^T \, B \, u_i,$$ hence $$\frac{v_m^T( \lambda_i \,B - A) \, u_i}{(\lambda_m-\lambda_i) \, x_m} \, u_m = - u_m \frac{v_m^T \, B \, u_i }{x_m}.$$

Also notice that if $$\lambda_i \neq 0$$ and $$\lambda_j \neq 0$$, then $$\lambda_m \, v_m^T \, B \, u_i = v_m^T \, A \, u_i = \lambda_i \, v_m^T \, B \, u_i,$$ hence $$v_m^T \, A \, u_i = v_m^T \, B \, u_i = 0 \quad \text{for i \neq m}.$$ This leads us to $$\frac{v_m^T( \lambda_i \,B - A) \, u_i}{(\lambda_m-\lambda_i) \, x_m} \, u_m = - u_m \frac{v_m^T \, B \, u_i }{x_m} = - \delta_{mi} \frac{v_m^T \, B \, u_m }{x_m} \, u_m$$

Btw. Be careful with the ordering of the eigenvalues:

RandomSeed[123];
n = 250;
A = RandomReal[{1*10^-7, 1*10^-6}, {n, n}];
B = RandomReal[{1*10^-12, 1*10^-11}, {n, n}];
X = RandomReal[{1*10^-12, 1*10^-11}, {n}];
{EWR, EVR} = Eigensystem[{A, B}];
{EWL, EVL} = Eigensystem[{Transpose[A], Transpose[B]}];

Max[Abs[EWR - EWL]]
perm = Flatten[Nearest[EWL -> "Index", EWR]];
EWR = EWR[[perm]];
EVR = EVR[[perm]];


1.99437

As you can see, due to small numerical errors, the orderings of left and right eigensystems do not coincide. This is easy to repair:

perm = Flatten[Nearest[EWL -> "Index", EWR]];
EWR = EWR[[perm]];
EVR = EVR[[perm]];
Max[Abs[EWR - EWL]/Abs[EWL]]


1.42593*10^-13

• Thank you so much. I completely missed that. This helps alot. Dec 3 '21 at 12:17
• Great! You're welcome! Dec 3 '21 at 12:23