I need to speed up this summation

Question

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}];

Thanks in advance.

Answer 1

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