1
$\begingroup$

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.

$\endgroup$
2
  • 2
    $\begingroup$ 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, $\endgroup$ Dec 3 '21 at 11:45
  • 1
    $\begingroup$ 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. $\endgroup$ Dec 3 '21 at 12:00
3
$\begingroup$

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.