ParDen: Surrogate Assisted Hyper-Parameter Optimisation for Portfolio Selection

1st T. L van Zyl1,
2nd M Woolway2,
3rd A Paskaramoorthy3
1Institute for Intelligent Systems

University of Johannesburg

Johannesburg, South Africa

Email<EMAIL_ADDRESS>2Faculty of Engineering and the Built Environment

University of Johannesburg

Johannesburg, South Africa

Email<EMAIL_ADDRESS>3School of Computer Science and Applied Mathematics

University of the Witwatersrand

Johannesburg, South Africa

Email<EMAIL_ADDRESS>
Portfolio optimisation is a multi-objective optimisation problem (MOP), where an investor aims to optimise the conflicting criteria of maximising a portfolio's expected return whilst minimising its risk and other costs. However, selecting a portfolio is a computationally expensive problem because of the cost associated with performing multiple evaluations on test data (“backtesting”) rather than solving the convex optimisation problem itself. In this research, we present ParDen, an algorithm for the inclusion of any discriminative or generative machine learning model as a surrogate to mitigate the computationally expensive backtest procedure. In addition, we compare the performance of alternative metaheuristic algorithms: NSGA-II, R-NSGA-II, NSGA-III, R-NSGA-III, U-NSGA-III, MO-CMA-ES, and COMO-CMA-ES. We measure performance using multi-objective performance indicators, including Generational Distance Plus, Inverted Generational Distance Plus and Hypervolume. We also consider meta-indicators, Success Rate and Average Executions to Success Rate, of the Hypervolume to provide more insight into the quality of solutions. Our results show that ParDen can reduce the number of evaluations required by almost a third while obtaining an improved Pareto front over the state-of-the-art for the problem of portfolio selection.

metaheuristics, genetic algorithms, surrogate modelling, multi-objective optimisation, portfolio selection, hyper-parameter optimisation

§ INTRODUCTION

Portfolio optimisation is a multi-objective optimisation problem (MOP), where an investor aims to optimise the conflicting criteria of maximising a portfolio's expected return whilst minimising its risk, and other costs [1]. The original mean-variance optimisation problem is specified by a bi-objective criterion, where risk is measured by the variance of the portfolio's return [2]. These objectives are combined linearly into a single optimisation problem:

\begin{eqnarray}
&\text{maximise} \quad &\mathbf{w}'\mathbf{\mu} - \frac{\gamma}{2}\mathbf{w}'\Sigma \mathbf{w} \\
&\text{subject to} \quad &\mathbf{w}'\mathbf{1} = 1
\end{eqnarray}

where $\mathbf{w}$ is the vector of portfolio weights, $\mathbf{\mu}$ is the vector of expected returns, $\Sigma$ is the covariance matrix of returns, and $\gamma$ is known as the risk-aversion parameter, which is used to trade-off the conflicting objectives. The resulting problem is a quadratic program that has an analytical solution. A particular choice of the risk aversion parameter yields a particular optimal solution. By varying the risk-aversion parameter and solving the optimisation problem, we can uncover the set of non-dominated solutions known as the Pareto frontier or, in the finance jargon, the efficient frontier.

In practice, portfolio selection problems are specified more extensively by incorporating a wide range of real-world considerations [1]. For example, the optimisation problem can include additional objectives for trading costs and holding costs, as well as a leverage constraint:

\begin{eqnarray}
\text{maximise} \quad &&\mathbf{w}'\mathbf{\mu} - \frac{\gamma}{2}\mathbf{w}'\Sigma \mathbf{w} - \gamma_{t}\phi^\text{trade}\left(\mathbf{w} - \mathbf{w}_0\right) \nonumber\\
&& \qquad - \gamma_h \phi^\text{hold}\left(\mathbf{w}\right) \\
\text{subject to} \quad &&\lVert\mathbf{w}\rVert_1 \leq L^\text{max}
\end{eqnarray}

Here, $\mathbf{w}_0$ is the initial portfolio, the difference $\mathbf{w} - \mathbf{w}_0$ denotes the trades required to attain portfolio $\mathbf{w}$, $\phi^\text{trade}$ and $\phi^\text{hold}\left(\mathbf{w}\right)$ are the trading and holding cost functions with corresponding trade-off parameters $\gamma_t$ and $\gamma_h$, and maximum leverage $L^{\text{max}}$.

Additional considerations can feature in the problem as hard constraints or otherwise can be included as additional objectives (i.e. soft constraints), depending on an investor's particular needs. In the latter case, these additional objectives are added to the mean-variance objective function and have associated trade-off parameters. The inclusion of hard and soft constraints implies that the mean-variance optimisation is no longer analytically tractable.

The use of metaheuristic algorithms for finding solutions to portfolio-optimisation problems has been covered extensively in the literature [3, 4, 5]. To motivate new algorithms, the convex sub-problem is contrived with appropriate constraints to make it non-convex. Nonetheless, most real-world considerations can be included in a convex specification of the optimisation problem and hence can be easily solved using Disciplined Convex Programming (DCP) [6, 7]. Here, a convex solver is used to produce an exact solution to a particular set of trade-off parameters. The Pareto frontier can be uncovered by repeatedly solving the problem for different choices of the trade-off parameters.

Motivated by hyper-parameter search, a standard practice in machine learning, Boyd et al. [1] considered the set of trade-off parameters as the hyper-parameters and the portfolio weights as the parameters in a hybrid approach. First, they specify the portfolio weights optimisation as a single-objective convex problem and use DCP to solve the corresponding exact solution. They then used a grid search over the hyper-parameters to uncover the subset of non-dominated solutions making up the Pareto-frontier. However, without prior knowledge of suitable trade-off parameters to evaluate, this process is typically computationally inefficient as many of the chosen trade-off parameters will result in sub-optimal solutions.

Nystrup [8] extends upon the work of Boyd et al. [1] by using a metaheuristic algorithm to generate the trade-off parameter set. They investigate using a single metaheuristic algorithm, the Multi-Objective Covariance-Matrix Adaptation Evolution Strategy (MO-CMA-ES), to generate new trade-off parameters to evaluate. By combining a metaheuristic and exact-solution approach, they can generate exact solutions on the frontier using far fewer evaluations than the previous grid search procedure. Other researchers have also considered the hybrid approach with similar success to both Boyd et al. [1] and Nsytrup [8], showing its utility over a purely metaheuristic approach when the sub-problem can be solved using convex optimisation [9, 10, 11, 12, 13].

In these scenarios, the computational intractability of evaluating all possible parameter combinations arises from the cost associated with performing multiple evaluations on test data (“backtesting”) rather than from solving the convex optimisation problem itself. In contrast, multi-objective algorithms can solve the Pareto frontier in a single stage by simultaneously solving both the trade-off parameters and the corresponding solutions. However, this gain in computational efficiency is at the expense of finding approximate solutions.

