# basic finance
#library(xts) # to manipulate time series of stock data
#library(pob) # book package with financial data
# data
#library(RcppRoll) # rolling methods
#library(TTR) # exponentially moving average
#library(rugarch) # time series models
#library(MARSS) # Kalman
#library(stochvol) # to compute volatility envelope
# graphs
library(clusterSim) # graph datasets
library(igraph) # for graph related functions
library(dbscan) # for kNN computation
library(glasso) # GLASSO implementation
library(spectralGraphTopology) # implementation of many graph learning functions
# plotting
#library(ggplot2) # for nice plots
library(cowplot)
#library(patchwork) # for combining plots
#library(scales)R Code for Portfolio Optimization
Chapter 4 – Financial Data: Graphs
Daniel P. Palomar (2025). Portfolio Optimization: Theory and Application. Cambridge University Press.
Last update: March 13, 2025
Loading packages
The following packages are used in the examples:
For convenience, we define the following wrapper function (based on the package igraph) to plot graphs:
library(igraph)
# define function to plot nice graphs
plot_graph <- function(W, L = NULL, main = "", layout = NULL, clusters, edges = TRUE, labels = FALSE) {
# if the Laplacian is passed instead, then use it to get W
if (!is.null(L)) {
W <- abs(L)
diag(W) <- 0
}
# default layout
if(is.null(layout)) layout <- layout_with_fr
# build graph
graph <- graph_from_adjacency_matrix(W, mode = "undirected", weighted = TRUE)
# colorify nodes and edges
V(graph)$cluster <- clusters
colors <- c("#706FD3", "#FF5252", "#33D9B2")
V(graph)$color <- colors[clusters] #V(graph)$color <- clusters
if (edges)
E(graph)$color <- apply(as.data.frame(get.edgelist(graph)), 1,
function(x) {
x <- as.numeric(x)
ifelse(V(graph)$cluster[x[1]] == V(graph)$cluster[x[2]],
colors[V(graph)$cluster[x[1]]], "grey")
})
else
E(graph)$color <- NA
E(graph)$weight <- E(graph)$weight / max(E(graph)$weight)
E(graph)$width <- 2*E(graph)$weight
plot(graph, layout = layout,
edge.width = E(graph)$width,
vertex.size = 3,
vertex.label = NA, #vertex.label.cex = 0.8, vertex.label.dist = 0.7, vertex.label.color = "black",
edge.curved = 0.25,
main = main)
}Learning graphs from similarity measures
First load some toy data:
library(clusterSim)
set.seed(42)
data2moons <- shapes.two.moon(100)
data3circles <- shapes.circles3(100)Let’s define the functions to learn the following graphs (in the form of the adjacency matrix \bm{W}) based on similarity measures:
- thresholded distance graph
- Gaussian graph
- thresholded Gaussian graph
- kNN graph
- feature correlation graph
- selftuned Gaussian graph
library(dbscan)
W_thresholded <- function(X, threshold = NULL) {
# compute Euclidean distance matrix
euclid_dist <- as.matrix(dist(X, method = "euclidean", diag = FALSE, upper = FALSE))
# choose threshold automatically if needed
if (is.null(threshold)) {
max_dist <- max(euclid_dist)
min_dist <- min(euclid_dist[euclid_dist != 0])
threshold <- min_dist + 0.15*(max_dist-min_dist)
}
# generate W
W <- ifelse(euclid_dist <= threshold, 1, 0)
diag(W) <- 0
return(W)
}
W_Gaussian <- function(X, sigma2 = 0.1) {
# compute Euclidean distance matrix
euclid_dist <- as.matrix(dist(X, method = "euclidean", diag = FALSE, upper = FALSE))
# choose sigma2 automatically if needed
if (is.null(sigma2))
sigma2 <- mean(euclid_dist^2)/10
# generate W
W <- exp(-(euclid_dist^2)/(2*sigma2))
diag(W) <- 0
return(W)
}
W_thresholded_Gaussian <- function(X, sigma2 = 0.1, threshold = NULL) {
# compute Euclidean distance matrix
euclid_dist <- as.matrix(dist(X, method = "euclidean", diag = FALSE, upper = FALSE))
# choose sigma2 automatically?
if (is.null(sigma2))
sigma2 <- mean(euclid_dist^2)/10
# choose threshold automatically if needed
if (is.null(threshold)) {
max_dist <- max(euclid_dist)
min_dist <- min(euclid_dist[euclid_dist != 0])
threshold <- min_dist + 0.15*(max_dist-min_dist)
}
# generate W
W <- exp(-(euclid_dist^2)/(2*sigma2))
W[euclid_dist > threshold] <- 0 # threshold
diag(W) <- 0
return(W)
}
W_knn <- function(X, k = NULL) {
p <- nrow(X)
# choose k automatically if needed
if (is.null(k))
k <- 1 + floor(log2(p))
# compute the knn distances and so
knn_nodes <- dbscan::kNN(X, k, sort = FALSE)$id
# generate W
W <- matrix(0, p, p)
for (i in 1:p)
W[knn_nodes[i, ], i] <- W[i, knn_nodes[i, ]] <- 1
colnames(W) <- rownames(W) <- rownames(X)
return(W)
}
W_feature_correlation <- function(X, threshold = 0.5) {
W <- X %*% t(X)
W[abs(W) < threshold*mean(diag(W))] <- 0
diag(W) <- 0
return(W)
}
W_selftuned_Gaussian <- function(X, k = NULL) {
p <- nrow(X)
# choose k automatically if needed
if (is.null(k))
k <- floor(1 + log2(p)) #0.5*log2(p)
# compute Euclidean distance matrix and knn distances
knn_dist <- dbscan::kNN(X, k, sort = TRUE)$dist[, k]
euclid_dist2 <- as.matrix(dist(X, method = "euclidean", diag = FALSE, upper = FALSE))^2
#euclid_dist2 <- diag(1/knn_dist) %*% euclid_dist2 %*% diag(1/knn_dist)
euclid_dist2 <- t((1/knn_dist) * t((1/knn_dist) * euclid_dist2))
# generate W
W <- exp(-euclid_dist2/2)
diag(W) <- 0
return(W)
}Graphs learned based on similarity measures for the two-moon dataset:
old_par <- par(mfcol = c(2, 3), mar = c(0.6, 0, 0.8, 0)) # b-l-t-r
# thresholded distance graph
plot_graph(W = W_thresholded(data2moons$data),
layout = data2moons$data, clusters = data2moons$clusters,
main = "Thresholded distance 0-1 graph")
# Gaussian graph
plot_graph(W = W_Gaussian(data2moons$data),
layout = data2moons$data, clusters = data2moons$clusters,
main = "Gaussian graph")
# thresholded Gaussian graph
plot_graph(W = W_thresholded_Gaussian(data2moons$data),
layout = data2moons$data, clusters = data2moons$clusters,
main = "Thresholded Gaussian graph")
# k-NN graph
plot_graph(W = W_knn(data2moons$data),
layout = data2moons$data, clusters = data2moons$clusters,
main = "k-NN 0-1 graph")
# feature correlation graph
plot_graph(W = W_feature_correlation(data2moons$data),
layout = data2moons$data, clusters = data2moons$clusters,
main = "Feature correlation graph")
# self-tuned Gaussian graph
plot_graph(W = W_selftuned_Gaussian(data2moons$data, k = 4),
layout = data2moons$data, clusters = data2moons$clusters,
main = "Self-tuned Gaussian graph")
par(old_par)
Graphs learned based on similarity measures for the three-circle dataset:
old_par <- par(mfcol = c(2, 3), mar = c(0.6, 0, 0.8, 0)) # b-l-t-r
# thresholded distance graph
plot_graph(W = W_thresholded(data3circles$data),
layout = data3circles$data, clusters = data3circles$clusters,
main = "Thresholded distance 0-1 graph")
# Gaussian graph
plot_graph(W = W_Gaussian(data3circles$data),
layout = data3circles$data, clusters = data3circles$clusters,
main = "Gaussian graph")
# thresholded Gaussian graph
plot_graph(W = W_thresholded_Gaussian(data3circles$data),
layout = data3circles$data, clusters = data3circles$clusters,
main = "Thresholded Gaussian graph")
# k-NN graph
plot_graph(W = W_knn(data3circles$data),
layout = data3circles$data, clusters = data3circles$clusters,
main = "k-NN 0-1 graph")
# feature correlation graph
plot_graph(W = W_feature_correlation(data3circles$data),
layout = data3circles$data, clusters = data3circles$clusters,
main = "Feature correlation graph")
# self-tuned Gaussian graph
plot_graph(W = W_selftuned_Gaussian(data3circles$data, k = 4),
layout = data3circles$data, clusters = data3circles$clusters,
main = "Self-tuned Gaussian graph")
par(old_par)
Learning graphs from smooth signals
We now learn the following graphs in the form of the adjacency matrix \bm{W} based on similarity measures:
Graph with fixed degrees (Nie et al., 2016): \[ \begin{array}{ll} \underset{\bm{W}}{\textm{minimize}} & \textm{Tr}(\bm{W}\bm{Z}) + \gamma\|\bm{W}\|_\textm{F}^2\\ \textm{subject to} & \textm{diag}(\bm{W})=\bm{0}, \quad \bm{W}=\bm{W}^\T \ge \bm{0}\\ & \bm{W}\bm{1} = \bm{1}. \end{array} \]
Graph with regularized degrees (Kalofolias, 2016): \[ \begin{array}{ll} \underset{\bm{W}}{\textm{minimize}} & \textm{Tr}(\bm{W}\bm{Z}) - \alpha\bm{1}^\T\textm{log}(\bm{W}\bm{1}) + \frac{\beta}{2}\|\bm{W}\|_\textm{F}^2\\ \textm{subject to} & \textm{diag}(\bm{W})=\bm{0}, \quad \bm{W}=\bm{W}^\T \ge \bm{0}, \end{array} \]
Clustered \(k\) component graph based on the rank-contrained Laplacian algorithm (Nie et al., 2016), which is implemented in the package
spectralGraphTopology: \[ \begin{array}{ll} \underset{\bm{W}}{\textm{minimize}} & \|\bm{W} - \bm{W}_0\|_\textm{F}^2\\ \textm{subject to} & \bm{W}\bm{1} = \bm{1}, \quad \textm{diag}(\bm{W})=\bm{0}, \quad \bm{W}=\bm{W}^\T \ge \bm{0}\\ & \textm{rank}\left(\mathcal{L}(\bm{W})\right) = p - k, \end{array} \]
W_smooth <- function(X, m = NULL) {
p <- nrow(X)
# choose m automatically if needed
if (is.null(m))
m <- ceiling(0.08*p)
# compute squared Euclidean distance matrix
Z <- as.matrix(dist(X, method = "euclidean", diag = FALSE, upper = FALSE))^2
diag(Z) <- 10*max(Z)
# generate W
W <- matrix(0, p, p)
for (i in 1:p) {
z <- Z[, i]
z_sorted <- sort(z)
w_ <- pmax(z_sorted[m+1] - z, 0)
W[, i] <- w_/sum(w_)
}
W <- 1/2*(W + t(W)) # symmetrize
colnames(W) <- rownames(W) <- rownames(X)
return(W)
}
W_kalofolias_cvx <- function(X, alpha = 1e-2, beta = 1e-4) {
p <- nrow(X)
# compute squared Euclidean distance matrix
Z <- as.matrix(dist(X, method = "euclidean", diag = FALSE, upper = FALSE))^2
# optimize
W <- Symmetric(p)
objective <- .5 * sum_entries(W * Z) - alpha * sum_entries(log(sum_entries(W, axis = 1))) + .5 * beta * sum_squares(W)
#objective_value <- .5 * sum(W * Z) - alpha * sum(log(rowSums(W))) + .5 * beta * norm(W, "F")^2
constraints <- list(W >= 0,
diag(W) == 0,
sum_entries(W, axis = 1) >= 1e-3)
problem <- Problem(Minimize(objective), constraints)
result <- solve(problem, solver = "SCS")
if (result$status != "optimal") message("Warning: CVX status: ", result$status)
W <- as.matrix(result$getValue(W))
W[W < 1e-7] <- 0
return(W)
}
W_kalofolias <- function(X, m = NULL, alpha_over_beta = 1e2) {
p <- nrow(X)
# choose m automatically if needed
if (is.null(m))
m <- ceiling(0.08*p)
# compute squared Euclidean distance matrix
Z <- as.matrix(dist(X, method = "euclidean", diag = FALSE, upper = FALSE))^2
diag(Z) <- 10*max(Z)
# generate W
W <- matrix(0, p, p)
for (i in 1:p) {
z <- Z[, i]
z_sorted <- sort(z)
w_ <- alpha_over_beta * pmax(1 - z/z_sorted[m+1], 0)
W[, i] <- w_/sqrt(sum(w_))
}
W <- 1/2*(W + t(W)) # symmetrize
colnames(W) <- rownames(W) <- rownames(X)
return(W)
}Graphs learned based on smooth signals for the two-moon dataset:
library(spectralGraphTopology)
old_par <- par(mfcol = c(1, 3), mar = c(0, 0, 0.8, 0)) # b-l-t-r
plot_graph(W = W_smooth(data2moons$data, m = 20),
layout = data2moons$data, clusters = data2moons$clusters,
main = "Graph with fixed degrees")
plot_graph(W = W_kalofolias(X = data2moons$data, m = 20),
layout = data2moons$data, clusters = data2moons$clusters,
main = "Graph with regularized degrees")
plot_graph(W = cluster_k_component_graph(data2moons$data, k = 2, m = 20)$adjacency,
layout = data2moons$data, clusters = data2moons$clusters,
main = "Low-rank graph")
par(old_par)
old_par <- par(mfcol = c(1, 3), mar = c(0, 0, 0.8, 0)) # b-l-t-r
plot_graph(W = W_smooth(data3circles$data, m = 20),
layout = data3circles$data, clusters = data3circles$clusters,
main = "Graph with fixed degrees")
plot_graph(W = W_kalofolias(X = data3circles$data, m = 20),
layout = data3circles$data, clusters = data3circles$clusters,
main = "Graph with regularized degrees")
plot_graph(W = cluster_k_component_graph(data3circles$data, k = 3, m = 20)$adjacency,
layout = data3circles$data, clusters = data3circles$clusters,
main = "Low-rank graph")
par(old_par)
Learning graphs from Gaussian graphical model networks
We start with learning the Laplacian matrix \bm{L} for the following graphs:
directly from the precision matrix
sparse precision matrix from the graphical lasso (GLASSO) (Banerjee et al., 2008; Friedman et al., 2008) \[ \begin{array}{ll} \underset{\bm{\Theta}\succ\bm{0}}{\textm{maximize}} & \textm{log det}(\bm{\Theta}) - \textm{Tr}(\bm{\Theta}\bm{S}) - \rho\|\bm{\Theta}\|_{1,\textm{off}}, \end{array} \]
the Laplacian-structured GLASSO (Egilmez et al., 2017; Lake and Tenenbaum, 2010; Zhao et al., 2019), which is implemented in the package
spectralGraphTopology\[ \begin{array}{ll} \underset{\bm{L}\succeq\bm{0}}{\textm{maximize}} & \textm{log gdet}(\bm{L}) - \textm{Tr}(\bm{L}\bm{S}) - \rho\|\bm{L}\|_{1,\textm{off}}\\ \textm{subject to} & \bm{L}\bm{1}=\bm{0}, \quad L_{ij} = L_{ji} \le 0, \quad \forall i\neq j, \end{array} \]
L_precision <- function(S, threshold = 0) {
L <- ginv(S)
diag_L <- diag(L) # safe diagonal elements
L <- ifelse(abs(L) <= threshold*mean(diag_L), 0, L)
diag(L) <- diag_L # recover diagonal elements
return(L)
}
L_GLASSO <- function(S, rho = 1e-5) {
L <- glasso::glasso(S, rho)$wi
colnames(L) <- rownames(L) <- rownames(S)
return(L)
}Graphs learned based on Gaussian graphical models for the two-moon dataset:
old_par <- par(mfcol = c(1, 3), mar = c(0, 0, 0.8, 0)) # b-l-t-r
plot_graph(L = L_precision(S = data2moons$data %*% t(data2moons$data), threshold = 1e-2),
layout = data2moons$data, clusters = data2moons$clusters,
main = "Partial correlation (precision matrix) graph")
plot_graph(L = L_GLASSO(S = data2moons$data %*% t(data2moons$data), rho = 1e-4),
layout = data2moons$data, clusters = data2moons$clusters,
main = "GLASSO graph")
plot_graph(L = learn_laplacian_gle_admm(S = data2moons$data %*% t(data2moons$data), maxiter = 200)$laplacian,
layout = data2moons$data, clusters = data2moons$clusters,
main = "Laplacian-structured GLASSO graph")
par(old_par)
Graphs learned based on Gaussian graphical models for the three-circle dataset:
old_par <- par(mfcol = c(1, 3), mar = c(0, 0, 0.8, 0)) # b-l-t-r
plot_graph(L = L_precision(S = data3circles$data %*% t(data3circles$data), threshold = 1e-2),
layout = data3circles$data, clusters = data3circles$clusters,
main = "Partial correlation (precision matrix) graph")
plot_graph(L = L_GLASSO(S = data3circles$data %*% t(data3circles$data), rho = 1e-4),
layout = data3circles$data, clusters = data3circles$clusters,
main = "GLASSO graph")
plot_graph(L = learn_laplacian_gle_admm(S = data3circles$data %*% t(data3circles$data), maxiter = 200)$laplacian,
layout = data3circles$data, clusters = data3circles$clusters,
main = "Laplacian-structured GLASSO graph")
par(old_par)
Learning \(k\)-component graphs
We learn \(k\)-component graphs based on:
the smooth signal graph formulation (Nie et al., 2016): \[ \begin{array}{ll} \underset{\bm{W}}{\textm{minimize}} & \textm{Tr}(\bm{W}\bm{Z}) + \gamma\|\bm{W}\|_\textm{F}^2\\ \textm{subject to} & \textm{diag}(\bm{W})=\bm{0}, \quad \bm{W}=\bm{W}^\T \ge \bm{0}\\ & \bm{W}\bm{1} = \bm{1}. \end{array} \]
the same graph further approximated by a low-rank one (Nie et al., 2016): \[ \begin{array}{ll} \underset{\bm{W}}{\textm{minimize}} & \|\bm{W} - \bm{W}_0\|_\textm{F}^2\\ \textm{subject to} & \bm{W}\bm{1} = \bm{1}, \quad \textm{diag}(\bm{W})=\bm{0}, \quad \bm{W}=\bm{W}^\T \ge \bm{0}\\ & \textm{rank}\left(\mathcal{L}(\bm{W})\right) = p - k, \end{array} \]
the joint low-rank smooth signal graph formulation: \[ \begin{array}{ll} \underset{\bm{L}\succeq\bm{0}, \bm{F}}{\textm{minimize}} & \textm{Tr}\left(\bm{X}\bm{L}\bm{X}^\T\right) + \frac{\gamma}{2}\|\bm{L}\|_\textm{F,off}^2 + \gamma \textm{Tr}\left(\bm{F}^\T\bm{L}\bm{F}\right)\\ \textm{subject to} & \bm{L}\bm{1}=\bm{0}, \quad L_{ij} = L_{ji} \le 0, \quad \forall i\neq j\\ & \textm{diag}(\bm{L})=\bm{1}\\ & \bm{F}^\T\bm{F}=\bm{I}. \end{array} \]
Graphs learned based on smooth signals for the three-circle dataset:
library(spectralGraphTopology)
W_smooth_3circle <- W_smooth(data3circles$data)
W_approx_lowrank_3circle <- cluster_k_component_graph(data3circles$data, k = 3, m = 7)$adjacency
W_lowrank_3circle <- learn_k_component_graph(S = crossprod(t(data3circles$data)),
k = 3,
beta = 1,
fix_beta = FALSE,
verbose = TRUE,
abstol = 1e-3)$adjacency
old_par <- par(mfcol = c(1, 3), mar = c(0, 0, 0.8, 0)) # b-l-t-r
plot_graph(W = W_smooth_3circle,
layout = data3circles$data, clusters = data3circles$clusters,
main = "Smooth-signal graph")
plot_graph(W = W_approx_lowrank_3circle,
layout = data3circles$data, clusters = data3circles$clusters,
main = "Low-rank two-step approximation")
plot_graph(W = W_lowrank_3circle,
layout = data3circles$data, clusters = data3circles$clusters,
main = "Low-rank smooth-signal graph")
par(old_par)
References
Banerjee, O., El Ghaoui, L., and d’Aspremont, A. (2008). Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data. Journal of Machine Learning Research (JMLR), 9, 485–516.
Egilmez, H. E., Pavez, E., and Ortega, A. (2017). Graph learning from data under Laplacian and structural constraints. IEEE Journal of Selected Topics in Signal Processing, 11(6), 825–841.
Friedman, J., Hastie, T., and Tibshirani, R. (2008). Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3), 432–441.
Kalofolias, V. (2016). How to learn a graph from smooth signals. In Proceedings of the international conference on artificial intelligence and statistics (AISTATS), pages 920–929.
Lake, B., and Tenenbaum, J. (2010). Discovering structure by learning sparse graphs. In Proceedings of the 33rd annual cognitive science conference.
Nie, F., Wang, X., Jordan, M., and Huang, H. (2016). The constrained Laplacian rank algorithm for graph-based clustering. In Proceedings of the AAAI conference on artificial intelligence, pages 1969–1976. Phoenix, Arizona, USA.
Zhao, L., Wang, Y., Kumar, S., and Palomar, D. P. (2019). Optimization algorithms for graph Laplacian estimation via ADMM and MM. IEEE Transactions on Signal Processing, 67(16), 4231–4244.