1

I have a under-determined equation system (Ax = b) with infinite solutions. My goal is to get n random solutions (i.e n vectors x) for this system. Is it possible to do it with R (if possible with base R functions) ?

With base::qr.solve() you can get 1 solution but it is always the same.

Here is a reproducible example to work with:

# Ex data illustrating this situation: Ax = b
A = matrix(
c(rep(1,11), 0.1803, 0.0071, 0.0063, 0.0201, 0.3333, 0.0043, 0.1573, 0.0007, 0.2439, 0.0072, 0.0011),
nrow = 2,
ncol = 11,
byrow = TRUE
)

b = as.matrix(c(1,0.1))

# currently not random x solutions
x = qr.solve(A, b)

Current result:

print(x)
           [,1]
 [1,] 0.5363741
 [2,] 0.4636259
 [3,] 0.0000000
 [4,] 0.0000000
 [5,] 0.0000000
 [6,] 0.0000000
 [7,] 0.0000000
 [8,] 0.0000000
 [9,] 0.0000000
[10,] 0.0000000
[11,] 0.0000000

Here are some related questions:

Paul
  • 2,850
  • 1
  • 12
  • 37

1 Answers1

2

According to the first link you provided you can get the general solution (using MASS): 

library(MASS)

Ag <- ginv(A)
Ag
xb <- Ag %*% b
xb
Aw <- diag(nrow=nrow(Ag)) - Ag %*% A
Aw

and thus you can create n solutions in the form:

w <- runif(11)
z <- xb + Aw %*% w

like so (example matrix of 100 solutions):

sapply(1:100, function(i) {
  w <- runif(11)
  z <- xb + Aw %*% w}
)
Jeroen Colin
  • 345
  • 1
  • 6
  • Many thanks for this tip, I am not familiar with these functions your answer helped a lot – Paul Nov 25 '19 at 10:43