# Automated Design using Neural Networks and Gradient Descent

Oliver Hennigh
Mexico
Guanajuato, Gto, Mexico
<EMAIL_ADDRESS>

###### Abstract

We propose a novel method that makes use of deep neural networks and gradient
decent to perform automated design on complex real world engineering tasks.
Our approach works by training a neural network to mimic the fitness function
of a design optimization task and then, using the differential nature of the
neural network, perform gradient decent to maximize the fitness. We
demonstrate this methods effectiveness by designing an optimized heat sink and
both 2D and 3D airfoils that maximize the lift drag ratio under steady state
flow conditions. We highlight that our method has two distinct benefits over
other automated design approaches. First, evaluating the neural networks
prediction of fitness can be orders of magnitude faster then simulating the
system of interest. Second, using gradient decent allows the design space to
be searched much more efficiently then other gradient free methods. These two
strengths work together to overcome some of the current shortcomings of
automated design.

##  1 Introduction

Automated Design is the process by which an object is designed by a computer
to meet or maximize some measurable objective. This is typically performed by
modeling the system and then exploring the space of designs to maximize some
desired property whether that be an automotive car styling with low drag or
power and cost efficient magnetic bearings (Ando et al., 2010) (Dyck &
Lowther, 1996) . A notable historic example of this is the 2006 NASA ST5
spacecraft antenna designed by an evolutionary algorithm to create the best
radiation pattern (Hornby et al., ). More recently, an extremely compact
broadband on-chip wavelength demultiplexer was design to split electromagnetic
waves with different frequencies (Piggott et al., 2015). While there have been
some significant successes in this field the dream of true automated is still
far from realized. The main challenges present are heavy computational
requirements for accurately modeling the physical system under investigation
and often exponentially large search spaces. These two problems negatively
complement each other making the computation requirements intractable for even
simple problems.

Our approach works to solve the current problems of automated design in two
ways. First, we learn a computationally efficient representation of the
physical system on a neural network. This trained network can be used to
evaluate the quality or fitness of the design several orders of magnitude
faster. Second, we use the differentiable nature of the trained network to get
a gradient on the parameter space when performing optimization. This allows
significantly more efficient optimization requiring far fewer iterations then
other gradient free methods such as genetic algorithms or simulated annealing.
These two strengths of our method overcome the present difficulties with
automated design and greatly accelerate optimization.

The first problem tackled in this work is designing a simple heat sink to
maximize the cooling of a heat source. The setup of our simulation is meant to
mimic the conditions seen with an aluminum heat sink on a computer processor.
We keep this optimization problem relatively simple and use this only as a
first test and introduction to the method. Our second test is on the
significantly more difficult task of designing both 2D and 3D airfoils with
high lift drag ratios under steady state flow conditions. This problem is of
tremendous importance in many engineering areas such as aeronautical,
aerospace and automotive engineering. Because this is a particularly
challenging problem and often times unintuitive for designers, there has been
considerable work using automated design to produce optimized designs. We
center much of the discussion in this paper around this problem because of its
difficulty and view this as a true test our method. While we only look at
these two problems in this work, we emphasize that the ideas behind our method
are applicable to a wide variety of automated design problems and present the
method with this in mind.

As we will go into more detail in later sections, in order to perform our
airfoil optimization we need a network that predicts the steady state flow
from an objects geometry. This problem has previously been tackled in Guo et
al. (2016) where they use a relatively simple network architecture. We found
that better perform could be obtained using some of the modern network
architecture developments and so, in addition to presenting our novel method
of design optimization, we also present this superior network for predicting
steady state fluid flow with a neural network.

##  2 Related Work

Because this work is somewhat multidisciplinary, we give background
information on the different areas. In particular, we provide a brief
discussion of other work related to emulating physics simulations with neural
networks as this is of key importance in our method. We also review some of
the prior work in automated design of airfoils because this is the main
problem used to test our method.

##  3 Speeding up Computational Physics with Neural Networks

In recent years, there has been incredible interest in applications of neural
networks to computational physics problems. One of the main pursuits being to
emulate the desired physics for less computation then the physics simulation.
Examples of this range from simulating 3D high energy particle showers seen in
Paganini et al. (2017) to solving the Schrödinger equation seen in Mills et
al. (2017). Computational Fluid Dynamics has gotten the most attention in this
regard because of its many uses in engineering as well as computer animation
(Tompson et al., 2016) (Hennigh, 2017). The prior work that is most related to
our own is Guo et al. (2016) where they train a neural network to predict the
steady state fluid flow from an objects geometry. Our method builds on this
idea and we use the same general approach for approximating the fluid flow but
with an improved architecture.

###  3.1 Automated Design Optimization of Airfoils

To date, there has been substantial work in automated aerodynamic design for
use in aeronautical and automotive applications (Ando et al., 2010) (Anderson
& Aftosmis, 2015). Airfoil optimization in particular has received a lot of
attention where the general methodology is to refine an airfoil geometry to
minimize drag (Drela, 1998) (Koziel & Leifsson, 2013). Roughly speaking, there
are two classes of optimization strategies used here. The first class being
gradient free methods like simulated annealing, genetic algorithms, and
particle swarm methods. A look at these methods and there applications to
airfoil optimization can be found in Mukesh et al. (2012). The other class
being gradient based methods such as steepest descent. Typically gradient
based methods can perform optimization in fewer steps then gradient free
methods however computing the gradient is often very costly. The simplest
approach in doing so is finite difference method however this requires
simulating the system a proportional number of times to the dimension of the
search space in order to approximate the gradient. This is infeasible if the
fluid simulation is computationally expensive and the search space is large.
Our approach can be viewed as a gradient based method but where the gradients
are coming from a neural network that is emulating the simulation.

In order to perform automated design of airfoils one needs to parameterize the
space of possible geometries. There are a variety of approaches in doing this
and a thorough list can be found in Salunke et al. (2014). In this work we use
the parameterization technique found in Lane & Marshall (2009) and Hilton
(2007) where the upper and lower surface are described by a polynomial and the
parameters are the coefficients of this polynomial.

##  4 Gradient Decent on Parameter Space

Figure 1: Illustration of Proposed Gradient Descent Method

An automated design optimization problem can be viewed in concrete terms as
maximizing some desired fitness function $F(x)$, where
$F:X\rightarrow\mathbb{R}$ for some space $X$ of design parameters.

$\max_{\forall x\in X}F(x)$ (1)

In most real world setting, evaluating the fitness function $F$ can be
computationally demanding as is the case with our fluid simulations. The first
aspect of our method is to replace $F$ with a computationally efficient neural
network $F_{net}$. This can offer considerable speed improvements as we will
discuss bellow. The second piece of our method is the observation that
$F_{net}$ is differentiable and can be used to obtain a usable gradient in the
direction of maximizing fitness. This is in contrast to $F$ where it may be
computationally infeasible to calculate the gradient and thus require other
search techniques such as simulated annealing or genetic algorithms. Using
this gradient allows faster optimization to be performed with fewer iterations
as we will demonstrate bellow. There are some details that need to be
addressed and to do so we go through the example problem of optimizing the fin
heights on a heat sink.

In our heat sink problem, $X$ contains 15 real valued parameters between 0 and
1. Each of these parameters correspond to the height of an aluminum fin on the
heat sink as seen in the figure 2. We also fix the amount of aluminum and
scale the total length of all the fins to meet this requirement. This presents
an interesting problem of determining the optimal length each fin should have
to maximize the cooling of the heat source. The simplest application of our
method is to use the 15 fin heights as inputs to a neural network that outputs
a single value corresponding to the temperature at the heat source. This
approach has the draw back that if you want to add another constraint to the
optimization like making the left side cooler then the right side you would
need to retrain the network. A solution to this problem is to have the network
again take in the fin parameters but output the full heat distribution of the
heat sink. This allows different quantities to be optimized but is still
limiting in that our network only runs on a single parameter setup. Our
solution to this problem is to train two networks. The first network,
$P^{heat}_{net}$, takes in the fin parameters and generates a binary image
corresponding to the geometry of the heat sink. We refer to this as the
parameterization network. The second network, $S^{heat}_{net}$, predicts the
steady state heat distribution from the geometry. Because the parameterization
network is performing an extremely simple task and training data can be
generating cheaply, we can quickly retrain $P^{heat}_{net}$ if we want to
change the parameter space. The network $S^{heat}_{net}$ is now learning the
more general task of predicting steady state heat flow on an arbitrary
geometry. The same approach is used for the steady state flow problem and an
illustration depicting this can be found in figure 1 . This approach allows
our network to be as versatile as possible while still allowing it to used on
many design optimization tasks.