In this research, we adopt the hybrid approach of Boyd et al. [1] and consider extensions to reduce the number of evaluations of the backtest. We first extend the research of Nystrup [8] by comparing his usage of MO-CMA-ES against the performance of the alternative metaheuristic algorithms: NSGA-II, R-NSGA-II, NSGA-III, R-NSGA-III U-NSGA-III and COMO-CMA-ES. We then build on this by considering the inclusion of a surrogate model to mitigate the computationally expensive backtest procedure. Our results include both the original grid-search of Boyd et al. [1] as well as a random search as baselines for comparison. In order to allow for quantitative comparison, we introduce several performance indicators, including Generational Distance Plus (GD+), Inverted Generational Distance Plus (IGD+) and Hypervolume (HV). We then consider several meta-indicators of the HV to provide more insight into the quality of solutions. In particular, we consider the Success Rate (SR), the Average Evaluations to SR (AESR) and Average Generations to SR (AGSR) [14].

§ BACKGROUND

§.§ Metaheuristic Multi-Objective Optimisation Methods

There exists a significant number of metaheuristic MOP algorithms in the literature [15]. Amongst these, certain algorithms have been shown to be successful in applications across several different domains [16, 17, 18, 19, 20]. Our objective is to apply various representative metaheuristic algorithms and demonstrate the utility of a surrogate assisted metaheuristic optimisation approach to solving portfolio optimisation problems of the nature described in this paper. To this end, we implement the metaheuristic MOP methods described below.

We use Latin Hyper-cube Sampling for random selection for all our algorithms, we use real Uniform Crossover (UX) for crossover, and we use real Polynomial Mutation (PM) for mutation. For real UX, we set the probability of crossover to 90%. For PM, we set the probability of mutation to 20%. For all techniques, we use a population size of 60 and offspring size of 30. For all techniques, we give a total evaluation budget of $510$.

§.§.§ MO-CMA-ES

In multi-objective CMAES (MO-CMA-ES), a population of individuals that adapt their search strategy as in the elitist CMA-ES is maintained. The elites are subjected to multi-objective selection pressure. The selection pressure originates from non-dominated sorting using the crowding-distance first and then hypervolume as the second criterion [21]. We set sigma – the initial step size of the complete system - to $0.1$. Other parameters to CMA-ES are left as their defaults.

§.§.§ NSGA-II

The algorithm follows the form of a genetic algorithm with modified mating and survival.  The individuals are chosen Pareto front-wise first. There will arise a situation where a front needs to be split since not all individuals can survive. When splitting the front, solutions are selected based on a Manhattan crowding distance in the objective space. However, we want to hold onto the extreme points from each generation, and as a result, they are assigned a crowding distance of infinity. Furthermore, to increase selection pressure, we use a binary tournament mating selection. Each individual is first compared by rank and then crowding distance [21].

§.§.§ R-NSGA-II

The algorithm follows NSGA-II with modified survival selection. However, unlike NSGA-II, in splitting the front, solutions are selected based on rank. This rank is calculated based on the euclidean distance to a set of supplied reference points. The solution closest to a reference point is assigned a rank of one, with each point taking its best rank.
After each reference point has selected the solution with the best rank for survival, all solutions within the epsilon distance of surviving solutions are “cleared”, meaning they can not be selected for survival until all fronts have been processed. If more solutions are still needed to be selected, the cleared points are considered. The free parameter epsilon has the effect that a smaller value results in a tighter set of solutions [21]. For our work, we choose an epsilon of 0.1 by manual tuning.

§.§.§ NSGA-III

The algorithm is based on reference directions instead of reference points. For survival, after non-dominated sorting, the algorithm has a modified mechanism for dealing with the front splitting. The algorithm fills up the underrepresented reference direction first. If the reference direction does not have any solution assigned, then the solution with the smallest perpendicular distance in the normalised objective space is assigned as the survivor. In the case that additional solutions for a reference direction are required, these are assigned randomly. Consequently, when this algorithm converges, each reference direction attempts to find a suitable non-dominated solution. We use Riesz s-Energy to generate reference directions [21].

§.§.§ U-NSGA-III

It has previously been shown that tournament selection performs better than random selection. Unlike NSGA-III, which selects parents randomly for mating, this algorithm uses tournament pressure as an improvement [21].  We use Riesz s-Energy to generate reference directions [21].

§.§.§ R-NSGA-III

The algorithm follows the general NSGA-III procedure with a modified survival selection operator. First, non-dominated sorting is done as in NSGA-III. Solutions are associated with aspiration points based on perpendicular distance, then solutions from the underrepresented reference direction are chosen. For this reason, when this algorithm converges, each reference line seeks to find a good representative non-dominated solution [21].

§.§.§ COMO-CMA-ES

The algorithm is a multi-objective version of the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) single objective optimiser. Each single-objective CMA-ES optimises an indicator function given $p-1$ fixed solutions. Dominated solutions try to minimise their distance to the empirical Pareto front defined by these $p-1$ solutions. Other parameters to CMA-ES are left as their defaults [21].

§.§ Surrogate Assisted Metaheuristic Multi-Objective Optimisation

Although surrogate methods have been extensively applied within a Bayesian framework aimed at single objective optimisation [22, 23, 24], the applications of surrogate assisted methods within a multi-objective setting are substantially less common. This is especially so for data-driven metaheuristic multi-objective surrogate assisted optimisation [25, 26, 27, 28]. To this end, we present our algorithm, ParDen, which is a Pareto driven multi-objective optimisation. ParDen is sufficiently general that it can be applied together with most metaheuristic multi-objective algorithms. It also allows us to move beyond using generative models as surrogates, allowing for the future integration of discriminative models such as deep neural networks and recurrent neural networks.

ParDen relies on the use of any supervised learning model to act as a surrogate drop-in replacement $\hat{f}(\cdot)$ to evaluate a computationally resource intensive real-world experiment or simulation $f_e(\cdot)$. In this research, the “backtesting” procedure, which runs for several minutes per a set of trade-off parameters, is the simulation under consideration.

The core idea of ParDen is to limit our evaluation of candidates using the simulation $f_e(\cdot)$ to only those that have a high probability of advancing the Pareto front (the non-dominated solution set $\mathcal{P}$). The likelihood that a candidate would advance the Pareto front is determined by two things: i) whether the surrogate model predicts that the candidate advances the Pareto front and ii) the extent to which our surrogate is correct. These two criteria allow us to balance exploration and exploitation of the surrogate. In particular, we attempt to estimate the probability that the surrogate is correct, and we use this probability as a threshold within a rejection sampling scheme to select candidates.

In order to determine the extent to which our surrogate is correct, we require a mechanism for assessing the probability $P[(x_i,y_i) \in \mathcal{P} | \hat{f}]$ that a pretender $(x_i, y_i)$ will be in the current non-dominated set $\mathcal{P}$ given the surrogate's prediction $\hat{y}_i = \hat{f}(\cdot)$. We employ $k$-fold cross-validation $\operatorname{CV}_k(\cdot)$ together with the non-dominating rank [29] of each candidate to estimate this probability as a non-dominated score ($\operatorname{NDScore}$).

