[d2 and d3] normalize_left in solve_sparse() is misleading if raise_on_mismatch=False and warning is thrown
afraser3 opened this issue · 1 comments
There's some behavior in solve_sparse() with left=True that might be too much of an edge case to be worth worrying about, and might be intended behavior (because it stems from disabling an error that would otherwise be raised, as described below), but I wanted to flag it because maybe it's worth rethinking this behavior. At the very least it stumped the heck out of me when I first encountered some issues that ultimately stemmed from this behavior.
Quick summary of what I propose: I think the user should be forbidden from setting both raise_on_mismatch=False and normalize_left=True, since the latter might be rendered completely pointless (and worse, misleading) if the left and right eigenmodes don't match.
Detailed explanation:
When using solve_sparse() with left=True, we first calculate the right eigenvalues/eigenvectors as usual, and then calculate the left eigenvalues/eigenvectors by solving for the right eigenvalues/eigenvectors of the conjugate-transpose system (here is where that's done in d2, but the same is done in d3 -- this method for calculating the left eigenvectors is demonstrated, e.g., in equation 6.28 of Keaton's thesis). We then check if the right eigenvalues are the same as the conjugate of the left eigenvalues, and if they aren't, raise an error if raise_on_mismatch=True (default), or merely a warning if raise_on_mismatch=False. Next, if no error is raised, and if normalize_left=True (default), the left eigenvectors are normalized: they are each multiplied by a scalar to assure that the modified left eigenvectors (the left eigenvectors multiplied by the mass matrix) are orthonormal to the right eigenvectors, rather than merely orthogonal (see equation 6.30 in Keaton's thesis).
My issue: the scalar that we multiply the left eigenvectors by when normalizing them is only the correct scalar if the left and right eigenvalues actually match. Let's say we have a system where we've solved for 2 left eigenmodes and 2 right eigenmodes, and the eigenvalues are indeed compatible except that the indexing is swapped. So, rather than eigenvalues[0] = np.conj(left_eigenvalues[0]) and eigenvalues[1] = np.conj(left_eigenvalues[1]), we have eigenvalues[0] = np.conj(left_eigenvalues[1]) and eigenvalues[1] = np.conj(left_eigenvalues[0]). If raise_on_mismatch=False, then a warning will be raised rather than an error, and the 0th (say) left eigenvector will be normalized by dividing it by modified_left_eigenvector.T.conj()[0] @ eigenvector[0] when in fact it should be divided by modified_left_eigenvector.T.conj()[0] @ eigenvector[1].
Maybe this is an inherent risk that the user assumes when setting raise_on_mismatch=False. But I wonder if better behavior might be to forbid the user to set both raise_on_mismatch=False and normalize_left=True, since the latter might be rendered completely pointless (and worse, misleading) if the left and right eigenmodes don't match.
In the specific case I'm dealing with, I had 300 left and 300 right eigenmodes, and I was able to confirm that the two sets of modes were indeed compatible but just indexed completely differently, and I was able to sort them manually and then do the normalization myself. But that's hard to replicate generally for reasons I'm happy to go into.