Up until now we have not discussed how to generate the data needed to train
these neural networks. Generating the data to train the parameterization
network is relatively simple. If the parameterization is known, we simply make
a set of parameter vectors and their corresponding geometries. In the case of
the heat sink, this is a set of examples composed of the 15 parameters and
there corresponding binary representation of the heat sink. Putting together a
dataset for $S^{heat}_{net}$ or $S^{flow}_{net}$ (fluid flow network) is
somewhat more complex. The simplest approach and the one used in this work is
to simulate the respective physics on objects drawn from the object design
space. For the heat sink problem this would entail a dataset of object
geometries and their corresponding steady state heat distributions. This
method has the disadvantage that the network only sees examples from the
current parameter search space and if it is changed the network may not be
able to accurately predict the physics. We argue this is not a significant
issue for two reasons. First, neural networks are very good at generalizing to
examples outside their train set. An example of this can be seen in Guo et al.
(2016) where the network is able to produce accurate fluid predictions on
vehicle cross sections even though it was only trained on simple polygons.
Second, it is easy to imagine a hybrid system where a network is trained on a
large set of diverse simulations and then fine tuned on the current desired
parameter space. For these reasons we feel that this approach of generating
simulation data is not significantly limiting and does not detract from the
generalizability of the approach.

###  4.1 Flow Prediction Network

In order for our method to work effectively we need a network to predict the
pressure and velocity field of the steady state flow from an objects geometry.
This is a difficult task because each point of flow is dependent on the
entirety of the geometry. This global information requirement is met in the
previous work (Guo et al., 2016) with a fully connected layer. This has
drawbacks because fully connected layers are often slow, difficult to train,
and parameter heavy. Our improved method keeps the entire network
convolutional and employs a U-network architecture seen in Ronneberger et al.
(2015) with gated residual blocks seen in Salimans et al. (2017). By making
the network deep and using many downsamples and upsamples we can provide
global information about the boundary when predicting each point of flow.
Keeping the network all convolutional also allows the spacial information to
be preserved. We found that the U-network style allowed us to train our
network on considerably smaller datasets then reported in the previous work.
The use of gated residual blocks also sped up training considerably. For input
into the network we use a simple binary representation of the geometry instead
of the Signed Distance Function representation used in the previous work as we
found no benefit in this added complexity. The steady state heat prediction
network uses the same basic network architecture and a complete description of
all networks including the parametrization networks can be found in the
appendix in figure 7.

##  5 Experiments

In the following sections we subject our method and model to a variety of
tests in order to see its performance.

###  5.1 Datasets

To train the parameterization networks we generate a set of 10,000 examples
for each system consisting of a parameter vector and their corresponding
geometry. An example of what a heat sink geometry looks like can be found in
figure 2. We use the parameterization found in Lane & Marshall (2009) for the
2D and 3D airfoils with 46 parameters that correspond to coefficients of a
polynomial describing the upper and lower surface of the foil. A complete
description of the parameterization can be found in the appendix.

The simulation datasets consists of 5,000, 5,000, and 2,500 training examples
for the heat sink simulation, 2D fluid simulation, and 3D fluid simulation
respectively. We use a 80-20 split in making the train and test sets. The
geometries used for the simulations are drawn from the distributions used in
the parameterization dataset. The heat simulations used a finite difference
solver and the fluid flow simulation used the Lattice Boltzmann method.

###  5.2 Training

We used the Adam optimizer for all networks (Kingma & Ba, 2014). For
$S^{heat}_{net}$ and $S^{flow}_{net}$ a learning rate of 1e-4 was used until
the loss plateaued and then the learning rate was dropped to 1e-5. Mean
Squared Error was used as the loss function however when training the flow
prediction network we scaled up the loss from the pressure field by a factor
of 10 to roughly match the magnitude of the velocity vector field. The
parameterization networks also used Mean Squared Error with a constant
learning rate of 1e-4. We found the parameterization networks trained
extremely quickly.

###  5.3 Gradient Decent Design Optimization Details

There are some complexities in how exactly the design parameters are optimized
that need explanation. The most naive approach is to scale the computed
gradient by some learning rate and add it to the design parameters. We found
this approach effective however it was prone to finding local optimum. We
found that adding momentum to the gradient reduced the chance of this and in
most cases accelerated optimization. We also found that adding a small amount
of noise too the parameters when computing gradients helped jump out of local
optima. We used momentum 0.9 and a learning rate of 0.05 and 0.001 for the
heat sink and airfoil problems respectively. The noise added to the parameters
used a Gaussian distribution with mean 0 and standard deviation 0.01.

If the above approach is used naively it can result in parameter values
outside of the original design space. To solve this problem we scale the input
to the parameterization network between 0 and 1 and use a hard sigmoid to
enforce this. This does not fix the problem completely though because if the
parameters being trained leave the range of -0.5 to 0.5, the gradient will be
zero and the parameter will be stuck at its current value. To prevent this we
simply add a small loss that pushes any parameters outside the -0.5 to 0.5
range back in.

###  5.4 Heat Sink Optimization

As discussed above, the heat sink optimization task is to find a set of fin
heights that maximally cool a constant heat source given a fixed total length
of the fins. The set up roughly corresponds to an aluminum heat sink placed on
a CPU where the heat source is treated as a continual addition of temperature.
There is no heat dissipation between the underside of the heat sink but all
other areas not on the heat sink are kept at a constant temperature. The
intuitive solution to this optimization problem is to place long fins near the
heat source and shorter fins farther away. Balancing this is a difficult task
though because changing the length of any fin has a global effect on how much
heat is dissipated by all the other fins.

After training our networks $P^{heat}_{net}$ and $S^{heat}_{net}$ we perform
our proposed gradient optimization on the 15 fin heights to minimize the
temperature at the source. In figure 2 we see the optimized heat sink and
observe that the design resembles what our intuition tells us. We also note
the extremely smooth optimization that occurs with only small bumps caused by
the addition of noise noted above. A natural question to ask is how this
compares to other search techniques. In order to answer these questions we use
simulated annealing to search designs and use the original heat diffusion
solver to evaluate their performance. In figure 2, we see that the optimized
heat sink design produced by the neural network closely resembles that
produced by simulated annealing. There are some minute differences however the
total effectiveness in cooling the system are almost identical. We also note
the iteration difference between the two methods. The gradient decent approach
required roughly 150 iterations to converge where as the simulated annealing
approach needed at least 800.

Figure 2: The top plot shows the optimization process and final design using
our gradient descent method. The bottom plot shows a comparison of our
optimization method to simulated annealing and the differences in final
designs. As we can see, the gradient descent method converges much faster and
finds roughly the same design.

###  5.5 Flow Prediction Accuracy

Before we move to our final test of designing 2D and 3D airfoils it is
important to know how accurately our model can predict steady state fluid
flow. We can also verify our claim of a superior network architecture over
previous work and show results indicating this. We omitted this discussion of
accuracy from the heat sink problem however a figure showing the accuracy in
predicting the heat at source can be found in figure 8 in the appendix.

The quantities of most interest in our predictions are the forces on the
object. These are the values being optimized so being able to predict them
accurately is of crucial importance. The forces are calculated from the
pressure field by doing a surface integral over the airfoil. This can be done
with any neural network library in a differentiable way by using a 3 by 3
transpose convolution on the boundary to determine the surface normals of the
object. Then multiplying this with the pressure field and summing to produce
the total force. Viscus forces are left out from this calculation as they are
relatively small for thin airfoils. In figure 3, we see that our model is very
accurate in predicting the forces. When comparing our network to the previous
model we see a clear increase in accuracy. We also visually inspect the flow
and see that the predicted flow is very sharp and doesn’t have any rough or
blurring artifacts.

