A Fully Analog Pipeline for Portfolio Optimization

Authors: james s cummins
natalia g berloff
Venue: NeurIPS_cc_2024_Workshop_MLNCP
Published: 2024-09-11 23:04:58

Abstract: Portfolio optimization is a ubiquitous problem in financial mathematics that relies on accurate estimates of covariance matrices for asset returns. However, estimates of pairwise covariance are notoriously poor and calculating time-sensitive optimal portfolios is energy-intensive for digital computers. We present an energy-efficient, fast, and fully analog pipeline for solving portfolio optimization problems that overcomes these limitations. The analog paradigm leverages the fundamental principles of physics to recover accurate optimal portfolios in a two-step process. Firstly, we utilize equilibrium propagation, an analog alternative to backpropagation, to train linear autoencoder neural networks to calculate low-rank covariance matrices. Then, analog continuous Hopfield networks output the minimum variance portfolio for a given desired expected return. The entire efficient frontier may then be recovered, and an optimal portfolio selected based on risk appetite.

Main paper: 
 introduction: 
 Portfolio optimization involves creating an investment portfolio that balances risk and return. The objective is to allocate assets optimally to maximize expected returns while minimizing risk. Naturally, this problem is of great interest to financial organizations and is pivotal in risk management. However, the problem, formulated by Markowitz's mean-variance model [1] , suffers from problems in practice. Namely, it is well known that estimates of pairwise covariance between assets are notoriously poor [2] . Data samples tend to include significant amounts of noise, distorting the underlying relationships between assets. To overcome this issue, factor models were introduced that vastly reduce the dimensionality [3] . Factor methods produce low-rank covariance matrices that retain only the largest eigenvalues and discard small eigenvalues associated with noise. Despite this development, the computation of optimal portfolios remains energy-intensive as the efficient frontier is mapped out in return-variance space. In high-frequency trading, this becomes a time-sensitive computation as assets are purchased and sold on microsecond timescales [4] , and portfolios must be regularly rebalanced so as not to exceed risk appetites. Much attention has been focused on portfolio optimization in the high-frequency domain [5, 6, 7] , including the use of evolutionary algorithms to update efficient frontiers [8] . While algorithmic enhancements provide incremental gains, the exploration of alternative hardware paradigms has the potential to drive significant advancements. By using fundamental principles such as minimizing entropy, energy, and dissipation [9] , or, perhaps, incorporating quantum phenomena like superposition and entanglement [10] , we can advance and surpass classical computations of these problems. At the forefront of this drive to alternate architectures is the integration of analog, physics-based algorithms and hardware, which involve translating complex optimization problems into universal spin Hamiltonians [11, 12, 13] . Indeed, the mean-variance portfolio optimization framework can be encoded into a Hamiltonian's coupling strengths with the physical system recovering the Hamiltonian's ground state, which corresponds to the optimal portfolio solution [14, 15] . Efficient mapping from the original problem description to spin Hamiltonian enables the problem to remain manageable despite increasing complexity [16] . 
mean-variance optimization: 
 We define µ i as the expected return of asset i, and [Σ] ij = σ ij = Cov(i, j) as the covariance between assets i and j. The decision variables are w i , the proportion of the total investment in asset i. For a universe of securities with n assets, the Markowitz mean-variance portfolio optimization problem isEQUATIONfor i = 1, . . . , n, and the condition w i ≥ 0 prohibits shorting [14] . The variance w T Σw quantifies the portfolio risk for positive semidefinite matrix Σ, while R is the desired expected return of the portfolio. µ and Σ are not known a priori and must be estimated from historical data. The efficient frontier is calculated by solving (1) for various R. The efficient frontier is the set of portfolios that minimize the risk for a given R. We illustrate a frontier in Appendix A for a toy model with n = 2 assets. It was recently suggested that portfolio optimization problems could be solved on analog spatial-photonic Ising machines for equal-weighted portfolios, that is, w i ∈ {0, 1/q} for q selected assets [15] . We go beyond this constraint by utilizing analog Hopfield networks and consider the quadratic continuous optimization problem (1) . In Section 3, we aim to recover the optimal asset weights w, given known expected returns µ and covariance matrix Σ. 
