QR factorization in R not giving correct answer with large sparse matrix

Question

I have a sparse matrix of type dgTMatrix built with the Matrix package in R and I'm trying to solve for the vector x in Ax=b using QR decomposition, but it is not working correctly.

For example below I am solving for a random A, b and you can see the method works with the norm of Ax-b = 0

A <- matrix(rnorm(9),ncol=3)
decomp <- qr(A)
b <- rnorm(3)
x <- qr.coef(decomp,b)
(A %*% x - b) %>% norm
[1] 3.885781e-16

My A and b are like this:

A: 173700 x 173700 sparse Matrix of class "dgTMatrix"``` 
b: 173700 x 1 sparse Matrix of class "dgCMatrix"```

When I take a QR decomposition of A I get the following:

qr_decomp <- qr(A, LAPACK = TRUE, tol = 1e-10)
qr_decomp
'MatrixFactorization' of Formal class 'sparseQR' [package "Matrix"] with 6 slots
  ..@ V   :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. ..@ i       : int [1:13350377] 0 173283 173284 173285 173286 1 2 3 2 3 ...
  .. .. ..@ p       : int [1:173701] 0 5 8 11 19 22 25 30 33 36 ...
  .. .. ..@ Dim     : int [1:2] 173700 173700
  .. .. ..@ Dimnames:List of 2
  .. .. .. ..$ : NULL
  .. .. .. ..$ : NULL
  .. .. ..@ x       : num [1:13350377] -9.27e+01 2.35e-03 -1.28e+02 2.35e-03 6.40e+01 ...
  .. .. ..@ factors : list()
  ..@ beta: num [1:173700] 6.88e-05 7.80e-14 1.24e-06 3.88e-05 7.80e-14 ...
  ..@ p   : int [1:173700] 34977 38066 38838 39610 38067 38839 34980 38069 38841 39613 ...
  ..@ R   :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. ..@ i       : int [1:24075359] 0 1 1 2 0 1 2 3 4 4 ...
  .. .. ..@ p       : int [1:173701] 0 1 2 4 8 9 11 12 13 15 ...
  .. .. ..@ Dim     : int [1:2] 173700 173700
  .. .. ..@ Dimnames:List of 2
  .. .. .. ..$ : NULL
  .. .. .. ..$ : NULL
  .. .. ..@ x       : num [1:24075359] 1.57e+02 3.58e+06 -1.26e+03 3.58e+06 1.92e-03 ...
  .. .. ..@ factors : list()
  ..@ q   : int [1:173700] 35749 38838 38066 36522 38839 38067 35752 38841 38069 36525 ...
  ..@ Dim : int [1:2] 173700 173700

According to the docs for sparseQR class,

For a sparse m×n (“long”: m≥n) rectangular matrix A, the sparse QR decomposition is either of the form PA=QR with a (row) permutation matrix P, (encoded in the p slot of the result) if the q slot is of length 0, or of the form PAP∗=QR with an extra (column) permutation matrix P∗ (encoded in the q slot).

Since qr_decomp@q exists, per the docs, this must be of the form PAP∗=QR

Anyway, the same procedure as above is not working in this case:

decomp <- qr(A)
x <- qr.coef(decomp, b)
(A %*% x - b) %>% norm
[1] 3.540814e+24

The above value is really far from zero.

Inspecting matrix A everything looks okay to me:

!(A@x %>% is.finite) %>% any
[1] FALSE
!(b@x %>% is.finite) %>% any
[1] FALSE

> A[1:10, 1:10]
10 x 10 sparse Matrix of class "dgTMatrix"
                                                                                                                                                 
 [1,] -5.366088e+08  5.301600e-02  .             .             .             .             .             .             .             .           
 [2,]  4.418000e-02 -5.366088e+08  4.418000e-02  .             .             .             .             .             .             .           
 [3,]  .             4.123467e-02 -5.366088e+08  4.123467e-02  .             .             .             .             .             .           
 [4,]  .             .             3.976200e-02 -5.366088e+08  3.976200e-02  .             .             .             .             .           
 [5,]  .             .             .             3.887840e-02 -5.366088e+08  3.887840e-02  .             .             .             .           
 [6,]  .             .             .             .             3.828933e-02 -5.366088e+08  3.828933e-02  .             .             .           
 [7,]  .             .             .             .             .             3.786857e-02 -5.366088e+08  3.786857e-02  .             .           
 [8,]  .             .             .             .             .             .             2.650800e-02 -5.366088e+08  1.251767e-02  .           
 [9,]  .             .             .             .             .             .             .             9.719600e-03 -5.366088e+08  9.719600e-03
[10,]  .             .             .             .             .             .             .             .             9.572333e-03 -5.366088e+08

However, there are large values present in the matrices:

max(b)
[1] 3.978441e+22
max(A)
[1] 3.979517e+22
min(b)
[1] 0
min(A)
[1] -7.958754e+22

Could it be the large values messing things up with some kind of round off error? Or am I missing some key info about QR not converging in some case?

Answer 1

tl;dr it does appear that large values will lead to problems; scaling the inputs appears to help.

An example that isn't problematic, even though I used Cauchy-distributed samples to fill in the matrix (which leads to a wide range of values):

set.seed(101)
library(Matrix)
## d <- 173700
d <- 1e3
## n <- 1e6
n <- 1e5
randfun <- rcauchy
As <- sparseMatrix(i=integer(0),
                   j=integer(0),
                   dims = c(d, d), repr = "T")
As[cbind(sample(d, size = n, replace = TRUE),
         sample(d, size = n, replace = TRUE))] <- randfun(n)
bs <- sparseMatrix(i = 1:d, j = rep(1, d), x = randfun(d))
qr_decomp <- qr(As, LAPACK = TRUE, tol = 1e-10)

decomp <- qr(As)
xs <- qr.coef(decomp, bs)
(As %*% xs - bs) |> norm()  ## 1.19e-09

If I use larger values for d (matrix/vector dimensions) and n (number of non-zero elements), either the computation takes a long time or I get a segmentation fault.

However, if I set randfun <- function(n) rcauchy(n, scale = 1e20) in the example above I do get large deviations.

Scaling the inputs seems to solve the problem, although I don't know if it will work for what you want to do ...

As2 <- As/1e20
bs2 <- bs/1e20
decomp2 <- qr(As2)
xs2 <- qr.coef(decomp2, bs2)
(As2 %*% xs2 - bs2) |> norm()