Figure 3: Comparison of steady state flow predicted by neural network and the
Lattice Boltzmann flow solver. The right plot compares our network
architecture (blue dots) with Guo et al. (2016) (green dots). As we can see,
our network predicts forces and the max velocity more accurately then the
other model.

###  5.6 Automated Design of 2D and 3D Airfoils

Figure 4: A 2D airfoil designed by our gradient descent method. The airfoil
works by producing a low pressure area above its surface.

Figure 5: A 3D airfoil designed by our gradient descent method. The cross
section of the wing is very similar to that seen in the 2D airfoil.

A conventional airfoil works by using a curved upper surface to create a low
pressure zone and produce lift. The most important quantity for an airfoil is
its lift drag ratio which in effect tells its efficiency. At different angles
with respect to the fluid flow (angles of attack) the airfoil will produce
different lift drag ratios. Roughly speaking, an airfoil should have a
increase in lift drag ratio as the angle of attack increases until a max value
is reached. For our optimization task, we maximize this lift drag ratio for an
airfoil at angles of attack ranging from -5 to 17.5 degrees. The gradient for
the airfoil is calculated 9 times at angles in this range and then combined to
produce one gradient update. This approach of multiple angle optimization is
common and can be found in Drela (1998). In figure 4 and 5 we see the
optimized designs produced for the 2D and 3D simulations. We see that our
method produces the expected shape and characteristic curve of lift drag ratio
versus angle of attack. We also simulated the optimized airfoil with the
Lattice Boltzmann solver and found that it performed similarly confirming that
optimized designs produced by our method translate well to the original
simulation.

We have seen that our method is quite effective at producing optimized designs
but it is worth investigating what the fitness space looks like. To do this we
selected a random airfoil and slowly changed one of its parameters to see the
effect on the lift drag ratio. A plot of this can be seen in figure 6. We
notice that while there are many local optima present, the change in lift drag
ratio is very smooth and produces a very clean gradient. We view this as the
reason our method optimizes so quickly. We found that local optima like the
ones seen in this plot did not pose a serious problem during the optimization
and when running multiple times with different starting designs the same basic
shape was found with similar total fitness. We believe this was a result of
both the momentum and addition of noise as well as optimizing multiple angles
of attack at once. Adding this multiple angle constraint limits the number of
possible designs and makes the chance of finding local optima smaller. We
leave a deeper investigation into the effect of local optima for future work.

Similar to the heat sink problem, we compare our gradient decent method to
simulated annealing. Unlike the heat sink problem though, performing simulated
annealing with the Lattice Boltzmann solver was too computationally demanding
and so we used our network to evaluate the lift drag ratio instead. We see
from the figure 6 that using the gradient accelerates the optimization and in
only 200 iterations it converges. In comparison, the simulated annealing
requires at least 1500 iterations to reach similar performance.

Figure 6: The left plot shows the change in lift drag ratio versus a change in
a single design parameter. We note that while there are many local optima, the
surface is very smooth and produces clean gradients. The right plot shows a
comparison of the gradient decent optimization to simulated annealing in the
2D airfoil problem for a variety of starting temperatures.

###  5.7 Comparison of Computation Times

The central purpose of our method is to accelerate the automated design
process and in this section we attempt to quantify this in real time. The most
important quantities are the time to perform a gradient update on the design
parameters and the time needed to perform a simulation. Using these values we
can give a very rough comparison of optimization using our method and other
gradient free methods that use the flow solver. We leave this section for the
airfoil design problems only.

The first quantity we look at is the raw speed of the fluid solver. We found
that our flow solver converged to steady state in an average of 37.8 seconds
for the 2D simulation and 163.8 seconds for the 3D simulation on a Nvidia 1080
GPU. We used the Sailfish library for these simulations as it performed faster
then every other non-proprietary Lattice Boltzmann based fluid flow library
(Januszewski & Kostur, 2014). In comparison to our neural network, performing
one gradient update required only 0.052 seconds for the 2D simulation and
0.711 seconds for the 3D simulation. A more complete list of values including
their relation to batch size can be found in the table 1 in the appendix.
Given that performing automated design on the 2D airfoil required roughly
1,500 iterations at 9 different angles, this represents a total computation
time of 141 hours. In comparison, our method only took 1.5 minutes to perform
its 200 iterations at the 9 angles of attack. While this does represent a
significant 5,000 times speed increase, we note that there are several methods
of accelerating Lattice Boltzmann steady state flow calculations not explored
in this work that under restricted conditions can give a significant speed
increase (Guo & Shu, 2013) (Bernaschi et al., 2002). We also note that there
are other applicable search methods such as genetic algorithms and particle
swarm methods that may be more sample efficient. With the understanding that
this comparison is somewhat rough, we view this result as strong evidence that
our novel method is able to overcome some of the current computational
limitations faced in automated design.

##  6 Conclusion

In this work we have presented a novel method for automated design and shown
its effectiveness on a variety of tasks. Our method makes use of neural
networks and gradient descent to provide powerful and fast optimization. There
are many directions for future work such as applying this method to new
domains like structural optimization and problems related to electromagnetism.
One area of particular interest is design optimization on airfoils in
turbulent time dependent flows. Another interesting area to explore is hybrid
approaches where the neural network method is used to generate a rough design
and then fine tuned with a high fidelity simulation.

#### Acknowledgments

This work was made possible through the http://aigrant.org created by Nat
Friedman and Daniel Gross. This work would not have be possible without this
very generous support.

## References

  * Anderson & Aftosmis (2015) George R Anderson and Michael J Aftosmis.  Adaptive shape parameterization for aerodynamic design.  _Nat. Aeronaut. Space Admin., Ames Res. Center, Moffett Field, CA, USA, NAS Tech. Rep. NAS-2015-02_ , 2015.
  * Ando et al. (2010) Kenichi Ando, Akio Takamura, and Isao Saito.  Automotive aerodynamic design exploration employing new optimization methodology based on cfd.  _SAE International Journal of Passenger Cars-Mechanical Systems_ , 3(2010-01-0513):398–406, 2010.
  * Bernaschi et al. (2002) Massimo Bernaschi, Sauro Succi, Hudong Chen, and Raoyang Zhang.  Computing steady state flows with an accelerated lattice boltzmann technique.  _International Journal of Modern Physics C_ , 13(05):675–687, 2002.
  * Drela (1998) Mark Drela.  Pros and cons of airfoil optimization.  _Frontiers of computational fluid dynamics_ , 1998, 1998.
  * Dyck & Lowther (1996) Derek N Dyck and David A Lowther.  Automated design of magnetic devices by optimizing material distribution.  _IEEE Transactions on Magnetics_ , 32(3):1188–1193, 1996.
  * Guo et al. (2016) Xiaoxiao Guo, Wei Li, and Francesco Iorio.  Convolutional neural networks for steady flow approximation.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pp. 481–490. ACM, 2016.
  * Guo & Shu (2013) Zhaoli Guo and Chang Shu.  _Lattice Boltzmann method and its applications in engineering_ , volume 3.  World Scientific, 2013.
  * Hennigh (2017) O. Hennigh.  Lat-Net: Compressing Lattice Boltzmann Flow Simulations using Deep Neural Networks.  _ArXiv e-prints_ , May 2017.
  * Hilton (2007) Reno Hilton.  A universal parametric geometry representation method–“cst”.  2007\.
  * (10) Gregory S Hornby, Al Globus, Derek S Linden, and Jason D Lohn.  Automated antenna design with evolutionary algorithms.
  * Januszewski & Kostur (2014) Michal Januszewski and Marcin Kostur.  Sailfish: A flexible multi-gpu implementation of the lattice boltzmann method.  _Computer Physics Communications_ , 185(9):2350–2368, 2014.
  * Kingma & Ba (2014) Diederik Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Koziel & Leifsson (2013) Slawomir Koziel and Leifur Leifsson.  Multi-level cfd-based airfoil shape optimization with automated low-fidelity model selection.  _Procedia Computer Science_ , 18:889–898, 2013.
  * Lane & Marshall (2009) Kevin A Lane and David D Marshall.  A surface parameterization method for airfoil optimization and high lift 2d geometries utilizing the cst methodology.  _AIAA Paper_ , 1461:5–8, 2009.
  * Mills et al. (2017) K Mills, M Spanner, and I Tamblyn.  Deep learning and the schr$\backslash$” odinger equation.  _arXiv preprint arXiv:1702.01361_ , 2017.
  * Mukesh et al. (2012) R Mukesh, K Lingadurai, and U Selvakumar.  Application of nontraditional optimization techniques for airfoil shape optimization.  _Modelling and Simulation in Engineering_ , 2012:46, 2012\.
  * Paganini et al. (2017) M. Paganini, L. de Oliveira, and B. Nachman.  CaloGAN: Simulating 3D High Energy Particle Showers in Multi-Layer Electromagnetic Calorimeters with Generative Adversarial Networks.  _ArXiv e-prints_ , May 2017.
  * Piggott et al. (2015) Alexander Y Piggott, Jesse Lu, Konstantinos G Lagoudakis, Jan Petykiewicz, Thomas M Babinec, and Jelena Vučković.  Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer.  _Nature Photonics_ , 9(6):374–377, 2015.
  * Ronneberger et al. (2015) Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  _CoRR_ , abs/1505.04597, 2015.  URL http://arxiv.org/abs/1505.04597.
  * Salimans et al. (2017) Tim Salimans, Andrej Karpathy, Xi Chen, and Diederik P Kingma.  Pixelcnn++: Improving the pixelcnn with discretized logistic mixture likelihood and other modifications.  _arXiv preprint arXiv:1701.05517_ , 2017.
  * Salunke et al. (2014) Nilesh P Salunke, Juned Ahamad RA, and SA Channiwala.  Airfoil parameterization techniques: A review.  _American Journal of Mechanical Engineering_ , 2(4):99–102, 2014.
  * Tompson et al. (2016) Jonathan Tompson, Kristofer Schlachter, Pablo Sprechmann, and Ken Perlin.  Accelerating eulerian fluid simulation with convolutional networks.  _arXiv preprint arXiv:1607.03597_ , 2016.

