diff --git a/dedalus/core/solvers.py b/dedalus/core/solvers.py
index 3de09ee3..ddd27155 100644
--- a/dedalus/core/solvers.py
+++ b/dedalus/core/solvers.py
@@ -212,8 +212,8 @@ def solve_dense(self, subproblem, rebuild_matrices=False, left=False, normalize_
         if left:
             self.eigenvalues, pre_left_evecs, pre_right_evecs = eig_output
             self.right_eigenvectors = self.eigenvectors = sp.pre_right @ pre_right_evecs
-            self.left_eigenvectors = sp.pre_left.conj().toarray() @ pre_left_evecs
-            self.modified_left_eigenvectors = (sp.M_min @ sp.pre_right_pinv).conj().toarray() @ pre_left_evecs
+            self.left_eigenvectors = sp.pre_left.conj().T @ pre_left_evecs
+            self.modified_left_eigenvectors = (sp.M_min @ sp.pre_right_pinv).conj().T @ pre_left_evecs
             if normalize_left:
                 norms = np.diag(pre_left_evecs.T.conj() @ sp.M_min @ pre_right_evecs)
                 self.left_eigenvectors /= np.conj(norms)
@@ -270,8 +270,8 @@ def solve_sparse(self, subproblem, N, target, rebuild_matrices=False, left=False
             # Note: this definition of "left eigenvectors" is consistent with the documentation for scipy.linalg.eig
             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.conj().toarray() @ pre_left_evecs
-            self.modified_left_eigenvectors = (sp.M_min @ sp.pre_right_pinv).conj().toarray() @ pre_left_evecs
+            self.left_eigenvectors = sp.pre_left.conj().T @ pre_left_evecs
+            self.modified_left_eigenvectors = (sp.M_min @ sp.pre_right_pinv).conj().T @ pre_left_evecs
             # Check that eigenvalues match
             if not np.allclose(self.eigenvalues, np.conjugate(self.left_eigenvalues)):
                 if raise_on_mismatch:
diff --git a/dedalus/libraries/matsolvers.py b/dedalus/libraries/matsolvers.py
index 811b4542..f301d4f2 100644
--- a/dedalus/libraries/matsolvers.py
+++ b/dedalus/libraries/matsolvers.py
@@ -137,7 +137,7 @@ def __init__(self, matrix, solver=None):
         if self.trans == "T":
             matrix = matrix.T
         elif self.trans == "H":
-            matrix = matrix.conj().toarray()
+            matrix = matrix.conj().T
         self.LU = spla.splu(matrix.tocsc(),
                             permc_spec=self.permc_spec,
                             diag_pivot_thresh=self.diag_pivot_thresh,
diff --git a/dedalus/tools/array.py b/dedalus/tools/array.py
index bf28cd73..ab9caf88 100644
--- a/dedalus/tools/array.py
+++ b/dedalus/tools/array.py
@@ -433,7 +433,7 @@ def matvec(x):
         # 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.conj().toarray().dot(x))
+            return solver.solve_H(B.conj().T.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)