diff --git a/pypfopt/cla.py b/pypfopt/cla.py
index bde53b19..595472cf 100644
--- a/pypfopt/cla.py
+++ b/pypfopt/cla.py
@@ -192,6 +192,45 @@ def _compute_lambda(self, covarF_inv, covarFB, meanF, wB, i, bi):
         res = float(res[0, 0])
         return res, bi
 
+    @staticmethod
+    def _invert_covar_free(covarF):
+        """
+        Invert the covariance submatrix of the currently free assets.
+
+        The critical line algorithm requires this submatrix to be invertible. A
+        singular covariance matrix – e.g. one containing perfectly correlated or
+        duplicated assets, or estimated from fewer observations than assets – can
+        make it singular, in which case numpy raises a bare ``LinAlgError`` from
+        deep inside the algorithm. Translate it into an actionable error instead.
+
+        Parameters
+        ----------
+        covarF : np.ndarray
+            covariance submatrix of the free assets
+
+        Raises
+        ------
+        ValueError
+            if the submatrix is singular
+
+        Returns
+        -------
+        np.ndarray
+            inverse of ``covarF``
+        """
+        try:
+            return np.linalg.inv(covarF)
+        except np.linalg.LinAlgError:
+            raise ValueError(
+                "The covariance matrix is singular, so the critical line algorithm "
+                "cannot invert it. This usually means some assets are perfectly "
+                "correlated or duplicated, or that the covariance matrix was "
+                "estimated from fewer observations than there are assets. Consider "
+                "removing redundant assets, or using a shrinkage estimator such as "
+                "risk_models.CovarianceShrinkage, which returns a positive definite "
+                "matrix."
+            ) from None
+
     def _get_matrices(self, f):
         # Slice covarF,covarFB,covarB,meanF,meanB,wF,wB
         covarF = self._reduce_matrix(self.cov_matrix, f, f)
@@ -339,7 +378,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_covar_free(covarF)
                 j = 0
                 for i in f:
                     lam, bi = self._compute_lambda(
@@ -354,7 +393,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_covar_free(covarF)
                     lam, bi = self._compute_lambda(
                         covarF_inv,
                         covarFB,
@@ -371,7 +410,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_covar_free(covarF)
                 meanF = np.zeros(meanF.shape)
             else:
                 # 4) decide lambda
@@ -383,7 +422,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_covar_free(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..e9412676 100644
--- a/tests/test_cla.py
+++ b/tests/test_cla.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 
-from pypfopt import risk_models
+from pypfopt import expected_returns, risk_models
 from pypfopt.cla import CLA
 from tests.utilities_for_tests import get_data, setup_cla
 
@@ -96,6 +96,42 @@ def test_cla_two_assets():
     assert CLA(mu, cov)
 
 
+def test_cla_singular_covariance_raises():
+    # Perfectly correlated assets give a positive semidefinite but singular
+    # covariance matrix, which CLA cannot invert. The user should get an
+    # actionable error rather than a bare numpy LinAlgError.
+    mu = np.array([0.1, 0.2])
+    cov = np.array([[0.01, 0.01], [0.01, 0.01]])
+
+    for method in ("max_sharpe", "min_volatility", "efficient_frontier"):
+        with pytest.raises(ValueError, match="cannot invert") as e:
+            getattr(CLA(mu, cov), method)()
+        # np.linalg.LinAlgError subclasses ValueError, so assert explicitly that
+        # the raw numpy error is not what surfaces to the user.
+        assert not isinstance(e.value, np.linalg.LinAlgError)
+
+
+def test_cla_singular_covariance_fixed_by_shrinkage():
+    # The remedy suggested by the error message should actually work: shrinkage
+    # returns a positive definite matrix, which CLA can handle.
+    df = get_data()
+    # duplicate a column to make the sample covariance singular
+    df = df.iloc[:, :3].copy()
+    df[df.columns[1]] = df[df.columns[0]]
+    mu = expected_returns.mean_historical_return(df)
+
+    sample_cov = risk_models.sample_cov(df)
+    assert np.min(np.linalg.eigvalsh(sample_cov)) < 1e-12
+    with pytest.raises(ValueError, match="cannot invert"):
+        CLA(mu, sample_cov).min_volatility()
+
+    shrunk_cov = risk_models.CovarianceShrinkage(df).ledoit_wolf()
+    assert np.min(np.linalg.eigvalsh(shrunk_cov)) > 0
+    w = CLA(mu, shrunk_cov).min_volatility()
+    assert len(w) == 3
+    np.testing.assert_almost_equal(sum(w.values()), 1)
+
+
 def test_cla_max_sharpe_semicovariance():
     df = get_data()
     cla = setup_cla()