##  7 Appendix

Figure 7: Complete network architectures. Figure 8: Difference between predicted heat at source of heat sink for our network and the original solver. Table 1: Complete list of computation times for different network evaluations with different batch sizes. Batch Size | 1 | 2 | 4 | 8 | 16
---|---|---|---|---|---
Flow Net $512^{2}$ | 0.020 sec | 0.017 sec | 0.016 sec | 0.015 sec | 0.015 sec
Param Net $512^{2}$ | 0.004 sec | 0.002 sec | 0.002 sec | 0.002 sec | 0.001 sec
Learn Step $512^{2}$ | 0.080 sec | 0.065 sec | 0.057 sec | 0.053 sec | 0.052 sec
Flow Net $144^{3}$ | 0.192 sec | 0.192 sec | 0.193 sec | 0.197 sec | Nan
Param Net $144^{3}$ | 0.025 sec | 0.024 sec | 0.025 sec | 0.025 sec | 0.025 sec
Learn Step $144^{3}$ | 0.911 sec | 0.711 sec | Nan | Nan | Nan

###  7.1 Airfoil Parameterization

The equation that parameterizes the upper and lower surface of the 2D airfoil
is

$S(x)=x^{n_{1}}(1-x)^{n_{2}}\sum_{i=0}^{N}A_{i}\frac{N}{i!(N-i)!}x^{i}(1-x)^{N-1}+hx$
(2)

The parameters present are $n_{1}$, $n_{2}$, $A_{i}$s, and $h$. We also add
the parameter $\theta$ that determines the angle of attack. In this work we
fixed $n_{1}$ to 0.5 and $n_{2}$ to 1.0 as this will produce a rounded head on
the airfoil. We also fix $h$ to zero making the tail the same height as the
head. Thus the trainable parameters are the 42 values corresponding to the
$A_{i}$s for the upper and lower surface. A illustration showing the
parameterization can be found in figure 9. The 3D airfoil has similar
parameterization.

$S(\phi,y)=\phi^{n_{1}}(1-\phi)^{n_{2}}\sum_{i=0}^{N_{x}}\sum_{j=0}^{N_{y}}A_{i}B_{j}\frac{N_{x}}{i!(N_{x}-i)!}\phi^{i}(1-\phi)^{N_{x}-1}\frac{N_{y}}{i!(N_{y}-i)!}y^{i}(1-y)^{N_{y}-1}+h\phi$
(3)

Where

$\phi=\frac{x-sy}{(l-0.5)y+0.5}$ (4)

This tells the height of the airfoil at a point $(x,y)$. The trainable
parameters here are $n_{1}$, $n_{2}$, $A_{i}$s, $B_{j}$s, $h$, $s$, and $l$.
Again, $n_{1}$, $n_{2}$, and $h$ are fixed to the values in the 2D case. We
also have 2 parameters for the angle $\theta$ and $\psi$ that determine the
rotation in the x and y direction. We keep $\psi$ at zero and only vary
$\theta$ at the desired angles during the optimization. The parameters $s$ and
$l$ correspond to the sweep present in the wing. This leaves the $A_{i}$s and
$B_{j}$s for optimization. We split the remaining 39 parameters equally so
that 13 values are used for $B_{i}$s and the remaining 26 are split between
the $A_{i}$s for the upper and lower surface. For a much more in depth look at
this parameterization, see Lane & Marshall (2009).

Figure 9: Example of a 2D airfoil.


# Learning Undirected Graphs in Financial Markets
††thanks: This work was supported by the Hong Kong GRF 16207019 research
grant.

José Vinícius de Miranda Cardoso  Department of Electronic and Computer
Engineering
The Hong Kong University of Science and Technology
Clear Water Bay, Hong Kong
<EMAIL_ADDRESS>Daniel P. Palomar  Department of Electronic
and Computer Engineering
Department of Industrial Engineering and Decision Analytics
The Hong Kong University of Science and Technology
Clear Water Bay, Hong Kong
<EMAIL_ADDRESS>
###### Abstract

We investigate the problem of learning undirected graphical models under
Laplacian structural constraints from the point of view of financial market
data. We show that Laplacian constraints have meaningful physical
interpretations related to the market index factor and to conditional
correlations between stocks. Those interpretations lead to a set of guidelines
that users should be aware of when estimating graphs in financial markets. In
addition, we propose algorithms to learn undirected graphs that account for
stylized facts and tasks intrinsic to financial data such as non-stationarity
and stock clustering.

###### Index Terms:

undirected graphical models, graph Laplacian, stock markets

##  I Introduction

Learning the structure of general graphical models is an NP-hard task [1]
whose importance is critical towards visualizing, understanding, and
leveraging the full potential contained in the data that live in such
structures. Learning graphs from data is a fundamental problem in the
statistical graph learning and signal processing fields [2, 3, 4], having a
direct impact on applications such as unsupervised learning, clustering, and
applied finance [5, 6, 7].

Nonetheless, most existing techniques for learning graphs are often unable to
impose a particular graph structure due to their inability to incorporate
prior information in the learning process. Moreover, most graph learning
frameworks are designed towards static networks [8, 9, 3, 4] which inherently
neglect dynamic time-domain variations in real data. As a consequence, they
usually lack practicality especially in non-stationary data regimes, which is
often the case in data from financial stock markets.

Motivated by practical applications such as clustering of stocks and
understanding their time-domain variations, we investigate the problem of
estimating graph matrices whose structure follow those of Laplacian matrices
of undirected weighted graphs in a financial context both for static and
dynamic graphs.

The main contributions of our paper are as follows:

  1. 1.

We for the first time provide natural interpretations for the Laplacian
constraints of graphs estimated from stock market data. This leads to
meaningful and intuitive guidelines on the data processing required prior to
learning graphs.

  2. 2.

We show that rank constraints alone, a practice often used by state-of-the-art
methods, are not sufficient to learn $k$-component graphs.

  3. 3.

We propose novel formulations to learn: i) $k$-component graphs and ii) time-
varying graphs.

  4. 4.

We develop simple trading strategies as a result of the estimated time-varying
graphs.

##  II Background and Related Work

