diff --git a/pypfopt/cla.py b/pypfopt/cla.py
index bde53b19..01b0c2a3 100644
--- a/pypfopt/cla.py
+++ b/pypfopt/cla.py
@@ -327,6 +327,15 @@ def _eval_sr(self, a, w0, w1):
         c = np.dot(np.dot(w.T, self.cov_matrix), w)[0, 0] ** 0.5
         return b / c
 
+    @staticmethod
+    def _invert_covariance_matrix(cov_matrix):
+        try:
+            return np.linalg.inv(cov_matrix)
+        except np.linalg.LinAlgError as exc:
+            raise ValueError(
+                "CLA requires a positive definite covariance matrix"
+            ) from exc
+
     def _solve(self):
         # Compute the turning points,free sets and weights
         f, w = self._init_algo()
@@ -339,7 +348,7 @@ def _solve(self):
             l_in = None
             if len(f) > 1:
                 covarF, covarFB, meanF, wB = self._get_matrices(f)
-                covarF_inv = np.linalg.inv(covarF)
+                covarF_inv = self._invert_covariance_matrix(covarF)
                 j = 0
                 for i in f:
                     lam, bi = self._compute_lambda(
@@ -354,7 +363,7 @@ def _solve(self):
                 b = self._get_b(f)
                 for i in b:
                     covarF, covarFB, meanF, wB = self._get_matrices(f + [i])
-                    covarF_inv = np.linalg.inv(covarF)
+                    covarF_inv = self._invert_covariance_matrix(covarF)
                     lam, bi = self._compute_lambda(
                         covarF_inv,
                         covarFB,
@@ -371,7 +380,7 @@ def _solve(self):
                 # 3) compute minimum variance solution
                 self.ls.append(0)
                 covarF, covarFB, meanF, wB = self._get_matrices(f)
-                covarF_inv = np.linalg.inv(covarF)
+                covarF_inv = self._invert_covariance_matrix(covarF)
                 meanF = np.zeros(meanF.shape)
             else:
                 # 4) decide lambda
@@ -383,7 +392,7 @@ def _solve(self):
                     self.ls.append(l_out)
                     f.append(i_out)
                 covarF, covarFB, meanF, wB = self._get_matrices(f)
-                covarF_inv = np.linalg.inv(covarF)
+                covarF_inv = self._invert_covariance_matrix(covarF)
             # 5) compute solution vector
             wF, g = self._compute_w(covarF_inv, covarFB, meanF, wB)
             for i in range(len(f)):
diff --git a/tests/test_cla.py b/tests/test_cla.py
index c5f69e23..2fd7300e 100644
--- a/tests/test_cla.py
+++ b/tests/test_cla.py
@@ -52,6 +52,16 @@ def test_cla_max_sharpe_short():
     assert sharpe > long_only_sharpe
 
 
+def test_cla_max_sharpe_rejects_singular_covariance():
+    mu = np.array([0.1, 0.2])
+    cov = np.array([[0.01, 0.01], [0.01, 0.01]])
+
+    cla = CLA(mu, cov)
+
+    with pytest.raises(ValueError, match="positive definite covariance matrix"):
+        cla.max_sharpe()
+
+
 def test_cla_custom_bounds():
     bounds = [(0.01, 0.13), (0.02, 0.11)] * 10
     cla = setup_cla(weight_bounds=bounds)