R Code for Portfolio Optimization
Chapter 4 – Financial Data: Graphs

Daniel P. Palomar (2025). Portfolio Optimization: Theory and Application. Cambridge University Press.

Last update: April 14, 2025

---

Loading packages

The following packages are used in the examples:

# 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

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)

Graphs learned based on smooth signals for the three-circle dataset.

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.