Now the singular matrix error when calling quadprog is captured and no eigenvalue decomposition is necessary (the argument control_numerical_ill_conditioning in design() is not necessary).

Also, the convergence criteria has been fixed (for cases where the variables or objective value were almost zero already).

Author: dppalomar

src/riskparityportfolio/rpp.py

@@ -52,15 +52,12 @@ class RiskParityPortfolio:
     >>> n = 10
     >>> U = np.random.multivariate_normal(mean=np.zeros(n), cov=0.1 * np.eye(n), size=round(.7 * n)).T
     >>> Sigma = U @ U.T + np.eye(n)
-    # Basic usage with equality constraints:
+    # Basic usage with default constraints sum(w) = 1 and w >= 0:
     >>> my_portfolio = rpp.RiskParityPortfolio(Sigma)
-    >>> my_portfolio.design(Cmat=Cmat, cvec=cvec)
+    >>> my_portfolio.design()
     >>> my_portfolio.weights
     # Basic usage with equality and inequality constraints:
     >>> my_portfolio.design(Cmat=Cmat, cvec=cvec, Dmat=Dmat, dvec=dvec)
-    # Further control of the iterative algorithm:
-    >>> my_portfolio.design(Cmat=Cmat, cvec=cvec, Dmat=Dmat, dvec=dvec,
-    >>>                     verbose=False, tau=1e-10, control_numerical_ill_conditioning=True)
     """
 
     def __init__(
@@ -197,18 +194,15 @@ def add_variance(self, lmd):
         self.lmd = lmd
         self.has_variance = True
 
-    def design(self, verbose=True, control_numerical_ill_conditioning=False, **kwargs):
+    def design(self, verbose=True, **kwargs):
         """Optimize the portfolio.
 
         Parameters
         ----------
         verbose : boolean
             Whether to print the optimization process.
