Repository: skfolio/skfolio

Issue:
Problem with max_weights in HERC
First, thank you for this library. I do appreciate that you have created skfolio and made it available for any interested investor. 

I like the idea of HERC, which resonates to my own experiences on the markets. I aim for reasonably high cumulative return while valuing a smooth return curve for the portfolio. 

However, using HERC as such, it seems quite often to lead to portfolio compositions which may be computationally correct, but very unbalanced.  As a remedy, I have tried to use weight constraints to hinder too large deviation from naïve (1/n) weighting in portfolio composition.

Here I have run into problems. Setting max_weights does not seem to have effect I thought it would. It does not seem to cap the weight of any single asset in the portfolio, as I’d like to see.

Using your example of Hierarchical Equal Risk Contribution – CdaR as a reference,
there are 16 assets in the portfolio.

If I have code like this


```
# The relevant parts only

weight_leeway = 0.05
max_weights = float(1/len(assets) + weight_leeway) # 0.1125
transaction_costs = float(0.004/(21*4))

assets = [
'BRK-B',
 'CME',
 'INVE-B.ST',
 'V',
 'NDA-FI.HE',
 'MUV2.DE',
 'MSFT',
 'NOVO-B.CO',
 'VVSM.DE',
 'TW',
 'FICO',
 'ALV.DE',
 'WRT1V.HE',
 'SAP.DE',
 'GOOGL',
 'BWXT']

# download the data from yf, startdate=2020-01-01 end_date=2024-10-21 and remove the blanks (NaN)

X = prices_to_returns(prices)
X_train, X_test = train_test_split(X, test_size=0.33, shuffle=False)

model3 = HierarchicalEqualRiskContribution(
    risk_measure=RiskMeasure.CDAR,
    distance_estimator=KendallDistance(absolute=True),
    max_weights=max_weights,
    transaction_costs=transaction_costs,
    portfolio_params=dict(name="HERC-CDaR-Ward-Kendal")
)
model3.fit(X_train)
model3.weights_
```

```
array([0.02376651, 0.10343303, 0.0746895 , 0.08083832, 0.07263537,
       0.06226895, 0.0201518 , 0.08157259, 0.032706  , 0.03112282,
       0.21489253, 0.01559069, 0.04922161, 0.10995555, 0.01663003,
       0.01052472])
```
```
import skfolio; print(skfolio.__version__)
-->
0.4.2
```

EqualRisk portfolio returns weights to be 0.0625. The Portfolio HERC-CDaR-Ward-Kendal gives 0.21489258 (here for Novo Nordisk, which has had quite some fluctuations this year), (HERC-CDaR-Ward-Pearson gives 0.308940 with weight_constraints similarly set) 

I would like to have it and any other asset capped at max 0.1125 in this test , but seem to get almost double the weight of a single asset into portfolio.

I have tried to look at the herc.py code to understand if I am doing it wrong or whether there is something simple to correct in the code, but my python skills are not at adequate level to find the culprit. 

I can see from the code that there is intention to use the weight constraints after the HERC weights have been calculated, but I cannot understand the exact working of this. (No, I am not a programmer) 

Of course, I could post-process the portfolio weights in a separate code set and to calculate my own returns for the portfolio, but this would mean to abandoning the skfolio stats (portfolio.summary), which I find otherwise appealing.

Have I misunderstood the max_weights use or is use of weight constraints still under development or is there something wrong in the code?


Best regards,

Olli

Hints from issue comments:
@HugoDelatte: Hello, thank you!

You are correct that HERC can be more unbalanced and concentrated compared to HRP.

As explained in this example: https://skfolio.org/auto_examples/5_clustering/plot_2_herc_cdar.html#linkage-methods, 
since HERC relies on the dendrogram structure (unlike HRP), the clustering can be significantly affected by the choice of linkage method and distance estimator.

The approach would be to first focus on the clustering (choosing the right linkage and distance), which will result in more balanced weights. As a last resort, you can indeed apply weight constraints.

Thanks for pointing out the bug. I'll investigate and push a fix this week.@HugoDelatte: @myllo issue is solved in v0.4.3.
You can find the details in the note section in:
https://skfolio.org/generated/skfolio.optimization.HierarchicalEqualRiskContribution.html

Thanks for spotting it.

Solution patch:
diff --git a/docs/api.rst b/docs/api.rst
index acdccd8d..fd84e48d 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -617,4 +617,5 @@ Functions
     stats.compute_optimal_n_clusters
     stats.rand_weights
     stats.rand_weights_dirichlet
+    stats.minimize_relative_weight_deviation
 