A graph is denoted as a triple
$\mathcal{G}=\left(\mathcal{V},\mathcal{E},\bm{W}\right)$, where
$\mathcal{V}=\left\\{1,2,\dots,p\right\\}$ is the vertex (or node) set,
$\mathcal{E}\displaystyle\subseteq\left\\{\left\\{u,v\right\\}:u,v\in\mathcal{V}\right\\}$
is the edge set, that is, a subset of the set of all possible unordered pairs
of $p$ nodes such that $\\{u,v\\}\in\mathcal{E}$ if and only if nodes $u$ and
$v$ are connected. We denote the number of elements in $\mathcal{E}$ by
$|\mathcal{E}|$. $\bm{W}\in\mathbb{R}_{+}^{p\times p}$ is the symmetric
weighted adjacency matrix that satisfies $W_{ii}=0,W_{ij}>0~{}\text{if and
only
if}~{}\\{i,j\\}\in\mathcal{E}~{}\text{and}~{}W_{ij}=0~{}\text{otherwise}$. The
graph Laplacian matrix $\bm{L}$ is defined as $\bm{L}\triangleq\bm{D}-\bm{W}$,
where $\bm{D}\triangleq\mathsf{Diag}(\bm{W}\mathbf{1})$ is the degree matrix.

An attractive improper Gaussian Markov Random Field (GMRF) [10] is denoted as
a $p$-dimensional, real-valued, Gaussian random variable $\bm{x}$ with mean
vector $\bm{\mu}$ and low-rank precision (inverse covariance) matrix
$\bm{\Xi}$. The data generating process is assumed to be a zero-mean,
attractive improper GMRF $\bm{x}\in\mathbb{R}^{p}$, such that $x_{i}$ is the
random variable generating a signal measured at node $i$, whose low-rank
precision matrix is modeled as a graph Laplacian matrix. Assume we are given
$n$ observations from $\bm{x}$, i.e., $\bm{X}\in\mathbb{R}^{n\times p}$,
$\bm{X}=\left[\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{p}\right]$,
$\bm{x}_{i}\in\mathbb{R}^{n}$. The Penalized Maximum Likelihood Estimator of
the precision matrix of $\bm{x}$, on the basis of the observed data $\bm{X}$,
is

$\begin{array}[]{ll}\underset{\bm{L}\succeq\mathbf{0}}{\mathsf{minimize}}&\mathsf{tr}\left(\bm{L}\bm{S}\right)-\mathsf{log~{}gdet}\left(\bm{L}\right)+h_{\bm{\alpha}}(\bm{L}),\\\
\mathsf{subject~{}to}&\bm{L}\mathbf{1}=\mathbf{0},~{}L_{ij}=L_{ji}\leq
0,\end{array}$ (1)

where $\bm{S}$ is a similarity matrix, such as sample covariance (or
correlation) matrix $\bm{S}\propto\bm{X}^{\top}\bm{X}$,
$\mathsf{gdet}(\bm{L})$ is the pseudo determinant of $\bm{L}$, i.e., the
product of its positive eigenvalues [11], and $h_{\bm{\alpha}}(\bm{L})$ is a
regularization function, with hyperparameter vector $\bm{\alpha}$, to promote
certain properties on $\bm{L}$ such as sparsity.

Problem (1) is a fundamental problem in the graph signal processing field that
has served as a cornerstone for many extensions, primarily those involving the
inclusion of structure onto $\bm{L}$ [3, 6]. Even though Problem (1) is
convex, assuming $h_{\bm{\alpha}}(\cdot)$ is convex, it is not adequate to be
solved by disciplined convex programming languages, such as $\mathsf{cvxpy}$
[12], due to scalability issues related to the computation of
$\mathsf{log~{}gdet}(\bm{L})$. Indeed, recently a few works have proposed
algorithms based on Majorization-Minimization (MM) [13], and Alternating
Direction Method of Multipliers (ADMM) [14] to solve Problem (1) in an
efficient and scalable fashion [3, 4].

On the other hand, due to nuisances involved in dealing with the term
$\mathsf{log~{}gdet}(\bm{L})$, several works focused on the assumption that
the underlying signals in a graph are smooth [9, 8]. In its simplest form,
learning a smooth graph from a data matrix $\bm{X}\in\mathbb{R}^{n\times p}$
is tantamount to finding an adjacency matrix $\bm{W}$ that minimizes the
Dirichlet energy:

$\begin{array}[]{ll}\underset{\bm{W}}{\mathsf{minimize}}&\frac{1}{2}\sum_{i,j}W_{ij}\left\lVert\bm{x}_{i}-\bm{x}_{j}\right\rVert^{2}_{2}\\\
\mathsf{subject~{}to}&W_{ij}=W_{ji}\geq
0,\mathsf{diag}(\bm{W})=\mathbf{0}.\end{array}$ (2)

In order for Problem (2) to be well-defined, several constraints have been
proposed in the literature. For instance, Kalofolias et al. [9] proposed a
convex formulation as follows

$\begin{array}[]{ll}\underset{\bm{W}}{\mathsf{minimize}}&\frac{1}{2}\mathsf{tr}\left({\bm{W}\bm{Z}}\right)-\alpha\mathbf{1}^{\top}\log(\bm{W}\mathbf{1})+\frac{\gamma}{2}\left\lVert\bm{W}\right\rVert^{2}_{\mathrm{F}},\\\
\mathsf{subject~{}to}&W_{ij}=W_{ji}\geq
0,\mathsf{diag}(\bm{W})=\mathbf{0},\end{array}$ (3)

where $Z_{ij}\triangleq\left\lVert\bm{x}_{i}-\bm{x}_{j}\right\rVert^{2}_{2}$.

Problem (3) is convex and can be solved via primal-dual, ADMM-like algorithms
[9]. It can be seen that the objective function in Problem (3) is actually an
approximation to that of Problem (1) where the $\mathsf{gdet}(\bm{L})$ term
has been upper bounded by $\prod_{i=1}^{p}L_{ii}$. Therefore, Problem (3) can
be thought of as an approximation of the penalized maximum likelihood
estimator.

A formulation to estimate $k$-component graphs based on the smooth signal
approach was proposed in [5]. They proposed a two-stage algorithm where it
first estimates a connected graph using, e.g., Problem (3) and then it
projects the graph onto the set of Laplacian matrices of dimension $p$ with
rank $p-k$, where $k$ is the given number of graph components.

Spectral constraints on the Laplacian matrix are an intuitive way to recover
$k$-component graphs as the multiplicity of its zero eigenvalue, i.e., the
nullity of $\bm{L}$, dictates the number of components of a graph. The first
proposed framework to impose structures on the estimated Laplacian matrix,
under the attractive improper GMRF model, was done by Kumar et al. [15, 6],
through the use of spectral constraints, as follows

$\begin{array}[]{ll}\underset{\bm{L},\bm{U},\bm{\lambda}}{\mathsf{minimize}}&\mathsf{tr}\left(\bm{L}\bm{S}\right)-\displaystyle\sum_{i=1}^{p-k}\mathsf{log}\left(\lambda_{i}\right)+\frac{\eta}{2}\left\lVert\bm{L}-\bm{U}\mathsf{Diag}(\bm{\lambda})\bm{U}^{\top}\right\rVert^{2}_{\mathrm{F}},\\\
\mathsf{subject~{}to}&\bm{L}\succeq\mathbf{0},~{}\bm{L}\mathbf{1}=\mathbf{0},~{}L_{ij}=L_{ji}\leq
0,\\\ &\bm{U}^{\top}\bm{U}=\bm{I},\bm{U}\in\mathbb{R}^{p\times(p-k)},\\\
&0<c_{1}<\lambda_{1}<\dots<\lambda_{p-k}<c_{2}.\end{array}$ (4)

Note that Problem (4) learns a $k$-component graph without the need for a two-
stage algorithm. However, a clear shortcoming of this formulation is that it
does not control the degrees of the nodes in the graph, which may result in a
trivial solution that contains isolated nodes, turning out not to be useful
for clustering tasks especially when applied to noisy data sets.

##  III Graph Laplacian Constraints Interpretation for Stock Signals

