diff --git a/dedalus/core/solvers.py b/dedalus/core/solvers.py
index c8f5dcbd..135b1933 100644
--- a/dedalus/core/solvers.py
+++ b/dedalus/core/solvers.py
@@ -250,15 +250,13 @@ def solve_sparse(self, subproblem, N, target, rebuild_matrices=False, left=False
         # Solve as sparse general eigenvalue problem
         A = sp.L_min
         B = - sp.M_min
-        # Solve for the right eigenvectors
-        self.eigenvalues, pre_right_evecs = scipy_sparse_eigs(A=A, B=B, N=N, target=target, matsolver=self.matsolver, **kw)
-        self.right_eigenvectors = self.eigenvectors = sp.pre_right @ pre_right_evecs
+        # Solve for the right (and optionally left) eigenvectors
+        eig_output = scipy_sparse_eigs(A=A, B=B, left=left, N=N, target=target, matsolver=self.matsolver, **kw)
+
         if left:
-            # Solve for the left eigenvectors
             # Note: this definition of "left eigenvectors" is consistent with the documentation for scipy.linalg.eig
-            self.left_eigenvalues, pre_left_evecs = scipy_sparse_eigs(A=A.getH(), B=B.getH(),
-                                                                      N=N, target=np.conjugate(target),
-                                                                      matsolver=self.matsolver, **kw)
+            self.eigenvalues, pre_right_evecs, self.left_eigenvalues, pre_left_evecs = eig_output
+            self.right_eigenvectors = self.eigenvectors = sp.pre_right @ pre_right_evecs
             self.left_eigenvectors = sp.pre_left.H @ pre_left_evecs
             self.modified_left_eigenvectors = (sp.M_min @ sp.pre_right_pinv).H @ pre_left_evecs
             # Check that eigenvalues match
@@ -278,6 +276,9 @@ def solve_sparse(self, subproblem, N, target, rebuild_matrices=False, left=False
                 norms = np.diag(pre_left_evecs.T.conj() @ sp.M_min @ pre_right_evecs)
                 self.left_eigenvectors /= np.conj(norms)
                 self.modified_left_eigenvectors /= np.conj(norms)
+        else:
+            self.eigenvalues, pre_right_evecs = eig_output
+            self.right_eigenvectors = self.eigenvectors = sp.pre_right @ pre_right_evecs
 
     def set_state(self, index, subsystem):
         """
diff --git a/dedalus/libraries/matsolvers.py b/dedalus/libraries/matsolvers.py
index 12e26f35..082e7e34 100644
--- a/dedalus/libraries/matsolvers.py
+++ b/dedalus/libraries/matsolvers.py
@@ -24,6 +24,9 @@ def __init__(self, matrix, solver=None):
     def solve(self, vector):
         pass
 
+    def solve_H(self,vector):
+        raise NotImplementedError("%s has not implemented 'solve_H' method" % type(self))
+
 
 @add_solver
 class DummySolver(SolverBase):
@@ -145,6 +148,14 @@ def __init__(self, matrix, solver=None):
     def solve(self, vector):
         return self.LU.solve(vector, trans=self.trans)
 
+    def solve_H(self,vector):
+        if self.trans == "N":
+            return self.LU.solve(vector, trans="H")
+        elif self.trans == "H":
+            return self.LU.solve(vector)
+        elif self.trans == "T":
+            return np.conj(self.LU.solve(np.conj(vector)))
+
 
 @add_solver
 class SuperluNaturalFactorized(_SuperluFactorizedBase):
diff --git a/dedalus/tools/array.py b/dedalus/tools/array.py
index e03a0fbd..0a196b7e 100644
--- a/dedalus/tools/array.py
+++ b/dedalus/tools/array.py
@@ -395,7 +395,7 @@ def drop_empty_rows(mat):
     return mat[nonempty_rows]
 
 
-def scipy_sparse_eigs(A, B, N, target, matsolver, **kw):
+def scipy_sparse_eigs(A, B, left, N, target, matsolver, **kw):
     """
     Perform targeted eigenmode search using the scipy/ARPACK sparse solver
     for the reformulated generalized eigenvalue problem
@@ -410,6 +410,8 @@ def scipy_sparse_eigs(A, B, N, target, matsolver, **kw):
         Sparse matrices for generalized eigenvalue problem
     N : int
         Number of eigenmodes to return
+    left: boolean
+        Whether to solve for the left eigenvectors or not
     target : complex
         Target σ for eigenvalue search
     matsolver : matrix solver class
@@ -427,7 +429,19 @@ def matvec(x):
     evals, evecs = spla.eigs(D, k=N, which='LM', sigma=None, **kw)
     # Rectify eigenvalues
     evals = 1 / evals + target
-    return evals, evecs
+    if left:
+        # Build sparse linear operator representing (A^H - conj(σ)B^H)^I B^H = C^-H B^H = left_D
+        # Note: left_D is not equal to D^H
+        def matvec_left(x):
+            return solver.solve_H(B.H.dot(x))
+        left_D = spla.LinearOperator(dtype=A.dtype, shape=A.shape, matvec=matvec_left)
+        # Solve using scipy sparse algorithm
+        left_evals, left_evecs = spla.eigs(left_D, k=N, which='LM', sigma=None, **kw)
+        # Rectify left eigenvalues
+        left_evals = 1 / left_evals + np.conj(target)
+        return evals, evecs, left_evals, left_evecs
+    else:
+        return evals, evecs
 
 
 def interleave_matrices(matrices):