我想你可以从一个具有正对角线的下三角形L开始构造一个正定矩阵M，例如M = L*L^T，使得如果x>0，则x^T*M*x = x^T*L*L^T*x = |L^T*x|^2应该是肯定的。
例如，实现方式：

n <- 9
L <- diag(n)
L[-1, 1] <- c(76:79, 81:84) / 100
v <- 1 - L[, 1]^2
for (k in 2:length(v)) {
  L[k, 2:k] <- sqrt((r <- runif(k - 1)) / sum(r) * v[k])
}
M <- tcrossprod(L)

我们可以看到

> M
      [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
 [1,] 1.00 0.7600000 0.7700000 0.7800000 0.7900000 0.8100000 0.8200000
 [2,] 0.76 1.0000000 0.8883572 0.7988221 0.7328434 0.7003986 0.7966054
 [3,] 0.77 0.8883572 1.0000000 0.8972778 0.7937081 0.8500339 0.8200438
 [4,] 0.78 0.7988221 0.8972778 1.0000000 0.9093460 0.8651427 0.8668928
 [5,] 0.79 0.7328434 0.7937081 0.9093460 1.0000000 0.8581035 0.9336535
 [6,] 0.81 0.7003986 0.8500339 0.8651427 0.8581035 1.0000000 0.9105348
 [7,] 0.82 0.7966054 0.8200438 0.8668928 0.9336535 0.9105348 1.0000000
 [8,] 0.83 0.6964236 0.7789939 0.7994304 0.8621129 0.9136291 0.9507678
 [9,] 0.84 0.7344271 0.8091267 0.7860827 0.8387232 0.9097534 0.9200662
           [,8]      [,9]
 [1,] 0.8300000 0.8400000
 [2,] 0.6964236 0.7344271
 [3,] 0.7789939 0.8091267
 [4,] 0.7994304 0.7860827
 [5,] 0.8621129 0.8387232
 [6,] 0.9136291 0.9097534
 [7,] 0.9507678 0.9200662
 [8,] 1.0000000 0.9652907
 [9,] 0.9652907 1.0000000

我们可以检查所有的特征值都大于0

> all(eigen(M)$values > 0)
[1] TRUE

可能有一些优雅的线性代数的作品，但蛮力似乎可以解决这个问题。

## set up matrix with required form
m <- matrix(NA, 9, 9)
diag(m) <- 1
m[1,-1] <- m[-1, 1] <- c(76:79, 81:84)/100

mk_matrix <- function(p) {
    ## fill in lower triangle
    m[col(m) > 1 & row(m) > col(m)] <- p
    ## symmetrize
    m[upper.tri(m)] <- t(m)[upper.tri(m)]
    m
}

objfun <- function(p = rnorm(8*7/2)) {
    r <- -1*det(mk_matrix(p))
    ## cat(r, "\n")
    r
}

objfun的正值表示负定矩阵。
似乎我们有50/50的机会通过填充随机的正常值来获得pos-def矩阵：

set.seed(101)
rr <- replicate(1000, objfun())
hist(rr, breaks = 50)

在任何情况下，我们都可以使用optim()和abstol，只要我们得到一个负值（==一个pos def结果），就立即停止。

set.seed(104) ## chosen to get a neg def starting matrix
optim(par = rnorm(8*7/2), 
      fn = objfun, control=list(abstol = 1e-8))

如果你真的不关心其他相关性，你可以直接在线性时间内提取样本（你可以对结果应用逐点线性变换而不改变相关性）。

correlations <- c(1, 0.76, 0.77, 0.78, 0.79, 0.81, 0.82, 0.83, 0.84)
n <- length(correlations)
samples <- correlations * rnorm(1) + sqrt(1 - correlations^2) * rnorm(n)