continuous hopfield network: 
 A continuous Hopfield network is a type of Hopfield neural network which has continuous states and dynamics [17] . It is an analog computational network for solving optimization problems. For a network of size n, the i-th network element at time t is described by a real input x i (t), and the network dynamics are governed byEQUATIONwhere v i = g(x i ) is a nonlinear activation function, p(t) is an annealing parameter, m i are the offset biases, and J ij are elements of the symmetric coupling matrix J. Should g(x) be a non-decreasing function, then the steady states of the continuous Hopfield network (2) are the minima of the Lyapunov functionEQUATIONWe choose the functional form of g(x), such that when p(t) → 0, the minima of E occur for v i ∈ [0, 1] and correspond to the minima of -v T Jv. Therefore, by setting J = -Σ, we can minimize the variance w T Σw of problem (1). To satisfy the constraints in problem (1) we introduce Lagrange multiplier-like scalars λ 1 , λ 2 and seek to minimize the expressionH = w T Σw + λ 1 (µ T w -R) 2 + λ 2 (1 T w -1) 2 .Therefore, after discarding constants, we seek to minimizeEQUATIONwhere J = -2Σ -2λ 1 µµ T -2λ 2 11 T , and m = 2Rλ 1 µ + 2λ 2 1. Equation ( 4) can be directly encoded into the Hopfield network (2), and if required, m can be absorbed into J by introducing an additional auxiliary spin. The non-decreasing monotonic function g(x) is chosen to be the logistic function g(x) = 1/[1 -exp(-x)] to limit possible values of v i such that 0 ≤ v i ≤ 1.We illustrate the Hopfield network dynamics in Appendix B for a randomly generated covariance matrix Σ and expected return vector µ. The energy minimization properties of Hopfield networks make them particularly suitable for solving combinatorial optimization problems. Further extensions have been proposed to increase convergence to optimal states in challenging optimization problems. For example, the first-order Eq. ( 2) can be momentum-enhanced and replaced with a second-order equation leading to Microsoft's analog iterative machine [18] or Toshiba's bifurcation machine [19] . 
low-rank approximation: 
 We now focus on calculating a low-rank approximation of the covariance matrix, which will be used in (1) . If x i ∈ R n are the i-th sample of asset returns over N total samples, and we assume that E[x] = 0, then the sample covariance matrix isS = 1 N N i=1 x i x T i .When the number of samples N is of the same magnitude as n, then the sample covariance matrix usually suffers a large estimation error [2, 20] . Many low-rank factor analysis techniques exist to improve the covariance matrix estimate. Here, we consider asset returns x as random variables that follow the modelEQUATIONwhere x ∈ R n is the observed data, A ∈ R n×r is a factor loading matrix, s ∈ R r is the vector of latent variables, and e ∈ R n is uncorrelated random noise, where r ≪ n [21] . Here, s represents macroeconomic factors like the growth rate of the GDP, unemployment, inflation etc. Under the assumption that s and e are uncorrelated, the covariance matrix is thenEQUATIONEQUATIONwhere P ≡ E[ss T ] ∈ R r×r has rank(P) ≤ r, rank(A) ≤ r, and Ψ is a diagonal matrix containing the variance of noise on its diagonal. Since rank(AB) ≤ min(rank(A), rank(B)), then rank(APA T ) ≤ r. Therefore, we have decomposed the covariance matrix Σ into a positive semidefinite low-rank matrix plus a positive semidefinite diagonal matrix. Defining M ≡ APA T , low-rank factor analysis concerns the estimation of M and Ψ. To calculate M and Ψ we solve the minimization problem minEQUATIONwhere || • || F denotes the Frobenius norm [22] . We present a common digital computing method in Appendix (C) for solving problem (8) based on principal component analysis (PCA). The eigendecomposition in PCA becomes computationally expensive as the data size grows. Alternatively, autoencoders -particularly when implemented using stochastic gradient descent -can handle larger datasets and higher-dimensional data more efficiently than PCA [23] . Additionally, when integrating dimensionality reduction as part of a larger neural network framework, an autoencoder can be easily embedded within the pipeline, whereas PCA would need to be applied as a separate pre-processing step [24] . 
