Belief Propagation Algorithm for Portfolio Optimization Problems

The typical behavior of optimal solutions to portfolio optimization problems with absolute deviation and expected shortfall models using replica analysis was pioneeringly estimated by S. Ciliberti et al. [Eur. Phys. B. 57, 175 (2007)]; however, they have not yet developed an approximate derivation method for finding the optimal portfolio with respect to a given return set. In this study, an approximation algorithm based on belief propagation for the portfolio optimization problem is presented using the Bethe free energy formalism, and the consistency of the numerical experimental results of the proposed algorithm with those of replica analysis is confirmed. Furthermore, the conjecture of H. Konno and H. Yamazaki, that the optimal solutions with the absolute deviation model and with the mean-variance model have the same typical behavior, is verified using replica analysis and the belief propagation algorithm.


Introduction
Portfolio optimization is one of the most fundamental frameworks of risk diversification management. Its theory was introduced by Markowitz in 1959 and is one of the most important areas being actively investigated in financial engineering [1][2][3]. In their theoretical research, Ciliberti and Mézard assessed the typical behavior of optimal solutions to portfolio optimization problems, in particular those described by the absolute deviation and expected shortfall models, using replica analysis, one of the most powerful approaches in disordered systems. With this approach, they showed that the phase transitions of these optimal solutions were nontrivial [1]. However, they did not develop an effective algorithm for finding the optimal portfolio with respect to a fixed return set. This requires a rapid algorithm for resolving the optimal portfolio problem with respect to a large enough in-sample set.
As a first step in such a research direction, we propose an algorithm based on belief propagation, which is well-known as one of the most prominent algorithms in probabilistic inference, to resolve the microscopic averages of the optimal solution in a feasible amount of time for a fixed return set. We also confirm whether the numerical experimental results of our novel algorithm are consistent with the ones of replica analysis. Furthermore, the conjecture of Konno and Yamazaki, that if the return at each period is independently and identically drawn from the normal probability distribution [2], the optimal portfolio of the mean-variance model is consistent with that of the absolute deviation model, is supported using replica analysis and belief propagation.

Model Setting and Proposed Algorithm
Let us define the model setting for our discussion. A portfolio of N assets and the return at period μ are represented by w = {w 1 , w 2 , Á Á Á, w N } T 2 R N and x μ = {x 1μ , x 2μ , Á Á Á, x Nμ } T 2 R N , respectively, where w k is the position of asset k, and we assume for simplicity that the mean of the return of asset k in period μ, x kμ , is zero. The notation T indicates matrix transposition. Given a return set for p periods as reference, the problem is to minimize the following cost function (i.e., Hamiltonian) for the portfolio: where R(u) represents a cost function, such as u 2 2 in the mean-variance model and juj in the absolute deviation model, respectively. Furthermore, since the budget is assumed to be finite, the following global constraint is set: One of our aims is to develop an effective general algorithm for solving this problem; in particular, our aim is an algorithm that works for all cost functions R(u) and all probability distributions of the returns.
As a basis for the proposed algorithm, following examples in statistical mechanics, we set the joint probability of portfolio w used in Eq (1) using finite inverse absolute temperature β as follows: where g(u) = e −βR(u) is the likelihood function and prior probability P 0 ðwÞ / expm P N k¼1 w k À N À Á Â Ã for sufficiently large N [4,5]. Notice that the partition function of this posterior probability is implicitly ignored in this analysis because intuitively it is possible to evaluate the firstand second-order moments of portfolio w k approximately without the partition function by the following procedure. An arbitrary test probability of portfolio is defined as follows: where the reducibility condition on beliefs b k (w k ) and b μ (w), must hold and w\w k denotes a subset of w from which w k is excluded. The Kullback-Liebler divergence (KLD) provides a useful guideline for deriving the belief propagation algorithm. However, since it is too complicated to directly assess KLD except in specific graphical models, we here approximate the Bethe free energy denoted as follows: where P 0k ðw k Þ / em w k is used. The purpose of this step is to derive the optimal portfolio using the beliefs b k (w k ) and b μ (w) that minimize the Bethe free energy under the reducibility condition of Eq (6). By adding the term (8), it is possible to minimize the Bethe free energy with respect to the beliefs to obtain Furthermore, for simplicity, we set as novel auxiliary functions, and then b k (w k ) and b μ (w) can be rewritten using 1 Moreover, applying the cumulant generating functions the first and second moments of w k have the compact forms m wk ¼ @ k ðy k Þ @y k ¼ @ m ðyÞ @y k and w wk ¼ This allows us to disregard the calculation of the partition function. Then, our proposed algorithm for sufficiently large N comprises the following: x km m um þw wk m wk ; ð17Þ x km m wk Àw um m um ; ð21Þ where [4][5][6][7]. In addition,w wk m wk andw um m um describe the Onsager reaction terms in the literature of spin glass theory (respectively [8,9]).
Four points should be noticed here. First, the calculation of this procedure is reduced from O(N 3 ) to O(N 2 ). For instance, in the case of the mean-variance model, although we are required to calculate the inverse matrix of the correlation matrix of return set XX T 2 M N×N , where return matrix X = {x 1 , Á Á Á, x p } 2 M N×p , in order to assess the optimal solution rigorously, it is well-known that this calculation is O(N 3 ). Moreover, fortunately it is found that in the case of the mean-variance model, this algorithm derives the exact optimal solution (see appendix A for details). Second, only Eqs (20) and (23) are dependent on the likelihood function g(u) = e −βR(u) , and the variables of index u are the only model dependent ones. Furthermore,m is determined by Eqs (2) and (16). Third, the randomness of return is not assumed to be sampled from specific distributions. Because it is plausible that the assumption on the Bethe free energy approximation works correctly if the return at each period is asymptotically not correlated with other returns. Lastly, we expect that in the limit as β ! 1, the estimate of the portfolio of asset k, m wk , asymptotically corresponds to the optimal portfolio with respect to the given return set.

