0

I am looking for a general way to find the solution (in R) to the determined parts of an under-determined linear equation system, like the following one

# Let there be an A * x = b equation system
# with three coefficients A,B,C
# where only the coefficient A has a unique solution (A=2)
A <- rbind(c(0,1,1), c(1,0,0), c(0,0,0))
colnames(A) <- LETTERS[1:3]
b <- c(1,2,0)

cbind(A,b)
#      A B C b   
# [1,] 0 1 1 1 
# [2,] 1 0 0 2
# [3,] 0 0 0 0

I would like to solve for the parameter (A) that is determined and receive nothing for the under-determined parts, in this case should be

A = 2

Importantly, I am searching for a general way to determine the unique solutions that is not specific to the above example.

I have tried playing around with the QR-decomposition qr.coef(qr(A),b), which only shows me that C has no solution, but lacks the information that B has none.

I also played around with the single value decomposition svd(A) but the decomposition d in the result of the latter just indicates that one of the three parameters has a solution.

I am sure I am missing something obvious here -- thanks a bunch for the help!

JBJ
  • 866
  • 9
  • 21
  • Thanks for this - will check it out. – JBJ Nov 19 '20 at 15:36
  • I don't think this scales. Ideally, we would like to avoid solving the largest linear equation system only to throw away the solutions that vary. Is there a more efficient solution? – JBJ Nov 19 '20 at 15:40
  • That sounds like a different question to me. How big are you talking about? – user2554330 Nov 19 '20 at 15:46
  • this is for a software package in which users can supply linear constraints as big as they need, i.e. hard to say, unfortunately. – JBJ Nov 23 '20 at 10:36

1 Answers1

1

The question here Finding all solutions of a non-square linear system with infinitely many solutions is similar, but your problem is a little easier. You could use the SVD approach given in https://stackoverflow.com/a/46774174/2554330, but this QR approach is probably better:

Find the QR decomposition of A using decomp <- qr(A). Look in particular at qr.R(decomp). Since we know that

qr.R(decomp) %*% x = t(qr.Q(decomp)) %*% b

is solved by solutions to your original system, we need to look at which rows of qr.R(decomp) have exactly one non-zero entry: those will be uniquely determined by the above equation. But once those entries are determined, you can determine others by removing corresponding columns (i.e. once the first entry is known, maybe it could be used to solve for the second entry). I think this code does everything you want:

A <- rbind(c(0,1,1), c(1,0,0), c(0,0,0))
colnames(A) <- LETTERS[1:3]
b <- c(1,2,0)

decomp <- qr(A)
R <- qr.R(decomp)
R
#>       A B C
#> [1,] -1 0 0
#> [2,]  0 1 1
#> [3,]  0 0 0

# Identify the rows with solutions:
unsolved <- seq_len(ncol(A))
repeat {
  solved <- unsolved[apply(R[, unsolved, drop = FALSE] != 0, 1,
                           function(row) if (sum(row) == 1) which(row) else NA)]
  if (all(is.na(solved))) break
  unsolved <- setdiff(unsolved, solved)
}
solved <- setdiff(seq_len(ncol(A)), unsolved)
colnames(A)[solved]
#> [1] "A"

# Find the solutions:
Q <- qr.Q(decomp)
Qtb <- t(Q) %*% b

solve(R[solved, solved, drop = FALSE], Qtb[solved])
#> A 
#> 2

Created on 2020-11-19 by the reprex package (v0.3.0)

user2554330
  • 37,248
  • 4
  • 43
  • 90
  • Hi, thank you - seems perfect (I don't understand the QR decomposition well enough to get there, sad to say). The iterative nature of the algorithm though leads to somehow weird results in over-determined inconsistent equation systems such as `A <- cbind(a=c(1,1,0,0), b=c(0,1,0,1), c=c(0,0,1,0))`with `b <- c(1,0,0,0)` – JBJ Nov 20 '20 at 09:58
  • I believe in this case you're seeing the least squares solution (i.e. `x` with smallest sum of squares difference between `A %*% x` and `b`). I don't know what it would mean in a system where some variables were omitted and some were over-determined: *maybe* the least squares solution in those variables?? – user2554330 Nov 21 '20 at 11:42
  • I think ideally, an error of the sort of "cannot solve inconsistent constraints" would be reasonable for over-determined inconsistent systems – JBJ Nov 23 '20 at 09:22