linear autoencoders: 
 A linear autoencoder is a classic neural network model for unsupervised learning that is trained to learn the identity function. The input and output layers have the same number of nodes, while the middle layer has fewer nodes. It aims to approximate the input through learning linear encodings and decodings between input and latent space. The encoderB ∈ R r×n maps input X = [x 1 , . . . , x N ] ∈ R n×N into a low-dimensional latent space [s 1 , . . . , s N ],and the decoder A ∈ R n×r maps [s 1 , . . . , s N ] back to the original representation X. We therefore recover the same model as in Eq. ( 5), and training the linear autoencoder becomes the minimization problem [25] minEQUATIONWe do not explicitly express the learnable biases in the network as these may be absorbed into the encoder B and decoder A by introducing an auxiliary row into X that is permanently clamped to values of 1. We illustrate the training of a linear autoencoder in Fig.(1)(a)-(c) with the backpropagation method and compare it to PCA in Fig.(A linear autoencoder is related to PCA. Indeed, under mild nondegeneracy conditions, any A at a local minimizer recovers the top rank-r eigenspace of XX T [26] . However, unlike actual PCA, the coordinates of the output of the middle layer in the network are correlated and are not sorted in descending order of variance [27] . Autoencoder neural networks typically use backpropagation to train the weights. However, backpropagation is energy-intensive and not biologically plausible.  
equilibrium propagation: 
 On dedicated analog hardware, equilibrium propagation is an energy-efficient alternative to backpropagation [28] . Therefore, in the supervised learning setting studied here, it may be used to train the weights of a linear autoencoder. Equilibrium propagation is an energy-based model because it relies on the concept of energy minimization to learn and make predictions. We consider the continuous Lyapunov function (3) with p(t) = 1 and m i = 0 for all i, where here the symmetric coupling weights J ij are to be learned, and nonlinear activation function g(x) need not be the same as in Section 3. Neurons x i are split in three sets: the input neurons, which are always clamped, the hidden neurons, and the output neurons. The discrepancy between the desired output y and the realized output x is measured by the cost functionEQUATIONwhich forms part of the total energy function F = E + βC. The clamping factor β ≥ 0 is a real-valued scalar that allows the output neurons to be weakly clamped [29] . The continuous-time dynamical system evolves according to the differential equation of motionEQUATIONwhich is formed of two parts. The first is the internal force induced by the internal Hopfield energy, given by Eq. ( 2) for all i, and the second, the external force, is induced by the cost function C asEQUATIONfor nodes in the output layer Y. Equilibrium propagation has two modes: the free phase and the weakly clamped phase. In the free phase β = 0 and only the inputs are clamped. The network then converges to a fixed point x * and the output units are read out. In the weakly clamped phase β > 0, which induces an external force that acts on the output units as in Eq. ( 12). This force nudges the outputs from their fixed point values in the direction of the target values y i . This perturbation propagates among the hidden neurons before a new fixed point x * β is found. Then, another weakly clamped phase is executed, this time with β → -β leading to the weakly clamped equilibrium x * -β . It was shown that the weakly clamped phase implements the propagation of error derivatives with respect to the synaptic weights [29] . In the limit β → 0, the update rule isEQUATIONwhich is a second-order approximation to the standard backpropagation derivative [30] . The process is iterated, at each step updating the weights J ij to minimize the loss function C. We choose activation functiong(x) = x if |x| ≤ c c • sgn(x)otherwise, (14) with constant c, so that under the condition that |x i | ≤ c for all i, Eq. ( 2) is linear and can thus represent a linear autoencoder. In this case, the output x of Eq. ( 2) is then the solution to the linear differential equation dx/dt = (J -I)x, and thereforeEQUATION)The constant c in Eq. ( 14) is chosen to be large enough such that after training, all neurons obey |x i | ≤ c, and we can associate the Hopfield network as a linear autoencoder. To achieve a steady state in Eq. ( 15), at least one eigenvalue of J -I should be zero, with all others having a negative real part.Proposition. We state, with proof given in Ref. [26] , that for any fixed n × r matrix A, Eq. ( 9) attains its minimum for B = (A T A) -1 A T .Lemma. The n × n matrix J -I, where J = AB, has at least one zero eigenvalue, with all others having negative real part.Proof. J = AB = A(A T A) -1 A T , and thereforeEQUATIONwhich shows that J is idempotent, that is J 2 = J. It follows that J is a projection operator on the column space C(J) along its null space N (J). The n eigenvalues λ i of J are either 0 or 1:λ i x i = Jx i = J 2 x i = λ i Jx i = λ 2 i x i , which implies λ i ∈ {0, 1}. By construction, J has rank at most r, and therefore there are at least n -r zero eigenvalues. It follows that there are between 1 and r nonzero eigenvalues of J, which must have value λ i = 1. Since J -I has eigenvalues µ i = λ i -1, then µ i ∈ {-1, 0}. Therefore, J -I has between 1 and r zero eigenvalues, with all others being equal to -1.The Lemma guarantees that should equilibrium propagation learn the weights that minimize Eq. ( 9), the corresponding Hopfield network will converge to a steady state. Yet, during training, this will, in general, not be the case, and positive eigenvalues of J -I will produce exponential growth in Eq. (15) . However, Eq. ( 15) only holds in the linear regime of the activation function (14) . Exponential growth is prevented by the symmetric clipping incorporated into the nonlinear activation function g(x) for neurons with |x i | > c.In the linear regime, the overall network dynamics is represented by the square matrix lim t→∞ exp{(J -I)t}, which for linear autoencoders we seek to decompose into its non-square constituent parts: encoder B and decoder A. We achieve this by treating the encoder and decoder as separate Hopfield networks, as shown in Fig . (1)(g), each with their own energy function. The encoder settles into an equilibrium representing the latent vector s without taking into account the decoder. s is then used as a fixed input to the decoder which then settles into its own equilibrium. The decoder then undergoes the weakly clamped phases, and its weights are updated according to Eq. ( 13). The encoder weights also need to be optimized to lower the reconstruction loss at the decoder output, which is achieved by settingEQUATIONin Eq. (11) , where x (dec) βis the weakly clamped decoder equilibrium state, and Eq. ( 18) only pertains to neurons in the encoder output layer Y. Equation ( 18) follows from the fact that it can be shown that equilibrium propagation also allows for finding the gradient of the loss with respect to the input [31] . We note that J, which contains the couplings of the continuous Hopfield network, is now a (n + r) × (n + r) matrix on account of the number of nodes in the encoder and decoder networks. Nonetheless, the factor loading matrix A can be recovered as the n × r block corresponding to the nodes of the decoder output layer. The equilibrium propagation training procedure is illustrated in  
