channel_capacity <- function(n, m, P, sum_x = 1) {
## x is probability distribution of the input signal X(t)
x <- Variable(n)
## y is the probability distribution of the output signal Y(t)
y <- P %*% x
## Compute c_j = sum_i p_{ij} log2(p_{ij})
## Use the convention 0*log(0) = 0 by handling zeros
P_log <- ifelse(P > 0, log2(P), 0)
c_vec <- colSums(P * P_log)
## Mutual information: I = c'x + sum(entr(y)) / log(2)
I <- t(c_vec) %*% x + sum_entries(entr(y)) / log(2)
## Channel capacity is maximized by maximising the mutual information
obj <- Maximize(I)
constraints <- list(sum(x) == sum_x, x >= 0)
## Form and solve problem
prob <- Problem(obj, constraints)
result <- psolve(prob)
list(status = status(prob),
capacity = result,
x = value(x))
}Channel Capacity
Introduction
This example is adapted from Boyd and Vandenberghe, Convex Optimization, exercise 4.57, pages 207–208.
Convex optimization can be used to find the channel capacity
In a discrete memoryless channel, the relation between the input and the output is given by the transition probability:
These transition probabilities form the channel transition matrix
Assume that
From Shannon’s theory, the channel capacity is given by the maximum possible mutual information
where
and
Given that
Due to the entropy function in CVXR, this can be written quite easily in DCP form.
Channel Capacity Function
We write a function that computes the channel capacity for a given transition matrix
Example
In this example we consider a communication channel with two possible inputs and outputs, so
Note that the columns of
n <- 2
m <- 2
P <- matrix(c(0.75, 0.25,
0.25, 0.75), nrow = m, ncol = n, byrow = TRUE)
result <- channel_capacity(n, m, P)
cat("Problem status:", result$status, "\n")
cat(sprintf("Optimal value of C = %.4g\n", result$capacity))
cat("Optimal variable x =\n")
print(round(result$x, 4))Problem status: optimal
Optimal value of C = 0.1887
Optimal variable x =
[,1]
[1,] 0.5
[2,] 0.5
The optimal input distribution is uniform (
A Larger Example
We can also consider a non-symmetric channel with
n2 <- 4
m2 <- 3
set.seed(42)
## Generate a random transition matrix (columns sum to 1)
P2 <- matrix(runif(m2 * n2), nrow = m2, ncol = n2)
P2 <- sweep(P2, 2, colSums(P2), "/")
result2 <- channel_capacity(n2, m2, P2)
cat("Problem status:", result2$status, "\n")
cat(sprintf("Optimal value of C = %.4g\n", result2$capacity))
cat("Optimal variable x =\n")
print(round(result2$x, 4))Problem status: optimal
Optimal value of C = 0.152
Optimal variable x =
[,1]
[1,] 0.4805
[2,] 0.0000
[3,] 0.5195
[4,] 0.0000
Session Info
R version 4.5.2 (2025-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.3
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/Los_Angeles
tzcode source: internal
attached base packages:
[1] stats graphics grDevices datasets utils methods base
other attached packages:
[1] CVXR_1.8.0.9207
loaded via a namespace (and not attached):
[1] slam_0.1-55 cli_3.6.5 knitr_1.51 ECOSolveR_0.6.1
[5] rlang_1.1.7 xfun_0.56 clarabel_0.11.2 otel_0.2.0
[9] gurobi_13.0-1 Rglpk_0.6-5.1 highs_1.12.0-3 cccp_0.3-3
[13] scs_3.2.7 S7_0.2.1 jsonlite_2.0.0 Rcplex_0.3-8
[17] backports_1.5.0 Rmosek_11.1.1 htmltools_0.5.9 gmp_0.7-5.1
[21] rmarkdown_2.30 piqp_0.6.2 grid_4.5.2 evaluate_1.0.5
[25] fastmap_1.2.0 yaml_2.3.12 compiler_4.5.2 codetools_0.2-20
[29] htmlwidgets_1.6.4 Rcpp_1.1.1 osqp_1.0.0 lattice_0.22-9
[33] digest_0.6.39 checkmate_2.3.4 Matrix_1.7-4 tools_4.5.2
References
- Boyd, S. and Vandenberghe, L. (2004). Convex Optimization. Cambridge University Press. Exercise 4.57, pages 207–208.
- Shannon, C.E. (1948). A Mathematical Theory of Communication. Bell System Technical Journal, 27, 379–423, 623–656.