2

R can solve underdetermined linear systems:

A = matrix((1:12)^2,3,4,T)
B = 1:3
qr(A)$rank  # 3
qr.solve(A, B)  # solutions will have one zero, not necessarily the same one
# 0.1875 -0.5000  0.3125  0.0000
solve(qr(A, LAPACK = TRUE), B)
# 0.08333333 -0.18750000  0.00000000  0.10416667

(It gives one solution among the infinity of solutions).

However, if the rank (here 2) is lower than the number of rows (here 3), it won't work:

A = matrix(c((1:8)^2,0,0,0,0),3,4,T)
B = c(1,2,0)
A
#      [,1] [,2] [,3] [,4]
# [1,]    1    4    9   16
# [2,]   25   36   49   64
# [3,]    0    0    0    0

qr.solve(A, B)  # Error in qr.solve(A, B) : singular matrix
solve(qr(A, LAPACK = TRUE), B)  # Error in qr.coef(a, b) : error code 3

but this system does have a solution!

I know that the general solution is to use SVD or generalized/pseudo inverse of A (see this question and its answers), but:

Is there a mean with solve or qr.solve to automatically reduce the system AX=B to an equivalent system CX=D of only rank(A) rows, for which qr.solve(C, D) would simply work out-of-the-box?

Example: 

C = matrix(c((1:8)^2),2,4,T)
D = c(1,2)
qr.solve(C, D)
# -0.437500  0.359375  0.000000  0.000000
Basj
  • 41,386
  • 99
  • 383
  • 673

1 Answers1

3

qr.coef along with qr seem to do the job:

(A <- matrix(c((1:8)^2, 0, 0, 0, 0), nrow = 3, ncol = 4, byrow = TRUE))
#     [,1] [,2] [,3] [,4]
# [1,]    1    4    9   16
# [2,]   25   36   49   64
# [3,]    0    0    0    0
(B <- c(1, 2, 0))
# [1] 1 2 0
(X0 <- qr.coef(qr(A), B))
# [1] -0.437500  0.359375        NA        NA
X0[is.na(X0)] <- 0
X0
# [1] -0.437500  0.359375  0.000000  0.000000
# Verification:
A %*% X0
#      [,1]
# [1,]    1
# [2,]    2
# [3,]    0

Second example:

(A<-matrix(c(1, 2, 0, 0, 1, 2, 0, 0, 1, 2, 1, 0), nrow = 3, ncol = 4, byrow = TRUE))
#      [,1] [,2] [,3] [,4]
# [1,]    1    2    0    0
# [2,]    1    2    0    0
# [3,]    1    2    1    0
(B<-c(1, 1, 2))
# [1] 1 1 2
qr.solve(A, B)
# Error in qr.solve(A, B) : singular matrix 'a' in solve
(X0 <- qr.coef(qr(A), B))
# [1]  1 NA  1 NA
X0[is.na(X0)] <- 0
X0
# [1] 1 0 1 0
A %*% X0
#      [,1]
# [1,]    1
# [2,]    1
# [3,]    2
Julius Vainora
  • 47,421
  • 9
  • 90
  • 102
  • Yes but here you have to manually set the lines to keep (1:2). What if the system has 100 lines and the rank is 70, i.e. 30 (maybe not-consecutive) rows have to be discarded? (Or maybe I missed something in your example). – Basj Oct 26 '18 at 13:28
  • To reformulate my question @JuliusVainora: how do you create a solution X0 such that A X0 = B, without specifiying which rows should be kept and which ones should be discarded? (the "1:2" in your code) – Basj Oct 26 '18 at 13:29
  • not paying too much attention, but does `w <- !is.na(coef); A[,w] %*% na.omit(coef)` work? – Ben Bolker Oct 26 '18 at 13:35
  • 1
    It's not necessary to know 1:2 in `qr.coef(qr(A), B)`, which is all your question is about. I only added a verification step. As to get rid of `NA`, you may use `coef[!is.na(coef)]`, or what @BenBolker suggests for the verification part. – Julius Vainora Oct 26 '18 at 13:40
  • @JuliusVainora What does `qr.coef(qr(A), B)` represent exactly? The question was about getting an equivalent system `CX=D` with only two rows in this example. How to go from `qr.coef(qr(A), B)` to the matrices `C`, `D`? I see we are quite close to it, but I just want to be sure to use the best method. – Basj Oct 26 '18 at 14:27
  • @Basj, truth is that I also don't know the details; looking at `qr.coef` leads to Fortran code. See the update, though. Looks like `qr.coef(qr(A), B)` does what you want (?): solves the system with your known methods after eliminating redundant rows, so that the rank coincides with the number of rows. Then it would seem like those `NA` could be just zeros, as in the output of `qr.solve`. I hope that helps. – Julius Vainora Oct 26 '18 at 15:21
  • Being curious @JuliusVainora: why the additional `(...)` around `X0 <- ...`? – Basj Oct 26 '18 at 15:51
  • `X0 <- ...` would be just an assignment, while adding parenthesis gives as an output the newly assigned value. That is, `(X <- 1)` works as `X <- 1; X`. – Julius Vainora Oct 26 '18 at 15:54
  • @JuliusVainora I have slightly modified the answer, and now it 100% works, I hope it's ok for you (I accepted the answer). I also added a second example. – Basj Oct 27 '18 at 07:03