Changes:

1) Documentation updated for RiskParityPortfolio(), SuccessiveConvexOptimizer(), and design()
2) Initial tau changed
3) Sign notation in Dmat changed for consistency (it does not affect anything from the outside).
4) Function design() now takes arguments verbose and control_numerical_ill_conditioning (which controls numerical issues).
5) Now the linear constraints can involve dummy variables apart from the portfolio weights.
6) Many more unit tests added.

Author: dppalomar

.gitignore

@@ -5,7 +5,8 @@ build/
 riskparityportfolio.egg-info/
 var/
 .DS_Store
-.coverage
+.idea/
+.coverage*
 docs/source/tutorials/.ipynb_checkpoints/*
 .eggs/
 .ipynb_checkpoints/*

riskparityportfolio/rpp.py

@@ -18,14 +18,22 @@ def __init__(self):
 
 
 class RiskParityPortfolio:
-    """Designs risk parity portfolios by solving the following optimization problem
+    """Designs risk parity portfolios by solving the following optimization problem:
 
-    minimize R(w) - alpha * mu.T * w + lambda * w.T Sigma w
-    subject to Cw = c, Dw <= d
+      minimize     R(w) - alpha * mu.T w + lambda * w.T Sigma w
+      subject to   Cw = c, Dw <= d
 
-    where R is a risk concentration function, and alpha and lambda are trade-off
+    where R(w) is a risk concentration function, and alpha and lambda are trade-off
     parameters for the expected return and the variance, respectively.
 
+    The risk concentration R(w) is computed as the squared l2-norm of the risk concentration vector,
+    which by default is obtained from RiskContribOverVarianceMinusBudget():
+
+      R(w) = sum_i (MR_i/sum(MR_i) - budget_i)^2
+
+    where MR_i = w_i * (Sigma @ w)_i are the marginal risks (sum(MR_i) = Var(w)), and
+          budget_i are the risk budgets (by default 1/n).
+
     Parameters
     ----------
     covariance : array, shape=(n, n)
@@ -36,6 +44,23 @@ class RiskParityPortfolio:
         weights of the portfolio
     risk_concentration : class
         any valid child class of RiskConcentrationFunction
+
+    Examples:
+    # Set up:
+    >>> import numpy as np
+    >>> import riskparityportfolio as rpp
+    >>> 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:
+    >>> my_portfolio = rpp.RiskParityPortfolio(Sigma)
+    >>> my_portfolio.design(Cmat=Cmat, cvec=cvec)
+    >>> 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__(
@@ -172,13 +197,18 @@ def add_variance(self, lmd):
         self.lmd = lmd
         self.has_variance = True
 
-    def design(self, **kwargs):
+    def design(self, verbose=True, control_numerical_ill_conditioning=False, **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.
+            Dictionary of parameters to be passed to SuccessiveConvexOptimizer().
         """
         self.sca = SuccessiveConvexOptimizer(self, **kwargs)
-        self.sca.solve()
+        self.sca.solve(verbose=verbose, control_numerical_ill_conditioning=control_numerical_ill_conditioning)

riskparityportfolio/sca.py

@@ -101,9 +101,25 @@ def maxiter(self, value):
 
 class SuccessiveConvexOptimizer:
     """
-    Successive Convex Approximation optimizer tailored for the risk parity problem.
+    Successive Convex Approximation optimizer tailored for the risk parity problem including the linear constraints:
+       Cmat @ w  = cvec
+       Dmat @ w <= dvec,
+    where matrices Cmat and Dmat have n columns (n being the number of assets). Based on the paper:
+
+    Feng, Y., and Palomar, D. P. (2015). SCRIP: Successive convex optimization methods for risk parity portfolios design.
+    IEEE Trans. Signal Processing, 63(19), 5285–5300.
+
+    By default, the constraints are set to sum(w) = 1 and w >= 0, i.e.,
+       Cmat = np.ones((1, n))
+       cvec = np.array([1.0])
+       Dmat = -np.eye(n)
+       dvec = np.zeros(n).
+
+    Notes:
+        1) If equality constraints are not needed, set Cmat = np.empty((0, n)) and cvec = [].
+        2) If the matrices Cmat and Dmat have more than n columns, it is assumed that the additional columns
+        (same number for both matrices) correspond to dummy variables (which do not appear in the objective function).
     """
-
     def __init__(
         self,
         portfolio,
@@ -119,28 +135,27 @@ def __init__(
         dvec=None,
     ):
         self.portfolio = portfolio
-        self.tau = tau or 0.05 * np.sum(np.diag(self.portfolio.covariance)) / (
-            2 * self.portfolio.number_of_assets
-        )
+        self.tau = tau or 1e-4  # 0.05 * np.trace(self.portfolio.covariance) / (2 * self.portfolio.number_of_assets)
         sca_validator = SuccessiveConvexOptimizerValidator()
         self.gamma = sca_validator.gamma = gamma
         self.zeta = sca_validator.zeta = zeta
         self.funtol = sca_validator.funtol = funtol
         self.wtol = sca_validator.wtol = wtol
         self.maxiter = sca_validator.maxiter = maxiter
         self.Cmat = Cmat
-        self.Dmat = Dmat
+        self.Dmat = Dmat  # Dmat @ w <= dvec
         self.cvec = cvec
         self.dvec = dvec
-        self.CCmat = np.vstack((self.Cmat, self.Dmat)).T
-        self.bvec = np.concatenate((self.cvec, self.dvec))
+        self.number_of_vars = self.Cmat.shape[1]
+        self.number_of_dummy_vars = self.number_of_vars - self.portfolio.number_of_assets
+        self.dummy_vars = np.zeros(self.number_of_dummy_vars)
+        self.CCmat = np.vstack((self.Cmat, -self.Dmat)).T  # CCmat.T @ w >= bvec
+        self.bvec = np.concatenate((self.cvec, -self.dvec))
         self.meq = self.Cmat.shape[0]
         self._funk = self.get_objective_function_value()
         self.objective_function = [self._funk]
         self._tauI = self.tau * np.eye(self.portfolio.number_of_assets)
-        self.Amat = (
-            self.portfolio.risk_concentration.jacobian_risk_concentration_vector()
-        )
+        self.Amat = self.portfolio.risk_concentration.jacobian_risk_concentration_vector()
         self.gvec = self.portfolio.risk_concentration.risk_concentration_vector
 
     @property
@@ -151,7 +166,7 @@ def Cmat(self):
     def Cmat(self, value):
         if value is None:
             self._Cmat = np.atleast_2d(np.ones(self.portfolio.number_of_assets))
-        elif np.atleast_2d(value).shape[1] == self.portfolio.number_of_assets:
+        elif np.atleast_2d(value).shape[1] >= self.portfolio.number_of_assets:
             self._Cmat = np.atleast_2d(value)
         else:
             raise ValueError(
@@ -166,9 +181,9 @@ def Dmat(self):
     @Dmat.setter
     def Dmat(self, value):
         if value is None:
-            self._Dmat = np.eye(self.portfolio.number_of_assets)
-        elif np.atleast_2d(value).shape[1] == self.portfolio.number_of_assets:
-            self._Dmat = -np.atleast_2d(value)
+            self._Dmat = -np.eye(self.portfolio.number_of_assets)
+        elif np.atleast_2d(value).shape[1] == self.Cmat.shape[1]:
+            self._Dmat = np.atleast_2d(value)
         else:
             raise ValueError(
                 "Dmat shape {} doesnt agree with the number of"
@@ -200,7 +215,7 @@ def dvec(self, value):
         if value is None:
             self._dvec = np.zeros(self.portfolio.number_of_assets)
         elif len(value) == self.Dmat.shape[0]:
-            self._dvec = -np.atleast_1d(value)
+            self._dvec = np.atleast_1d(value)
         else:
             raise ValueError(
                 "dvec shape {} doesnt agree with Dmat shape"
@@ -215,19 +230,43 @@ def get_objective_function_value(self):
             obj += self.portfolio.lmd * self.portfolio.volatility ** 2
         return obj
 
-    def iterate(self):
+    def iterate(self, verbose=True, control_numerical_ill_conditioning=False):
         wk = self.portfolio.weights
         g = self.gvec(wk)
         A = np.ascontiguousarray(self.Amat(wk))
         At = np.transpose(A)
         Q = 2 * At @ A + self._tauI
-        q = 2 * np.matmul(At, g) - np.matmul(Q, wk)
+        q = 2 * np.matmul(At, g) - Q @ wk  # np.matmul() is necessary here since g is not a numpy array
         if self.portfolio.has_variance:
             Q += self.portfolio.lmd * self.portfolio.covariance
         if self.portfolio.has_mean_return:
             q -= self.portfolio.alpha * self.portfolio.mean
+        if self.number_of_dummy_vars > 0:
+            Q = np.vstack([np.hstack([Q, np.zeros((self.portfolio.number_of_assets, self.number_of_dummy_vars))]),
+                           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
+                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]
-        self.portfolio.weights = wk + self.gamma * (w_hat - wk)
+        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 = (
@@ -238,19 +277,29 @@ def iterate(self):
             np.abs(self._funk - fun_next)
             <= 0.5 * self.funtol * (np.abs(self._funk) + np.abs(fun_next))
         ).all()
-        if has_w_converged or has_fun_converged:
+        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))
+            ).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:
             return False
         self.gamma = self.gamma * (1 - self.zeta * self.gamma)
         self._funk = fun_next
         return True
 
-    def solve(self):
+    def solve(self, verbose=True, control_numerical_ill_conditioning=False):
         i = 0
-        with tqdm(total=self.maxiter) as pbar:
-            while self.iterate() and i < self.maxiter:
-                i += 1
-                pbar.update()
-
+        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):
+                break
+            i += 1
 
 def project_line_and_box(weights, lower_bound, upper_bound):
     def objective_function(variable, weights):