R：如何从4D超椭圆外的几何space采样？

Question

原标题(in-precise): 如何找到所有满足4-未知non-linear不等式的向量？

题目是求一个space C，使得C中任意一个4x1的向量x满足：

t(x) %*% t(Q) %*% Q %*% x > a,

其中Q是我已知的100 x 4矩阵，a是正常数

我试图从 ellipse、rootSolve 和 bvpSolve 等包中找到解决方案。但是我找不到合适的解决方案。

任何想法或解决方案将不胜感激。

Answer 1

备注：下面的算法也可用于从4D超椭圆的表面/流形或其内部space采样。

我已经更改了你的问题标题。不可能列出所有解决方案，尽管 space 具有简单的数学表示。我们可以从这样的 space.

中获取最好的样本

椭圆到球体的变换

这里是一些基于 Cholesky 分解的数学。或者，考虑对称 Eigen 分解，我在 Obtain vertices of the ellipse on an ellipse covariance plot (created by car::ellipse) 上对这两者进行了演示/比较，并提供了很好的几何图形。

既然Q是已知的，那么R也是已知的。以下 R 代码得到 R:

R <- chol(crossprod(Q))

y来自半径大于sqrt(a)的超球体。如果我们可以从这样的space中采样y ，然后我们可以通过求解三角系统将其映射到 x：

x <- backsolve(R, y)

y

的抽样

我们可以使用n-sphere coordinates来参数化这样的space。对于 4D space，我们有：

以下 R 函数从这样的 space 中采样 n y-vectors。由于浮点数的有限表示，我们不能有无穷大的半径，但最多 .Machine$double.xmax。但如果我们想要更受限的半径，我们也使用可选参数 rmax。

ry <- function (n, rmin, rmax = NA) {
  if (is.na(rmax)) rmax <- .Machine$double.xmax
  if (rmin > rmax) stop("larger `rmax` expected!")
  r <- runif(n, rmin, rmax)
  phi1 <- runif(n, 0, pi)
  phi2 <- runif(n, 0, 2 * pi)
  phi3 <- runif(n, 0, 2 * pi)
  matrix(c(r * cos(phi1), 
           r * sin(phi1) * cos(phi2),
           r * sin(phi1) * sin(phi2) * cos(phi3),
           r * sin(phi1) * sin(phi2) * sin(phi3)),
         nrow = 4L, byrow = TRUE, dimnames = list(paste0("y", 1:4), NULL))
  }

尝试一些例子：

## radius between 4 and 10
set.seed(0); ry(5, 4, 10)
#         [,1]        [,2]       [,3]       [,4]      [,5]
#y1  7.5594886 -5.31049687 -6.1388372 -3.5991830 -3.728597
#y2  5.1402945  0.47936481  0.4799181 -2.5085948 -6.480402
#y3  0.2614002 -1.68833263 -0.1950092 -5.9975328 -4.213166
#y4 -2.0859078  0.02440839 -0.9452077  0.3052708  3.954674

## radius between 4 and "inf"
set.seed(0); ry(5, 4)
#             [,1]           [,2]           [,3]           [,4]           [,5]
#y1  1.299100e+308 -4.531902e+307 -6.588856e+307 -4.983772e+307 -6.442420e+307
#y2  8.833607e+307  4.090831e+306  5.150993e+306 -3.473640e+307 -1.119710e+308
#y3  4.492167e+306 -1.440799e+307 -2.093047e+306 -8.304756e+307 -7.279678e+307
#y4 -3.584637e+307  2.082977e+305 -1.014498e+307  4.227070e+306  6.833046e+307

我选择使用每一列而不是每一行作为样本，以便以后更容易进行矩阵计算。

将y转换为x

现在假设我们有

set.seed(0); Q <- matrix(runif(10 * 4), 10L, 4L)

我们得到R

R <- chol(crossprod(Q))
#         [,1]     [,2]      [,3]      [,4]
#[1,] 2.176848 1.420882 1.2517326 1.4481875
#[2,] 0.000000 1.077816 0.1045581 0.4646328
#[3,] 0.000000 0.000000 1.2284251 0.3961126
#[4,] 0.000000 0.000000 0.0000000 0.9019771

假设你有a = 4，那么我们将y映射到x:

a <- 4
set.seed(0); y <- ry(5, sqrt(a), 10)  ## we set an `rmax` here
x <- backsolve(R, y)  ## each column is a sample of `x`
#           [,1]        [,2]       [,3]       [,4]      [,5]
#[1,]  0.7403534 -1.49866534 -2.2359350  2.0269516  2.948561
#[2,]  5.5481682  0.41827601  0.7024109 -1.7606720 -7.288339
#[3,]  0.9373905 -1.01984708  0.1430660 -4.4180688 -4.749419
#[4,] -2.2616584  0.01995357 -0.8367956  0.2995693  4.299268

正在检查

我们可以检查上面的x是否满足我们的要求。

z <- Q %*% x
ax <- colSums(x ^ 2)  ## value of `diag(x'Q'Qx)`
#[1] 84.15453 17.00795 24.77044 43.33361 85.85250

都大于4。

R：如何从4D超椭圆外的几何space采样？

R: How to sample from the geometric space outside a 4D hyper ellipse?

r

ellipse