Sampling sorted vectors whose multiplication with a design matrix is also sorted

Question

I want to sample the elements of vectors a, b, c of length C in a way that both the elements of each vector a, b, c as well as the elements of D %*% A1 are sorted. Here, A1 is the matrix that results from row-binding the vectors a, b and c. D is a design matrix (see code).

I cannot use brute force until I find a solution because the real problem involves more than three vectors (a, b, c, d, ...) and potentially higher values of C. Therefore, brute force takes too long to find a solution.

Minimal example:

# Size of the vectors
C <- 3

# Sample sorted vectors
a <- sort(runif(C, 0, 1))
b <- sort(runif(C, 0, 1))
c <- sort(runif(C, 0, 1))

# Row-bind the vectors a, b, c to matrix A1
A1 <- rbind(a, b, c)

# Create a design matrix D
D <- matrix(c(1, 1, 0, -1, 0, 1, 0, -1, -1), nrow = 3)

# Multiply
A2 <- D %*% A1

# Is the result as desired?
all(t(apply(A2, 1, sort)) == A2)

General Example

Code to create the 'larger cases' mentioned in the bounty.

# Size of the vectors
C <- 4

# Number of vectors
nvec <- 100

# Create matrix A1 with sorted rows
A1 <- vector("list", nvec)
for (i in 1:nvec) {
  A1[[i]] <- sort(runif(C, 0, 1))
}
A1 <- do.call(rbind, A1)

# Create a design matrix D
DPosition <- combn(nvec, 2)
D <- matrix(0, ncol(DPosition), nvec)
for (i in 1:nrow(D)) {
  D[i, DPosition[, i]] <- c(1, -1)
}

# Multiply
A2 <- D %*% A1

# Is the result as desired?
all(t(apply(A2, 1, sort)) == A2)

Do you have an idea how I could achieve this (in R)?

EDIT: Added a more general example to clarify how the problem scales.

Answer 1

Update

Here is an update for f, such that it returns a list with both A1 and D with given C and nr, and A2 is not negligibe

f <- function(C, nr) {
  res <- matrix(nrow = nr, ncol = C)
  res[1, ] <- sort(runif(C, 1 - 1 / nr, 1))
  for (k in 2:nr) {
    d <- diff(res[k - 1, ])
    v <- c(runif(1, 0, d[1]), d)
    res[k, ] <- 1 - k / nr + cumsum(v - runif(1, 0, v[1]))
  }

  DPosition <- combn(nr, 2)
  D <- matrix(0, ncol(DPosition), nr)
  for (i in 1:nrow(D)) {
    D[i, DPosition[, i]] <- c(1, -1)
  }

  list(A1 = res, D = D)
}

for example

> C <- 4

> nr <- 6

> (r <- f(C, nr))
$A1
             [,1]         [,2]       [,3]       [,4]
[1,] 0.9074367087 0.9136961138 0.92941219 0.94173666
[2,] 0.6696354939 0.6726631943 0.68514756 0.69424033
[3,] 0.5017587194 0.5035297048 0.51475736 0.52259342
[4,] 0.3333633168 0.3347584666 0.34561028 0.35307051
[5,] 0.1669981771 0.1675170401 0.17749257 0.18407651
[6,] 0.0002043442 0.0006086951 0.01046971 0.01693914

$D
      [,1] [,2] [,3] [,4] [,5] [,6]
 [1,]    1   -1    0    0    0    0
 [2,]    1    0   -1    0    0    0
 [3,]    1    0    0   -1    0    0
 [4,]    1    0    0    0   -1    0
 [5,]    1    0    0    0    0   -1
 [6,]    0    1   -1    0    0    0
 [7,]    0    1    0   -1    0    0
 [8,]    0    1    0    0   -1    0
 [9,]    0    1    0    0    0   -1
[10,]    0    0    1   -1    0    0
[11,]    0    0    1    0   -1    0
[12,]    0    0    1    0    0   -1
[13,]    0    0    0    1   -1    0
[14,]    0    0    0    1    0   -1
[15,]    0    0    0    0    1   -1


> all(apply(r$A1, 1, Negate(is.unsorted)))
[1] TRUE

> (A2 <- with(r, D %*% A1))
           [,1]      [,2]      [,3]      [,4]
 [1,] 0.2378012 0.2410329 0.2442646 0.2474963
 [2,] 0.4056780 0.4101664 0.4146548 0.4191432
 [3,] 0.5740734 0.5789376 0.5838019 0.5886662
 [4,] 0.7404385 0.7461791 0.7519196 0.7576602
 [5,] 0.9072324 0.9130874 0.9189425 0.9247975
 [6,] 0.1678768 0.1691335 0.1703902 0.1716469
 [7,] 0.3362722 0.3379047 0.3395373 0.3411698
 [8,] 0.5026373 0.5051462 0.5076550 0.5101638
 [9,] 0.6694311 0.6720545 0.6746778 0.6773012
