Introduction

In this vignette, we consider optimal tensor transport (OTT), which is an extension of OT to be able to handle tensors of any order by learning possibly multiple transport plans.

Here, we reproduce the experiments in the original paper (Kerdoncuff 2022). For the details of the methodology, see the original paper.

library("otTensor")

.show_matrix <- function(mat, main = ""){
    mat_rev <- apply(mat, 2, rev)
    mat_rev <- t(mat_rev)

    row_index <- 1:ncol(mat_rev)
    col_index <- 1:nrow(mat_rev)

    # grayscale
    image(mat_rev, col = gray((0:255)/255), xaxt = "n", yaxt = "n",
        xlab = "", ylab = "", axes = FALSE, main = main)
}

OTT_1 (OT)

D <- 1 A <- 1 Is <- c(4) Ks <- c(7) f <- c(1) arrX <- array(rep(0, prod(Is)), Is) arrY <- array(rep(0, prod(Ks)), Ks)

for (i1 in 1:Is[1]) { arrX[i1] <- i1 } for (k1 in 1:Ks[1]) { arrY[k1] <- k1 }

ps <- list() for (a in 1:A) { ds <- which(f == a) d <- ds[1] length_of_p_a <- dim(arrX)[d] ps[[a]] <- rep(0.01, length_of_p_a); ps[[a]][c(1, 3)] <- 1 ps[[a]] <- ps[[a]] / sum(ps[[a]]) } qs <- list() for (a in 1:A) { ds <- which(f == a) d <- ds[1] length_of_q_a <- dim(arrY)[d] qs[[a]] <- rep(1, length_of_q_a); qs[[a]][c(2, 3)] <- 0 qs[[a]] <- qs[[a]] / sum(qs[[a]]) }

X <- as.tensor(arrX) Y <- as.tensor(arrY)

out <- OTT(X = X, Y = Y, D = D, A = A, Is = Is, Ks = Ks, f = f, ps=ps, qs=qs, num.sample=1000, loss = function (x, y) {abs(x - y)}, num.iter=200, epsilon=1e-10)

options(repr.plot.width=6, repr.plot.height=10) par(mfrow=c(3, 2)) plot(ps[[1]], type=“h”, col=“red”, ylim=c(0, 1), main=“ps[[1]]”) plot(qs[[1]], type=“h”, col=“red”, ylim=c(0, 1), main=“qs[[1]]”) .show_matrix(out$Ts[[1]], main=“Ts[[1]]”) plot(arrX, type=“h”, col=“black”, main=“arrX”) plot(arrY, type=“h”, col=“black”, main=“arrY”)

OTT_12 (Co-OT)

D <- 2 A <- 2 Is <- c(4, 5) Ks <- c(7, 8) f <- c(1, 2) arrX <- array(rep(0, prod(Is)), Is) arrY <- array(rep(0, prod(Ks)), Ks)

for (i1 in 1:Is[1]) { for (i2 in 1:Is[2]) { arrX[i1, i2] <- i1 + i2 } } for (k1 in 1:Ks[1]) { for (k2 in 1:Ks[2]) { arrY[k1, k2] <- k1 + k2 } }

ps <- list() for (a in 1:A) { ds <- which(f == a) d <- ds[1] length_of_p_a <- dim(arrX)[d] ps[[a]] <- rep(0.01, length_of_p_a); ps[[a]][c(1, 3)] <- 1 ps[[a]] <- ps[[a]] / sum(ps[[a]]) } qs <- list() for (a in 1:A) { ds <- which(f == a) d <- ds[1] length_of_q_a <- dim(arrY)[d] qs[[a]] <- rep(1, length_of_q_a); qs[[a]][c(2, 3)] <- 0 qs[[a]] <- qs[[a]] / sum(qs[[a]]) }

X <- as.tensor(arrX) Y <- as.tensor(arrY)

out <- OTT(X = X, Y = Y, D = D, A = A, Is = Is, Ks = Ks, f = f, ps=ps, qs=qs, num.sample=1000, loss = function (x, y) {abs(x - y)}, num.iter=200, epsilon=1e-10)

options(repr.plot.width=6, repr.plot.height=10) par(mfrow=c(3, 2)) plot(ps[[1]], type=“h”, col=“red”, ylim=c(0, 1), main=“ps[[1]]”) plot(qs[[1]], type=“h”, col=“red”, ylim=c(0, 1), main=“qs[[1]]”) .show_matrix(out$Ts[[1]], main=“Ts[[1]]”) .show_matrix(arrX, main=“arrX”) .show_matrix(arrY, main=“arrY”)

par(mfrow=c(3, 2)) plot(ps[[2]], type=“h”, col=“red”, ylim=c(0, 1), main=“ps[[2]]”) plot(qs[[2]], type=“h”, col=“red”, ylim=c(0, 1), main=“qs[[2]]”) .show_matrix(out$Ts[[2]], main=“Ts[[2]]”) .show_matrix(arrX, main=“arrX”) .show_matrix(arrY, main=“arrY”)