results: 
 We collect real data samples x i ∈ R n from stock returns of a selection of n = 100 stocks in the S&P 500 index. We restrict ourselves to only N = 100 observations such that the sample covariance matrix has a tendency to contain significant noise. Two continuous Hopfield networks, structured as the encoder and decoder parts of a linear autoencoder, are trained using equilibrium propagation. The latent variables [s 1 , . . . , s N ] are calculated as the subset Y of steady-state solutions of the encoder network, while the factor loading matrix A is the n × r block of the decoder matrix representation lim t→∞ exp{(J -I)t} corresponding to its output layer Y. In practice, we cannot take the limit to infinity, and instead, we use a suitably large value of t such that exp{(J -I)t} changes minimally from t to t + 1. We depict the full-rank sample covariance matrix and the equilibrium propagation-based low-rank approximation in Figs. Analog Hopfield networks can be implemented as electronic circuits [32] and photonic neural networks [33] . Photonic systems operate on picosecond to femtosecond timescales as high bandwidth signals flow through a single optical waveguide. Consequently, such implementations can have dense connectivity while maintaining fast convergence times. However, physical analog platforms are subject to noise sensitivity, thermal effects, and non-idealities in circuit components which can degrade performance. In addition, real-world portfolio optimization problems often involve complex constraints such as transaction costs, market liquidity, regulatory requirements, and cardinality constraints. While some of these can be readily incorporated into the objective function (4), for example, an ℓ 1 -norm can enforce sparsity to satisfy a cardinality constraint, others take more complex forms. To address the limitations, a hybrid approach that combines analog Hopfield networks with digital computing could be explored.The efficient frontier represents the set of optimal portfolios that offer the highest expected return for a given level of risk, or the lowest risk for a given return. It is derived by plotting the risk-return profiles of various portfolios, with the frontier itself being the curve where no portfolio exists with both a higher return and lower risk. For the toy model in Fig. (3) with n = 2 assets, the entire set of possible portfolios lies on a single hyperbola. The Sharpe ratio S r for each portfolio in the efficient frontier. We now see that the purple circle is indeed the portfolio that maximizes the Sharpe ratio.The quadratic continuous optimization problem (4) is solved using continuous Hopfield network (2) . The trajectories of the neurons x i (t) are illustrated in Fig. ConclusionsThis paper introduces a fully analog pipeline for portfolio optimization problems. Starting with raw data samples, the proposed pipeline leverages the energy-efficient analog operation of continuous Hopfield networks to calculate optimal portfolio weights. The analog pipeline distinguishes itself from traditional digital methods by its speed and scalability, with applications in time-sensitive domains such as high-frequency trading. At the heart of the pipeline are continuous Hopfield networks, used in two separate applications: autoencoder neural networks and minimum variance portfolios. By shifting to analog architectures, we reduce the reliance on binary logic operations typical of digital systems, paving the way for a more energy-efficient approach to computation. This efficiency can reduce power consumption in data centers and other computing environments, addressing the growing energy demands of digital computing. Specifically, companies can reduce their energy consumption while optimizing large portfolios as part of their risk management processes. 
 
 Appendix:  