diff --git a/src/skfolio/optimization/cluster/hierarchical/_base.py b/src/skfolio/optimization/cluster/hierarchical/_base.py
index cfe5cc78..a8d40eca 100644
--- a/src/skfolio/optimization/cluster/hierarchical/_base.py
+++ b/src/skfolio/optimization/cluster/hierarchical/_base.py
@@ -52,8 +52,6 @@ class BaseHierarchicalOptimization(BaseOptimization, ABC):
             * ENTROPIC_RISK_MEASURE
             * FOURTH_CENTRAL_MOMENT
             * FOURTH_LOWER_PARTIAL_MOMENT
-            * SKEW
-            * KURTOSIS
 
         The default is `RiskMeasure.VARIANCE`.
 
@@ -80,12 +78,12 @@ class BaseHierarchicalOptimization(BaseOptimization, ABC):
 
     min_weights : float | dict[str, float] | array-like of shape (n_assets, ), default=0.0
         Minimum assets weights (weights lower bounds). Negative weights are not allowed.
-        If a float is provided, it is applied to each asset. `None` is equivalent to
-        `-np.Inf` (no lower bound). If a dictionary is provided, its (key/value) pair
-        must be the (asset name/asset minium weight) and the input `X` of the `fit`
-        methods must be a DataFrame with the assets names in columns. When using a
-        dictionary, assets values that are not provided are assigned a minimum weight
-        of `0.0`. The default is 0.0 (no short selling).
+        If a float is provided, it is applied to each asset.
+        If a dictionary is provided, its (key/value) pair must be the
+        (asset name/asset minium weight) and the input `X` of the `fit` methods must be
+        a DataFrame with the assets names in columns.
+        When using a dictionary, assets values that are not provided are assigned a
+        minimum weight of `0.0`. The default is 0.0 (no short selling).
 
         Example:
 
@@ -96,12 +94,12 @@ class BaseHierarchicalOptimization(BaseOptimization, ABC):
 
     max_weights : float | dict[str, float] | array-like of shape (n_assets, ), default=1.0
         Maximum assets weights (weights upper bounds). Weights above 1.0 are not
-        allowed. If a float is provided, it is applied to each asset. `None` is
-        equivalent to `+np.Inf` (no upper bound). If a dictionary is provided, its
-        (key/value) pair must be the (asset name/asset maximum weight) and the input `X`
-        of the `fit` method must be a DataFrame with the assets names in columns. When
-        using a dictionary, assets values that are not provided are assigned a minimum
-        weight of `1.0`. The default is 1.0 (each asset is below 100%).
+        allowed. If a float is provided, it is applied to each asset.
+        If a dictionary is provided, its (key/value) pair must be the
+        (asset name/asset maximum weight) and the input `X` of the `fit` method must be
+        a DataFrame with the assets names in columns.
+        When using a dictionary, assets values that are not provided are assigned a
+        minimum weight of `1.0`. The default is 1.0 (each asset is below 100%).
 
         Example:
 