OTT_11 (GW)

D <- 2 A <- 1 Is <- c(4, 4) Ks <- c(6, 6) f <- c(1, 1) arrX <- array(rep(0, prod(Is)), Is) arrY <- array(rep(0, prod(Ks)), Ks)

for (i1 in 1:Is[1]) { for (i2 in 1:Is[2]) { arrX[i1, i2] <- i1 + i2 } } for (k1 in 1:Ks[1]) { for (k2 in 1:Ks[2]) { arrY[k1, k2] <- k1 + k2 } }

ps <- list() for (a in 1:A) { ds <- which(f == a) d <- ds[1] length_of_p_a <- dim(arrX)[d] ps[[a]] <- rep(0.01, length_of_p_a); ps[[a]][c(1, 3)] <- 1 ps[[a]] <- ps[[a]] / sum(ps[[a]]) } qs <- list() for (a in 1:A) { ds <- which(f == a) d <- ds[1] length_of_q_a <- dim(arrY)[d] qs[[a]] <- rep(1, length_of_q_a); qs[[a]][c(2, 3)] <- 0 qs[[a]] <- qs[[a]] / sum(qs[[a]]) }

X <- as.tensor(arrX) Y <- as.tensor(arrY)

out <- OTT(X = X, Y = Y, D = D, A = A, Is = Is, Ks = Ks, f = f, ps=ps, qs=qs, num.sample=1000, loss = function (x, y) {abs(x - y)}, num.iter=200, epsilon=1e-10)

options(repr.plot.width=6, repr.plot.height=10) par(mfrow=c(3, 2)) plot(ps[[1]], type=“h”, col=“red”, ylim=c(0, 1), main=“ps[[1]]”) plot(qs[[1]], type=“h”, col=“red”, ylim=c(0, 1), main=“qs[[1]]”) .show_matrix(out$Ts[[1]], main=“Ts[[1]]”) .show_matrix(arrX, main=“arrX”) .show_matrix(arrY, main=“arrY”)

OTT_111 (triplets)

D <- 3 A <- 1 Is <- c(4, 4, 4) Ks <- c(6, 6, 6) f <- c(1, 1, 1) arrX <- array(rep(0, prod(Is)), Is) arrY <- array(rep(0, prod(Ks)), Ks)

for (i1 in 1:Is[1]) { for (i2 in 1:Is[2]) { for (i3 in 1:Is[3]) { arrX[i1, i2, i3] <- i1 + i2 + i3 } } } for (k1 in 1:Ks[1]) { for (k2 in 1:Ks[2]) { for (k3 in 1:Ks[3]) { arrY[k1, k2, k3] <- k1 + k2 + k3 } } }

ps <- list() for (a in 1:A) { ds <- which(f == a) d <- ds[1] length_of_p_a <- dim(arrX)[d] ps[[a]] <- rep(0.01, length_of_p_a); ps[[a]][c(1, 3)] <- 1 ps[[a]] <- ps[[a]] / sum(ps[[a]]) } qs <- list() for (a in 1:A) { ds <- which(f == a) d <- ds[1] length_of_q_a <- dim(arrY)[d] qs[[a]] <- rep(1, length_of_q_a); qs[[a]][c(2, 3)] <- 0 qs[[a]] <- qs[[a]] / sum(qs[[a]]) }

X <- as.tensor(arrX) Y <- as.tensor(arrY)

out <- OTT(X = X, Y = Y, D = D, A = A, Is = Is, Ks = Ks, f = f, ps=ps, qs=qs, num.sample=1000, loss = function (x, y) {abs(x - y)}, num.iter=200, epsilon=1e-10)

options(repr.plot.width=6, repr.plot.height=10) par(mfrow=c(3, 2)) plot(ps[[1]], type=“h”, col=“red”, ylim=c(0, 1), main=“ps[[1]]”) plot(qs[[1]], type=“h”, col=“red”, ylim=c(0, 1), main=“qs[[1]]”) .show_matrix(out$Ts[[1]], main=“Ts[[1]]”) .show_matrix(arrX[,,1], main=“arrX[,,1]”) .show_matrix(arrY[,,1], main=“arrY[,,1]”)

OTT_123 (triCo-OT)

D <- 3 A <- 3 Is <- c(4, 5, 6) Ks <- c(7, 8, 9) f <- c(1, 2, 3) arrX <- array(rep(0, prod(Is)), Is) arrY <- array(rep(0, prod(Ks)), Ks)

for (i1 in 1:Is[1]) { for (i2 in 1:Is[2]) { for (i3 in 1:Is[3]) { arrX[i1, i2, i3] <- i1 + i2 + i3 } } } for (k1 in 1:Ks[1]) { for (k2 in 1:Ks[2]) { for (k3 in 1:Ks[3]) { arrY[k1, k2, k3] <- k1 + k2 + k3 } } }

