diff --git a/kwant/physics/boundstate.py b/kwant/physics/boundstate.py
index e948cd94a716a1c64c5ad3eb41e16a6b5ca3ca8a..9419be19cceeca7538a3c5551296b8d10b2fdd1c 100644
--- a/kwant/physics/boundstate.py
+++ b/kwant/physics/boundstate.py
@@ -50,9 +50,12 @@ else:
             return self.solve(x.astype(self.dtype))
 
     def sparse_diag(matrix, k, sigma, **kwargs):
-        if sigma != 0:
-            matrix = matrix - sigma * sp.identity(matrix.shape[0])
-        return sp.linalg.eigsh(matrix, k, sigma=sigma, OPinv=_LuInv(matrix), **kwargs)
+        if sigma == 0:
+            shifted_matrix = matrix
+        else:
+            shifted_matrix = matrix - sigma * sp.identity(matrix.shape[0])
+        return sp.linalg.eigsh(matrix, k, sigma=sigma, OPinv=_LuInv(shifted_matrix),
+                               **kwargs)
 
 
 def find_boundstates(