Fix modified left eigenvectors to take account of new eigenvalue formulation

[csskene]

EVPs in Dedalus are of the form $L_\textrm{min} PR^H x = - \lambda M\textrm{min} P_R^H . x$, and are now solved without the preconditioner for $y= PR^H x$. The left eigenvectors are correct (as the transposed eigenvalue system can be written $L\textrm{min}^H x^\dagger = - \lambda^* M\textrm{min}^H . x^\dagger$ and the preconditioner drops out), but the modified left eigenvectors need to be altered to take account of the right preconditioner in the mass matrix $M= - M\textrm{min} P_R^H$

This pull request modifies _build_modified_left_eigenvectors so that the preconditioner in $M$ is taken account of when forming the modified left eigenvectors.

[csskene]

Apologies, I have submitted this too soon. Let me test it a bit more, I think I may be missing something.

[csskene]

I've gone through it carefully this morning and believe the fix is needed. I just explained it in a confusing way which was not quite correct. The full story is that the matrices $L\textrm{min}$ and $M\textrm{min}$ now contain the right preconditioner. This affects the modified left eigenvectors as outlined in my original post which is the current behaviour of Dedalus, i.e. $L_\textrm{min} PR^H x = - \lambda M\textrm{min} P_R^H . x$ is how eigenvalue problems are currently solved. The inclusion of $PR$ in the definition of $M\textrm{min}$ means that it no longer cancels out in the definition of the left modified eigenvectors, hence the need for this fix.

I hope this makes more sense!

[kburns]

I'm trying to go through the details here, and I think before merging we should add some documentation that describes exactly what we mean by left and modified left eigenvectors. In particular, we probably should be un-preconditioning the left eigenvectors themselves.

As a start, the right eigenvectors we return satisfy

$\lambda M X + L X = 0$

Adding preconditioners:

$\lambda (P_L M P_R) P_R^+ X + (P_L L P_R) P_R^+ X = 0$

This is

$\lambda \tilde{M} \tilde{X} + \tilde{L} \tilde{X} = 0$ where $\tilde{M} = P_L M P_R$, $\tilde{L} = P_L L P_R$ and $\tilde{X} = P_R^+ X$.

This is the "right" problem we actually solve for $(\lambda, \tilde{X})$, from which we recover $X = P_R \tilde{X}$, which is what we set as solver.eigenvectors.

Applying the same to the left eigenvectors, following scipy's notation

$\lambda Y^ M + Y^ L = 0$

$\lambda Y^ P_L^+ (P_L M P_R) + Y^ P_L^+ (P_L L P_R) = 0$

$\lambda \tilde{Y}^ \tilde{M} + \tilde{Y}^ \tilde{L} = 0$ where $\tilde{Y}^ = Y^ P_L^+$.

This is the "left" problem scipy solves for $(\lambda, \tilde{Y})$. So we should probably set solver.left_eigenvectors to be $Y = P_L^* \tilde{Y}$.

[kburns]

For orthogonality and normalizations we have

$Y_i^ L X_j = - \lambda_i Y_i^ M X_j = - \lambda_j Y_i^* M X_j$

$(\lambda_i - \lambda_j) Y_i^* M X_j = 0$

so $Y_i^* M X_j = 0$ if $\lambda_i \neq \lambda_j$.

We can define orthonormal duals such that $Z_i^* Xj = \delta{ij}$ by taking

$$Z_i^ = \frac{Y_i^ M}{Y_i^ M X_i} = \frac{\tilde{Y}_i^ \tilde{M} P_R^+}{\tilde{Y}_i^* \tilde{M} \tilde{X}_i}$$

or

$$Z_i = \frac{P_R^{+ } \tilde{M}^ \tilde{Y}_i}{(\tilde{Y}_i^ \tilde{M} \tilde{X}_i)^}$$

In our special case where the preconditioners are semi-unitary, we just have $P_R^{+* } = P_R$. I think this then basically matches what is in the code, up to a minus sign that I'm not quite sure about. Does this seem right @csskene and @jsoishi?

[csskene]

This all seems correct to me and I'm happy to help with adding documentation. I guess essentially at the moment the adjoints are with respect to $\tilde{L}$ rather than $L$, so the left pre-conditioner is incorporated in the adjoint. Switching to adjoints with respect to $L$ may make it easier to use them so is probably a good change.

I think the minus sign is because the modified left-eigenvectors are defined with the mass matrix as $-M$, but it drops out in the normalisation so the left-eigenvectors are not changed.

[kburns]

Ok I think this new commit fixes things. It adds some docs explaining all the above (but without preconditioners, for simplicity), and chances the code to match, so the stored vectors are always in the un-preconditioned state.

[csskene]

Thanks! I've checked it with some of my codes (cartesian and spherical shell) and it works as expected. The new documentation looks good too.

[kburns]

Great, thanks for testing!