ps <- list() for (a in 1:A) { ds <- which(f == a) d <- ds[1] length_of_p_a <- dim(arrX)[d] ps[[a]] <- rep(0.01, length_of_p_a); ps[[a]][c(1, 3)] <- 1 ps[[a]] <- ps[[a]] / sum(ps[[a]]) } qs <- list() for (a in 1:A) { ds <- which(f == a) d <- ds[1] length_of_q_a <- dim(arrY)[d] qs[[a]] <- rep(1, length_of_q_a); qs[[a]][c(2, 3)] <- 0 qs[[a]] <- qs[[a]] / sum(qs[[a]]) }

X <- as.tensor(arrX) Y <- as.tensor(arrY)

out <- OTT(X = X, Y = Y, D = D, A = A, Is = Is, Ks = Ks, f = f, ps=ps, qs=qs, num.sample=1000, loss = function (x, y) {abs(x - y)}, num.iter=200, epsilon=1e-10)

options(repr.plot.width=6, repr.plot.height=10) par(mfrow=c(3, 2)) plot(ps[[1]], type=“h”, col=“red”, ylim=c(0, 1), main=“ps[[1]]”) plot(qs[[1]], type=“h”, col=“red”, ylim=c(0, 1), main=“qs[[1]]”) .show_matrix(out$Ts[[1]], main=“Ts[[1]]”) .show_matrix(arrX[,,1], main=“arrX[,,1]”) .show_matrix(arrY[,,1], main=“arrY[,,1]”)

par(mfrow=c(3, 2)) plot(ps[[2]], type=“h”, col=“red”, ylim=c(0, 1), main=“ps[[2]]”) plot(qs[[2]], type=“h”, col=“red”, ylim=c(0, 1), main=“qs[[2]]”) .show_matrix(out$Ts[[2]], main=“Ts[[2]]”) .show_matrix(arrX[,,2], main=“arrX[,,2]”) .show_matrix(arrY[,,2], main=“arrY[,,2]”)

par(mfrow=c(3, 2)) plot(ps[[3]], type=“h”, col=“red”, ylim=c(0, 1), main=“ps[[3]]”) plot(qs[[3]], type=“h”, col=“red”, ylim=c(0, 1), main=“qs[[3]]”) .show_matrix(out$Ts[[3]], main=“Ts[[3]]”) .show_matrix(arrX[,,3], main=“arrX[,,3]”) .show_matrix(arrY[,,3], main=“arrY[,,3]”)

OTT_112 (GW Collection)

D <- 3 A <- 2 Is <- c(4, 4, 5) Ks <- c(6, 6, 7) f <- c(1, 1, 2) arrX <- array(rep(0, prod(Is)), Is) arrY <- array(rep(0, prod(Ks)), Ks)

for (i1 in 1:Is[1]) { for (i2 in 1:Is[2]) { for (i3 in 1:Is[3]) { arrX[i1, i2, i3] <- i1 + i2 + i3 } } } for (k1 in 1:Ks[1]) { for (k2 in 1:Ks[2]) { for (k3 in 1:Ks[3]) { arrY[k1, k2, k3] <- k1 + k2 + k3 } } }

ps <- list() for (a in 1:A) { ds <- which(f == a) d <- ds[1] length_of_p_a <- dim(arrX)[d] ps[[a]] <- rep(0.01, length_of_p_a); ps[[a]][c(1, 3)] <- 1 ps[[a]] <- ps[[a]] / sum(ps[[a]]) } qs <- list() for (a in 1:A) { ds <- which(f == a) d <- ds[1] length_of_q_a <- dim(arrY)[d] qs[[a]] <- rep(1, length_of_q_a); qs[[a]][c(2, 3)] <- 0 qs[[a]] <- qs[[a]] / sum(qs[[a]]) }

X <- as.tensor(arrX) Y <- as.tensor(arrY)

out <- OTT(X = X, Y = Y, D = D, A = A, Is = Is, Ks = Ks, f = f, ps=ps, qs=qs, num.sample=1000, loss = function (x, y) {abs(x - y)}, num.iter=200, epsilon=1e-10)

options(repr.plot.width=6, repr.plot.height=10) par(mfrow=c(3, 2)) plot(ps[[1]], type=“h”, col=“red”, ylim=c(0, 1), main=“ps[[1]]”) plot(qs[[1]], type=“h”, col=“red”, ylim=c(0, 1), main=“qs[[1]]”) .show_matrix(out$Ts[[1]], main=“Ts[[1]]”) .show_matrix(arrX[,,1], main=“arrX[,,1]”) .show_matrix(arrY[,,1], main=“arrY[,,1]”)

par(mfrow=c(3, 2)) plot(ps[[2]], type=“h”, col=“red”, ylim=c(0, 1), main=“ps[[2]]”) plot(qs[[2]], type=“h”, col=“red”, ylim=c(0, 1), main=“qs[[2]]”) .show_matrix(out$Ts[[2]], main=“Ts[[2]]”) .show_matrix(arrX[,,2], main=“arrX[,,2]”) .show_matrix(arrY[,,2], main=“arrY[,,2]”)