@@ -388,57 +386,6 @@ def _convert_weights_bounds(self, n_assets: int) -> tuple[np.ndarray, np.ndarray
 
         return min_weights, max_weights
 
-    @staticmethod
-    def _apply_weight_constraints_to_alpha(
-        alpha: float,
-        max_weights: np.ndarray,
-        min_weights: np.ndarray,
-        weights: np.ndarray,
-        left_cluster: np.ndarray,
-        right_cluster: np.ndarray,
-    ) -> float:
-        """Apply weight constraints to the alpha multiplication factor of the
-        Hierarchical Tree Clustering algorithm.
-
-        Parameters
-        ----------
-        alpha : float
-            The alpha multiplication factor of the Hierarchical Tree Clustering
-            algorithm.
-
-         min_weights : ndarray of shape (n_assets,)
-            The weight lower bound 1D array.
-
-        max_weights : ndarray of shape (n_assets,)
-            The weight upper bound 1D array.
-
-        weights : np.ndarray of shape (n_assets,)
-            The assets weights.
-
-        left_cluster : ndarray of shape (n_left_cluster,)
-            Indices of the left cluster weights.
-
-        right_cluster : ndarray of shape (n_right_cluster,)
-            Indices of the right cluster weights.
-
-        Returns
-        -------
-        value : float
-            The transformed alpha incorporating the weight constraints.
-        """
-        alpha = min(
-            np.sum(max_weights[left_cluster]) / weights[left_cluster[0]],
-            max(np.sum(min_weights[left_cluster]) / weights[left_cluster[0]], alpha),
-        )
-        alpha = 1 - min(
-            np.sum(max_weights[right_cluster]) / weights[right_cluster[0]],
-            max(
-                np.sum(min_weights[right_cluster]) / weights[right_cluster[0]],
-                1 - alpha,
-            ),
-        )
-        return alpha
-
     def get_metadata_routing(self):
         # noinspection PyTypeChecker
         router = (
diff --git a/src/skfolio/optimization/cluster/hierarchical/_herc.py b/src/skfolio/optimization/cluster/hierarchical/_herc.py
index a25a01e2..39964461 100644
--- a/src/skfolio/optimization/cluster/hierarchical/_herc.py
+++ b/src/skfolio/optimization/cluster/hierarchical/_herc.py
@@ -3,8 +3,7 @@
 # Copyright (c) 2023
 # Author: Hugo Delatte <delatte.hugo@gmail.com>
 # License: BSD 3 clause
-# The risk measure generalization and constraint features are derived
-# from Riskfolio-Lib, Copyright (c) 2020-2023, Dany Cajas, Licensed under BSD 3 clause.
+# Weight constraints is a novel implementation, see docstring for more details.
 
 import numpy as np
 import numpy.typing as npt
@@ -20,6 +19,7 @@
     BaseHierarchicalOptimization,
 )
 from skfolio.prior import BasePrior, EmpiricalPrior
+from skfolio.utils.stats import minimize_relative_weight_deviation
 from skfolio.utils.tools import check_estimator
 
 
@@ -45,6 +45,32 @@ class HierarchicalEqualRiskContribution(BaseHierarchicalOptimization):
         which is more stable and has better properties than the single-linkage
         method [4]_.
 
+        Also, the initial paper does not provide an algorithm for handling weight
+        constraints, and no standard solution currently exists.
+        In contrast to HRP (Hierarchical Risk Parity), where weight constraints
+        can be applied to the split factor at each bisection step, HERC
+        (Hierarchical Equal Risk Contribution) cannot incorporate weight constraints
+        during the intermediate steps of the allocation. Therefore, in HERC, the
+        weight constraints must be enforced after the top-down allocation has been
+        completed.
+        In skfolio, we minimize the relative deviation of the final weights from
+        the initial weights. This is formulated as a convex optimization problem:
+
+        .. math::
+            \begin{cases}
+            \begin{aligned}
+            &\min_{w} & & \Vert \frac{w - w_{init}}{w_{init}} \Vert_{2}^{2} \\
+            &\text{s.t.} & & \sum_{i=1}^{N} w_{i} = 1 \\
+            & & & w_{min} \leq w_i \leq w_{max}, \quad \forall i
+            \end{aligned}
+            \end{cases}
+
+        The reason for minimizing the relative deviation (as opposed to the absolute
+        deviation) is that we want to limit the impact on the risk contribution of
+        each asset. Since HERC allocates inversely to risk, adjusting the weights
+        based on relative deviation ensures that the assets' risk contributions
+        remain proportionally consistent with the initial allocation.
+
     Parameters
     ----------
     risk_measure : RiskMeasure or ExtraRiskMeasure, default=RiskMeasure.VARIANCE
@@ -70,8 +96,6 @@ class HierarchicalEqualRiskContribution(BaseHierarchicalOptimization):
             * ENTROPIC_RISK_MEASURE
             * FOURTH_CENTRAL_MOMENT
             * FOURTH_LOWER_PARTIAL_MOMENT
-            * SKEW
-            * KURTOSIS
 
         The default is `RiskMeasure.VARIANCE`.
 
@@ -98,12 +122,12 @@ class HierarchicalEqualRiskContribution(BaseHierarchicalOptimization):
 
     min_weights : float | dict[str, float] | array-like of shape (n_assets, ), default=0.0
         Minimum assets weights (weights lower bounds). Negative weights are not allowed.
-        If a float is provided, it is applied to each asset. `None` is equivalent to
-        `-np.Inf` (no lower bound). If a dictionary is provided, its (key/value) pair
-        must be the (asset name/asset minium weight) and the input `X` of the `fit`
-        methods must be a DataFrame with the assets names in columns. When using a
-        dictionary, assets values that are not provided are assigned a minimum weight
-        of `0.0`. The default is 0.0 (no short selling).
+        If a float is provided, it is applied to each asset.
+        If a dictionary is provided, its (key/value) pair must be the
+        (asset name/asset minium weight) and the input `X` of the `fit` methods must be
+        a DataFrame with the assets names in columns.
+        When using a dictionary, assets values that are not provided are assigned a
+        minimum weight of `0.0`. The default is 0.0 (no short selling).
 
         Example:
 
@@ -114,12 +138,12 @@ class HierarchicalEqualRiskContribution(BaseHierarchicalOptimization):
 
     max_weights : float | dict[str, float] | array-like of shape (n_assets, ), default=1.0
         Maximum assets weights (weights upper bounds). Weights above 1.0 are not
-        allowed. If a float is provided, it is applied to each asset. `None` is
-        equivalent to `+np.Inf` (no upper bound). If a dictionary is provided, its
-        (key/value) pair must be the (asset name/asset maximum weight) and the input `X`
-        of the `fit` method must be a DataFrame with the assets names in columns. When
-        using a dictionary, assets values that are not provided are assigned a minimum
-        weight of `1.0`. The default is 1.0 (each asset is below 100%).
+        allowed. If a float is provided, it is applied to each asset.
+        If a dictionary is provided, its (key/value) pair must be the
+        (asset name/asset maximum weight) and the input `X` of the `fit` method must be
+        a DataFrame with the assets names in columns.
+        When using a dictionary, assets values that are not provided are assigned a
+        minimum weight of `1.0`. The default is 1.0 (each asset is below 100%).
 
         Example:
 
@@ -208,6 +232,19 @@ class HierarchicalEqualRiskContribution(BaseHierarchicalOptimization):
         `management_fees`, `previous_weights` and `risk_free_rate` are copied from the
         optimization model and passed to the portfolio.
 
+    solver : str, default="CLARABEL"
+        The solver used for the weights constraints optimization. The default is
+        "CLARABEL" which is written in Rust and has better numerical stability and
+        performance than ECOS and SCS.
+        For more details about available solvers, check the CVXPY documentation:
+        https://www.cvxpy.org/tutorial/advanced/index.html#choosing-a-solver
+
+    solver_params : dict, optional
+        Solver parameters. For example, `solver_params=dict(verbose=True)`.
+        The default (`None`) is to use the CVXPY default.
+        For more details about solver arguments, check the CVXPY documentation:
+        https://www.cvxpy.org/tutorial/advanced/index.html#setting-solver-options
+
     Attributes
     ----------
     weights_ : ndarray of shape (n_assets,)
@@ -251,6 +288,8 @@ def __init__(
         hierarchical_clustering_estimator: HierarchicalClustering | None = None,
         min_weights: skt.MultiInput | None = 0.0,
         max_weights: skt.MultiInput | None = 1.0,
+        solver: str = "CLARABEL",
+        solver_params: dict | None = None,
         transaction_costs: skt.MultiInput = 0.0,
         management_fees: skt.MultiInput = 0.0,
         previous_weights: skt.MultiInput | None = None,
@@ -268,6 +307,8 @@ def __init__(
             previous_weights=previous_weights,
             portfolio_params=portfolio_params,
         )
+        self.solver = solver
+        self.solver_params = solver_params
 
     def fit(
         self, X: npt.ArrayLike, y: None = None, **fit_params
@@ -301,6 +342,13 @@ def fit(
             raise TypeError(
                 "`risk_measure` must be of type `RiskMeasure` or `ExtraRiskMeasure`"
             )
+
+        if self.risk_measure in [ExtraRiskMeasure.SKEW, ExtraRiskMeasure.KURTOSIS]:
+            # Because Skew and Kurtosis can take negative values
+            raise ValueError(
+                f"risk_measure {self.risk_measure} currently not supported" f"in HERC"
+            )
+
         self.prior_estimator_ = check_estimator(
             self.prior_estimator,
             default=EmpiricalPrior(),
@@ -393,21 +441,12 @@ def _recurse(node):
 
             left_cluster = np.array(left_cluster)
             right_cluster = np.array(right_cluster)
+
             left_risk = np.sum(cluster_risks[left_cluster])
             right_risk = np.sum(cluster_risks[right_cluster])
 
             alpha = 1 - left_risk / (left_risk + right_risk)
 
-            # Weights constraints
-            alpha = self._apply_weight_constraints_to_alpha(
-                alpha=alpha,
-                weights=weights,
-                max_weights=max_weights,
-                min_weights=min_weights,
-                left_cluster=left_cluster,
-                right_cluster=right_cluster,
-            )
-
             clusters_weights[left_cluster] *= alpha
             clusters_weights[right_cluster] *= 1 - alpha
 
@@ -421,5 +460,15 @@ def _recurse(node):
         for i, cluster_ids in enumerate(clusters):
             weights[cluster_ids] *= clusters_weights[i]
 
+        # Apply weights constraints
+        weights = minimize_relative_weight_deviation(
+            weights=weights,
+            min_weights=min_weights,
+            max_weights=max_weights,
+            solver=self.solver,
+            solver_params=self.solver_params,
+        )
+
         self.weights_ = weights
+
         return self
diff --git a/src/skfolio/optimization/cluster/hierarchical/_hrp.py b/src/skfolio/optimization/cluster/hierarchical/_hrp.py
index 0f20b66b..04309eb4 100644
--- a/src/skfolio/optimization/cluster/hierarchical/_hrp.py
+++ b/src/skfolio/optimization/cluster/hierarchical/_hrp.py
@@ -72,8 +72,6 @@ class HierarchicalRiskParity(BaseHierarchicalOptimization):
             * ENTROPIC_RISK_MEASURE
             * FOURTH_CENTRAL_MOMENT
             * FOURTH_LOWER_PARTIAL_MOMENT
-            * SKEW
-            * KURTOSIS
 
         The default is `RiskMeasure.VARIANCE`.
 
@@ -100,9 +98,9 @@ class HierarchicalRiskParity(BaseHierarchicalOptimization):
 
     min_weights : float | dict[str, float] | array-like of shape (n_assets, ), default=0.0
         Minimum assets weights (weights lower bounds). Negative weights are not allowed.
-        If a float is provided, it is applied to each asset. `None` is equivalent to
-        `-np.Inf` (no lower bound). If a dictionary is provided, its (key/value) pair
-        must be the (asset name/asset minium weight) and the input `X` of the `fit`
+        If a float is provided, it is applied to each asset.
+        If a dictionary is provided, its (key/value) pair must be the
+        (asset name/asset minium weight) and the input `X` of the `fit`
         methods must be a DataFrame with the assets names in columns. When using a
         dictionary, assets values that are not provided are assigned a minimum weight
         of `0.0`. The default is 0.0 (no short selling).
@@ -116,12 +114,12 @@ class HierarchicalRiskParity(BaseHierarchicalOptimization):
 
     max_weights : float | dict[str, float] | array-like of shape (n_assets, ), default=1.0
         Maximum assets weights (weights upper bounds). Weights above 1.0 are not
-        allowed. If a float is provided, it is applied to each asset. `None` is
-        equivalent to `+np.Inf` (no upper bound). If a dictionary is provided, its
-        (key/value) pair must be the (asset name/asset maximum weight) and the input `X`
-        of the `fit` method must be a DataFrame with the assets names in columns. When
-        using a dictionary, assets values that are not provided are assigned a minimum
-        weight of `1.0`. The default is 1.0 (each asset is below 100%).
+        allowed. If a float is provided, it is applied to each asset.
+        If a dictionary is provided, its (key/value) pair must be the
+        (asset name/asset maximum weight) and the input `X` of the `fit` method must
+        be a DataFrame with the assets names in columns.
+        When using a dictionary, assets values that are not provided are assigned a
+        minimum weight of `1.0`. The default is 1.0 (each asset is below 100%).
 
         Example:
 
@@ -296,6 +294,13 @@ def fit(
             raise TypeError(
                 "`risk_measure` must be of type `RiskMeasure` or `ExtraRiskMeasure`"
             )
+
+        if self.risk_measure in [ExtraRiskMeasure.SKEW, ExtraRiskMeasure.KURTOSIS]:
+            # Because Skew and Kurtosis can take negative values
+            raise ValueError(
+                f"risk_measure {self.risk_measure} currently not supported" f"in HRP"
+            )
+
         self.prior_estimator_ = check_estimator(
             self.prior_estimator,
             default=EmpiricalPrior(),
@@ -365,7 +370,7 @@ def fit(
                 left_cluster, right_cluster = clusters_ids
                 alpha = 1 - left_risk / (left_risk + right_risk)
                 # Weights constraints
-                alpha = self._apply_weight_constraints_to_alpha(
+                alpha = _apply_weight_constraints_to_split_factor(
                     alpha=alpha,
                     weights=weights,
                     max_weights=max_weights,
@@ -379,3 +384,54 @@ def fit(
 
         self.weights_ = weights
         return self
+
+
+def _apply_weight_constraints_to_split_factor(
+    alpha: float,
+    max_weights: np.ndarray,
+    min_weights: np.ndarray,
+    weights: np.ndarray,
+    left_cluster: np.ndarray,
+    right_cluster: np.ndarray,
+) -> float:
+    """
+    Apply weight constraints to the split factor alpha of the ,Hierarchical Tree
+    Clustering algorithm.
+
+    Parameters
+    ----------
+    alpha : float
+        The split factor alpha of the Hierarchical Tree Clustering algorithm.
+
+    min_weights : ndarray of shape (n_assets,)
+        The weight lower bound 1D array.
+
+    max_weights : ndarray of shape (n_assets,)
+        The weight upper bound 1D array.
+
+    weights : np.ndarray of shape (n_assets,)
+        The assets weights.
+
+    left_cluster : ndarray of shape (n_left_cluster,)
+        Indices of the left cluster weights.
+
+    right_cluster : ndarray of shape (n_right_cluster,)
+        Indices of the right cluster weights.
+
+    Returns
+    -------
+    value : float
+        The transformed split factor alpha incorporating the weight constraints.
+    """
+    alpha = min(
+        np.sum(max_weights[left_cluster]) / weights[left_cluster[0]],
+        max(np.sum(min_weights[left_cluster]) / weights[left_cluster[0]], alpha),
+    )
+    alpha = 1 - min(
+        np.sum(max_weights[right_cluster]) / weights[right_cluster[0]],
+        max(
+            np.sum(min_weights[right_cluster]) / weights[right_cluster[0]],
+            1 - alpha,
+        ),
+    )
+    return alpha
diff --git a/src/skfolio/utils/stats.py b/src/skfolio/utils/stats.py
index 95e22a17..5102513f 100644
--- a/src/skfolio/utils/stats.py
+++ b/src/skfolio/utils/stats.py
@@ -10,9 +10,11 @@
 # Statsmodels, Copyright (C) 2006, Jonathan E. Taylor, Licensed under BSD 3 clause.
 from enum import auto
 
+import cvxpy as cp
 import numpy as np
 import scipy.cluster.hierarchy as sch
 import scipy.optimize as sco
+import scipy.sparse.linalg as scl
 import scipy.spatial.distance as scd
 import scipy.special as scs
 from scipy.sparse import csr_matrix
@@ -34,6 +36,7 @@
     "compute_optimal_n_clusters",
     "rand_weights",
     "rand_weights_dirichlet",
+    "minimize_relative_weight_deviation",
 ]
 
 
@@ -488,3 +491,87 @@ def compute_optimal_n_clusters(distance: np.ndarray, linkage_matrix: np.ndarray)
     # k=0 represents one cluster
     k = np.argmax(gaps) + 2
     return k
+
+
+def minimize_relative_weight_deviation(
+    weights: np.ndarray,
+    min_weights: np.ndarray,
+    max_weights: np.ndarray,
+    solver: str = "CLARABEL",
+    solver_params: dict | None = None,
+) -> np.ndarray:
+    r"""
+    Apply weight constraints to an initial array of weights by minimizing the relative
+    weight deviation of the final weights from the initial weights.
+
+    .. math::
+            \begin{cases}
+            \begin{aligned}
+            &\min_{w} & & \Vert \frac{w - w_{init}}{w_{init}} \Vert_{2}^{2} \\
+            &\text{s.t.} & & \sum_{i=1}^{N} w_{i} = 1 \\
+            & & & w_{min} \leq w_i \leq w_{max}, \quad \forall i
+            \end{aligned}
+            \end{cases}
+
+    Parameters
+    ----------
+    weights : ndarray of shape (n_assets,)
+        Initial weights.
+
+    min_weights : ndarray of shape (n_assets,)
+        Minimum assets weights (weights lower bounds).
+
+    max_weights : ndarray of shape (n_assets,)
+        Maximum assets weights (weights upper bounds).
+
+    solver : str, default="CLARABEL"
+        The solver to use. The default is "CLARABEL" which is written in Rust and has
+        better numerical stability and performance than ECOS and SCS.
+        For more details about available solvers, check the CVXPY documentation:
+        https://www.cvxpy.org/tutorial/advanced/index.html#choosing-a-solver
+
+    solver_params : dict, optional
+        Solver parameters. For example, `solver_params=dict(verbose=True)`.
+        The default (`None`) is to use the CVXPY default.
+        For more details about solver arguments, check the CVXPY documentation:
+        https://www.cvxpy.org/tutorial/advanced/index.html#setting-solver-options
+    """
+    if not (weights.shape == min_weights.shape == max_weights.shape):
+        raise ValueError("`min_weights` and `max_weights` must have same size")
+
+    if np.any(weights < 0):
+        raise ValueError("Initial weights must be strictly positive")
+
+    if not np.isclose(np.sum(weights), 1.0):
+        raise ValueError("Initial weights must sum to one")
+
+    if np.any(max_weights < min_weights):
+        raise ValueError("`min_weights` must be lower or equal to `max_weights`")
+
+    if np.all((weights >= min_weights) & (weights <= max_weights)):
+        return weights
+
+    if solver_params is None:
+        solver_params = {}
+
+    n = len(weights)
+    w = cp.Variable(n)
+
+    objective = cp.Minimize(cp.norm(w / weights - 1))
+    constraints = [cp.sum(w) == 1, w >= min_weights, w <= max_weights]
+    problem = cp.Problem(objective, constraints)
+
+    try:
+        problem.solve(solver=solver, **solver_params)
+
+        if w.value is None:
+            raise cp.SolverError("No solution found")
+
+    except (cp.SolverError, scl.ArpackNoConvergence):
+        raise cp.SolverError(
+            f"Solver '{solver}' failed. Try another"
+            " solver, or solve with solver_params=dict(verbose=True) for more"
+            " information"
+        ) from None
+
+    return w.value


Test patch:
diff --git a/tests/test_optimization/test_cluster/test_hierarchical/test_herc.py b/tests/test_optimization/test_cluster/test_hierarchical/test_herc.py
index 996471a2..71beea39 100644
--- a/tests/test_optimization/test_cluster/test_hierarchical/test_herc.py
+++ b/tests/test_optimization/test_cluster/test_hierarchical/test_herc.py
@@ -11,7 +11,11 @@
 
 @pytest.fixture(
     scope="module",
-    params=list(RiskMeasure) + list(ExtraRiskMeasure),
+    params=[
+        x
+        for x in [*RiskMeasure, *ExtraRiskMeasure]
+        if x not in [ExtraRiskMeasure.SKEW, ExtraRiskMeasure.KURTOSIS]
+    ],
 )
 def risk_measure(request):
     return request.param
@@ -149,3 +153,76 @@ def test_metadata_routing(X_medium, implied_vol_medium):
 
     # noinspection PyUnresolvedReferences
     assert model.prior_estimator_.covariance_estimator_.r2_scores_.shape == (20,)
+
+
+def test_herc_weight_constraints(X):
+    model = HierarchicalEqualRiskContribution(
+        risk_measure=RiskMeasure.STANDARD_DEVIATION,
+    )
+    model.fit(X)
+    np.testing.assert_almost_equal(model.weights_[0], 0.02584786935368332)
+    np.testing.assert_almost_equal(model.weights_[-1], 0.0421175871045949)
+    np.testing.assert_almost_equal(sum(model.weights_), 1.0)
+
+    # Min Weights
+    model.set_params(min_weights={"AAPL": 0.05, "XOM": 0.08})
+    model.fit(X)
+    np.testing.assert_almost_equal(model.weights_[0], 0.05)
+    np.testing.assert_almost_equal(model.weights_[-1], 0.08)
+    np.testing.assert_almost_equal(sum(model.weights_), 1.0)
+
+    model.set_params(min_weights=0.05)
+    model.fit(X)
+    np.testing.assert_almost_equal(model.weights_, np.ones(20) * 0.05)
+    np.testing.assert_almost_equal(sum(model.weights_), 1.0)
+
+    # Max Weights
+    model.set_params(min_weights=0)
+    model.set_params(max_weights={"AAPL": 0.01, "XOM": 0.03})
+    model.fit(X)
+    np.testing.assert_almost_equal(model.weights_[0], 0.01)
+    np.testing.assert_almost_equal(model.weights_[-1], 0.03)
+    np.testing.assert_almost_equal(sum(model.weights_), 1.0)
+
+    model.set_params(max_weights=0.05)
+    model.fit(X)
+    np.testing.assert_almost_equal(model.weights_, np.ones(20) * 0.05)
+    np.testing.assert_almost_equal(sum(model.weights_), 1.0)
+
+    # Both
+    model.set_params(min_weights={"AAPL": 0.05}, max_weights={"XOM": 0.03})
+    model.fit(X)
+    np.testing.assert_almost_equal(model.weights_[0], 0.05)
+    np.testing.assert_almost_equal(model.weights_[-1], 0.03)
+    np.testing.assert_almost_equal(sum(model.weights_), 1.0)
+
+    model.set_params(min_weights=0.03, max_weights=0.06)
+    model.fit(X)
+    np.testing.assert_almost_equal(
+        model.weights_,
+        np.array(
+            [
+                0.03218923,
+                0.05585784,
+                0.06,
+                0.03,
+                0.05336249,
+                0.06,
+                0.04001747,
+                0.05895233,
+                0.06,
+                0.06,
+                0.03674032,
+                0.04708938,
+                0.03484055,
+                0.06,
+                0.04486381,
+                0.06,
+                0.05002204,
+                0.03796618,
+                0.05914401,
+                0.05895433,
+            ]
+        ),
+    )
+    np.testing.assert_almost_equal(sum(model.weights_), 1.0)
diff --git a/tests/test_optimization/test_cluster/test_hierarchical/test_hrp.py b/tests/test_optimization/test_cluster/test_hierarchical/test_hrp.py
index e4f53496..59b10efc 100644
--- a/tests/test_optimization/test_cluster/test_hierarchical/test_hrp.py
+++ b/tests/test_optimization/test_cluster/test_hierarchical/test_hrp.py
@@ -16,7 +16,11 @@ def small_X(X):
 
 @pytest.fixture(
     scope="module",
-    params=list(RiskMeasure) + list(ExtraRiskMeasure),
+    params=[
+        x
+        for x in [*RiskMeasure, *ExtraRiskMeasure]
+        if x not in [ExtraRiskMeasure.SKEW, ExtraRiskMeasure.KURTOSIS]
+    ],
 )
 def risk_measure(request):
     return request.param
@@ -173,3 +177,76 @@ def test_metadata_routing(X_medium, implied_vol_medium):
 
     # noinspection PyUnresolvedReferences
     assert model.prior_estimator_.covariance_estimator_.r2_scores_.shape == (20,)
+
+
+def test_hrp_weight_constraints(X):
+    model = HierarchicalRiskParity(
+        risk_measure=RiskMeasure.STANDARD_DEVIATION,
+    )
+    model.fit(X)
+    np.testing.assert_almost_equal(model.weights_[0], 0.030624328591088226)
+    np.testing.assert_almost_equal(model.weights_[-1], 0.05811358822991056)
+    np.testing.assert_almost_equal(sum(model.weights_), 1.0)
+
+    # Min Weights
+    model.set_params(min_weights={"AAPL": 0.05, "XOM": 0.08})
+    model.fit(X)
+    np.testing.assert_almost_equal(model.weights_[0], 0.05)
+    np.testing.assert_almost_equal(model.weights_[-1], 0.08)
+    np.testing.assert_almost_equal(sum(model.weights_), 1.0)
+
+    model.set_params(min_weights=0.05)
+    model.fit(X)
+    np.testing.assert_almost_equal(model.weights_, np.ones(20) * 0.05)
+    np.testing.assert_almost_equal(sum(model.weights_), 1.0)
+
+    # Max Weights
+    model.set_params(min_weights=0)
+    model.set_params(max_weights={"AAPL": 0.01, "XOM": 0.03})
+    model.fit(X)
+    np.testing.assert_almost_equal(model.weights_[0], 0.01)
+    np.testing.assert_almost_equal(model.weights_[-1], 0.03)
+    np.testing.assert_almost_equal(sum(model.weights_), 1.0)
+
+    model.set_params(max_weights=0.05)
+    model.fit(X)
+    np.testing.assert_almost_equal(model.weights_, np.ones(20) * 0.05)
+    np.testing.assert_almost_equal(sum(model.weights_), 1.0)
+
+    # Both
+    model.set_params(min_weights={"AAPL": 0.05}, max_weights={"XOM": 0.03})
+    model.fit(X)
+    np.testing.assert_almost_equal(model.weights_[0], 0.05)
+    np.testing.assert_almost_equal(model.weights_[-1], 0.03)
+    np.testing.assert_almost_equal(sum(model.weights_), 1.0)
+
+    model.set_params(min_weights=0.03, max_weights=0.06)
+    model.fit(X)
+    np.testing.assert_almost_equal(
+        model.weights_,
+        np.array(
+            [
+                0.03,
+                0.03,
+                0.03,
+                0.04712633,
+                0.05191067,
+                0.06,
+                0.06,
+                0.06,
+                0.03,
+                0.06,
+                0.06,
+                0.06,
+                0.04068081,
+                0.06,
+                0.05825588,
+                0.06,
+                0.03115588,
+                0.05087043,
+                0.06,
+                0.06,
+            ]
+        ),
+    )
+    np.testing.assert_almost_equal(sum(model.weights_), 1.0)
diff --git a/tests/test_utils/test_stats.py b/tests/test_utils/test_stats.py
index f2a6f7ba..4d5253f3 100644
--- a/tests/test_utils/test_stats.py
+++ b/tests/test_utils/test_stats.py
@@ -1,3 +1,4 @@
+import cvxpy as cp
 import numpy as np
 import pytest
 import scipy.cluster.hierarchy as sch
@@ -17,6 +18,7 @@
     cov_nearest,
     cov_to_corr,
     is_cholesky_dec,
+    minimize_relative_weight_deviation,
     n_bins_freedman,
     n_bins_knuth,
     rand_weights,
@@ -69,6 +71,21 @@ def linkage_matrix(distance):
     return linkage_matrix
 
 
+@pytest.fixture()
+def weights():
+    weights = np.array(
+        [
+            0.15686274156862742,
+            0.09803922098039221,
+            0.07843137078431371,
+            0.12500000125000002,
+            0.20833333208333335,
+            0.3333333333333333,
+        ]
+    )
+    return weights
+
+
 def test_n_bins_freedman(returns):
     n_bins = n_bins_freedman(returns)
     assert n_bins == 329
@@ -478,3 +495,28 @@ def test_raise_value_error_if_input_covariance_matrix_not_symmetric(self):
         cov = np.array([[1, 2], [3, 4]])
         with pytest.raises(ValueError):
             cov_nearest(cov)
+
+
+class TestMinimizeRelativeWeightDeviation:
+    def tests_no_constraints(self, weights):
+        new_weights = minimize_relative_weight_deviation(
+            weights=weights, min_weights=np.zeros(6), max_weights=np.ones(6)
+        )
+
+        np.testing.assert_array_almost_equal(weights, new_weights)
+
+    def tests_constraints(self, weights):
+        new_weights = minimize_relative_weight_deviation(
+            weights=weights, min_weights=np.ones(6) * 0.1, max_weights=np.ones(6) * 0.18
+        )
+
+        assert np.isclose(sum(new_weights), 1.0)
+        np.testing.assert_array_almost_equal(
+            new_weights, np.array([0.18, 0.16116994, 0.11883005, 0.18, 0.18, 0.18])
+        )
+
+    def test_non_feasible(self, weights):
+        with pytest.raises(cp.SolverError):
+            _ = minimize_relative_weight_deviation(
+                weights=weights, min_weights=np.zeros(6), max_weights=np.ones(6) * 0.1
+            )