Reviews: review: title: title: A use-case of continuous Hopfield networks and equilibrium propagation for portfolio optimization
main_review: main_review: This paper is about solving the Markowitz mean-variance portfolio optimization problem using algorithms compatible (in principle) with analog hardware, specifically: continuous Hopfield networks, and equilibrium propagation.

The portfolio optimization problem is formulated as a constrained optimization problem: minimizing the risk for a given level of expected returns. The challenge is that the covariance matrix and expected returns are not explicitly known ; these must be estimated using samples from historical data. The authors estimate the covariance matrix using a low-rank approximation method, which is known to yield better results than naive estimation methods. Specifically, they use an autoencoder, consisting of a linear encoder and a linear decoder. The encoder and decoder take the form of Continuous Hopfield Networks, trained using Equilibrium Propagation.

Once the encoder and decoder are trained, a new continuous Hopfield network is built, representing the low-rank approximation of the covariance matrix, and whose ground state is an approximate solution of the original portfolio optimization problem.



**Presentation**

In Section 2, it could be useful to remind the reader that one of the challenges of Problem (1) is that the vector of expected returns mu and the covariance matrix Sigma are not known. We must estimate their values using samples from historical data.

The analytical expression of J in Eq.3 is (re)defined several times: line 60, line 64, and again line 151, which makes it a little difficult to follow.
I am wondering if Section 3 (Continuous Hopfield Networks) would be better placed between Section 5 (Linear Autoencoders) and Section 6 (Equilibrium Propagation), in which case you might consider removing the first two definitions of J (line 60 and line 64).
Alternatively, it could be clarified that Sections 2 and 3 are about “finding the optimal w given Sigma and mu”, whereas Sections 4, 5, 6 are about “finding good approximates of Sigma and mu”. (see also my third question below)