Numerical Experimental Results
In order to confirm the effectiveness of our method, the numerical experimental results of the proposed algorithm and those of the replica analysis for the case of the Markowitz model are shown in Figs 1 and 2, where x kμ are independently and identically drawn from the normal distribution with mean and variance 0 and 1, respectively. The numerical experimental result of belief propagation is assessed from 10 2 samples of the number of assets N = 500 and is denoted by error bars and the result of replica analysis is denoted by a solid line. Both findings indicate that the two approaches are consistent with each other. With regard to the conjecture of Konno and Yamazaki, the variables in Eqs (20) and (23), in the case of the mean-variance model and the absolute deviation model w um ¼ À @m um @h um ; ð27Þ ffiffiffiffiffi ffi 2p p u À Á À1 e À u 2 2 in the case of u ) 1, m um ' À h um w um and w um ' 1 w um are estimated; that is, this finding indicates that the conjecture of Konno and Yamazaki is valid in part in the sense of the belief propagation approach. See appendices for details.

Summary
In conclusion, we have discussed an effective algorithm for finding the optimal solution of the portfolio optimization problem with respect to an arbitrary cost function according to Ciliberti and Mézard [1]. With loss of generality, applying the likelihood function g(u) defined by the cost function R(u) dependent on the risk diversification problem, we proposed a novel approximation derivation method based on one of the most powerful estimation methods in probabilistic inference. In addition, since two types of Onsager reaction terms are derived in Eqs (17) and (21), our algorithm provides the Thouless, Anderson, and Palmer approach rather than the mean-field approximation in the literature of spin glass theory. One advantage of our algorithm is that it rapidly converges by excluding the effect of self-response. In order to confirm the effectiveness of the proposed approach, we have described the case of the mean-variance model. Furthermore, we have shown that the conjecture of Konno and Yamazaki is supported by employing both approaches developed in cross-disciplinary research involving statistical mechanics and information sciences. In future work, we will assess the properties of R(u) and the randomness of return that make solving the portfolio optimization problem using belief propagation possible.
We here confirm the exactness of the proposed belief propagation algorithm for the case of the Markowitz model. Our discussion is restricted to α > 1 for simplicity. From Eqs (21), (24), and (25), we obtain m u ¼ À b ffiffi ffi