Graphical representations of data are increasingly important tools in
financial signal processing to uncover hidden relationships between variables
[7]. In stock markets, one is generally interested in learning about
conditional dependencies among stocks and how to leverage this information
into practical scenarios such as portfolio design and crisis forecasting.

Mathematically, we would like to estimate a precision matrix $\bm{L}$ that
enjoys the following two key properties

(P1)

$\bm{L}\mathbf{1}=\mathbf{0}$,

(P2)

${L}_{ij}={L}_{ji}\leq 0~{}\forall~{}i\neq j$.

The first property states that the Laplacian matrix $\bm{L}$ is singular and
its null space contains the $\mathbf{1}$ vector. That means that any signal
sampled from $\bm{L}$, say $\bm{x}$, is constrained to a subspace of rank
$p-1$ satisfying $\bm{x}^{\top}\mathbf{1}=0$. In practice, (P1) implies that
signals living in a graph $\mathcal{G}$ have zero graph-mean.

Property (P2) together with (P1) implies that $\bm{L}$ is positive
semidefinite. The fact that the off diagonal entries are non-positive means
that the Laplacian matrix only represents non-negative conditional
dependencies111The correlation between any two pair of nodes conditioned on
the rest of the graph is given as $-\frac{L_{ij}}{\sqrt{L_{ii}L_{jj}}}$.. This
assumption is often met for stock data as assets are typically positively
dependent [16, 17].

These two properties along with efficient learning frameworks make the
Laplacian-based graphical model a natural candidate for learning graphs of
stock data. As a consequence of using the Laplacian model, we propose the
following guidelines when estimating Laplacian matrices with stock market
data.

  * •

Correlation vs Covariance: Both the GMRF and smooth signal approaches rely on
the term $\mathsf{tr}{(\bm{S}\bm{L})}\propto\mathsf{tr}(\bm{W}\bm{Z})$, where
$\bm{S}$ is the sample covariance matrix. From the definition of $\bm{Z}$ (3),
we observe that two perfectly correlated stocks but with large Euclidean
distances would appear very different on the graph. Hence, we should use the
sample correlation matrix
$\bar{\bm{S}}=\mathsf{Diag}(\bm{S})^{-1/2}\bm{S}\mathsf{Diag}(\bm{S})^{-1/2}$
in case we want two highly correlated stocks to have a strong graph connection
regardless of their variances.

  * •

Removing the market trend: A widely used and tested model for the returns of
the stocks includes explicitly the dependency on the market factor:
$\bm{x}_{t}=\bm{\beta}x_{\mathsf{mkt},t}+\bm{\epsilon}_{t}$, where
$\bm{\beta}$ denotes the market loadings, $x_{\mathsf{mkt},t}$ denotes the
market index, and $\bm{\epsilon}_{t}$ is the residual idiosyncratic component
with covariance matrix $\bm{\Psi}$. Since all the stocks are heavily dominated
by the market index $x_{\mathsf{mkt},t}$, it may be convenient to remove that
component if we seek to explore the structure of the residual cross-dependency
among the stocks, $\bm{\epsilon}_{t}$. Thus, an alternative to using the full
covariance matrix $\bm{\Sigma}$ is to use the covariance matrix $\bm{\Psi}$ of
the idiosyncratic component. However, since $\bm{\beta}\approx\mathbf{1}$, it
turns out that the market factor is automatically removed in the squared
distance matrix $\bm{Z}$:

$Z_{ij}\triangleq\|\bm{x}^{(i)}-\bm{x}^{(j)}\|^{2}_{2}\approx\|\bm{\epsilon}^{(i)}-\bm{\epsilon}^{(j)}\|^{2}_{2}.$
(5)

Even more interestingly, if one first normalizes each stock, whose variances
are $\mathbb{V}(\bm{x}^{(i)})\approx\beta^{2}_{i}$, we have
$\bar{\bm{x}}=\mathbf{1}\bar{x}_{\mathsf{mkt}}+\bar{\bm{\epsilon}}_{t}$, which
implies that the market factor is automatically removed in the squared
distance matrix.

  * •

Degree control: Enforcing a rank smaller than $p-1$ for the Laplacian matrix
will generate a $k$-component graph, which is one desired goal. However, one
may get the undesired result of having isolated nodes. The way to avoid
isolated nodes is by controlling the degrees of the nodes. The GMRF
formulation has the natural penalty term $\mathsf{log~{}gdet}(\bm{L})$ in the
objective, but that does not help in controlling the degrees of the nodes.
Instead, some of the graph learning formulations from smooth signals include
degree control via the constraint $\bm{W}\mathbf{1}=\mathbf{1}$, which fixes
the degrees of all the nodes to $1$. The regularization term
$\mathbf{1}^{\top}\mathsf{log}(\bm{W}\mathbf{1})$ also avoids the trivial
solution of any degree equals $0$. Hence, any graph learning formulation that
enforces a $k$-component graph (or low-rank Laplacian matrix) should also
control the degrees of the nodes to avoid a trivial solution with isolated
nodes.

##  IV Proposed Formulations

In this section, we propose graph learning formulations to account for: 1)
$k$-component structures and 2) non-stationarity of financial stock market
data.

###  IV-A $k$-component graphs: GMRF formulation

We propose the following formulation to learn a $k$-component graph:

$\begin{array}[]{ll}\underset{\bm{L}\succeq\mathbf{0}}{\mathsf{minimize}}&\mathsf{tr}\left(\bm{L}\bm{S}\right)-\mathsf{log~{}gdet}\left(\bm{L}\right)\\\
\mathsf{subject~{}to}&\bm{L}\mathbf{1}=\mathbf{0},~{}L_{ij}=L_{ji}\leq
0,~{}\forall~{}i\neq j,\\\
&\mathsf{diag}(\bm{L})=\mathbf{1},\mathsf{rank}(\bm{L})=p-k.\end{array}$ (6)

Problem (6) is highly non-convex and non-differentiable due to the constraint
$\mathsf{rank}(\bm{L})=p-k$, which is equivalently to saying that the sum of
the $k$ smallest eigenvalues of $\bm{L}$ is equal to zero, i.e.,
$\sum_{i=1}^{k}\lambda_{i}\left(\bm{L}\right)=0$ [5] (assuming eigenvalues in
increasing order). By Fan’s theorem [18],

$\sum_{i=1}^{k}\lambda_{i}\left(\bm{L}\right)=\underset{\bm{V}\in\mathbb{R}^{p\times
k},\bm{V}^{\top}\bm{V}=\bm{I}}{\mathsf{minimize}}\mathsf{tr}\left(\bm{V}^{\top}\bm{L}\bm{V}\right).$
(7)

Thus, a relaxed version of Problem (6) becomes

$\begin{array}[]{ll}\underset{\bm{L}\succeq\mathbf{0},\bm{V}\in\mathbb{R}^{p\times
k}}{\mathsf{minimize}}&\mathsf{tr}\left(\bm{L}\bm{S}\right)-\mathsf{log~{}gdet}\left(\bm{L}\right)+\eta\mathsf{tr}\left(\bm{V}^{\top}\bm{L}\bm{V}\right),\\\
\mathsf{subject~{}to}&\bm{L}\mathbf{1}=\mathbf{0},~{}L_{ij}=L_{ji}\leq
0,~{}\forall~{}i\neq j,\\\
&\mathsf{diag}(\bm{L})=\mathbf{1},\bm{V}^{\top}\bm{V}=\bm{I}.\end{array}$ (8)

Even though still non-convex, but now differentiable, Problem (8) can be
solved in an alternating fashion. More precisely, for a given $\bm{L}$, say
$\bm{L}^{l}$, we have the following subproblem for $\bm{V}$:

$\begin{array}[]{ll}\underset{\bm{V}\in\mathbb{R}^{p\times
k}}{\mathsf{minimize}}&\mathsf{tr}\left(\bm{V}^{\top}\bm{L}^{l}\bm{V}\right),\\\
\mathsf{subject~{}to}&\bm{V}^{\top}\bm{V}=\bm{I},\end{array}$ (9)

whose solution is given by the $k$ eigenvectors associated with the $k$
smallest eigenvalues of $\bm{L}^{l}$ [19].

For a fixed value of $\bm{V}$, say $\bm{V}^{l}$, we have the following
subproblem for $\bm{L}$:

$\begin{array}[]{ll}\underset{\bm{L}\succeq\mathbf{0}}{\mathsf{minimize}}&\mathsf{tr}\left(\bm{L}\left(\bm{S}+\eta\bm{V}^{l}{\bm{V}^{l}}^{\top}\right)\right)-\mathsf{log~{}gdet}\left(\bm{L}\right)\\\
\mathsf{subject~{}to}&\bm{L}\mathbf{1}=\mathbf{0},L_{ij}=L_{ji}\leq
0,\mathsf{diag}(\bm{L})=\mathbf{1}.\end{array}$ (10)

Problem (10) is convex and can be solved efficiently by, e.g., the methods
developed in [3, 4, 6]. Algorithm 1 summarizes the implementation to solve
Problem (8).

Data: Similarity matrix $\bm{S}$, initial estimate $\bm{L}^{0}$, rank
constraint hyperparameter $\eta>0$.

Result: Laplacian estimation $\bm{L}$

1 while  _not converged_ do

2 update $\bm{V}^{l+1}$ by solving Problem (9) fixing $\bm{L}$ at $\bm{L}^{l}$

3 update $\bm{L}^{l+1}$ by solving Problem (10) fixing $\bm{V}$ at
$\bm{V}^{l+1}$

4

5 end while

Algorithm 1 GMRF $k$-component graph learning

###  IV-B Time-varying graphs

Most graph learning frameworks are designed towards static graphs, which
inherently neglect time-domain variations in real data. As a result, they
usually lack practicality for finance especially during nonstationary regimes,
e.g., an economic crisis.

A time-varying graph is defined as a sequence of graphs stacked over time,
i.e.,
$\left\\{\mathcal{G}_{t}\right\\}_{t=1}^{T}=\left\\{\mathcal{V}_{t},\mathcal{E}_{t},\bm{W}_{t}\right\\}_{t=1}^{T}$.
We assume the node set of each graph to be the same, i.e.,
$\mathcal{V}_{t}=\left\\{1,2,\dots,p\right\\},\forall~{}t=1,2,\dots,T$.

The data generating process is as follows. Assume that for every graph
$\mathcal{G}_{t}$ we associate an attractive improper GMRF
$\bm{x}_{t}\sim\mathcal{N}\left(\mathbf{0},\bm{L}^{\dagger}_{t}\right)$, where
$\bm{L}_{t}$ is the precision matrix of the $t$-th GMRF that is assumed to
have graph Laplacian structure. Further, suppose we are given $n_{t}$
observations from $\bm{x}_{t}$, i.e., $\bm{X}_{t}\in\mathbb{R}^{n_{t}\times
p}$, then we propose the following optimization program to learn the Laplacian
matrices $\left\\{\bm{L}_{t}\right\\}_{t=1}^{T}$ on the basis of
$\left\\{\bm{X}_{t}\right\\}^{T}_{t=1}$, in particular $\bm{L}_{t}$ is
obtained from

$\begin{array}[]{ll}\underset{\bm{L}_{1},\bm{L}_{2},\dots,\bm{L}_{t}}{\mathsf{minimize}}&\sum_{\tilde{t}=1}^{t}n_{\tilde{t}}\left[\mathsf{tr}\left(\bm{S}_{\tilde{t}}\bm{L}_{\tilde{t}}\right)-\mathsf{log~{}gdet}\left(\bm{L}_{\tilde{t}}\right)\right]\\\
&+\delta\sum_{\tilde{t}=2}^{t}d\left(\bm{L}_{\tilde{t}},\bm{L}_{\tilde{t}-1}\right),\\\
\mathsf{subject~{}to}&\left\\{\bm{L}_{\tilde{t}}\succeq\mathbf{0},\bm{L}_{\tilde{t}}\mathbf{1}=\mathbf{0},(\bm{L}_{\tilde{t}})_{ij}=(\bm{L}_{\tilde{t}})_{ji}\leq
0\right\\}_{\tilde{t}=1}^{t},\end{array}$ (11)

where $\\{\bm{S}_{t}\\}_{t=1}^{T}$ is a sequence of similarity matrices,
$d:\mathbb{R}^{p\times p}\times\mathbb{R}^{p\times
p}\rightarrow\mathbb{R}_{+}$ is a distance function that measures the
similarity between $\bm{L}_{t}$ and $\bm{L}_{t-1}$ in order to impose time
consistency, e.g.,
$d(\bm{L}_{t},\bm{L}_{t-1})\triangleq\|\bm{L}_{t}-\bm{L}_{t-1}\|^{2}_{\text{F}}$,
and $\delta\in\mathbb{R}_{++}$ is its corresponding hyperparameter. The
solution to Problem (11) is obtained on a rolling window basis, which is
tantamount to a causal estimator $\hat{\bm{L}}_{t|t},t=1,\dots,T$. In other
words, to estimate, e.g., $\bm{L}_{3}$, we only use information up to and
including time $t=3$, i.e., $\\{\bm{S}_{t}\\}_{t=1}^{3}$.

The time-varying graph learning formulations proposed in [20] and [21] are not
adequate to our particular scenario because those formulations are based on
the smooth signal approach, rather than the statistical GMRF model, and solve
a dynamic graph conditioned on all the $T$ chunks of observations, i.e.,
$\bm{\hat{L}}_{t\mid T},t=1,\dots,T$, which inevitably introduces look-ahead
biases.

Algorithm 2 summarizes the implementation to solve Problem (11).

Data: $\bm{S}_{1},\dots,\bm{S}_{T}$, $\\{n_{t}\\}^{T}_{t=1}$

Result: Causal Graph Laplacians estimates $\\{\bm{\hat{L}}_{t\mid
t}\\}^{T}_{t=1}$

1 for  _$t=1:T$_ do

2 $\triangleright$ estimate $\\{\bm{L}_{\tilde{t}}\\}_{\tilde{t}=1}^{t}$ using
data $\\{\bm{S}_{\tilde{t}}\\}^{t}_{\tilde{t}=1}$.

3 $\triangleright$ store $\bm{L}_{t}$ as $\bm{\hat{L}}_{t}$.

4

5 end for

return $\\{\bm{\hat{L}}\\}_{t=1}^{T}$.

Algorithm 2 Time-varying graph learning

Remark: for practical programming language implementation of the proposed
algorithms, one can take into account that the Laplacian matrix $\bm{L}$ is
symmetric and its diagonal elements are mappings of the off-diagonal ones.
Therefore, only elements above (or below) the diagonal of $\bm{L}$ need to be
updated. This approach is similar to that used by Kalofolias et al. [9] and
Kumar et al. [6].

##  V Numerical Experiments

In the experiments that follow, we use log-returns time-series data from
stocks belonging to three sectors of the S&P500 index, namely
$\mathsf{Industrials}$, $\mathsf{Consumer~{}Staples}$, and $\mathsf{Energy}$.
We collect price data of 130 stocks from Jan. 1st 2016 to Jan. 1st 2019 from
Yahoo! Finance, which represents 753 days worth of data.

###  V-A Effects of market factor and data scaling

In order to measure the effects of the market factor and the data scaling, we
learn four graphs using the two-stage algorithm proposed by [5].

Fig. 1a shows that removing the market and using the covariance matrix as
input leads to a graph with many possibly spurious connections (grey edges).
Fig. 1b shows that using the sample correlation matrix introduces an
improvement, but there exist still many spurious dependencies. Fig. 1c shows
that not removing the market presents an improvement, but the use of the
covariance matrix leads to many possibly fake connections. Fig. 1d combines
the correlation matrix as input and not removing the market, which clearly
shows a meaningful representation of a graph from stocks belonging to three
distinct sectors.

(a) Market removed and no data scaling.

(b) Market removed and data scaled.

(c) Market not removed and no data scaling.

(d) Market not removed and data scaling.

Figure 1: Graphs estimated with different data preprocessing. Panel (d) gives
the best results with scaled data and no market removed.

###  V-B Effect of degree control