[10,] 0.1683954 0.1687712 0.1691471 0.1695229
[11,] 0.3347605 0.3360127 0.3372648 0.3385169
[12,] 0.5015544 0.5029210 0.5042876 0.5056543
[13,] 0.1663651 0.1672414 0.1681177 0.1689940
[14,] 0.3331590 0.3341498 0.3351406 0.3361314
[15,] 0.1667938 0.1669083 0.1670229 0.1671374

> all(apply(A2, 1, Negate(is.unsorted)))
[1] TRUE

---

Previous Answer

I guess you need some artificial interventions before sampling in a really "random" manner. Here is an example you can try

f <- function(C, nr) {
  res <- matrix(nrow = nr, ncol = C)
  res[1, ] <- sort(runif(C, 0, 1))
  for (k in 2:nr) {
    d <- diff(res[k - 1, ])
    v <- c(runif(1, 0, d[1]), d)
    res[k, ] <- cumsum(v - runif(1, 0, v[1]))
  }
  res
}

such that

> (A1 <- f(3, 3))
          [,1]      [,2]      [,3]
[1,] 0.3800352 0.7774452 0.9919061
[2,] 0.1579321 0.5128165 0.6847518
[3,] 0.0979778 0.4387943 0.5966617

> all(t(apply(A1, 1, sort)) == A1)
[1] TRUE

> (A2 <- D %*% A1)
          [,1]       [,2]       [,3]
[1,] 0.2221031 0.26462868 0.30715428
[2,] 0.2820574 0.33865089 0.39524441
[3,] 0.0599543 0.07402221 0.08809012

> all(t(apply(A2, 1, sort)) == A2)
[1] TRUE

Answer 2

Not a complete answer, but too long for a comment.

The goal as I understand it is to sample n vectors v₁, v₂, …, v_n ∈ (0, 1)^d having sorted entries such that, for all 1 ≤ i < j ≤ n, the vector v_i − v_j has sorted entries.

Define the difference map Δ : R^d → R^d−1 as

(x₁, x₂, …, x_d) ↦ (x₂ − x₁, x₃ − x₂, …, x_d − x_d−1).

A vector x ∈ R^d is sorted if and only if Δx ≥ 0. In particular, the vector v_i − v_j is sorted if and only if Δ(v_i − v_j) ≥ 0. Since Δ is linear, this condition is equivalent to Δv_i ≥ Δv_j, or, expanding the universal quantification and including the proposition that v_n is sorted,

Δv₁ ≥ Δv₂ ≥ … ≥ Δv_n ≥ 0.

Thomas’s (first) proposal is

• Sample a sorted vector v₁ ∈ (0, 1)^d uniformly at random (i.e., draw a random vector and sort it, since the degenerate cases are theoretically measure zero).
• For i from 2 to n, sample y_i ∈ (0, (Δv_i−1)₁) × (0, (Δv_i−1)₂) × … × (0, (Δv_i−1)_d−1) uniformly at random, then sample v_i ∈ (0, 1)^d such that Δv_i = y_i uniformly at random (i.e., sample the first element from the appropriate range and compute the cumulative sums).

I’ve said uniformly at random a lot, but the overall result is not uniform random. The entries of Δv_n are disproportionately likely to be small compared to brute force rejection sampling.

I suspect that it would be an improvement to sample, independently for each j ∈ {1, 2, …, d}, uniform random (y_n)_j ≤ (y_n−1)_j ≤ … ≤ (y₂)_j ≤ (Δv₁)_j. But this is still not right, since (intuitively) the more we shrink the differences, the more options we have for the first elements.

---

I can’t seem to find a nice way to sample v₁, v₂, …, v_n uniformly at random. I did find a way to sample Δv₁, Δv₂, …, Δv_n, by grabbing n(d−1) + 1 exponentially distributed random values, normalizing out the extra degree of freedom, and computing row and column cumulative sums (cf. sampling from the Dirichlet distribution). Since you don’t care about the sampler being uniform, we can just sample uniformly from the v₁, v₂, …, v_n that yield those deltas. In (probably bad) R:

C <- 4
nvec <- 20
x <- rexp(nvec * (C - 1))
x <- x/(sum(x) + rexp(1))
A1 <- matrix(0, nvec, C)
A1[1:nvec, 2:C] <- matrix(x, nvec, C - 1)
A1 <- apply(A1, 2, cumsum)
A1 <- apply(A1, 2, rev)
A1 <- t(apply(A1, 1, cumsum))
A1 <- A1 + (1 - A1[1:nvec, C]) * runif(nvec)