diff --git a/pypfopt/cla.py b/pypfopt/cla.py
index bde53b19..95ff995f 100644
--- a/pypfopt/cla.py
+++ b/pypfopt/cla.py
@@ -192,6 +192,16 @@ def _compute_lambda(self, covarF_inv, covarFB, meanF, wB, i, bi):
         res = float(res[0, 0])
         return res, bi
 
+    @staticmethod
+    def _invert(covarF):
+        try:
+            return np.linalg.inv(covarF)
+        except np.linalg.LinAlgError as e:
+            raise ValueError(
+                "CLA requires an invertible covariance matrix. "
+                "The supplied covariance matrix appears to be singular."
+            ) from e
+
     def _get_matrices(self, f):
         # Slice covarF,covarFB,covarB,meanF,meanB,wF,wB
         covarF = self._reduce_matrix(self.cov_matrix, f, f)
@@ -339,7 +349,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(covarF)
                 j = 0
                 for i in f:
                     lam, bi = self._compute_lambda(
@@ -354,7 +364,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(covarF)
                     lam, bi = self._compute_lambda(
                         covarF_inv,
                         covarFB,
@@ -371,7 +381,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(covarF)
                 meanF = np.zeros(meanF.shape)
             else:
                 # 4) decide lambda
@@ -383,7 +393,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(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..f319fba6 100644
--- a/tests/test_cla.py
+++ b/tests/test_cla.py
@@ -96,6 +96,16 @@ def test_cla_two_assets():
     assert CLA(mu, cov)
 
 
+def test_cla_singular_covariance_raises_clear_error():
+    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="invertible covariance matrix"):
+        cla.max_sharpe()
+
+
 def test_cla_max_sharpe_semicovariance():
     df = get_data()
     cla = setup_cla()