R：如何从 4D 超椭圆外的几何空间中采样？答案

【问题标题】：R: How to sample from the geometric space outside a 4D hyper ellipse?R：如何从 4D 超椭圆外的几何空间中采样？
【发布时间】：2016-11-03 23:28:54
【问题描述】：

原标题（精确）：如何找到满足 4-unknown 非线性不等式的所有向量？

问题是找一个空间C，使得C中任意4x1向量x满足：

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

其中Q 是我已经知道的100 x 4 矩阵，a 是一个正常数。

我尝试从ellipse、rootSolve 和bvpSolve 等软件包中找到解决方案。但我无法找到合适的解决方案。

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

【问题讨论】：

标签： r ellipse

【解决方案1】：

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

我已更改您的问题标题。不可能列出所有解，尽管这样的空间具有简单的数学表示。我们可以从这样的空间中获得最好的样本。

从椭圆到球体的变换

这里有一些基于 Cholesky 分解的数学。或者，考虑对称特征分解，我在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) 的超球面。 如果我们可以从这样的空间中采样y，我们可以通过求解三角系统将其映射到x ：

x <- backsolve(R, y)

y的抽样

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

以下 R 函数从此类空间中采样 n y 向量。由于浮点数的有限表示，我们不能有无限半径，但最好是.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。

【讨论】：