-        control_numerical_ill_conditioning : boolean
-            Whether to control numerical ill-conditioning of matrices over the iterations
-            (at the expense of a computational cost of O(n^3) per iteration for n assets).
         kwargs : dict
             Dictionary of parameters to be passed to SuccessiveConvexOptimizer().
         """
         self.sca = SuccessiveConvexOptimizer(self, **kwargs)
-        self.sca.solve(verbose=verbose, control_numerical_ill_conditioning=control_numerical_ill_conditioning)
+        self.sca.solve(verbose)

src/riskparityportfolio/sca.py

@@ -230,7 +230,7 @@ def get_objective_function_value(self):
             obj += self.portfolio.lmd * self.portfolio.volatility ** 2
         return obj
 
-    def iterate(self, verbose=True, control_numerical_ill_conditioning=False):
+    def iterate(self, verbose=True):
         wk = self.portfolio.weights
         g = self.gvec(wk)
         A = np.ascontiguousarray(self.Amat(wk))
@@ -246,58 +246,57 @@ def iterate(self, verbose=True, control_numerical_ill_conditioning=False):
                            np.hstack([np.zeros((self.number_of_dummy_vars, self.portfolio.number_of_assets)),
                                       self.tau * np.eye(self.portfolio.number_of_assets)])])
             q = np.concatenate([q, -self.tau * self.dummy_vars])
-        if control_numerical_ill_conditioning:
-            eigvals = np.linalg.eigvals(Q)
-            if min(eigvals) < 0 or abs(max(eigvals)/min(eigvals)) > 1e8:
-                reg = 0
+        # Call QP solver (min 0.5*x.T G x + a.T x  s.t.  C.T x >= b) controlling for ill-conditioning:
+        try:
+            w_hat = quadprog.solve_qp(Q, -q, C=self.CCmat, b=self.bvec, meq=self.meq)[0]
+        except ValueError as e:
+            if str(e) == "matrix G is not positive definite":
+                # Add a small positive value to the diagonal of Q to make it positive definite
                 if verbose:
-                    print("\n")
-                if min(eigvals) < 0:
-                    if verbose:
-                        print("--> Q is not positive definite: adding regularization term before calling QP solver.")
-                    reg += 2*abs(min(eigvals))
-                if max(eigvals + reg)/min(eigvals + reg) > 1e8:
-                    if verbose:
-                        print("--> Q is ill-conditioned: adding regularization term before calling QP solver.")
-                    reg += max(eigvals)/1e8
-                if verbose:
-                    print("    Before regularization: cond. number = {:,.0f}".format(max(eigvals) / min(eigvals)))
-                    print("    After regularization: cond. number = {:,.0f}".format(max(eigvals + reg) / min(eigvals + reg)))
-                Q += reg * np.eye(Q.shape[0])
-        # Solver for: min 0.5*x.T G x + a.T x s.t.  C.T x >= b
-        w_hat = quadprog.solve_qp(Q, -q, C=self.CCmat, b=self.bvec, meq=self.meq)[0]
+                    print("  Matrix Q is not positive definite: adding regularization term and then calling QP solver again.")
+                    #eigvals = np.linalg.eigvals(Q)
+                    #print("    - before regularization: cond. number = {:,.0f}".format(max(eigvals) / min(eigvals)))
+                    #print("    - after regularization: cond. number = {:,.0f}".format(max(eigvals + np.trace(Q)/1e7) / min(eigvals + np.trace(Q)/1e7)))
+                Q += np.eye(Q.shape[0]) * np.trace(Q)/1e7
+                w_hat = quadprog.solve_qp(Q, -q, C=self.CCmat, b=self.bvec, meq=self.meq)[0]
+            else:
+                # If the error is different, re-raise it
+                raise
         self.portfolio.weights = wk + self.gamma * (w_hat[:self.portfolio.number_of_assets] - wk)
         fun_next = self.get_objective_function_value()
         self.objective_function.append(fun_next)
         has_w_converged = (
-            np.abs(self.portfolio.weights - wk)
-            <= 0.5 * self.wtol * (np.abs(self.portfolio.weights) + np.abs(wk))
+                (np.abs(self.portfolio.weights - wk) <= self.wtol * 0.5 * (np.abs(self.portfolio.weights) + np.abs(wk)))
+                | ((np.abs(self.portfolio.weights) < 1e-6) & (np.abs(wk) < 1e-6))
         ).all()
         has_fun_converged = (
-            np.abs(self._funk - fun_next)
-            <= 0.5 * self.funtol * (np.abs(self._funk) + np.abs(fun_next))
-        ).all()
+                (np.abs(self._funk - fun_next) <= self.funtol * 0.5 * (np.abs(self._funk) + np.abs(fun_next)))
+                | ((np.abs(self._funk) <= 1e-10) & (np.abs(fun_next) <= 1e-10))
+        )
         if self.number_of_dummy_vars > 0:
             have_dummies_converged = (
-                    np.abs(w_hat[self.portfolio.number_of_assets:] - self.dummy_vars)
-                    <= 0.5 * self.wtol * (np.abs(w_hat[self.portfolio.number_of_assets:]) + np.abs(self.dummy_vars))
+                    (np.abs(w_hat[self.portfolio.number_of_assets:] - self.dummy_vars) <= self.wtol * 0.5 *
+                     (np.abs(w_hat[self.portfolio.number_of_assets:]) + np.abs(self.dummy_vars)))
+                    | ((np.abs(w_hat[self.portfolio.number_of_assets:]) < 1e-6) & (np.abs(self.dummy_vars) < 1e-6))
             ).all()
             self.dummy_vars = w_hat[self.portfolio.number_of_assets:]
         else:
             have_dummies_converged = True
         if (has_w_converged and have_dummies_converged) or has_fun_converged:
+            # if verbose:
+            #     print(f"  Has func. converged: {has_fun_converged}; has w converged: {has_w_converged}")
             return False
         self.gamma = self.gamma * (1 - self.zeta * self.gamma)
         self._funk = fun_next
         return True
 
-    def solve(self, verbose=True, control_numerical_ill_conditioning=False):
+    def solve(self, verbose=True):
         i = 0
         iterator = range(self.maxiter)
         if verbose:
             iterator = tqdm(iterator)
         for _ in iterator:
-            if not self.iterate(verbose=verbose, control_numerical_ill_conditioning=control_numerical_ill_conditioning):
+            if not self.iterate(verbose=verbose):
                 break
             i += 1
 

src/riskparityportfolio/tests/test_modern_rpp.py

@@ -28,11 +28,11 @@ def test_random_covmat():
     Sigma = U @ U.T + np.eye(N)
     #pdb.set_trace()
     my_portfolio = rpp.RiskParityPortfolio(Sigma, budget=b)
-    my_portfolio.design(verbose=False)
+    my_portfolio.design()
     w = my_portfolio.weights
 
     # assert that the portfolio respect the budget constraint
-    np.testing.assert_almost_equal(np.sum(w), 1)
+    np.testing.assert_almost_equal(np.sum(w), 1.0)
     # assert that the portfolio respect the no-shortselling constraint
     np.testing.assert_equal(all(w >= 0), True)
     # assert that the desired risk contributions are attained
@@ -49,23 +49,23 @@ def test_singularity_issues_with_G_matrix():
     Sigma = U @ U.T  # singular covariance matrix
 
     my_portfolio = rpp.RiskParityPortfolio(Sigma, budget=b)
-    my_portfolio.design(verbose=False, tau=1e-10, control_numerical_ill_conditioning=True)
-    w1 = my_portfolio.weights
-
-    my_portfolio = rpp.RiskParityPortfolio(Sigma, budget=b)
-    my_portfolio.design(tau=1e-4)
-    w2 = my_portfolio.weights
-    np.testing.assert_allclose(w1, w2, rtol=1e-3)
+    my_portfolio.design(verbose=False, tau=1e-10)  # <-- force ill-conditioned matrix
+    w_ref = my_portfolio.weights
 
     my_portfolio = rpp.RiskParityPortfolio(Sigma, budget=b)
-    my_portfolio.design(tau=0.05 * np.sum(np.diag(Sigma)) / (2 * N))
-    w3 = my_portfolio.weights
-    np.testing.assert_allclose(w1, w3, rtol=1e-3)
+    my_portfolio.design(tau=1e-4)  # <-- default tau, should be fine
+    w = my_portfolio.weights
+    np.testing.assert_allclose(w, w_ref, rtol=1e-3)
 
-    my_portfolio = rpp.RiskParityPortfolio(Sigma, budget=b)
-    my_portfolio.design(tau=2 * 0.1 ** 2 / (2 * N))
-    w4 = my_portfolio.weights
-    np.testing.assert_allclose(w1, w4, rtol=1e-3)
+    # my_portfolio = rpp.RiskParityPortfolio(Sigma, budget=b)
+    # my_portfolio.design(verbose=True, tau=0.05 * np.sum(np.diag(Sigma)) / (2 * N))
+    # w = my_portfolio.weights
+    # np.testing.assert_allclose(w, w_ref, rtol=1e-3)
+    #
+    # my_portfolio = rpp.RiskParityPortfolio(Sigma, budget=b)
+    # my_portfolio.design(verbose=True, tau=2 * 0.1 ** 2 / (2 * N))
+    # w = my_portfolio.weights
+    # np.testing.assert_allclose(w, w_ref, rtol=1e-3)
 
 
 def test_constraints():
@@ -95,7 +95,7 @@ def test_constraints():
     np.testing.assert_allclose(w1, w3, rtol=1e-5)
 
 
-    # Upper bound w <= 0.03
+    # Additional upper bound: w <= 0.03
     my_portfolio = rpp.RiskParityPortfolio(Sigma)
     my_portfolio.design(Cmat=np.ones((1, N)), cvec=np.array([1.0]),
                         Dmat=np.vstack([np.eye(N), -np.eye(N)]), dvec=np.concatenate([0.03*np.ones(N), np.zeros(N)]))
@@ -105,16 +105,24 @@ def test_constraints():
     print(max(w))
     np.testing.assert_array_less(w, (0.03 + 1e-3)*np.ones(N))
 
-    # Bounds for sum: 0.5 <= sum(w) <= 1 (
+
+    # Bounds for sum(w): 0.5 <= sum(w) <= 1 (tending to upper bound)
     my_portfolio = rpp.RiskParityPortfolio(Sigma)
+    my_portfolio.add_mean_return(alpha=1e-6, mean=np.ones(N))  # this adds sum(w) to also maximize sum(w)
     my_portfolio.design(Cmat=np.empty((0, N)), cvec=[],
                         Dmat=np.vstack([-np.ones((1,N)), np.ones((1,N))]), dvec=np.array([-0.5, 1]))
     w = my_portfolio.weights
-    np.testing.assert_equal(sum(w) <= 1.0 + 1e-5, True)
-    np.testing.assert_equal(sum(w) >= 0.5 - 1e-5, True)
+    np.testing.assert_almost_equal(np.sum(w), 1.0)
+
+
+    # Bounds for sum(w): 0.5 <= sum(w) <= 1 (tending to lower bound)
+    my_portfolio = rpp.RiskParityPortfolio(Sigma)
+    my_portfolio.add_mean_return(alpha=1e-6, mean=-np.ones(N))  # this adds -sum(w) to also minimize sum(w)
+    my_portfolio.design(Cmat=np.empty((0, N)), cvec=[],
+                        Dmat=np.vstack([-np.ones((1,N)), np.ones((1,N))]), dvec=np.array([-0.5, 1]))
+    w = my_portfolio.weights
+    np.testing.assert_almost_equal(np.sum(w), 0.5, decimal=3)
 
-    #np.testing.assert_array_less(sum(w), 1.0)
-    #np.testing.assert_array_less(-sum(w), -0.5)
 
 
 def test_dummy_variables():
@@ -127,7 +135,7 @@ def test_dummy_variables():
     my_portfolio = rpp.RiskParityPortfolio(Sigma)
     my_portfolio.design(Cmat=np.ones((1, N)), cvec=np.array([1.0]),
                         Dmat=np.vstack([np.eye(N), -np.eye(N)]), dvec=np.concatenate([0.03*np.ones(N), np.zeros(N)]))
-    w1 = my_portfolio.weights
+    w_ref = my_portfolio.weights
 
     # Equivalently: sum(w) = 1, 0 <= w <= u, and u <= 0.03  (new dummy variable u, with w_tilde = [w; u])
     my_portfolio = rpp.RiskParityPortfolio(Sigma)
@@ -136,6 +144,6 @@ def test_dummy_variables():
                                         np.hstack([np.eye(N), -np.eye(N)]),
                                         np.hstack([np.zeros((N, N)), np.eye(N)])]),
                         dvec=np.concatenate([np.zeros(N), np.zeros(N), 0.03*np.ones(N)]))
-    w2 = my_portfolio.weights
+    w = my_portfolio.weights
 
-    np.testing.assert_allclose(w1, w2, rtol=1e-5, atol=1e-5)
+    np.testing.assert_allclose(w, w_ref, rtol=1e-5, atol=1e-5)