To illustrate the importance of controlling the nodes degrees while learning
$k$-component graphs, we conduct a comparison between the algorithm proposed
in [6] and Algorithm 1 on the basis of the sample correlation matrix. Fig. 2
shows the estimated financial stock networks with $k=3$. It is clear that
$\mathsf{SGL}$ [6] (benchmark) is unable to account for the trivial solution
with isolated nodes, whereas the proposed algorithm returns a graph with a
meaningful representation.

(a) $\mathsf{SGL}$ [6].

(b) Algorithm 1 (proposed).

Figure 2: Effect of lack of degree control on the optimization formulation.
Benchmark algorithm $\mathsf{SGL}$, panel (a), fails to obtain a meaningful
solution, showing that spectral constraints alone are not sufficient to
recover a non-trivial $k$-component graph. The proposed algorithm in panel (b)
shows a meaningful estimated structure. We fixed the $\beta=10$ for both
methods.

###  V-C Time-varying Experiment

We consider an experiment with data from FAAMUNG companies (Facebook, Apple,
Amazon, Microsoft, Uber, Netflix, and Google) from June 1st 2019 to May 1st
2020, totalling 230 days worth of data, which includes the most recent
economic crisis due to the pandemic associated with COVID-19. For each period
of 30 days we estimate graphs on a rolling window fashion, shifting the window
one day at a time, on the basis of the sample correlation matrices
$\bm{\bar{S}}_{1},\dots,\bm{\bar{S}}_{200}$. At the end, we estimate 200
graphs, i.e., $\bm{\hat{L}}_{1},\dots,\bm{\hat{L}}_{200}$.

Based on this 200 estimated graphs we compute the algebraic connectivity,
i.e., the second smallest eigenvalue of $\bm{L}_{t}$, as an indicator of
variation of the graph. We use this indicator to acquire insights on possible
trends of the stock market. Note that other indicators could be used in
practice such as the spectral radius,
$\bm{\lambda}_{\mathsf{max}}(\bm{\hat{L}}_{t})$, and the time consistency,
$\|\bm{\hat{L}}_{t}-\bm{\hat{L}}_{t-1}\|^{2}_{\text{F}}$.

Fig. 3a shows the S&P500 log-price, where the impact of the COVID-19 pandemic
is clear around March 2020. Fig. 3b depicts the algebraic connectivity
indicator computed from each estimated graph along the time axis. Fig. 3c
shows the evolution of the graph network at certain dates. It is clear both
from the indicator (Fig. 3b) and the network visualization (Fig. 3c) that
around September 2019 the market has changed significantly. That is consistent
with news involving the impeachment inquiry of US President Donald J. Trump.
From the middle of March 2020 to the beginning of May the market saw its
largest drop since the financial crisis in 2008. This can also be noticed
through the indicator and the network visualization.

(a)

(b)

(c)

Figure 3: Estimated time-varying indicators and network visualization of
FAAMUNG companies during 2019-2020. Panel (a) shows the log-price of the
S&P500 index. Panel (b) show the algebraic connectivity indicator for
$\delta=100$. Panel (c) shows the estimated networks with $\delta=100$ for
several dates. It can be noticed that there is an increase in conditional
correlation among the stocks during times of economic crisis.

###  V-D Trading Application

By leveraging the proposed time-varying graph learning algorithm, we perform
an experiment comparing two simple trading strategies: (S1) uniformly invest a
unit of budget during the whole period; (S2) uniformly invest a unit of budget
according to whether or not the algebraic connectivity of the graph falls
below a fixed threshold of $\tau=1.0$. Fig. 4 shows the cumulative sum of the
profits and losses (PnL) over time for (S1) and (S2). It can be observed that
(S2) outperforms (S1) by smartly entering/exiting the market based on the
algebraic connectivity of the estimated graphs.

Figure 4: Uniform investment strategies on FAAMUNG companies. The cumulative
PnL of strategy (S2), which enters/exits the market according to the knowledge
of algebraic connectivity of the estimated graph, outperforms that of strategy
(S1), which continuously invests over the whole time period.

## References

  * [1] A. Anandkumar, V. Y. F. Tan, F. Huang, and A. S. Willsky, “High-dimensional Gaussian graphical model selection: Walk summability and local separation criterion,” _Journal of Machine Learning Research_ , vol. 13, no. 1, pp. 2293–2337, 2012.
  * [2] B. M. Lake and J. B. Tenenbaum, “Discovering structure by learning sparse graph,” in _Proceedings of the 33rd Annual Cognitive Science Conference_ , 2010.
  * [3] H. E. Egilmez, E. Pavez, and A. Ortega, “Graph learning from data under Laplacian and structural constraints,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 11, no. 6, pp. 825–841, 2017.
  * [4] L. Zhao, Y. Wang, S. Kumar, and D. P. Palomar, “Optimization algorithms for graph laplacian estimation via admm and mm,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 16, pp. 4231–4244, 2019.
  * [5] F. Nie, X. Wang, M. I. Jordan, and H. Huang, “The constrained Laplacian rank algorithm for graph-based clustering,” in _Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence_ , ser. AAAI’16. AAAI Press, 2016, pp. 1969–1976.
  * [6] S. Kumar, J. Ying, J. V. de M. Cardoso, and D. Palomar, “A unified framework for structured graph learning via spectral constraints,” _Journal of Machine Learning Research_ , vol. 21, pp. 1–60, 2020.
  * [7] G. Marti, F. Nielsen, M. Bińkowski, and P. Donnat, “A review of two decades of correlations, hierarchies, networks and clustering in financial markets,” in _arXiv: 1703.00485_ , 2017.
  * [8] X. Dong, D. Thanou, P. Frossard, and P. Vandergheynst, “Learning Laplacian matrix in smooth graph signal representations,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 23, pp. 6160–6173, 2016.
  * [9] V. Kalofolias, “How to learn a graph from smooth signals,” in _Proceedings of the 19th International Conference on Artificial Intelligence and Statistics_ , ser. Proceedings of Machine Learning Research, vol. 51. Cadiz, Spain: PMLR, 2016, pp. 920–929.
  * [10] H. Rue and L. Held, _Gaussian Markov Random Fields: Theory And Applications (Monographs on Statistics and Applied Probability)_. Chapman & Hall/CRC, 2005.
  * [11] O. Knill, “Cauchy–binet for pseudo-determinants,” _Linear Algebra and its Applications_ , vol. 459, pp. 522 – 547, 2014.
  * [12] S. Diamond and S. Boyd, “CVXPY: A Python-embedded modeling language for convex optimization,” _Journal of Machine Learning Research_ , vol. 17, no. 83, pp. 1–5, 2016.
  * [13] Y. Sun, P. Babu, and D. P. Palomar, “Majorization-minimization algorithms in signal processing, communications, and machine learning,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 3, pp. 794–816, 2017\.
  * [14] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” _Foundations and Trends in Machine Learning_ , vol. 3, no. 1, pp. 1–122, 2011.
  * [15] S. Kumar, J. Ying, J. V. de M. Cardoso, and D. P. Palomar, “Structured graph learning via laplacian spectral constraints,” in _Advances in Neural Information Processing Systems (NeurIPS)_ , 2019.
  * [16] R. Agrawal, U. Roy, and C. Uhler, “Covariance Matrix Estimation under Total Positivity for Portfolio Selection,” _arXiv e-prints_ , Sep. 2019\.
  * [17] Y. Wang, U. Roy, and C. Uhler, “Learning high-dimensional gaussian graphical models under total positivity without adjustment of tuning parameters,” in _Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics_ , vol. 108, 2020, pp. 2698–2708.
  * [18] K. Fan, “On a theorem of Weyl concerning eigenvalues of linear transformations i,” _Proceedings of the National Academy of Sciences_ , vol. 35, no. 11, pp. 652–655, 1949.
  * [19] P.-A. Absil, R. Mahony, and R. Sepulchre, _Optimization Algorithms on Matrix Manifolds_. Princeton, NJ: Princeton University Press, 2008.
  * [20] V. Kalofolias, A. Loukas, D. Thanou, and P. Frossard, “Learning time varying graphs,” in _IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2017, pp. 2826–2830.
  * [21] A. O. Koki Yamada, Yuichi Tanaka, “Time-varying graph learning with constraints on graph temporal variation,” in _arXiv: 2001.03346_ , 2020.