Original title (in-precise): How to find all vectors that satisfy a 4-unknown non-linear inequalities?
The question is to find a space C, so that any 4x1 vector x in C satisfies:
t(x) %*% t(Q) %*% Q %*% x > a,
in which Q is a 100 x 4 matrix I already know, and a is a positive constant.
I tried to find the solution from packages such as ellipse, rootSolve, and bvpSolve. But I can't reach a suitable solution.
Any idea or solution will be sincerely appreciated.
Remark: the algorithm below can be used to sample from the surface / manifold of a 4D hyper ellipse, or its inner space, too.
I have changed your question title. It is impossible to list all solutions, although such space has simple mathematical representation. We can at our best sample from such space.
Transformation from Ellipse to Sphere
Here is some mathematics based on Cholesky factorization. Alternatively, consider symmetric Eigen decomposition, and I have a demonstration / comparison between these two at Obtain vertices of the ellipse on an ellipse covariance plot (created by car::ellipse), with nice figures for the geometry.
Since Q is known, R is known. The following R code gets R:
R <- chol(crossprod(Q))
y is from the hypersphere with radius greater than sqrt(a). If we can sample y from such space, we can then map it to x by solving a triangular system:
x <- backsolve(R, y)
Sampling of y
We can use n-sphere coordinates to parametrize such space. For 4D space, we have:
The following R function samples n y-vectors from such space. Because of finite representation of floating point numbers, we can't have infinity radius but .Machine$double.xmax at the best. But we also use an optional argument rmax, if we want a more restricted radius.
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))
}
Try some examples:
## 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
I have chosen to use each column than each row as a sample, to make matrix computation easier later.
Transforming y to x
Now assume we have
set.seed(0); Q <- matrix(runif(10 * 4), 10L, 4L)
We get 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
Suppose you have a = 4, then we map y to 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
Checking
We can check that the above x does satisfy our requirement.
z <- Q %*% x
ax <- colSums(x ^ 2) ## value of `diag(x'Q'Qx)`
#[1] 84.15453 17.00795 24.77044 43.33361 85.85250
They are all greater than 4.