**Questions**

The original constrained optimization problem (Eq 1) is approximated by an unconstrained optimization problem using a soft penalty on the constraints (Eq 4). How does this affect the solution?

Similarly, in the unconstrained optimization problem, outputs are normalized during postprocessing. It is not directly clear if the v (asset allocation) obtained this way is a solution to the original problem. How does this method affect the solution?

If I understood correctly, there are not 2, but 3 different continuous Hopfield networks: one for the encoder, one for the decoder, and a third one whose coupling strengths correspond to the covariance matrix of the original portfolio optimization problem (which is used for inference, once the encoder and decoder have been trained). Is this correct?



**Minor remarks**

For the function E (Eq.3) to be a Lyapunov function of the dynamics (Eq.2), I believe that the nonlinear activation function g must be increasing. It is mentioned at line 66 that g is non-decreasing, but it could be explained earlier (at line 56) that this is a requirement for the E to be a Lyapunov function.

If I understood correctly, the method presented in Appendix C is a common method to solve the problem on a digital computer, but is not applicable on analog systems. If so, this could be clarified.

Ref [17] should probably be replaced by:
Hopfield, John J. "Neurons with graded response have collective computational properties like those of two-state neurons." Proceedings of the national academy of sciences 81.10 (1984): 3088-3092.
title: title: An interesting use case of Equilibrium Propagation, however very difficult to read
main_review: main_review: **Summary.** The paper proposes to apply Equilibrium Propagation (EP), an algorithm suited for analog hardware, to portfolio optimization. Namely, given $n$ assets $X_k$ and the covariance matrix $\Sigma:=\mathbb{E}[X \cdot X^\top]$ associated to these assets and $w_i$ the proportion of asset $i$ inside the portfolio, the optimal asset allocation should minimize $E:=w^\top \cdot \Sigma \cdot w$ under positivity and convexity constraint on $w$ and ensuring a return $R$ (assuming some value for the assets given by another vector $\mu$). The subtlety though is that the empirical covariance matrix $\Sigma$ may not be used itself, but a *low-rank* approximation thereof. This can be done by using a *linear* autoencoder $A\cdot B$ where $A$ and $B$ have to be learned such that $\ell := || X - A \cdot B \cdot X ||_F^2$ is minimized. Therefore, the goal is to find $A, B$ which minimize $\ell$ *under the constraint that* the energy $E$ is minimized. The approach is demonstrated on 100 observations of 100 stocks taken from the S\&P 500 index.

**Strengths**. It is an interesting use case for analog computing that is worth spreading at the workshop.

**Weakenesses**.

- I still have a hard time understanding well how the problem is framed how the application of EP reads here. Eq. 1 is very misleading because it wrongly conveys that  $E:=w^\top \cdot \Sigma \cdot w$ is the outer objective, but it actually turns out (if I understood correctly?) to be the *inner* objective, i.e. the constraint to be satisfied *while minimizing the reconstruction error* $\ell := || X - A \cdot B \cdot X ||_F^2$. There are lots of different notations that are introduced for the free variables, for the various parameters such that I really fail to understand eventually what are precisely the variables of the *outer* optimization problem ($\ell$ as defined above) and **how the energy function $E$ is precisely defined in terms of the variables appearing inside** $\ell$.

- "it can be shown that equilibrium propagation also allows for finding the gradient of the loss with respect to the input" should be attributed to [1], which precedes [28] by two years. 

[1] Scellier, B. (2021). A deep learning theory for neural networks grounded in physics. arXiv preprint arXiv:2103.09985
[28] Tom Van Der Meersch, Johannes Deleu, and Thomas Demeester. Training a hopfield variational autoencoder with equilibrium propagation. In Associative Memory & Hopfield Networks in 2023, 2023.