To calculate a $\operatorname{NDScore}$ we fit surrogates $\hat{f}_{_{X_T,Y_T}}(\cdot)$ using the training splits $(X_T,Y_T)$ with the multi-objective $Y_T$ as the target. For each point in the validation split $(X_V,Y_V)$, we assign its non-dominating rank in the split as its class label $r$. We use the surrogates to estimate $\hat{Y}_V=\hat{f}_{_{X_T,Y_T}}(X_V)$, the objective values of the validation split. We then assign the non-dominating rank of the estimated values $\hat{Y}_V$ as the predicted class labels $\hat{r}$. We can use any classification metric $\mathcal{E}(\hat{r}, r)$ constrained to sum to one, to compare $r$ with $\hat{r}$:

\begin{equation}
\operatorname{NDScore}(\cdot) = \operatorname*{aggregation}_{\substack{(X_T,Y_T),(X_V,Y_V)\\ \in \operatorname{CV}_k(X)}} \mathcal{E}\left(\hat{r}, r\right) \\
\end{equation}

\begin{equation}
\hat{r} = \lceil Y_V \rceil^{\uparrow} \quad \mathrm{and} \quad
r = \lceil\hat{f}_{_{X_T,Y_T}}(X_V) \rceil^{\uparrow}
\end{equation}

where $\lceil\cdot\rceil^\uparrow$ assigns the non-dominating rank to each vector. Since we are using $k$-fold cross-validation to estimate this value and we are interested in the worst case, we use the minimum across all $k$ runs as our aggregation. We note that accuracy could be replaced with a more balanced metric such as $F_1$-Score and that using the mean instead of minimum might be a better aggregation over $CV_k$. We did not investigate this further.

ParDen Algorithm
Metaheuristic ($M$), Terminating Criterion ($TC$), Simulation ($f_e$), Surrogate Class ($\mathcal{H}$), Loss ($L$)
Approximate Optimal Frontier ($\mathcal{P}$)
Daily Trading Data ($D$)

$M, f_e, \ldots$
$X_W \sim M$                *[r]Generate warm start candidates
$Y_W=f_e(X_W,D)$              *[r]Actual simulated fitness
$\mathcal{G}=\{(X_W,Y_W)\}$ *[r]Set ground-truth
$\mathcal{P} = \lceil \mathcal{G} \rceil$ *[r]Compute non-dominated set
$\operatorname{NDScore} = 0$  *[r]Set non-dominated score
*[r]Train surrogate on $\mathcal{G}$ with loss $L$
$\hat{f} = \operatorname{argmin}_{h \in \mathcal{H}} \mathbb{E}[L(h,\mathcal{G})]$
$\hat{f}, \mathcal{G}, \operatorname{NDScore}, \mathcal{P}$

$M, f_e, \ldots$
$X_C \sim M$                *[r]Generate new candidates
*[r]Estimate candidates' fitness with surrogate
$\hat{Y}_C=\hat{f}(X_C, D)$
*[r]Join candidates to non-dominated set
$\mathcal{P}_C = \lceil \{\mathcal{C} \cup \mathcal{P}\} \rceil$  where $\mathcal{C}=\{(X_C,\hat{Y}_C)\}$
*[r]Pretenders are non-dominated candidates
($X_S, \hat{Y}_S) = \{\mathcal{P}_C \cup \mathcal{C}\}$
*[r]Acceptance sampling with $\operatorname{NDScore}$ as threshold to add additional pretenders
$X_S \quad\operatorname{append}\quad X_C$ if $\operatorname{NDScore} \le x \sim [0,1]$
$Y_S=f_e(X_S,D)$      *[r]Actual pretenders' fitness
$\mathcal{G} = \lbrace\mathcal{G}\cup (X_S,Y_S) \rbrace$ *[r]Add to ground-truth
$\mathcal{P} = \lceil \mathcal{G} \rceil$ *[r]Update non-dominated set
*[r]Update non-dominated score
$\operatorname{NDScore} = \operatorname{NDScore}(\hat{f}, \mathcal{G})$
*[r]Train new surrogate on $\mathcal{G}$ with loss $L$
$\hat{f} = \operatorname*{argmin}_{h \in \mathcal{H}} \mathbb{E}[L(h,\mathcal{G})]$
$\hat{f}, \mathcal{G}, \operatorname{NDScore}, \mathcal{P}$

($M, f_e$)
$TC$ = False

($M, f_e, \ldots$)

$X \sim I$: sample $X$ from candidate generator $I$
$\lceil \mathcal{X} \rceil$: Return non-dominated set of $\mathcal{X}$

startstop = [rectangle, rounded corners, minimum width=2cm, minimum height=1cm,text centered, draw=ashgrey, drop shadow, fill=ashgrey!50,]
io = [trapezium, trapezium left angle=70, trapezium right angle=110, minimum width=2cm, minimum height=1cm, text centered, draw=black, drop shadow]
process = [rectangle, minimum width=4.5cm, minimum height=1cm, text centered, text width=4cm, draw=orange, drop shadow, fill = orange!20]
process2 = [rectangle, minimum width=2cm, minimum height=1cm, text centered, text width=2cm, draw=black, drop shadow, fill = white]
process3 = [rectangle, minimum width=4.0cm, minimum height=1cm, text centered, text width=4.0cm, draw=orange, drop shadow, fill = orange!20]
process4 = [rectangle, minimum width=2.0cm, minimum height=1cm, text centered, text width=2.0cm, draw=orange, drop shadow, fill = orange!20]
decision = [diamond, minimum width=3cm, minimum height=1cm, text centered, draw=darkred, drop shadow, fill = darkred!20]
arrow = [thick,->,>=stealth]
Tank/.style=draw=black,fill=white,thick,rectangle,rounded corners=10pt,minimum width=0.75cm,minimum height=1.5cm,text width=0.75cm,align=center, drop shadow
surround = [fill=gray!15,thick,draw=black,rounded corners=2mm]

line = [draw, -stealth]

[node distance=1cm,every node/.style=fill=white, align=center]

(start) [startstop,xshift=3cm] Start;

(pro1) [process, below of =start, yshift=-1cm] $X_W \sim M$;
(pro2) [process, below of =pro1, yshift=-1cm] $Y_W=f_e(X_W,D)$;
(pro3) [process, below of =pro2, yshift=-1cm] $\mathcal{G}=\{(X_W,Y_W)\}$;
(pro4) [process, below of =pro3, yshift=-1cm] $\mathcal{P} = \lceil \mathcal{G} \rceil$;
(pro5) [process, below of =pro4, yshift=-1cm] $\operatorname{NDScore} = 0$;
(pro6) [process, below of =pro5, yshift=-1cm] $\hat{f} = \operatorname{argmin}_{h \in \mathcal{H}} \mathbb{E}[L(h,\mathcal{G})]$;

[surround] (background) [inner sep=0.5cm, fit = (pro1) (pro6), label=[rotate=270, anchor=south]right:] ;

(dec2) [decision, below of=pro5, yshift=-5cm] $TC$ = False;