B. Replica Analysis
According to Ciliberti and Mézard, replica symmetry solution of the portfolio optimization problem, where x kμ is independently and identically distributed with N(0,1), is derived as the following extremum: where Z ¼ P w P 0 ðwÞ is the partition function and the notation [Á Á Á] q denotes the quenched average over the return set. Moreover, the quenched overlap parameters become and q otherwise by employing replica indices a, b = 1, 2, Á Á Á, n and the assumption of replica symmetry. Furthermore, for large N and p, α = p/ N * O(1) remains finite and plays an important role as a control parameter with respect to phase transition phenomena. If gðuÞ ¼ e À b 2 u 2 , then q ¼ 1 À 1 a À Á À1 and χ = (β(α − 1)) −1 can be exactly calculated in the case α > 1 and q ! 1, and χ ! 1 otherwise. This analytical finding is also verified in by the following. It is well known that the eigenvalue distribution of the correlation matrix C = 1 N XX T in the limit of N ! 1 is asymptotically close to the Marčhenko-Pastur law rðlÞ ¼ ½1 À a þ dðlÞ þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ½lÀl À þ ½l þ Àl þ p 2pl with l AE ¼ 1 AE ffiffi ffi a p ð Þ 2 and [u] + = max{u,0} [10].
Therefore, q ¼ h 1 l 2 ih 1 l i À2 and one degree of the cost function ε ¼ lim N!1 3 if α > 1 and approach infinity otherwise follows directly. This is consistent with the findings of replica analysis [5].
In general, the order parameters are derived as follows: From Eqs (29) and (30), is obtained. In the limit of sufficiently large β of g(u) = e −βjuj , if we assess Z ' À ffi ffi q p w and d ' q w 2 asymptotically, then the conjecture of Konno and Yamazaki is confirmed as correct in the sense of replica analysis.

C. The Conjecture of Konno and Yamazaki
This conjecture is related to the assessment of an annealed system in the context of spin glass theory. If the return at period μ, x μ , is independently and identically drawn from N(0, S), where S 2 M N×N is variance-covariance matrix and w is fixed, the novel variable z ¼ is distributed as N(0, s 2 (w)) with s 2 ðwÞ ¼ 1 N w T P w 2 R. With respect to fixed w, employing one degree of the cost function of the annealed optimization problem DuR usðwÞ ð Þ, which becomes ε MV ðwÞ ¼ a 2 s 2 ðwÞ in the case of the mean-variance model and ε AD ðwÞ ¼ 2a ffiffiffi ffi 2p p j sðwÞ j in the absolute deviation model. This implies that the optimal portfolios of the annealed situations of the two models are consistent with each other. Note that one degree of the cost function in the case of the annealed portfolio problem with the expected shortfall model, ε ES ðwÞ ¼ min v!0 a vg þ H v sðwÞ n o with γ > 0 can also be assessed. If sðwÞ 1 ffiffiffi ffi 2p p g , then this optimal solution is identical to those of the previous mentioned models. This finding, that is, arg min w T e¼N ε MV ðwÞ ¼ arg min w T e¼N ε AD ðwÞ, is one part of the contributions reported by Konno and Yamazaki. However, they optimistically assumed w MV = w AD with respect to a given return set X without any mathematical proof, using w k x km : ð35Þ As explained above, arg min w T e¼N ε MV ðwÞ ¼ arg min w T e¼N ε AD ðwÞ with respect to the annealed optimization problem strictly holds; however, w MV = w AD is not always satisfied. For example, in the simple case of N = p = 2 for the two returns x 1 = {a, c} T and x 2 = {b, d} T , their assumption w MV = w AD does not hold, except under specific special situations.
Although this is apparently contradictory to these obtained findings from both approaches, it is necessary to recognize that the relation w MV = w AD with a fixed return set is equivalent to the sufficient condition q MV = q AD , where q MV ¼ lim N!1 1 N w T MV w MV Â Ã q and q AD ¼ lim N!1 1 N ½w T AD w AD q are quenched averages of overlap parameters. Moreover, although w MV = w AD does not hold in general, it is expected that the inner product w T MV w AD jw MV jjw AD j is approximately 1 because 1 N P m>v x km x jv ! 0 in the case of sufficiently large N.
Supporting Information S1 Data. Figure 1 and figure 2 from S1 Data. S1 Data is the result from numerical experiment of S1 File and S2 File. (TXT) S1 File. S1 File is main program of C language.
(TXT) S2 File. S2 File is sub program of C language. (TXT)