-## Citation
-If you made use of this software please consider citing:
-- ## Citation
+## Citation
Please cite:
- [A Javaheri](https://javaheriamirhossein.github.io/), [JVM Cardoso](https://mirca.github.io) and
From c9c9c403f2c40f7a2a8fa19443b65b122cbb060e Mon Sep 17 00:00:00 2001
From: Amirhossein Javaheri
<46240571+javaheriamirhossein@users.noreply.github.com>
Date: Mon, 26 May 2025 00:54:28 +0800
Subject: [PATCH 07/25] Create learn_kcomp_heavytail_graph_balanced.R
---
R/learn_kcomp_heavytail_graph_balanced.R | 230 +++++++++++++++++++++++
1 file changed, 230 insertions(+)
create mode 100644 R/learn_kcomp_heavytail_graph_balanced.R
diff --git a/R/learn_kcomp_heavytail_graph_balanced.R b/R/learn_kcomp_heavytail_graph_balanced.R
new file mode 100644
index 0000000..7697dd9
--- /dev/null
+++ b/R/learn_kcomp_heavytail_graph_balanced.R
@@ -0,0 +1,230 @@
+library(spectralGraphTopology)
+
+#' @title Laplacian matrix of a balanced k-component graph with heavy-tailed data
+#'
+#' Computes the Laplacian matrix of a balanced graph on the basis of an observed data matrix,
+#' where we assume the data to be Student-t distributed.
+#'
+#' @param X an n x p data matrix, where n is the number of observations and p is
+#' the number of nodes in the graph.
+#' @param k the number of components of the graph.
+#' @param heavy_type a string which selects the statistical distribution of the data .
+#' Valid values are "gaussian" or "student".
+#' @param nu the degrees of freedom of the Student-t distribution.
+#' Must be a real number greater than 2.
+#' @param w0 initial vector of graph weights. Either a vector of length p(p-1)/2 or
+#' a string indicating the method to compute an initial value.
+#' @param beta hyperparameter that controls the regularization to obtain a
+#' k-component graph
+#' @param update_beta whether to update beta during the optimization.
+#' @param alpha spare regularization term parameter
+#' @param t upperbound scale parameter
+#' @param d the nodes' degrees. Either a vector or a single value.
+#' @param rho ADMM hyperparameter.
+#' @param update_rho whether or not to update rho during the optimization.
+#' @param maxiter maximum number of iterations.
+#' @param reltol relative tolerance as a convergence criteria.
+#' @param verbose whether or not to show a progress bar during the iterations.
+#' @export
+#' @import spectralGraphTopology
+learn_kcomp_heavytail_graph_balanced <- function(X,
+ alpha = 0.1,
+ k = 1,
+ heavy_type = "gaussian",
+ nu = NULL,
+ w0 = "naive",
+ d = 1,
+ t = 1.2,
+ beta = 1e-8,
+ update_beta = TRUE,
+ early_stopping = FALSE,
+ rho = 1,
+ update_rho = FALSE,
+ maxiter = 10000,
+ reltol = 1e-5,
+ verbose = TRUE,
+ record_objective = FALSE) {
+
+ X <- scale(as.matrix(X))
+ # number of nodes
+ p <- ncol(X)
+
+ t <- t*d/sqrt(p-1)
+ # number of observations
+ n <- nrow(X)
+ LstarSq <- vector(mode = "list", length = n)
+ for (i in 1:n)
+ LstarSq[[i]] <- Lstar(X[i, ] %*% t(X[i, ])) / n
+ # w-initialization
+ if (assertthat::is.string(w0)) {
+ w <- spectralGraphTopology:::w_init(w0, MASS::ginv(cor(X)))
+ A0 <- A(w)
+ A0 <- A0 / rowSums(A0)
+ w <- spectralGraphTopology:::Ainv(A0)
+ }
+ else {
+ w <-w0
+ }
+ # Theta-initilization
+ Lw <- L(w)
+ Aw <- A(w)
+ Theta <- Lw
+ U <- eigen(Lw, symmetric = TRUE)$vectors[, (p - k + 1):p]
+ Y <- matrix(0, p, p)
+ y <- rep(0, p)
+ # Z-initialization
+ Lambda <- matrix(0, p, p)
+ Z <- A0
+
+ # ADMM constants
+ mu <- 2
+ tau <- 2
+ # residual vectors
+ primal_lap_residual <- c()
+ primal_deg_residual <- c()
+ dual_residual <- c()
+ # augmented lagrangian vector
+ lagrangian <- c()
+ beta_seq <- c()
+ if (verbose)
+ pb <- progress::progress_bar$new(format = "<:bar> :current/:total eta: :eta",
+ total = maxiter, clear = FALSE, width = 80)
+ elapsed_time <- c()
+ start_time <- proc.time()[3]
+ for (i in 1:maxiter) {
+
+ for (j in 1:1){
+ # update w
+ LstarLw <- Lstar(Lw)
+ DstarDw <- Dstar(diag(Lw))
+ LstarSweighted <- rep(0, .5*p*(p-1))
+ if (heavy_type == "student") {
+ for (q in 1:n)
+ LstarSweighted <- LstarSweighted + LstarSq[[q]] * compute_student_weights(w, LstarSq[[q]], p, nu)
+ } else if (heavy_type == "gaussian") {
+ for (q in 1:n)
+ LstarSweighted <- LstarSweighted + LstarSq[[q]]
+ }
+ grad <- LstarSweighted + Lstar(beta * crossprod(t(U)) - Y - rho * Theta) + Dstar(y - rho * d) + rho * (LstarLw + DstarDw)
+ grad <- grad + Astar( rho * (Aw-Z) + Lambda )
+ eta <- 1 / (2*rho * (2*p - 1)+ 2*rho)
+ wi <- w - eta * grad
+ thr <- sqrt(2*alpha *eta )
+ wi[wi< thr] <- 0
+ Lwi <- L(wi)
+ Awi <- A(wi)
+ }
+
+ # Update Z
+ Z <- Awi + Lambda/rho
+ for (i in 1:p){
+ norm_i <- norm(Z[i,], type="2")
+ if (norm_i>t) {
+ Z[i,] <- Z[i,]/norm_i *t
+ }
+ }
+
+
+ # update U
+ U <- eigen(Lwi, symmetric = TRUE)$vectors[, (p - k + 1):p]
+
+ # update Theta
+ eig <- eigen(rho * Lwi - Y, symmetric = TRUE)
+ V <- eig$vectors[,1:(p-k)]
+ gamma <- eig$values[1:(p-k)]
+ Thetai <- V %*% diag((gamma + sqrt(gamma^2 + 4 * rho)) / (2 * rho)) %*% t(V)
+
+ # update Y
+ R1 <- Thetai - Lwi
+ Y <- Y + rho * R1
+
+ # update y
+ R2 <- diag(Lwi) - d
+ y <- y + rho * R2
+
+ R3 = Awi - Z
+ Lambda <- Lambda + rho * R3
+ # compute primal, dual residuals, & lagrangian
+ primal_lap_residual <- c(primal_lap_residual, norm(R1, "F"))
+ primal_deg_residual <- c(primal_deg_residual, norm(R2, "2"))
+ dual_residual <- c(dual_residual, rho*norm(Lstar(Theta - Thetai), "2"))
+ lagrangian <- c(lagrangian, compute_augmented_lagrangian_kcomp_ht_balanced(wi, LstarSq, Thetai, U, Y, y, d, heavy_type, n, p, k, rho, beta, nu, alpha, Z, Lambda))
+
+ # update rho
+ if (update_rho) {
+ eig_vals <- spectralGraphTopology:::eigval_sym(Theta)
+ n_zero_eigenvalues <- sum(eig_vals < 1e-9)
+ if (k < n_zero_eigenvalues)
+ rho <- .5 * rho
+ else if (k > n_zero_eigenvalues)
+ rho <- 2 * rho
+ else {
+ if (early_stopping) {
+ has_converged <- TRUE
+ break
+ }
+ }
+ }
+ if (update_beta) {
+ eig_vals <- spectralGraphTopology:::eigval_sym(L(wi))
+ n_zero_eigenvalues <- sum(eig_vals < 1e-9)
+ if (k < n_zero_eigenvalues)
+ beta <- .5 * beta
+ else if (k > n_zero_eigenvalues)
+ beta <- 2 * beta
+ else {
+ if (early_stopping) {
+ has_converged <- TRUE
+ break
+ }
+ }
+ beta_seq <- c(beta_seq, beta)
+ }
+ if (verbose)
+ pb$tick()
+
+ elapsed_time <- c(elapsed_time, proc.time()[3] - start_time)
+ has_converged <- (norm(Lwi - Lw, 'F') / norm(Lw, 'F') < reltol) && (i > 1)
+ if (has_converged)
+ break
+ w <- wi
+ Lw <- Lwi
+ Aw <- Awi
+
+ Theta <- Thetai
+ }
+ results <- list(laplacian = L(wi), adjacency = A(wi), theta = Thetai, maxiter = i,
+ convergence = has_converged, beta_seq = beta_seq,
+ primal_lap_residual = primal_lap_residual,
+ primal_deg_residual = primal_deg_residual,
+ dual_residual = dual_residual,
+ lagrangian = lagrangian,
+ elapsed_time = elapsed_time)
+ return(results)
+}
+
+compute_augmented_lagrangian_kcomp_ht_balanced <- function(w, LstarSq, Theta, U, Y, y, d, heavy_type, n, p, k, rho, beta, nu, alpha, Z, Lambda) {
+ eig <- eigen(Theta, symmetric = TRUE, only.values = TRUE)$values[1:(p-k)]
+ Lw <- L(w)
+ Aw <- A(w)
+ Dw <- diag(Lw)
+ u_func <- 0
+ if (heavy_type == "student") {
+ for (q in 1:n)
+ u_func <- u_func + (p + nu) * log(1 + n * sum(w * LstarSq[[q]]) / nu)
+ } else if (heavy_type == "gaussian"){
+ for (q in 1:n)
+ u_func <- u_func + sum(n * w * LstarSq[[q]])
+ }
+ u_func <- u_func / n
+ return(u_func - sum(log(eig)) + sum(y * (Dw - d)) + sum(diag(Y %*% (Theta - Lw)))
+ + .5 * rho * (norm(Dw - d, "2")^2 + norm(Lw - Theta, "F")^2) + beta * sum(w * Lstar(crossprod(t(U))))
+ + .5 * rho * norm(Aw - Z, "2")^2 + sum(diag(Lambda %*% (Aw - Z)))
+ + alpha * sum(w>0) )
+}
+
+hardThresh <- function(v, thr){
+
+ return( v * (abs(v) > thr) )
+
+}
From 908addd721d880ed04255de0d822a76df4c79cba Mon Sep 17 00:00:00 2001
From: Amirhossein Javaheri
<46240571+javaheriamirhossein@users.noreply.github.com>
Date: Mon, 26 May 2025 00:55:59 +0800
Subject: [PATCH 08/25] Delete README.html
---
README.html | 568 ----------------------------------------------------
1 file changed, 568 deletions(-)
delete mode 100644 README.html
diff --git a/README.html b/README.html
deleted file mode 100644
index d0d7438..0000000
--- a/README.html
+++ /dev/null
@@ -1,568 +0,0 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-fingraph README
- -
-
-
-
-
-
-fingraph
-
This repo contains ADMM implementations to estimate weighted undirected graphs (Markov random fields) under Student-t assumptions with applications to financial markets.
-
-
-Installation
-fingraph depends on the development version of spectralGraphTopology, which can be installed as:
-> devtools::install_github("convexfi/spectralGraphTopology")The stable version of fingraph can be installed directly from CRAN:
-> install.packages("fingraph")
-
-Microsoft Windows
-On MS Windows environments, make sure to install the most recent version of Rtools.
-
-Usage
-
-
-Learning a graph of cryptocurrencies
-library(igraph)
-library(fingraph)
-library(fitHeavyTail)
-library(xts)
-set.seed(123)
-
-# load crypto prices into an xts table
-crypto_prices <- readRDS("examples/crypto/crypto-prices.rds")
-colnames(crypto_prices)
-#> [1] "BTC" "ETH" "USDT" "BNB" "USDC" "XRP"
-#> [7] "ADA" "HEX" "DOGE" "SOL" "MATIC" "DOT"
-#> [13] "TRX" "LTC" "BUSD" "SHIB" "AVAX" "DAI"
-#> [19] "LEO" "LINK" "ATOM" "UNI7083" "XMR" "OKB"
-#> [25] "ETC" "TON11419" "XLM" "BCH" "ICP" "CNX"
-#> [31] "TUSD" "FIL" "HBAR" "CRO" "LDO" "NEAR"
-#> [37] "VET" "QNT" "ALGO" "USDP" "FTM" "GRT6719"
-
-# compute log-returns
-log_returns <- diff(log(crypto_prices), na.pad = FALSE)
-
-# estimate a weighted, undirected graph (markov random field)
-graph_mrf <- learn_kcomp_heavytail_graph(scale(log_returns),
- k = 8,
- heavy_type = "student",
- nu = fit_mvt(scale(log_returns))$nu,
- verbose = FALSE)
-
-# plot network
-net <- graph_from_adjacency_matrix(graph_mrf$adjacency,
- mode = "undirected",
- weighted = TRUE)
-cfg <- cluster_fast_greedy(as.undirected(net))
-la_kcomp <- layout_nicely(net)
-V(net)$label.cex = 1
-plot(cfg, net, vertex.label = colnames(crypto_prices),
- layout = la_kcomp,
- vertex.size = 4.5,
- col = "black",
- edge.color = c("#686de0"),
- vertex.label.family = "Helvetica",
- vertex.label.color = "black",
- vertex.label.dist = 1.25,
- vertex.shape = "circle",
- edge.width = 20*E(net)$weight,
- edge.curved = 0.1)
-
-Learning a network of S&P500 stocks
-library(xts)
-library(igraph)
-library(fingraph)
-library(fitHeavyTail)
-library(readr)
-set.seed(123)
-
-# load table w/ stocks and their sectors
-SP500 <- read_csv("examples/stocks/SP500-sectors.csv")
-
-# load stock prices into an xts table
-stock_prices <- readRDS("examples/stocks/stock-data-2014-2018.rds")
-colnames(stock_prices)
-#> [1] "AEE" "AEP" "AES" "AIV" "AMT" "ARE" "ATO" "ATVI" "AVB"
-#> [10] "AWK" "BXP" "CBRE" "CCI" "CHTR" "CMCSA" "CMS" "CNP" "CTL"
-#> [19] "D" "DIS" "DISCA" "DISCK" "DISH" "DLR" "DRE" "DTE" "DUK"
-#> [28] "EA" "ED" "EIX" "EQIX" "EQR" "ES" "ESS" "ETR" "EVRG"
-#> [37] "EXC" "EXR" "FB" "FE" "FRT" "GOOG" "GOOGL" "HST" "IPG"
-#> [46] "IRM" "KIM" "LNT" "LYV" "MAA" "NEE" "NFLX" "NI" "NRG"
-#> [55] "NWS" "NWSA" "O" "OMC" "PEAK" "PEG" "PLD" "PNW" "PPL"
-#> [64] "PSA" "REG" "SBAC" "SLG" "SO" "SPG" "SRE" "T" "TMUS"
-#> [73] "TTWO" "TWTR" "UDR" "VNO" "VTR" "VZ" "WEC" "WELL" "WY"
-#> [82] "XEL"
-
-# compute log-returns
-log_returns <- diff(log(stock_prices), na.pad = FALSE)
-
-# estimate a weighted, undirected graph (markov random field)
-graph_mrf <- learn_kcomp_heavytail_graph(scale(log_returns),
- rho = 10,
- k = 3,
- heavy_type = "student",
- nu = fit_mvt(scale(log_returns))$nu,
- verbose = FALSE)
-#> Warning in tclass.xts(x): index does not have a 'tclass' attribute
-
-#> Warning in tclass.xts(x): index does not have a 'tclass' attribute
-
-# map stock names and sectors
-stock_sectors <- c(SP500$GICS.Sector[SP500$Symbol %in% colnames(stock_prices)])
-stock_sectors_index <- as.numeric(as.factor(stock_sectors))
-
-# plot network
-net <- graph_from_adjacency_matrix(graph_mrf$adjacency,
- mode = "undirected",
- weighted = TRUE)
-la_kcomp <- layout_nicely(net)
-V(net)$label.cex = 1
-colors <- c("#FD7272", "#55E6C1", "#25CCF7")
-V(net)$color <- colors[stock_sectors_index]
-V(net)$type <- stock_sectors_index
-V(net)$cluster <- stock_sectors_index
-E(net)$color <- apply(as.data.frame(get.edgelist(net)), 1,
- function(x) ifelse(V(net)$cluster[x[1]] == V(net)$cluster[x[2]],
- colors[V(net)$cluster[x[1]]], 'grey'))
-plot(net, vertex.label = colnames(stock_prices),
- layout = la_kcomp,
- vertex.size = 4.5,
- vertex.label.family = "Helvetica",
- vertex.label.dist = 1.25,
- vertex.label.color = "black",
- vertex.shape = "circle",
- edge.width = 20*E(net)$weight,
- edge.curved = 0.1)
-
-
-Citation
-If you made use of this software please consider citing:
--
-
- Cardoso JVM, Ying J, Palomar DP (2021). Graphical Models in Heavy-Tailed Markets. Advances in Neural Information Processing Systems (NeurIPS’21). -