(main1) [process3, left of =pro1, xshift=-6cm] $X_C \sim M$;
(main2) [process3, below of =main1, yshift=-1cm] $\hat{Y}_C=\hat{f}(X_C, D)$;
(main3) [process3, below of =main2, yshift=-1cm] $\mathcal{P}_C = \lceil \{\mathcal{C} \cup \mathcal{P}\} \rceil$;
(main4) [process3, below of =main3, yshift=-1cm] ($X_S, \hat{Y}_S) = \{\mathcal{P}_C \cup \mathcal{C}\}$;
(main5) [process3, below of =main4, yshift=-1cm] $X_S \operatorname{append} X_C$ if $\operatorname{NDScore} \le x \sim [0,1]$;
(main6) [process3, below of =main5, yshift=-1cm] $Y_S=f_e(X_S,D)$;
(main7) [process3, below of =main6, yshift=-1cm] $\mathcal{G} = \lbrace\mathcal{G}\cup (X_S,Y_S) \rbrace$;
(main8) [process3, below of =main7, yshift=-1cm] $\mathcal{P} = \lceil \mathcal{G} \rceil$;
(main9) [process3, below of =main8, yshift=-1cm] $\operatorname{NDScore} = \operatorname{NDScore}(\hat{f}, \mathcal{G})$;
(main10) [process3, below of =main9, yshift=-1cm] $\hat{f} = \operatorname*{argmin}_{h \in \mathcal{H}} \mathbb{E}[L(h,\mathcal{G})]$;
[surround] (background2) [inner sep=0.5cm, fit = (main1) (main10), label=[rotate=90, anchor=south]left:] ;

(out1) [process4, right of =dec2, xshift=2cm, yshift=-2cm] $\mathcal{P}$;

[arrow] (start) – (pro1);
[arrow] (pro1) – (pro2);
[arrow] (pro2) – (pro3);
[arrow] (pro3) – (pro4);
[arrow] (pro4) – (pro5);
[arrow] (pro5) – (pro6);

[arrow] (pro6) – (dec2);
[arrow] (dec2.west) node[below, xshift=-0.5cm, yshift=-0.25cm]No |- ++ (-2.0cm, -0cm) |- (main1);

[arrow] (main10) |- ++(0, -15mm) -| (dec2.south);

[arrow] (dec2.east) node[below, xshift=0.5cm, yshift=0.75cm]Yes -| (out1);

[arrow] (main1) – (main2);
[arrow] (main2) – (main3);
[arrow] (main3) – (main4);
[arrow] (main4) – (main5);
[arrow] (main5) – (main6);
[arrow] (main6) – (main7);
[arrow] (main7) – (main8);
[arrow] (main8) – (main9);
[arrow] (main9) – (main10);

Flowchart for the ParDen Algorithm

§ METHODOLOGY

In order to demonstrate the efficacy of the ParDen algorithm, we first evaluate several MOP EAs on our problem to see which of them performs best and provide a baseline. We then test ParDen in the two best performing EAs: MO-CMA-ES and NSGA-II. Since we cannot clearly ascertain a single best EA, we instead select the best one with a fixed number of points in the front and one with a dynamic number of points in the front. Finally, we compare all of these results using the following quality metrics on the HV score.

Success Rate (SR): Success is defined as a required quality criterion that needs to be achieved, i.e. a solution within 99% or 95% of the known optimal value. The success rate as a metric is defined as the percentage of algorithm executions in which the success criterion is met.

Average Evaluations to a Solution (AES): Average number of evaluations of the simulation required to reach a solution at a given SR, i.e. the AES at a SR of 99%, would be the average number of evaluations it took to reach there. The particular runs of the EA which do not meet the SR criterion are ignored.

Average Generations to Success Rate (AGSR): Average number of generations of the simulation required to reach the solution at a given SR, i.e. the AGSR at a SR of 99%, would be the average number of generations required to reach there.

§ RESULTS AND DISCUSSION

Performance Indicators

Method                                          GD$+$          IGD$+$         HV

Random $6510$                                   -              -              1313.0

Grid Search [1]                .5941           .5703         1259.1

NSGA-II                                         .1137           .1585         1308.4

R-NSGA-II                                       .0259      .2443         1309.7

NSGA-III                                        .1062           .2320         1305.0

R-NSGA-III                                      .1215           .2279         1304.3

U-NSGA-III                                      .1314           .2567         1302.7

COMO-CMA-ES                                     .0413          1.775          1141.1

MO-CMA-ES [8]      .1002           .0965         1317.4

$\mathcal{H}$:MO-CMA-ES*                        .0912          .0951     1320.0

$\mathcal{H}$:NSGA-II*                          .0999           .1141         1316.1

* ParDen algorithm applied

We note from Table <ref> that all of the EAs perform similarly when comparing them using hypervolume. The difficulty with using hypervolume as a metric is that it is sensitive to the number of points in the Pareto front; thus, methods with a higher number of points have an advantage when compared via hypervolume. Next, we turn our attention to GD+ and IGD+ but are cognizant that neither of these is being measured with respect to the true Pareto front but instead measured against the Random $6510$ points as a reference. Fortunately, both GD+ and IGD+ are semi Pareto respecting, meaning they do not negatively impact the score of methods that outperform the reference Pareto front. In this respect, we note that although COMO-CMA-ES has one of the best GD+ scores, this has come at the cost of excessive crowing as seen in Figure <ref> and reinforced by the IGD+ score. Finally, we note that the surrogate assisted  $\mathcal{H}$:MO-CMA-ES has the best performing IGD+ score and the second-best GD+. Combined with the fact that it has the best HV, this indicates that our ParDen algorithm can outperform all other methods.

Quality Indicators

Method                                          SR$@99$     SR$@95$     AESR$@99$     AESR$@95$     AGSR$@99$

Grid Search [1]                0.0         100.0       -              510.0         -

NSGA-II                                         100.0       100.0       327.0          63.0     9.9

R-NSGA-II                                       100.0       100.0       297.0          63.0     8.9

NSGA-III                                        100.0       100.0       405.0          81.0          12.5

R-NSGA-III                                      60.0        100.0       415.3          75.2          12.8

U-NSGA-III                                      70.0        100.0       338.6          96.0          10.3

COMO-CMA-ES                                     0.0         0.0         -              -             -

MO-CMA-ES [8]      100.0       100.0       333.0          99.0          10.1

$\mathcal{H}$:MO-CMA-ES*                        100.0       100.0       247.7     75.9          10.7

$\mathcal{H}$:NSGA-II*                          100.0       100.0       290.2          120.0         10.7

* ParDen algorithm applied

[group style=
group name=myplot,
group size= 2 by 5,height=4.75cm,width=8.75cm]
[ylabel=Return %,
grid style=gray!30, loosely dashed,
legend entries=Random 6510,Grid Search [1],
legend style=at=(.9875,0.025),anchor=south east,
legend cell align=left,
axis background/.style=fill=gray!7]

+[only marks,mark size = 1.75, mark options=
draw = black!90,
fill = black!70,
fill opacity=0.6,
draw opacity=0.3,
col sep=comma,
+[only marks,mark=square*, mark size = 1.75,mark options=
draw = Set1-B!90,
fill = Set1-B!70,
col sep=comma,
legend entries=Random 6510,NSGA-II,
legend style=at=(.9875,0.025),anchor=south east,
legend cell align=left,
grid style=gray!30, loosely dashed,
axis background/.style=fill=gray!7

+[only marks,mark size = 1.75, mark options=
draw = black!90,
fill = black!70,
fill opacity=0.6,
draw opacity=0.3,
col sep=comma,
+[only marks,mark=diamond*, mark size = 1.75,mark options=
draw = Dark2-A!90,
fill = Dark2-A!70,
col sep=comma,
[ylabel=Return %,
legend entries=Random 6510,R-NSGA-II,
legend style=at=(.9875,0.025),anchor=south east,
legend cell align=left,
grid style=gray!30, loosely dashed,
axis background/.style=fill=gray!7]
+[only marks,mark size = 1.75, mark options=
draw = black!90,
fill = black!70,
fill opacity=0.6,
draw opacity=0.3,
col sep=comma,
+[only marks,mark=triangle*, mark size = 1.75,mark options=
draw = Set1-A!90,
fill = Set1-A!70,
col sep=comma,
legend entries=Random 6510,NSGA-III,
legend style=at=(.9875,0.025),anchor=south east,
legend cell align=left,
grid style=gray!30, loosely dashed,
axis background/.style=fill=gray!7]
+[only marks,mark size = 1.75, mark options=
draw = black!90,
fill = black!70,
fill opacity=0.6,
draw opacity=0.3,
col sep=comma,
+[only marks,mark=oplus*, mark size = 1.75,mark options=
draw = Dark2-B!90,
fill = Dark2-B!70,
col sep=comma,
[ylabel=Return %,
legend entries=Random 6510,R-NSGA-III,
legend style=at=(.9875,0.025),anchor=south east,
legend cell align=left,
grid style=gray!30, loosely dashed,
axis background/.style=fill=gray!7]
+[only marks,mark size = 1.75, mark options=
draw = black!90,
fill = black!70,
fill opacity=0.6,
draw opacity=0.3,
col sep=comma,
+[only marks,mark=halfsquare left*, mark size = 1.75,mark options=
draw = Dark2-D!90,
fill = Dark2-D!70,
col sep=comma,
legend entries=Random 6510,U-NSGA-III,
legend style=at=(.9875,0.025),anchor=south east,
legend cell align=left,
grid style=gray!30, loosely dashed,
axis background/.style=fill=gray!7]
+[only marks,mark size = 1.75, mark options=
draw = black!90,
fill = black!70,
fill opacity=0.6,
draw opacity=0.3,
col sep=comma,
+[only marks,mark=halfsquare right*, mark size = 1.75,mark options=
draw = Set1-B!90,
fill = Set1-B!70,
col sep=comma,
[ylabel=Return %,
legend entries=Random 6510,COMO-CMA-ES,
legend style=at=(.9875,0.025),anchor=south east,
legend cell align=left,
grid style=gray!30, loosely dashed,
axis background/.style=fill=gray!7]
+[only marks,mark size = 1.75, mark options=
draw = black!90,
fill = black!70,
fill opacity=0.6,
draw opacity=0.3,
col sep=comma,
+[only marks,mark=otimes*, mark size = 1.75,mark options=
draw = Dark2-F!90,
fill = Dark2-F!70,
col sep=comma,
legend entries=Random 6510,MO-CMA-ES [8],
legend style=at=(.9875,0.025),anchor=south east,
legend cell align=left,
grid style=gray!30, loosely dashed,
axis background/.style=fill=gray!7]
+[only marks,mark size = 1.75, mark options=
draw = black!90,
fill = black!70,
fill opacity=0.6,
draw opacity=0.3,
col sep=comma,
+[only marks,mark=halfcircle*, mark size = 1.75,mark options=
draw = Dark2-G!90,
fill = Dark2-G!70,
col sep=comma,
[xlabel=Risk %,
ylabel=Return %,
legend entries=Random 6510,$\mathcal{H}$:MO-CMA-ES,
legend style=at=(.9875,0.025),anchor=south east,
legend cell align=left,
grid style=gray!30, loosely dashed,
axis background/.style=fill=gray!7

+[only marks,mark size = 1.75, mark options=
draw = black!90,
fill = black!70,
fill opacity=0.6,
draw opacity=0.3,
col sep=comma,
+[only marks,mark=halfdiamond*, mark size = 1.75,mark options=
draw = Set1-A!90,
fill = Set1-A!70,
col sep=comma,
[xlabel=Risk %,
legend entries=Random 6510,$\mathcal{H}$:NSGA-II,
legend style=at=(.9875,0.025),anchor=south east,
legend cell align=left,
grid style=gray!30, loosely dashed,
axis background/.style=fill=gray!7]
+[only marks,mark size = 1.75, mark options=
draw = black!90,
fill = black!70,
fill opacity=0.6,
draw opacity=0.3,
col sep=comma,
+[only marks,mark=pentagon*, mark size = 1.75,mark options=
draw = Set1-C!90,
fill = Set1-C!70,
col sep=comma,

Optimal Frontiers for each Method in Table <ref>

every non boxed x axis/.append style=x axis line style=-

xmin  = -1,  xmax = 24.5,
ymin  =  1180,
ymax = 1330.0,
ylabel=Hypervolume (HV),
xlabel=Generation Number,
axis lines = left,
axis background/.style=fill=gray!7, draw=gray,
[Dark2-A, ultra thick]
table [x = it, y = val, col sep=comma] convergence_hvs.csv;

axis y line=right,
axis x line=none,

ymin = -0.0,
ymax = 1.3,
xmin = -1,
xmax = 24.5,

ylabel= GD+/IGDI+,
legend cell align=left,
legend style=at=(0.99,0.35),
anchor=south east,legend columns=1,

/pgfplots/refstyle=hvHypervolume (HV)

[Set1-B, ultra thick]
table [ x = it, y = val, col sep=comma] convergence_gdps.csv ;
Generational Distance Plus (GD+)

[Set1-B, ultra thick, dashed]
table [ x = it, y = val, col sep=comma] convergence_igdps.csv ;
Inverted Generational Distance Plus (IGD+)

Convergence Plots of HV, GD+ and IGD+ Metrics for $\mathcal{H}$:MO-CMA-ES over Successive Generations

every non boxed x axis/.append style=x axis line style=-stealth

xlabel=Risk %,
ylabel=Return %,
axis lines = left,
legend cell align=left,
legend style=at=(0.99,0.05),
anchor=south east,legend columns=1,
axis background/.style=fill=gray!7, draw=gray

+[only marks, mark=halfdiamond*, color=Set1-A, opacity=1] table[col sep=comma,x=y22,y=x22]convergence_pfs.csv;
+[only marks, mark=halfdiamond*, color=Set1-A,opacity=0.1+0.01*#1] table[col sep=comma,x=y#1,y=x#1]

Descending Opacity Illustrating the Advancing Optimal Frontier of $\mathcal{H}$:MO-CMA-ES from Generation to Generation

Although the ParDen algorithm can achieve superior performance with respect to HV, GD+, and IGD+, of interest is the quality of the algorithm in arriving at that solution. From Table <ref>, we note that for both surrogate assisted methods, $\mathcal{H}$:MO-CMA-ES and $\mathcal{H}$:NSGA-II, the number of evaluations required to reach a solution within 99% of random $6510$ reference HV is noticeably lower. In fact, $\mathcal{H}$:MO-CMA-ES requires almost one third or approximately $100$ fewer evaluations. Turing to Figure <ref> and Figure <ref> we see the rapid initial gains in HV, GD+ and IGD+ for $\mathcal{H}$:MO-CMA-ES. Again, these results lend support to the efficacy of the ParDen algorithm, especially in cases where there is a realisable cost associated with each execution of the simulation.

In previous work, Nydstrup et al. [8] mentions that one need only consider the number of evaluations. However, this is not entirely correct as EAs improve their solutions from generation to generation. As a result the total evaluations of the simulation must be considered with respect to incremental improvements between generations. To this end, we consider the number of generations required as a proxy for Efficiency. It is meaningful to make this assumption as all evaluations within a generation could be done in parallel and the time taken for an evaluation of the backtest is orders of magnitude larger than any of the operations performed by optimisation algorithms. Inspecting the Average number of Generations to a SR$@99$ (AGSR$@99$) we note that, although the surrogates assisted models can achieve better results with fewer evaluations of the simulation, they still require the same number of generations. This observation is significant since it highlights the need for further investigations into algorithms that are able to reduce the number of generations required to reach optimal solutions.

§ CONCLUSION

In this paper, we have extended the existing research by applying a number of additional EA algorithms to the constrained multi-objective portfolio optimisation problem. We further built upon this work by introducing ParDen; a Surrogate assisted algorithm that, together with an NDScore, allows any generative or discriminative machine learning model to be used. Our results show the superiority of our ParDen framework across several performance and quality metrics. Our research highlights the need for further investigation into algorithms better suited to problems with computationally resource intensive objective evaluations.

[1]
S. Boyd, E. Busseti, S. Diamond, R. N. Kahn, K. Koh, P. Nystrup, J. Speth
et al., “Multi-period trading via convex optimization,”
Foundations and Trends® in Optimization, vol. 3,
no. 1, pp. 1–76, 2017.

[2]
H. Markowitz, “Portfolio selection,” The Journal of Finance, vol. 7,
no. 1, pp. 77–91, 1952.

[3]
J. Doering, R. Kizys, A. A. Juan, À. Fitó, and O. Polat,
“Metaheuristics for rich portfolio optimisation and risk management: Current
state and future trends,” Operations Research Perspectives, vol. 6,
p. 100121, 2019.

[4]
N. Loukeris, D. Donelly, A. Khuman, and Y. Peng, “A numerical evaluation of
meta-heuristic techniques in portfolio optimisation,” Operational
Research, vol. 9, no. 1, pp. 81–103, 2009.

[5]
E. Fernandez, C. Gomez, G. Rivera, and L. Cruz-Reyes, “Hybrid metaheuristic
approach for handling many objectives and decisions on partial support in
project portfolio optimisation,” Information Sciences, vol. 315, pp.
102–122, 2015.

[6]
M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex
programming, version 2.1,” <http://cvxr.com/cvx>, Mar. 2014.

[7]
M. Grant, S. Boyd, and Y. Ye, “Disciplined convex programming,” in
Global optimization.1em plus 0.5em minus 0.4emSpringer,
2006, pp. 155–210.

[8]
P. Nystrup, E. Lindström, and H. Madsen, “Hyperparameter optimization for
portfolio selection,” The Journal of Financial Data Science, vol. 2,
no. 3, pp. 40–54, 2020.

[9]
J. Branke, B. Scheckenbach, M. Stein, K. Deb, and H. Schmeck, “Portfolio
optimization with an envelope-based multi-objective evolutionary algorithm,”
European Journal of Operational Research, vol. 199, no. 3, pp.
684–693, 2009.

[10]
R. Ruiz-Torrubiano and A. Suarez, “Hybrid approaches and dimensionality
reduction for portfolio selection with cardinality constraints,” IEEE
Computational Intelligence Magazine, vol. 5, no. 2, pp. 92–107, 2010.

[11]
R. Ruiz-Torrubiano and A. Suárez, “A memetic algorithm for
cardinality-constrained portfolio optimization with transaction costs,”
Applied Soft Computing, vol. 36, pp. 125–142, 2015.

[12]
A. Baykasoğlu, M. G. Yunusoglu, and F. B. Özsoydan, “A grasp based
solution approach to solve cardinality constrained portfolio optimization
problems,” Computers & Industrial Engineering, vol. 90, pp.
339–351, 2015.

[13]
R. Qi and G. G. Yen, “Hybrid bi-objective portfolio optimization with
pre-selection strategy,” Information Sciences, vol. 417, pp.
401–419, 2017.

[14]
A. Eiben and J. Smith, Introduction to Evolutionary Computing.1em
plus 0.5em minus 0.4emSpringer, 2015.

[15]
Q. Liu, X. Li, H. Liu, and Z. Guo, “Multi-objective metaheuristics for
discrete optimization problems: A review of the state-of-the-art,”
Applied Soft Computing, p. 106382, 2020.

[16]
K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist
multiobjective genetic algorithm: Nsga-ii,” IEEE transactions on
evolutionary computation, vol. 6, no. 2, pp. 182–197, 2002.

[17]
Q. Zhang and H. Li, “Moea/d: A multiobjective evolutionary algorithm based on
decomposition,” IEEE Transactions on evolutionary computation,
vol. 11, no. 6, pp. 712–731, 2007.

[18]
C. Igel, N. Hansen, and S. Roth, “Covariance matrix adaptation for
multi-objective optimization,” Evolutionary computation, vol. 15,
no. 1, pp. 1–28, 2007.

[19]
M. Helbig and A. P. Engelbrecht, “Population-based metaheuristics for
continuous boundary-constrained dynamic multi-objective optimisation
problems,” Swarm and Evolutionary computation, vol. 14, pp. 31–47,

[20]
H. AlSattar, A. Zaidan, B. Zaidan, M. A. Bakar, R. Mohammed, O. Albahri,
M. Alsalem, and A. Albahri, “Mogsabat: A metaheuristic hybrid algorithm for
solving multi-objective optimisation problems,” Neural Computing and
Applications, vol. 32, no. 8, pp. 3101–3115, 2020.

[21]
J. Blank and K. Deb, “pymoo: Multi-objective optimization in python,”
IEEE Access, vol. 8, pp. 89 497–89 509, 2020.

[22]
M. Ben Salem, O. Roustant, F. Gamboa, and L. Tomaso, “Universal prediction
distribution for surrogate models,” SIAM/ASA Journal on Uncertainty
Quantification, vol. 5, no. 1, pp. 1086–1109, 2017.

[23]
Z. Zhou, Y. S. Ong, P. B. Nair, A. J. Keane, and K. Y. Lum, “Combining global
and local surrogate models to accelerate evolutionary optimization,”
IEEE Transactions on Systems, Man, and Cybernetics, Part C
(Applications and Reviews), vol. 37, no. 1, pp. 66–76, 2006.

[24]
X. Wan, J. F. Pekny, and G. V. Reklaitis, “Simulation-based optimization with
surrogate models—application to supply chain management,” Computers
& chemical engineering, vol. 29, no. 6, pp. 1317–1328, 2005.

[25]
H. Wang and Y. Jin, “A random forest-assisted evolutionary algorithm for
data-driven constrained multiobjective combinatorial optimization of trauma
systems,” IEEE transactions on cybernetics, vol. 50, no. 2, pp.
536–549, 2018.

[26]
C. Yang, J. Ding, Y. Jin, and T. Chai, “Offline data-driven multiobjective
optimization: Knowledge transfer between surrogates and generation of final
solutions,” IEEE Transactions on Evolutionary Computation, vol. 24,
no. 3, pp. 409–423, 2019.

[27]
T. Chugh, N. Chakraborti, K. Sindhya, and Y. Jin, “A data-driven
surrogate-assisted evolutionary algorithm applied to a many-objective blast
furnace optimization problem,” Materials and Manufacturing Processes,
vol. 32, no. 10, pp. 1172–1178, 2017.

[28]
L. Stander, M. Woolway, and T. van Zyl, “Data-driven evolutionary optimisation
for the design parameters of a chemical process: A case study,” in
2020 IEEE 23rd International Conference on Information Fusion
(FUSION).1em plus 0.5em minus 0.4emIEEE, 2020, pp. 1–8.

[29]
S. T. Selvi, S. Baskar, and S. Rajasekar, “Application of evolutionary
algorithm for multiobjective transformer design optimization,” in
Classical and Recent Aspects of Power System Optimization.1em
plus 0.5em minus 0.4emElsevier, 2018, pp. 463–504.


$\displaystyle=-16\pi^{2}\alpha_{k}\rho\mu_{\rho,n}\sum_{l=1}^{m}\sum_{i=1,2}\frac{\tau_{l}s_{i,n}^{(l)}}{\sqrt{8}}\frac{\partial^{2}\mathcal{F}}{\partial
x_{i}^{(l)}\partial x_{j}^{(k)}}(\xi_{1},\ldots,\xi_{m})+o(\rho),$

Hence

$\displaystyle-\frac{16}{3\sqrt{8}\rho}\pi\frac{\alpha_{k}}{\tau_{k}}s_{j,n}^{(k)}(1+o(1))(1-\mu_{\rho,n})$
$\displaystyle=-16\pi^{2}\alpha_{k}\rho\mu_{\rho,n}\sum_{l=1}^{m}\sum_{i=1,2}\frac{\tau_{l}s_{i,n}^{(l)}}{\sqrt{8}}\frac{\partial^{2}\mathcal{F}}{\partial
x_{i}^{(l)}\partial x_{j}^{(k)}}(\xi_{1},\ldots,\xi_{m})+o(\rho).$

Having set

$\eta_{n}:=\frac{\sum_{l=1}^{m}\sum_{i=1,2}\tau_{l}\tau_{k}s_{i,n}^{(l)}\frac{\partial^{2}\mathcal{F}}{\partial
x_{i}^{(l)}\partial x_{j}^{(k)}}(\xi_{1},\ldots,\xi_{m})}{s_{j,n}^{(k)}},$

we get

$s_{j,n}^{(k)}\eta_{n}=\sum_{l=1}^{m}\sum_{i=1,2}\tau_{l}\tau_{k}s_{i,n}^{(l)}\frac{\partial^{2}\mathcal{F}}{\partial
x_{i}^{(l)}\partial x_{j}^{(k)}}(\xi_{1},\ldots,\xi_{m}),$

so that
$s_{n}=\left(s_{1,n}^{(1)},\ldots,s_{1,n}^{(m)},s_{2,n}^{(1)},\ldots,s_{2,n}^{(m)}\right)$
and $\eta_{n}$ are respectively the eigenvectors and the eigenvalues of the
matrix $D(\mathrm{Hess}\,\mathcal{F})D$. We can therefore conclude that

$(1-\mu_{\rho,n})=3\pi\rho^{2}\eta_{n}+o(\rho^{2}).$

∎

Finally we can prove the following

###### Proof of Theorem 2.

.

Proposition 23 proves point 1 and point 3, while Theorem 17 jointly with point
(i) of Proposition 20 proves point 2. ∎

And

###### Proof of Theorem 3.

.

With similar calculus of those in Proposition 20 we can prove that for
$3m+1<j\leqslant 4m$, $\mu_{\rho,j}$ still goes to one. In the proof of point
(iii) in Proposition 20 we showed that for $3m+1<j\leqslant 4m$,
$t^{(k)}_{j}\neq 0$, so Lemma 14 holds and in such a case at the end of
Theorem 17 we proved (4.35), therefore point 1 and point 3 are proved.
Finally, considering Proposition 12 and the fact that the functions
$\frac{8-|x|^{2}}{8+|x|^{2}}$ and $\frac{x_{i}}{8+|x|^{2}}$, for $i=1,2$, are
orthogonal and that eigenfunctions relative to different eigenspaces have to
be orthogonal we get point 2. ∎

We are now ready to prove the main result.

###### Proof of Theorem 4.

.

For $\rho$ sufficiently small, by Proposition 10 we have that $\mu_{\rho,j}<1$
for every $j=1,\ldots,m$. By Proposition 20 we have that $\mu_{\rho,j}>1$ for
$j>3m$ and for every $m+1\leqslant j\leqslant 3m$ we have that
$\mu_{\rho,j}\leqslant 1+C\rho^{2}$, so by Theorem 17 we get
$t_{j}=\textbf{0}$. Therefore, since the hypotheses of Proposition 23 are
satisfied, we get

$\mu_{\rho,j}=1-3\pi\rho^{2}\eta_{j}(1+o(1)),$

where the $\eta_{j}$’s, for $j=1,\ldots,2m$, are the eigenvalues of the matrix
$D(\mathrm{Hess}\,\mathcal{F})D$. Then, to any positive eigenvalue $\eta_{j}$
of the matrix $D(\mathrm{Hess}\,\mathcal{F})D$ there corresponds a negative
$\mu_{\rho,j}$. Since $D$ is a diagonal, positive definite matrix, the
signature of the matrix $D(\mathrm{Hess}\,\mathcal{F})D$ is equal to the
signature of the matrix $\mathrm{Hess}\,\mathcal{F}$, and hence the theorem is
proved. ∎

We deduce the following corollary.

###### Proof of Corollary 5.

.

It is enough to notice that $\mathrm{Hess}\,\mathcal{F}$ is a $2m\times 2m$
matrix, and so $0\leqslant\mathcal{M}(\mathrm{Hess}\,\mathcal{F})\leqslant
2m$. ∎

## References

  * [1] D. Bartolucci and A. Pistoia.  Existence and qualitative properties of concentrating solutions for the sinh-Poisson equation.  IMA J. Appl. Math., 72(6):706–729, 2007.
  * [2] T. Bartsch, K.-C. Chang, and Z.-Q. Wang.  On the morse indices of sign changing solutions of nonlinear elliptic problems.  Mathematische Zeitschrift, 233(4):655–677, 2000.
  * [3] T. Bartsch, A. Pistoia, and T. Weth.  N-vortex equilibria for ideal fluids in bounded planar domains and new nodal solutions of the sinh-poisson and the lane-emden-fowler equations.  Communications in Mathematical Physics, 297(3):653–686, 2010.
  * [4] E. Caglioti, P.-L. Lions, C. Marchioro, and M. Pulvirenti.  A special class of stationary flows for two-dimensional Euler equations: a statistical mechanics description.  Communications in Mathematical physics, 143(3):501–525, 1992.
  * [5] E. Caglioti, P.-L. Lions, C. Marchioro, and M. Pulvirenti.  A special class of stationary flows for two-dimensional Euler equations: a statistical mechanics description. Part II.  Communications in Mathematical Physics, 174(2):229–260, 1995.
  * [6] A. J. Chorin.  Vorticity and turbulence, volume 103.  Springer Science & Business Media, 2013.
  * [7] A. J. Chorin, J. E. Marsden, and J. E. Marsden.  A mathematical introduction to fluid mechanics, volume 3.  Springer, 1990.
  * [8] K. Chow, N. Ko, R. Leung, and S. Tang.  Inviscid two dimensional vortex dynamics and a soliton expansion of the sinh-Poisson equation.  Physics of fluids, 10(5):1111–1119, 1998.
  * [9] K. Chow, S. Tsang, and C. Mak.  Another exact solution for two-dimensional, inviscid sinh Poisson vortex arrays.  Physics of Fluids, 15(8):2437–2440, 2003.
  * [10] T. Dauxois, S. Fauve, and L. Tuckerman.  Stability of periodic arrays of vortices.  Physics of Fluids, 8(2):487–495, 1996.
  * [11] F. De Marchis, I. Ianni, and F. Pacella.  Exact Morse index computation for nodal radial solutions of Lane–Emden problems.  Mathematische Annalen, 367(1-2):185–227, 2017.
  * [12] K. El Mehdi and M. Grossi.  Asymptotic estimates and qualitative properties of an elliptic problem in dimension two.  Advanced Nonlinear Studies, 4(1):15–36, 2004.
  * [13] P. Esposito, M. Grossi, and A. Pistoia.  On the existence of blowing-up solutions for a mean field equation.  In Annales de l’IHP Analyse non linéaire, volume 22, pages 227–257, 2005.
  * [14] F. Gladiali, M. Grossi, H. Ohtsuka, and T. Suzuki.  Morse indices of multiple blow-up solutions to the two-dimensional Gel’fand problem.  Comm. Partial Differential Equations, 39(11):2028–2063, 2014.
  * [15] M. Grossi, H. Ohtsuka, T. Suzuki, et al.  Asymptotic non-degeneracy of the multiple blow-up solutions to the Gel’fand problem in two space dimensions.  Adv. Differential Equations, 16(1/2):145–164, 2011.
  * [16] M. Grossi and A. Pistoia.  Multiple blow-up phenomena for the sinh-poisson equation.  Archive for Rational Mechanics and Analysis, 209(1):287–320, 2013\.
  * [17] D. Gurarie and K. Chow.  Vortex arrays for sinh-Poisson equation of two-dimensional fluids: Equilibria and stability.  Physics of Fluids, 16(9):3296–3305, 2004.
  * [18] J. Jost, G. Wang, D. Ye, and C. Zhou.  The blow up analysis of solutions of the elliptic sinh-gordon equation.  Calculus of Variations and Partial Differential Equations, 31(2):263–276, 2008.
  * [19] G. Joyce and D. Montgomery.  Negative temperature states for the two-dimensional guiding-centre plasma.  Journal of Plasma Physics, 10(1):107–121, 1973.
  * [20] S. Julien, J.-M. Chomaz, and J.-C. Lasheras.  Three-dimensional stability of periodic arrays of counter-rotating vortices.  Physics of Fluids, 14(2):732–743, 2002.
  * [21] G. R. Kirchhoff.  Vorlesungen über mathematische physik: mechanik, volume 1.  Teubner, 1876.
  * [22] A. J. Majda and M. Holen.  Dissipation, topography, and statistical theories for large-scale coherent structure.  Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 50(12):1183–1234, 1997.
  * [23] R. Mallier and S. Maslowe.  A row of counter-rotating vortices.  Physics of Fluids A: Fluid Dynamics, 5(4):1074–1075, 1993.
  * [24] C. Marchioro and M. Pulvirenti.  Mathematical theory of incompressible nonviscous fluids, volume 96.  Springer Science & Business Media, 2012.
  * [25] W. Matthaeus, W. Stribling, D. Martinez, S. Oughton, and D. Montgomery.  Decaying, two-dimensional, Navier-Stokes turbulence at very long times.  Physica D: Nonlinear Phenomena, 51(1-3):531–538, 1991.
  * [26] W. H. Matthaeus, W. T. Stribling, D. Martinez, S. Oughton, and D. Montgomery.  Selective decay and coherent vortices in two-dimensional incompressible turbulence.  Physical review letters, 66(21):2731–2734, 1991.
  * [27] D. Montgomery and G. Joyce.  Statistical mechanics of “negative temperature” states.  The Physics of Fluids, 17(6):1139–1145, 1974.
  * [28] D. Montgomery, W. H. Matthaeus, W. Stribling, D. Martinez, and S. Oughton.  Relaxation in two dimensions and the “sinh-Poisson” equation.  Physics of Fluids A: Fluid Dynamics, 4(1):3–6, 1992.
  * [29] D. Montgomery, W. H. Matthaeus, W. Stribling, D. Martinez, and S. Oughton.  Relaxation in two dimensions and the “sinh-Poisson”equation.  Physics of Fluids A: Fluid Dynamics, 4(1):3–6, 1992.
  * [30] H. Ohtsuka.  To what extent can the Hamiltonian of vortices illustrate the mean field of equilibrium vortices? (Modern approach and developments to Onsager’s theory on statistical vortices).  RIMS Kôkyûroku Bessatsu, 1798:1–17, 2012.
  * [31] D. Schecter, D. Dubin, K. Fine, and C. Driscoll.  Vortex crystals from 2D Euler flow: Experiment and simulation.  Physics of Fluids, 11(4):905–914, 1999.
  * [32] F. Spineanu and M. Vlad.  Stationary vortical flows in two-dimensional plasma and in planetary atmospheres.  Physical review letters, 94(23):235003, 2005.
  * [33] J. Spruck.  The elliptic sinh gordon equation and the construction of toroidal soap bubbles.  In Calculus of variations and partial differential equations, pages 275–301. Springer, 1988.
  * [34] Z. Yin, D. Montgomery, and H. Clercx.  Alternative statistical-mechanical descriptions of decaying two-dimensional turbulence in terms of “patches” and “points”.  Physics of Fluids, 15(7):1937–1953, 2003.